Nonlinear model order reduction for flexible multibody dynamics: a modal derivatives approach

An effective reduction technique is presented for flexible multibody systems, for which the elastic deflection could not be considered small. We consider here the planar beam systems undergoing large elastic rotations, in the floating frame description. The proposed method enriches the classical linear reduction basis with modal derivatives stemming from the derivative of the eigenvalue problem. Furthermore, the Craig–Bampton method is applied to couple the different reduced components. Based on the linear projection, the configuration-dependent internal force can be expressed as cubic polynomials in the reduced coordinates. Coefficients of these polynomials can be precomputed for efficient runtime evaluation. The numerical results show that the modal derivatives are essential for the correct approximation of the nonlinear elastic deflection with respect to the body reference. The proposed reduction method constitutes a natural and effective extension of the classical linear modal reduction in the floating frame.

reference (CFR), is the ability to naturally allow Model Order Reduction (MOR): the local generalized coordinates can be expressed as a linear combination of a small number of modes.
The traditional FFR approach combined with linear elastic finite element (FE) models was illustrated by Shabana [2]. This formulation has been used to numerous problems featuring large rigid body displacement but small deflections. However, nonlinear effects due to elastic geometric nonlinearities are not incorporated in this method, and cannot be ignored in many FMBS applications. In [3,4], the classical geometric stiffness, obtained from an expression for the strain energy that includes only some higher-order terms of the strain tensor, was included in the motion equations. This approximation ignores the foreshortening displacement, and may lead to diverging solutions in applications involving large deflections and large axial forces [5]. Mayo [6] extended this formulation and obtained additional geometric stiffness matrix and nonlinear elastic force vectors. The inclusion of this effect improved both the axial and transverse response.
The FE discretization of the elastic bodies in FFR introduces a large number of degrees of freedom (dofs), and the simulation of the multibody system becomes computationally expensive, especially when the internal forces are nonlinear. Therefore, an essential step in the modeling of FMBS is the reduction of the elastic dofs. While the traditional linear MOR methods have been widely applied in FFR formulation [2,3], some other non-modal model reduction techniques have also been used in large scale industrial models in the last few years [7][8][9]. However, efficient reduction techniques with elastic geometrical nonlinearities in FFR formulation still remain a relevant research topic, given the broad range of applications tackled by FMBS. One proposed approach, under the name of ad-hoc modes, is to specifically select some axial vibration modes, in addition to low-frequency bending and torsion modes, and include them in the reduction basis in order to properly account for the nonlinear membrane response [10,11]. However, the frequencies of fundamental axial modes are normally much higher than the ones associated to bending modes. The extraction of such modes is difficult and expensive, and therefore not practical for realistic applications. Along this line, Holm-Jørgensen [12] extended the truncated modal basis with a quasi-static correction, assuming that the high-frequency elastic modes only cause quasi-static displacements. Similar approaches were previously proposed by Schwertassek [13,14]. A set of suitable quasi-comparison function, by combining eigenfunctions and static modes, is used in the FFR formulation applying Ritz method. This method was further applied to the deployment of a solar panel array by Wallrapp [15].
Higher-order modes, also known as modal derivatives (MDs), have been proved to be an efficient approach to enrich the modal basis and represent the effects of nonlinearity [16][17][18]. This method has been successfully applied in the inertial frame description to solve nonlinear problems without large rigid body motion. Interesting applications are also found in computer graphics and haptics [19,20]. Recently, this method has also been extended from the planar beam element to a general three-dimensional shell element implementation in an inertial frame [21].
In this paper, a reduction method based on the enrichment of reduction basis constituted of vibration modes with modal derivatives is presented. The FFR formulation is considered here. The elastic nonlinearity is modeled by employing the full quadratic Green strain expression, and the MDs can be interpreted as a static correction of the selected vibration modes that correctly represent the nonlinear forces with respect to the body reference. The proposed technique is implemented with the Craig-Bampton method on the floating frame [22,23], to give an exact compatibility at boundaries. The effectiveness of the proposed method will be illustrated through several numerical examples.

Floating frame of reference formulation
The proposed method is illustrated by planar beam systems featuring multiple components. The following assumptions have been made: 1. We consider only Euler-Bernoulli beam theory; 2. The material nonlinearities are not taken into consideration; 3. Damping is neglected.

Kinematic description
In the FFR description, the absolute motion of an arbitrary point on element P j of body S i is described as the superposition of the motion of the body coordinate X i Y i and the position of the points with respect to the body reference, as shown in Fig. 1. The definition of body coordinate X i Y i is not unique, we adopt here the nodal fixed frame [24]: the origin of the body coordinate is fixed to one node of the body S i , and the O i X i axis connects the origin to the end node. The position vector r ij on body S i can be defined in FFR formulation as where R i represents the position of origin of body system X i Y i with respect to global system XY , A i is the transformation matrix from X i Y i to XY , N ij are the FE shape functions, e ij 0 is the nodal coordinate vector in the undeformed state and q ij f is the vector of elastic deformation at the nodal points. In the presented 2D framework, the transformation matrix A i depends only on the angle θ i . The kinetic energy T ij for element P j on body S i is defined as where the velocityṙ ij is given bẏ andȦ i =θ i e Z × R i , e Z is the unit vector associated with the global axis Z. Global quantities (therefore not equipped with index j ) are obtained with standard FE assembly. The equations of motion for the entire multibody system can be derived from Lagrange's equations where T and U are, respectively, the kinetic energy and strain energy for the entire system; C q is the constraint Jacobian matrix, obtained from the constraint equation C(q) = 0; λ is the vector of Lagrange multipliers, representing the constraint forces; Q e is the vector of externally applied forces; t indicates time; q is the vector of the body generalized coordinates, formed as q T = [q 1 T . . . q N T ]. Here, N is the number of bodies forming the system. Furthermore, q i can be divided as where q i r ∈ R 3 refers to the displacement and orientation of the body coordinate S i , q i f ∈ R n refers to the elastic displacement in the body coordinate S i , and n is the corresponding number of elastic dofs. For the remainder of this paper, we drop the superscript i for the sake of clarity.
The inertia coupling between the body coordinates q r and elastic coordinates q f leads to inertia terms which are configuration and velocity dependent where M is the configuration-dependent mass matrix, and Q v is the quadratic velocity vector. The details of the derivation can be found in Appendix A.1.

Nonlinear strain expression
For a planar Euler-Bernoulli beam, the strain energy U could be written as [25] where E is Young's modulus; A is the cross-sectional area of the beam; L is the length, and ε xx is the axial strain. The axial strain of an Euler-Bernoulli beam is given by Green-Lagrange strain expression [26] where u and v are, respectively, the axial and transverse displacement at any point of the cross-section. Furthermore, the displacement field can be typically specified by the following linearized form where subscription (·) 0 indicates the position at the centroid of the cross-section, and y is the transverse coordinate. We adopt here the moderate rotations, von Kármán kinematic, which states that axial deformations and curvature are small compared to the bending rotation Therefore, the quadratic strain expression applied in this work is suitable for the model with small to moderate deflections [27] ε xx = ∂u ∂x In case of isotropic linear elastic material, the strain energy could be expressed as where I is the moment of inertia. The classical linear strain energy expression, which is commonly applied in the FFR, only contains first and second order integrals in (12). The last two integrals in (12) are cubic and quartic functions of displacements, respectively. Once the FE discretization and assembly are applied, the elastic forces here can be directly generated from the differentiation of strain energy U as where the internal force vector Q nl is a third-order polynomial function of the displacements q f , and it can be divided into two contributions: the first term contains the classical linear internal forces, while the second term Q f contains the higher order geometrically nonlinear terms.

Equations of motion
The equations of motion could be obtained by simply substituting (6) and (13) into (4), and written as or more compactly, where f = Q e + Q v . The explicit dependency on q is here dropped for clarity. Equation (15) can be conveniently written in a partitioned form in terms of a coupled set of reference and elastic coordinates: In multibody dynamics, the constraints are often differentiated twice with respect to time and incorporated in the inertial terms [2,3]. Then, additional regularization for constraint equations is required to ensure that the constraints are satisfied on the displacement and velocity level. This constraint regularization could be avoided by solving the original constraints (usually nonlinear) together with the equations of motion Since the constraints will generally introduce infinite stiffness into the system, it is necessary to use unconditionally stable time integration schemes [12]. Usually, the constraints are acting only on specific boundary nodes. It is therefore convenient to further partition the elastic dofs as where the subscripts (·) b and (·) i replace (·) f to represent boundary and internal dofs. q b ∈ R n b , q i ∈ R n i and n = n b + n i . The constraint Jacobian matrix can then be written as where the zero block matrix 0 qi indicates that no constraints are acting on the internal dofs. Equation (16) can then be partitioned in the same fashion as

Nonlinear model order reduction
In the classical application of modal analysis in FFR [3], the elastic displacement field q f is represented as a linear combination of mode shapes. In this section, we discuss the extension of the well-known Craig-Bampton method for the effective reduction of flexible multibody systems characterized by elastic geometric nonlinearities.

Craig-Bampton method
In order not to introduce any error in the constraints during the modal transformation, we adopt here the well known Craig-Bampton method [28] for the reduction of the internal elastic dofs. The reduction can be written as where constraint modes ∈ R n i ×n b are used to account for local effects at boundaries. In the linear modal analysis, the fixed interface modes only contains m vibration modes (VMs) of the system when constrained at the interface (i.e., q b = 0), solution of the eigenvalue problem where ω j is the j th eigenfrequency and φ j is the associated VM. In (20), η ∈ R m is a vector of modal coordinates. The reduction will be achieved by forming only with m n i internal VMs. In practice, the reduction is performed on each component of the multibody system independently.

Augmented reduction basis with modal derivatives
Although Craig-Bampton method has been successfully applied in the FFR formulation in [22], the reduction basis in (20) is usually linearized around a reference equilibrium position.
In geometrically nonlinear systems, linear VMs usually fail to accurately reproduce the motion because of their inability to account for bending-stretching coupling caused by finite deflections. Typically, low-frequency bending dominated VMs must be accompanied by axial modes to accommodate for such effects. Axial VMs can in principle be calculated. However, their extraction is expensive since they are typically associated to much higher frequencies with respect to the bending modes. In addition, for more complex and realistic components, the distinction between axial and bending/twisting dominated modes can be difficult to establish.
To overcome this difficulty, modal derivatives (MDs) stemming from VMs can be appended to the existing linear reduction basis , as shown in [16,17]. When the internal displacement q i cannot be considered small, we first assume a nonlinear mapping between the modal coordinate vector η and internal displacement vector q i , which gives where X(η) is a configuration-dependent matrix of modes. The mapping (23) can be expanded in Taylor series around the equilibrium position where the derivatives of the displacement vector q i with respect to the modal amplitudes η j around the equilibrium position can be computed from (24) as and where (26) represents the linear VM φ j (i.e., calculated at the equilibrium), and (27) gives the corresponding MD, denoted here as θ jk , which represents how X j changes because of an imposed perturbation in the shape X k , at equilibrium. Equation (25) can therefore be written as The MDs can be computed analytically by differentiating the linearized eigenvalue problem with respect to the modal amplitudes: It has already been shown that the inertia related terms in (29) can be neglected [16][17][18]. In this case, the MDs could be interpreted as a static correction, and are calculated by solving the linear system Since the system is fixed at its boundary nodes by applying a nodal-fixed frame in FFR, the internal stiffness matrix K ii is nonsingular. By neglecting all the inertia terms in (29), it can be proved that MDs are symmetric, i.e., θ jk = θ kj . Therefore, given m VMs, r = m(m+1)/2 MDs can be calculated.
In order to properly capture the contribution of the nonlinearity, we now augment the reduction basis for the internal dofs with a set of MDs collected in the matrix , by loosing the quadratic mapping (28) between VMs and MDs and adding additional modal coordinates ξ , Note that, although the number of MDs r grows quadratically with the number of chosen VMs m, it is possible to use simple selection criteria that indicate the most significant k MDs for the given analysis [29]. The reduction basis for an arbitrary body S i is therefore written as Once the reduction basis (32) has been derived, the equations of motion (19) can be evaluated and projected to obtain a model of greatly reduced dimensions. Unfortunately, this procedure is inefficient since the cost for the evaluation of the nonlinear terms (i.e., inertial and elastic forces) scales with the size of the nodal coordinates: the nonlinear terms have to be computed in nodal coordinates and then projected on the reduced subspace. In our case, the nonlinear terms are written directly in terms of the modal coordinates. This is discussed in detail in the next section.

Precomputing polynomial coefficients
The adopted kinematic model yields a multivariate third-order polynomial expression of the nonlinear elastic forces Q f . A generic component Q I f can be written as where Q 1 ∈ R n×n×n , Q 2 ∈ R n×n×n×n are constant third order and fourth order tensor coefficients with generic components Q I ij 1 and Q I ijl 2 . We adopted here the Einstein summation convention over the repeated indexes in the superscript.
Consequently, with the modal transformation of the form q f = Vγ , where and the reduced internal forces Q f = V T Q f (Vγ ) is a multivariate cubic polynomial in reduced coordinates γ , where Q 1 ∈ R v×v×v , Q 2 ∈ R v×v×v×v are constant third order and forth order tensor coefficients in the reduced coordinates with v = n b + m + k. Similarly, each component of the reduced tangent stiffness matrix K nl is also a multivariate quadratic polynomial in γ , with constant tensor coefficients K 3 ∈ R v×v×v , K 4 ∈ R v×v×v×v . The constant linear stiffness can be expressed as The tensors Q 1 , Q 2 , K 3 , K 4 and K L can be precomputed offline for efficient runtime evaluation.

Reduced equations
The reduced equations are obtained via a classical Galerkin projection, i.e., the residual obtained by introducing (32) in (16) is projected onto the same reduced basis used for the approximation of the displacements. As discussed in the previous section, it is possible to precompute the nonlinear terms to directly obtain modal terms rather than performing the full evaluation and projection. The resulting reduced equations are written as where Note that C qf only contains nonzero terms on the column corresponding to the boundary dofs q b . Because of the inertia coupling between the reference motion of the body and the elastic deformation of elements, also the inertia related terms M rf , M f r and f f will be configuration-dependent. Similar to the elastic terms, the inertial terms can also be directly expressed in modal coordinates. The detailed formulation of the reduced mass matrix and quadratic velocity vector is reported in Appendix A.2.
The detailed MOR procedure proposed is illustrated in Fig. 2 on a flexible slider-crank mechanism. The end nodes in both crank and connecting rod are treated as boundary nodes: their corresponding dofs will not be transformed to modal coordinates. The material and geometrical properties are taken from the numerical example in Sect. 5.3. From Fig. 2 we can see that the first few VMs feature bending displacements only, while the corresponding MDs, which describe the second-order nonlinearities, exhibit only longitudinal displacement.

Computational cost estimation for reduced time integration
In this section, we present a complexity analysis to highlight the benefits of the proposed method in terms of computational cost. The number of dofs N 1 = 3 + n b + n i relative to the full analysis is reduced to N 2 = 3 + n b + m + k where m + k n i . In the present work, the numerical solution of the second order differential equations resulting from the reduced order model is performed with the implicit Newmark method, with time integration parameters β = 0.25 and α = 0.5. A similar scheme has already been successfully applied in FFR with full models [30], and now it will be applied to reduced models.
If we express the computational cost in terms of total number of dofs in the full and reduced models, i.e., N 1 and N 2 , respectively, then a detailed comparison can be summarized in Table 1 for a single iteration step.
For a fair estimation, the computational cost involved in this table also includes the time spent on the evaluation of the reduction basis and polynomial coefficients, sometimes also referred as offline cost (steps 1 and 2). The two most time-consuming operations are tangent stiffness and force assembly (steps 4 and 5), as well as the solution for displacement increment (step 6), which involves the factorization of large sparse symmetric matrices.
In this work, the effort in steps 4 and 5 is reduced by directly expressing the nonlinear force vector as a function of modal coordinates, with precomputing polynomial coefficients as shown in Sect. 3.3.
In addition, the reduced model here results in a much smaller system of equations to be solved for the displacement increment (step 6) during every corrector step in the time integration scheme. In Table 1 Table 1 clearly indicates the potential time savings if the number of reduced dofs N 2 grows much more slowly than the number of full dofs N 1 .

Numerical examples
Three numerical examples are presented in this section to show the effectiveness of the proposed approach. The reduced solutions in the FFR are not only compared with the corresponding full solutions, but also with the response obtained by using a Corotational Frame of Reference (CFR) formulation. The CFR is a more general and expensive framework that is able to deal with an arbitrary large elastic deflection. This allows assessing that the magnitudes of the elastic deflections are within the validity of the adopted approximated von Kármán kinematic model. We considered here the CFR formulation discussed in [31] as a reference, where cubic shape functions are used to derive both the inertia and elastic terms. Recently, this formulation has also been successfully extended to the dynamic analysis of 3D flexible beam with good accuracy [32].

Test 1: rotating beam
The dynamic analysis of a rotating beam, which has been used as a benchmark in many papers dealing with flexible beams and geometric nonlinearities [33][34][35][36], is here presented. The geometry of the beam and the corresponding material properties are shown in Fig. 3. An imposed end rotation θ s is applied as The beam reaches steady state motion after T s seconds and then rotates at a constant angular velocity ω s . Transient responses are computed for T s = 15 s, and for angular velocity ω s = 6 rad/s.
In order to make a comprehensive comparison, this example is analyzed by different approaches: (i) the CFR is taken as reference solution; (ii) the nonlinear floating frame is mentioned as the full analysis (denoted as NLFFR); (iii) the linear floating frame is performed by neglecting the nonlinear internal forces Q f in equation of motion (denoted as LFFR); (iv) the MDs based reduction solution is applied in a nonlinear floating frame, where the modal basis is composed of the first 10 VMs and 10 MDs corresponding to the first 4 VMs (denoted as 10 VMs + 10 MDs); (v) a reduced solution obtained with the first 20 VMs (denoted as 20 VMs) is also calculated. The two components of the tip displacement vector are shown in Fig. 4.
The results obtained with the NLFFR are in very good agreement with the CFR solution and clearly differ from the LFFR response, to confirm that the adopted kinematic model is adequate for the problem at hand; see Fig. 4. Furthermore, the proposed modal derivatives based reduction method with 10 VMs and 10 MDs shows a good agreement with the refer- To offer further insight, some other reduced models with different modal basis are compared with the full analysis (NLFFR), as shown in Fig. 5. The comparison is performed between different reduction basis of the same size. The tip displacements components are shown. In addition, the root mean square (RMS) error RMS relative to the entire displacement vector is shown, calculated as where u i , w i and u i , w i are the horizontal and vertical components of the node displacement from the full and reduced models, respectively. The reduced basis with 10 VMs enriched by 6 MDs (indicated as 10 VMs + 6 MDs) yields much better results when compared to the reduced solution obtained with 16 VMs. In addition, by increasing the reduction basis with two additional MDs (10 VMs + 8 MDs), the results show a marked improvement, while the results are almost unchanged if two additional VMs are added in the 18 VMs cases. In some related work [10,11], high-frequency axial vibration modes (AVMs) are specifically added to the basis in order to properly account for the nonlinear bending-stretching coupling. It is therefore interesting to compare this approach with the method proposed here. We form the reduction basis with first 10 VMs and the first 8 AVMs. The eigenfrequencies associated to the AVMs are reported in Table 2. Note that the frequency range of interest has to be extended from 145 Hz to almost 2100 Hz. In addition, the obtained results do not match the accuracy given by the 10 VMs + 8 MDs case. Therefore, we can conclude that in order to reproduce the nonlinear behavior, MDs based reduction basis provides better accuracy than an equal size reduction constructed with AVMs.

Test 2: swinging rubber bar
We consider here a slender rubber beam connected to a fixed inertial frame via a joint and exposed to a constant gravitational force. The beam is initially at rest in horizontal position. The material and geometric properties of the beam are taken from [37,38]: E = 5 × 10 6 N/m 2 , ρ = 1.1 × 10 3 kg/m 3 , L = 1.0 m. The cross-section is circular with a radius of 5 mm. The analysis is performed for a time interval of 1 second with a fixed time step T of 0.001 seconds.
The dynamic response is shown in Fig. 6. The full response is compared with reduced basis formed with the first 5 VMs plus 10 MDs, as well as 15 VMs, respectively. As shown in Fig. 6, a reduction basis containing 5 VMs plus 10 MDs clearly outperforms a basis of the same size formed with 15 VMs only.
The RMS error RMS of four reduced models with different modal basis is shown in Fig. 7. The number of VMs is fixed and equal to 5, while the number of MDs is increased to include the contribution of the first 2, 3 and 4 VMs, respectively. The error rapidly decreases as more MDs associated to the low frequency VMs are included in the basis.

Test 3: flexible slider-crank mechanism
In the third numerical example, a flexible slider-crank mechanism is analyzed. The system is depicted in Fig. 8 together with the material and geometric properties. The connection rod is constrained to the crank at the left end by a joint, and is fixed in vertical direction at the right end. A constant angular acceleration is prescribed to the crank such that the rotation θ s reaches 270 degrees after 5 seconds, as shown in Fig. 9(a). At this point, the angular rotation of the crank is locked and the system starts to exhibit elastic oscillations in the geometrically nonlinear range.
In the FFR model, two fixed nodal frames are attached to the crank and rod, respectively. The two end nodes of both crank and rod are set as the boundary nodes, and therefore 6 constraint modes are present in the reduced model for each substructure.
The elastic deflection of the middle point of both the crank and rod is shown in Fig. 9. The results obtained with the CFR have been set as a reference to check the accuracy of the full analysis in FFR. The reduction basis is formed with 5 VMs and 10 MDs (relative to the first 4 VMs) for both substructures, and the obtained results are compared with the case when only first 15 VMs are used.
The response of the system during the first 5 seconds (i.e., before the rotation locking) exhibits mainly rigid motion, and therefore both of the two reduced models show accurate results. After the locking, the elastic deflections are large and the reduction basis featuring only VMs is not able to reproduce the full solution, while the reduced basis formed with VMs and MDs yields very good results, as can be seen in Fig. 9.
The RMS error between the full and reduced solutions is computed also for this example in order to gain additional insight on the properties of the reduction basis, see Fig. 10. It can be noticed that, if not enough MDs are included in the reduction basis, the reduced solution is of poor accuracy. This could be attributed to the fact that the impulsive load generated by the rotation locking triggers the large response of the first few VMs, and therefore their interactions, which is given by the corresponding MDs, must be properly included in the basis.

Computational efficiency
In this subsection, the computational efficiency is compared between the nonlinear modal analysis with augmented basis (VMs + MDs) and the full nonlinear FEM analysis (NLFFR). The performance is measured in terms of computational time for Test 1: rotating beam. All the simulations are performed with MATLAB ® R2012b, on an Intel ® Core TM i5-3470 @ 3.2 GHz and 16 GB RAM machine. Table 3 compares the CPU time required for Test 1. The simulations are performed for 20 s with a time step t = 0.02 s. As shown in Table 3, the computational time in the case of nonlinear modal analysis mainly depends on the size of selected model reduction bases, and only slightly rises with the increase of element number owing to the time spend on modal transformation. Therefore, a substantial CPU time reduction can be achieved for application characterized by a large number of elements, which is typically the case for realistic applications.

Conclusions
We propose a model order reduction method for the dynamic analysis of flexible multibody systems featuring large overall motion and nonlinear elastic deflection, which can be described with the von Kármán kinematic assumption in the body reference. The equations of motion are written in the floating frame of reference (FFR) for each flexible component. This enables the description of the elastic motion with an internal reduction basis complemented with interface modes, which allow the connection between the bodies via constraints. The internal basis of linear vibration modes (VMs) is enriched with modal derivatives (MDs), which describe the essential nonlinear contributions for the elastic geometric nonlinearities induced by the VMs.
The MDs can be obtained through a sensitivity analysis of the eigenvalue problem associated to the internal VMs for each component, with respect to the modal amplitudes. The reduced basis is formed by enriching a classical Craig-Bampton reduction with MDs relative to the VMs considered for each component. Subsequently, the reduced equations of motion are obtained with a Galerkin projection.
The reduction in computational time is obtained by the substantial reduction in size of the governing equations of motion. This results in the solution of a much smaller system during the time integration via the implicit scheme. Moreover, the polynomial form of the nonlinear elastic forces allows the offline computation of the reduced nonlinear terms directly in modal coordinates.
The presented numerical examples highlight the superior performance of the proposed approach with respect to classical reduction with VMs only. It is worth noting that the adhoc inclusion of axial vibration modes in the reduction basis does not provide results as accurate as those obtained with the proposed approach. The proposed method extends the common practice of linear modal analysis to systematically tackle geometrically nonlinear multibody problems. As a consequence, the proposed technique does not require expensive sampling of the full solution to form the reduction basis, as necessary in several model order reduction techniques based on the proper orthogonal decomposition.
The proposed approach does not pose additional conceptual difficulties for the extension into three dimensional cases. In this case, it is recommended to resort on the mean-axis definition for the floating frame of reference. This choice could be more practical than the nodal fixed frame, as the best choice of the specified fixed nodes in the nodal fixed frame is not straightforward in three dimensional cases. Also, the MDs have been successfully calculated for three dimensional models [29,39] in the inertia frame, to solve nonlinear problems without large rigid body motion. It is worth mentioning that the benefits of the proposed technique will be even more apparent for three dimensional problems that would feature, in general, larger finite element meshes and therefore provide larger gains from the adopted modal approach.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: The mass matrix and quadratic velocity vector
This Appendix contains a detailed formulation of the mass matrix and quadratic velocity vector. Appendix A.1 gives all the quantities in nodal formulation in Eq. (6), while Appendix A.2 refers to the reduced coordinates (i.e., modal) formulation in Eq. (39).

A.1 Precomputed offline quantities in nodal coordinates
The configuration-dependent mass matrix M and quadratic velocity vector Q v can be written in a partitioned form as and In the following, we define the components of (A.