A constraint elimination-based harmonic balance method for multibody systems

● A constraint elimination–based harmonic balance (CE-HB) method is proposed ● Lagrange multipliers are eliminated to reduce unknowns and improve convergence ● Two-step matrix reshaping reveals FFT-based Jacobian structure for nonlinear terms ● The proposed method reduces total computational cost by about 40% The harmonic balance (HB) semi-analytical method has shown strong potential for periodic response analysis of nonlinear dynamical systems. However, its direct application to multibody dynamics (MBD) equations formulated as differential-algebraic equations (DAEs) with Lagrange multipliers leads to increased computational costs, convergence difficulties, and strong sensitivity to initial conditions. In this study, a constraint elimination–based harmonic balance (CE-HB) method is proposed for the periodic solution of DAE-based MBD systems. The primary contribution lies in the analytical elimination of Lagrange multipliers prior to the HB formulation, which transforms the original DAEs into an equivalent ordinary differential equation (ODE) system. This transformation reduces the number of unknowns, avoids convergence issues associated with Lagrange multipliers, and enables the direct application of Floquet–Lyapunov stability analysis and arc-length continuation techniques within a unified ODE framework. For the Newton–Raphson iteration, a fast Fourier transform-based Jacobian construction is employed, and a two-step matrix reshaping procedure is presented to clarify and regularize the Jacobian formulation for nonlinear terms, which improves both numerical consistency and implementation transparency. Several numerical examples are provided to validate the proposed method. The results show excellent agreement with Runge–Kutta time integration and capture both stable and unstable periodic solutions. Compared with the direct method that considers differential and algebraic equations simultaneously, the proposed CE-HB method achieves a reduction in computational costs of over 40%. Moreover, it demonstrates improved convergence robustness, particularly with respect to the choice of initial conditions.