# Parallel Direct Solution of Flexible Multibody Systems Based on Block Gaussian Elimination
## 2 Flexible Multibody Formulation with the Craig-Bampton Method
### Time Integration Procedures of Flexible Multibody Systems
#### 2.1.2 Time Integration Scheme  
To solve the algebraic-differential Equation (3), we need to apply a time integration to convert the general velocity \(\boldsymbol{\xi}\) and acceleration \(\ddot{\boldsymbol{\xi}}\) into displacement \(\boldsymbol{\xi}\) at the time integrationsteps. Since multibody problems are typically stiff, meaning that the system responses contain both low and high frequencies simultaneously, the stiff integrator--backward difference formulas (BDFs) is utilized in this paper.  
For simplicity, we use a constant step size \(h\) and the first-order BDFs to demonstrate the integration process. To avoid confusion, we also employ a right superscript within parentheses \(\square^{(n)}\) to denote variable values at the discrete time \(t=nh\), such as \(\boldsymbol{\xi}^{(n)}\triangleq\boldsymbol{\xi}\,|_{t=nh}\) and \(\boldsymbol{\lambda}^{(n)}\triangleq\boldsymbol{\lambda}\,|_{t=nh}\). Assuming the previous \(n-1\) time steps are known and the \(n^{\text{th}}\) time step needs to be computed, the general velocity and acceleration can be expressed in terms of displacement as  
\[\left\{\begin{array}{l}\dot{\boldsymbol{\xi}}^{(n)}=\frac{\xi^{(n)}-\xi^{(n -1)}}{h}\triangleq f_{\text{v}}(\boldsymbol{\xi}^{(n)}),\\ \bar{\boldsymbol{\xi}}^{(n)}=\frac{\dot{\xi}^{(n)}-\dot{\xi}^{(n-1)}}{h}=\frac {\xi^{(n)}-2\xi^{(n-1)}+\dot{\xi}^{(n-2)}}{h^{2}}\triangleq f_{\text{a}}( \boldsymbol{\xi}^{(n)}),\end{array}\right. \tag{4}\]  
wherein \(f_{\text{v}}\) represents the velocity's algebraic function and \(f_{\text{a}}\) represents the acceleration's algebraic function.  
Readers may notice that Formula (4) is not applicable when \(n=1\), where \(\xi^{(-1)}\) is requested in forming \(\bar{\boldsymbol{\xi}}^{(1)}\). A possible solution is to estimate \(\bar{\boldsymbol{\xi}}^{(-1)}\) using the Taylor series as  
\[\bar{\boldsymbol{\xi}}^{(-1)}=\bar{\boldsymbol{\xi}}^{(0)}-\dot{\bar{ \boldsymbol{\xi}}}^{(0)}h+\bar{\bar{\boldsymbol{\xi}}}^{(0)}\frac{h^{2}}{2}+o( h^{2}), \tag{5}\]  
wherein \(\bar{\boldsymbol{\xi}}^{(0)}\) and \(\dot{\bar{\boldsymbol{\xi}}}^{(0)}\) are given by the initial conditions and \(\bar{\bar{\boldsymbol{\xi}}}^{(0)}\) can be calculated by solving the system Equation (3) at time \(t=0\).  
By substituting Equation (4) into Equation (3), nonlinear algebraic equations at the time step \(t=nh\) with unknown variables \(\boldsymbol{\xi}^{(n)}\) and \(\boldsymbol{\lambda}^{(n)}\) are obtained as  
\[\left\{\begin{array}{l}\boldsymbol{F}_{\text{iner}}\Big{(}\boldsymbol{\xi} ^{(n)},f_{\text{v}}(\boldsymbol{\xi}^{(n)}),f_{\text{a}}(\boldsymbol{\xi}^{( n)})\Big{)}+\boldsymbol{F}_{\text{elas}}\Big{(}\boldsymbol{\xi}^{(n)}\Big{)}= \boldsymbol{\tau}\Big{(}\boldsymbol{\xi}^{(n)},f_{\text{v}}(\boldsymbol{ \xi}_{nh})\Big{)}-c\bar{\boldsymbol{\xi}}^{\text{I}}(\boldsymbol{\xi}^{(n)}) \boldsymbol{\lambda}^{(n)},\\ \boldsymbol{c}(\boldsymbol{\xi}^{(n)})=\boldsymbol{0},\end{array}\right. \tag{6}\]  
wherein the nonlinearity may arise from the expression of \(\boldsymbol{F}_{\text{iner}}\), \(\boldsymbol{F}_{\text{elas}}\), \(\boldsymbol{\tau}\), and \(\boldsymbol{c}\).