An enhanced parametric nonlinear reduced order model for imperfect structures using Neumann expansion
Jacopo Marconi, Paolo Tiso, Davide E. Quadrelli, Francesco Braghin
AAn enhanced parametric nonlinear reduced order model for imperfectstructures using Neumann expansion
Jacopo Marconi a , Paolo Tiso b, ∗ , Davide E. Quadrelli a , Francesco Braghin a a Department of Mechanical Engineering, Politecnico di Milano, Via La Masa 1, 20156 Milan, Italy b Institute for Mechanical Systems, ETH Z¨urich, Leonhardstrasse 21, 8092 Z¨urich, Switzerland
Abstract
We present an enhanced version of the parametric nonlinear reduced order model for shape imperfections instructural dynamics we studied in a previous work [1]. The model is computed intrusively and with no training using information about the nominal geometry of the structure and some user-defined displacement fieldsrepresenting shape defects , i.e. small deviations from the nominal geometry parametrized by their respectiveamplitudes. The linear superposition of these artificial displacements describe the defected geometry and canbe embedded in the strain formulation in such a way that, in the end, nonlinear internal elastic forces can beexpressed as a polynomial function of both these defect fields and the actual displacement field. This way, atensorial representation of the internal forces can be obtained and, owning the reduction in size of the modelgiven by a Galerkin projection, high simulation speed-ups can be achieved. We show that by adopting arigorous deformation framework we are able to achieve better accuracy as compared to the previous work. Inparticular, exploiting
Neumann expansion in the definition of the Green-Lagrange strain tensor, we show thatour previous model is a lower order approximation with respect to the one we present now. Two numericalexamples of a clamped beam and a MEMS gyroscope finally demonstrate the benefits of the method in termsof speed and increased accuracy.
Keywords:
Nonlinear Modeling, Reduced Order Models, Parametric, Geometric Nonlinearities, Defects
1. Introduction
The Finite Element (FE) method has long been a fundamental analysis and design tool in many areas ofscience and engineering. In structural mechanics it is almost mandatory to use FE models to investigate thebehavior of complex systems, which often have many geometric details that would be difficult to handle withalternative approaches, such as lumped parameter or analytical models [2]. However, large FE simulationswould often require considerable computational resources and time, so in some cases designers may prefer toperform real experiments rather than numerical ones. On the one hand, this need for fast and affordable FEsimulations has given rise to numerical techniques to improve computational efficiency: domain decomposition ∗ Corresponding author
Email address: [email protected] (Paolo Tiso)
Preprint submitted to arXiv February 4, 2021 a r X i v : . [ c s . C E ] F e b nd substructuring [3, 4] and FE Tearing and Interconnecting (FETI, [5]) are just a few examples. On theother hand, model order reduction methods have emerged, consisting in the construction of a Reduced OrderModel (ROM), whose number of degrees of freedom (dofs) is much smaller than that of the reference FullOrder Model (FOM). The use of linear ROMs also in industrial contexts is nowadays well established asthe theory underlying them. Guyan reduction [6] and modal analysis [7] are two well-known examples inmechanical statics and dynamics, respectively, where FOM’s static deformations and Vibration Modes (VMs,also known as eigenmodes or natural modes of the linear system) are used to construct a Reduced Basis (RB)that projects the governing equations onto a lower dimension subspace. Linear ROMs were also successfullycoupled with substructuring techniques in the Craig-Bampton and Rubin methods [8, 9], which are availablein many commercial software.For nonlinear FE studies, where the demand for reduction is dire, many solutions have been proposedover the last decades, but none of them seems to have prevailed over the others, as each of them offerscertain advantages, requires certain costs and/or targets specific problems. Overall, however, the literature ismature enough to provide the analyst with many different options in several practical applications, rangingfrom bolted joints [10], gears [11], contacts [12, 13], friction [14] and viscoplasticity [15] to flexible multi-bodydynamics with geometric nonlinearities [16] and substructuring [17].Nonlinear ROMs can be classified according to (i) whether they are RB-projection based or not, (ii)whether they are data- or model-driven and (iii) their (non-)intrusiveness. In the following we considermostly projection approaches, as the one adopted in this work; alternatively, one could resort to differentstrategies, such as normal form theory or Spectral Submanifolds. The most recent contributions in this senseinclude [18] and [19, 20]. In (ii), for data-driven ROMs we usually refer to ROMs constructed using previousFOM simulation data (or experimental data, [21]), as opposed to model-driven methods that rely on someintrinsic properties of the model itself for ROM construction, such as modal approaches [22, 23, 24, 25]. Asfor intrusiveness, we usually denote a ROM as non-intrusive [26] if it can be used with routines and solversof commercial FE software and, conversely, as intrusive a method requiring dedicated routines. Specifically,intrusive methods require access and manipulation to element-level quantities, as for instance nonlineargeneralized forces and jacobians. Other distinctions can be made in terms of the types of nonlinearities thata given model can handle and the way nonlinear functions are evaluated [27]. All these differences ultimatelyaffect the two phases that all ROMs have in common: the offline phase, in which the ROM is constructed,and the online phase, in which the simulation responses are retrieved. As the main goal of ROMs is toreduce computational effort and time, a key aspect to keep in mind when choosing a method is the overheadcost to pay in the offline phase; in the case of data-driven methods, this cost can be as high as the costassociated to the solution of the FOM [12]. Generally speaking then, data-driven methods (usually basedon Proper Orthogonal Decomposition, or POD, strategies [28]) are used in scenarios where the high costassociated to the data generation can be amortized: typically, this is the case of multi-query analysis. In thissense, although not as versatile and generally applicable as data-driven POD-based approaches, model-driven2trategies in structural dynamics are desirable, for no FOM simulation is required a priori. Rayleigh-Ritzprocedures [29], dual modes [26] and Modal Derivatives (MDs) [30, 31, 32] are some popular examples.One way to mitigate the offline overhead costs of all the aforementioned methods, but especially thedata-driven ones, is to resort to (nonlinear) parametric ROMs, (NL-)pROMs. Also in this context, theliterature on linear systems is quite well developed and consolidated. An extensive survey and comparison ofthese methods can be found in [33, 34]. The reduction of nonlinear parametric Partial Differential Equations(PDEs) is instead still an active research topic, which has attracted increasing interest in various disciplinesover the years. Interestingly, the vast majority of nonlinear parametric model order reduction methods isdata-driven, POD-based. Some recent examples include non-intrusive interpolation methods for evaluatingnonlinear functions with hypersurfaces [35, 36] and use of Gaussian Processes and machine learning for errorevaluation and refinement of the pROM [37] or interpolation on the Grassman manifold via tangent spaces[38]. Alternatively, many of these methods approximate the nonlinear function using hyper-reduction methodsas the Discrete Empirical Interpolation Method (DEIM) [39, 40] to speed up the evaluation, and in this senseonline basis selection and adaptive algorithms were studied [41, 42]. However, as mentioned above, POD(and DEIM) needs a number of FOM simulations to construct the ROM. For this reason, [43] implementeda Multi-Fidelity strategy in which the parametric dependence was reconstructed using a large number oflow-fidelity models and a minimal number of high-fidelity evaluations. Other approaches exploit machinelearning to construct an input-output relationship, with convolutional neural networks [44] and autoencoders[45], which require the training of a network, again, using preexisting data. Note that most of the abovemethods lead to pROMs that are only evaluated in the online phase, i.e. no simulation is actually performed ,but the solutions at the known parameter locations are “interpolated” to obtain the result.Although model-driven NL-pROMs seem to be less popular, they offer the undeniable advantage of beingsimulation-free, thus considerably cutting down the offline costs. Interesting recent examples are looselybased on the extension of methods for linear systems, such as the Non-Linear Moment Matching (NLMM)scheme [46, 47, 48]. In Ref. [49], a non-parametric ROM is constructed with NLMM and DEIM for eachparameter instance sampled from the parameter space. These models are then “adjusted” onto a commonsubspace where they are interpolated to produce the pROM. This strategy, however, requires the solution ofa set of nonlinear algebraic equations on the FOM at different time instances, for different signal generators,and at each point on the parameter grid. For large systems, the computational effort could still be significant,although lower than that of POD methods.In this paper we propose a NL-pROM for geometric nonlinearities and parametrized shape defects tostudy the behavior of imperfect structures. This is motivated by the fact that, as it is observed in manyengineering applications, even small imperfections can significantly change characteristics and performances By simulation we refer to the solution of a set of equations describing a system in any kind of analysis setting (e.g. in timeor frequency domain).
3f a system, as for instance in the case of MEMS devices [50, 51] and mistuning of gas turbine blades [14].Other ROMs have already been developed in this sense [52, 53], but limited to localized defects. Regardinggeometric nonlinearities, we recall that in the case of continuum finite elements with linear elastic constitutivelaw and Total Lagrangian formulation, as in our study, the nonlinear elastic forces are a polynomials whichcan be represented using tensors, so that qualitatively the FOM governing equations write M ¨ u F + C d ˙ u F + K F u F + K F : ( u F ⊗ u F ) + K F ... ( u F ⊗ u F ⊗ u F ) = f ext ( t ) (1)where M , C d ∈ R n × n are the mass and damping matrices, u F , ˙ u F , ¨ u F ∈ R n the displacement, velocity andacceleration vectors, and f ext ( t ) ∈ R n an external forcing, being n the FOM number of dofs. K F ∈ R n × n , K F ∈ R n × n × n and K F ∈ R n × n × n × n are the stiffness tensors for the linear, quadratic and cubic elasticinternal forces.Conceptually, the method retrace the one we presented in [1], but it is based on a different deformationscheme (of which our earlier work resulted to be a sub-case). An overview of the individual steps of themethod is shown in Fig. 1. The user defines as input data the nominal structure (in terms of geometry,material properties and FE mesh) and a number m of displacement fields representing the shape defects , whichare intended as small deviations from the nominal geometry (Fig. 1a). These can be known analytically, fromexperimental measurements or previous simulations, and finally they can be discretized with displacementfield vectors U i and collected in a matrix U = [ U , ..., U m ]. Each defect can then be leveraged in amplitude bythe parameter vector ξ = [ ξ , ..., ξ m ] T (Fig. 1b) so that the final defected geometry represented by the modelis given by the global defect displacement field u d = U ξ , i.e. a linear superposition of the selected defects(Fig. 1c). With this information about the nominal structure and shape defects, we assemble the RB usinga modal approach with VMs, MDs and Defect Sensitivities (DSs). We then construct the reduced stiffnesstensors, once and for all , projecting the element-level tensors with the selected RB (Fig. 1d). In this way,linear, quadratic and cubic elastic forces can be evaluated directly with respect to the reduced coordinates and shape defect magnitudes without switching between the full and reduced order space when evaluatingthe nonlinear function. Our strategy can then be classified as model-driven (simulation-free). Finally, in theonline phase, the simulation is performed with the reduced governing equations (Fig. 1e). Notice that themodel is used to run a simulation, not to evaluate a solution as in interpolation-like techniques: as such,different forcing terms and also different analysis types (e.g. transient, frequency response) are possible.All of this is possible thanks to the modified definition of the Green-Lagrange strain tensor we use.Specifically, our strain tensor embeds two subsequent transformations: (i) the one from nominal to defectedgeometry (which, at the end, will be parametrized), and (ii) the one from the defected configuration to thedeformed/final one. The deformation produced by the latter is the one we measure, so no strain/stresses Due to memory limitations, third and fourth order stiffness tensors cannot be computed for the FOM, but they can beconstructed in reduced form directly operating at element level [27]. ⊗ denotes the outer product, : and ... the double and triple contraction operations. x (a) (b) (c)(d) (e) Figure 1: Overview of the proposed method, schematically illustrated for a pinned beam. are introduced by the presence of the defect in (i); however, the deformation of (ii) will depend on (i). Theformulation we obtain however contains rational terms which cannot be used for a tensorial representation(which can describe polynomials only). Given the assumed small entity of the shape defects, we advocatethe use of a Neumann expansion to approximate the Green-Lagrange tensor, obtaining again a polynomialform. Applying standard FE procedures, we finally get to the expression of the reduced elastic internal forces,which will parametrically depend on the defect amplitudes ξ . In this framework, we show that the modelin [1] (whose deformation formulation was based on [54]) corresponds to a lower order Neumann expansionwith integrals evaluated on the nominal volume, and that the higher order approximation we propose hereleads to better accuracy and to a larger applicability range.The work is organized as follows: the modified strain formulation is given in Section 2 and approximatedusing Neumann expansion in Section 3. In Section 4 the FE discretization is developed and then used inSection 5 to construct the reduced order stiffness tensors. The choice and computation of the RB is describedin Section 6. Finally, numerical studies in Sections 7 and 8 demonstrate the effectiveness of the proposedapproach on a 2D FE clamped beam and on a MEMS gyroscope and computational times are discussed.
2. Strain formulation: a two-steps deformation approach
Let us consider the scheme depicted in Fig. 2. A nominal body of coordinates x = { x , y , z } undergoesa first deformation described by the map Φ , which brings the body in a new configuration with coordinates x d = { x d , y d , z d } = Φ ( x ). The displacement corresponding to this operation is u d = { u d , v d , w d } = x d − x . We will refer to this second configuration as the defected configuration . As it will be detailed later,in our method u d will be a user-defined displacement field representing a shape defect which, superimposedto the nominal geometry, defines the configuration with respect to which we measure deformation. Let us5 igure 2: Scheme for the considered deformation setting. A nominal structure, of coordinates x , undergoes a deformationdescribed by the transformation map Φ . The structure is now in the deformed configuration (coordinates x d ). A secondtransformation Φ and the displacement u describe the deformation from the defected configuration to the final one. now consider a second deformation, described by the map Φ , from the defected configuration to the finalone, with coordinates x ( t ) = { x ( t ) , y ( t ) , z ( t ) } = Φ ( x d , t ). We will refer to the latter as to the deformed orfinal configuration , whose displacement is given by u = { u, v, w } = x − x d .Considering the infinitesimal line segment d x in the nominal geometry, we can define the line segmentsd x d and d x in the defected and deformed configurations asd x d = F d x , (2a)d x = F d x d = F F d x , (2b)where the deformation gradients F and F are given by F = ∇ x d = ∂ x d ∂ x = I + ∇ u d = I + ∂ u d ∂ x = I + D d , (3a) F = ∇ d x = ∂ x ∂ x d = I + ∇ d u = I + ∂ u ∂ x d = I + D . (3b)and where D d and D are the displacement derivative matrices of the first and second transformations,respectively. Using the chain rule, we can also define D = ∂ u ∂ x = ∂ u ∂ x d ∂ x d ∂ x = D F , (4)so that D = DF − can be referred to the nominal coordinates.Using Eqs. (2)–(4), the stretch between deformed and defected configurations writes S = d x T d x − d x Td d x d = d x T F T ( F T F − I ) F d x = d x T ( D + D T + D T D + D Td D + D T D d )d x . (5)6easuring the deformation with respect to the defected configuration, the second order Green-Lagrangestrain tensor E is defined as S = 2d x Td E d x d = 2d x T F T E F d x , (6)which, rearranged, leads to E = 12 F − T ( D + D T + D T D + D Td D + D T D d ) F − . (7)Looking at Eqs. (5) and (7), it can be easily verified that E correctly satisfies the minimum requirementsfor a strain measure to vanish under a rigid body translation ( F = I ) and/or rotation ( F T F = R T R = I ,with R orthonormal rotation matrix), for any F . Eq. (7) is indeed an exact expression for the strains fromdefected to final configuration. Notice however that in this form all the quantities are computed with respectto the nominal coordinates x .
3. Strain approximations
The introduced strain measure, being referred to the nominal geometry only, paves the way for the pre-computation of the stiffness tensors, as it will be shown in the following sections. However, as mentioned inthe introduction, a tensorial formulation can be applied only when the internal forces display a polynomialdependence on the displacements, which in the present case include both u d and u . The inverse of thedeformation gradient F in Eq. (7) entails a rational dependence on u d , and therefore needs some attention.Let us consider the following known result: Neumann expansion. If P is a square matrix and the Neumann series (cid:80) + ∞ n =0 P n is convergent, we have that( I − P ) − = + ∞ (cid:88) n =0 P n (8)A spectral norm ε = (cid:107) P (cid:107) < sufficient condition for the convergence of the Neumann series. Moreover,it can be shown [55] that truncating the sum to order N the norm error is bounded as δ N = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ( I − P ) − − N (cid:88) n =0 P n (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ≤ ε N +1 − ε = δ lim . (9)Letting P = − D d , we can expand F − using the Neumann series as F − = ( I + D d ) − ≈ N (cid:88) n =0 ( − D d ) n . (10) The spectral norm of a matrix A is defined as the square root of the largest singular value of A ∗ A , being A ∗ the conjugatetranspose of A , that is: (cid:107) A (cid:107) = (cid:112) λ max ( A ∗ A ) (cid:107) D d (cid:107) (cid:28)
1) the series is guaranteed to converge. Moreover,we can truncate the sum in Eq. (10) to N = 1, obtaining: E ,N = 12 ( I − D d ) T ( D + D T + D T D + D Td D + D T D d )( I − D d ) , (11)which, solving the product, can be rewritten as: E ,N = 12 (cid:104) (cid:16) D + D T + D T D + (cid:8)(cid:8)(cid:8) D Td D + (cid:24)(cid:24)(cid:24)(cid:24) D T D d (cid:17) + − (cid:16) (cid:8)(cid:8)(cid:8) D Td D + D Td D T + D Td D T D + D Td D Td D + (cid:24)(cid:24)(cid:24)(cid:24)(cid:24) D Td D T D d (cid:17) + − (cid:16) DD d + (cid:24)(cid:24)(cid:24)(cid:24) D T D d + D T DD d + (cid:24)(cid:24)(cid:24)(cid:24) D Td DD d + D T D d D d (cid:17) ++ (cid:16) (cid:24)(cid:24)(cid:24)(cid:24) D Td DD d + (cid:24)(cid:24)(cid:24)(cid:24)(cid:24) D Td D T D d + D Td D T DD d + D Td D Td DD d + D Td D T D d D d (cid:17) (cid:105) . (12)Finally, neglecting the terms O ( D d ), i.e. assuming that the first transformation Φ is linear, Eq. (12) reduces: E ,N = 12 (cid:16) D + D T + D T D − D Td D T − DD d − D Td D T D − D T DD d (cid:17) (13)The modified Green-Lagrange strain tensor E ,N is a polynomial function of the derivatives of the displace-ment fields u and u d , and can be thus used to compute a ROM using tensors. Remark 1 (on Budiansky approximation).
The strain formulation in [54], used by Budiansky to studybuckling in presence of defects, was obtained by subtracting the strain that a defect would produce onthe nominal structure from the strain of the deformed structure measured with respect to the nominalconfiguration. It can be shown that truncating the Neummann series to the zero-th order (i.e. setting N = 0,so that F − = I ) and using Eq. (7) and (10), the strain writes: E ,N = 12 (cid:16) D + D T + D T D + D Td D + D T D d (cid:17) , (14)which is the same strain tensor we adopted in [1] following Budiansky’s approximation.
4. Finite Element formulation
In this section we derive the elastic internal forces (at element-level) for the FE discretization of the fullorder model based on the strain as defined in Eq. (13). We remark that this full model represents just anapproximation of the reference full order model FOM-d (where the defect is embedded directly in the meshby shifting the position of the nodes). Although not offering any direct advantage over FOM-d, this fullmodel will allow us to compute the parametric ROM, as it will be explained in Section 5.First, it is convenient to switch to Voigt notation. Let θ = { u ,x u ,y u ,z v ,x v ,y v ,z w ,x w ,y w ,z } T be thevectorized form of D and θ d = { u d,x u d,y u d,z v d,x v d,y v d,z w d,x w d,y w d,z } T be the vectorized form of D d . For the derivatives of the displacement components, we use the notation u ,x = ∂u∂x and u d,x = ∂u d ∂x (similar definition for v , w and the other spatial coordinates). E ,N = (cid:18) H + 12 A ( θ ) + A ( θ d ) + A ( θ d ) A ( θ ) (cid:19) θ (15)where H = , A ( θ ) = u ,x v ,x w ,x u ,y v ,y w ,y
00 0 u ,z v ,z w ,z u ,y u ,x v ,y v ,x w ,y w ,x u ,z u ,x v ,z v ,x w ,z w ,x u ,z u ,y v ,z v ,y w ,z w ,y , s.t.: D + D T + D T D ←→ (2 H + A ( θ )) θ , (16) A ( θ d ) = ( − u d,x v d,x w d,x u d,y v d,y w d,y u d,z v d,z w d,z u d,y v d,y w d,y u d,x v d,x w d,x u d,z v d,z w d,z u d,x v d,x w d,x u d,z v d,z w d,z u d,y v d,y w d,y , s.t.: − D Td D T − DD d ←→ A ( θ d ) θ , (17) A ( θ d ) = (cid:18) − (cid:19) u d,x v d,x w d,x
00 2 v d,y u d,y w d,y w d,z u d,z v d,z u d,y v d,x u d,x + v d,y w d,y w d,x u d,z w d,x v d,z u d,x + w d,z v d,x v d,z w d,y u d,z u d,y v d,y + w d,z , s.t.: − D Td D T D − D T DD d ←→ A ( θ d ) A ( θ ) θ . (18)Let us now define u e and u ed as the nodal displacement vectors of a (continuum) finite element. Calling G the shape function derivatives matrix, such that θ = Gu e and θ d = Gu ed , and exploiting the property bywhich A ( θ ) δ θ = A ( δ θ ) θ , the virtual variation of the strain in Eq. (15) writes δ E ,N = ( H + A + A + 2 A A ) G δ u e = B N δ u e , (19)where B is the strain-displacements matrix and where we dropped the explicit dependencies on θ d and θ toease the notation. The virtual work of internal forces is given by W int = (cid:90) V d δ E T ,N S d V d = ( δ u e ) T (cid:90) V d B TN S d V d , (20)9here S = CE is the Piola-Kirchhoff stress in Voigt notation, being C the linear elastic constitutive matrix,and where V d is the volume of the defected configuration. The expression for the internal forces f int followsfrom the virtual work: f int = (cid:90) V d B TN CE ,N d V d . (21)Finally, the tangent stiffness matrix can be computed as usual taking the virtual variation of the internalforces (see Appendix B). Equations (15) and (21) can be used to perform tests and/or simulations ofthe full model and to compare the results to the corresponding FOM-d in order to assess the quality ofthe approximation before the reduction of the model. In the next section, the DpROM derived from thisformulation is presented.
5. Tensorial representation of internal forces
Under the hypothesis that for small defects V d ≈ V o , Eq. (21) in full can be written as f int = (cid:90) V o G T ( H + A + A + 2 A A ) T C (cid:18) H + 12 A + A + A A (cid:19) Gu e d V o . (22)To “extract” the displacement vectors u and u d from matrices A , A and A , we can write: A = L · θ = L · ( Gu e ) , (23a) A = L · θ d = L · ( Gu ed ) , (23b) A A = ( L · θ d ) · θ = ( L · ( Gu ed )) · ( Gu e ) , (23c)where L , L ∈ R × × and L ∈ R × × × are constant sparse matrices (see Appendix A).Dropping for convenience the integral operation over the volume V o (implicitly assumed), we can separatethe contributions in Eq. (22) as f = G T (cid:16) H T CH + H T CA + A T CH + A T CA (cid:17) Gu e , (24a) f = G T (cid:18) H T CA + A T CH + 12 A T CA + A T CA +2 A T A T CH + H T CA A + 2 A T A T CA + A T CA A (cid:17) Gu e , (24b) f = G T (cid:18) A T CA + A T CA A + A T A T CA + 2 A T A T CA A (cid:19) Gu e , (24c)where f , f and f are the linear, quadratic and cubic terms in the displacement u , respectively. These canbe recasted in tensorial form as f = [ K n + K d · u ed + K dd : ( u ed ⊗ u ed )] · u = K ( u ed ) · u e , (25a) f = [ K n + K d · u ed + K dd : ( u ed ⊗ u ed )] : ( u ⊗ u ) = K ( u ed ) : ( u e ⊗ u e ) , (25b) f = [ K n + K d · u ed + K dd : ( u ed ⊗ u ed )] ... ( u ⊗ u ⊗ u ) = K ( u ed ) ... ( u e ⊗ u e ⊗ u e ) . (25c)10 .2. Reduced tensors and internal forces We now derive the reduced internal forces and tensors via Galerkin projection. Let us assume the followingreduction for the displacement vectors u e and u ed : u e ≈ V e η , with V e ∈ R n e × m , η ∈ R m , (26a) u ed = U e ξ , with U e ∈ R n e × m d , ξ ∈ R m d (26b)being V e and U e the partitions of the RBs ( V ∈ R n × m and U ∈ R n × m d , with n number of dofs of the fullorder model) relative to the element, η and ξ the reduced coordinates. m (cid:28) n is thus the number of vectorsincluded in the RB V while m d represents the number of the assumed shape defects, collected column-wisein U . Introducing Γ = GV e and Υ = GU e , we can directly define the reduced order tensors using Einstein’snotation as Q n IJ = Γ iI H ji C jk H kl Γ kJ , (27a) Q d IJK = Γ iI (cid:0) H ji C jk L kla Υ aK + L jia Υ aK C jk H kl (cid:1) Γ lJ , (27b) Q dd IJKL = Γ iI L jia Υ aK C jk L klb Υ bL Γ lJ , (27c) Q n IJK = Γ iI (cid:18) H ji C jk L kla Γ aK + L jia Γ aK C jk H kl (cid:19) Γ lJ , (27d) Q d IJKL = Γ iI (cid:18) L jia Υ aL C jk L klb Γ bK + L jia Γ aK C jk L klb Υ bL ++ 2 L jiab Υ bL Γ aK C jk H kl + H ji C jk L klab Υ bL Γ aK (cid:17) Γ lJ , (27e) Q dd IJKLM = Γ iI (cid:0) L jiab Υ bL Γ aK C jk L klc Υ cM + L jia Υ aL C jk L klbc Υ cM Γ bK (cid:1) Γ lJ , (27f) Q n IJKL = 12 Γ iI L jia Γ aK C jk L klb Γ bL Γ lJ , (27g) Q d IJKLM = Γ iI (cid:0) L jia Γ aK C jk L klbc Υ cM Γ bL + L jiab Υ bM Γ aK C jk L klc Γ lL (cid:1) Γ lJ , (27h) Q dd IJKLMN = 2 Γ iI L jiab Υ bM Γ aL C jk L klcd Υ dN Γ cK Γ lJ . (27i)where, for convenience, tensor dimensions of size m are denoted by capital letter subscripts, dimensions ofsize m d by underlined capital letter ones. So, for example, Q d ∈ R m × m × m d . The global reduced tensors ofthe full structure can then be computed directly summing up the element contributions, a procedure whichis highly parallelizable. Reduced (global) internal forces can therefore be defined as f int,r = f r + f r + f r ,where f r = [ Q n + Q d · ξ + Q dd : ( ξ ⊗ ξ )] · η = Q ( ξ ) · η , (28a) f r = [ Q n + Q d · ξ + Q dd : ( ξ ⊗ ξ )] : ( η ⊗ η ) = Q ( ξ ) : ( η ⊗ η ) , (28b) f r = [ Q n + Q d · ξ + Q dd : ( ξ ⊗ ξ )] ... ( η ⊗ η ⊗ η ) = Q ( ξ ) ... ( η ⊗ η ⊗ η ) , (28c)and the reduced tangent stiffness matrix Q t can be written simply as Q t IJ = Q IJ + ( Q IJj + Q IjJ ) η j + ( Q IJij + Q IiJj + Q IijJ ) η i η j . (29)11inally, the equations of motion for the reduced system write: M r ¨ η ( t ) + C r ˙ η ( t ) + f int,r ( η ( t ) , ξ ) = f ext,r ( t ) (30)where M r = V T MV and C r = V T C d V are the reduced mass and damping matrices, f ext,r ( t ) = V T f ext ( t )the reduced external forces acting on the system. Notice that, in accordance to the hypothesis made for theinternal forces, also these matrices must be integrated over the nominal volume V o . Remark 2 (on tensor computation) . Equations (27) give directly the stiffness tensors in reduced form, andthis is in general highly desirable as their integration over the element volume takes multiple evaluations (e.g.through Gauss quadrature). Since the computational complexity highly depends on the number of dofs ofthe tensor, it is preferable to integrate directly the reduced ones as long as the number of reduced coordinates m is lower than the number of element’s dofs n e (e.g. n e = 60 for a serendipity hexahedron with quadraticshape functions). In case m > n e , it is computationally more efficient to compute the element tensors first(using Eqs (27) and replacing both Γ and Υ with G ) for Gauss integration, then project the element tensorsusing V and U accordingly. A similar reasoning can be done for m d , but under the very likely hypothesisthat m d (cid:28) n e it results almost always convenient to adopt the reduced form. The tensors in Eqs. (28) can be computed once and for all integrating on the nominal volume V o , under thesaid hypothesis that the defect does not change the defected volume V d significantly. When this hypothesiscannot be made, one can adopt the following approximation. Let Q e be the generic expression of a tensor tobe integrated over the volume V d (element level). We can compute the final reduced tensor Q as: Q = N el (cid:88) e =1 (cid:90) V d Q e d V d = N el (cid:88) e =1 (cid:90) V o Q e det( F ) d V o . (31)where N el is the total number of elements. The determinant of F can now be approximated retaining onlyfirst order terms. To the purpose of illustration, let us consider the following 2D example where the globaldefect is given by the linear superposition of two shape defects, that is: u d ( x , ξ ) = u d v d = f (1) u ( x ) f (2) u ( x ) f (1) v ( x ) f (2) v ( x ) ξ ξ = (cid:104) f (1) , f (2) (cid:105) ξ . (32)where we denote with f ( i ) = [ f ( i ) u , f ( i ) v ] T the vector of the functions describing the i-th shape-defect for thex-displacement u d and the y-displacement v d , respectively. We can approximate the determinant of F asdet( F ) = 1 + u d,x + v d,y + u d,x v d,y − u d,y v d,x = 1 + ξ (cid:16) f (1) u,x + f (1) v,y (cid:17) + ξ (cid:16) f (2) u,x + f (2) v,y (cid:17) + ξ (cid:16) f (1) u,x f (1) v,y − f (1) u,y f (1) v,x (cid:17) ++ ξ ξ (cid:16) f (1) u,x f (2) v,y + f (2) u,x f (1) v,y − f (1) u,y f (2) v,x − f (2) u,y f (1) v,x (cid:17) + ξ (cid:16) f (2) u,x f (2) v,y − f (2) u,y f (2) v,x (cid:17) ≈ ξ (cid:16) f (1) u,x + f (1) v,y (cid:17) + ξ (cid:16) f (2) u,x + f (2) v,y (cid:17) m d defects, we can writedet( F ) ≈ m d (cid:88) i =1 ξ i (cid:16) div f ( i ) (cid:17) (33)so that Eq. (31) can be approximated as: Q ≈ N el (cid:88) e =1 (cid:90) V o Q e d V + m d (cid:88) i =1 (cid:32) ξ i N el (cid:88) e =1 (cid:90) V o Q e (cid:16) div f ( i ) (cid:17) d V o (cid:33) = Q (cid:48) + m d (cid:88) i =1 ξ i Q (cid:48)(cid:48) i (34)where Q (cid:48) is the tensor evaluated on the nominal volume and Q (cid:48)(cid:48) i is the contribution of the i-th defect, whichcan be computed again once and for all offline, referring to the nominal volume. The additional compu-tational burden to compute Q (cid:48)(cid:48) i grows less than linearly with the number of defects, since in a quadratureintegration scheme we can use the Q e evaluated at integration points (using Eqs. (27)) both for Q (cid:48) i and for Q (cid:48)(cid:48) i . The additional computations therefore involve only scalar by tensor multiplications and tensor sums, sothat most of the added computational time is merely due to memory access management. Notice that onecould also compute all the additional tensors needed to describe det( F ) with no approximation (even toughthis is in most cases unnecessary, for h.o.t. do not improve accuracy significantly). However, the first orderapproximation we presented has the advantage to introduce only one new term for every additional defect.Finally, if adopting this correction one should use the mass matrix computed over the defected volume V d in accordance with the present formulation. This way a new mass matrix M d must be computed for eachnew parameter realization, however, being the matrix constant during the analysis, this additional cost isnegligible. Remark 3 (on computational efficiency) . The corrective terms Q (cid:48)(cid:48) i in Eq. (34) are null for an isochorictransformation between nominal and defected domain (det( F )=1). In practice, one can set up a procedureto avoid the computation of these terms to speed up the construction of the reduced tensors.
6. Reduction Basis
To construct the system described so far it is necessary to select the bases V and U . The latter is simply acollection of user-defined displacement vectors, each representing one specific defect, so that the final defectedshape is given by a linear superposition (see Eq. (26b)). The (properly said) RB is V , whose choice maynot be trivial, as it must correctly represent the system response over a range of parameters without thepossibility to be changed (since a change would require to recompute the stiffness tensors). As previouslydone in [1], our choice is to use a modal-based approach including VMs, MDs and Defect-Sensitivities (DSs)[56] in the RB, as this solution offers a way to construct a basis in a direct way, that is without convolutedbasis selection strategies, the need of computing all (or an excessively high number of) eigenvectors or theneed for previously computed simulations. We remark however that, in principle, one could use also otherRBs, as long as they are valid over the parameter space.13 able 1: Acronyms for the different models considered in the numerical studies. Model Description
FOM-d
Full Order Model with defect included by shifting the mesh nodes from the nominal configuration(no approximation). It is the reference model.
ROM-d
Reduced Order Model computed from FOM-d. Its reduction basis comprises VMs and MDs.
DpROM-N0
Defect parametric Reduced Order Model, based on the 0-th order Neumann expansion (see Eq. (14)).Its reduction basis comprises VMs, MDs and DSs. Tensors are up to the 4-th order (see [1]).
DpROM-N1
Defect parametric Reduced Order Model, based on the 1-st order Neumann expansion (see Eq. (13)).Its reduction basis comprises VMs, MDs and DSs. Tensors are up to the 6-th order (see Eqs. (27)).
DpROM-N1t
Truncated version of DpROM-N1, where terms of order O ( D d D ) in Eq. (13) have been neglected.As a consequence, Q d , Q dd , Q dd and the last two terms in Eq. (27e) are null. -v (suffix) as the corresponding DpROM, but with the volume integration correction described in Section 5.3. Let us consider the following eigenvalue problem (cid:0) K t − ω i M (cid:1) Φ i = (35)where K t = K t ( u , u d ) is the tangent stiffness matrix, M the mass matrix, ω i the i-th eigenfrequency and Φ i the corresponding eigenvector. Static Modal Derivatives θ ij (MDs) are computed neglecting the mass term,by taking the derivative of Eq. (35) with respect to η j and evaluating the resulting expression at equilibrium(i.e. η j = 0) and for ξ = : θ ij = ∂ Φ i ∂η j (cid:12)(cid:12)(cid:12)(cid:12) = − K − ∂ K t ( Φ j η j , ) ∂η j (cid:12)(cid:12)(cid:12)(cid:12) Φ i . (36)Defect Sensitivities (DSs) Ξ i,j can be obtained following a similar procedure, differentiating with respectto ξ j : Ξ i,j = ∂ Φ i ∂ξ j (cid:12)(cid:12)(cid:12)(cid:12) = − K − ∂ K t ( , U j ξ j ) ∂ξ j (cid:12)(cid:12)(cid:12)(cid:12) Φ i . (37)Expressions for the tangent stiffness derivatives are given in Appendix B. Remark 4 (on higher derivatives) . Given the higher accuracy of the model, larger defect magnitudes canbe considered as compared to [1]. To fully exploit the increased applicability range, a richer RB might benecessary, reason why we here introduce second Defect-Sensitivities (DS2s) and MDs-Sensitivities (MDSs).Let us take the derivative of Eq. (36) with respect to the k-th defect amplitude ξ k . We define the MDS θ ij,k as: θ ij,k = ∂ θ ij ∂ξ k (cid:12)(cid:12)(cid:12)(cid:12) = − K − (cid:18) ∂ K t ∂ξ k (cid:12)(cid:12)(cid:12)(cid:12) θ ij + ∂ K t ∂η j ∂ξ k (cid:12)(cid:12)(cid:12)(cid:12) Φ i + ∂ K t ∂η j (cid:12)(cid:12)(cid:12)(cid:12) Ξ i,k (cid:19) . (38)Notice that θ ij,k (cid:54) = θ ji,k . In the same manner, the second Defect Sensitivities (DS2s) with respect to ξ k write: Ξ i,jk = ∂ Ξ i,j ∂ξ k (cid:12)(cid:12)(cid:12)(cid:12) = − K − (cid:18) ∂ K t ∂ξ k (cid:12)(cid:12)(cid:12)(cid:12) Ξ i,j + ∂ K t ∂ξ j ∂ξ k (cid:12)(cid:12)(cid:12)(cid:12) Φ i + ∂ K t ∂ξ j (cid:12)(cid:12)(cid:12)(cid:12) Ξ i,k (cid:19) . (39)14 igure 3: (a) Model-I: Nominal mesh and defected mesh with ξ = 1 (and a x5 scale factor). (b) Frequency Responses andbackbone curves for different defect amplitudes ξ using the Harmonic Balance method (7 harmonics) for: ROM-d ( ),DpROM-N1 ( ), DpROM-N0 ( , only backbones). Backbones have been computed also for the FOM-d ( ) using the shooting method for validation. The vertical displacement of the mid-span of the beam is shown (first harmonic coefficient of theFourier series, normalized over the beam thickness). For each plot, the detuning parameter σ is referred to the correspondingFOM-d first eigenfrequency. The bottom-right figure collects the backbone curves for comparison. It is evident that the blind inclusion of DS2s and/or MDSs in the RB would add an unacceptable numberof unknowns, especially when considering MDSs. Depending on the type of the analysis (linear/nonlinear),on the kind of the defect (i.e. affecting the linear or the nonlinear dynamics) and on the entity of the defectitself (large/small), one can decide whether to include some vectors or not. Pre-selection strategies to reducethe basis size, as the one presented in [57] and [32], are beyond the scopes of this work and are not treatedhereafter.
7. Numerical tests – I
We consider now a FE model of an aluminum beam, of length l x = 2 m , thickness t y = 50 mm andwidth w z = 0 . m , clamped at both ends. We use a 2D plain strain model, with a mesh of 80 quadrilateralelements with quadratic shape functions (630 dofs). A Rayleigh damping matrix C d = α M d + β K isintroduced by imposing a quality factor Q = Q = 100 on the first and second modes of the linear system15 able 2: Average computational times. The data of the two DpROMs, being very similar, are clustered togeter. For the FR,the times refer to the higher forcing ( F = 4 kN). ROM-d, DpROM-N1/0 and FOM-d have 20, 25 (due to the 5 additional DSs)and 630 dofs, respectively. Notice that, being the size of the FOM-d very small, no significant conclusions in terms of speedupscan be drawn from this data (refer to the next section for a detailed discussion). Model–I
ROM-d DpROM-N1/0 FOM-dHarmonic Balance (HB)
Frequency Response
649 s 673 s –
Backbone
237 s 273 s –
Shooting
Backbone
31 min 35 min 83 h 18 min
ROM construction α = 3 . , β = 6 . · − ). A nodal load F is applied to the mid-span of the beam (with F = 1 kNand F = 4 kN for the forced responses). A shape defect defined as the vertical translation of the nodes v d ( x, ξ ) = ξt y sin( π/l x x ) is imposed, deforming the nominal geometry of the straight beam into a shallowarch. Notice that this kind of defect represents an isochoric transformation (see Section 5.3).Table 1 reports the acronyms used for the models of this and the next numerical studies. For ξ = { , . , . , . , } , backbones and Frequency Responses (FR) are computed for ROM-d, DpROM-N1 andDpROM-N0, constructed using 5 VMs, 15 MDs, and 5 DSs (only for DpROMs). The Harmonic Balance (HB)method was used (with 7 harmonics) using the NLvib Matlab tool [58] (slightly modified to adapt the directuse of tensors) and our in-house Matlab FE code. To validate the results of the ROMs, the Shooting Methodis used to find the backbones of the corresponding FOM-d. Results are shown in Fig. 3. Computationswere carried out in Matlab 2020a on a local machine equipped with an Intel(R) Xeon(R) Silver 4214 [email protected] GHz and 256 GB RAM @2666 MHz. Tensors were built in a Julia sub-routine, called by the mainMatlab code, which uses the TensorOperations package [59] for the tensor contraction. At present, the tensorconstruction is implemented serially, therefore leaving space for possible future speed-ups exploiting parallelcomputing, as remarked earlier.As it can be observed, the shift from hardening to softening behavior is well captured by all the models,with a minor loss of accuracy of the DpROMs as ξ increases. In particular, DpROM-N0 shows a significantfrequency offset of the first linear eigenfrequency which remains constant throughout the backbone curve(the same happens for the FRs, but we omitted to plot them for the sake of figures clarity). The main goalof the present test was to assess the accuracy of the method verifying the results against the FOM and overa range of frequencies. However, computational times are collected in Table 2 for completeness. Run timesfor the shooting method with the (Dp)ROMs are included for comparison. These figures, however, must betaken just as an indication, first because of the difference between FOM and ROMs in terms of convergenceduring continuation (ROMs are less likely to incur into numerical artifacts) and, secondly, because speed andconvergence of this kind of analysis is highly sensitive to several parameters and finding the best combination16 a) (b) Figure 4: (a) Model–II, a MEMS gyroscope. Drive and sense direction are indicated by arrows. (b) Meshed model, with 14,920quadratic hexahedra, for a total of 87,767 nodes and 261,495 dofs. Approximate dimension are given. by trial-and-error usually leads to sub-optimal performances. Last but not least, the size of the FOM in thiscase is too low to really appreciate the savings in terms of ROM construction.
8. Numerical tests – II
The last example we present is a prototype MEMS mono-axial gyroscope, shown in Fig. 4a. The deviceconsists in a mass suspended by four S-shaped springs, connected to the ground on the bottom of the anchors.It is a monolithic piece, produced via Deep Reactive Ion Etching (DRIE), a process which removes materialfrom a plane silicon wafer to obtain the final geometry. The etching procedure is the main cause of productionshape defects of MEMS devices, as it will be detailed later. In operative conditions, the mass is kept in motionby comb finger electrodes at the natural frequency of the drive mode (i.e. a mode featuring motion mainly inthe x-direction), so that in presence of an external angular rate Ω (along the y-axis) a vertical displacement w sense arises due to Coriolis effect along the z-axis (sense). The latter is then converted into an electricalsignal through the parallel plate electrode placed on the ground below the mass, providing the measurefor the angular rate. In general, a defect or a combination of them may create a coupling between the x-and z-axes so that the drive motion generates an additional out-of-plane displacement which superimposesto the Coriolis displacement to be measured. This is usually referred to as quadrature error since, beingproportional to the drive displacement, it is in phase quadrature with the Coriolis signal, proportional tothe drive velocity. Though it is possible to tell apart the two contributions, this is highly undesirable as itrequires dedicated, over-sized electronics to accommodate the larger displacements. Ultimately, this resultsin higher power consumption. 17 a) (b) Figure 5: (a) First defect shape, U : constant wall angle α y (only one beam is shown). The colormap indicates x-displacement.(b) Second defect shape, U : tapering of the suspension beams, 3D and top views (only one beam is shown). The colormapindicates y-displacement (absolute value). The FE model is shown in Fig. 4b and describes in detail the geometry and mesh of the device, counting14,920 quadratic hexahedral elements for a total of 261,495 dofs. For the present study we selected twotypical defects occurring in the production of MEMS devices, namely the wall angle (shown in Fig. 5a) anda restriction of the cross section of the beams (Fig. 5b). The first is generated by the fact that the plasmabeam of the DRIE process might be not perfectly perpendicular to the working plane, while the second onetypically comes from an overexposure to the chemical attacks ( over-etching ). In the spirit of our method, wecan describe the global defects as the superposition of these two displacement fields (see Eq. (26b)), letting U = [ U , U ] with the associated amplitude parameter vector ξ = [ ξ , ξ ] T . The wall angle shape defect U = [ u d , v d , w d ] T is given by u d ( ξ , z ) = ξ tan( α y ) z, (40)and v d = w d = 0. The tapering of the beams U = [ u d , v d , w d ] T is defined as v d ( ξ , x, y ) = ξ sin (cid:18) πL b ( x − x off ) (cid:19) ( y − y mid ) 2 W b , (41)and u d = w d = 0, where L b and W b are the length and the width of the beam, x off an offset depending onthe location of each beam and y mid is the y-coordinate corresponding to the middle line of each beam. Toease the interpretation of the amplitude parameters, in the following ξ is reported in degrees to representthe physical wall angle coming from the product ξ tan( α y ) in Eq. (40), while ξ is reported as a percentageof the beam thickness.We compute the FR of the MEMS gyroscope using the NLvib Matlab tool and our in-house MatlabFE code. We used a reduction basis with 3 VMs, the corresponding 6 MDs and 3 DSs per defect (onlyfor the DpROMs), and we used the HB method with H = 5 harmonics (with N s = 3 H + 1 time samplesper period, the minimum number of samples by which no sampling error is introduced in the harmonics up18he the H-th order when considering polynomial nonlinearities up to the third order [60], as in our case).Given the size of the model, we take as reference the results of ROM-d, as it would be prohibitively timeand memory consuming to compute the frequency response for FOM-d. Apart from the practical issues, wejustify this choice considering on the one side the good results obtained for lower dimensional models (as theone presented in the previous section), and on the other side considering that, ultimately, our DpROMs willbe at best as good as ROM-d, which is not parametric and not approximated in its formulation.The frequency response was obtained forcing the system in the center of the suspended mass with anodal load directed along the x-direction, with amplitude F = 0 . µ N, and using a Rayleigh damping matrixwith α = 105 and β = 0. Figure 6 reports the FRs around the first eigenfrequency of the system for thex-displacement u (drive direction) and the z-displacement w (sense direction) for all the combinations of ξ = { ◦ , . ◦ , ◦ } and ξ = { , . , , . , } . For the present study, we considered also a truncated version of DpROM-N1 (named DpROM-N1t) wherewe considered negligible in Eq. (11) the strains of order O ( D d D ). This further assumption allows us todrop all the terms multiplying A or, equivalently, L , so that Q d , Q dd , Q dd and the last two termsin Eq. (27e) are null. Although introducing a new approximation, DpROM-N1t is significantly cheaper toconstruct and, as it will be shown, does not introduce any appreciable accuracy loss in our studies. For eachof the presented pROMs then, we test the same models with the volume-integration correction described inSection 5.3. We will address to these models appending the suffix “-v” to the name of the model itself (e.g.DpROM-N1t-v). To recap, the results for a total of 6 models are reported in the following: DpROM-N0(-v),DpROM-N1t(-v) and DpROM-N1(-v). Again, see Table 1 for a quick reference. Considering first the effect of the wall angle defect only, it is apparent how DpROM-N0 performancesquickly degrade as soon as the parameter ξ is increased. This can be seen both in the error on the lineareigenfrequency and especially in the overestimated w -response, approximately one order of magnitude higherthan the reference. This may be due to the fact that the S-shaped beams are specifically designed tominimize the cross-coupling between the drive (x-) and sense (z-) axes created by the wall angle, so thatthe w -response is so small (2 orders of magnitude lower than the u -response) that it cannot be accuratelycaptured by DpROM-N0. The same observations can be made for DpROM-N0-v, as the wall angle defect byitself represent an isochoric transformation. The responses of all the other tested DpROMs show instead aperfect match with the reference when ξ = 0%.If the tapering defect only is considered (i.e. with ξ = 0 ◦ ), we observe that DpROM-N0/N1t/N1 havesimilar responses, with an error on the eigenfrequency that translates the whole response by an approximatelyconstant ∆ f . The models with the volume-integration correction were then tested. If on the one hand19 a) Drive response.(b) Sense response (out of plane). Figure 6: Frequency Responses and backbone curves for different defect amplitudes ξ (shown in degrees) and ξ (reportedas percentage of the beam width) using the Harmonic Balance method (5 harmonics) for: ROM-d ( ), DpROM-N0 ( ),DpROM-N1t ( ), DpROM-N1 ( ), DpROM-N0-v ( ), DpROM-N1t-v ( ) and DpROM-N1-v ( ). The displace-ments u and w of the center of the mass is shown (first harmonic coefficient of the Fourier series). For each plot, the percentdetuning parameter σ is referred to the corresponding FOM-d first eigenfrequency. DpROM-N0-v still presents relevant errors, on the other hand DpROM-N1t-v and DpROM-N1v show veryaccurate results in the full range of the tested ξ . Such an improvement was expected, as the volume changedby this defect affects the suspension beams dimensions, to which the eigenfrequencies of the system are verysensitive.For the remaining cases, the trends observed for the parameters ξ and ξ individually mix. Notice thatlooking at some results (e.g. u -response for ξ = 1 ◦ , ξ = 0 . igure 7: Transient response of the center of the suspended mass of the gyroscope for: ROM-d ( ), DpROM-N0-v ( ),DpROM-N1t-v ( ) and DpROM-N1-v ( ). The forcing is harmonic with period T . Case with ξ = 1 ◦ , ξ = 2%. results than DpROM-N1. This is however just a coincidence, as for DpROM-N0 the first defect shifts thefirst eigenfrequency to lower frequencies while the second defect to higher frequencies, so that the two errorsin this case cancel out. Indeed, when the volume correction is used in DpROM-N0-v, only the first effect isobserved, and the frequencies are shifted to the left.In Fig. 7 we also show the transient response of the forced node for ROM-d and DpROMs-v (case with ξ = 1 ◦ , ξ = 0 . f (as it is usually the casefor MEMS gyroscopes) with a harmonic forcing, taking 100 samples per period and for a time span equalto 10 times T = 1 /f , with F = 50 µ N. The integration was carried out in Matlab with our in-house code,using a Newmark integration scheme. Looking at the responses along the three axes, we observe that thethree DpROMs yield correct results but for DpROM-N0-v along the sense z-direction ( w component). Also,considering the z-response, we can see that DpROM-N1 is slightly better DpROM-N1t, fact that was notvery visible in the FRs. Table 3 reports the average time for the FR analyses and for the construction of the different models. Tocompare in terms of time ROM-d and the DpROMs, it is convenient to consider the variable costs (T var ),i.e. the ones that have to be sustained for each new parameter realization, and the overhead costs (T oh ), i.e.the ones sustained once and for all independently from the number of realizations. In the case of ROM-dwe have that T var = T constr + T sim , being T constr the time to construct the model (i.e. RB and tensorscomputation) and T sim the time for one simulation, while T oh = 0. For ROM-d indeed, there are no commonoverhead costs, but a new model must be constructed for each new realization of the parameters. In the caseof DpROMs instead, we have that T pvar = T psim and T poh = T pconstr (we use the superscript “p” to distinguishthe parametric models from ROM-d). For the parametric models we have in fact to pay upfront the costof model construction, which is generally more expensive than the one for ROM-d, but thereafter only thesimulation cost must be sustained for each new case. The first trivial conclusion is then that there exist anumber ¯ N of parameter realization above which DpROMs become convenient, that is: N ≥ ¯ N = (cid:24) T poh T var − T pvar (cid:25) . (42)21 able 3: Average computational times for the FR with the HB method and for the construction of each ROM (comprising thetime to compute VMs, MDs, DSs and the tensors). ROM-d counts 9 dofs, while DpROMs 15 dofs. Apart from the constructioncosts, items for the DpROMs are clustered and averaged. Notice that the construction time for ROM-d is sustained for eachnew parameter realization and contributes to the variable costs. DpROM names are reported by their respective suffix. Model–II
ROM-d N0 N1t N1 N0-v N1t-v N1-vROM construction
209 s 335 s 333 s 816 s 353 s 357 s 1,063 s
HB FR
Transient analysis
Overhead cost – 335 s 333 s 816 s 353 s 357 s 1,063 s
Variable cost (FR/transient)
Speedup (FR/transient) – / – 4 . × / 410 . × For ¯ N to be positive and finite, it follows that T var > T pvar ←→ T constr + T sim > T psim . (43)From Eq. (43) it can be seen how the convenience of the parametric model over the non-parametric onedepends on the relative weight between the simulation and construction times of the latter and the simulationtime of the former, as it can be observed in Table 3 looking at the different speedups for the FR and transientanalyses.That said, it is clearly difficult to draw general and definitive conclusions on the benefits of the twosolutions, ROM-d and DpROMs, time-wise. In the experience of the authors, transient analysis offer thebest gains, as simulation speed is very high, grows almost linearly with the simulation time span, and isless sensitive to the number of dofs than other kind of analysis, as the ones requiring continuation methods.When continuation is required, one could potentially find greater benefits in using a model with a low numberof dofs, so that ROM-d could actually become the best choice. We remark however that for ROM-d we haveto take into account also the construction cost as a variable cost, and that for large FE models the solecomputation of structural eigenmodes can already take several minutes, making this cost very high.
9. Conclusions
We presented a ROM for geometric nonlinearities that can parametrically describe a shape imperfectionwith respect to the nominal (blueprint) design, named for brevity DpROM. The imperfection is given bythe superposition of user-defined defect shapes, whose amplitudes are parameters of the model and can bechanged without reconstructing the model itself. This result has been made possible thanks to a polynomial Speedups are computed considering the variable costs only, with respect to ROM-d.
Acknowledgements
The authors want to express their gratitude to the AMS division at STMicroelectronics in Cornaredo(MI, Italy), for their support in this work. 23 ppendix A. L , L and L matrices We report in Tables A.4 and A.5 the expressions for the matrices L , L and L defined in Eqs. 23. Table A.4: Elements L (1) ijk , L (2) ijk and L (3) ijkl of the sparse 3 × × L , L and of the sparse 3 × × × L ,respectively, in the 2D case (name subscripts are moved to superscripts to avoid confusion with the indexes). L (1)111 = 1, L (1)321 = 1, L (1)312 = 1, L (1)222 = 1, L (1)123 = 1, L (1)343 = 1, L (1)324 = 1, L (1)244 = 1. L (2)111 = 1, L (2)331 = 1, L (2)312 = 1, L (2)232 = 1, L (2)123 = 1, L (2)343 = 1, L (2)324 = 1, L (2)244 = 1. L (3)1111 = 1, L (3)3211 = , L (3)3121 = , L (3)1331 = 1, L (3)3431 = , L (3)3341 = , L (3)3112 = 1, L (3)2212 = , L (3)2122 = , L (3)3332 = 1, L (3)2432 = , L (3)2342 = , L (3)1213 = , L (3)1123 = , L (3)3223 = 1, L (3)1433 = , L (3)1343 = , L (3)3443 = 1, L (3)3214 = , L (3)3124 = , L (3)2224 = 1, L (3)3434 = , L (3)3344 = , L (3)2444 = 1. Table A.5: Elements L (1) ijk , L (2) ijk and L (3) ijkl of the sparse 6 × × L , L and of the sparse 6 × × × L ,respectively, in the 3D case (name subscripts are moved to superscripts to avoid confusion with the indexes). L (1)111 = 1, L (1)421 = 1, L (1)531 = 1, L (1)412 = 1, L (1)222 = 1, L (1)632 = 1, L (1)513 = 1, L (1)623 = 1, L (1)333 = 1, L (1)144 = 1, L (1)454 = 1, L (1)564 = 1, L (1)445 = 1, L (1)255 = 1, L (1)665 = 1, L (1)546 = 1, L (1)656 = 1, L (1)366 = 1, L (1)177 = 1, L (1)487 = 1, L (1)597 = 1, L (1)478 = 1, L (1)288 = 1, L (1)698 = 1, L (1)579 = 1, L (1)689 = 1, L (1)399 = 1. L (2)111 = 1, L (2)441 = 1, L (2)571 = 1, L (2)412 = 1, L (2)242 = 1, L (2)672 = 1, L (2)513 = 1, L (2)643 = 1, L (2)373 = 1, L (2)124 = 1, L (2)454 = 1, L (2)584 = 1, L (2)425 = 1, L (2)255 = 1, L (2)685 = 1, L (2)526 = 1, L (2)656 = 1, L (2)386 = 1, L (2)137 = 1, L (2)467 = 1, L (2)597 = 1, L (2)438 = 1, L (2)268 = 1, L (2)698 = 1, L (2)539 = 1, L (2)669 = 1, L (2)399 = 1. L (3)1111 = 1, L (3)4211 = , L (3)5311 = , L (3)4121 = , L (3)5131 = , L (3)1441 = 1, L (3)4541 = , L (3)5641 = , L (3)4451 = , L (3)5461 = , L (3)1771 = 1, L (3)4871 = , L (3)5971 = , L (3)4781 = , L (3)5791 = , L (3)4112 = 1, L (3)2212 = , L (3)6312 = , L (3)2122 = , L (3)6132 = , L (3)4442 = 1, L (3)2542 = , L (3)6642 = , L (3)2452 = , L (3)6462 = , L (3)4772 = 1, L (3)2872 = , L (3)6972 = , L (3)2782 = , L (3)6792 = , L (3)5113 = 1, L (3)6213 = , L (3)3313 = , L (3)6123 = , L (3)3133 = , L (3)5443 = 1, L (3)6543 = , L (3)3643 = , L (3)6453 = , L (3)3463 = , L (3)5773 = 1, L (3)6873 = , L (3)3973 = , L (3)6783 = , L (3)3793 = , L (3)1214 = , L (3)1124 = , L (3)4224 = 1, L (3)5324 = , L (3)5234 = , L (3)1544 = , L (3)1454 = , L (3)4554 = 1, L (3)5654 = , L (3)5564 = , L (3)1874 = , L (3)1784 = , L (3)4884 = 1, L (3)5984 = , L (3)5894 = , L (3)4215 = , L (3)4125 = , L (3)2225 = 1, L (3)6325 = , L (3)6235 = , L (3)4545 = , L (3)4455 = , L (3)2555 = 1, L (3)6655 = , L (3)6565 = , L (3)4875 = , L (3)4785 = , L (3)2885 = 1, L (3)6985 = , L (3)6895 = , L (3)5216 = , L (3)5126 = , L (3)6226 = 1, L (3)3326 = , L (3)3236 = , L (3)5546 = , L (3)5456 = , L (3)6556 = 1, L (3)3656 = , L (3)3566 = , L (3)5876 = , L (3)5786 = , L (3)6886 = 1, L (3)3986 = , L (3)3896 = , L (3)1317 = , L (3)4327 = , L (3)1137 = , L (3)4237 = , L (3)5337 = 1, L (3)1647 = , L (3)4657 = , L (3)1467 = , L (3)4567 = , L (3)5667 = 1, L (3)1977 = , L (3)4987 = , L (3)1797 = , L (3)4897 = , L (3)5997 = 1, L (3)4318 = , L (3)2328 = , L (3)4138 = , L (3)2238 = , L (3)6338 = 1, L (3)4648 = , L (3)2658 = , L (3)4468 = , L (3)2568 = , L (3)6668 = 1, L (3)4978 = , L (3)2988 = , L (3)4798 = , L (3)2898 = , L (3)6998 = 1, L (3)5319 = , L (3)6329 = , L (3)5139 = , L (3)6239 = , L (3)3339 = 1, L (3)5649 = , L (3)6659 = , L (3)5469 = , L (3)6569 = , L (3)3669 = 1, L (3)5979 = , L (3)6989 = , L (3)5799 = , L (3)6899 = , L (3)3999 = 1. ppendix B. Tangent stiffness matrix derivatives The virtual variation wrt u of the internal elastic forces as defined in Eq. (22) writes δ f int = (cid:90) V o (cid:20) G T ( H + A + A + 2 A A ) T C (cid:18) H + 12 A + A + A A (cid:19) G δ u ++ G T ( H + A + A + 2 A A ) T C (cid:18) δ A + A δ A (cid:19) Gu ++ G T ( δ A + 2 A δ A ) T C (cid:18) H + 12 A + A + A A (cid:19) Gu (cid:21) d V o . (B.1)Recalling that A δ θ = δ A θ , we can write δ f int = (cid:90) V o (cid:104) G T ( H + A + A + 2 A A ) T C ( H + A + A + 2 A A ) G δ u ++ G T δ A T (cid:16) I + 2 A T (cid:17) C (cid:18) H + 12 A + A + A A (cid:19) Gu (cid:124) (cid:123)(cid:122) (cid:125) N ( u , u d ) (cid:105) d V o , (B.2)where the second term on the right-hand side can be rewritten to put in evidence the displacement virtualvariation δ u as δf (cid:48)(cid:48) I = (cid:90) V o G iI L jik G kl δu l N j d V o , (B.3)where Einstein notation was used for convenience. The tangent stiffness matrix therefore writes: K t = K (cid:48) + K (cid:48)(cid:48) (B.4)where K (cid:48) = (cid:90) V o G T ( H + A + A + 2 A A ) T C ( H + A + A + 2 A A ) G d V o , (B.5) K (cid:48)(cid:48) IJ = (cid:90) V o G iI L jik G kJ N j d V o . (B.6)Substituting u = Φ i η i and u d = U j ξ j in Eq. (B.4), taking the derivative wrt either η i and/or ξ j andevaluating the resulting expressions at equilibrium and with zero defect amplitudes, as required by Eqs.(36)–(39), we can write the derivatives of K t as: ∂ K t ∂η j (cid:12)(cid:12)(cid:12)(cid:12) = (cid:90) V o (cid:104) G T (cid:16) H T CA j + A Tj CH (cid:17) G + G · [( L · G ) · ( CHGΦ i )] (cid:105) d V o , (B.7a) ∂ K t ∂ξ j (cid:12)(cid:12)(cid:12)(cid:12) = (cid:90) V o G T (cid:16) H T CA j + A Tj CH (cid:17) G d V o , (B.7b) ∂ K t ∂η j ∂ξ k (cid:12)(cid:12)(cid:12)(cid:12) = (cid:90) V o (cid:104) G T (cid:16) A Tk CA j + A Tj CA k + 2 A Tj A Tk CH + 2 H T CA k A j (cid:17) G ++ G · (cid:104) ( L · G ) · (cid:16)(cid:16) CA + 2 A T CH (cid:17) GΦ i (cid:17)(cid:105) (cid:105) d V o , (B.7c) ∂ K t ∂ξ j ∂ξ k (cid:12)(cid:12)(cid:12)(cid:12) = (cid:90) V o G T (cid:16) A Tj CA k + A Tk CA j (cid:17) G d V o , (B.7d)where, recalling that η i and ξ j are scalars, we used A ( GΦ i η i ) = A i η i and A ( GU j ξ j ) = A j ξ j (same for A ) to avoid a cumbersome notation, and where · ij denotes the contraction of the i-th dimension of the firstterm with the j-th dimension of the second term. 25 eferences [1] J. Marconi, P. Tiso, F. Braghin, A nonlinear reduced order model with parametrized shape defects,Computer Methods in Applied Mechanics and Engineering 360 (2020) 112785.[2] T. Belytschko, W. K. Liu, B. Moran, K. Elkhodary, Nonlinear finite elements for continua and structures,John wiley & sons, 2013.[3] A. Toselli, O. B. Widlund, Domain Decomposition Methods — Algorithms and Theory, volume 34 of Springer Series in Computational Mathematics , Springer Berlin Heidelberg, Berlin, Heidelberg, 2005.doi: .[4] D. D. Klerk, D. J. Rixen, S. N. Voormeeren, General Framework for Dynamic Substructuring: History,Review and Classification of Techniques, AIAA Journal 46 (2008) 1169–1181.[5] C. Farhat, F.-X. Roux, A method of finite element tearing and interconnecting and its parallel solutionalgorithm, International Journal for Numerical Methods in Engineering 32 (1991) 1205–1227.[6] R. J. Guyan, Reduction of stiffness and mass matrices, AIAA Journal 3 (1965) 380–380.[7] J. He, Z.-F. Fu, Modal Analysis, Elsevier, 2001. doi: .[8] R. Craig, M. Bampton, Coupling of Substructures for Dynamic Analyses 6 (1968) 1313–1319.[9] S. Rubin, Improved Component-Mode Representation for Structural Dynamic Analysis, AIAA Journal13 (1975) 995–1006.[10] F. Pichler, W. Witteveen, P. Fischer, Reduced-Order Modeling of Preloaded Bolted Structures inMultibody Systems by the Use of Trial Vector Derivatives, Journal of Computational and NonlinearDynamics 12 (2017) 051032.[11] B. Blockmans, T. Tamarozzi, F. Naets, W. Desmet, A nonlinear parametric model reduction methodfor efficient gear contact simulations, International Journal for Numerical Methods in Engineering 102(2015) 1162–1191.[12] M. Balajewicz, D. Amsallem, C. Farhat, Projection-based model reduction for contact problems, Inter-national Journal for Numerical Methods in Engineering 106 (2015) 644–663.[13] M. G´eradin, D. J. Rixen, A ‘nodeless’ dual superelement formulation for structural and multibodydynamics application to reduction of contact problems, International Journal for Numerical Methods inEngineering 106 (2016) 773–798.[14] S. Mehrdad Pourkiaee, S. Zucca, A Reduced Order Model for Nonlinear Dynamics of Mistuned BladedDisks with Shroud Friction Contacts, Journal of Engineering for Gas Turbines and Power 141 (2019)1–13. 2615] F. Ghavamian, P. Tiso, A. Simone, POD–DEIM model order reduction for strain softening viscoplastic-ity, Computer Methods in Applied Mechanics and Engineering 317 (2017) 458–479.[16] L. Wu, P. Tiso, K. Tatsis, E. Chatzi, F. van Keulen, A modal derivatives enhanced Rubin substructuringmethod for geometrically nonlinear multibody systems, Multibody System Dynamics 45 (2019) 57–85.[17] L. Wu, P. Tiso, F. van Keulen, Interface Reduction with Multilevel Craig–Bampton Substructuring forComponent Mode Synthesis, AIAA Journal (2018) 1–15.[18] A. Vizzaccaro, L. Salles, C. Touz´e, Comparison of nonlinear mappings for reduced-order modellingof vibrating structures: normal form theory and quadratic manifold method with modal derivatives,Nonlinear Dynamics (2020).[19] S. Jain, P. Tiso, G. Haller, Exact nonlinear model reduction for a von K´arm´an beam: Slow-fast decom-position and spectral submanifolds, Journal of Sound and Vibration 423 (2018) 195–211.[20] S. Ponsioen, S. Jain, G. Haller, Model reduction to spectral submanifolds and forced-response calculationin high-dimensional mechanical systems, Journal of Sound and Vibration 488 (2020) 115640.[21] R. Perez, G. Bartram, T. Beberniss, R. Wiebe, S. M. Spottswood, Calibration of aero-structural reducedorder models using full-field experimental measurements, Mechanical Systems and Signal Processing 86(2017) 49–65.[22] C. Touz´e, M. Vidrascu, D. Chapelle, Direct finite element computation of non-linear modal couplingcoefficients for reduced-order shell models, Computational Mechanics 54 (2014) 567–580.[23] J. J. Hollkamp, R. W. Gordon, Reduced-order models for nonlinear response prediction: Implicitcondensation and expansion, Journal of Sound and Vibration 318 (2008) 1139–1153.[24] R. J. Kuether, B. J. Deaner, J. J. Hollkamp, M. S. Allen, Evaluation of geometrically nonlinear reduced-order models with nonlinear normal modes, AIAA Journal 52 (2015) 3273–3285.[25] M. Amabili, Reduced-order models for nonlinear vibrations, based on natural modes: the case of thecircular cylindrical shell, Philosophical Transactions of the Royal Society A: Mathematical, Physicaland Engineering Sciences 371 (2013) 20120474.[26] M. P. Mignolet, A. Przekop, S. A. Rizzi, S. M. Spottswood, A review of indirect/non-intrusive reducedorder modeling of nonlinear geometric structures, Journal of Sound and Vibration 332 (2013) 2437–2460.[27] S. Jain, Model Order Reduction for Non-linear Structural Dynamics, 2015. URL: http://resolver.tudelft.nl/uuid:cb1d7058-2cfa-439a-bb2f-22a6b0e5bb2a .2728] K. Lu, Y. Jin, Y. Chen, Y. Yang, L. Hou, Z. Zhang, Z. Li, C. Fu, Review for order reduction basedon proper orthogonal decomposition and outlooks of applications in mechanical systems, MechanicalSystems and Signal Processing 123 (2019) 264–297.[29] A. K. Noor, J. M. Peterst, Reduced Basis Technique for Nonlinear Analysis of Structures, AIAA Journal18 (1980) 455–462.[30] S. R. Idelsohn, A. Cardona, A reduction method for nonlinear structural dynamic analysis, ComputerMethods in Applied Mechanics and Engineering 49 (1985) 253–279.[31] C. S. Sombroek, P. Tiso, L. Renson, G. Kerschen, Numerical computation of nonlinear normal modesin a modal derivative subspace, Computers and Structures 195 (2018) 34–46.[32] S. Jain, P. Tiso, J. B. Rutzmoser, D. J. Rixen, A quadratic manifold for model order reduction ofnonlinear structural dynamics, Computers and Structures 188 (2017) 80–94.[33] P. Benner, S. Gugercin, K. Willcox, A Survey of Projection-Based Model Reduction Methods forParametric Dynamical Systems, SIAM Review 57 (2015) 483–531.[34] U. Baur Peter Benner Bernard Haasdonk Christian Himpe Immanuel Maier Mario Ohlberger, F. Dy-namik Komplexer, Comparison of methods for parametric model order reduction of instationary prob-lems, 2017. URL: .[35] D. Xiao, F. Fang, A. G. Buchan, C. C. Pain, I. M. Navon, A. Muggeridge, Non-intrusive reduced ordermodelling of the Navier-Stokes equations, Computer Methods in Applied Mechanics and Engineering293 (2015) 522–541.[36] D. Xiao, F. Fang, C. C. Pain, I. M. Navon, A parameterized non-intrusive reduced order model anderror analysis for general time-dependent nonlinear partial differential equations and its applications,Computer Methods in Applied Mechanics and Engineering 317 (2017) 868–889.[37] D. Xiao, Error estimation of the parametric non-intrusive reduced order model using machine learning,Computer Methods in Applied Mechanics and Engineering 355 (2019) 513–534.[38] R. Zimmermann, Manifold interpolation and model reduction (2019) 1–36.[39] M. Barrault, Y. Maday, N. C. Nguyen, A. T. Patera, An ‘empirical interpolation’ method: applicationto efficient reduced-basis discretization of partial differential equations, Comptes Rendus Mathematique339 (2004) 667–672.[40] S. Chaturantabut, D. C. Sorensen, Nonlinear Model Reduction via Discrete Empirical Interpolation,SIAM Journal on Scientific Computing 32 (2010) 2737–2764.2841] P. Phalippou, S. Bouabdallah, P. Breitkopf, P. Villon, M. Zarroug, ‘On-the-fly’ snapshots selectionfor Proper Orthogonal Decomposition with application to nonlinear dynamics, Computer Methods inApplied Mechanics and Engineering 367 (2020) 113120.[42] H. Cho, S. J. Shin, H. Kim, M. Cho, Enhanced model-order reduction approach via online adaptationfor parametrized nonlinear structural problems, Computational Mechanics 65 (2020) 331–353.[43] M. Kast, M. Guo, J. S. Hesthaven, A non-intrusive multifidelity method for the reduced order modelingof nonlinear problems, Computer Methods in Applied Mechanics and Engineering 364 (2020) 112947.[44] J. S. Hesthaven, S. Ubbiali, Non-intrusive reduced order modeling of nonlinear problems using neuralnetworks, Journal of Computational Physics 363 (2018) 55–78.[45] R. Maulik, B. Lusch, P. Balaprakash, Reduced-order modeling of advection-dominated systems withrecurrent neural networks and convolutional autoencoders, arXiv (2020).[46] A. Astolfi, Model reduction by moment matching for nonlinear systems, in: 2008 47thIEEE Conference on Decision and Control, volume 2015-Febru, IEEE, 2008, pp. 4873–4878.URL: http://ieeexplore.ieee.org/document/7039956/http://ieeexplore.ieee.org/document/4738791/ . doi: .[47] A. Astolfi, Model Reduction by Moment Matching for Linear and Nonlinear Systems, IEEE Transactionson Automatic Control 55 (2010) 2321–2336.[48] T. C. Ionescu, A. Astolfi, Nonlinear moment matching-based model order reduction, IEEE Transactionson Automatic Control 61 (2016) 2837–2847.[49] D. Rafiq, M. A. Bazaz, A framework for parametric reduction in large-scale nonlinear dynamical systems,Nonlinear Dynamics 102 (2020) 1897–1908.[50] C. Acar, A. Shkel, MEMS Vibratory Gyroscopes, Springer, 2008.[51] M. Izadi, F. Braghin, D. Giannini, D. Milani, F. Resta, M. F. Brunetto, L. G. Falorni, G. Gattere,L. Guerinoni, C. Valzasina, A comprehensive model of beams’ anisoelasticity in MEMS gyroscopes,with focus on the effect of axial non-vertical etching, 5th IEEE International Symposium on InertialSensors and Systems, INERTIAL 2018 - Proceedings (2018) 1–4.[52] X. Q. Wang, G. P. Phlipot, R. A. Perez, M. P. Mignolet, Locally enhanced reduced order modelingfor the nonlinear geometric response of structures with defects, International Journal of Non-LinearMechanics 101 (2018) 1–7.[53] X. Q. Wang, P. J. O’Hara, M. P. Mignolet, J. J. Hollkamp, Reduced Order Modeling with LocalEnrichment for the Nonlinear Geometric Response of a Cracked Panel, AIAA Journal 57 (2018) 421–436. 2954] B. Budiansky, Dynamic Buckling of Elastic Structures: Criteria and Estimates, in: Proceedings ofan International Conference Held at Northwestern University, Evanston, Illinois, Pergamon Press Ltd,1967. URL: http://linkinghub.elsevier.com/retrieve/pii/B9781483198217500107 . doi: .[55] X. Wang, S. Cen, C. Li, Generalized neumann expansion and its application in stochastic finite elementmethods, Mathematical Problems in Engineering 2013 (2013).[56] A. Hay, J. Borggaard, I. Akhtar, D. Pelletier, Reduced-order models for parameter dependent geometriesbased on shape sensitivity analysis, Journal of Computational Physics 229 (2010) 1327–1352.[57] P. Tiso, Optimal second order reduction basis selection for nonlinear transient analysis, in:Proceedings of the 29th IMAC A Conference on Structural Dynamics 2011, volume 3, 2011,pp. 27–39. URL: http://link.springer.com/10.1007/978-1-4419-9299-4http://link.springer.com/10.1007/978-1-4419-9299-4{_}3 . doi: .[58] M. Krack, J. Gross, Harmonic Balance for Nonlinear Vibration Problems, 2019. doi: .[59] Jutho, getzdan, S. Lyon, M. Protter, M. P. S, Leo, J. Garrison, F. Otto, E. Saba, D. Iouchtchenko,A. Privett, A. Morley, Jutho/tensoroperations.jl: v1.1.0, 2019. URL: https://doi.org/10.5281/zenodo.3245497 . doi:10.5281/zenodo.3245497