# 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.3 Equation Linearization and Assembly  
The nonlinear algebraic Equation (6) can then be linearized and solved iteratively, for instance, using the Newton iteration, as shown in our previous work [15]. To start the Newton iteration, an initial guess close enough is needed for \(\boldsymbol{\xi}^{(n)}\). With the first-order constant step BDFs, the initial guess of the differential variable \(\boldsymbol{\xi}^{(n)}\) is simply given by the same order extrapolation formula as  
\[\boldsymbol{\xi}^{(n,0)}=2\boldsymbol{\xi}^{(n-1)}-\boldsymbol{\xi}^{(n-2)}, \tag{7}\]  
wherein the two numbers in the right superscript, enclosed by parenthesis \(\square^{(n,i)}\), represent the results of the \(j^{\text{th}}\) iteration for the \(n^{\text{th}}\) time step. Consequently, \(\boldsymbol{\xi}^{(n,0)}\) represents the initial value in the iteration for the \(n^{\text{th}}\) time step. The Lagrange multipliers are not differentiable with respect to time, and the initial guess is set to zero as  
\[\boldsymbol{\lambda}^{(n,0)}=\boldsymbol{0}. \tag{8}\]  
We rewrite the nonlinear Equation (6) on one side to obtain their residual form:  
\[\left\{ \begin{array}{l}\boldsymbol{b}(\boldsymbol{\xi}^{(n)})\triangleq \boldsymbol{F}_{\text{iner}}\Big{(}\boldsymbol{\xi}^{(n)},f_{\text{v}}( \boldsymbol{\xi}^{(n)}),f_{\text{a}}(\boldsymbol{\xi}^{(n)})\Big{)}+F_{\text{ elas}}\Big{(}\boldsymbol{\xi}^{(n)}\Big{)}-\boldsymbol{\tau}(\boldsymbol{\xi}_{nh},f_{ \text{v}}(\boldsymbol{\xi}_{nh}))+c\boldsymbol{\xi}^{\text{I}}_{\boldsymbol{ \xi}}(\boldsymbol{\xi}^{(n)})\boldsymbol{\lambda}^{(n)}=\boldsymbol{0}\,\\ \boldsymbol{c}(\boldsymbol{\xi}^{(n)})=\boldsymbol{0}\,\end{array} \right. \tag{9}\]and the Newton iteration is given by  
\[ \begin{bmatrix}\frac{\partial\mathbf{b}}{\partial\xi^{(n)}}&\left(\frac{\partial\mathbf{c} }{\partial\xi^{(n)}}\right)^{\mathrm{T}}\end{bmatrix}\begin{bmatrix}\Delta\xi^{ (n,j)}\\ \Delta\lambda^{(n,j)}\end{bmatrix} =- \begin{bmatrix}\mathbf{b}(\xi^{(n,j)})\\ \mathbf{c}(\xi^{(n,j)})\end{bmatrix} , \tag{10}\]  
wherein  
\[\frac{\partial\mathbf{b}}{\partial\xi^{(n)}}=\frac{\partial\mathbf{F}_{\mathrm{iner}} }{\partial\xi}+\frac{\partial\mathbf{F}_{\mathrm{iner}}}{\partial\dot{\xi}}\frac{ \partial f_{\mathrm{v}}}{\partial\xi^{(n)}}+\frac{\partial\mathbf{F}_{\mathrm{iner }}}{\partial\dot{\xi}}\frac{\partial f_{\mathrm{a}}}{\partial\xi^{(n)}}+\frac {\partial\mathbf{F}_{\mathrm{elas}}}{\partial\xi}-\frac{\partial\mathbf{\tau}}{ \partial\xi}\frac{\partial f_{\mathrm{v}}}{\partial\xi^{(n)}}+\frac{\partial \Big{(}\bar{c}_{\xi}^{\intercal}\lambda^{(n)}\Big{)}}{\partial\xi}. \tag{11}\]  
Typically, the second-order differential term \(\frac{\partial\big{(}\bar{c}_{\xi}^{\intercal}\lambda^{(n)}\big{)}}{\partial\xi}\) is omitted in Formula (11) to ensure the stability of the Newton iteration because of the quick variation in the variable \(\lambda^{(n)}\). Utilizing the time discrete Formula (4), we have  
\[\frac{\partial f_{\mathrm{v}}}{\partial\xi^{(n)}}=\frac{1}{h},\ \frac{\partial f_{\mathrm{a}}}{\partial\xi^{(n)}}=\frac{1}{h^{2}}, \tag{12}\]  
and denote the linearized general mass, damping, and stiffness matrices of the bodies as  
\[\mathbf{M}\triangleq\frac{\partial\mathbf{F}_{\mathrm{iner}}}{\partial\dot{\xi}},\mathbf{ C}\triangleq\frac{\partial\mathbf{F}_{\mathrm{iner}}}{\partial\dot{\xi}},\mathbf{K} \triangleq\frac{\partial\mathbf{F}_{\mathrm{elas}}}{\partial\dot{\xi}}+\frac{ \partial\mathbf{F}_{\mathrm{iner}}}{\partial\dot{\xi}}, \tag{13}\]  
The Newton iteration (10) for the nonlinear Equation (9) is rewritten as  
\[\left[ \begin{array}{cc}\frac{M}{h^{2}}+\frac{C}{h}+\mathbf{K}-\left(\frac{ \partial\mathbf{\tau}}{\partial\xi}+\frac{1}{h}\cdot\frac{\partial\mathbf{\tau}}{ \partial\dot{\xi}}\right)&\frac{\partial\mathbf{c}^{\intercal}}{\partial\xi}\\ \frac{\partial c}{\partial\xi}&\mathbf{0}\end{array} \right]^{(n,j)}\left[\Delta \xi^{(n,j)}_{\mathbf{\xi}}\right]=- \begin{bmatrix}\mathbf{b}\\ \mathbf{c}\end{bmatrix}