New Insights into a Three-Sub-Step Composite Method and Its Performance on Multibody Systems

: This paper develops a new implicit solution procedure for multibody systems based on a three-sub-step composite method, named TTBIF (trapezoidal–trapezoidal backward interpolation formula). The TTBIF is second-order accurate, and the effective stiffness matrices of the ﬁrst two sub-steps are the same. In this work, the algorithmic parameters of the TTBIF are further optimized to minimize its local truncation error. Theoretical analysis shows that for both undamped and damped systems, this optimized TTBIF is unconditionally stable, controllably dissipative, third-order accurate, and has no overshoots. Additionally, the effective stiffness matrices of all three sub-steps are the same, leading to the effective stiffness matrix being factorized only once in a step for linear systems. Then, the implementation procedure of the present optimized TTBIF for multibody systems is presented, in which the position constraint equation is strictly satisfied. The advantages in accuracy, stability, and energy conservation of the optimized TTBIF are validated by some benchmark multibody dynamic problems.


Introduction
For solving dynamic problems after spatial discretization, time integration methods are a powerful numerical tool. The first one was originally developed for structural dynamic systems and was then generalized to multibody systems. Structural dynamic systems are governed by ordinary differential equations (ODEs), whereas multibody systems are typically formulated as a set of differential-algebraic equations (DAEs) [1]. The introduction of constraint equations results in the numerical treatment of DAEs being more challenging than that of ODEs [2]. By differentiating position constraints with respect to time, the DAEs can be transformed into ODEs; hence, time integration methods developed for structural dynamic systems can be directly used to solve multibody systems. However, this procedure is costly [2] and produces numerical drifts [3][4][5]. Therefore, in the simulation of multibody systems, the direct discretization of DAEs has gained more attention.
Generally, conventional time integration methods [6,7] are divided into two categories, explicit and implicit methods [8]. Explicit methods [8][9][10] are not suitable for solving DAEs directly because they cannot satisfy constraint equations at the position level. In contrast, implicit methods together with an iterative procedure (e.g., Newton-Raphson method) are applicable. Therefore, the well-known trapezoidal rule (TR), also denominated Newmark's constant average acceleration method [8], the backward differentiation formula (BDF) [11], the generalized-α method [12,13], and the Hilber-Hughes-Taylor-α method (HHT-α) [14] have been employed in the solutions of DAEs. The non-dissipative TR can keep the energy of systems, but it cannot filter out unwanted or spurious high-frequency information [15]. The asymptotic annihilating (or L-stable) BDF is particularly useful for stiff problems, whereas it cannot automatically start [11]. The controllably dissipative (or A-stable) generalized-α method and HHT-α method are more popular in multibody systems [16][17][18][19][20]. However, they both use the time weighted residual representation of motion equations, resulting in their acceleration being first-order accurate [21].
To our knowledge, there are many excellent implicit methods for structural dynamics, such as the single-step parameters methods [21,22], energy momentum methods [23][24][25], linear multistep methods [26][27][28][29], and composite methods [30][31][32][33][34][35][36][37][38][39][40]. Using them to simulate multibody systems seems natural. Among these methods, the self-starting composite methods possessing advantages in accuracy, efficiency, dissipation, and stability have become more attractive in recent decades. The composite methods divide a time step into several sub-steps, in which different methods can be adopted. The L-stable Bathe method [30,31] proposed in 2005 is a representative composite method. In the Bathe method, the TR is used in the first sub-step to ensure low-frequency accuracy, and the second sub-step adopts the BDF to provide numerical dissipation in the high-frequency range. Adopting the same idea, some more accurate composite methods [32][33][34][35] with L-stability, such as the TR-TR-BDF (TTBDF) [32], the Wen method [33], and the TR-TR-TR-BDF (TTTBDF) [34], were proposed. For flexibly controlling the amount of dissipation, the controllably dissipative composite methods [36][37][38][39][40] based on the TR and the backward interpolation formula (BIF) were developed, such as the ρ ∞ -Bathe method [37], the Kim method [38], and the TR-TR-BIF (TTBIF) [39]. Among them, the three-sub-step TTBIF [39] proposed by the present authors has been proved to have excellent performance in linear and nonlinear structural dynamics, and it is superior to the ρ ∞ -Bathe method and the Kim method in some aspects. Additionally, the TTBIF has been employed in seismic response analysis by other researchers [41,42], further showing its superiority.
In this context, the three-sub-step TTBIF [39] is further optimized and generalized to multibody systems here. The algorithmic parameters which can minimize local truncation errors are found in this work, and the numerical properties, including spectral characteristics, overshoot characteristics, and convergent rates of the present optimized TTBIF for both undamped and damped systems are deliberately investigated. Then, the solution procedure for multibody systems of the present TTBIF is given, and its validation is demonstrated by some benchmark multibody dynamic problems. This paper is organized as follows. In Section 2, the governing equation of the multibody system, the updating formulations, and the optimized parameters of the TTBIF are given. The solution procedure of the present TTBIF for multibody systems is provided in Section 3. Then, the numerical properties of the present TTBIF are analyzed in Section 4. The numerical experiments are presented in Section 5. Finally, the conclusions are drawn in Section 6.

Formulation
This paper develops an effective implicit solution procedure for multibody systems. In this section, we first describe the governing equations of multibody systems, and then provide the updating equations and optimized parameters of the present TTBIF.

Governing Equations
The governing equations of multibody systems, in the absence of nonlinear hysteretic forces [43][44][45], can be described as where M is the system's mass matrix; q is the generalized coordinate vector, and the over dot represents its derivative with respect to time t; Q is the generalized force vector; Φ is the position constraint equation; Φ q (q, t) = ∂ Φ(q, t)/∂q is the constraint's Jacobian matrix; λ is the Lagrange multiplier vector. Equation (1) can be directly solved by the implicit methods, such as the TR, BDF, and generalized-α method. By differentiating the position constraints Φ with respect to time twice, the constraint equation given in Equation (1) Then, Equation (1) can be reformulated as where ,t = ∂ /∂t, ,q = ∂ /∂q, ,qt = ∂ 2 /∂q∂t, and ,tt =∂ 2 /∂t 2 . Most time integration methods developed for structural dynamic systems can be used to solve Equation (4). However, transforming DAEs (1) into ODEs (4) causes violations of the constraints in position and velocity, resulting in unacceptable global results. For this reason, the TTBIF given in Section 2.2 is directly used to deal with DAEs (1) in this paper.

Updating Equations of the TTBIF
The TTBIF consists of three sub-steps, [t, t + γ 1 ∆t], [t + γ 1 ∆t, t + γ 2 ∆t], and [t + γ 2 ∆t, t + ∆t], in which γ 1 and γ 2 , respectively, stand for the ratios of the first sub-step size and second sub-step size to the entire step size ∆t. In the TTBIF, the first two sub-steps adopt TR for preserving as much low-frequency content as possible. The updating equations of the first-sub-step are from which the state vectors, q t+γ 1 ∆t , . q t+γ 1 ∆t , .. q t+γ 1 ∆t , and λ t+γ 1 ∆t , at the moment t + γ 1 ∆t can be solved, and they are known information in the solution of the second sub-step. The displacement, velocity, acceleration, and Lagrange multiplier at the moment t + γ 2 ∆t are updated by The four-point BIF is employed in the last sub-step, and its updating equations have the forms One can find from Equations (5)-(7) that the TTBIF has six algorithmic parameters, θ 0 , θ 1 , θ 2 , θ 3 , γ 1 , and γ 2 ; and the values of them are closely related to the numerical properties of the TTBIF. These parameters were designed in the previous work [39] for satisfying the basic requirements, including second-order accuracy, unconditional stability, and controllable dissipation, and they were written as the functions of ρ ∞ (the spectral radius at the frequency limit [8]), as follows: where Here are some remarks for the above algorithmic parameters. Equation (8) ensures that the TTBIF is unconditionally stable for undamped systems. The TTBIF has the ability to control the degree of numerical dissipation when Equation (9) holds. When Equations (10) and (11) hold, the TTBIF is second-order accurate in displacement, velocity, and acceleration. Additionally, the value of γ 1 is required to be located in the range given in Equation (15) to ensure that the value of c 3 is positive (or avoids the value of θ 0 is imaginary), and the value of γ 2 is assumed to be twice that of γ 1 , leading to the effective stiffness matrices of the first two sub-steps being the same.

Further Optimization of the TTBIF
From Equations (8)- (16), one can find that the numerical properties of the TTBIF closely depend on the value of γ 1 . In this section, the value of γ 1 is optimized to minimize the local truncation errors [8] which can be used to evaluate the numerical accuracy of a time integration method. It can be seen in Equations (5)-(7) that the difference formulas of displacement and velocity are the same in the TTBIF; therefore, Equations (5)-(7) can be reformulated as where y T = [q T . q T ]. Generally, the numerical properties of time integration methods are evaluated by the single degree-of-freedom equation .. q + 2ξω . q + ω 2 q = 0 (ξ and ω stand for the damping ratio and natural frequency, respectively) [8], which is common and sufficient owing to the mode superposition principle. The second-order equation q + ω 2 q = 0 can be equivalently rewritten as the following first-order equation [27]: Thus, the above equivalent model (18) is used in the optimization of algorithmic parameters to simplify the analysis. Applying Equation (17) to Equation (18) yields the following recursive equation: where A represents the amplification factor, and has the form as The detailed derivations of recursion Equations (19) and (20) can be found in Appendix A. For discussing the accuracy of the TTBIF, the local truncation error σ [8] is introduced here.
If σ = O(τ n+1 ), the method is said to be nth-order accurate, which requires that up to nth-derivatives of A at τ = 0 are all equal to 1; i.e., It can be found from Equation (20) together with Equations (10) and (11) that Thus, the TTBIF is at least second-order accurate, and when A (3) (0) = 1 holds, the TTBIF can obtain third-order accuracy for linear systems. The changes of A (3) (0) with γ 1 are plotted in Figure 1, where in one can find that: (2) for the case of , the TTBIF can be third-order accurate, and the corresponding optimal γ 1 can be determined by A (3) (0) = 1; (3) additionally, we found that γ 1 −2θ 3 = 0 holds when the optimal γ 1 obtained by ϕ(γ 1 , ρ ∞ ) = 0 is adopted; then, the effective stiffness matrices of all three sub-steps are the same. This is important for linear ODEs, leading to the effective stiffness matrix being factorized only once in three sub-steps, as in single-step methods; refer to Appendix B.
The expressions of the optimal γ 1 are complicated, so its numerical values obtained by the bisection method are listed in Tables 1 and 2, ensuring that (γ 1 −2θ 3 ) is almost zero, as γ 1 −2θ 3 ≈ 10 −16 . Change of (A (1)

Implementation
For avoiding the numerical drifts, the multibody systems governed by DAEs (1) are directly solved by the TTBIF, and its procedure of one loop t → t + γ 1 ∆t → t + γ 2 ∆t → t+ ∆t is shown in Figure 2. In calculations, the increments x = [∆ .. q T ∆λ T ] T are calculated by the Newton-Raphson technique, as follows: where t = t + γ 1 ∆t, t + γ 2 ∆t, or t + ∆t, and the vector G = [G T 1 G T 2 ] T has the forms where α = γ 1 ∆t/2 (or θ 3 ∆t). Then, the Jacobian matrix J can be derived from Equation (25), as From Equations (24)-(26), one can find that the introduction of α can ensure that all blocks of J are O(∆t 0 ), avoiding the Jacobian matrix J becoming ill-conditioned from ∆t being too small. Newton-Raphson iteration stops if G 1 | t < ε (the iteration threshold error) is satisfied; then, we have ..

Numerical Properties
The spectral characteristic, overshoot characteristic, and convergence rate of the TTBIF for undamped and damped systems are analyzed in this section. Hereinafter, the TTBIFs using the optimal γ 1 provided in Tables 1 and 2 are called TTBIFa and TTBIFbn, respectively, in which n represents accuracy order. Consider the following single degree-of-freedom test equation [8] .. (27) where ξ and ω stand for the damping ratio and natural frequency, respectively, and f (t) is the external force. The above test equation (Equation (27)) is used to analyze the numerical properties of the TTBIF below.

Spectral Characteristics
The numerical stability, dissipation capability, and low-frequency accuracy of a time integration method can be evaluated by spectral characteristics [8]. Applying Equations (5)- (7) to the test equation (Equation (27)) can yield the following recursive equations: where A represents the amplification matrix, and its elements refer to Appendix C. Then, the characteristic polynomial of A can be derived from Equation (28), as where A 1 , A 2 , and A 3 represent the trace, the sum of 2 × 2 principle minors, and the determinant of A, respectively. For a convergent time integration method, three non-zero characteristic roots of A satisfy |λ 3 | ≤ |λ 1,2 | ≤ 1, where λ 1,2 are a pair of conjugate complex numbers, referred to as the principal roots, and λ 3 is the spurious root. The definition of spectral radius [8] is The numerical damping ratio [8] and the period elongation ratio [8] can be derived from Equation (29), which are employed to evaluate the amplitude and phase accuracy. Their definitions are numerical damping ratio : ξ = − ln(ρ) 2ω∆t (31) period elongation ratio : where ω∆t = arctan(b/a), and a and b are real and imaginary parts of the principal roots λ 1,2 , respectively. Figure 3 shows the variations of TTBIFa's spectral radius vs. Ω = ω∆t, and the results of the TTBIFbn are plotted in Figure 4. It can be seen that: (1) the TTBIFa and the TTBIFbn are unconditionally stable for undamped and damped systems; (2) ξ only affects ρ in the middle-frequency range; (3) if ξ = 0, the range corresponding to ρ(Ω) → 1 of the TTBIFa is wider than that of the TTBIFbn; (4) the size of γ 1 has little effect on ρ of the TTBIFbn-refer to Figure 4a-c; (5) if 0 < ξ ≤ 1, the TTBIFa and the TTBIFbn nearly have the same ρ for 0 < Ω ≤ 1.  In addition, the spectral radius curves shown in Figures 3 and 4 can help the users to choose a reasonable time step size. To stably solve nonlinear systems, the dissipative TTBIF (ρ ∞ < 1) is recommended. For dynamic problems with strong nonlinearity, the smaller ρ ∞ is recommended, and the larger ρ ∞ can be used for dynamic problems with weak nonlinearity. Once the value of ρ ∞ is selected, the time step size can be determined by ρ(ω cr ∆t) = 0.95, improving stability and keeping low-frequency accuracy. Herein, the ω cr represents the maximum frequency that needs to be kept, which can be calculated by the mass matrix and Jacobin matrix at the initial moment.
In the following, the accuracy of the TTBIF is compared to those of other time integration methods, including the generalized-α method [12], the SS2 method [26], and the Bathe method [31]. The changes of ρ with δΩ = ω∆t/n are plotted in Figure 5, in which δΩ = ω∆t/n can ensure that the accuracy comparison is conducted under the same cost for all methods. It follows that for 0 ≤ ρ ∞ < 1, the range corresponding to ρ(Ω) → 1 of the TTBIFa is the widest, whereas that of the TTBIFb is the narrowest.
To compare the spectral characteristics of different methods in the presence of damping ratio, the numerical spectral radius and the exact spectral radius ρ exact = exp(−ξΩ) [46] are shown in Figure 6. It can be concluded that: (1) if ξ = 0, the generalized-α method may be unstable, as ρ(Ω) > 1; (2) if ρ ∞ = 1, the physical damping has no effect on the ρ of the SS2 method, and if 0 ≤ ρ ∞ < 1, the ρ of the SS2 method agrees well with the exact one in the range of 0 < Ω ≤ 0.3; (3) the ρ of the Bathe method matches well with the exact one in the range of 0 < Ω ≤ 1; (4) in the range of 0 < Ω ≤ 4, the ρ of the TTBIFa agrees well with the exact one, and the ρ of the TTBIFb is consistent with the exact one in the range of 0 < Ω ≤1.  It can be found that for both undamped and damped systems, the TTBIFa's spectral radius is more accurate than those of other methods. Then, the amplitude and phase accuracy are further compared, as shown in Figures 7 and 8, in which one can see that: (1) the TTBIFa exhibits desirable amplitude and phase accuracy for δΩ ∈ (0, 1]; (2) the TTBIFb2 has considerable amplitude and phase errors in the low-frequency range, but the third-order accurate TTBIFb3 has the highest phase accuracy in the range of δΩ ∈ (0, 0.1]. Hence, the use of TTBIFb2 is not recommended owing to its larger numerical errors, and the TTBIFb3 can be applied to damped systems considering that it can quickly damp out the unwanted information and it has higher phase accuracy. Similar phenomena can be observed for other TTBIFb2 and TTBIFb3 under different ρ ∞ ; thus, they are not plotted herein for brevity.

Overshoot Characteristics
The overshooting phenomenon may occur in the first several time steps. For a convergent method, there is no overshoot as Ω → 0, so only the case of Ω → ∞ needs to be considered. The analysis of overshooting should take into account the effect of physical damping. With physical damping, first-order overshooting components enter into several well-known time integration methods [47] which were previously thought to exhibit zeroorder overshooting. Therefore, this section discusses the overshoot characteristics of the TTBIF for undamped (ξ = 0) and damped (ξ = 1) cases. If Ω → ∞, the recursive schemes of the TTBIF at the first step (t = 0) can be derived from Equation (28); i.e., The above equations disclose that when Ω → ∞ with a large ∆t, displacement, velocity, and acceleration of the TTBIF in the first step decay to constants for arbitrary conditions. It can be seen in Equations (33)-(35) that the overshoot characteristics of the TTBIF are independent of the value of γ 1 , so this section discusses the TTBIFa only; refer to  The test equation (Equation (27)), two initial conditions are considered, including two situations, q(0) = 0, . q(0) = 1 and q(0) = 1, . q(0) = 0. One can see in Figures 9-12 that the TTBIF does not have the overshoots for these two initial conditions, and the physical damping has no influence on the overshoot characteristics of the TTBIF.

Convergence Rates
The decrease rate of the absolute error as the step size decreases at a certain moment is defined as the convergence rate, which is usually employed to evaluate the accuracy order of a time integration method by displacement, velocity, and acceleration. In the test equation (Equation (27)), we assume that ξ = 2 √ 5, ω = √ 5, and f (t) = sin2t; and the exact solution is x(t) = exp(−2t)(cost + 2sint) − (8cos2t − sin2t)/65. The absolute errors of the TTBIF at time t = 1 are shown in Figure 13, in which one can see that: (1) the TTBIFa method is strictly second-order accurate for displacement, velocity, and acceleration, and its computational accuracy can be improved with an increase in ρ ∞ ; (2) the TTBIFb3 method can be third-order accurate for displacement, velocity, and acceleration; (3) the generalized-α method (G-α) is first-order accurate for acceleration when ρ ∞ < 1. The theoretical analysis given in this section illustrates that: (1) the low-frequency responses of the TTBIFa are more accurate, and it has the ability to damp out high-frequency modes; therefore, the TTBIFa are applicable to conservative systems and stiff systems; (2) the low-frequency response accuracy of the TTBIFb2 is low, so the TTBIFb2 is not practical; (3) the TTBIFb3 has strong numerical dissipation and lower numerical dispersion; hence, the TTBIFb3 is more suitable for damped systems.
Here, we provide a summary of the above time integration methods, and their accuracy, type of numerical dissipation, low-frequency accuracy, and scope of application are provided in Table 3, in which one can find that: (1) the generalized-α method [12], the Bathe method [31], and the TTBIF have been applied to multibody dynamic systems governed by DAEs; (2) among these methods, only TTBIF can achieve third-order accuracy for displacement, velocity, and acceleration, and it has numerical dissipation; (3) the ρ ∞ -Bathe method [37] and the Kim method [38], developed for structural dynamic systems governed by ODEs, have the same accuracy in the low-frequency range as the Bathe method [31] for the dissipative case ρ ∞ = 0; (4) from the values of ξ and PE, one can see that among the second-order accurate methods, the TTBIFa's low-frequency accuracy, including amplitude and phase, is the highest, meaning that it can give more accurate predictions when applied to dynamic problems, including structural dynamics and the multibody systems. Table 3. Numerical characteristics of some advanced time integration methods.

Numerical Experiments
To validate the desirable performances of the TTBIF, several benchmark multibody dynamic experiments are presented in this section. The generalized-α method [12], the SS2 method [26], and the Bathe method [31], are also considered. To ensure that the accuracy comparison was conducted with the same computational cost, the time step size relations of these methods were ∆t(Generalized-α method) = ∆t(SS2) = ∆t(Bathe)/2 = ∆t(TTBIF)/3. Therefore, in all examples, only the size of the generalized-α method is provided, and the reference solutions were obtained by the generalized-α method using a small time step size.

Flexible Beam
To show the advantages of the TTBIFa in stability, energy-conservation, and accuracy, this example considers a flexible beam model [48], as shown in Figure 14. The model parameters are: the length l = 3 m, the radius r = 0.02 m, the modulus of elasticity E = 2 × 10 6 Pa, and the density ρ = 7200 kg/m 3 . In the original configuration, the beam is horizontal and has zero initial velocity. The beam is discretized by 20 elements. For damping out numerical oscillations induced by the spatial discretization, the ρ ∞ = 0 is used in all methods. The cantilever beam under the effect of gravity is considered first, in which the time step size of the generalized-α method is 0.001 s. In this case, the beam is assumed to fall under the effect of gravity, as g = 9.81 m/s 2 . The configurations of the cantilever beam at different moments are plotted in Figure 15, and the results are provided by the TTBIFa. One can see that this model represents an extreme case of a large deformation problem that involves high inertia forces. The free falling beam is a conservative system, and its mechanical energy (or total energy), consisting of the kinetic energy, the strain energy, and the potential energy, should be zero. Figure 16 shows the variations of the energy with time. It can be seen that (1) the generalized-α method loses stability, and the other methods can give convergent results; (2) among the convergent methods, with the increase in time, the total energies of the SS2 method and the Bathe method show considerable attenuation, whereas the TTBIFa preserves energy well.  Then, the simply supported beam is considered to compare the accuracies of these methods, and the size of the generalized-α method is also 0.001 s. In the second case, the flexible beam is accelerated by increasing the gravity constant to g = 30 m/s 2 . The numerical results and the absolute errors at mid-point of the beam are plotted in Figures 17-19. It can be clearly seen that the TTBIFa's accuracy is the highest. Last, the accuracies of these methods in solving the problem of a flexible beam under the effect of vertical concentrated force f = 100 N, as shown in Figure 14c, are investigated, in which the size of the generalizedα method is assumed to be 0.001 s. The displacements of these methods for the case (x = l) and the case (x = l/2) are shown in Figures 20 and 21, respectively. It can be seen that the TTBIFa has a considerable advantage in the aspect of accuracy.

Slider-Pendulum
As shown in Figure 22, the second example considers a slider-pendulum model [49] to demonstrate the accuracy and dissipation of the TTBIFa in dealing with rigid multibody systems. The slider is constrained by the spring, and one end of the pendulum is hinged to the center of mass of the slider. The system parameters are m 1 = m 2 = 1 kg, L = 1 m, J 2 = 1/12 kg·m 2 , g = 9.81 m/s 2 , k = 1 N/m, and 10 16 N/m for the compliant and stiff systems. ρ ∞ = 0 is used in all methods, and the time step size of the generalized-α method is 0.06 s. Figure 22. Slider-pendulum model.
The compliant case (k = 1 N/m) is considered first, in which the slider is motivated by the horizontal velocity v 0 = 1 m/s. Figures 23-25 show that the TTBIFa performs best in the simulations. The results of the two single-step methods, the generalized-α method and the SS2 method, become more inaccurate with time. Additionally, the performances of these methods in dealing with numerical drifts are also considered here. The position constraint equation of this model has the form The variations of Φ 1 and Φ 2 vs. time are plotted in Figure 26. It can be observed that the TTBIFa behaves best in eliminating the drift phenomenon.
Then, the stiff case (k = 10 16 N/m) is considered, in which the pendulum is motivated by the horizontal velocity v 0 = 1 m/s. The time histories of x 1 and θ are shown in Figure 27, from which one can see that all methods can effectively filter out the high-frequency information induced by the stiff constraint. Additionally, among these methods, the calculation accuracy of the TTBIFa is the highest.

Moving Cable
The last example considers a moving cable model, as shown in Figure 28, to test the performances of the TTBIFa and the TTBIFb3 in solving damped problems. The cable with spatial moving and axial material flow is part of the arresting-cable system [50]. The absolute nodal coordinate formulation (ANCF) in the framework of the arbitrary Lagrange-Euler (ALE) description is adopted to deal with the geometric nonlinearity and variable length of the moving cable elements. The model has two flexible ALE cable elements, and 12 degrees of freedom. Partial spatial movements or material flows of these nodes are given and realized by the kinematic time-varying constraint Φ, as follows: The ρ ∞ = 0 is used in the generalized-α method, the SS2 method, and the TTBIFa to quickly filter out the numerical oscillations introduced by time-varying Φ. In addition, the third-order accurate TTBIFb3 with ρ ∞ = 0.7 and γ 1 = 1.64 is also considered here. In this example, the time step size of the generalized-α method is set to 0.01 s. The kinematic characteristics of mid-node 2 in the horizontal and vertical directions are shown in Figures 29-31, in which one can see that the TTBIFb3 is the most accurate, followed by the TTBIFa. In addition, it can be found from Figure 31 that the acceleration calculated by the generalized-α method has larger errors compared with other methods.

Conclusions
The authors constructed a three-sub-step composite time integration method with controllable numerical dissipation in 2020, which was named as TTBIF [39]. This work presented a deep understanding of the TTBIF, and the algorithmic parameters of the TTBIF was further optimized, so two new optimized schemes of the TTBIF were developed and applied to the simulations of multibody systems.
The main contributions of this study can be summarized as follows: (1) The algorithmic parameter of the TTBIF (γ 1 ) was optimized in this work to minimize local truncation error, yielding two new optimized schemes, the TTBIFa and the TTB-IFbn (n = 2,3). The second-order accurate TTBIFa can accurately keep low-frequency information, and all sub-steps share the same effective stiffness matrix. The TTB-IFb3 is third-order accurate, and possesses lower period errors than the second-order methods. Since the TTBIFb3 also has high dissipation capability, the TTBIFb3 is more suitable for damped systems. (2) The effects of algorithmic parameters on stability, dissipation, dispersion, overshoot, and convergence rate for damped and undamped systems were obtained in this work. The investigation is an important theoretical supplement for the TTBIF [39]. (3) The calculation formulation of the TTBIF was extended for the dynamic analysis of multibody systems, which can be borrowed for other implicit multi-sub-step methods, such as the TTBDF [32], the Wen method [33], and the TSSBN method [51]. (4) The numerical experiments demonstrated that the second-order accurate, controllably dissipative, unconditionally stable TTBIFa exhibits advantages in accuracy, dissipation, stability, and energy-conservation, so it seems to be a good candidate for the response analysis of multibody systems.

Conflicts of Interest:
The authors declare that they have no conflict of interest.
where A represents the TTBIF's amplification factor.