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

: Coordinate reduction has been widely used for efﬁcient simulation of ﬂexible multibody dynamics. To achieve the reduction of ﬂexible 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 ﬂexible 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 ﬂexible 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 veriﬁed through numerical examples. It is also evaluated that efﬁcient and accurate ﬂexible multibody dynamics simulation is available with the reduced ﬂexible body, generated by the proposed algorithm.


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.

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.

Floating Frame of Reference Formulation
The Lagrange equation for the i-th body in the FMBD system can be written as: where q i is the generalized coordinate vector including the reference (rigid) and flexible coordinates. T i is the kinetic energy. C T q i is the constraint Jacobian matrix. λ i is the vector of Lagrange multipliers. Q i is the vector of generalized forces. The generalized forces Q i acting on the i-th body can be derived using the virtual work δW i d of the elastic force and δW i e of the external forces as: where δW i d and δW i e are defined as: where K i is the stiffness matrix and Q i e is the external force vector. Substituting Equation (3) into Equation (2), the generalized force Q i can be expressed as: The terms on the left side in Equation (1) without the third term can be written as [1]: where M i is the mass matrix and Q i v is the quadratic velocity vector defined as: Using Equations (4) and (5) in Equation (1), the following EOMs neglecting the damping are obtained: and it can be specifically written in a partitioned form of the reference and flexible coordinates as: where subscripts R and F denote the quantities with respect to the reference and flexible coordinates, respectively.

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: and its eigenvalue problem is: where χ i n and ϕ i n are the n-th eigenvalue and its corresponding eigenvector, respectively. Here, N i F is the number of DOFs of the flexible body.
The mass and stiffness matrices, M i F and K i F , 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: The numbers of interior and boundary DOFs are denoted as N i s and N i b , respectively. In the CB method, the fixed-interface normal mode Φ i s and the interface constraint mode Ψ 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: where κ i n and φ i n denote the eigenvalue and eigenvector of the interior part, respectively. The interface constraint mode is then defined as [7]: The transformation matrix T i 0 can be constructed using Equations (12a) and (13), then the generalized coordinate q i F can be written with the hybrid coordinate p i F as: The generalized coordinate vector of the interior q i s is expressed as: in which the subscripts d and r denote the dominant and residual terms, respectively. The number of the dominant modes Φ i d and the residual modes Φ i r is denoted as N i d and N i r , respectively (N i . Neglecting the residual modes, the original generalized coordinate vector q i F can be approximated using the reduced transformation matrix T i 0 as: Using Galerkin projection with the transformation matrix T i 0 , Equation (9) can be reduced as: and the corresponding eigenvalue problem is: The number of DOFs for the reduced flexible body N i depends on the number of dominant modes N d and the number of boundary DOFs . The number of dominant modes is much smaller than the number of residual modes, and N i can be dramatically reduced compared to the original DOFs N i F . From the eigenvector Ξ i n , the original eigenvector ϕ i n in Equation (10) is approximated as:

Floating Frame of Reference Formulation with the Reduced Flexible Bodies
The generalized coordinate vector q i can be approximated using the transformation matrix T i 0 as: The EOMs of a reduced FMBD system can be obtained by substituting Equation (20) into Equation (8) and premultiplying by B i T as follows: where the submatrices and vectors in Equation (21) are expressed as: 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 i d may be the simplest and the most popular remedy to get more accurate χ i n .

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: The value of ζ i n implies a difference between χ i n in Equation (18) and χ i n 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 ζ i n requires a huge burden because the reference eigenvalues, χ i n , 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 ζ i n without the reference eigenvalues (χ i n ) [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 ϕ i n T by the eigenvalue problem of the original flexible body in Equation (10), the following scalar equation could be obtained: The original eigenvector ϕ i n is then decomposed to the approximated eigenvector and the error term as: The approximated eigenvector ϕ i n of the flexible body is defined in Equation (19). Considering the residual modal compensation, ϕ i n could be more accurately approximated [4,10] as: where ω is an angular frequency. F i rs is a flexibility matrix including the residual modes that are defined in Equation (15). Λ i r is a diagonal matrix of the residual eigenvalues with respect to the residual eigenvectors Φ i r . Equation (26a) shows that the approximated eigenvector in Equation (19) can be refined with the additional transformation matrix T i r , 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: Substituting Equation (27) into Equation (24), the original eigenvector ϕ i n can be represented as: If the approximated eigenvector is close to the exact eigenvector in Equation (25), the error term (δϕ i n ) becomes small, and then, all the terms with respect to δϕ i n could be neglected to approximate Equation (28). In addition, using Equation (26a), Equation (28) is written as: Using the mass orthonormality and stiffness orthogonality conditions of the reduced flexible body (29) is expressed as: where ω 2 is χ i n , we have: The left side of Equation (31) is the relative eigenvalue error ζ i n defined in Equation (23). Therefore, the right side of Equation (31) is a direct approximation of ζ i n . In the same manner as Equation (25), χ i n is decomposed as: In this derivation, we assumed that the approximated eigenvector is close to the exact eigenvector in Equation (25). This also implies that δχ i n is very small, and then, χ i n can be approximated as χ i n . Replacing χ i n with χ i n in the right-hand side of Equation (31), we get: and the component matrices in Equation (33) are defined by: Using the definition of the constraint mode Ψ i in Equation (13) and the orthogonality of the eigenvector matrices Φ i d and Φ i r yields: Substituting Equations (35a), (35b), and (35c) into Equations (34a) and (34b), Equation (33) can be rewritten using Equations (34a) and (34c) as: 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), Λ i r 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: where κ i N d +1 is the first residual eigenvalue and κ i n is among the dominant eigenvalues (n ≤ N d < N d + 1). Therefore, the inequality in Equation (37) From Equations (37) and (38), we have: (36) and using the inequality in Equation (39), the new error estimator (ζ i n ) at the upper bound of the relative eigenvalue error (ζ i n ) can be derived as: ζ i n can be written as: and considering component matrices' operation, the new error estimatorζ i n is simplified as: For the lumped mass case (M i c = 0) [30], the error estimatorζ i n is defined as: The lumped mass matrix M i s is a diagonal matrix, then Equation (43) can be efficiently computed compared with Equation (42). The components (z i ) nm of the matrix Z i can be rewritten by the following tensor notation as: In addition, as we described, neglecting T i r T M i f T i r , the previous error estimator [23] was derived as: 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 χ i n . 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.

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 i d ) is proposed to make a reliable RFlex. The goal of the process is to find N i d satisfying the specified relative eigenvalue error tolerance e i within the desired mode range N i m . The initial input variables are the error tolerance e i , the initial number of dominant modes N i init , and the mode range N i m . The size of the initial RFlex is then . 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, · · · , N i m . 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 i m -th modes are less then the desired error tolerance e i . When the dominant modes are added N i iter times, the size of RFlex becomes (N i init . 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 i m modes will be used to solve the EOMs in Equation (21) in the FMBD analysis. The coordinate p i F is then transformed with the approximated mode shape matrixΦ i solving Equation (18) for the RFlex with N i m dominant modes as: where a i F is a modal coordinate. Then, the coordinates in Equation (21) can be defined as: UsingB, the EOMs in Equation (21) are additionally reduced as: and the component matrices and vectors are: When Equation (48) is solved for the end time t end of the FMBD analysis, the solutions t q i R and t a i F are obtained, where the prefix superscript t denotes the time instant for t = t 0 , t 1 , · · · , t end . To restore the original generalized coordinate t q i from the solutions, the following equations that premultiply the transformation matrices B i andB i by the solutions should be employed:

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-α 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.

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 ρ 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 i init 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 i m 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.   (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)  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 i m = 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 i d ) 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 Figures 5-7, respectively. For the position x, y, and z in Figure 5, the absolute error |δx| of the accelerations are calculated as |δx| = |x re f − x e i |, where the superscripts re f 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 Figures 9 and 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.

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 Figures 11 and 12. Free-free boundary conditions are applied for the flexible connecting rod. The Young's modulus E, Poisson's ratio ν, and density ρ 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 i init is five, and the selected mode range N i m is 10. We select node A, which represents the maximum stress in Figure 11 (middle) to compare the dynamic responses.  The proposed algorithm is applied for the desired mode range N i m = 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 a x and 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 i d ) 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 |δa| of the accelerations calculated as |δa| = |a re f − 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.

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.  Acknowledgments: This research was supported by FunctionBay, Inc.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Enhanced Transformation Matrix
Equation (14) can be explicitly presented as: Using the transformation matrix T i 0 for the projection in Equation (9), we have: F rs can be simply computed by subtracting the dominant flexibility matrix from the full flexibility matrix, which were already computed in the original CB reduction formulation. Using Equation (A7) in Equation (A5), q s can be also approximated as: Considering the newly derived q i s , an enhanced transformation matrix T i 1 could be derived. Consequently, q i F can be approximated with high fidelity as: This clearly shows that T i 1 is a refined transformation matrix with a residual modal effect. Therefore, using T i 1 , one can expect a better representation of the modal behavior of the original flexible bodies than with the conventional CB transformation matrix T i 0 .