A modal derivatives enhanced Rubin substructuring method for geometrically nonlinear multibody systems

This paper presents a novel model order reduction technique for 3D flexible multibody systems featuring nonlinear elastic behavior. We adopt the mean-axis floating frame approach in combination with an enhanced Rubin substructuring technique for the construction of the reduction basis. The standard Rubin reduction basis is augmented with the modal derivatives of both free-interface vibration modes and attachment modes to consider the bending–stretching coupling effects for each flexible body. The mean-axis frame generally yields relative displacements and rotations of smaller magnitude when compared to the one obtained by the nodal-fixed floating frame. This positively impacts the accuracy of the reduction basis. Also, when equipped with modal derivatives, the Rubin method better considers the geometric nonlinearities than the Craig–Bampton method, as it comprises vibration modes and modal derivatives featuring free motion of the interface. The nonlinear coupling between free-interface modes and attachment modes is also considered. Numerical tests confirm that the proposed method is more accurate than Craig–Bampton’s, a nodal fixed floating frame counterpart originally proposed in Wu and Tiso (Multibody Syst. Dyn. 36(4): 405–425, [2016]), and produces significant speed-ups. However, the offline cost is increased because the mean-axis formulation produces operators with decreased sparsity patterns.


Introduction
The simulation of flexible multibody systems (FMBS) often relies on finite element (FE) discretization of flexible components, which are then embedded into a floating frame of reference (FFR) formulation [2,3]. The FFR represents the position of each body as a superposition of two components: (i) the motion of the reference frame which follows the overall rigid body motion of the flexible body; (ii) the relative motion of the flexible body with respect to the reference frame. The floating frame approach is usually preferred to the description of the multibody motion with respect to the inertial frame (see, for instance, [4]) as it naturally distinguishes the elastic deformation from the rigid body motion. The resulting models often comprise a large number of degrees of freedom (DoFs), which render time integration schemes extremely costly. A relevant example of unaffordable computational burden could be found in the simulation of large-scale offshore wind turbines. To assess their fatigue life, thousands of load cases need to be simulated, resulting in disproportionally large computation times. At present, this can be achieved only by relying on extremely simplified beam models that reduce the computational cost to a bearable level. Such models do not inherit the complexity of the actual three-dimensional model of the blade, and, as a result, the complex dynamic behavior may not be appropriately represented. For this reason, many model order reduction (MOR) strategies for three-dimensional FMBS have been proposed in the past. These techniques are based on classic modal truncation [5,6] or singular value decomposition (SVD) based MOR techniques as in [7][8][9]. In [10], a global modal parametrization based MOR method is proposed, where the motion of the FMBS is described in terms of configuration dependent modes. Using this reduction method, the nonlinear holonomic constraints are naturally satisfied without the adoption of Lagrange Multipliers. However, in most of the MOR techniques, the elastic behavior is assumed to be linear. As discussed in [5], the linear MOR with FFR formulation is only suitable for structures featuring large rigid body motions but small relative displacements with respect to the reference frame, as well as slow rotational speeds. For FMBS featuring high rigid body rotation rates, the centrifugal force is of great significance, and therefore, the centrifugal stiffening effect and foreshortening effect have to be considered.
For many FMBS applications involving finite but moderate relative rotations with respect to the reference frame, neglecting geometrical nonlinearities may lead to incorrect and even diverging solutions [11,12]. In [13], the geometrical nonlinearities are introduced in the equations of motion. As a result, the internal force vector and tangent matrix need to be recomputed for every iteration within each time step, therefore significantly impacting the computational cost. It is then a must to extend the linear MOR methods to the geometrically nonlinear regime for three-dimensional FMBS.
When one substructure of the FMBS features geometrically nonlinear behavior, dominant low-frequency modes are not sufficient for adequately representing the relative motion with respect to the reference frame. Typically, large slender structures exhibit coupling between bending and axial displacements when excited in the nonlinear regime. The corresponding bending-stretching coupling could be in principle provided by adding membrane-dominant (usually high-frequency) modes to the bending-dominant (typically low-frequency) modes based reduction basis. For flat structures, where each vibration mode exhibits purely bending or membrane displacement, such membrane modes can be easily identified and added to the reduced-order basis (ROB). The inclusion of these so-called ad hoc modes has been applied in the FFR formulation in [14,15]. However, for more complex geometries, the extraction of such modes is (i) challenging, as it is not straightforward to identify membrane-dominated modes, and (ii) expensive, as several modes need to be extracted.
In previous work [1], the linear Craig-Bampton (CB) substructuring basis [16] was enriched with modal derivatives (MDs) [17,18] corresponding to low-frequency fixedinterface modes. The augmented ROB was capable of capturing both the rigid body motions and the nonlinear relative displacement of the FMBS effectively. The nonlinear MOR technique was applied for nodal-fixed frame reference [19], which is the most straightforward implementation of the FFR formulation. In this case, the reference frame is attached to specified nodes of the moving body. However, for complex structures, e.g., discretized with shell and solid elements, it is difficult to determine the optimal node whereon the reference frame should be attached. This arbitrary definition of the nodal-fixed frame results in significantly different relative displacements and rotations with respect to the reference frame [19], and ultimately degrades the accuracy if the relative displacement and rotations are too large.
The use of mean-axis frame [20], which alleviates the need for the reference frame to be attached to a specified node of the structure, aims at minimizing the relative kinetic energy with respect to the reference frame. As a result, the largest relative displacement and rotation observed from a mean-axis frame will be smaller than the largest one observed when standing at the origin of the nodal-fixed frame, as underlined in [19]. This is especially relevant when one assumes geometrical nonlinearities based on the von Kármán kinematic assumption, which is suitable for small strains and moderate rotations [21] with respect to the reference frame. Since the MDs are obtained from a truncated Taylor expansion of the nonlinear static equilibrium around the reference position [22,23] and are not updated during the time integration, the accuracy of using MDs will be determined by how far the structure departs from the equilibrium position. Therefore, the use of MDs further supports the argument of using the mean-axis formulation.
In this paper, the standard Rubin substructuring technique [24] is enhanced with MDs and then implemented on the mean-axis frame formulation for the construction of reduced-order models (ROMs) for the FMBS featuring moderate relative displacements and rotations with respect to the reference frame. Each body is reduced by forming the ROB with attachment modes, free-interface modes, and corresponding MDs. The Rubin method fits the mean-axis formulation more naturally than the CB method when applied to the geometrically nonlinear problem, for two reasons. First, the Rubin method is based on a truncated set of free-interface vibration modes, which naturally describe the elastic deformation of the component with respect to the reference frame (i.e., free-interface deformation with respect to the reference frame as in mean-axis frame formulation). Second, the nonlinear behavior occurring at the interface is better represented by MDs of both free-interface modes and attachment modes (related to the Rubin method) than by MDs of fixed interface modes coming from the CB method. In [25], the inclusion of only the MDs relative to rigid body modes (i.e., vibration modes of zero frequency) in the ROB significantly increases the accuracy. In our approach, the MDs relative to rigid body modes are avoided since the rigid body motion has already been described by the reference frame motion. Therefore, a ROB of very limited size can be achieved.
This paper is organized as follows. Section 2 describes the FFR description featuring geometric nonlinearities. The nodal-fixed and mean-axis frame are applied to the FFR formulation in Sect. 3. The assembled EoMs of all FMBS, as well as the holonomic joint constraints, are presented in Sect. 4. The nonlinear MOR method based on the enhanced Rubin method is proposed in Sect. 5. Section 6 shows numerical examples to assess the accuracy of the present formulation, especially emphasizing the improvements with respect to [1]. Finally, conclusions are given in Sect. 7.

Equations of motion in floating frame of reference
In the FFR formulation, we describe the absolute motion of an arbitrary point P j,(s) on the j th finite element of the sth body as the superposition of the motion of the reference frame O (s) X (s) Y (s) Z (s) and the position of the point with respect to the reference frame, as shown in Fig. 1. The position vector r j,(s) ∈ R 3 of the point P j,(s) is defined as where R (s) ∈ R 3 represents the position of origin of the reference frame O (s) X (s) Y (s) Z (s) with respect to global frame OXY Z, u j,(s) ∈ R 3 is the relative nodal position of P j,(s) with respect to the reference frame, and A (s) ∈ R 3×3 is the transformation matrix from the reference frame O (s) X (s) Y (s) Z (s) to the global frame OXY Z. The matrix of shape functions in the reference frame is indicated by N j,(s) ∈ R 3×ne , where n e is the number of DoFs per element, q j,(s) 0 ∈ R ne is the vector of nodal coordinates in the undeformed state and q j,(s) f ∈ R ne is the vector of relative DoFs of the j th element. For the remainder of this paper, we drop the superscript (s) for the sake of clarify, unless it is necessary to distinguish between different bodies in the FMBS.
The rotation matrix A is defined as Here, the notation col(...) indicates the column stacking of vectors or scalar quantities. The unit vector along the rotation axis is given by v = col(v 1 , v 2 , v 3 ) and θ is the rotation angle. The axis of rotation along v and the rotation angle θ are defined for each flexible body separately. The absolute velocityṙ j is given bẏ where the argument of functional dependency is enclosed in braces. The second term in (3) is rewritten asȦ in order to isolate the velocity termsθ . The matrix B j is thus a function of θ , q j f and q j 0 . The dependency of B j on q j 0 will not be explicitly expressed since q j 0 is constant for each FE model.
The kinetic energy T j for the j th element can be evaluated by where q j = col(R, θ, q j f ) and ρ j is the density of the element material. The mass matrix M j is configuration-dependent (its formulation is given in detail in Appendix A). The kinetic energy of the sth body can be determined by summing up the kinetic energy T j of all its elements. The mass matrix M of the sth body is obtained by standard FE assembly. The vector q = col(R, θ , q f ) indicates the generalized coordinates of a single flexible body, where q f ∈ R n refers to the total relative DoFs in the reference frame and n is the number of relative DoFs.
The equations of motion (EoMs) for each flexible body can be derived from Lagrange's equations as d dt

∂T ∂q
where T and U are the kinetic energy and strain energy, respectively; g is the vector of externally applied generalized loads. At this stage, the flexible bodies (if more than one) have not been assembled and the prescribed motions of the flexible bodies have not been imposed yet. The EoMs (6) can be rewritten in matrix form as where Q is the quadratic velocity vector, which includes the effect of apparent forces (such as centrifugal force and Coriolis force), and f is the nonlinear elastic force vector. The quadratic velocity vector Q results from the inertia coupling between the motion col(R, θ ) of the reference frame and the relative motion q f . The derivation of Q is given in Appendix A.
In this work, we adopt the von Kármán kinematic assumption for geometric nonlinearities, which is suitable for small strains and moderate rotations [21]. The elastic force f can be directly derived from the differentiation of the strain energy and may be written as a third-order polynomial function of the relative DoFs q f . Equation (7) can be conveniently written in a partitioned form that highlights the coupling between col(R, θ ) and q f as ⎡ where the explicit dependency on q is dropped for clarity. The subscripts R , θ and f indicate the partitions corresponding to R, θ and q f , respectively. In this work, we rigidize the interface by rigidly linking each interface set with a reference virtual node, and expressing all relative interface DoFs q bp at the pth interface set through 6 DoFs q vp ∈ R 6 (3 translational DoFs and 3 rotational DoFs) of the virtual node, as illustrated in Fig. 2. The 6 DoFs q vp represent the relative translations and rotations of the pth virtual node with respect to the reference frame. The rigid body constraints are commonly applied at the interface, when the flexible bodies are connected through rigid joints.
To be specific, we split the relative DoFs q f ∈ R n in the sth body into the sets of relative interface DoFs q b ∈ R n b and relative internal DoFs q i ∈ R n i as where the interface DoFs q b have been further divided into different interface sets from q b 1 ∈ R n b 1 to q bw ∈ R n bw , and w is the number of interface sets. It holds that n b 1 + · · · + n bw = n b . The transformation from DoFs q b of all interface nodes to the DoFs q v of all virtual nodes can be written as where L v ∈ R n b ×nv is the transformation matrix of the entire interface DoFs, and n v is the number of DoFs for all virtual nodes. Matrix L v is formed according to the position of each interface node. The detailed expression of L v is given in Appendix B.
It should be noticed that the FE discretized components, without imposed constraints, allow relative rigid body motion of the flexible bodies with respect to the body reference frame. In the FFR formulation, however, the rigid body motion has already been described by the translation and rotation of the reference frame. To define a unique displacement field, we need to eliminate redundant DoFs, by imposing a set of reference constraints. This is discussed in the next section.

Floating frame definition
We now briefly summarize the nodal-fixed definition [19] and the mean-axis definition [20] of the FFR, together with the embedding technique utilized to impose the constraints introduced by the mean-axis frame definition.

Nodal-fixed frame
The nodal-fixed frame is commonly applied since its definition is straightforward. In this work, the origin of the reference frame is attached to a specified node of the moving body, i.e., no relative translations and rotations of the attached node with respect to the reference frame are allowed. Clearly, the choice of this attached node is not unique. Here, we choose the virtual node of the kth interface set. Mathematically, this is simply done by fixing the 6 relative DoFs q v k of the corresponding virtual node with respect to the reference frame as which corresponds to 6 scalar constraints. For illustration, we show the kinematic description of a crank-shaft system for the nodal-fixed floating frame in Fig. 3. The gray mesh denotes the rigid body motion of each body, defined by the position and orientation of the reference frame. When "standing" at the origin of the reference frame, one observes the relative displacements and rotations of the body as the flexible body is clamped at the origin of the frame.

Mean-axis frame
If an approximated kinematic model for only moderate relative rotations formulation is adopted, as the case of the von Kármán model in the present work, one should try to keep the relative rotations with respect to the reference frame as small as possible. Since the magnitude of the relative displacements and rotations with respect to the nodal-fixed frame largely depends on the choice of the attachment node, the mean-axis floating frame is a more clever choice. Unlike the nodal-fixed frame, the mean-axis frame imposes constraint condition as a function of relative DoFs at each body impartially. The basic idea is to locate the reference frame in such a way that the relative kinetic energy is minimum with respect to an observer stationed on the reference frame. The relative kinetic energy of the j th element in the sth flexible body is defined as [26] According to (3), the relative velocityu j f of an arbitrary point in the j th element is rewritten by statingu Therefore, the relative kinetic energy T r can be expressed as IfṘ andθ are to satisfy the mean-axis condition, the kinetic energy T r should be minimum. As discussed, for instance, in [26], we first rewrite the time derivatives of the Euler parameters as a function of the angular velocity vector ω, i.e., Then, the minimum for T r can be found by posing Equation (16) yields 6 constraint equations to satisfy the mean-axis condition. In [26], the mean-axis constraint equations are further simplified and finally linearized as a function ofq f , which is the time derivative of the relative DoFs with respect to reference frame. The approximated mean-axis condition is expressed as where S ∈ R 6×n is a matrix of constant parameters, usually referred to as inertia integrals. The detailed derivation and linearization from (16) to (17) are given in Appendix C.
In order to express the mean-axis condition in terms of q f , Eq. (17) is integrated in time to obtain Sq f = 0. (18) By applying the mean-axis frame condition (18), the flexible body can no longer undergo rigid body motion with respect to the reference frame. For illustration, the kinematic description of a crank-shaft system for mean-axis frame is also given in Fig. 3. In the mean-axis frame, the relative displacement and rotation (green mesh) of the body exhibit a interfacefree vibration with respect to the reference frame (gray mesh). Generally, the relative displacement and rotation observed from a mean-axis frame will be smaller than their counterparts observed from the nodal-fixed frame, as discussed in [19].

Embedding of mean-axis and interface constraints
While enforcing Eq. (11) for nodal-fixed frame is straightforward, the treatment of Eq. (18) requires more attention, since the constraint conditions are expressed as an explicit form of all relative DoFs q f . By noticing that the mean-axis frame only yields linear constraints, we apply the so-called embedding technique [5] to obtain a minimum number of equations expressed in terms of independent coordinates. As mentioned in [27], the process of imposing all the reference conditions is actually equivalent to static condensation, where the slave variables are eliminated. We can define the generalized DoFs vector q g as where the virtual, boundary and internal DoFs are further split into independent (master) and dependent (slave) sets of coordinates, denoted by the subscript m and s , respectively. Note that all interface DoFs q b are set as slave DoFs q s,b and the DoFs of the virtual nodes are set as master DoFs q m,v since q s,b are determined by q m,v . The rigid interface condition in (10) and mean-axis frame constraint in (18) can be written together as where D ∈ R (n b +6)×(nv +n b +n i ) is the Jacobian matrix of all constraint conditions with respect to the generalized DoFs q g in the mean-axis frame, and the S matrix has been partitioned accordingly. Equation (20) can also be written as where q s = col(q s,b , q s,i ) and q m = col(q m,v , q m,i ), (22) and the matrices D s and D m contain the columns of D corresponding to slave DoFs q s ∈ R n b +6 and master DoFs q m ∈ R nm , respectively, and n m are the number of master DoFs. The generalized DoFs q g can then be written as a function of q m as where H m is the generalized condensation matrix. Finally, according to (23), the relative DoFs q f can be directly written as a function of the master DoFs q m as where H f m contains the rows of H m corresponding to q f . By substituting (24) into (8) and performing a Galerkin projection, we can obtain the EoMs as ⎡ where In (25), the EoMs are expressed in terms of only the master DoFs q m . The˜ refers to quantities relative to a flexible body constrained on the mean-axis frame. The constraint condition in (20) will be identically satisfied. This procedure is referred as the embedding technique in [5]. This embedding technique is not as computationally efficient as using Lagrange multipliers, since the condensed tangent stiffness and mass matrices are characterized by a worse sparsity pattern as compared to the unconstrained, full counterparts. However, it is strongly preferred when applying MOR for the relative DoFs, as any mode extracted from Eq. (25) and used to form the reduction basis would satisfy the mean axis and rigid interface condition exactly.

Flexible multibody equations
Holonomic joint constraints are applied to connect neighboring bodies and/or impose prescribed motion through virtual nodes. For instance, a rigid connection between j th and kth bodies can be imposed as where q 0,v and q m,v are the initial position and relative DoFs of the connecting virtual nodes, respectively, and N v here is equal to the Boolean matrix of selecting the translation DoFs.
It is emphasized here that the joint constraints are not imposed at the internal DoFs q m,i . The constraint Jacobian matrix can thus be written as where the block 0 i reflects the fact that the constraints are not imposed on the internal DoFs q m,i and C m = [ ∂C ∂qm,v 0 i ]. The joint constraints (usually nonlinear) are included with Lagrange multipliers λ as where H is the number of bodies in the FMBS.

Enhanced Rubin substructuring method
The inertial terms of (28) are configuration dependent and therefore need to be updated at every time step during time integration. Likewise, when geometric nonlinearities have to be considered, also the internal force vector f m is configuration dependent. The computational cost of large size nonlinear FMBS using the FFR reference may thus become significant, and MOR is required. The idea is to reduce the size of q m for each body by expressing them as a combination of modes computed after suppressing the rigid body motion of the floating frame, see [1,5,28]. The EoMs for each body in (28) can thus be simplified by suppressing col (R, θ ) as By fixing the rigid body motion of the frame, the quadratic velocity terms Q m and the couplings between rigid body motion and relative displacement vanish. The last term in (29) represents the connecting forces which are imposed only at the virtual node DoFs. The last term can be further rewritten as where B ∈ R nv ×nm is the local-Boolean matrix that selects the interface DoFs q v from q m , and g v ∈ R nv is the interface force imposed at the virtual node. Furthermore, we linearized (29) as where is the linear stiffness matrix after the constraint embedding, and K ff is the linear sparse stiffness matrix of the body. In general, K mm features decreased sparsity patterns as compared to K ff .

Augmented Rubin reduction bases with modal derivatives
In this section, we extend the standard Rubin substructuring method by augmenting the associated reduction basis with MDs to properly consider geometric nonlinear effects. The ROBs are established for each body separately.
The MDs were first proposed in [17,18] for a single structure not undergoing rigid body motion, by differentiating the eigenvalue problem associated to the free vibration with respect to the modal amplitude. The computation of MDs is discussed in detail in [30]. The methods in [30] require an explicit form of the internal nonlinear forces. Alternatively, Slaats [29] proposed the use of finite difference, which allows the computation of the MDs by means of standard nonlinear FE programs, as this method does not require access to the nonlinear forces and Jacobians. Related to this property, we applied the simplified definition of MDs by neglecting these inertia related terms. This technique is usually addressed as the definition without mass consideration, or more recently, static MDs [30]. A more theoretical grounding of the validity of MDs is given in [31].
When the inertial terms are neglected, Eq. (29) becomes We assume here that the external load g m can be written as a stiffness-scaled linear superposition of the free-interface modes (FVMs) as where the FVMs Φ can be obtained by solving the linear eigenvalue problem associated to (31): where ν j is the j th eigenfrequency and φ j is the corresponding FVM. Generally, a truncated set of the first r m FVMs is selected in the reduction basis Φ ∈ R nm×rm based on the frequency range of interest. The reduction will be achieved by letting r m n m . Note that Φ does not contain any rigid body motion since the system has already been fully constrained, see Sect. 3. By substituting (33) into (32), we obtain a static nonlinear problem where static response q m is determined by the modal amplitude ζ ∈ R rm+nv . Instead of finding a solution of (35), we assume that q m is C 2 differentiable with respect to the modal parameter ζ and we expand q m into a Taylor series around equilibrium position as where ζ j is the j th component of the modal parameter ζ . In order to find the derivatives in (36) we differentiate both sides of (35) with respect to ζ , and evaluate them around the equilibrium position as where K −1 mm exists as rigid body motions are suppressed. This procedure distinguishes from the standard Rubin method, where a pseudo-inverse matrix needs to be computed due to the presence of rigid body modes. The matrix Ψ = K −1 mm B T , Ψ ∈ R nm×nv includes the so-called attachment modes (AMs). The AMs represent deformations due to the unit generalized force at one interface DoF and zero to all other interface DoFs. Therefore, the connecting interface force vector g v represents the modal amplitudes of the AMs Ψ . Expression (37) can be compactly written as where X is the standard Rubin ROB. In order to calculate the second-order derivatives of q m with respect to ζ , we differentiate (35) twice to get Evaluating (39) around the equilibrium position gives where the right-hand side of (40) is a null vector, since the external load and interface forces are assumed to be a linear superposition of the modal parameters ζ . The second derivatives of the nonlinear response with respect to the modal amplitudes ∂ 2 qm ∂ζ j ∂ζ k 0 are the MDs, computed from (40) as and it holds that ϑ jk = ϑ kj , see [23]. Note that d 2 fm d(qm) 2 0 X j is the directional derivative of internal tangent stiffness matrix K mm with respect to the modal amplitudes ζ j of mode X j , i.e., Equation (42) certifies that MDs can be computed numerically by the finite difference method, as proposed by [29]. This numerical approach recomputes the configuration dependent stiffness matrix when the structure perturbed from the equilibrium position in the direction of one mode. The perturbation should be small enough for the accuracy of finite difference, but it cannot be too small in order not to incur in round-off errors. In this paper, the von Kármán kinematic model is applied. Since the internal force vector and stiffness matrix can be explicitly expressed as a polynomial function of the DoFs, the MDs can be computed analytically.
Having defined the AMs, FVMs and corresponding MDs, q m can be approximated by the second-order Taylor expansion in (36) as Fig. 4 The nonlinear ROBs for a crank-shaft system. The first body (left) and its corresponding enhanced CB ROB is illustrated using the nodal-fixed frame, while the second body (right) and its corresponding enhanced Rubin ROB is illustrated using the mean-axis frame which constitutes a quadratic manifold for q m in the ζ space. In this work, we will not directly apply (43), as done, for instance, in [30,32]. Instead, the MDs will be included in the ROB as additional independent modes to reproduce geometric nonlinearities. 1 The relative DoFs vector q m with respect to the floating frame is then given by where Θ is the matrix containing the vectors of independent MDs, and ξ is the modal coordinates vector associated to the MDs in Θ. Since the MDs are calculated around the equilibrium position and will not be updated during the time integration, the accuracy of using the modal transformation in (44) will be determined by how far the structure departs from reference position. For illustration, some representative modes of the ROB for a crank-shaft system are shown in Fig. 4. For the first body, the ROB is constructed as done in [1], i.e., a nodal-fixed frame is applied, and the enhanced CB method is applied. The origin of the reference frame is attached to the interface B 1 . Therefore, the fixed-interface modes Φ (1) , MDs Θ (1) are fixed at interface B 1 and B 2 . The compatibility with neighboring bodies is considered by the constraint modes (CMs) Ψ (1) . The mean-axis frame is utilized to describe the motion of the second body. The corresponding FVMs Φ (2) , MDs Θ (2) and AMs Ψ (2) exhibit free motion at the interface sets B 2 and B 3 . It should be noted that all modes are obtained after the meanaxis constraints are included. Therefore, the FVMs in the ROB of the mean-axis frame model contain no rigid body motions. While these low-frequency FVMs of the flat plate models contain bending-dominant vibrations, the corresponding MDs exhibit membrane-dominant vibrations to properly model the nonlinear effects.
In order to assemble the reduced components, Eq. (44) should contain only modal amplitudes. A further transformation is applied to replace the force vector g v with the interface DoFs vector q m,v . The interface partition of (44) is where Ψ v , Φ v and Θ v are the rows of the AMs, FVMs and MDs associated to the interface DoFs, respectively. The interface force vector g v can thus be expressed as By substituting (46) into (44) and recalling (22), the interface DoFs q m,v can be retained in the final coordinates transformation as

Reduced equation of motion
The final reduced EoMs for the sth substructure can be obtained by substituting (47) into (28) and performing a Galerkin projection as ⎡ All nonlinear terms in (48) can be directly expressed as a function of the modal coordinates col(R (s) , θ (s) , γ (s) f ) by a tensorial form as in [1]. The nonlinear force vector f (s) f , whose update is the most time consuming operation during each iteration of the time integration, can be directly expressed as are constant quadratic, cubic and quartic tensors, respectively, and r (s) g is the number of modes in the enhanced ROB V (s) f . The tensors 2 W (s) , 3 W (s) and 4 W (s) can be calculated offline, once the reduction basis of each flexible body is determined.
The system reduced EoMs can be obtained by assembling the contribution from each body and by appending the constraint conditions as It is emphasized here that the constraint equation C = 0 is not projected onto the reduced basis, and therefore, exact interface compatibility has been guaranteed.
In this work, we use the implicit Newmark scheme for the time integration of (50) by setting the integration parameters γ = 1 2 and β = 1 4 . The artificial damping coefficient α is set to zero for all the presented examples. The constraint equation is treated as discussed in [33], where the Lagrange multipliers have been set as additional DoFs. Substantial computational cost reduction can be achieved, in comparison to full analysis, thanks to the reduction in size and the efficient treatment of the nonlinear terms (49). The computational efficiency of applying the implicit Newmark scheme has been discussed in [1], and will not be repeated here.

Numerical examples
In this section, two numerical examples are presented to assess the performance of the proposed reduction method. All the models contain elastic bodies meshed with triangular FE shell elements [34] featuring 3 nodes per element and 6 DoFs per node. The von Kármán kinematic model is adopted. The detailed formulation of the shell element can be found in [34]. In our discussion, the following labeling is utilized to refer to the solutions obtained from different approaches, namely:

Model 1: rotating beam
We consider here the dynamic analysis of a rotating beam, which has been used as a benchmark in many papers dealing with flexible beams and geometric nonlinearities [35][36][37][38]. In all the previous publications, the system shown in Fig. 5 was meshed with planar beam elements. The geometry of the beam and material properties are also illustrated in Fig. 5. An imposed end rotation θ s with respect to OX axis is applied as where T s = 15 s and ω s = 6 rad/s. The nonlinear responses of the tip displacement components W , U and V (see Fig. 5(b)) obtained with different methods are compared. When the global position vector at the tip node r j is obtained from Eq. (1), the tip displacement W , U and V can be calculated for this example as (51) The vector R presents the position of virtual node of the joint, and the superscript t denotes the tip node with nodal position col(0, 10, 0) at the initial state. Notice that since the beam only rotates about the x-axis, it holds that W = 0.
The tip displacement components U and V are shown in Fig. 6. The nonlinear response obtained from the corotational frame of reference (CFR) featuring planar beam element is also included as a reference solution. The CFR [39] is a more general and expensive framework that is able to deal with arbitrary large elastic displacement. This reference solution is denoted as CFR-HFM-NL. From Fig. 6 we can observe a good agreement between the time response of the CFR-HFM-NL, NFR-HFM-NL and MFR-HFM-NL. The good agreement confirms that the adopted von Kármán kinematic assumption is adequate for this numerical example.
In addition, the accuracy of all ROMs is compared to the HFMs. For illustration, the number of modes for different ROMs is listed in Table 1. Since the origin of the NFR is  fixed at the hinge, no interface modes (i.e., constraint modes) are included in the NFR-ECB-NL. The MDs enhanced substructuring method for both nodal-fixed frame and mean-axis frame (i.e., NFR-ECB-NL and MFR-ERubin-NL) can achieve a good approximation of the reference solution. On the contrary, the standard substructuring techniques (without MDs) for both nodal-fixed and mean-axis frame (i.e., NFR-CB-NL and MFR-Rubin-NL) fail to reproduce the full response, even though as many as 50 modes are included in the ROBs.
To further compare the effectiveness of the ROMs between nodal-fixed frame and meanaxis frame, the root-mean square (RMS) error, defined as is plotted in Fig. 7(a) for ω s = 6 rad/s, and the relative error (RE), defined as is plotted in Fig. 7(b) for ω s ranging from 6 to 12 rad/s. Here, r and r are the global position vectors obtained from the HFM and ROMs, respectively. The subscripts X , Y , Z indicate directional components along the OX, OY , OZ axis. As can be observed in Fig. 7(a), although both the ROMs with nodal-fixed frame and mean-axis frame can reproduce the full solution with satisfactory accuracy (RMS errors below 1.5 × 10 −4 ), the MFR-ERubin-NL is more accurate (RMS error below 2 × 10 −5 ) than NFR-ECB-NL, for an ROM of the same size. As can be seen in Fig. 7(b), when the rotational velocity parameter ω s increases from

Model 2: 5 MW/61.5 m wind turbine blade
We consider here a more complex example, namely a 61.5 m long blade of the NREL 5 MW reference wind turbine, which is originally presented in [40]. This model is constructed by assuming constant thickness and homogeneous material. The effective material properties and geometrical parameters are shown in Fig. 8. Rayleigh damping [41] is adopted: a modal damping factor of 2% for the first two modes is used to determine the Rayleigh coefficients. The aerodynamic loads experienced by the blade are calculated using the Blade Element Momentum (BEM) theory, which constitutes a broadly adopted industrial practice for design and analysis of wind turbines. The aerodynamic loads (i.e., normal force F N and tangential force F T per section) are computed as discussed in [42]. For this example, they result in prescribed, time-varying nodal forces at the leading edge of the blade. No aero-elastic interaction is here considered. For illustration, the normal force F N and tangential force F T , calculated as in [42], along the leading edge nodes with length L of 10, 33 and 50 m, respectively, are shown in Fig. 9.
The blade is assumed to rotate around the x-axis with a constant speed Ω = 8 rad/s and a physical time of 100 s is simulated. For the time integration, we use a fixed time step of 0.02 s, with updating of the tangential operator at each iteration within one time step, with  Fig. 10. The computation of the tip displacement has been shown in (51). A clear difference between the linear and nonlinear tip displacement W , U and V can be observed, confirming that the blade vibrates in the nonlinear regime. On the other hand, the relative displacements obtained from MFR-HFM-NL and NFR-HFM-NL are in good agreement.
The different ROBs for the nonlinear analysis, for both nodal-fixed frame and mean-axis frame, are listed in Table 2. The reduced nonlinear responses obtained from different ROMs are shown in Fig. 11. As can be observed, the ROMs with MDs (i.e., MFR-ERubin-NL and NFR-ECB-NL) yield much better approximations than their counterparts without MDs (i.e., MFR-Rubin-NL and NFR-CB-NL), even though many more modes are included in these latter ones. This confirms the efficiency of using MDs to account for the nonlinear effect. Furthermore, the enhanced Rubin basis performs better than the enhanced CB counterpart in the nodal-fixed frame, even though the same number of MDs are included. The better accuracy is clearly highlighted in Fig. 11(d), where the RMS errors of the overall displacement field are compared.
The computational time is compared between the ROMs enriched with MDs (NFR-ECB-NL and MFR-ERubin-NL) and the HFMs (MFR-HFM-NL and NFR-HFM-NL). All simu-  The computational cost of the HFM is denoted by t full while the one of ROMs has been divided into three different components: (i) the construction of the reduction basis, which is regarded as offline cost (t off 1 ); (ii) the calculation of all the higher-order tensors as in (49), also included in the offline cost (t off 2 ), and (iii) the time required for time integration, which constitutes the online cost (t on ). Clearly, the HFM does not bear any offline cost. The computational time is shown in Table 3.  In order to compare the time gains given by the ROMs, a speed-up factor S is then computed as where C off and C on are weight factors for the offline and online stages, respectively. The offline calculation cost is neglected by setting C off = 0, C on = 1. The so obtained speed up factor, denoted as S 1 , is justified when the same ROM is used for many different load cases. Alternatively, one can set equal weights to offline and online costs, i.e., S 2 : C off = 0.5, C on = 0.5. This covers the limit case in which the ROM is used only once. In the context of wind turbines, a considerable number of loading conditions is prescribed by the design standard IEC 61400-1 [43], resulting in a minimum acceptable number of simulations in the order of 1880 [44]. This number grows rapidly when one or more of the parameters associated with the site environmental conditions lies outside the range of IEC reference conditions and may quickly result in up to 3,200,000 simulations [45]. The significance of online cost becomes even more pronounced in the context of digital twin technology, which constitutes the state-of-the-art approach for lifecycle management of wind turbine structures. Such technology consists in combining a virtual system model, i.e., digital twin of the wind turbine, with operational sensor data so as to afford real-time assessment of the structural condition of wind turbines, making thus S 1 a decisive factor in comparison to S 2 . It can be observed that t full and t off 1 of the mean-axis frame are larger than their counterparts of the nodal-fixed frame, since the stiffness and mass matrices M mm and K mm feature a worse sparsity pattern due to the condensation of the mean axis constraints. Therefore, the eigenvalue analysis in Eq. (34), the calculation of MDs in Eq. (41), as well as tangent operator calculation in Newmark time integration are more expensive than their correspondents in the nodal fixed frame. The offline cost t off 2 is similar for the mean-axis frame and nodal-fixed frame ROMs, as t off 2 is mainly determined by the size of ROBs, i.e., the number of modes included in the reduction basis. On the contrary, t on in the mean-axis frame is much smaller than its counterpart in the nodal-fixed frame. This is due to the fact that the MFR-ERubin-NL requires fewer iterations for a given time step because of smaller relative DoFs q f , although NFR-ECB-NL and MFR-ERubin-NL contain the same number of modes in the ROB and their corresponding computational time per iteration is similar.

Conclusions
This paper presents a model-order reduction technique for flexible multibody systems featuring geometrically nonlinear elastic behavior. The overall motion of each body is described with the mean-axis floating frame of reference. The relative displacements of each body are then represented by a basis obtained by enhancing the standard Rubin substructuring basis with modal derivatives computed for both free vibration modes and attachment modes. This allows to accurately capture the geometrically nonlinear elastic behavior of the deformable body. When compared with a previous contribution [1], where a modal derivatives-enhanced Craig-Bampton substructuring method is applied in the nodal-fixed floating frame, the present approach offers a better representation of the nonlinearity at the interface, since the coupling between attachment modes and free-interface modes is considered.
For the reduced-order model, the modal derivatives essentially represent second-order terms of the Taylor expansion of the displacements from the undeformed configuration. As such, it is essential to minimize the relative displacements and rotations with respect to the reference frame. The mean-axis formulation indeed provides generally smaller relative displacements and rotations than their counterpart in the nodal-fixed frame, thus improving the accuracy of the reduced-order model, as shown in the numerical examples.
The method provides significant computational gains when tested on the simulation of a flexible wind turbine blade featuring about 24000 degrees of freedom. The necessary offline cost for the computation of reduction basis is higher than the one proposed obtained in [1], since the projected matrices feature worse sparsity pattern due to the embedding of the mean-axis constraints. However, the online speed-up is better for the chosen test as compared to the one achieved in [1] for a reduced-order model of equal size. This is due to the fewer iterations within a time step required for convergence, as the relative displacements with respect to the reference frame are smaller. This technique is particularly useful when several load conditions need to be simulated so that the offline cost can be amortized. Notice that N j kl is a skew symmetric matrix. Therefore, assumption (67) can be justified on the basis that velocityq j f occur approximately collinear with displacement q j f , and hence the cross product vanishes.
With this assumption, the mean-axis condition can be linearized as and Equations (68) and (69) can be further expressed, with a classic FE assembly, as where S 1 ∈ R 3×n and S 2 ∈ R 3×n are the precomputed inertia integrals, which can be assembled over the elementary component S j 1 ∈ R 3×ne and S j 2 ∈ R 3×ne . The mean-axis condition in (70) implies that the mean-axis condition can be written in terms of onlyq f that have been defined locally w.r.t. the body axis.