Variable Step Implicit Numerical Integration of Stiff Multibody Systems

This paper presents a variable step size implicit numerical integration algorithm for dynamic analysis of stiff multibody systems. Stiff problems are very common in real world applications, and their numerical treatment by means of explicit integration is cumbersome or infeasible. Until recently, implicit numerical integration of the equations of motion of stiff mechanical systems has been problematic. Prior attempts to solve this problem stopped at the level of theoretical considerations, or resulted in numerical methods that lacked generality and robustness. A better theoretical understanding of the process of state space reduction, along with the ability to efficiently generate the needed derivative information for implicit integration, has allowed for the development of robust and general implicit numerical methods. The benefit of variable step size implicit integration is shown using a stiff double pendulum model. An experimental general purpose code that performs implicit integration is capable to successfully treat a large class of rigid body models. 1 The DAE of Multibody Dynamics Dynamic analysis of a constrained mechanical system model requires the numerical integration of a set of differential-algebraic equations (DAE). This type of problem is known to exhibit numerical difficulties, and it is prone to intense computation effort. There are many specific methods for the numerical solution of DAE of multibody dynamics, some of them outlined in greater detail by Potra (1994) and Haug et al. (1997, b). In this paper, the numerical solution of the differential-algebraic problem is obtained by reducing it to a set of ordinary differential equations (ODE). The reduction is based on state space parametrization, the evolution of the mechanical system being described by a number of independent variables equal to the number of degrees of freedom of the mechanical system. The existing theory for the numerical solution of ODE can then be employed for the solution of the resulting state space ordinary differential equations (SSODE). Once the solution of the SSODE is available, the homeomorphism induced by the parametrization considered is used to recover the solution of the differential-algebraic problem. If q = [q , q , , q ] 1 2 n T L is the vector of generalized coordinates used to model a mechanical system subject to the constraints in Eq. 1, the parametrization considered in this paper is defined by a subset of the generalized coordinates. They are called independent generalized coordinates, and are denoted by v ∈R , where ndof n m = − . Φ( ) [ ( ), ( ), , ( ) q q q q 0 ≡ = Φ Φ Φ 1 2 m T ] K (1) Details of the partitioning of q into independent and dependent vectors v and u ∈R, respectively, are provided by Haug et al. (1997,b). Differentiating Eq. 1 with respect to time yields the kinematic velocity equation,