An enhanced parametric nonlinear reduced order model for imperfect structures using Neumann expansion

We present an enhanced version of the parametric nonlinear reduced order model for shape imperfections in structural 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 fields representing shape defects, i.e. small deviations from the nominal geometry parametrized by their respective amplitudes. The linear superposition of these artificial displacements describe the defected geometry and can be embedded in the strain formulation in such a way that, in the end, nonlinear internal elastic forces can be expressed as a polynomial function of both these defect fields and the actual displacement field. This way, a tensorial representation of the internal forces can be obtained and, owning the reduction in size of the model given by a Galerkin projection, high simulation speed-ups can be achieved. We show that by adopting a rigorous deformation framework we are able to achieve better accuracy as compared to the previous work. In particular, exploiting Neumann expansion in the definition of the Green-Lagrange strain tensor, we show that our previous model is a lower order approximation with respect to the one we present now. Two numerical examples of a clamped beam and a MEMS gyroscope finally demonstrate the benefits of the method in terms of speed and increased accuracy.


Introduction
The Finite Element (FE) method has long been a fundamental analysis and design tool in many areas of science and engineering. In structural mechanics it is almost mandatory to use FE models to investigate the behavior of complex systems, which often have many geometric details that would be difficult to handle with alternative approaches, such as lumped parameter or analytical models [2]. However, large FE simulations would often require considerable computational resources and time, so in some cases designers may prefer to perform real experiments rather than numerical ones. On the one hand, this need for fast and affordable FE simulations has given rise to numerical techniques to improve computational efficiency: domain decomposition and substructuring [3,4] and FE Tearing and Interconnecting (FETI, [5]) are just a few examples. On the other hand, model order reduction methods have emerged, consisting in the construction of a Reduced Order Model (ROM), whose number of degrees of freedom (dofs) is much smaller than that of the reference Full Order Model (FOM). The use of linear ROMs also in industrial contexts is nowadays well established as the theory underlying them. Guyan reduction [6] and modal analysis [7] are two well-known examples in mechanical 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 successfully coupled with substructuring techniques in the Craig-Bampton and Rubin methods [8,9], which are available in many commercial software.
For nonlinear FE studies, where the demand for reduction is dire, many solutions have been proposed over the last decades, but none of them seems to have prevailed over the others, as each of them offers certain advantages, requires certain costs and/or targets specific problems. Overall, however, the literature is mature enough to provide the analyst with many different options in several practical applications, ranging from bolted joints [10], gears [11], contacts [12,13], friction [14] and viscoplasticity [15] to flexible multi-body dynamics 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 consider mostly projection approaches, as the one adopted in this work; alternatively, one could resort to different strategies, such as normal form theory or Spectral Submanifolds. The most recent contributions in this sense include [18] and [19,20]. In (ii), for data-driven ROMs we usually refer to ROMs constructed using previous FOM simulation data (or experimental data, [21]), as opposed to model-driven methods that rely on some intrinsic properties of the model itself for ROM construction, such as modal approaches [22,23,24,25]. As for intrusiveness, we usually denote a ROM as non-intrusive [26] if it can be used with routines and solvers of 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 nonlinear generalized forces and jacobians. Other distinctions can be made in terms of the types of nonlinearities that a given model can handle and the way nonlinear functions are evaluated [27]. All these differences ultimately affect 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 to reduce computational effort and time, a key aspect to keep in mind when choosing a method is the overhead cost to pay in the offline phase; in the case of data-driven methods, this cost can be as high as the cost associated to the solution of the FOM [12]. Generally speaking then, data-driven methods (usually based on Proper Orthogonal Decomposition, or POD, strategies [28]) are used in scenarios where the high cost associated to the data generation can be amortized: typically, this is the case of multi-query analysis. In this sense, although not as versatile and generally applicable as data-driven POD-based approaches, model-driven strategies in structural dynamics are desirable, for no FOM simulation is required a priori. Rayleigh-Ritz procedures [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 the data-driven ones, is to resort to (nonlinear) parametric ROMs, (NL-)pROMs. Also in this context, the literature on linear systems is quite well developed and consolidated. An extensive survey and comparison of these 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 disciplines over the years. Interestingly, the vast majority of nonlinear parametric model order reduction methods is data-driven, POD-based. Some recent examples include non-intrusive interpolation methods for evaluating nonlinear functions with hypersurfaces [35,36] and use of Gaussian Processes and machine learning for error evaluation 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 methods as the Discrete Empirical Interpolation Method (DEIM) [39,40] to speed up the evaluation, and in this sense online 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] implemented a Multi-Fidelity strategy in which the parametric dependence was reconstructed using a large number of low-fidelity models and a minimal number of high-fidelity evaluations. Other approaches exploit machine learning 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 above methods lead to pROMs that are only evaluated in the online phase, i.e. no simulation is actually performed 1 , 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 being simulation-free, thus considerably cutting down the offline costs. Interesting recent examples are loosely based 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 each parameter instance sampled from the parameter space. These models are then "adjusted" onto a common subspace where they are interpolated to produce the pROM. This strategy, however, requires the solution of a 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 to study the behavior of imperfect structures. This is motivated by the fact that, as it is observed in many engineering applications, even small imperfections can significantly change characteristics and performances of 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. Regarding geometric nonlinearities, we recall that in the case of continuum finite elements with linear elastic constitutive law and Total Lagrangian formulation, as in our study, the nonlinear elastic forces are a polynomials which can be represented using tensors, so that qualitatively 2 the FOM governing equations write 3 where M, C d ∈ R n×n are the mass and damping matrices, u F ,u F ,ü F ∈ R n the displacement, velocity and acceleration vectors, and f ext (t) ∈ R n an external forcing, being n the FOM number of dofs.
F ∈ R n×n×n×n are the stiffness tensors for the linear, quadratic and cubic elastic internal forces.
Conceptually, the method retrace the one we presented in [1], but it is based on a different deformation scheme (of which our earlier work resulted to be a sub-case). An overview of the individual steps of the method 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, which are intended as small deviations from the nominal geometry (Fig. 1a). These can be known analytically, from experimental measurements or previous simulations, and finally they can be discretized with displacement field vectors U i and collected in a matrix U = [U 1 , ..., U m ]. Each defect can then be leveraged in amplitude by the parameter vector ξ = [ξ 1 , ..., ξ m ] T (Fig. 1b) so that the final defected geometry represented by the model is 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 using a modal approach with VMs, MDs and Defect Sensitivities (DSs). We then construct the reduced stiffness tensors, 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 evaluating the nonlinear function. Our strategy can then be classified as model-driven (simulation-free). Finally, in the online phase, the simulation is performed with the reduced governing equations (Fig. 1e). Notice that the model 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 defected geometry (which, at the end, will be parametrized), and (ii) the one from the defected configuration to the deformed/final one. The deformation produced by the latter is the one we measure, so no strain/stresses 2 Due to memory limitations, third and fourth order stiffness tensors cannot be computed for the FOM, but they can be constructed in reduced form directly operating at element level [27]. 3 ⊗ denotes the outer product, : and . . . the double and triple contraction operations. are introduced by the presence of the defect in (i); however, the deformation of (ii) will depend on (i). The formulation 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 advocate the use of a Neumann expansion to approximate the Green-Lagrange tensor, obtaining again a polynomial form. 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 model in [1] (whose deformation formulation was based on [54]) corresponds to a lower order Neumann expansion with integrals evaluated on the nominal volume, and that the higher order approximation we propose here leads 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 approximated using Neumann expansion in Section 3. In Section 4 the FE discretization is developed and then used in Section 5 to construct the reduced order stiffness tensors. The choice and computation of the RB is described in Section 6. Finally, numerical studies in Sections 7 and 8 demonstrate the effectiveness of the proposed approach on a 2D FE clamped beam and on a MEMS gyroscope and computational times are discussed.

Strain formulation: a two-steps deformation approach
Let us consider the scheme depicted in Fig. 2. A nominal body of coordinates x 0 = {x 0 , y 0 , z 0 } undergoes a first deformation described by the map Φ 1 , which brings the body in a new configuration with coordinates 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, superimposed to the nominal geometry, defines the configuration with respect to which we measure deformation. Let us now consider a second deformation, described by the map Φ 2 , from the defected configuration to the final one, with coordinates We will refer to the latter as to the deformed or final configuration, whose displacement is given by Considering the infinitesimal line segment dx 0 in the nominal geometry, we can define the line segments dx d and dx in the defected and deformed configurations as where the deformation gradients F 1 and F 2 are given by and where D d and D 2 are the displacement derivative matrices of the first and second transformations, respectively. Using the chain rule, we can also define so that D 2 = DF −1 1 can be referred to the nominal coordinates. Using Eqs. (2)-(4), the stretch between deformed and defected configurations writes Measuring the deformation with respect to the defected configuration, the second order Green-Lagrange strain tensor E 2 is defined as which, rearranged, leads to Looking at Eqs. (5) and (7), it can be easily verified that E 2 correctly satisfies the minimum requirements for a strain measure to vanish under a rigid body translation (F 2 = I) and/or rotation (F T 2 F 2 = R T R = I, with R orthonormal rotation matrix), for any F 1 . Eq. (7) is indeed an exact expression for the strains from defected to final configuration. Notice however that in this form all the quantities are computed with respect to the nominal coordinates x 0 .

Strain approximations
The introduced strain measure, being referred to the nominal geometry only, paves the way for the precomputation of the stiffness tensors, as it will be shown in the following sections. However, as mentioned in the introduction, a tensorial formulation can be applied only when the internal forces display a polynomial dependence on the displacements, which in the present case include both u d and u. The inverse of the deformation gradient F 1 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 +∞ n=0 P n is convergent, we have that A spectral norm 4 ε = P 2 < 1 is a 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 Letting P = −D d , we can expand F −1 1 using the Neumann series as 4 The spectral norm of a matrix A 1 is defined as the square root of the largest singular value of A * 1 A 1 , being A * 1 the conjugate transpose of A 1 , that is: Under the assumption of small defects (i.e. D d 2 1) the series is guaranteed to converge. Moreover, we can truncate the sum in Eq. (10) to N = 1, obtaining: which, solving the product, can be rewritten as: Finally, neglecting the terms O(D 2 d ), i.e. assuming that the first transformation Φ 1 is linear, Eq. (12) reduces: The modified Green-Lagrange strain tensor E 2,N 1 is a polynomial function of the derivatives of the displacement 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 study buckling in presence of defects, was obtained by subtracting the strain that a defect would produce on the nominal structure from the strain of the deformed structure measured with respect to the nominal configuration. It can be shown that truncating the Neummann series to the zero-th order (i.e. setting N = 0, so that F −1 1 = I) and using Eq. (7) and (10), the strain writes: which is the same strain tensor we adopted in [1] following Budiansky's approximation.

Finite Element formulation
In this section we derive the elastic internal forces (at element-level) for the FE discretization of the full order model based on the strain as defined in Eq. (13). We remark that this full model represents just an approximation of the reference full order model FOM-d (where the defect is embedded directly in the mesh by shifting the position of the nodes). Although not offering any direct advantage over FOM-d, this full model 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.
Eq. (13) rewrites: Let us now define u e and u e d as the nodal displacement vectors of a (continuum) finite element. Calling G the shape function derivatives matrix, such that θ = Gu e and θ d = Gu e d , and exploiting the property by which A 1 (θ)δθ = A 1 (δθ)θ, the virtual variation of the strain in Eq. (15) writes where B is the strain-displacements matrix and where we dropped the explicit dependencies on θ d and θ to ease the notation. The virtual work of internal forces is given by where S 1 = CE 1 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 follows from the virtual work: Finally, the tangent stiffness matrix can be computed as usual taking the virtual variation of the internal forces (see Appendix B). Equations (15) and (21) can be used to perform tests and/or simulations of the full model and to compare the results to the corresponding FOM-d in order to assess the quality of the approximation before the reduction of the model. In the next section, the DpROM derived from this formulation is presented.

Element-level tensors
Under the hypothesis that for small defects V d ≈ V o , Eq. (21) in full can be written as To "extract" the displacement vectors u and u d from matrices A 1 , A 2 and A 3 , we can write: where L 1 , L 2 ∈ R 6×9×9 and L 3 ∈ R 6×9×9×9 are constant sparse matrices (see Appendix A).
Dropping for convenience the integral operation over the volume V o (implicitly assumed), we can separate the contributions in Eq. (22) as where f 1 , f 2 and f 3 are the linear, quadratic and cubic terms in the displacement u, respectively. These can be recasted in tensorial form as

Reduced tensors and internal forces
We now derive the reduced internal forces and tensors via Galerkin projection. Let us assume the following reduction for the displacement vectors u e and u e d : 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 full order model) relative to the element, η and ξ the reduced coordinates. m n is thus the number of vectors included in the RB V while m d represents the number of the assumed shape defects, collected column-wise in U. Introducing Γ = GV e and Υ = GU e , we can directly define the reduced order tensors using Einstein's notation as where, for convenience, tensor dimensions of size m are denoted by capital letter subscripts, dimensions of size m d by underlined capital letter ones. So, for example, Q 3d ∈ R m×m×m d . The global reduced tensors of the full structure can then be computed directly summing up the element contributions, a procedure which is highly parallelizable. Reduced (global) internal forces can therefore be defined as and the reduced tangent stiffness matrix Q t can be written simply as Finally, the equations of motion for the reduced system write: 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 the internal 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, and this 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 of the 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 quadratic shape 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 tensors using V and U accordingly. A similar reasoning can be done for m d , but under the very likely hypothesis that m d n e it results almost always convenient to adopt the reduced form.

Volume integration
The tensors in Eqs. (28) can be computed once and for all integrating on the nominal volume V o , under the said hypothesis that the defect does not change the defected volume V d significantly. When this hypothesis cannot be made, one can adopt the following approximation. Let Q e be the generic expression of a tensor to be integrated over the volume V d (element level). We can compute the final reduced tensor Q as: where N el is the total number of elements. The determinant of F 1 can now be approximated retaining only first order terms. To the purpose of illustration, let us consider the following 2D example where the global defect is given by the linear superposition of two shape defects, that is: where we denote with T the vector of the functions describing the i-th shape-defect for the x-displacement u d and the y-displacement v d , respectively. We can approximate the determinant of F 1 as v,y where we neglected higher order terms. Generalizing this result for m d defects, we can write ξ i div f (i) (33) so that Eq. (31) can be approximated as: where Q is the tensor evaluated on the nominal volume and Q i is the contribution of the i-th defect, which can be computed again once and for all offline, referring to the nominal volume. The additional computational burden to compute Q i grows less than linearly with the number of defects, since in a quadrature integration scheme we can use the Q e evaluated at integration points (using Eqs. (27)  to avoid the computation of these terms to speed up the construction of the reduced tensors.

Reduction Basis
To construct the system described so far it is necessary to select the bases V and U. The latter is simply a collection of user-defined displacement vectors, each representing one specific defect, so that the final defected shape is given by a linear superposition (see Eq. (26b)). The (properly said) RB is V, whose choice may not be trivial, as it must correctly represent the system response over a range of parameters without the possibility to be changed (since a change would require to recompute the stiffness tensors). As previously done 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 convoluted basis selection strategies, the need of computing all (or an excessively high number of) eigenvectors or the need for previously computed simulations. We remark however that, in principle, one could use also other RBs, as long as they are valid over the parameter space.
Its reduction basis comprises VMs, MDs and DSs. Tensors are up to the 4-th order (see [1]).
Its reduction basis comprises VMs, MDs and DSs. Tensors are up to the 6-th order (see Eqs. (27)). As a consequence, Q 5d , Q 5dd , Q 6dd and the last two terms in Eq. (27e) are null.

DpROM-N1t
-v (suffix) as the corresponding DpROM, but with the volume integration correction described in Section 5.3.
Let us consider the following eigenvalue problem Defect Sensitivities (DSs) Ξ i,j can be obtained following a similar procedure, differentiating with respect to ξ j : 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 can be considered as compared to [1]. To fully exploit the increased applicability range, a richer RB might be necessary, 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: Notice that θ ij,k = θ ji,k . In the same manner, the second Defect Sensitivities (DS2s) with respect to ξ k write: It is evident that the blind inclusion of DS2s and/or MDSs in the RB would add an unacceptable number of 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 defect itself (large/small), one can decide whether to include some vectors or not. Pre-selection strategies to reduce the basis size, as the one presented in [57] and [32], are beyond the scopes of this work and are not treated hereafter.

Numerical tests -I
We consider now a FE model of an aluminum beam, of length l x = 2 m, thickness t y = 50 mm and width w z = 0.2 m, clamped at both ends. We use a 2D plain strain model, with a mesh of 80 quadrilateral elements with quadratic shape functions (630 dofs). A Rayleigh damping matrix C d = αM d + βK 0 is introduced by imposing a quality factor Q 1 = Q 2 = 100 on the first and second modes of the linear system    [59] for the tensor contraction. At present, the tensor construction is implemented serially, therefore leaving space for possible future speed-ups exploiting parallel computing, 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 significant frequency 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 goal of the present test was to assess the accuracy of the method verifying the results against the FOM and over a range of frequencies. However, computational times are collected in Table 2 for completeness. Run times for the shooting method with the (Dp)ROMs are included for comparison. These figures, however, must be taken just as an indication, first because of the difference between FOM and ROMs in terms of convergence during continuation (ROMs are less likely to incur into numerical artifacts) and, secondly, because speed and convergence of this kind of analysis is highly sensitive to several parameters and finding the best combination by trial-and-error usually leads to sub-optimal performances. Last but not least, the size of the FOM in this case is too low to really appreciate the savings in terms of ROM construction.

MEMS gyroscope
The last example we present is a prototype MEMS mono-axial gyroscope, shown in Fig. 4a. The device consists 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 material from a plane silicon wafer to obtain the final geometry. The etching procedure is the main cause of production shape defects of MEMS devices, as it will be detailed later. In operative conditions, the mass is kept in motion by comb finger electrodes at the natural frequency of the drive mode (i.e. a mode featuring motion mainly in the 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 electrical signal through the parallel plate electrode placed on the ground below the mass, providing the measure for the angular rate. In general, a defect or a combination of them may create a coupling between the xand z-axes so that the drive motion generates an additional out-of-plane displacement which superimposes to the Coriolis displacement to be measured. This is usually referred to as quadrature error since, being proportional to the drive displacement, it is in phase quadrature with the Coriolis signal, proportional to the drive velocity. Though it is possible to tell apart the two contributions, this is highly undesirable as it requires dedicated, over-sized electronics to accommodate the larger displacements. Ultimately, this results in higher power consumption.

FE model, defects and simulation details
The FE model is shown in Fig. 4b and describes in detail the geometry and mesh of the device, counting 14,920 quadratic hexahedral elements for a total of 261,495 dofs. For the present study we selected two typical defects occurring in the production of MEMS devices, namely the wall angle (shown in Fig. 5a) and a restriction of the cross section of the beams (Fig. 5b). The first is generated by the fact that the plasma beam of the DRIE process might be not perfectly perpendicular to the working plane, while the second one typically comes from an overexposure to the chemical attacks (over-etching). In the spirit of our method, we can describe the global defects as the superposition of these two displacement fields (see Eq. (26b)), letting U = [U 1 , U 2 ] with the associated amplitude parameter vector ξ = [ξ 1 , ξ 2 ] T . The wall angle shape defect and v d1 = w d1 = 0. The tapering of the beams and u d2 = w d2 = 0, where L b and W b are the length and the width of the beam, x of f an offset depending on the location of each beam and y mid is the y-coordinate corresponding to the middle line of each beam. To ease the interpretation of the amplitude parameters, in the following ξ 1 is reported in degrees to represent the physical wall angle coming from the product ξ 1 tan(α y ) in Eq. (40), while ξ 2 is reported as a percentage of the beam thickness.
We compute the FR of the MEMS gyroscope using the NLvib Matlab tool and our in-house Matlab FE code. We used a reduction basis with 3 VMs, the corresponding 6 MDs and 3 DSs per defect (only for the DpROMs), and we used the HB method with H = 5 harmonics (with N s = 3H + 1 time samples per period, the minimum number of samples by which no sampling error is introduced in the harmonics up the 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 time and memory consuming to compute the frequency response for FOM-d. Apart from the practical issues, we justify this choice considering on the one side the good results obtained for lower dimensional models (as the one presented in the previous section), and on the other side considering that, ultimately, our DpROMs will be 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 a nodal load directed along the x-direction, with amplitude F = 0.4 µN, and using a Rayleigh damping matrix with α = 105 and β = 0. Figure 6 reports the FRs around the first eigenfrequency of the system for the x-displacement u (drive direction) and the z-displacement w (sense direction) for all the combinations of , 0.5%, 1%, 1.5%, 2%}.

Tested models
For the present study, we considered also a truncated version of DpROM-N1 (named DpROM-N1t) where we considered negligible in Eq. (11) the strains of order O(D d D 2 ). This further assumption allows us to drop all the terms multiplying A 3 or, equivalently, L 3 , so that Q 5d , Q 5dd , Q 6dd and the last two terms in Eq. (27e) are null. Although introducing a new approximation, DpROM-N1t is significantly cheaper to construct and, as it will be shown, does not introduce any appreciable accuracy loss in our studies. For each of the presented pROMs then, we test the same models with the volume-integration correction described in Section 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.

Results
Considering first the effect of the wall angle defect only, it is apparent how DpROM-N0 performances quickly degrade as soon as the parameter ξ 1 is increased. This can be seen both in the error on the linear eigenfrequency and especially in the overestimated w-response, approximately one order of magnitude higher than the reference. This may be due to the fact that the S-shaped beams are specifically designed to minimize the cross-coupling between the drive (x-) and sense (z-) axes created by the wall angle, so that the w-response is so small (2 orders of magnitude lower than the u-response) that it cannot be accurately  For the remaining cases, the trends observed for the parameters ξ 1 and ξ 2 individually mix. Notice that looking at some results (e.g. u-response for ξ 1 = 1 • , ξ 2 = 0.5%), it may seem that DpROM-N0 gives better results than DpROM-N1. This is however just a coincidence, as for DpROM-N0 the first defect shifts the first eigenfrequency to lower frequencies while the second defect to higher frequencies, so that the two errors in this case cancel out. Indeed, when the volume correction is used in DpROM-N0-v, only the first effect is observed, 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 = 1 • , ξ 2 = 0.5%). Each model is forced at its own first resonance frequency f 0 (as it is usually the case for MEMS gyroscopes) with a harmonic forcing, taking 100 samples per period and for a time span equal to 10 times T 0 = 1/f 0 , 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 the three 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 not very visible in the FRs. Table 3  of DpROMs instead, we have that T p var = T p sim and T p oh = T p constr (we use the superscript "p" to distinguish the parametric models from ROM-d). For the parametric models we have in fact to pay upfront the cost of model construction, which is generally more expensive than the one for ROM-d, but thereafter only the simulation cost must be sustained for each new case. The first trivial conclusion is then that there exist a numberN of parameter realization above which DpROMs become convenient, that is:

Computational times
From Eq. (43) it can be seen how the convenience of the parametric model over the non-parametric one depends on the relative weight between the simulation and construction times of the latter and the simulation time of the former, as it can be observed in Table 3 looking at the different speedups 6 for the FR and transient analyses.
That said, it is clearly difficult to draw general and definitive conclusions on the benefits of the two solutions, ROM-d and DpROMs, time-wise. In the experience of the authors, transient analysis offer the best gains, as simulation speed is very high, grows almost linearly with the simulation time span, and is less 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 number of dofs, so that ROM-d could actually become the best choice. We remark however that for ROM-d we have to take into account also the construction cost as a variable cost, and that for large FE models the sole computation of structural eigenmodes can already take several minutes, making this cost very high.

Conclusions
We presented a ROM for geometric nonlinearities that can parametrically describe a shape imperfection with respect to the nominal (blueprint) design, named for brevity DpROM. The imperfection is given by the superposition of user-defined defect shapes, whose amplitudes are parameters of the model and can be changed without reconstructing the model itself. This result has been made possible thanks to a polynomial 6 Speedups are computed considering the variable costs only, with respect to ROM-d.
representation of the internal forces resulting from a two-step deformation process (which brings the nominal geometry into the defected one and then into the deformed one) and from the approximation of the strains obtained by a Neumann expansion. The latter allowed to eliminate rational expressions under the hypothesis of small defects, so that the elastic internal forces are written as simple polynomials both with respect to the displacement field representing the defect and with respect to the actual displacement field. Using a Galerkin projection and a modal-based approach for selecting the RB, the reduced internal forces have been recast in tensorial form, where the linear, quadratic and cubic stiffness tensors are found to be functions of a parameter vector collecting the amplitudes of the defects imposed on the structure. Within this framework we tested different versions of the DpROM for different degrees of approximation. In particular, we have shown that the model we had previously developed using Budiansky's approach corresponds to the 0th-order expansion of our model, without volume integration correction (i.e. DpROM-N0). Finally, in the numerical studies we showed that the higher order approximation DpROM-N1 effectively leads to more accurate results and that for volume-changing defects a large improvement can be achieved by approximating the tensor integral over the real volume of the defective geometry (DpROM-N1-v). The truncated version DpROM-N1t(-v) was also presented, which has almost the same accuracy as its complete counterpart, but without the need to construct tensors with dimensionality higher than four. The computational costs were then critically discussed, taking into account different types of analysis. In particular, we showed that in transient studies we can usually expect very high speedups from the parametric models. In the case of FR analysis, which we used to assess the quality of the solutions over a range of frequencies as an alternative to multiple time analyses, the gains will be more contained. In this context, to reduce the dofs of both the parametric and non-parametric ROMs and make FR analysis faster and thus closer to transient analysis in terms of time and speedups, we think that an a priori selection of the RB vectors and hyperreduction strategies would actually be very beneficial, and they can actually constitute the spur for future investigation.
Appendix A. L 1 , L 2 and L 3 matrices We report in Tables A.4 and A.5 the expressions for the matrices L 1 , L 2 and L 3 defined in Eqs. 23.
ijk and L (3) ijkl of the sparse 3 × 4 × 4 matrices L 1 , L 2 and of the sparse 3 × 4 × 4 × 4 matrix L 3 , respectively, in the 2D case (name subscripts are moved to superscripts to avoid confusion with the indexes).
ijk and L ijkl of the sparse 6 × 9 × 9 matrices L 1 , L 2 and of the sparse 6 × 9 × 9 × 9 matrix L 3 , respectively, in the 3D case (name subscripts are moved to superscripts to avoid confusion with the indexes).

Appendix B. Tangent stiffness matrix derivatives
The virtual variation wrt u of the internal elastic forces as defined in Eq. (22) writes (B.1) Recalling that A 1 δθ = δA 1 θ, we can write where the second term on the right-hand side can be rewritten to put in evidence the displacement virtual variation δu as where Einstein notation was used for convenience. The tangent stiffness matrix therefore writes: Substituting u = Φ i η i and u d = U j ξ j in Eq. (B.4), taking the derivative wrt either η i and/or ξ j and evaluating 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: where, recalling that η i and ξ j are scalars, we used A 1 (GΦ i η i ) = A 1i η i and A 2 (GU j ξ j ) = A 2j ξ j (same for A 3 ) to avoid a cumbersome notation, and where · ij denotes the contraction of the i-th dimension of the first term with the j-th dimension of the second term.