Iterative Coordinate Reduction Algorithm of Flexible Multibody Dynamics Using a Posteriori Eigenvalue Error Estimation

Abstract

Coordinate reduction has been widely used for efficient simulation of flexible multibody dynamics. To achieve the reduction of flexible bodies with reasonable accuracy, the appropriate number of dominant modes used for the reduction process must be selected. To handle this issue, an iterative coordinate reduction strategy is introduced. In the iteration step, more dominant modes of flexible bodies are selected than the ones in the previous step. Among the various methods, the conventional frequency cut-off rule is here considered. As a stop criterion, a novel a posteriori error estimator that can evaluate the relative eigenvalue error between full and reduced flexible bodies is proposed. Through the estimated relative eigenvalue error obtained, the number of dominant modes is automatically selected to satisfy the error tolerance up to the desired mode range. The applicability to the automation process is verified through numerical examples. It is also evaluated that efficient and accurate flexible multibody dynamics simulation is available with the reduced flexible body, generated by the proposed algorithm.

1. Introduction

Flexible multibody dynamics (FMBD) simulation is a popular way to consider system dynamics with elastic deformation. The FMBD approach might be powerful because it can provide stress and strain information unlike with rigid body dynamics, but it requires much heavier computational cost than multibody dynamics (MBD) does. The computational cost comes from flexible bodies that are finite element (FE) models with large degrees of freedom (DOFs). Coordinate reduction techniques of the flexible bodies are then essential to effective FMBD simulation [1,2,3,4,5,6]. Among the various methods, component mode synthesis (CMS) [7,8,9,10,11,12], a projection technique with an eigenbasis, may be the most popular method in linear vibration with FE analysis. In the CMS methods, the residual eigenvectors that are less important for representing the dynamic characteristics of the original FE model are truncated. The constraint modes are then used to define the interface boundary condition of the reduced model [7,10,12], or the attachment modes may be employed instead of the constraint modes [9,13]. Those depend on the type of CMS method. The small number of selected eigenvectors, known as dominant eigenvectors, are only used in the projection of the original FE model, which is a key to the model reduction process. The following two questions then arise [2]: (a) Which eigenvectors are dominant? (b) How many eigenvectors should be used to reduce FE models accurately? With this motivation, many researchers have studied mode selection and error estimation methods for use with the CMS methods.

The rule of thumb of mode selection may be the frequency cut-off method in the mode superposition manner, but in the CMS methods, the effective interface mass (EIM) methods by Kammer and Tiller [14,15] may be the first trial for the mode selection method. The EIM methods were developed for the Craig–Bampton (CB) method, which is the most popular CMS method. The optimal model reduction (OMR) methods based on the Dirichlet-to-Neumann (DTN) map have been studied under the substructuring strategy [16,17]. Those may be the first mode selections considering the coupling effect between substructures. The CMS${}_{\chi }$ method is an extension of the OMR method using the moment-matching approach instead of the DTN map [18]. Kim et al. inherited the CMS${}_{\chi }$ method without the relatively rough assumption in the derivation process, which is called the CMS${}_{\sigma }$ method [19]. Unlike the above approaches, the eigenvector-based modal contribution (EMC) technique is an a posteriori data-driven method that computes the substructural modal contribution by considering the relation between the original and reduced matrices [20]. Therefore, the EMC method can provide better mode selection performance than a priori methods can, but requires relatively huge computational cost. Those issues were well examined by Kim et al. [20].

Ranking of the dominant modes can be arranged using the mode selection method. We can then make a selection of the dominant modes, and then, the modes could be added to the reduced matrices in the iterative way for better accuracy. However, more modes lead to larger reduced system matrices, requiring a compromise between computation cost and accuracy. The use of mode selection only may not be enough to resolve this issue. It can be treated using the error evaluation methods of the reduced matrices. Among the various methods, the error estimators of the CB method and the Guyan reduction proposed by Kim et al. may be the most highly accurate [21,22,23,24]. Those approaches are direct approximations of the relative eigenvalue errors of the CB and Guyan reductions. The direct eigenvalue error estimation is then available without knowing the reference eigenvalue.

Combining the error estimators and the above mode selection methods can provide an automated model reduction algorithm excluding empirical judgment [25]. This work was accelerated toward that idea. An iterative coordinate reduction algorithm of the FMBD analysis is then proposed. The floating frame of reference formulation (FFRF) [1,26] and its CB-based coordinate reduction are here considered so that those may be the standard approaches in various commercial FMBD software. In the proposed iterative algorithm, the conventional frequency cut-off rule is considered to select more simply the important flexible modes, and the reduced matrices could be initially constructed. The accuracy of the eigenvalues approximated by obtaining the reduced eigenvalue problem can be evaluated using the new CB error estimation method, which is a more generalized formulation than in the previous work [22]. If the accuracy level is less than the desired tolerance, the flexible modes are iteratively added to the reduced matrices constructed in the previous step. In this manner, reduced matrices with the required accuracy and a reasonable dimensions could be obtained.

In the following section, the FFRF approach with the CB coordinate reduction is reviewed, and the CB error estimator is presented in Section 3. The iterative coordinate reduction algorithm is described in Section 4. The performances and characteristics of the proposed algorithm are investigated through numerical examples in Section 5. The conclusions are given in Section 6.

2. Equations of Motion for Flexible Multibody Systems

The floating frame of reference formulation (FFRF) [1,26] is widely used to describe the behavior of multibody dynamics systems with flexible bodies. Two types of coordinate systems are combined in FFRF; one type is the reference coordinates used to describe the translational and rotational motion of a body. The other type involves the flexible coordinates used for the deformation of a body. To derive equations of motion (EOMs), the i-th body in an FMBD system is considered.

2.1. Floating Frame of Reference Formulation

The Lagrange equation for the i-th body in the FMBD system can be written as:

$$\begin{array}{c}\hfill \frac{d}{dt}{\left(\frac{\partial T^i}{\partial {\dot{\mathbf{q}}}^i}\right)}^T-{\left(\frac{\partial T^i}{\partial {\mathbf{q}}^i}\right)}^T+{\mathbf{C}}_{{\mathbf{q}}^i}^T{\lambda }^i={\mathbf{Q}}^i,\end{array}$$

(1)

where ${\mathbf{q}}^i$ is the generalized coordinate vector including the reference (rigid) and flexible coordinates. $T^i$ is the kinetic energy. ${\mathbf{C}}_{{\mathbf{q}}^i}^T$ is the constraint Jacobian matrix. ${\lambda }^i$ is the vector of Lagrange multipliers. ${\mathbf{Q}}^i$ is the vector of generalized forces. The generalized forces ${\mathbf{Q}}^i$ acting on the i-th body can be derived using the virtual work $\delta W_d^i$ of the elastic force and $\delta W_e^i$ of the external forces as:

$$\begin{array}{c}\hfill \delta W^i=\delta W_d^i+\delta W_e^i,\end{array}$$

(2)

where $\delta W_d^i$ and $\delta W_e^i$ are defined as:

$$\begin{array}{c}\hfill \delta W_d^i=-{{\mathbf{q}}^i}^T{\mathbf{K}}^i\delta {\mathbf{q}}^i,\delta W_e^i={{\mathbf{Q}}_e^i}^T\delta {\mathbf{q}}^i,\end{array}$$

(3)

where ${\mathbf{K}}^i$ is the stiffness matrix and ${\mathbf{Q}}_e^i$ is the external force vector. Substituting Equation (3) into Equation (2), the generalized force ${\mathbf{Q}}^i$ can be expressed as:

$$\begin{array}{c}\hfill {\mathbf{Q}}^i=-{\mathbf{K}}^i{\mathbf{q}}^i+{\mathbf{Q}}_e^i.\end{array}$$

(4)

The terms on the left side in Equation (1) without the third term can be written as [1]:

$$\begin{array}{c}\hfill \frac{d}{dt}{\left(\frac{\partial T^i}{\partial {\dot{\mathbf{q}}}^i}\right)}^T-{\left(\frac{\partial T^i}{\partial {\mathbf{q}}^i}\right)}^T={\mathbf{M}}^i{\ddot{\mathbf{q}}}^i-{\mathbf{Q}}_v^i,\end{array}$$

(5)

where ${\mathbf{M}}^i$ is the mass matrix and ${\mathbf{Q}}_v^i$ is the quadratic velocity vector defined as:

$$\begin{array}{c}\hfill {\mathbf{Q}}_v^i=-{\dot{\mathbf{M}}}^i{\dot{\mathbf{q}}}^i+{\left[\frac{\partial }{\partial {\mathbf{q}}^i}\left(\frac{1}{2}{{\dot{\mathbf{q}}}^i}^T{\mathbf{M}}^i{\dot{\mathbf{q}}}^i\right)\right]}^T.\end{array}$$

(6)

Using Equations (4) and (5) in Equation (1), the following EOMs neglecting the damping are obtained:

$$\begin{array}{c}\hfill {\mathbf{M}}^i{\ddot{\mathbf{q}}}^i+{\mathbf{K}}^i{\mathbf{q}}^i+{\mathbf{C}}_{{\mathbf{q}}^i}^T{\lambda }^i={\mathbf{Q}}_e^i+{\mathbf{Q}}_v^i,\end{array}$$

(7)

and it can be specifically written in a partitioned form of the reference and flexible coordinates as:

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{\mathbf{M}}_R^i& {\mathbf{M}}_{FR}^i\\ {\mathbf{M}}_{RF}^i& {\mathbf{M}}_F^i\end{array}\right]\left[\begin{array}{c}{\ddot{\mathbf{q}}}_R^i\\ {\ddot{\mathbf{q}}}_F^i\end{array}\right]+\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{K}}_F^i\end{array}\right]\left[\begin{array}{c}{\mathbf{q}}_R^i\\ {\mathbf{q}}_F^i\end{array}\right]+\left[\begin{array}{c}{\mathbf{C}}_{{\mathbf{q}}_R^i}^T\\ {\mathbf{C}}_{{\mathbf{q}}_F^i}^T\end{array}\right]{\lambda }^i=\left[\begin{array}{c}{\left({\mathbf{Q}}_e^i\right)}_R\\ {\left({\mathbf{Q}}_e^i\right)}_F\end{array}\right]+\left[\begin{array}{c}{\left({\mathbf{Q}}_v^i\right)}_R\\ {\left({\mathbf{Q}}_v^i\right)}_F\end{array}\right],\end{array}$$

(8)

where subscripts R and F denote the quantities with respect to the reference and flexible coordinates, respectively.

2.2. 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},$$

(9)

and its eigenvalue problem is:

$${\mathbf{K}}_F^i{\mathit{\phi }}_n^i={\chi }_n^i{\mathbf{M}}_F^i{\mathit{\phi }}_n^{i}n=1,2,\cdots ,N_F^i,$$

(10)

where ${\chi }_n^i$ and ${\mathit{\phi }}_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=\left[\begin{array}{cc}{\mathbf{M}}_s^i& {\mathbf{M}}_c^i\\ {{\mathbf{M}}_c^i}^T& {\mathbf{M}}_b^i\end{array}\right],{\mathbf{K}}_F^i=\left[\begin{array}{cc}{\mathbf{K}}_s^i& {\mathbf{K}}_c^i\\ {{\mathbf{K}}_c^i}^T& {\mathbf{K}}_b^i\end{array}\right],{\mathbf{q}}_F^i=\left[\begin{array}{c}{\mathbf{q}}_s^i\\ {\mathbf{q}}_b^i\end{array}\right].$$

(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 ${\mathsf{\Phi }}_s^i$ and the interface constraint mode ${\mathsf{\Psi }}^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:

$$\begin{array}{c}\hfill \left[{\mathbf{K}}_s^i-{\kappa }_n^i{\mathbf{M}}_s^i\right]{\mathit{\varphi }}_n^i=0,n=1,2,\cdots ,N_s^i,\end{array}$$

(12a)

$$\begin{array}{c}\hfill {\mathsf{\Phi }}_s^i=\left[{\mathit{\varphi }}_1^i{\mathit{\varphi }}_2^i\cdots {\mathit{\varphi }}_{N_s^i}^i\right],\end{array}$$

(12b)

where ${\kappa }_n^i$ and ${\mathit{\varphi }}_n^i$ denote the eigenvalue and eigenvector of the interior part, respectively. The interface constraint mode is then defined as [7]:

$$\begin{array}{c}\hfill {\mathsf{\Psi }}^i=-{\left[{\mathbf{K}}_s^i\right]}^{-1}{\mathbf{K}}_c^i.\end{array}$$

(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:

$$\begin{array}{c}\hfill {\mathbf{q}}_F^i={\mathbf{T}}_0^i{\mathbf{p}}_F^i,{\mathbf{T}}_0^i=\left[\begin{array}{cc}{\mathsf{\Phi }}_s^i& {\mathsf{\Psi }}^i\\ \mathbf{0}& \mathbf{I}\end{array}\right],{\mathbf{p}}_F^i=\left[\begin{array}{c}{\mathbf{p}}_s^i\\ {\mathbf{q}}_b^i\end{array}\right].\end{array}$$

(14)

The generalized coordinate vector of the interior ${\mathbf{q}}_s^i$ is expressed as:

$$\begin{array}{c}\hfill {\mathbf{q}}_s^i={\mathsf{\Phi }}_s^i{\mathbf{p}}_s^i+{\mathsf{\Psi }}^i{\mathbf{q}}_b^i,{\mathsf{\Phi }}_s^i=\left[\begin{array}{cc}{\mathsf{\Phi }}_d^i& {\mathsf{\Phi }}_r^i\end{array}\right],{\mathbf{p}}_s^i=\left[\begin{array}{c}{\mathbf{p}}_d^i\\ {\mathbf{p}}_r^i\end{array}\right],\end{array}$$

(15)

in which the subscripts d and r denote the dominant and residual terms, respectively. The number of the dominant modes ${\mathsf{\Phi }}_d^i$ and the residual modes ${\mathsf{\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$).

Neglecting the residual modes, the original generalized coordinate vector ${\mathbf{q}}_F^i$ can be approximated using the reduced transformation matrix ${\overline{\mathbf{T}}}_0^i$ as:

$$\begin{array}{c}\hfill {\mathbf{q}}_F^i\approx {\overline{\mathbf{q}}}_F^i={\overline{\mathbf{T}}}_0^i{\overline{\mathbf{p}}}_F^i,\end{array}$$

(16a)

$$\begin{array}{c}\hfill {\overline{\mathbf{T}}}_0^i=\left[\begin{array}{cc}{\mathsf{\Phi }}_d^i& {\mathsf{\Psi }}^i\\ \mathbf{0}& \mathbf{I}\end{array}\right],{\overline{\mathbf{p}}}_F^i=\left[\begin{array}{c}{\mathbf{p}}_d^i\\ {\mathbf{q}}_b^i\end{array}\right].\end{array}$$

(16b)

Using Galerkin projection with the transformation matrix ${\overline{\mathbf{T}}}_0^i$, Equation (9) can be reduced as:

$$\begin{array}{c}\hfill {\overline{\mathbf{M}}}_F^i{\ddot{\overline{\mathbf{p}}}}_F^i+{\overline{\mathbf{K}}}_F^i{\overline{\mathbf{p}}}_F^i=\mathbf{0},\end{array}$$

(17a)

$$\begin{array}{c}\hfill {\overline{\mathbf{M}}}_F^i={{\overline{\mathbf{T}}}_0^i}^T{\mathbf{M}}_F^i{\overline{\mathbf{T}}}_0^i,{\overline{\mathbf{K}}}_F^i={{\overline{\mathbf{T}}}_0^i}^T{\mathbf{K}}_F^i{\overline{\mathbf{T}}}_0^i,\end{array}$$

(17b)

and the corresponding eigenvalue problem is:

$$\begin{array}{c}\hfill {\overline{\mathbf{K}}}_F^i{\mathsf{\Xi }}_n^i={\overline{\chi }}_n^i{\overline{\mathbf{M}}}_F^i{\mathsf{\Xi }}_n^{i}n=1,2,\cdots ,{\overline{N}}^i.\end{array}$$

(18)

The number of DOFs for the reduced flexible body ${\overline{N}}^i$ depends on the number of dominant modes $N_d$ and the number of boundary DOFs $N_b$$({\overline{N}}^i=N_d^i+N_b^i)$. The number of dominant modes is much smaller than the number of residual modes, and ${\overline{N}}^i$ can be dramatically reduced compared to the original DOFs $N_F^i$. From the eigenvector ${\mathsf{\Xi }}_n^i$, the original eigenvector ${\mathit{\phi }}_n^i$ in Equation (10) is approximated as:

$$\begin{array}{c}\hfill {\mathit{\phi }}_n^i\approx {\overline{\mathit{\phi }}}_n^i={\overline{\mathbf{T}}}_0^i{\mathsf{\Xi }}_n^i.\end{array}$$

(19)

2.3. Floating Frame of Reference Formulation with the Reduced Flexible Bodies

The generalized coordinate vector ${\mathbf{q}}^i$ can be approximated using the transformation matrix ${\overline{\mathbf{T}}}_0^i$ as:

$$\begin{array}{c}\hfill {\mathbf{q}}^i\approx {\overline{\mathbf{q}}}^i={\mathbf{B}}^i\left[\begin{array}{c}{\mathbf{q}}_R^i\\ {\overline{\mathbf{p}}}_F^i\end{array}\right],{\mathbf{B}}^i=\left[\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{0}& {\overline{\mathbf{T}}}_0^i\end{array}\right].\end{array}$$

(20)

The EOMs of a reduced FMBD system can be obtained by substituting Equation (20) into Equation (8) and premultiplying by ${{\mathbf{B}}^i}^T$ as follows:

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{\mathbf{M}}_R^i& {\overline{\mathbf{M}}}_{RF}^i\\ {\overline{\mathbf{M}}}_{FR}^i& {\overline{\mathbf{M}}}_F^i\end{array}\right]\left[\begin{array}{c}{\ddot{\mathbf{q}}}_R^i\\ {\ddot{\overline{\mathbf{p}}}}_F^i\end{array}\right]+\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\overline{\mathbf{K}}}_F^i\end{array}\right]\left[\begin{array}{c}{\mathbf{q}}_R^i\\ {\overline{\mathbf{p}}}_F^i\end{array}\right]+\left[\begin{array}{c}{\mathbf{C}}_{{\mathbf{q}}_R^i}^T\\ {\overline{\mathbf{C}}}_{{\overline{\mathbf{p}}}_F^i}^T\end{array}\right]{\lambda }^i=\left[\begin{array}{c}{\left({\mathbf{Q}}_e^i\right)}_R\\ {\left({\overline{\mathbf{Q}}}_e^i\right)}_F\end{array}\right]+\left[\begin{array}{c}{\left({\mathbf{Q}}_v^i\right)}_R\\ {\left({\overline{\mathbf{Q}}}_v^i\right)}_F\end{array}\right],\end{array}$$

(21)

where the submatrices and vectors in Equation (21) are expressed as:

$$\begin{array}{cc}\hfill & {\overline{\mathbf{M}}}_{RF}^i={{\overline{\mathbf{M}}}_{FR}^i}^T={\mathbf{M}}_{RF}^i{\overline{\mathbf{T}}}_0^i,{\overline{\mathbf{C}}}_{{\overline{\mathbf{p}}}_F^i}={\mathbf{C}}_{{\mathbf{q}}_F^i}{\overline{\mathbf{T}}}_0^i,\hfill \\ \hfill & {\left({\overline{\mathbf{Q}}}_e^i\right)}_F=\overline{\mathbf{T}}{{}_0^i}^T{\left({\mathbf{Q}}_e^i\right)}_F,{\left({\overline{\mathbf{Q}}}_v^i\right)}_F=\overline{\mathbf{T}}{{}_0^i}^T{\left({\mathbf{Q}}_v^i\right)}_F.\hfill \end{array}$$

(22)

It is clear that the accuracy of the reduced EOMs depends on approximation of the flexible bodies. Various approaches such as mode selection methods [14,15,16,17,19], modal compensation [4,5,10,11,12,29], etc., have been studied to achieve better approximation of the reduced flexible bodies in structural dynamics, but increasing $N_d^i$ may be the simplest and the most popular remedy to get more accurate ${\overline{\chi }}_n^i$.

3. Error Estimation for a Reduced Flexible Body

The following relative eigenvalue error has been widely used to evaluate the accuracy of a reduced flexible body:

$$\begin{array}{c}\hfill {\zeta }_n^i=\frac{{\overline{\chi }}_n^i-{\chi }_n^i}{{\chi }_n^i}=\frac{{\overline{\chi }}_n^i}{{\chi }_n^i}-1,n=1,2,\cdots ,{\overline{N}}^i.\end{array}$$

(23)

The value of ${\zeta }_n^i$ implies a difference between ${\overline{\chi }}_n^i$ in Equation (18) and ${\chi }_n^i$ in Equation (10), and thus, it can provide information about how many modes are needed for accurate reduction of flexible bodies. However, the computation of ${\zeta }_n^i$ requires a huge burden because the reference eigenvalues, ${\chi }_n^i$, are obtained by solving the eigenvalue problems of the original flexible bodies. To overcome this issue, Kim et al. has proposed a posteriori error estimation techniques that offer a way to directly estimate ${\zeta }_n^i$ without the reference eigenvalues (${\chi }_n^i$) [22,23]. The previous error estimators provide accurate eigenvalue error prediction, but the estimated errors are generally less than the exact errors [23]. This feature may cause excessive and inaccurate model reduction. In this work, a more reasonable error estimator toward the upper bound of the exact error is proposed from a conservative point of view.

Premultiplying ${{\mathit{\phi }}_n^i}^T$ by the eigenvalue problem of the original flexible body in Equation (10), the following scalar equation could be obtained:

$$\begin{array}{c}\hfill \frac{1}{{\chi }_n^i}{{\mathit{\phi }}_n^i}^T{\mathbf{K}}_F^i{\mathit{\phi }}_n^i={{\mathit{\phi }}_n^i}^T{\mathbf{M}}_F^i{\mathit{\phi }}_n^i,n=1,2,\cdots ,N_F^i.\end{array}$$

(24)

The original eigenvector ${\mathit{\phi }}_n^i$ is then decomposed to the approximated eigenvector and the error term as:

$$\begin{array}{c}\hfill {\mathit{\phi }}_n^i={\overline{\mathit{\phi }}}_n^i+\delta {\mathit{\phi }}_n^i.\end{array}$$

(25)

The approximated eigenvector ${\overline{\mathit{\phi }}}_n^i$ of the flexible body is defined in Equation (19). Considering the residual modal compensation, ${\overline{\mathit{\phi }}}_n^i$ could be more accurately approximated [4,10] as:

$$\begin{array}{c}\hfill {\overline{\mathit{\phi }}}_n^i={\overline{\mathbf{T}}}_1^i{\mathsf{\Xi }}_n^i,{\overline{\mathbf{T}}}_1^i={\overline{\mathbf{T}}}_0^i+{\omega }^2{\overline{\mathbf{T}}}_r^i,{\overline{\mathbf{T}}}_r^i=\left[\begin{array}{cc}\mathbf{0}& {\mathbf{F}}_{rs}^i{\mathbf{A}}^i\\ \mathbf{0}& \mathbf{0}\end{array}\right],\end{array}$$

(26a)

$$\begin{array}{c}\hfill {\mathbf{A}}^i={\mathbf{M}}_s^i{\mathsf{\Psi }}^i+{\mathbf{M}}_c^i,{\mathbf{F}}_{rs}^i={\mathsf{\Phi }}_r^i{\left[{\mathsf{\Lambda }}_r^i\right]}^{-1}{{\mathsf{\Phi }}_r^i}^T,\end{array}$$

(26b)

where $\omega$ is an angular frequency. ${\mathbf{F}}_{rs}^i$ is a flexibility matrix including the residual modes that are defined in Equation (15). ${\mathsf{\Lambda }}_r^i$ is a diagonal matrix of the residual eigenvalues with respect to the residual eigenvectors ${\mathsf{\Phi }}_r^i$. Equation (26a) shows that the approximated eigenvector in Equation (19) can be refined with the additional transformation matrix ${\overline{\mathbf{T}}}_r^i$, which provides better representation of the flexible mode shapes [4,10]. The derivation details are presented in Appendix A.

Substituting Equation (26a) into Equation (25), we have:

$$\begin{array}{c}\hfill {\mathit{\phi }}_n^i={\overline{\mathbf{T}}}_1^i{\mathsf{\Xi }}_n^i+\delta {\mathit{\phi }}_n^i.\end{array}$$

(27)

Substituting Equation (27) into Equation (24), the original eigenvector ${\mathit{\phi }}_n^i$ can be represented as:

$$\begin{array}{c}\hfill \frac{1}{{\chi }_n^i}{[{\overline{\mathbf{T}}}_1^i{\mathsf{\Xi }}_n^i+\delta {\mathit{\phi }}_n^i]}^T{\mathbf{K}}_F^i[{\overline{\mathbf{T}}}_1^i{\mathsf{\Xi }}_n^i+\delta {\mathit{\phi }}_n^i]={[{\overline{\mathbf{T}}}_1^i{\mathsf{\Xi }}_n^i+\delta {\mathit{\phi }}_n^i]}^T{\mathbf{M}}_F^i[{\overline{\mathbf{T}}}_1^i{\mathsf{\Xi }}_n^i+\delta {\mathit{\phi }}_n^i].\end{array}$$

(28)

If the approximated eigenvector is close to the exact eigenvector in Equation (25), the error term ($\delta {\mathit{\phi }}_n^i$) becomes small, and then, all the terms with respect to $\delta {\mathit{\phi }}_n^i$ could be neglected to approximate Equation (28). In addition, using Equation (26a), Equation (28) is written as:

$$\begin{array}{c}\hfill \frac{1}{{\chi }_n^i}{{\mathsf{\Xi }}_n^i}^T{[{\overline{\mathbf{T}}}_0^i+{\omega }^2{\overline{\mathbf{T}}}_r^i]}^T{\mathbf{K}}_F^i[{\overline{\mathbf{T}}}_0^i+{\omega }^2{\overline{\mathbf{T}}}_r^i]{\mathsf{\Xi }}_n^i\approx {{\mathsf{\Xi }}_n^i}^T{[{\overline{\mathbf{T}}}_0^i+{\omega }^2{\overline{\mathbf{T}}}_r^i]}^T{\mathbf{M}}_F^i[{\overline{\mathbf{T}}}_0^i+{\omega }^2{\overline{\mathbf{T}}}_r^i]{\mathsf{\Xi }}_n^i.\end{array}$$

(29)

Using the mass orthonormality and stiffness orthogonality conditions of the reduced flexible body (${{\mathsf{\Xi }}_n^i}^T{\overline{\mathbf{K}}}_F^i{\mathsf{\Xi }}_n^i={\overline{\chi }}_n^i$ and ${{\mathsf{\Xi }}_n^i}^T{\overline{\mathbf{M}}}_F^i{\mathsf{\Xi }}_n^i=1$), Equation (29) is expressed as:

$$\begin{array}{c}\hfill \frac{{\overline{\chi }}_n^i}{{\chi }_n^i}-1\approx {{\mathsf{\Xi }}_n^i}^T\left[{\omega }^2{{\overline{\mathbf{T}}}_0^i}^T{\mathbf{M}}_f^i{\overline{\mathbf{T}}}_r^i+{\omega }^2{{\overline{\mathbf{T}}}_r^i}^T{\mathbf{M}}_f^i{\overline{\mathbf{T}}}_0^i+{\omega }^4{{\overline{\mathbf{T}}}_r^i}^T{\mathbf{M}}_f^i{\overline{\mathbf{T}}}_r^i\right]{\mathsf{\Xi }}_n^i\\ \hfill -\frac{1}{{\chi }_n^i}{{\mathsf{\Xi }}_n^i}^T\left[{\omega }^2{{\overline{\mathbf{T}}}_0^i}^T{\mathbf{K}}_f^i{\overline{\mathbf{T}}}_r^i+{\omega }^2{{\overline{\mathbf{T}}}_r^i}^T{\mathbf{K}}_f^i{\overline{\mathbf{T}}}_0^i+{\omega }^4{{\overline{\mathbf{T}}}_r^i}^T{\mathbf{K}}_f^i{\overline{\mathbf{T}}}_r^i\right]{\mathsf{\Xi }}_n^i.\end{array}$$

(30)

where ${\omega }^2$ is ${\chi }_n^i$, we have:

$$\begin{array}{c}\hfill \frac{{\overline{\chi }}_n^i}{{\chi }_n^i}-1\approx 2{{\mathsf{\Xi }}_n^i}^T{{\overline{\mathbf{T}}}_0^i}^T\left[{\chi }_n^i{\mathbf{M}}_F^i-{\mathbf{K}}_F^i\right]{\overline{\mathbf{T}}}_r^i{\mathsf{\Xi }}_n^i+{{\mathsf{\Xi }}_n^i}^T{{\overline{\mathbf{T}}}_r^i}^T\left[{\left[{\chi }_n^i\right]}^2{\mathbf{M}}_F^i-{\chi }_n^i{\mathbf{K}}_F^i\right]{\overline{\mathbf{T}}}_r^i{\mathsf{\Xi }}_n^i.\end{array}$$

(31)

The left side of Equation (31) is the relative eigenvalue error ${\zeta }_n^i$ defined in Equation (23). Therefore, the right side of Equation (31) is a direct approximation of ${\zeta }_n^i$.

In the same manner as Equation (25), ${\chi }_n^i$ is decomposed as:

$$\begin{array}{c}\hfill {\chi }_n^i={\overline{\chi }}_n^i+\delta {\chi }_n^i.\end{array}$$

(32)

In this derivation, we assumed that the approximated eigenvector is close to the exact eigenvector in Equation (25). This also implies that $\delta {\chi }_n^i$ is very small, and then, ${\chi }_n^i$ can be approximated as ${\overline{\chi }}_n^i$. Replacing ${\chi }_n^i$ with ${\overline{\chi }}_n^i$ in the right-hand side of Equation (31), we get:

$$\begin{array}{c}\hfill {\zeta }_n^i={{\mathsf{\Xi }}_n^i}^T\left[2{\overline{\chi }}_n^i{{\overline{\mathbf{T}}}_0^i}^T{\mathbf{M}}_F^i{\overline{\mathbf{T}}}_r^i-2{{\overline{\mathbf{T}}}_0^i}^T{\mathbf{K}}_F^i{\overline{\mathbf{T}}}_r^i+{\left[{\overline{\chi }}_n^i\right]}^2{{\overline{\mathbf{T}}}_r^i}^T{\mathbf{M}}_F^i{\overline{\mathbf{T}}}_r^i-{\overline{\chi }}_n^i{{\overline{\mathbf{T}}}_r^i}^T{\mathbf{K}}_F^i{\overline{\mathbf{T}}}_r^i\right]{\mathsf{\Xi }}_n^i,\end{array}$$

(33)

and the component matrices in Equation (33) are defined by:

$$\begin{array}{c}\hfill {{\overline{\mathbf{T}}}_0^i}^T{\mathbf{M}}_F^i{\overline{\mathbf{T}}}_r^i=\left[\begin{array}{cc}\mathbf{0}& {\mathbf{X}}_a^i\\ \mathbf{0}& {\mathbf{X}}_b^i\end{array}\right],{{\overline{\mathbf{T}}}_r^i}^T{\mathbf{K}}_F^i{\overline{\mathbf{T}}}_r^i=\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{X}}_b^i\end{array}\right],{\mathbf{X}}_a^i={{\mathsf{\Phi }}_d^i}^T{\mathbf{M}}_s^i{\mathbf{F}}_{rs}^i{\mathbf{A}}^i,{\mathbf{X}}_b^i={{\mathbf{A}}^i}^T{\mathbf{F}}_{rs}^i{\mathbf{A}}^i,\end{array}$$

(34a)

$$\begin{array}{c}\hfill {{\overline{\mathbf{T}}}_0^i}^T{\mathbf{K}}_F^i{\overline{\mathbf{T}}}_r^i=\left[\begin{array}{cc}\mathbf{0}& {\mathbf{Y}}_a^i\\ \mathbf{0}& {\mathbf{Y}}_b^i\end{array}\right],{\mathbf{Y}}_a^i={{\mathsf{\Phi }}_d^i}^T{\mathbf{K}}_s^i{\mathbf{F}}_{rs}^i{\mathbf{A}}^i,{\mathbf{Y}}_b^i=[{{\mathsf{\Psi }}^i}^T{\mathbf{K}}_s^i+{{\mathbf{K}}^i}_c^T]{\mathbf{F}}_{rs}^i{\mathbf{A}}^i,\end{array}$$

(34b)

$$\begin{array}{c}\hfill {{\overline{\mathbf{T}}}_r^i}^T{\mathbf{M}}_F^i{\overline{\mathbf{T}}}_r^i=\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {{\mathbf{A}}^i}^T{\mathsf{\Phi }}_r^i{\left[{\mathsf{\Lambda }}_r^i\right]}^{-2}{{\mathsf{\Phi }}_r^i}^T{\mathbf{A}}^i\end{array}\right].\end{array}$$

(34c)

Using the definition of the constraint mode ${\mathsf{\Psi }}^i$ in Equation (13) and the orthogonality of the eigenvector matrices ${\mathsf{\Phi }}_d^i$ and ${\mathsf{\Phi }}_r^i$ yields:

$$\begin{array}{cc}\hfill {\mathbf{X}}_a^i& ={\mathsf{\Phi }}_d^i{\mathbf{M}}_s^i{\mathsf{\Phi }}_r^i{\left[{\mathsf{\Lambda }}_r^i\right]}^{-1}{{\mathsf{\Phi }}_r^i}^T{\mathbf{A}}^i=\mathbf{0},\hfill \end{array}$$

(35a)

$$\begin{array}{cc}\hfill {\mathbf{Y}}_a^i& ={\mathsf{\Phi }}_d^i{\mathbf{K}}_s^i{\mathsf{\Phi }}_r^i{\left[{\mathsf{\Lambda }}_r^i\right]}^{-1}{{\mathsf{\Phi }}_r^i}^T{\mathbf{A}}^i=\mathbf{0},\hfill \end{array}$$

(35b)

$$\begin{array}{cc}\hfill {\mathbf{Y}}_b^i& =\left[-{{\mathbf{K}}^i}_c^T{\left[{\mathbf{K}}_s^i\right]}^{-1}{\mathbf{K}}_s^i+{{\mathbf{K}}^i}_c^T\right]{\mathbf{F}}_{rs}^i{\mathbf{A}}^i=\mathbf{0}.\hfill \end{array}$$

(35c)

Substituting Equations (35a), (35b), and (35c) into Equations (34a) and (34b), Equation (33) can be rewritten using Equations (34a) and (34c) as:

$$\begin{array}{c}\hfill {\zeta }_n^i={{\mathsf{\Xi }}_n^i}^T\left({\overline{\chi }}_n^i\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{X}}_b^i\end{array}\right]+{\left[{\overline{\chi }}_n^i\right]}^2\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {{\mathbf{A}}^i}^T{\mathsf{\Phi }}_r^i{\left[{\mathsf{\Lambda }}_r^i\right]}^{-2}{{\mathsf{\Phi }}_r^i}^T{\mathbf{A}}^i\end{array}\right]\right){\mathsf{\Xi }}_n^i.\end{array}$$

(36)

In the previous work [23], the second term of the right-hand side of Equation (36) was neglected for computational efficiency of the error estimator, but it causes the predicted error to be at the low bound of the exact error. To handle this problem, a numerical treatment is considered in this work. In Equation (34c), ${\mathsf{\Lambda }}_r^i$ is a diagonal matrix that includes residual eigenvalues that are larger than the maximum dominant eigenvalue. When the eigenvalues are sorted in ascending order, the following inequality equation with the matrix norm is obtained:

$$\begin{array}{c}\hfill \frac{1}{{\kappa }_{N_d+1}^i}\parallel {\mathbf{F}}_{rs}^i\parallel <\frac{1}{{\kappa }_n^i}\parallel {\mathbf{F}}_{rs}^i\parallel ,\end{array}$$

(37)

where ${\kappa }_{N_d+1}^i$ is the first residual eigenvalue and ${\kappa }_n^i$ is among the dominant eigenvalues ($n\le N_d<N_d+1$). Therefore, the inequality in Equation (37) is clear with ${\kappa }_n^i<{\kappa }_{N_d+1}^i$.

We also found that $\parallel {\left[{\mathsf{\Lambda }}_r^i\right]}^{-1}\parallel$ becomes smaller than the case when all residual eigenvalues located at the diagonal terms of ${\mathsf{\Lambda }}_r^i$ are ${\kappa }_{N_d+1}^i$. Assuming that all of the diagonal components are ${\kappa }_{N_d+1}^i$ of ${\left[{\mathsf{\Lambda }}_r^i\right]}^{-1}$, it can be written as $(1/{\kappa }_{N_d+1}^i)\mathbf{I}$. We then have the following equation:

$$\begin{array}{c}\hfill \parallel {\mathsf{\Phi }}_r^i{\left[{\mathsf{\Lambda }}_r^i\right]}^{-2}{{\mathsf{\Phi }}_r^i}^T\parallel \le \frac{1}{{\kappa }_{N_{d+1}}^i}\parallel {\mathbf{F}}_{rs}^i\parallel .\end{array}$$

(38)

From Equations (37) and (38), we have:

$$\begin{array}{c}\hfill \parallel {\mathsf{\Phi }}_r^i{\left[{\mathsf{\Lambda }}_r^i\right]}^{-2}{{\mathsf{\Phi }}_r^i}^T\parallel \le \frac{1}{{\kappa }_{N_{d+1}}^i}\parallel {\mathbf{F}}_{rs}^i\parallel <\frac{1}{{\kappa }_n^i}\parallel {\mathbf{F}}_{rs}^i\parallel .\end{array}$$

(39)

Replacing ${\mathsf{\Phi }}_r^i{\left[{\mathsf{\Lambda }}_r^i\right]}^{-2}{{\mathsf{\Phi }}_r^i}^T$ with $(1/{\kappa }_n^i){\mathbf{F}}_{rs}^i$ in Equation (36) and using the inequality in Equation (39), the new error estimator $({\tilde{\zeta }}_n^i)$ at the upper bound of the relative eigenvalue error $\left({\zeta }_n^i\right)$ can be derived as:

$$\begin{array}{c}\hfill {\zeta }_n^i\le {\tilde{\zeta }}_n^i={{\mathsf{\Xi }}_n^i}^T\left({\overline{\chi }}_n^i\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{X}}_b^i\end{array}\right]+\frac{{\left[{\overline{\chi }}_n^i\right]}^2}{{\kappa }_n^i}\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{X}}_b^i\end{array}\right]\right){\mathsf{\Xi }}_n^i.\end{array}$$

(40)

${\tilde{\zeta }}_n^i$ can be written as:

$$\begin{array}{c}\hfill {\tilde{\zeta }}_n^i=({\overline{\chi }}_n^i+\frac{{{\overline{\chi }}_n^i}^2}{{\kappa }_n^i}){{\mathsf{\Xi }}_n^i}^T\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{X}}_b^i\end{array}\right]{\mathsf{\Xi }}_n^i,{\mathsf{\Xi }}_n^i=\left[\begin{array}{c}{\left({\mathsf{\Xi }}_n^i\right)}_s\\ {\left({\mathsf{\Xi }}_n^i\right)}_b\end{array}\right],\end{array}$$

(41)

and considering component matrices’ operation, the new error estimator ${\tilde{\zeta }}_n^i$ is simplified as:

$$\begin{array}{c}\hfill {\tilde{\zeta }}_n^i=({\overline{\chi }}_n^i+\frac{{{\overline{\chi }}_n^i}^2}{{\kappa }_n^i}){\left({\mathsf{\Xi }}_n^i\right)}_b^T{{\mathbf{A}}^i}^T{\mathbf{F}}_{rs}^i{\mathbf{A}}^i{\left({\mathsf{\Xi }}_n^i\right)}_b.\end{array}$$

(42)

For the lumped mass case (${\mathbf{M}}_c^i=\mathbf{0}$) [30], the error estimator ${\tilde{\zeta }}_n^i$ is defined as:

$$\begin{array}{c}\hfill {\tilde{\zeta }}_n^i=({\overline{\chi }}_n^i+\frac{{{\overline{\chi }}_n^i}^2}{{\kappa }_n^i}){\left({\mathsf{\Xi }}_n^i\right)}_b^T{{\mathbf{Z}}^i}^T{\mathbf{F}}_{rs}^i{\mathbf{Z}}^i{\left({\mathsf{\Xi }}_n^i\right)}_b,{\mathbf{Z}}^i={\mathbf{M}}_s^i{\mathsf{\Psi }}^i.\end{array}$$

(43)

The lumped mass matrix ${\mathbf{M}}_s^i$ is a diagonal matrix, then Equation (43) can be efficiently computed compared with Equation (42). The components ${\left(z^i\right)}_{nm}$ of the matrix ${\mathbf{Z}}^i$ can be rewritten by the following tensor notation as:

$$\begin{array}{c}\hfill {\left(z^i\right)}_{nm}={\left(m_s^i\right)}_{nn}{\left({\psi }^i\right)}_{nm},n=1,2,\cdots N_s^i,m=1,2,\cdots N_b^i.\end{array}$$

(44)

In addition, as we described, neglecting ${{\overline{\mathbf{T}}}_r^i}^T{\mathbf{M}}_f^i{\overline{\mathbf{T}}}_r^i$, the previous error estimator [23] was derived as:

$$\begin{array}{c}\hfill {\eta }_n^i={\overline{\chi }}_n^i{\left({\mathsf{\Xi }}_n^i\right)}_b^T{{\mathbf{Z}}^i}^T{\mathbf{F}}_{rs}^i{\mathbf{Z}}^i{\left({\mathsf{\Xi }}_n^i\right)}_b.\end{array}$$

(45)

Therefore, the previous and proposed error estimators may become lower and upper bounds of the exact relative eigenvalue errors, respectively. This could be investigated numerically through example problems. The proposed error estimators in Equations (42) and (43) only need the eigensolutions obtained from the reduced flexible models in Equation (18), and thus, those could be efficiently computed. It should be also noted that the error estimation performance highly depends on how close the reduced eigensolution is to the reference eigensolution because the error estimator is derived with the following two assumptions: a small eigenvector error and the use of ${\overline{\chi }}_n^i$. In practice, when the relative eigenvalue error is ≤10%, the error estimator works well [21,24,25].

It is a better way to evaluate the model reliability through the eigenvector error rather than the eigenvalue error. The modal assurance criterion (MAC) is then a typical way to compare the reference and approximated eigenvectors. However, it should be noted that plotting MAC might require the reference eigenvectors. To avoid this problem, the proposed technique indirectly estimates the eigenvalue errors without using the reference eigensolutions. In addition, the comparison of the eigenvalues, which are scalar, is much less expensive than the eigenvector operation. The trend of the accuracy of the eigenvalues follows the one of the eigenvectors since the eigenvalues are computed by using eigenvectors in the typical eigenvalue solvers such as subspace iteration and Lanczos algorithms. However, since eigenvectors are less convergent than eigenvalues, and we here use $e^i$ as 0.1% in the numerical examples for a conservative point of view.

4. Iterative Reduction Process

The reduced flexible body (RFlex) can decrease the computational cost compared to the full-ordered flexible body (FFlex), but it has a disadvantage in that the accuracy is relatively low because the residual modes are truncated. Therefore, it is necessary to construct a reliable RFlex when performing effective FMBD analysis. Using the frequency cut-off rule and the a posteriori error estimator derived in Equation (43), the iterative coordinate reduction process with a decision about the number of dominant modes ($N_d^i$) is proposed to make a reliable RFlex. The goal of the process is to find $N_d^i$ satisfying the specified relative eigenvalue error tolerance $e^i$ within the desired mode range $N_m^i$. The initial input variables are the error tolerance $e^i$, the initial number of dominant modes $N_{init}^i$, and the mode range $N_m^i$. The size of the initial RFlex is then $(N_{init}^i+N_b^i)\times (N_{init}^i+N_b^i)$. After constructing the reduced flexible model, its eigenvalue problem of the RFlex in Equation (18) should be solved to construct the diagonal matrices of the reduced flexible model, and then, those could also be used to estimate its relative eigenvalue error using Equation (43) for the n-th mode, where $n=1,2,\cdots ,N_m^i$. Using the estimated error, it is determined whether the relative eigenvalue error of a certain mode is greater than the error tolerance or not. If there is a mode with an error larger than the error tolerance, the dominant modes are added in the transformation matrix to create a new RFlex. If the n-th mode in the previous step did not satisfy the error tolerance, the iteration process starts from the n-th mode when determining the eigenvalue error of the newly created RFlex. Because adding dominant modes means improving the accuracy of RFlex, the eigenvalues before the n-th mode naturally satisfy the error tolerance. The iteration process stops when the estimated errors from the first to $N_m^i$-th modes are less then the desired error tolerance $e^i$. When the dominant modes are added $N_{iter}^i$ times, the size of RFlex becomes $(N_{init}^i+N_{iter}^i+N_b^i)\times (N_{init}^i+N_{iter}^i+N_b^i)$. Through this process, the size of RFlex that achieves the desired accuracy can be selected. The schematic diagram of the process is shown in Figure 1.

After the iterative reduction process, the reliable $N_m^i$ modes will be used to solve the EOMs in Equation (21) in the FMBD analysis. The coordinate ${\overline{\mathbf{p}}}_F^i$ is then transformed with the approximated mode shape matrix ${\tilde{\mathsf{\Phi }}}^i$ solving Equation (18) for the RFlex with $N_m^i$ dominant modes as:

$$\begin{array}{c}\hfill {\overline{\mathbf{p}}}_F^i={\tilde{\mathsf{\Phi }}}^i{\mathbf{a}}_F^i,{\tilde{\mathsf{\Phi }}}^i=\left[{\mathsf{\Xi }}_1^i,{\mathsf{\Xi }}_2^i,\cdots ,{\mathsf{\Xi }}_{N_m^i}^i\right],\end{array}$$

(46)

where ${\mathbf{a}}_F^i$ is a modal coordinate. Then, the coordinates in Equation (21) can be defined as:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{\mathbf{q}}_R^i\\ {\overline{\mathbf{p}}}_F^i\end{array}\right]={\tilde{\mathbf{B}}}^i\left[\begin{array}{c}{\mathbf{q}}_R^i\\ {\mathbf{a}}_F^i\end{array}\right],{\tilde{\mathbf{B}}}^i=\left[\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{0}& {\tilde{\mathsf{\Phi }}}^i\end{array}\right].\end{array}$$

(47)

Using $\tilde{\mathbf{B}}$, the EOMs in Equation (21) are additionally reduced as:

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{\mathbf{M}}_R^i& {\tilde{\mathbf{M}}}_{RF}^i\\ {\tilde{\mathbf{M}}}_{FR}^i& \mathbf{I}\end{array}\right]\left[\begin{array}{c}{\ddot{\mathbf{q}}}_R^i\\ {\ddot{\mathbf{a}}}_F^i\end{array}\right]+\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\tilde{\mathsf{\Lambda }}}_F^i\end{array}\right]\left[\begin{array}{c}{\mathbf{q}}_R^i\\ {\mathbf{a}}_F^i\end{array}\right]+\left[\begin{array}{c}{\mathbf{C}}_{{\mathbf{q}}_R^i}^T\\ {\tilde{\mathbf{C}}}_{{\mathbf{a}}_F^i}^T\end{array}\right]{\lambda }^i=\left[\begin{array}{c}{\left({\mathbf{Q}}_e^i\right)}_R\\ {\left({\tilde{\mathbf{Q}}}_e^i\right)}_F\end{array}\right]+\left[\begin{array}{c}{\left({\mathbf{Q}}_v^i\right)}_r\\ {\left({\tilde{\mathbf{Q}}}_v^i\right)}_F\end{array}\right],\end{array}$$

(48)

and the component matrices and vectors are:

$$\begin{array}{cc}\hfill & {\tilde{\mathbf{M}}}_{RF}^i={{\tilde{\mathbf{M}}}_{FR}^i}^T={\overline{\mathbf{M}}}_{RF}^i{\tilde{\mathsf{\Phi }}}^i,{\tilde{\mathbf{C}}}_{{\mathbf{a}}_F^i}={\overline{\mathbf{C}}}_{{\mathbf{p}}_F^i}{\tilde{\mathsf{\Phi }}}^i,{\tilde{\mathsf{\Lambda }}}_F=\mathrm{diag}\left({\overline{\chi }}_1^i,\cdots ,{\overline{\chi }}_{N_m}^i\right).\hfill \\ \hfill & {\left({\tilde{\mathbf{Q}}}_e^i\right)}_F={{\tilde{\mathsf{\Phi }}}^i}^T{\left({\overline{\mathbf{Q}}}_e^i\right)}_F,{\left({\tilde{\mathbf{Q}}}_v^i\right)}_F={{\tilde{\mathsf{\Phi }}}^i}^T{\left({\overline{\mathbf{Q}}}_v^i\right)}_F.\hfill \end{array}$$

(49)

When Equation (48) is solved for the end time $t_{end}$ of the FMBD analysis, the solutions ${}^t{\mathbf{q}}_R^i$ and ${}^t{\mathbf{a}}_F^i$ are obtained, where the prefix superscript ${}^t$ denotes the time instant for $t=t_0,t_1,\cdots ,t_{end}$. To restore the original generalized coordinate ${}^t{\mathbf{q}}^i$ from the solutions, the following equations that premultiply the transformation matrices ${\mathbf{B}}^i$ and ${\tilde{\mathbf{B}}}^i$ by the solutions should be employed:

$$\begin{array}{cc}\hfill {}^t{\mathbf{q}}^i& =\left[\begin{array}{c}{}^t{\mathbf{q}}_R^i\\ {}^t{\mathbf{q}}_F^i\end{array}\right]\approx {\mathbf{B}}^i{\tilde{\mathbf{B}}}^i\left[\begin{array}{c}{}^t{\mathbf{q}}_R^i\\ {}^t{\mathbf{a}}_F^i\end{array}\right]=\left[\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{0}& {\overline{\mathbf{T}}}_0^i\end{array}\right]\left[\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{0}& {\tilde{\mathsf{\Phi }}}^i\end{array}\right]\left[\begin{array}{c}{}^t{\mathbf{q}}_R^i\\ {}^t{\mathbf{a}}_F^i\end{array}\right].\hfill \end{array}$$

(50)

5. Numerical Examples

In this section, we report our investigation of the performance of the RFlex generation iterative process using two numerical examples. In this study, a quadruple pendulum problem and a slider-crank mechanism are considered. The FMBD simulation results were obtained using RecurDyn [29], and the external subroutine program for the iterative algorithm was developed for an RFlex generation. The generalized-$\alpha$ method [31] is used for a time integrator, where the initial step size and the maximum time step were ${10}^{-6}$ and ${10}^{-2}$, respectively. The error tolerance of the Newton–Raphson method is 5 ×${10}^{-3}$. The damping is neglected in this study.

5.1. Quadruple Pendulum Problem

A quadruple pendulum model consisting of four cylinders is considered. The pendulum consists of three rigid bodies and one flexible body, and they are connected by spherical joints sequentially. Body 1 is connected to the ground, and one side of the flexible body is modeled as a free end. The model is described in Figure 2. When the simulation starts, the bodies fall freely in the $-y$ direction due to gravitational acceleration g = 9.806 m/s${}^2$, and the chaotic behavior appears. Here, the simulation end time is seven seconds, and the number of steps is 500. The length and the radius of the cylinders are 0.3 m and 0.015 m, respectively. The distance from a spherical joint to a cylinder is 0.05 m. For the flexible body modeling, the free-free boundary condition is applied. An ideal Young’s modulus, which is half of steel’s, is 100 GPa. Poisson’s ratio v is 0.285, and the density $\rho$ is 7850 kg/m${}^3$.

The flexible body is modeled using 3528 tetrahedron elements. The average and the minimum size of the elements are 0.0075 m and 0.0056 m, respectively. The initial number of dominant modes $N_{init}^i$ is five, and the multi-point constraint (MPC) [29] node on the right side of the flexible body is selected as the boundary node. The mode range $N_m^i$ selected is 30, and the desired error tolerance is considered to be 0.1%. The dynamic responses of node A in Figure 2 (right), which were selected arbitrarily, are compared in the following numerical studies.

Figure 3 is a relative eigenvalue error comparison of the previous [23] and the new error estimators in Equation (43). Two cases with 50 and 100 dominant modes are considered here. The previous error estimator well describes the relative eigenvalue errors with high accuracy, but it generally tends to be at the low bounds of the exact error. This may lead to excessive coordinate reduction. Otherwise, the new error estimator tends to be at the upper bounds of the exact error. In this work, in a conservative manner, the new error estimator is considered, but both error estimators could be used as the lower and upper bounds of the exact relative eigenvalue error. This is because neglecting the additional term of the new error estimator directly leads to the previous one, as shown in Equations (42) and (45).

To evaluate the performance of the proposed algorithm, the estimated and exact relative eigenvalue errors change as the dominant mode is added for the n-th mode (shown in Figure 4). Among the mode range, some modes ($n=1,6,10,15,16,17,19,26,$ and 30) are only considered in Figure 4 for better readability. We found that the accuracy of lower modes shows fast convergence adding dominant modes than the accuracy of higher modes in Figure 4, which follows a typical tendency of the mode superposition. To satisfy the error tolerance $e^i=0.1\%$ within the mode range of $N_m^i$ = 30,446 dominant modes are selected in the algorithm. It can be confirmed that the exact relative eigenvalue error of the RFlex using the error estimator satisfies the specified error tolerance for the n-th mode. The sufficiently accurate RFlex that satisfies the error tolerance $e^i$ for the desired mode range can be obtained by using the newly derived error estimator at the upper bounds.

Here, dynamic responses are compared for three RFlex models that have different error tolerances $e^i=10\%,1\%$, and $0.1\%$. The selected number of dominant modes ($N_d^i$) is 35 and 58 for the error tolerances $e^i=10\%$ and $1\%$, respectively. Although the number of dominant modes selected is different, the size of the reduced flexible bodies for the error tolerances are identical to the one we described in Equation (48). The position, velocity, and acceleration of node A in Figure 2 are plotted in Figure 5, Figure 6 and Figure 7, respectively. For the position $x,y$, and z in Figure 5, the absolute error $\left|\delta \mathbf{x}\right|$ of the accelerations are calculated as $\left|\delta \mathbf{x}\right|=$$|{\mathbf{x}}^{ref}-{\mathbf{x}}^{e^i}|$, where the superscripts $ref$ and $e^i$ denote the solutions from the reference model and the RFlex model generated using the error tolerance $e^i$, respectively. It can be seen that the smaller error tolerance leads the smaller the absolute error of the dynamic response. It can be also found in the von Mises stress and strain in Figure 8.

We compared the absolute error of the von Mises stress and strain for the specific time instant in Figure 9 and Figure 10, where the selected time instant are $t=6.04$ s and $t=6.8$, respectively. The snapshots of the contours are selected at the arbitrary time. The absolute error of von Mises stress described in Figure 9 shows the area with a large error becoming smaller when decreasing error tolerance. In the case of the von Mises strain, comparable errors occur when the error tolerance is 10% and 1%, but the RFlex of $e^i=0.1\%$ gives a similar strain result to the reference model in the overall area.

5.2. Slider-Crank Mechanism

Let us consider the slider-crank system shown in Figure 11 (left) with the flexible connecting rod (right). When torque is imposed on the crank, the slider reciprocates up and down inside the cylinder. The crank is constrained to the ground and the connecting rod with a revolute joint, and the connecting rod is connected to the pin. The inner cylindrical surfaces of the connecting rod are modeled using MPC (see Figure 11, middle), and the nodes on the center of the cylindrical surfaces are the connecting points of the joints. The revolute joint is also used to pin and the slider moves only in the x-direction due to the translational joint with the cylinder, where the cylinder is fixed to the ground. The dimensions of the bodies are described in Figure 11 and Figure 12. Free-free boundary conditions are applied for the flexible connecting rod. The Young’s modulus E, Poisson’s ratio $\nu$, and density $\rho$ of the flexible connecting rod are $E=70$ GPa, $v=0.33$, and 2710 kg/m${}^3$, respectively, which are the material properties of aluminum alloy 6061-T6. The total mass of the system is measured from the CAD kernel of RecurDyn [29] based on the geometries as 5.407 kg. The number of tetrahedron elements used to model the flexible connecting rod is 4380. The average and the minimum size of elements are 0.0036 m and 0.001 m, respectively. The driving torque T shown in Figure 11 (right) is applied to the crank, which increases to 1 Nm for up to 0.1 s and remains until the simulation end time 0.5 s. The initial number of the dominant mode $N_{init}^i$ is five, and the selected mode range $N_m^i$ is 10. We select node A, which represents the maximum stress in Figure 11 (middle) to compare the dynamic responses.

The relative eigenvalue errors from the previous [23] and the new error estimator in Equation (43) are compared in Figure 13. In this example, the cases of the number of dominant modes are 30 and 60. The results show that the new error estimator is at the upper bounds in all mode ranges in both cases.

The proposed algorithm is applied for the desired mode range $N_m^i=10$ and the error tolerance $e^i=0.1\%$. For Modes 4, 5, and 6, the changes of the relative eigenvalue error and its estimation when the number of dominant modes is added are shown in Figure 14.

The time transient accelerations ${\mathbf{a}}_x$ and ${\mathbf{a}}_y$ of node A for the error tolerances $e^i=10\%,$$1\%$, and $0.1\%$ are shown in Figure 15. The selected number of dominant modes ($N_d^i$) for the error tolerances $10\%$, $1\%$, and $0.1\%$ are 8, 25, and 105, respectively, but the reduced flexible models with $e^i=10\%,$$1\%$, and $0.1\%$ are identical to 10 by 10 matrices to provide a fair comparison study. The results on the left column of Figure 15 shows the absolute error $\left|\delta \mathbf{a}\right|$ of the accelerations calculated as $\left|\delta \mathbf{a}\right|=$$|{\mathbf{a}}^{ref}-{\mathbf{a}}^{e^i}|$. It can be seen that the errors of the model with $e^i=10\%$ and $1\%$ are bigger than in the case of $e^i=10\%$, whereas in the case of $e^i=0.1\%$, the error is close to zero for the overall simulation time.

The von Mises stress and strain are compared in Figure 16 in the same manner. The absolute errors of the von Mises stress and strain show that the RFlex model with $e^i=0.1\%$ leads to a more accurate solution for the reference model than the models with the error tolerances $e^i=10\%$ and $e^i=1\%.$ In order to confirm the trend in the global area, the absolute error contours for von Mises stress and von Mises strain were compared for a specific time instant. Figure 17 shows the absolute stress error for different errors $e^i$ at $t=0.4495$ s, and Figure 18 shows the strain error at $t=0.207$ s. We found that local errors of the connecting rod, which has a relatively complicate geometry, become more significant than the first simple cylinder problem. The snapshots of the contours are selected at the arbitrary time.

It should also be noted that the decisions of the desired error tolerance and the target mode range, which are problem dependent parameters, are not covered in the the proposed iterative algorithm.

6. Conclusions

In this work, an iterative coordinate reduction algorithm for fast FMBD simulation is introduced. In this proposed algorithm, the dominant flexible modes within the initially determined number are selected based on the well-known frequency cut-off rule. Using the selected dominant modes, the flexible model is reduced, and then, a novel a posteriori error estimator is employed to evaluate its reliability. The error estimator newly derived here provides the upper bound of the relative eigenvalue error. If the estimated error does not satisfy the desired error tolerance, the above process is iterated with the increase in the number of dominant modes. The iterative algorithm with the theoretical mode selection and error estimation criteria can overcome the lack of consistent coordinate reduction in the empirical approaches. The proposed algorithm is developed for the CB-based coordinate reduction of the FFRF, which has been widely used for the FMBD analysis. However, the main idea could be employed for various coordinate reduction techniques with those mode selection and error estimation methods. It should also be noted that the proposed algorithm is based on the linear elastic motion of the flexible bodies and inherits the independent coordinate reduction of the flexible bodies without considering the nonlinearity of the total system dynamics. Therefore, the accuracy of the FMBD simulation using the proposed algorithm may not be enough in highly nonlinear problems. To handle this issue, developing robust mode selection and error-estimation techniques based on nonlinear system dynamics is a prerequisite.