## 1 Introduction  
Multibody system dynamics are widely used in aerospace and mechanical engineering, as they can predict the dynamic response of structures and their components which undergo large deformations during long-period simulation times. In contrast to the conventional rigid-body dynamics, which neglect elastic deformation, multibody system dynamics focus on developing nonlinear transient solutions, in which large rigid-body motion and elastic deformation occur simultaneously. Thus, it can predict system dynamics with unprecedented accuracy [1]. Geometric nonlinearity needs to be addressed in solving transient problems of flexible structures with finite rotations and large displacements. Besides rigid body motion, elastic deformation often appears during analysis. An accurate measurement on the strain is required as strain energy produced by the elastic deformation is usually not negligible. Therefore, the partition of total deformation into elastic deformation and rigid body motion needs to be studied. Due to the complexity of the nature of the problem, together with the various complex shapes of structures, analytical solutions to flexible multibody dynamics are usually not available. Application engineers must use numerical methods to explore reliable solutions. The combination of the Finite Element Method (FEM) with time integration algorithms is the predominant solution [2]. When using FEM to solve a flexible multibody system, kinematic equations of flexible bodies and formulation of global dynamic equations are the two challenges [3]. Three approaches to model flexible bodies can be identified according to the selection of the coordinate system: floating frame formulation, inertial frame formulation, and co-rotational frame formulation. Despite the equivalence in the dynamic governing equation derivation, the different selections of coordinate systems show remarkable differences in computational performance. In this paper, the co-rotational method is adopted. The success of the co-rotational method lies in the explicit incorporation of geometric nonlinearity stemming from rigid-body motions. Compared to the other two approaches, the co-rotational method sets up a local coordinate system in each element at the beginning of analysis, and they are consistently updated during the analysis. The local coordinate system needs a proper selection so that the element deformation at each step is decomposed into a rigid-body motion and an elastic deformation. The elimination of rigid-body motion is achieved by a transformation matrix. The transformation matrix is formulated to connect the local and the global coordinate systems. When constructing element tangent stiffness matrices and internal force vectors, the co-rotational method naturally eliminates the parasitic coupling effects induced by rigid-body motion through location and orientation updates of the local system. Thus, a balance between computational efficiency and numerical accuracy is achieved when formulating the element kinematic equations. In contrast to static problems, in the study of flexible multibody systems, a robust time stepping algorithm is demanded to accurately simulate the transient response of the flexible bodies. The generalized-\(\alpha\) method and high-order backward differentiation formulas [4] are widely used as time discretization schemes due to their robustness and excellent performance in linear transient analyses. However, these algorithms exhibit numerical instability in nonlinear problems [5]. Energy conservation is a criterion in stability assessment for nonlinear time-history integration methods. Following this principle, researchers have pioneered approaches to derive governing equations based on energy-conserving principles since the 1960s. In 1978, [PERSON] et al. [6] developed the Constraint Energy Method by incorporating Lagrange multipliers into the trapezoidal rule, yet this algorithm fails to preserve momentum conservation.