Abstract: The flexible multibody system dynamics issues have high geometrically nonlinear, which is not only due to the large deformation of the flexible components but also caused by the extensive rigid body motion. Recently, the modeling method based on local frame of Lie group (LFLG) has been verified that it can be combined with various modeling methods to eliminate geometric nonlinearity caused by rigid body motion. At the same time, the geometric nonlinearity caused by large deformation will gradually weaken with spatial discrete encryption. Therefore, LFLG can eliminate the geometric nonlinearity of the components in the multi-flexible system, and the generalized inertial forces and internal forces as well as their Jacobian matrices are invariable under the arbitrary rigid body motion. However, due to the nonlinearity of constraints and loads in multibody systems, whether LFLG method can improve the overall Jacobian matrix reuse efficiency in practical applications is still to be discussed. In addition, the time integral strategy of variable order and variable step size is usually used to solve the dynamic equation of multibody flexible systems. The change of algorithm order and step size will also lead to Jacobian change of multibody systems, which intensifies the Jacobian reuse difficulty of the system. In order to analyze the numerical characteristics of the LFLG method objectively in practical simulation, this paper firstly presents a 3D corotational Timoshenko beam element based on the Lie group local frame, which can reuse the element algorithm of small deformation finite element method to the maximum, and reduce the difficulty of element development compared with geometrically exact beam formulation and absolute nodal coordinate formulation. Secondly, the calculation flow of backward difference formula (BDF) with variable step size and order as well as generalized α integrator with variable step size of LFLG method is set up, and the geometric nonlinear characteristics of the element elastic force and damping force are analyzed for small deformation and large deformation conditions, and the Jacobian reuse efficiency is compared. Finally, the numerical characteristics of LFLG method and traditional modeling method in the global frame are analyzed.