## 2 Equations of Motion for Flexible Multibody Systems
### Coordinate Reduction of Flexible Bodies  
In this section, we review the FMBD formulation including the coordinate reduction of flexible bodies [2; 27; 28]. The popular Craig-Bampton (CB) method [7] is here considered to reduce the degrees of freedom (DOFs) of a flexible body. The equations of motion of an undamped flexible body under free vibration can be expressed as:  
\[\mathbf{M}_{F}^{i}\ddot{\mathbf{q}}_{F}^{i}+\mathbf{K}_{F}^{i}\mathbf{q}_{F}^ {i}=\mathbf{0}, \tag{9}\]and its eigenvalue problem is:  
\[\mathbf{K}_{F}^{i}\boldsymbol{\varphi}_{n}^{i}=\chi_{n}^{i}\mathbf{M}_{F}^{i} \boldsymbol{\varphi}_{n}^{i}\quad n=1,2,\cdots,N_{F}^{i}, \tag{10}\]  
where \(\chi_{n}^{i}\) and \(\boldsymbol{\varphi}_{n}^{i}\) are the \(n\)-th eigenvalue and its corresponding eigenvector, respectively. Here, \(N_{F}^{i}\) is the number of DOFs of the flexible body. The mass and stiffness matrices, \(\mathbf{M}_{F}^{i}\) and \(\mathbf{K}_{F}^{i}\), can be partitioned to the interior DOFs and the boundary DOFs denoted using subscripts \(s\) and \(b\). We consider the connecting points of joints, external forces, and the boundary conditions as boundary DOFs for the flexible body. Then, we have:  
\[\mathbf{M}_{F}^{i}= \begin{bmatrix}\mathbf{M}_{s}^{i}&\mathbf{M}_{c}^{i}\\ \mathbf{M}_{c}^{i\,T}&\mathbf{M}_{b}^{i}\end{bmatrix} ,\;\mathbf{K}_{F}^{i}= \begin{bmatrix}\mathbf{K}_{s}^{i}&\mathbf{K}_{c}^{i}\\ \mathbf{K}_{c}^{i\,T}&\mathbf{K}_{b}^{i}\end{bmatrix} ,\;\mathbf{q}_{F}^{i}= \begin{bmatrix}\mathbf{q}_{s}^{i}\\ \mathbf{q}_{b}^{i}\end{bmatrix} . \tag{11}\]  
The numbers of interior and boundary DOFs are denoted as \(N_{s}^{i}\) and \(N_{b}^{i}\), respectively. In the CB method, the fixed-interface normal mode \(\boldsymbol{\Phi}_{s}^{i}\) and the interface constraint mode \(\boldsymbol{\
abla}^{i}\) are used to reduce and assemble the flexible body coordinates. The fixed-interface normal mode can be obtained through the eigenvalue problem of the mass and stiffness matrices with respect to the interior part as:  
\[\left[\mathbf{K}_{s}^{i}-\kappa_{n}^{i}\mathbf{M}_{s}^{i}\right]\boldsymbol{ \phi}_{n}^{i}=0,\;\;n=1,2,\cdots,N_{s}^{i}, \tag{12a}\] \[\boldsymbol{\Phi}_{s}^{i}=\left[\boldsymbol{\phi}_{1}^{i}\;\boldsymbol{\phi}_ {2}^{i}\;\cdots\;\boldsymbol{\phi}_{N_{s}^{i}}^{i}\right], \tag{12b}\]  
where \(\kappa_{n}^{i}\) and \(\boldsymbol{\phi}_{n}^{i}\) denote the eigenvalue and eigenvector of the interior part, respectively. The interface constraint mode is then defined as [7]:  
\[\boldsymbol{\
abla}^{i}=-\left[\mathbf{K}_{s}^{i}\right]^{-1}\mathbf{K}_{c}^ {i}. \tag{13}\]  
The transformation matrix \(\mathbf{T}_{0}^{i}\) can be constructed using Equations (12a) and (13), then the generalized coordinate \(\mathbf{q}_{F}^{i}\) can be written with the hybrid coordinate \(\mathbf{p}_{F}^{i}\) as:  
\[\mathbf{q}_{F}^{i}=\mathbf{T}_{0}^{i}\mathbf{p}_{F}^{i},\;\;\mathbf{T}_{0}^{i }= \begin{bmatrix}\boldsymbol{\Phi}_{s}^{i}&\boldsymbol{\
abla}^{i}\\ \boldsymbol{0}&\mathbf{I}\end{bmatrix} ,\;\;\mathbf{p}_{F}^{i}= \begin{bmatrix} \mathbf{p}_{s}^{i}\\ \mathbf{q}_{b}^{i}\end{bmatrix} . \tag{14}\]  
The generalized coordinate vector of the interior \(\mathbf{q}_{s}^{i}\) is expressed as:  
\[\mathbf{q}_{s}^{i}=\boldsymbol{\Phi}_{s}^{i}\mathbf{p}_{s}^{i}+\boldsymbol{ \
abla}^{i}\mathbf{q}_{b}^{i},\;\boldsymbol{\Phi}_{s}^{i}=\left[\boldsymbol{ \Phi}_{d}^{i}\quad\boldsymbol{\Phi}_{r}^{i}\right],\;\mathbf{p}_{s}^{i}= \begin{bmatrix}\mathbf{p}_{d}^{i}\\ \mathbf{p}_{r}^{i}\end{bmatrix} , \tag{15}\]  
in which the subscripts \(d\) and \(r\) denote the dominant and residual terms, respectively. The number of the dominant modes \(\boldsymbol{\Phi}_{d}^{i}\) and the residual modes \(\boldsymbol{\Phi}_{r}^{i}\) is denoted as \(N_{d}^{i}\) and \(N_{r}^{i}\), respectively (\(N_{s}^{i}=N_{d}^{i}+N_{r}^{i}\), \(N_{d}^{i}\ll N_{r}^{i}\)).