Component-level proper orthogonal decomposition for flexible multibody systems

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 361 (2020) 112690 www.elsevier.com/locate/cma

Component-level proper orthogonal decomposition for flexible multibody systems Yunsen Houa , Cheng Liua,b ,∗, Haiyan Hua a

MOE Key Laboratory of Dynamics and Control of Flight Vehicle, School of Aerospace Engineering, Beijing Institute of Technology, Beijing 100081, China b State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, China Received 19 July 2019; received in revised form 30 September 2019; accepted 7 October 2019 Available online 13 January 2020

Abstract Nonlinear finite element methods such as isogeometric analysis and absolute nodal coordinate formulation can be applied to modeling the flexible multibody systems (FMBS) subject to both large rotation and deformation. However, the computational expense is prohibitive for complex systems, especially for a parameterized FMBS. Proper orthogonal decomposition (POD), nonmodal model order reduction (MOR) technique, has been widely applied to nonlinear problems in order to reduce the computational expense of simulating the dynamic responses of FMBSs. This paper presents a component-level POD method, which applies the singular value decomposition process to snapshots of each component of an FMBS and utilizes constraint equations to satisfy the connecting conditions among the components. Furthermore, the paper provides a systematic study of the parameterized MOR (PMOR) for an FMBS, including the reduction of the constraint equations, energy conserving mesh sampling and weighting (ECSW), and the greedy-POD method. The ECSW method aims to further reducing the computational cost of matrix multiplications for evaluating the reduced stiffness matrix. Greedy-POD is applied to get a robust set of reduced order bases with parametric changes in simulation. Finally, three numerical examples validate the feasibility of the componentlevel POD as well as the whole PMOR process for an FMBS, thereby demonstrating that the component-level POD performs better numerical simulations in accuracy computational costs than the classical POD. c 2019 Published by Elsevier B.V. ⃝ Keywords: Component-level greedy-POD; Isogeometric analysis; Flexible multibody system; Parametric model order reduction

1. Introduction The studies on flexible multibody systems (FMBSs) have received great attention from various fields, such as aerospace engineering, vehicle engineering and robotics engineering over the past two decades. Based on an assumption of small deformation, the floating frame of reference formulation (FFRF) [1] defines a floating frame that follows the rigid body motion of each flexible component, such that the small deformation can be separated from the large rotation. Nonlinear finite element methods based on continuum mechanics are employed for systems that couple large rotations and large deformations, such as isogeometric analysis (IGA) [2], absolute nodal coordinate ∗ Corresponding author at: MOE Key Laboratory of Dynamics and Control of Flight Vehicle, School of Aerospace Engineering, Beijing

Institute of Technology, Beijing 100081, China. E-mail addresses: [email protected] (Y. Hou), [email protected] (C. Liu), [email protected] (H. Hu). https://doi.org/10.1016/j.cma.2019.112690 c 2019 Published by Elsevier B.V. 0045-7825/⃝

2

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

formulation (ANCF) [3], and the geometrically exact formulation [4,5]. The dynamic response of a complex system is governed by a set of large-scale nonlinear differential–algebraic equations, and the numerical simulation of the response involves solving large-scale systems of equations, which requires a large amount of computational resources. In order to alleviate the computational burden, several techniques of model order reduction (MOR) have been developed to replace the full-order model (FOM) of a high-dimensional system with a reduced-order model (ROM), but without great losses in accuracy. The component mode synthesis proposed by Craig and Bampton [6] is one of the most widely used MOR approaches in structural dynamics. In the FFRF, the small deformation can be expressed in terms of a few mode shapes obtained numerically or experimentally [7]. Thus, the FFRF can be regarded as a type of MOR technique. Br¨uls et al. [8] proposed the global modal parameterization (GMP), which is a system-level MOR technique for FMBS described via FFRF. This approach was developed further in subsequent studies [9,10]. However, the centrifugal and Coriolis forces still appear in the reduced kinematic equations. In addition, Eberhard et al. [11,12] developed the non-modal technique of MOR based on Krylov-subspaces and singular value decomposition (SVD) in the floating frame description. Furthermore, Eberhard et al. [13,14] developed the time-dependent parametric technique of MOR based on the matrix interpolation for moving load problems and material removal problems in FMBSs. However, these techniques are not suitable for FMBSs with geometric nonlinearity due to large deformation. Wu and Tiso [15] extended the component mode synthesis and applied the modal derivatives [16] to a few planar beam systems undergoing large elastic deflection in the floating frame description. The proposition of the quadratic manifold in further studies [16–19] provided a solid theoretical foundation for this technique which has been widely applied to structural dynamics. Snapshot-based methods, such as proper orthogonal decomposition (POD) [20] and reduced basis method [21– 25], are alternative approaches to MOR. They have been increasingly popular for both linear and nonlinear problems because of their great adaptability. These methods of MOR are among the most important ones for the constrained optimization of partial differential equations [22,24,26–29] in optimal design and optimal control, as well as the solution of inverse problems [26]. POD is a statistical method that aims to finding the low-dimensional subspace, within which a high-dimensional collection of samples is distributed. It is also called principal component analysis (PCA) [30] in machine learning, where it is employed as a method of data compression (an application of PCA). Owing to its simplicity and robustness in dealing with real-world data, POD can be applied to not only the MOR with geometric nonlinearity [31], but also to those with material nonlinearities, such as elasto-plasticity [32,33], hyper-elasticity and visco-elasticity [34], and visco-plasticity [35]. Furthermore, POD has been widely used in the field of computational fluid dynamics (CFD) [36–42]. As a type of statistical method, POD requires a priori simulation of the FOM in order to obtain the snapshot collections extracted from the dynamic response which involves information regarding the location of the lowdimensional subspace. The orthonormal bases of the subspace, which are also called the reduced order bases (ROBs), can be obtained by applying SVD to the snapshot collection. Thus, the MOR technique is only suitable for repeated simulations, e.g., the dynamic simulations of a system with various parameters. Therefore, the method of parametric model order reduction (PMOR) [31,37,43,44] was proposed. It aims to finding the set of ROBs appropriate for the dynamic response of a system in terms of a specific parameter without computing the FOM, and the ROM enables one to obtain the response directly and rapidly. Amsallem and Farhat [37] proposed an interpolation method based on Grassmann manifold interpolation. It assumes that the relationship between the parameters and the ROB sets follows a continuous mapping, such that the ROB set for a new parameter can be obtained by the special Grassmann manifold interpolation with the ROB sets for pre-computed parameters. A more suitable approach for PMOR is the greedy-POD method [22,24,25,43,45] which assumes that the dynamic responses associated with all parameters fall into a single low-dimensional subspace. This approach iteratively determines the configuration that leads to the maximal error under the current ROBs, solves the FOM for that configuration and updates the ROBs with the corresponding FOM response. Paul-Dubois-Taine and Amsallem [43] applied Kriging interpolation [46] to find the maximal error in each iteration. This technique essentially reduces the times required for computing the ROMs under the current ROBs, thereby eliminating most of the time needed for the training process. However, calculating the exact error between the ROMs and FOMs requires the dynamic responses of the FOMs which are not available in practice. Fortunately, the norm of the residual of the FOM is often highly correlated with the exact error, and thus it is generally employed as a convenient indicator of the exact error to reflect the consistency between the FOM and ROM [42,43,47,48] because of its much lower computational cost and rigorous error bound for the exact error.

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

3

POD can help to essentially reduce the dimensions of a large-scale system but it fails to allow the efficient simulation of nonlinear models because the nonlinear terms in the original system need to be evaluated and projected onto the subspace in each time step. Unfortunately, this process scales to the size of the FOM, thereby making it high costly in terms of the computational resources required. This problem has been addressed by developing hyper reduction [49] techniques, such as gappy POD [41], missing point estimation [38], the Gauss–Newton with approximated tensors method [42], and discrete empirical interpolation method [34,35,39,50], where the latter is now the most popular method for nonlinear first-order parabolic problems in CFD. However, it cannot be applied to nonlinear second-order hyperbolic problems in structural dynamics because of its numerical instability [49] and certainly not to FMBS. Instead, Farhat et al. [49,51,52] developed the energy conserving mesh sampling and weighting (ECSW) method which is suitable for structural dynamics with no numerical instability, because it preserves the Lagrangian structure associated with Hamilton’s principle of a system. The cubature scheme first proposed by Kim et al. [53] allowed the development of this method. POD and hyper reduction techniques have facilitated much progress in structural dynamics and CFD, but their application to FMBS is still in its infancy. Stadlmayr et al. [54,55] investigated the reduction of the constraint equations and their redundancy. In the previous study of Luo [31], an attempt was made to apply the classical POD to FMBS (described by ANCF) with coupled large rotations and deformations, which could be regarded as a systemlevel MOR technique. In this study, it was found that the numerical performance can be improved (better accuracy and lower computational costs) if one applies the SVD process to each component of an FMBS rather than to the whole system, and utilizes constraint equations to satisfy the connected conditions among the components. This technique can help to reduce the number of scalar multiplications when evaluating the reduced stiffness matrix in each Newton–Raphson iteration, thereby further improving the computational efficiency. It was inspired by previous studies [32,33,56] applying the domain decomposition technique to POD, where it was referred to as the DD-POD method in [33]. In this study, this technique is interpreted as component-level POD whereas the classical POD in the previous study of Luo [31] is interpreted as system-level POD. This study aims to demonstrating the benefits of the component-level POD in terms of its computational cost and accuracy, as well as presenting a systematic study of the PMOR for an FMBS, including the reduction of the constraint equations, application of the ECSW technique, and the greedy-POD method based on Kriging interpolation. The remainder of this paper is organized as follows. In Section 2, a brief introduction is given to model an FMBS subject to both large rotation and large deformation. Then, the component-level POD process is proposed in Section 3 for an FMBS, including the reduction of the constraint equations and the hyper reduction technique (ECSW). In Section 4, the greedy algorithm is discussed for parametric POD. In Section 5, three numerical examples are provided to validate the effectiveness of the proposed framework for the PMOR process with a parameterized FMBS. Finally, some concluding remarks are made in Section 6. 2. Brief introduction to modeling an FMBS with large deformation The kinematic equation for an FMBS can be expressed in a scalar form as the equilibrium of virtual work ˙ δq) + m (q, ¨ δq) + l (q, q; ˙ δq) + λ · Φq (q, t; δq) = 0 a (q, q; , δλ · Φ (q, t) = 0

(1)

where q, δq, and δλ are the generalized coordinates, the variation in the generalized coordinates, and the variation ˙ δq), m (q, ¨ δq), and l (q, q; ˙ δq) are respectively the virtual in the Lagrange multipliers, respectively. The a (q, q; work of the internal force, inertial force, and external force. The λ, Φ(q, t), and Φq (q, t; δq) denote separately the Lagrange multipliers, constraint equations, and direction derivative of the constraint equations w.r.t. q. In this work, the IGA and ANCF are chosen to describe configurations of an FMBS, making the kinematic equation free from any gyroscopic or centrifugal term. Since δq and δλ are arbitrary in each equation, the vector form of Eq. (1) can be denoted as ˙ + Mq¨ + Q (q, q) ˙ + ΦTq λ = 0 F (q, q) , (2) Φ (q, t) = 0 ˙ and Q(q, q) ˙ are the internal force and external force, respectively, and M is the total mass matrix. By where F(q, q) using the generalized-α algorithm [57,58] and the Newton–Raphson iteration, the discretization and linearization

4

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

of Eq. (2) leads to the linear algebraic equations as follows [ ][ ] [ ] ˆ + γˆ C + K) ΦTq ∆q c(βM R =− Φ , ∆λ Φq 0 R = c(F + Mq¨ + Q) + ΦTq λ

(3)

ˆ γˆ are the parameters determined by the spectral radius ρ∞ of the generalized-α algorithm [57,58], where λ = cλ, β, R is the residual of the equilibrium equation, K and C are separately the Jacobian of the residual R w.r.t. q (the tangential stiffness matrix) and q˙ (the damping matrix), and c is a scaling factor used to reduce the condition number for the total matrix. In practice, c is taken as h2 β, where h is the time step while β is the other parameter determined by ρ∞ . For the sake of brevity, it is assumed that F is independent of q˙ and the Jacobians of Q w.r.t q˙ and q are omitted in this work. If they do have a significant impact, the extension is direct. This study focuses on the FMBSs involving slender beam (cable) components and thin plate (shell) components. Isogeometric analysis (IGA) [2] has proved a kind of powerful tools in CAD geometry for modeling those flexible components. So did the gradient-deficient finite elements of ANCF proposed by Shabana [59,60]. IGA utilizes the basis function generated from B-splines and non-uniform rational B-splines (NURBS) for both exact geometric description and finite approximation. As shown in [2], the univariate B-spline basis functions of p-degree denoted by Ni, p (ξ ) corresponding to a non-decreasing knot vector Ξ = {ξ1 , ξ2 , . . . , ξn+ p+1 } are defined as the following recursive formulae { 1 if ξi ≤ ξ ≤ ξi+1 Ni,0 (ξ ) = for p = 0, (4) 0 otherwise and Ni, p (ξ ) =

ξi+ p+1 − ξ ξ − ξi Ni, p−1 (ξ ) + Ni+1, p−1 (ξ ) ξi+ p − ξi ξi+ p+1 − ξi+1

The NURBS basis is defined as Ni, p (ξ ) wi p R i = ∑n , k=1 Nk, p (ξ ) wk

for p ≥ 1.

(5)

(6)

where wi is the weight corresponding to the ith B-spline basis function. The two-dimensional NURBS basis can be defined by the tensor product of the one-dimensional ones. For the parameter ξ and η, and the associated polynomial degrees p and q Ni, p (ξ ) M j,q (η) wi, j ∑m , k=1 l=1 Nk, p (ξ ) Ml,q (η) wk,l

p,q Ri, j (ξ, η) = ∑n

(7)

where M j,q (η) is the B-spline basis defined in the η direction. The continuity of the basis function µ (i.e., Cµ continuous) is determined by the degree of the function p and the maximal multiplicity of the interior knots k, i.e., µ = p − k. Further details are referred to [2]. For modeling slender beam (cable) components and thin plate (shell) components, ANCF is quite similar to IGA since they both are rotation-free methods of nonlinear finite elements [61]. In fact, there exists a linear transformation, which defines the relationship between the B-spline surface (curve) control points and the ANCF generalized coordinates including position and gradient vectors [62]. The shape functions of the beam elements and shell elements of ANCF are the 3rd Hermite interpolations, as presented in previous studies [60] and [59,63], respectively. It should be noted that the application of the componentlevel POD method discussed in Section 3.4 is not limited to the finite elements of IGA (or ANCF), but also suitable for other nonlinear finite elements, such as those of geometrically exact formulation [4,5]. The derivations of the mass matrix, internal force, and the corresponding Jacobian for a slender beam (cable) element or a thin plate (shell) element were presented in the previous study of Liu [64,65], and thus the details are omitted for brevity. It should be noted that the derived mass matrix is constant, which is a great advantage for the POD method because the reduced mass matrix in the ROM is also constant without gyroscopic and centrifugal terms, and free from updating during simulations. In fact, when the system of concern is a large complex one, the set of linear algebraic equations in Eq. (3) is extremely large. Solving these linear equations is the main computational burden in the entire simulation process. Furthermore, it is common to conduct many reruns for a large complex system to be parameterized, e.g., in the

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

5

Fig. 1. FMBS comprising m components (Body I, II, III comprise Component m-2) where the joints connecting all of the components are shown. Each component has its own ROBs, modal coordinates, and constraint equations. Other constraint equations of the system are referred to as the joints connecting components.

design optimizations [27] and uncertainty quantifications [47]. In order to reduce the dimensionality of the linear equations to accelerate the simulation process, it is natural to employ the MOR techniques, i.e., the POD (parametric POD) approach, as introduced in the following section. 3. Component-level POD for an FMBS This section presents the framework of component-level POD for an FMBS. Sections 3.1 and 3.2 explain the basic ideas involved, including how to obtain a set of ROBs for an FMBS and how to produce the kinematic equation of a ROM, i.e., the reduced kinematic equation for an FMBS. Then, Section 3.3 gives the reduction of the constraint equations for an FMBS, which does not appear in structural dynamics. The last subsection focuses on the computational cost of the POD method with ECSW introduced as a hyper reduction technique. Before describing the component-level POD framework for an FMBS, it is necessary to note that the word “component” can refer to a single body of an FMBS but also to a substructure of an FMBS. As shown in Fig. 1, an FMBS is considered to comprise m components with various joints connecting all of the components, where Body I, II, III comprise Component m−2. 3.1. ROBs for an FMBS The procedure employed to obtain the snapshots of the FMBS is analogous to the framework used in structural dynamics [20]. The original snapshot collection matrix can be obtained as [ ] S0 = qt1 qt2 · · · qtn ∈ R N ×n , (8) where qt1 , qt2 , . . . , qtn ∈ R N are the vectors of the generalized coordinates of the FOM at n moments t1 , t2 , . . . , tn , and N is the total number of degrees of freedom (DOFs) for the FOM. In practice, it is usual to select the n moments evenly from the total simulation time and to modify the original snapshot matrix to obtain a centralized snapshot matrix denoted by ] [ ] [ S = qt1 − q qt2 − q · · · qtn − q = q˜ 1 q˜ 2 · · · q˜ n ∈ R N ×n n ∑ , (9) q= qti /n i=1

where q ∈ R N is the vector of the generalized, averaged coordinates.

6

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

Based on a specific decomposition, the snapshot matrix S can be rewritten as ⎡ ⎤ ⎡ ⎤ S1 q˜ 11 q˜ 12 · · · q˜ 1n ⎥ ⎢ ⎥ ⎢ ⎢ S2 ⎥ ⎢ q˜ 2 q˜ 2 · · · q˜ 2 ⎥ ⎢ ⎥ ⎢ 1 n· ⎥ 2 S=⎢ . ⎥=⎢ . , .. . . .. ⎥ ⎢ . ⎥ ⎢ . ⎥ . . . ⎦ ⎣ . ⎦ ⎣ . Sm

q˜ m 1

q˜ m 2

···

(10)

q˜ m n b

where m is the number of components in the FMBS, q˜ ib ∈ R N is the ith (i = 1, 2, . . . , n) snapshot of the bth b (b = 1, 2, . . . , m) component, Sb ∈ R N ×n is the snapshot matrix of the bth component, N b is the number of DOFs for the bth component. The superscripts denote the values of the variables for each component. In addition, the average vector q can be decomposed as follows [ ]T q = q1T q2T · · · qmT (11) A set of ROBs for the bth component can be defined in terms of the following optimization problem [20]: find [ ] vectors Ub = ub1 ub2 · · · urbb , r b ≤ n that satisfy ⎛ b ⎞ n r ∑ ∑   ⟨ ⟩ 2 ⎝ Ub = argmax q˜ ib , ubj ⎠ , with u j 2 = 1, (12) i=1

j=1

where ⟨·, ·⟩ denotes the inner product of two vectors, and ∥·∥ indicates the 2-norm of a vector. In this manner, all of the vectors in the snapshot matrix of the bth component Sb can be maximally represented by linear combinations of the ROBs. The solution to Eq. (12) can be obtained by SVD to Sb as ( )diag bT b b Sb = Ub0 Λb0 V0 , U0 ∈ R N ×n , Λb0 ∈ Rn , Vb0 ∈ Rn×n , (13) ( b )diag [ ]T where Λ0 is a diagonal matrix reformulated via the singular value vector Λb0 = λb1 λb2 · · · λbn with descending sorting, and Ub0 and Vb0 are two orthonormal matrices where their column vectors are the left b vectors and [singular ] right ones, respectively. Next, the singular value vector is assumed to be truncated as Λ = b b b T λ1 λ2 · · · λr b with an error ε defined as / n n ∑ ( b )2 ∑ ( b )2 ε= λi λi . (14) i=r b +1

i=1

In this study, the number for ROBs of a component needs to be determined from Eq. (14) with a prescribed error. The truncation error ε is assumed to be extremely small to guarantee the fidelity of the ROM simulations. Correspondingly, Ub0 and Vb0 also need to be truncated as follows ) ( ) ( (15) Ub = Ub0 :, 1 : r b , Vb = Vb0 :, 1 : r b Hence, the centralized snapshot matrix can be approximated as ( )diag bT b b b b b b Sb ≈ Ub Λb V , U ∈ R N ×r , Λb ∈ Rr , Vb ∈ Rr ×r .

(16)

In this manner, the low-dimensional subspace (or affine subspace) of the dynamic response of the bth component can be spanned via the vectors in Ub and qb , i.e. { } Vb = span Ub ⊕ qb . (17) Therefore, in the ROM, any generalized coordinates of the bth component qb can be assumed to fall into the subspace, i.e., ∀qb ∈ Vb . 3.2. Reduced kinematic equation for an FMBS Without loss of generality, this subsection (and Section 3.3) deals with an FMBS comprising two components to obtain the reduced kinematic equation. If qb (b = 1, 2) is restricted to the subspace Vb , then it can be represented

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

7

by a linear combination of ub1 , ub2 , . . . , urbb plus qb b

qb = q˜ b + qb =

r ∑

ξib uib + qb = Ub ξb + qb .

(18)

i=1

Thus, δqb , the variation in qb , can be derived as ( ) δqb = δ Ub ξb + qb = Ub δξb .

(19)

All of the constraint equations are assumed to belong to two classes: a class comprising the internal constraint equations for Component 1 and Component 2 denoted by Φ1 and Φ2 , respectively, and the other class comprising the interacting constraint equations between Component 1 and Component 2 denoted by Φ12 , as shown in Fig. 1. Now, the reduced kinematic equation for the FMBS can be denoted as ⎧ ( ( ) ) ) ( ( ) ⎪ a 1 U1 ξ1 + q1 ; U1 δξ1 + m 1 U1 ξ¨ 1 , U1 δξ1 − l 1 U1 δξ1 + λ1 · Φ1q1 U1 ξ1 + q1 , t; U1 δξ1 ⎪ ⎪ ⎪ (a) ⎪ ) ( 1 1 ⎪ ⎪ 1 2 12 2 2 1 1 12 ⎪ q , U ξ + q , t; U δξ = 0 +λ · Φ U ξ + ⎪ 1 ⎪ q ⎪ ( ) ⎪ ( ) ) ( ) ( ⎪ ⎪a 2 U2 ξ2 + q2 ; U2 δξ2 + m 2 U2 ξ¨ 2 , U2 δξ2 − l 2 U2 δξ2 + λ2 · Φ2 U2 ξ2 + q2 , t; U2 δξ2 ⎪ ⎨ q2 (b) ) ( 1 1 , (20) 1 2 12 12 2 2 2 2 +λ · Φ q , U ξ + q , t; U δξ = 0 U ξ + ⎪ ⎪ q2 ⎪ ⎪ ( ) ⎪ ⎪ δλ1 · Φ1 U1 ξ1 + q1 , t = 0 (c) ⎪ ⎪ ⎪ ) ( 2 2 ⎪ ⎪ 2 2 2 ⎪ δλ · Φ U ξ + q , t = 0 (d) ⎪ ⎪ ⎪ ( ) ⎩ 12 1 2 δλ · Φ12 U1 ξ1 + q , U2 ξ2 + q , t = 0 (e) where λ1 , λ2 , and λ12 are the Lagrange multipliers corresponding to Φ1 , Φ2 , and Φ12 , respectively. Eqs. (20)(a) and (b) are the equilibrium of the virtual work for Component 1 and Component 2, respectively. Here, the last terms in Eqs. (20)(a) and (b) denote the action of the interacting constraint force between the two components. Since δξ1 , δξ2 , δλ1 , δλ2 , and δλ12 are arbitrary in each equation, the vector forms of Eqs. (20) are as follows ( ( ) ) U1T R1 = U1T F1 U1 ξ1 + q1 + M1 U1 ξ¨ 1 − Q1 + Φ1q1 T λ1 + Φ12T λ12 = 0 q1 ) ( ( ) T 12 . (21) λ =0 U2T R2 = U2T F2 U2 ξ2 + q2 + M2 U2 ξ¨ 2 − Q2 + Φ2q2 T λ2 + Φ12 q2 ) ( ) ( ) ( Φ1 U1 ξ1 + q1 , t = 0, Φ2 U2 ξ2 + q2 , t = 0, Φ12 U1 ξ1 + q1 , U2 ξ2 + q2 , t = 0 The linearization of Eqs. (21) leads to the linear algebraic equations denoted as ⎡ ⎤ 1 1 1T 1 T 1T 12T ˆ 0 U Φq1 0 U Φq1 ⎥ ⎡ 1 ⎤ ⎢c(βMr + Kr ) ⎡ 1T 1 ⎤ ⎢ ⎥ ∆ξ U R 2 2 2T 2 T 2T 12 T ⎥ ⎢ 2 ˆ 0 c(βMr + Kr ) 0 U Φq2 U Φq2 ⎥ ⎢ ∆ξ ⎥ ⎢ ⎢U2T R2 ⎥ ⎥ ⎢ ⎢ ⎢ ⎥ ⎥⎢ 1 ⎥ ⎢ Φ1 U1 ⎥ ⎢ ∆λ ⎥ = − ⎢ Φ1 ⎥ , 0 0 0 0 ⎢ ⎢ ⎥ ⎥⎢ 2 ⎥ q1 ⎢ ⎥ ⎣ ∆λ ⎦ ⎣ Φ2 ⎦ ⎢ ⎥ 2 2 0 Φq2 U 0 0 0 ⎢ ⎥ 12 Φ12 ⎣ ⎦ ∆λ Φ12 U1 Φ12 U2 0 0 0 q1 q2 b

b

b

b

(22)

where Mrb = UbT Mb Ub ∈ Rr ×r and Krb = UbT Kb Ub ∈ Rr ×r are the reduced mass matrix and reduced stiffness matrix for the bth component, respectively. In particular, if a system has a single body only, or multiple bodies but it is regarded as a single component artificially, then the reduced kinematic equation for the system will only consist of Eqs. (20)(a) and (c) without λ12 · Φ12 q1 , which is exactly the same as the system-level POD described in the previous study of Luo [31]. In general, the number of ROBs for all the components is much smaller than the DOFs for an FMBS, i.e., r1 + r2 ≪ N . Thus, the scale of the reduced tangent matrix in Eq. (22) is much smaller than the tangent matrix in Eq. (3). Solving the linear algebraic equations in Eq. (22) will be much more efficient, which is why the computational time will be decreased greatly by computing the ROM for an FMBS. However, matrix multiplication is an extremely

8

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

time-consuming process when evaluating the reduced stiffness matrix. Details of how to handle this problem are discussed in Section 3.4. 3.3. Reduction of the constraint equations Unlike the problems of structural dynamics, this study must deal with the reduction of the constraint equations for an FMBS. According to Eq. (22), the Jacobian of the reduced constraint equations reads ⎤ ⎡ 1 1 U 0 Φ [ 1 ] ⎢ q1 ⎥ U ⎥ ⎢ Φq U = Φq 2 = ⎢ 0 (23) Φ2q2 U2 ⎥ . U ⎦ ⎣ Φ12 U1 Φ12 U2 q1 q2 The so-called “reduction of the constraint equations” is simply a right multiplication to the Jacobian of the constraint equations with the ROBs. However, it should be noted that the redundancy in Φq U is different from the redundancy in Φq . Now, Section 3.3.1 explains how to remove the redundancy in Φq U. 3.3.1. Removing the redundancy in Φq U Previous studies [31,54] suggested that the entries in each row of Φq U should be checked at the beginning of the simulation. If an all-zero row is detected, it is necessary to remove the corresponding constraint equation and the Lagrange multiplier. A more generalized method [55] applies the other SVD process to Φq U to obtain the linearly independent rows of Φq U. However, there is no clear explanation why the all-zero rows appear in Φq U. For numerical implementation, the linear independent rows in Φq U can be obtained via the Gauss elimination method. To analyze the redundancy in Φq U, all of the constraint equations are assigned again to three classes, i.e., the time-independent linear constraint equations such as the case of a spherical joint, the time-dependent constraint equations such as control inputs, and the nonlinear constraint equations, such as a pair of revolution joints and a pair of gear joints. As for the time-independent linear constraint equations in Φb (b = 1, 2) denoted by Φlb ∈ Rml , which can be expressed as ( ) Φlb qb ≡ Ab qb + cb = 0, (24) b

where Ab ∈ Rml ×N is a constant Jacobian of the time-independent linear constraint equations, cb ∈ Rml is a constant vector, and ml is the number of time-independent linear constraint equations. It should be noted that qb is the average of all the snapshots of Component b, it must satisfy Ab qb + cb = 0 ⇒ cb = −Ab qb .

(25)

Therefore, Eq. (24) can be rewritten as ( ) ( ) Φlb qb ≡ Ab qb − qb = Ab q˜ b = 0.

(26) b

b

b

Hence, the constant linear constraint equations can be regarded as a restriction of q˜ = q − q to the kernel subspace of the linear transformation Ab . In addition, from the procedure for obtaining the ROBs, one has { } ( ) { } (27) span Ub ⊆ span q˜ b1 , q˜ b2 , . . . , q˜ bn ⊆ ker Ab . ) ( ) ( ) ( b Therefore, ∀qb = Ub ξb + qb , ξb ∈ Rr , satisfies Φlb qb = Φlb Ub ξb + qb ≡ 0. The gradient of Φlb qb w.r.t. ξ reads ( ) ( b) ∂Φlb Ub ξb + qb ∂Φlb ∂qb = = Φl qb Ub = Ab Ub ≡ 0, (28) ∂ξb ∂qb ∂ξb which is why the all-zero rows appear in Φq U. Thus, it is easy to conclude that all of the constant linear constraint equations will disappear after reduction. Now, if all of the time-independent linear constraint equations in Φ12 are taken into account, one has ( ) ( ) A1 U1 ξ1 + q1 + A2 U2 ξ2 + q2 + c12 = D1 ξ1 + D2 ξ2 + d12 = 0, (29)

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690 12

9

12

b

where Db = Ab Ub ∈ Rm ×r (b = 1, 2) and d12 = A1 q1 + A2 q2 + c12 ∈ Rm are constant matrices and a constant vector, respectively. m12 is the number of time-independent linear constraint equations in Φ12 . The number of ROBs for one component must be greater than m12 , i.e., r1 > m12 or r2 > m12 ; otherwise, there will be no solution to Eq. (29). Of course, it is possible to obtain a solution via the least square fitting, but the solution may violate a constraint. An alternative approach for dealing with the time-independent linear constraint equations in Φ12 involves removing them and reducing some modal coordinates of ξ1 or ξ2 if r1 > m12 or r2 > m12 holds, respectively. 1 12 Without loss of generality, it is assumed that r1 > m12 and two Boolean matrices comprising Rd ∈ Rr ×m and 1 1 12 R f ∈ Rr ×(r −m ) are introduced to pick up some columns from D1 and also to pick up some entries from ξ1 that satisfy ( )( ) ( )( ) D1 ξ1 = D1 Rd RTd ξ1 + D1 R f RTf ξ1 = D1d ξ1d + D1f ξ1f , (30) 12

12

12

1

12

where D1d = D1 Rd ∈ Rm ×m and D1f = D1 R f ∈ Rm ×(r −m ) are formed by the columns from D1 , and 1 12 12 ξ1d = RTd ξ1 ∈ Rm and ξ1f = RTf ξ1 ∈ R(r −m ) are formed from the entries of ξ1 . After applying Eqs. (30), (29) can be rewritten as D1d ξ1d + D1f ξ1f + D2 ξ2 + d12 = 0.

(31)

In order to reduce the modal coordinates ξ1d , it is necessary to make D1d invertible such that ( )−1 1 1 ( 1 )−1 2 2 ( 1 )−1 12 ξ1d = − D1d D f ξ f − Dd D ξ − Dd d .

(32)

Thus, ξ1d can be removed from ξ1 to reduce the DOF of the ROM by [ ] [ ( )−1 1 ( )−1 2 ] [ 1 ] [ ( )−1 12 ] ξf −Rd D1d D f + R f −Rd D1d D ξ1 −Rd D1d d = + 2 2 , ξ 0 0 Ir 2 ξ

(33)

= Ty + d where ( )−1 1 ( )−1 2 ] −Rd D1d D f + R f −Rd D1d D T= 0 Ir 2 [ ( ( ) [ ( )T ( ) ] )T T , d = − Rd D1d −1 d12 y = ξ1f ξ2 [

01×r 2

]T .

After applying Eqs. (33) to (22), the linearized reduced kinematic equation can be rewritten as [ ][ ] [ ] ˆ r t + Kr t ) (Φr t )Ty ∆y R c(βM = − rt , ∆λ Φ 0 (Φr t )y rt rt

(34)

(35)

where T

Mr t = T

[

Mr1 0

] [ 1 0 T Kr T, Kr t = T Mr2 0

The Jacobian of Φr t w.r.t. y reads ⎡( ) ⎤ Φ1n q1 U1 0 ⎢ ⎥ ( 2) 0 Φn q2 U2 ⎥ (Φr t )y = ⎢ ⎣( ⎦ T, ) ( 12 ) 12 1 2 Φn q1 U Φn q2 U

] [ 1T 1 ] 0 T U R T, Rr t = T . Kr2 U2T R2

(36)

(37)

where the subscript n designates the remaining linear independent rows of the Jacobians of the constraint equations. In Eq. (35), Φr t and ∆λr t are ⎡ ⎤ ⎡ ⎤ ∆λ1n Φ1n ⎢ ⎥ ⎢ ⎥ ∆λr t = ⎣ ∆λ2n ⎦ , and Φr t = ⎣ Φ2n ⎦ , (38) 12 ∆λ12 Φ n n respectively, where •n designates the entries in ∆λ or Φ corresponding to the independent rows of Jacobians of the constraint equations. It should be noted that an ill-conditioned D1d (chosen from D1 ) makes it difficult for the

10

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

Newton–Raphson iteration method to converge in numerical implementations. The Gauss elimination method can be again employed to select the smallest condition number for D1d . However, the analogous characteristics are not suitable for the other two types of constraint equations, i.e., the time-dependent constraint equations and the nonlinear constraint equations. In most cases, the Jacobians for these two types of constraint equations discussed above are linear-independent, i.e., irredundant, after reduction. Nevertheless, some scenarios exist where redundancy appears from the two types of constraint equations discussed above. In summary, in order to remove the redundancy in Φq U, it is possible to firstly remove all of the timeindependent linear constraint equations in Φ1 and Φ2 , as well as in Φ12 and removing ξ1d from ξ1 in passing, and afterwards to apply the Gauss elimination method to the Jacobians for the remaining two types of constraint equations to remove the redundancy. 3.3.2. “Gappy” norm of the residual In the greedy-POD framework for PMOR described in Section 4, the information about the norm of the residual in Eq. (3) is essential for greedy iteration because it plays a role as an error indicator in the iteration process. However, some constraint equations in the original model are removed after reduction. Thus, the corresponding Lagrange multipliers are also removed. As a result, the constraint forces that have been removed and the residual of the FOM can no longer be obtained accurately in the ROM. To this end, a “gappy” norm is defined to quantify the size of the residual, which is based on the study of gappy POD (a hyper reduction technique) [41]. It defines a mask vector g to describe for a particular vector where entry is relevant (or irrelevant) to the redundant constraint equations after reduction. It is defined as: gi = 0 (or 1) if the ith entry is relevant (or irrelevant) to the redundant constraint equations. Then, the “gappy” inner product is defined as ⟨u, v⟩g =

N ∑

gi u i vi ,

(39)

i=1

and the induced norm is ∥v∥2g = ⟨v, v⟩g . The constraint forces are circumvented by computing the “gappy” norm of the residual for the FOM. In fact, the “gappy” norm is also a necessity for structural dynamics because the constraint forces corresponding to the boundary condition in the ROM of a system also cannot be obtained. 3.4. Hyper reduction technique and component-level POD This subsection focuses on the computational cost of matrix multiplication in Eq. (22) during the simulation of a ROM. As for finite element method of NURBS (or ANCF), a major advantage is that the mass matrix is a constant matrix in the kinematic equation. Thus, the reduced mass matrix is also a constant matrix and it can be evaluated in a preliminary process. However, due to the nonlinearity of F with respect to q in Eq. (2), updating the reduced stiffness matrix of a ROM accounts for most of the time, including evaluating the stiffness matrix of the associated FOM and matrix multiplication in Eq. (22), but especially the latter. In addition, the time spent on calculating Φq U can be ignored compared with matrix multiplication when evaluating the reduced stiffness matrix. It is better to note the sparsity of the original stiffness matrix of each component Kb , where w is used to denote the maximal number of non-zero elements in each row of Kb , i.e., the bandwidth. The computational complexity involved when calculating the term Kb Ub is of the order O(N b wrb ) and that for the term UbT (Kb Ub ) is of the order O(N b (rb )2 ). Therefore, the computational complexity of matrix multiplication is proportional to the number of DOFs for each component and to the SQUARE of the number of ROBs required for each component. Thus, improving the computational efficiency depends on reducing the number of DOFs N b or reducing the number of ROBs rb when evaluating Krb . 3.4.1. ECSW ECSW is a type of hyper reduction technique for finding a reduced mesh comprising some of the finite elements of the original mesh for each component. The reduced stiffness matrix Krb can be evaluated based only on the reduced mesh, which helps to reduce the number of DOFs N b . A brief introduction to the ECSW method is provided in the following.

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

11

In this subsection, the superscript b used to denote the values of the variables for each component are omitted for the sake of brevity. In fact, the virtual work done by the internal force can be regarded as the sum of all the virtual work done by every single element e of each component. Hence, the virtual work reads δWint = δqT F (q) =

K ∑

( ) δqeT Fe qe ,

(40)

e=1

where δWint denotes the total virtual work, K denotes the number of finite elements in each component, Fe and qe are the internal force and generalized coordinate of the eth element, respectively, and the superscript e denotes the e values of the variables for each finite element (instead of each component). The Boolean matrix Le ∈ R N ×N maps the global coordinates (forces) to the element coordinates (forces), then one has qe = Le (Uξ + q) = (Le U)ξ + Le q = Ue ξ + qe , δqe = Ue δξ, Fe = Le F

(41)

where Ue = Le U. After applying Eqs. (41), (40) can be rewritten as δWint =

K ∑

K ∑ ( ) ) ( δξT UeT Fe Ue ξ + qe = δξT UeT Fe Ue ξ + qe .

e=1

(42)

e=1

ECSW aims to finding an approximation δ W˜ int for δWint via a sum over some chosen finite elements instead of all finite elements. All of the chosen finite elements comprise the reduce mesh [51]. To compensate for the missing contributions of the unselected elements, the virtual work of every selected element is weighted by a positive weighting factor ωe > 0. The approximation for δWint then reads δWint

≈ δ W˜ int = δξT

Kr ∑

) ( ˜ ωe˜ UeT Fe˜ Ue˜ ξ + qe˜ ,

(43)

e=1 ˜

where Kr denotes the number of finite elements in the reduced mesh, i.e., the number of all the chosen finite elements, and e˜ is the index of a finite element in the reduced mesh. By using Eq. (43), the approximation for the reduced internal force can be written as Kr ( ) ∑ ˜ Fr ≈ ωe˜ UeT Fe˜ Ue˜ ξ + qe˜ . (44) e=1 ˜

The approximation for the reduced stiffness matrix can be written as Kr ≈

Kr ∑

) ( ˜ ωe˜ UeT Ke˜ Ue˜ ξ + qe˜ Ue˜ .

(45)

e=1 ˜

It is natural to require that the approximations of Eqs. (44) and (45) give the desired level of accuracy, i.e., an error within a given tolerance η, with the smallest possible number of finite elements in the reduced mesh. Finding the reduced mesh and the corresponding positive weight factors ω ∈ R Kr is an NP-hard problem, but an approximate solution can be obtained with a sparse non-negative least square (sNNLS) solver, as described in [51]. In contrast to other hyper reduction methods, the ECSW based on the principle of virtual work can preserve the Lagrangian structure associated with Hamilton’s principle of a system, thereby guaranteeing the numerical stability during the simulation of a ROM. The proof was provided in [49]. After obtaining the reduced mesh, it is easy to determine all the DOFs associated with the reduced mesh, which refers to as the reduced DOFs. The other Boolean matrix Lr ∈ R Nr ×N is defined to pick up the reduced DOFs from the total DOFs, where Nr denotes the number of reduced DOFs for each component. The approximations for the reduced internal force and reduced stiffness matrix can be rewritten as ( [ ( )]) Kr e˜ e˜ e˜ e˜ Fr ≈ (Lr U)T Lr Ae=1 ω F U ξ + q ˜ [ ( [ ( )]) ] , (46) Kr e˜ T e˜ e˜ e˜ T Kr ≈ (Lr U) Lr Ae=1 ω K U ξ + q L (L U) r r ˜

12

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

where A denotes the assemble operator based on the original mesh. In general, the number of reduced DOFs is much less than the number of total DOFs, i.e., Nr < N . Thus, the computational complexity when computing Kr should be of the order O(Nr wr + Nr r2 ), which is a much lower computational cost. However, numerical tests indicate that the number of finite elements in the reduced mesh is strongly related to the number of ROBs, where the sNNLS has more finite elements when each component needs more ROBs. This is analogous to the conclusion reached by Kim [53] when using the cubature scheme for hyper reduction. In addition, the sNNLS used for finding the reduced mesh can be extremely time-consuming when many finite elements are required to form the reduced mesh. The sNNLS is an offline preliminary process required for running a ROM, but its time cost should be valued in the framework of the parametric POD. 3.4.2. Advantages of component-level POD In the system-level POD, reducing r to improve the computational efficiency also decreases the accuracy of the simulation result for a ROM. However, component-level POD can reduce the number of ROBs required for each component in the ROM at the same level of simulation accuracy. In fact, the working of the componentlevel POD depends on the correlations among the snapshot collections of the components. For example, a double pendulum comprising two identical beams where the joint connecting the two beams is a spherical joint, the snapshot collections for the two beams will be almost irrelevant and the component-level POD can essentially reduce the number of ROBs required for each beam. If the joint is a weld joint, the snapshot collection for one of the beams will be highly correlated with that for the other, and although the component-level POD does not work, one can treat the two beams as a single component. For a complex FMBS, determining the bodies that comprise a “component” is an arbitrary problem that requires further research. To consider the computational complexity of the system-level POD and the component-level POD, it is assumed that an FMBS comprises m identical components, where the number of DOFs and ROBs required for each component are N b and rb , respectively. Furthermore, it is assumed further the optimal case such that there are no correlations between two arbitrary components. Thus, an identical number of ROBs will give the same level of accuracy for the system-level POD and component-level POD. The number of DOFs and ROBs required for the whole system are mN b and mrb . It is possible to show that the computational complexities of the matrix multiplications are O(N b rb m (w + rb )) and O(N b rb m (mw + m2 rb )) for the component-level POD and system-level POD, respectively. The ratio of the computational complexity between these two methods reads R=

w + rb . mw + m 2r b

(47)

In particular, when rb ≫ w, it is natural to expect that the theoretical speed-up ratio will reach m2 if the componentlevel POD is applied. Furthermore, when the ECSW is combined with the component-level POD, the computational complexity of the matrix multiplications will be reduced to O(m Nrb (wr b + (r b )2 )), where Nrb is the number of reduced DOFs for each component. Thus, one can reduce the number of DOFs N b as well as decreasing the number of ROBs rb required for each component when evaluating Krb . In addition, the cost for sNNLS will also be reduced if ECSW is applied to the component-level POD compared with the system-level POD. In fact, the number of finite elements required for the reduced mesh in each component m rb will be reduced correspondingly as rb decreases. The complexity of the matrix multiplications in sNNLS is of the order O(nrmr (mr + 1)(2mr + 1)/6) according to the algorithm given in [51], where mr is the number of the finite elements required in the reduced mesh, n is the number of snapshots, and r is the number of ROBs. More details are provided in the numerical examples, as follows. 4. Parametric POD for an FMBS This section presents the systematic MOR algorithm for a parametric FMBS. Here, the parametric POD determines how to find a proper set of ROBs for a changing parameter within a parametric domain D. The greedyPOD method is the most popular PMOR technique because of its high robustness with respect to parametric changes. This method gives a set of global ROBs suitable for all parameters under the assumption that the dynamic responses associated with all parameters fall into a single low-dimensional subspace. The greedy algorithm is employed for training and updating the ROBs in an iteration process. In particular, during each iteration, it finds the parameter

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

13

under the current ROBs that leads to the maximal error, solves the FOM associated with the parameter, and updates the ROBs with the FOM response. The subsequent paragraphs discuss three technical issues of this method. The first issue is how to measure the relative errors. The relative displacement error between the responses of the FOM and associated ROM corresponding to a parameter µ is defined as ∫T ∥q (µ) − (q + Uξ (µ))∥2 dt e (µ) = 0 . (48) ∫T ∥q dt (µ)∥ 2 0 In practice, however, the FOM responses are unknown since they are not computed. Fortunately, Eq. (3) gives rise to the residual R in the FOM (there is no residual for Φ because the constraint equations are satisfied) when computing the ROM response. It is convenient to use the residual R as an error indicator defined as ∫ T ∥R (µ)∥g dt, I (µ) = |∥R (µ)∥|g = (49) 0

with the “gappy” norm. The reason why the “gappy” norm is used (see Section 3.3) is that several constraint forces are not computable due to the disappearance of some Lagrange multipliers. The second issue is how to find the parameter that leads to the maximal error indicator. For numerical implementations, the parametric domain D is represented by a set of parameter candidates. However, if the number of candidate parameters is small, on the one hand these parameters are unlikely to contain the point in the parameter domain corresponding to the largest error indicator, thereby leading to a suboptimal parameter during each iteration. On the other hand, it will very likely lead to overfitting if these parameters are used for fitting a surface (or curve) over the parametric domain to find the largest error indicator [45]. In addition, if the candidate parameters are sufficient to obtain an optimal parameter, then computation of the ROM responses for all the candidate parameters makes the greedy-POD algorithm computationally inefficient. Therefore, Kriging interpolation [46] is employed to overcome this difficulty in order to find the maximal error indicator. It is summarized as the following four steps. 1. Start with an original sample parameter point collection and the corresponding error indicators. 2. Construct a surrogate model to predict the error indicators for all the parameter candidates. 3. Find some parameters most likely to be the maximum and add these parameters to the current sample collection. 4. Repeat steps 2 and 3 until the parameters found in step 3 are already in the sample collection. Applying Kriging interpolation to the greedy-POD method can essentially avoid the “blindness” when finding the parameter that actually leads to the maximal error indicator, thereby greatly reducing the computational effort because most of the ROM responses associated with the candidate parameters do not need to be computed. In addition, Kriging interpolation is much cheaper than computing a ROM response. Thus, the candidate parameters can be sufficient to guarantee a relatively optimal solution for the parameter that leads to the maximal error indicator. In the third issue, an increasing number of snapshots of the FOM responses is added to the snapshot matrix as the process iterates in the third step of the greedy algorithm. Thus, the application of SVD process to all snapshots will rapidly become prohibitive because the computational complexity of the SVD process to a matrix M ∈ Ra×b is O(ab2 ). To this end, the less costly approach proposed in [43] can be adopted. The snapshot matrix S in Eq. (9) is replaced by a surrogate with less columns ∑iter [ ] (50) S˜ = U1 (Λ1 )diag U2 (Λ2 )diag · · · Uiter (Λiter )diag ∈ R N × i=1 ri , where Ui and Λi (i = 1, . . . , iter) are the truncated left singular vectors and singular value vector of the ith iteration of the ∑ snapshot matrix, respectively, iter denotes the number of iterations, and ri is the number of elements of Λi . b ˜b Since iter i=1 ri is much smaller than iter × n, applying the SVD process to S will be much faster than that to S . The truncated error of the approximation to S was given in [43]. Algorithm 1 shows an updated version of the greedy algorithm with hyper reduction and component-level POD in detail. Here, Steps 7–31 comprise the outer loop as the iteration process for the greedy algorithm. The iteration stop criterion can be set to a maximal number of iterations (step 7), or to stop when the parameter that leads to the maximal error is found (steps 20–22), or when the maximal error indicator is smaller than a tolerance (steps 24–26). Steps 10–18 comprise the inner loop for Kriging interpolation. The iteration stop criterion can be set to a maximal number of iterations (step 10), or when the parameter that most likely leads to the maximal error indicator is found (steps 13–15). Steps 28–30 are for updating the ROBs and the weighting factors associated with each component according to the component-level POD described in Section 3.

14

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

15

Fig. 2. A rotating cantilever beam.

5. Numerical examples This section presents three numerical examples to demonstrate the feasibility of the PMOR process for FMBSs based on a FORTRAN compiler with an Intel® CoreTM i7-8550U @1.80 GHz and 20 GB RAM. The first example validates the effectiveness of the POD for geometrically nonlinear problems and illustrates the need of the “gappy” norm to compute the residual norm. The second example shows that the component-level POD can reduce the computational cost while retaining the same level of accuracy. The last example demonstrates the improved accuracy at the same level of computational complexity obtained via component-level POD and validates the effectiveness of the greedy-POD algorithm when applied to PMOR for FMBS. The second and the third examples also give the comparison among the FOM and four reduced models comprising the system-level ROM, system-level ROM with ECSW, component-level ROM, and component-level ROM with ECSW. 5.1. Example 1: A rotating cantilever beam As shown in Fig. 2, this example is a planer rotating cantilever beam with the left end hinged to the ground [66]. The length L, Young’s modulus E, Poisson ratio ν, the diameter d of the circular cross-section, and material density ρ of the beam were set to 1 m, 1 × 109 Pa, 0.3, 0.01 m, and 1800 kg/m3 , respectively. The angular displacement of the left end was given as ⎧ [ ( )2 ( )] 2 t T 2π t ω ⎪ s s ⎪ ⎪ + cos − 1 , 0 ≤ t ≤ Ts ⎨ Ts 2 2π Ts , (51) θ (t) = ( ) ⎪ Ts ⎪ ⎪ ⎩ωs t − , t > Ts 2 where Ts = 1.5 s and ωs = 2π rad/s. The total simulation time was set to 3 s, the time step h and the numerical damping parameter for the generalized-α method ρ∞ were set to 5 × 10−4 s and 0.8, respectively. A fairly fine mesh comprising 100 C2 -continuous cubic B-spline beam elements was used, thereby resulting in a system with 404 DOFs. 5.1.1. Comparison between POD and FFRF At first, a simple comparison was made between the POD method and FFRF. In this example, the generalized coordinates were collected every 0.005 s to form the original snapshot matrix. The truncated error ε for the singular value vector in Eq. (14) was set to smaller than 10−12 . The tolerance η for the ECSW method was set to 0.001. For comparison, 10 ROBs were used to build the ROM and were used in the FFRF. Fig. 3 shows the beam deflections of the right end obtained via the FOM, ROM, ROM with ECSW, and FFRF. The maximal beam deflection reached 0.3176 m according to the result obtained via NURBS (FOM), thereby indicating very large deformation. As expected, FFRF resulted in wrong simulation results. However, the POD method provided a fine approximation even when the ECSW technique was applied. The FFRF cannot capture the geometric nonlinearity caused by the large deformation of an FMBS, but POD is suitable for highly nonlinear systems. Thus, POD is an appropriate MOR method for geometrically nonlinear systems. 5.1.2. The “gappy” norm Now, further consideration was focused on the constraint equations. In fact, this example has two timeindependent linear constraint equations corresponding to the spherical joint and a time-dependent linear constraint corresponding to the angle displacement. In this example, the coordinates of the left end were numbered as the first

16

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

Fig. 3. Beam deflections of the right end from the FOM, ROM, ROM with ECSW, and NURBS.

two DOFs of the generalized coordinates, the two time-independent linear constraint equations can be written as [ ]T Φl (q) = q1 q2 = 02×1 . In addition, since the first two entries of each snapshot were equal to zero, the ROBs can be obtained as ]T [ (52) U = 010×2 UT , where U represents the remaining entries of the 10 ROB vectors. The Jacobian of Φl w.r.t. ξ was (Φl )ξ = (Φl )q U = [ ] I2×2 02×402 U = 02×10 . Thus, the two time-independent linear constraint equations disappeared in the ROM, thereby eliminating the associated Lagrange multipliers. As a result, the constraint force for the spherical joint was no longer available. In the ROM, two types of norms were obtained for the residual of the kinematic equation, as shown in Fig. 4, where one was the traditional two-norm and the other was the “gappy” norm. In the FOM, the norm of the constraint force corresponding to the spherical joint was computed for comparison, as shown in Fig. 5. Clearly, the constraint force made the dominant contribution to the residual when the two-norm was taken, which could not be regarded as an error indicator. 5.2. Example 2: A spatial double pendulum This example aims at demonstrating the advantages of the component-level POD compared with the system-level POD. As shown in Fig. 6, the second example is a double pendulum, where the left end of beam I was hinged to the ground with a spherical joint at the origin A, and beam II was connected with beam I at point B with a spherical joint. The length L, Young’s modulus E, Poisson ratio ν, the diameter d of the circular cross-section, and the material density ρ of each beam were set to 1 m, 1 × 109 Pa, 0.3, 0.01 m, and 1800 kg/m3 , respectively. Beam I and beam II were initially rest along the x-axis and y-axis, respectively. The system started falling down under the action of gravity along the negative direction of the z-axis, and the acceleration due to gravity was set to 9.81 m/s2 . Each beam was meshed via 100 C1 -continuous cubic B-spline beam elements, thereby resulting in a system with 1212 DOFs and six constraint equations (time-independent and linear). The total simulation time was set to 3 s, the time step h and spectra radius ρ∞ for the generalized-α method were set to 5 × 10−4 s and 0.8, respectively. 5.2.1. Advantages of component-level POD: computational cost Similar to Example 1, the generalized coordinates were collected every 0.005 s to form the original snapshot matrix. The truncated error ε of the singular value vector was set smaller than 10−12 for both system-level ROM and component-level ROM, and the tolerance η was set to 10−3 for the ECSW method in both system-level ROM and component-level ROM. The FOM and four reduced models were computed. The configurations at six instants in the FOM results are presented in Fig. 7. The z-displacements of point C obtained from the five models are shown in Fig. 8. Clearly, the simulation results obtained from four reduced models were highly consistent with the results produced by the FOM. Furthermore, the relative errors er between the FOM and four reduced models are shown in Fig. 9. The relative error between the FOM and an ROM at the instant t is defined as er = ∥q (t) − (q + Uξ (t))∥ / ∥q (t)∥ .

(53)

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

17

Fig. 4. Two norms for the residual of the kinematic equation in the ROM.

Fig. 5. Norm of the constraint force corresponding to the spherical joint in the FOM.

Fig. 6. A double pendulum.

The relative errors were slightly smaller with the two component-level ROMs (with or without ECSW) than those obtained from the two system-level ROMs. However, the two component-level ROMs reduced the computational cost in matrix multiplication and sNNLS essentially compared to the two system-level ROMs. Based on the analysis given in Section 3.4.2, the number of scalar multiplications among matrix multiplication can be estimated for computing Kr , before simulating the reduced models. Table 1 lists the number r of ROBs, the number N of DOFs, the number Nr of reduced DOFs, the bandwidth w , the number p of scalar multiplications, and

18

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

Fig. 7. Configurations of the FOM at six instants.

Fig. 8. The z-displacements of point C obtained from five models.

Fig. 9. Errors between the FOM and the reduced models.

the theoretical speed-up ratios s between the system-level ROM and the other three reduced models. The systemlevel ROM needed 59 ROBs and had to deal with all finite elements to evaluate Kr . Using ECSW, the number of finite elements of concern was reduced from 200 to 36, thereby decreasing the number of DOFs from 1212 to 432. The speed-up ratio for the number of scalar multiplications reached 2.81. In the component-level ROM, the

19

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690 Table 1 Comparison of the number of scalar multiplications among the matrix multiplications for computing Kr for the four reduced models. Reduced model

r

K(Kr )

N(Nr )

w

P

s

System-level ROM

–

59

200

1212

18

5 506 116

–

System-level ROM with ECSW

–

59

36

432

18

1 962 576

2.81

Component-level ROM

Beam I Beam II

33 32

100 100

606 606

18 18

1 019 898 969 600

Component-level ROM with ECSW

Beam I Beam II

33 32

18 18

216 216

18 18

363 528 345 600

1 989 498

2.77

709 128

7.76

Table 2 Comparison on the time required to compute Kr in each Newton–Raphson iteration for the four reduced models. Reduced model System-level ROM System-level ROM with ECSW Component-level ROM Component-level ROM with ECSW

Time required to compute Kr in each iteration (s) 2.348 6.606 8.861 2.656

Ratio

10−3

– 3.55 2.65 8.84

× × 10−4 × 10−4 × 10−4

Table 3 Comparison of the time required for sNNLS in ROM with the ECSW model and component-level ROM with the ECSW model. Reduced model

Time for sNNLS (s)

Ratio –

ROM with ECSW

–

1.1964

Component-level ROM with ECSW

Beam I Beam II

0.1685 0.1344

0.3029

3.95

number of ROBs for each beam was nearly a half of that in the system-level ROM (32 and 33 vs. 59). Thus, the number of scalar multiplications among matrix multiplications was greatly reduced, with a speed-up ratio of 2.77. Furthermore, when the ECSW was used, the number of finite elements of concern decreased (100 and 100 vs. 18 and 18). So did the number of DOFs for matrix multiplications (606 and 606 vs. 216 and 216). The speed-up ratio reached 7.76 with this reduced model because of a greater reduction in the number of scalar multiplications. The times required to compute Kr in each Newton–Raphson iteration are listed in Table 2, which shows that the practical speed-up ratio for the component-level POD based on the matrix multiplications was 2.65, and the combination of the ECSW and component-level POD could achieve a practical speed-up ratio of 8.84. The test results were consistent with the analysis of the matrix multiplications. This gives an assertion that both ECSW and component-level POD could help to reduce the computational cost, and the combination of these two methods reached even greater reductions, thereby decreasing the simulation time for large systems. In addition, in the component-level ROM, the number of finite elements of the reduced mesh in each beam was a half of the number of finite elements in the system-level ROM (18 and 18 vs. 36), mainly because the number of ROBs for each beam was nearly a half of that in the ROM. Table 3 shows the times required to run the sNNLS in the system-level ROM with ECSW and the component-level ROM with ECSW. The test results show that the practical speed-up ratio was 3.95. This indicates that if the ECSW is combined with the component-level POD, the computational cost can be greatly reduced for the sNNLS process. 5.3. Example 3: A crank–slider mechanism As shown in Fig. 10, the left end of the crank was fixed to the ground through revolute joint A with a diameter of 0.03 m. The linkage was connected with the crank via revolute joint B with a diameter of 0.02 m. The slider was attached to the linkage via revolute joint C with a diameter of 0.02 m. The slider could only move along the x direction. Other required geometric sizes are shown in Fig. 10. The thickness t, material density ρ, Young’s modulus E, and Poisson ratio ν for the crank and linkage were set to 0.01 m, 2700 kg/m3 , 70 GPa, and 0.3, respectively.

20

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

Fig. 10. A crank–slider mechanism.

All revolute joints were assumed to be massless, but the slider was assumed to be 0.5 kg. The torque M = 15 N m was applied to the crank. In the initial configuration, the crank and linkage were placed along the x direction. Both crank and linkage were densely meshed with gradient-deficient plate elements of ANCF [63], thereby resulting in 1130 finite elements and 18252 DOFs for the crank, and 1900 finite elements and 30 189 DOFs for the linkage, with a total of 3030 finite elements, 48 453 DOFs, and 548 nonlinear and 2 linear constraint equations for the whole system. The total simulation time was set as 0.2 s. The time step h and the numerical damping parameter for the generalized-α method ρ∞ were set to 1 × 10−5 s and 0.8, respectively. 5.3.1. Advantage of the component-level POD: accuracy First, a comparison was made for the accuracy of the component-level POD and system-level POD. The generalized coordinates at every 3 × 10−4 s during the first 0.15 s and at every 5 × 10−5 s in the last 0.05 s were collected to form the original snapshot matrix. The truncated error ε of the singular value vector was set smaller than 10−18 for the ROM and component-level ROM, and the tolerance η was set to 10−6 for the ECSW in the ROM and component-level ROM. Similar to Example 2, the FOM and four reduced models were computed. In Table 4, the first five rows list the number of DOFs, the bandwidth, the number of ROBs, the number of finite elements in the reduced meshes, and the number of reduced DOFs for all five models. The following five rows list the times required to compute the generalized forces and Jacobians, for solving the linear equations, matrix multiplications, simulating the whole models, and running the sNNLS for ECSW. The last row shows the maximal relative error between the reduced models and the FOM. It can be seen that the two ROMs without ECSW cannot save computational effort because of the matrix multiplications, and the speed-up ratios of running the two ROMs with ECSW are about 4 compared with running the FOM. Given the same truncated error ε = 10−18 , although the component-level ROM needed more ROBs (20 and 45, vs. 46), the time required for matrix multiplications greatly decreased (2378.3 s vs. 3111.8 s) without ECSW, and increased slightly (645.6 s vs. 609.7 s) with ECSW. However, the component-level ROM was much more accurate than the ROM without ECSW (6.0574e−7 vs. 2.4335e−5) or with ECSW (6.2388e−7 vs. 2.4477e−5). The component-level POD improved the accuracy of the reduced simulation results by about 40 times at the same computational cost. In addition, the time required for the extremely time-consuming sNNLS decreased (1082.2 s and 85.4 s, vs. 3632.7 s) for the component-level POD instead of the system-level POD. The x-displacements of point C from the five models are shown in Fig. 11, which indicates that the maximal rotational speed reached about 6000 rpm. The relative errors between the FOM and the reduced models are presented in Fig. 12. They show that the component-level ROM was much more accurate than the system-level ROM throughout the whole simulation process because the component-level ROM used more ROBs, i.e., more DOFs, to compute the dynamic response. 5.3.2. Greedy algorithm for PMOR The greedy algorithm was used in this example in order to obtain the robust ROBs with parametric changes. According to Algorithm 1 in Section 4, the iteration process for Kriging interpolation (Kriging iteration) was

21

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690 Table 4 Comparison between the POD and component-level POD in terms of the computational efficiency and accuracy. Model

FOM

ROM

Component-level ROM Crank

DOFs Bandwidth ROBs Elem. in Red. mesh Red. DOFs

48 453 99 – – –

48 453 99 46 – –

Force and Jacobian (s) Solve linear Eq. (s) Mat. Multi. (s) Total time (s) sNNLS (s)

3843.7 2482.7 – 6819.8 –

3652.9 7.9 3111.8 6879.1 –

Max. relative error

–

2.4335e−5

ROM with ECSW

Linkage

18 258

30 195

Component-level ROM with ECSW Crank

– –

48 453 99 46 541 17 013

3810.4 11.9 2378.3 6633.8 –

783.2 10.5 609.7 1555.5 3632.7

85.4

2.4477e−5

6.2388e−7

99 20

45

6.0574e−7

Linkage

18 258

30 195 99

20 208 6669

45 344 10 584 756.4 11.6 645.6 1681.9 1082.2

Fig. 11. The x-displacements of point C from five models.

Fig. 12. Relative errors between the FOM and reduced models.

nested into the iteration of the greedy algorithm (greedy iteration). In this example, the Young’s modulus E or the loading torque M were assumed to change among 101 candidate parameters distributed evenly in each range and the maximal number of sample parameters for Kriging iteration Nsmax was set to 15. In order to compare the greedy-POD with the fast-to-slow ROM mentioned in [31], the relative errors of 26 parameters distributed evenly in each range were computed. The truncated error ε for SVD was 10−18 and the tolerance η for the ECSW was 10−6 .

22

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

Fig. 13. Relative errors for 26 parameters under twice greedy iterations and the fast-to-slow ROM with respect to changes in Young’s modulus.

Fig. 14. Relative errors for 26 parameters under twice greedy iterations and the fast-to-slow ROM with respect to changes in the loading torque.

(a) Changing Young’s modulus E. The initial sample parameters of Young’s modulus of the crank and linkage were selected as {40, 47.5, 55, 62.5, 70} GPa from the range of [40, 70] GPa for the Kriging iteration, and E = 55 GPa was set as the initial sample parameter for the greedy iteration. During the greedy process, it underwent 2 runs of the FOM, 12 runs of the component-level ROM with ECSW, and 2 sNNLS processes to obtain the final reduced model. Fig. 13 shows the relative errors obtained from Eq. (48) for the 26 parameters under twice greedy iterations and the fast-to-slow ROM. (b) Changing loading torque M. The loading moment M was set in the range of [5,15] N m. Values of {5, 7.5, 10, 12.5, 15} N m were selected as the initial sample parameters for the Kriging iteration. M = 10.5 N m was selected as the initial sample parameter for the greedy iteration. During the greedy process, it underwent 2 runs of the FOM, 21 runs of the component-level ROM with ECSW, and 2 sNNLS processes to obtain the final reduced model. Fig. 14 shows the relative errors for the 26 parameters under twice greedy iterations and the fast-to-slow ROM. The numerical simulations show that after greedy iterations, the ROBs obtained for PMOR were more robust and accurate as the parameters changed. Furthermore, the greedy algorithm performed much better than the fast-to-slow ROM method, which failed to obtain robust ROBs. Given the final ROBs obtained by the greedy algorithm, it was possible to rapidly determine the dynamic response corresponding to any Young’s modulus (or loading moment) in each range, and it was not necessary to simulate the original FOM.

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

23

6. Conclusions In this paper, a systematic study is made on the parametric MOR to simulate the dynamic responses of flexible multibody systems modeled via IGA and ANCF. The dynamic equations for the reduced order model are constructed by projecting the original dynamic equations for an FBMS onto the low-dimensional subspace spanned by the ROBs produced via SVD with the snapshot collection for each component, that is, the component-level POD. It is emphasized that the redundancies of the reduced constraint equations and the original constraint equations are different. Thus, the Gauss elimination method is used to eliminate the redundant reduced constraint equations from the ROM of an FMBS. The “gappy” norm of the residual is employed to quantify the size of the residual vector because the residual of the FOM cannot be computed due to the disappearance of the Lagrange multipliers corresponding to the constraint equations removed from the original problem. Furthermore, the POD fails to facilitate efficient simulations because the matrix multiplications are extremely time consuming when evaluating the reduced stiffness matrix. Thus, the ECSW is used as a type of hyper reduction technique. In addition, the component-level POD is used to SVD to each component of a system and joins them up based on the constraint equations among components, thereby reducing the computational cost of the matrix multiplications. Three numerical examples validate the effectiveness of the ECSW and component-level POD and demonstrate that the combination of both can further reduce the matrix multiplications. In addition, the componentlevel POD can improve the accuracy of the reduced simulation result at the same level of computational complexity. Finally, the greedy algorithm with Kriging interpolation is applied to MOR when the system is parameterized. The last example demonstrates that the ROBs obtained by using the greedy algorithm were robust to parametric changes. Further studies may focus on how to optimally select the components for an FMBS. Acknowledgments This work was supported in part by National Natural Science Foundations of China under Grants 11672034 and 11832005. References [1] A.A. Shabana, R. Schwertassek, Equivalence of the floating frame of reference approach and finite element formulations, Int. J. Nonlinear Mech. 33 (1998) 417–432, http://dx.doi.org/10.1016/S0020-7462(97)00024-3. [2] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (2005) 4135–4195, http://dx.doi.org/10.1016/j.cma.2004.10.008. [3] A.A. Shabana, H.A. Hussien, J.L. Escalona, Application of the absolute nodal coordinate formulation to large rotation and large deformation problems, J. Mech. Des. 120 (1998) 188–195, http://dx.doi.org/10.1115/1.2826958. [4] J.C. Simo, L. Vu-Quoc, On the dynamics in space of rods undergoing large motions - A geometrically exact approach, Comput. Methods Appl. Mech. Engrg. 66 (1988) 125–161, http://dx.doi.org/10.1016/0045-7825(88)90073-4. [5] J.C. Simo, D.D. Fox, On a stress resultant geometrically exact shell model. Part I: Formulation and optimal parametrization, Comput. Methods Appl. Mech. Engrg. 72 (1989) 267–304, http://dx.doi.org/10.1016/0045-7825(89)90002-9. [6] R.R. Craig JR, M.C.C. Bampton, Coupling of substructures for dynamic analyses, AIAA J. 6 (1968) 1313–1319, http://dx.doi.org/10. 2514/3.4741. [7] T.M. Wasfy, A.K. Noor, Computational strategies for flexible multibody systems, Appl. Mech. Rev. 56 (2003) 553–613, http: //dx.doi.org/10.1115/1.1590354. [8] O. Brüls, P. Duysinx, J.C. Golinval, The global modal parameterization for non-linear model-order reduction in flexible multibody dynamics, Internat. J. Numer. Methods Engrg. 69 (2007) 948–977, http://dx.doi.org/10.1002/nme.1795. [9] F. Naets, G.H.K. Heirman, D. Vandepitte, W. Desmet, Inertial force term approximations for the use of Global Modal Parameterization for planar mechanisms, Internat. J. Numer. Methods Engrg. 85 (2011) 518–536, http://dx.doi.org/10.1002/nme.2984. [10] G.H.K. Heirman, F. Naets, W. Desmet, A system-level model reduction technique for the efficient simulation of flexible multibody systems, Internat. J. Numer. Methods Engrg. 85 (2011) 330–354, http://dx.doi.org/10.1002/nme.2971. [11] J. Fehr, P. Eberhard, Simulation process of flexible multibody systems with non-modal model order reduction techniques, Multibody Syst. Dyn. 25 (2010) 313–334, http://dx.doi.org/10.1007/s11044-010-9238-3. [12] M. Fischer, P. Eberhard, Linear model reduction of large scale industrial models in elastic multibody dynamics, Multibody Syst. Dyn. 31 (2014) 27–46, http://dx.doi.org/10.1007/s11044-013-9347-x. [13] M. Baumann, D. Hamann, P. Eberhard, Time-dependent parametric model order reduction for material removal simulations, in: P. Benner, M. Ohlberger, A. Patera, G. Rozza, K. Urban (Eds.), Model Reduct. Parametr. Syst., Springer, Cham, 2017, pp. 491–504, http://dx.doi.org/10.1007/978-3-319-58786-8_30. [14] M. Fischer, P. Eberhard, Application of parametric model reduction with matrix interpolation for simulation of moving loads in elastic multibody systems, Adv. Comput. Math. 41 (2015) 1049–1072, http://dx.doi.org/10.1007/s10444-014-9379-7.

24

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

[15] L. Wu, P. Tiso, Nonlinear model order reduction for flexible multibody dynamics: a modal derivatives approach, Multibody Syst. Dyn. 36 (2016) 405–425, http://dx.doi.org/10.1007/s11044-015-9476-5. [16] S.R. Idelsohn, A. Cardona, A reduction method for nonlinear structural dynamic analysis, Comput. Methods Appl. Mech. Engrg. 49 (1985) 253–279, http://dx.doi.org/10.1016/0045-7825(85)90125-2. [17] S. Jain, P. Tiso, J.B. Rutzmoser, D.J. Rixen, A quadratic manifold for model order reduction of nonlinear structural dynamics, Comput. Struct. 188 (2017) 80–94, http://dx.doi.org/10.1016/j.compstruc.2017.04.005. [18] J.B. Rutzmoser, D.J. Rixen, P. Tiso, S. Jain, Generalization of quadratic manifolds for reduced order modeling of nonlinear structural dynamics, Comput. Struct. 192 (2017) 196–209, http://dx.doi.org/10.1016/j.compstruc.2017.06.003. [19] C.S.M. Sombroek, P. Tiso, L. Renson, G. Kerschen, Numerical computation of nonlinear normal modes in a modal derivative subspace, Comput. Struct. 195 (2018) 34–46, http://dx.doi.org/10.1016/j.compstruc.2017.08.016. [20] G. Kerschen, J.C. Golinval, A.F. Vakakis, L.A. Bergman, The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: An overview, Nonlinear Dynam. 41 (2005) 147–169, http://dx.doi.org/10.1007/s11071-0052803-2. [21] D.B.P. Huynh, F. Pichi, G. Rozza, Reduced basis approximation and a posteriori error estimation: Applications to elasticity problems in several parametric settings, in: D.A. Di Pietro, A. Ern, L. Formaggia (Eds.), Numer. Methods PDEs, Springer, Cham, 2018, pp. 203–247, http://dx.doi.org/10.1007/978-3-319-94676-4_8. [22] M. Kärcher, S. Boyaval, M.A. Grepl, K. Veroy, Reduced basis approximation and a posteriori error bounds for 4D-Var data assimilation, Optim. Eng. 19 (2018) 663–695, http://dx.doi.org/10.1007/s11081-018-9389-2. [23] I. Martini, B. Haasdonk, G. Rozza, Certified reduced basis approximation for the coupling of viscous and inviscid parametrized flow models, J. Sci. Comput. 74 (2018) 197–219, http://dx.doi.org/10.1007/s10915-017-0430-y. [24] M.A. Dihlmann, B. Haasdonk, Certified PDE-constrained parameter optimization using reduced basis surrogate models for evolution problems, Comput. Optim. Appl. 60 (2015) 753–787, http://dx.doi.org/10.1007/s10589-014-9697-1. [25] N.C. Nguyen, G. Rozza, D.B.P. Huynh, A.T. Patera, Reduced basis approximation and a posteriori error estimation for parametrized parabolic PDEs: Application to real-time Bayesian parameter estimation, in: L. Biegler, G. Biros, O. Ghattas, M. Heinkenschloss, D. Keyes, B. Mallick, Y. Marzouk, L. Tenorio, B. van B. Waanders, K. Willcox (Eds.), Large-Scale Inverse Probl. Quantif. Uncertain., John Wiley & Sons, Ltd, Chichester, UK, 2010, pp. 151–177, http://dx.doi.org/10.1002/9780470685853.ch8. [26] M. Ulbrich, B. van Bloemen Waanders, An introduction to partial differential equations constrained optimization, Optim. Eng. 19 (2018) 515–520, http://dx.doi.org/10.1007/s11081-018-9398-1. [27] D. Amsallem, M. Zahr, Y. Choi, C. Farhat, Design optimization using hyper-reduced-order models, Struct. Multidiscip. Optim. 51 (2015) 919–940, http://dx.doi.org/10.1007/s00158-014-1183-y. [28] Y. Yue, K. Meerbergen, Accelerating optimization of parametric linear systems by model order reduction, SIAM J. Optim. 23 (2013) 1344–1370, http://dx.doi.org/10.1137/120869171. [29] P. LeGresley, J. Alonso, Airfoil design optimization using reduced order models based on proper orthogonal decomposition, in: Fluids 2000 Conf. Exhib., American Institute of Aeronautics and Astronautics, 2000, http://dx.doi.org/10.2514/6.2000-2545. [30] H. Hotelling, Analysis of a complex of statistical variables into principal components, J. Educ. Psychol. 24 (1933) 417–441, http://dx.doi.org/10.1037/h0071325. [31] K. Luo, H. Hu, C. Liu, Q. Tian, Model order reduction for dynamic simulation of a flexible multibody system via absolute nodal coordinate formulation, Comput. Methods Appl. Mech. Engrg. 324 (2017) 573–594, http://dx.doi.org/10.1016/j.cma.2017.06.029. [32] A. Radermacher, S. Reese, Model reduction in elastoplasticity: Proper orthogonal decomposition combined with adaptive sub-structuring, Comput. Mech. 54 (2014) 677–687, http://dx.doi.org/10.1007/s00466-014-1020-6. [33] A. Corigliano, M. Dossi, S. Mariani, Model order reduction and domain decomposition strategies for the solution of the dynamic elastic–plastic structural problem, Comput. Methods Appl. Mech. Engrg. 290 (2015) 127–155, http://dx.doi.org/10.1016/j.cma.2015.02. 021. [34] A. Radermacher, S. Reese, POD-based model reduction with empirical interpolation applied to nonlinear elasticity, Internat. J. Numer. Methods Engrg. 107 (2016) 477–495, http://dx.doi.org/10.1002/nme.5177. [35] F. Ghavamian, P. Tiso, A. Simone, POD–DEIM Model order reduction for strain softening viscoplasticity, Comput. Methods Appl. Mech. Engrg. 317 (2017) 458–479, http://dx.doi.org/10.1016/j.cma.2016.11.025. [36] T. Lieu, C. Farhat, M. Lesoinne, Reduced-order fluid/structure modeling of a complete aircraft configuration, Comput. Methods Appl. Mech. Engrg. 195 (2006) 5730–5742, http://dx.doi.org/10.1016/j.cma.2005.08.026. [37] D. Amsallem, C. Farhat, Interpolation method for adapting reduced-order models and application to aeroelasticity, AIAA J. 46 (2008) 1803–1813, http://dx.doi.org/10.2514/1.35374. [38] P. Astrid, S. Weiland, K. Willcox, T. Backx, Missing point estimation in models described by proper orthogonal decomposition, IEEE Trans. Automat. Control 53 (2008) 2237–2251, http://dx.doi.org/10.1109/TAC.2008.2006102. [39] D. Amsallem, M.J. Zahr, K. Washabaugh, Fast local reduced basis updates for the efficient reduction of nonlinear systems with hyper-reduction, Adv. Comput. Math. 41 (2015) 1187–1230, http://dx.doi.org/10.1007/s10444-015-9409-0. [40] K. Washabaugh, D. Amsallem, M. Zahr, C. Farhat, Nonlinear model reduction for CFD problems using local reduced-order bases, in: 42nd AIAA Fluid Dyn. Conf. Exhib., American Institute of Aeronautics and Astronautics, 2012, http://dx.doi.org/10.2514/6.2012-2686. [41] K. Willcox, Unsteady flow sensing and estimation via the gappy proper orthogonal decomposition, Comput. Fluids 35 (2006) 208–226, http://dx.doi.org/10.1016/j.compfluid.2004.11.006. [42] K. Carlberg, C. Farhat, J. Cortial, D. Amsallem, The GNAT method for nonlinear model reduction: Effective implementation and application to computational fluid dynamics and turbulent flows, J. Comput. Phys. 242 (2013) 623–647, http://dx.doi.org/10.1016/j.jcp. 2013.02.028.

Y. Hou, C. Liu and H. Hu / Computer Methods in Applied Mechanics and Engineering 361 (2020) 112690

25

[43] A. Paul-Dubois-Taine, D. Amsallem, An adaptive and efficient greedy procedure for the optimal training of parametric reduced-order models, Internat. J. Numer. Methods Engrg. 102 (2015) 1262–1292, http://dx.doi.org/10.1002/nme.4759. [44] P. Benner, S. Gugercin, K. Willcox, A survey of projection-based model reduction methods for parametric dynamical systems, SIAM Rev. 57 (2015) 483–531, http://dx.doi.org/10.1137/130932715. [45] J. Fehr, M. Fischer, B. Haasdonk, P. Eberhard, Greedy-based approximation of frequency-weighted Gramian matrices for model reduction in multibody dynamics, ZAMM - J. Appl. Math. Mech. / Z. Angew. Math. Mech. 93 (2013) 501–519, http://dx.doi.org/10. 1002/zamm.201200014. [46] D.R. Jones, A taxonomy of global optimization methods based on response surfaces, J. Global Optim. 21 (2001) 345–383, http://dx.doi.org/10.1023/A:1012771025575. [47] T. Bui-Thanh, K. Willcox, O. Ghattas, Parametric reduced-order models for probabilistic analysis of unsteady aerodynamic applications, AIAA J. 46 (2008) 2520–2529, http://dx.doi.org/10.2514/1.35850. [48] T. Bui-Thanh, K. Willcox, O. Ghattas, Model reduction for large-scale systems with high-dimensional parametric input space, SIAM J. Sci. Comput. 30 (2008) 3270–3288, http://dx.doi.org/10.1137/070694855. [49] C. Farhat, T. Chapman, P. Avery, Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models, Internat. J. Numer. Methods Engrg. 102 (2015) 1077–1110, http://dx.doi.org/10.1002/nme.4820. [50] S. Chaturantabut, D.C. Sorensen, Nonlinear model reduction via discrete empirical interpolation, SIAM J. Sci. Comput. 32 (2010) 2737–2764, http://dx.doi.org/10.1137/090766498. [51] C. Farhat, P. Avery, T. Chapman, J. Cortial, Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency, Internat. J. Numer. Methods Engrg. 98 (2014) 625–662, http://dx.doi.org/10.1002/nme.4668. [52] T. Chapman, P. Avery, P. Collins, C. Farhat, Accelerated mesh sampling for the hyper reduction of nonlinear computational models, Internat. J. Numer. Methods Engrg. 109 (2017) 1623–1654, http://dx.doi.org/10.1002/nme.5332. [53] T. Kim, D.L. James, Skipping steps in deformable simulation with online model reduction, ACM Trans. Graph. 28 (2009) 1, http://dx.doi.org/10.1145/1618452.1618469. [54] D. Stadlmayr, W. Witteveen, W. Steiner, Reduction of physical and constraint degrees-of-freedom of redundant formulated multibody systems, J. Comput. Nonlinear Dyn. 11 (2016) 1–9, http://dx.doi.org/10.1115/1.4031553. [55] D. Stadlmayr, W. Witteveen, W. Steiner, A generalized constraint reduction method for reduced order MBS models, Multibody Syst. Dyn. 41 (2017) 259–274, http://dx.doi.org/10.1007/s11044-016-9557-0. [56] P. Kerfriden, O. Goury, T. Rabczuk, S.P.A. Bordas, A partitioned model order reduction approach to rationalise computational expenses in nonlinear fracture mechanics, Comput. Methods Appl. Mech. Engrg. 256 (2013) 169–188, http://dx.doi.org/10.1016/j.cma.2012.12.004. [57] S. Erlicher, L. Bonaventura, O.S. Bursi, The analysis of the Generalized -α method for non-linear dynamic problems, Comput. Mech. 28 (2002) 83–104, http://dx.doi.org/10.1007/s00466-001-0273-z. [58] M. Arnold, O. Brüls, Convergence of the generalized-α scheme for constrained mechanical systems, Multibody Syst. Dyn. 18 (2007) 185–202, http://dx.doi.org/10.1007/s11044-007-9084-0. [59] A.M. Mikkola, A.A. Shabana, A non-incremental finite element procedure for the analysis of large deformation of plates and shells in mechanical system applications, Multibody Syst. Dyn. 9 (2003) 283–309, http://dx.doi.org/10.1023/A:1022950912782. [60] J. Gerstmayr, A.A. Shabana, Analysis of thin beams and cables using the absolute nodal co-ordinate formulation, Nonlinear Dynam. 45 (2006) 109–130, http://dx.doi.org/10.1007/s11071-006-1856-1. [61] A. Borkovi´c, S. Kovaˇcevi´c, G. Radenkovi´c, S. Milovanovi´c, D. Majstorovi´c, Rotation-free isogeometric dynamic analysis of an arbitrarily curved plane Bernoulli–Euler beam, Eng. Struct. 181 (2019) 192–215, http://dx.doi.org/10.1016/j.engstruct.2018.12.003. [62] A. Mikkola, A.A. Shabana, C. Sanchez-Rebollo, J.R. Jimenez-Octavio, Comparison between ANCF and B-spline surfaces, Multibody Syst. Dyn. 30 (2013) 119–138, http://dx.doi.org/10.1007/s11044-013-9353-z. [63] K. Dufva, A.A. Shabana, Analysis of thin plate structures using the absolute nodal coordinate formulation, Proc. Inst. Mech. Eng. K J. Multi-Body Dyn. 219 (2005) 345–355, http://dx.doi.org/10.1243/146441905X50678. [64] C. Liu, Q. Tian, D. Yan, H. Hu, Dynamic analysis of membrane systems undergoing overall motions, large deformations and wrinkles via thin shell elements of ANCF, Comput. Methods Appl. Mech. Engrg. 258 (2013) 81–95, http://dx.doi.org/10.1016/j.cma.2013.02.006. [65] C. Liu, Q. Tian, H. Hu, New spatial curved beam and cylindrical shell elements of gradient-deficient Absolute Nodal Coordinate Formulation, Nonlinear Dynam. 70 (2012) 1903–1918, http://dx.doi.org/10.1007/s11071-012-0582-0. [66] S.-C. Wu, E.J. Haug, Geometric non-linear substructuring for dynamics of flexible mechanical systems, Internat. J. Numer. Methods Engrg. 26 (1988) 2211–2226, http://dx.doi.org/10.1002/nme.1620261006.