Numerical Solution of Stiff Multibody Dynamic Systems Based on Kinematic Derivatives

The Hermite–Obreshkov–Padé (HOP) method of numerical integration is applicable to stiff systems of differential equations, where the linearization has large range of eigenvalues. A practical implementation of HOP requires the ability to determine high-order time derivatives of the system variables. In the case of a constrained multibody dynamical system, the power series solution for the kinematic differential equation is the foundation for an algorithmic differentiation (AD) procedure determining those derivatives. The AD procedure is extended in this paper to determine rates of change in the time derivatives with respect to variation in the position and velocity state variables of the multibody system. The coefficients of this variation form the Jacobian matrix required for Newton–Raphson iteration. That procedure solves the implicit relations for the state variables at the end of each integration time step. The resulting numerical method is applied to the rotation of a dynamically unbalanced constant-velocity (CV) shaft coupling, where the deflection angle of the output shaft is constrained to low levels by springs of high rate and damping.