Model Order Reduction based on Direct Normal Form: Application to Large Finite Element MEMS Structures Featuring Internal Resonance

Dimensionality reduction in mechanical vibratory systems poses challenges for distributed structures including geometric nonlinearities, mainly because of the lack of invariance of the linear subspaces. A reduction method based on direct normal form computation for large finite element (FE) models is here detailed. The main advantage resides in operating directly from the physical space, hence avoiding the computation of the complete eigenfunctions spectrum. Explicit solutions are given, thus enabling a fully non-intrusive version of the reduction method. The reduced dynamics is obtained from the normal form of the geometrically nonlinear mechanical problem, free of non-resonant monomials, and truncated to the selected master coordinates, thus making a direct link with the parametrisation of invariant manifolds. The method is fully expressed with a complex-valued formalism by detailing the homological equations in a systematic manner, and the link with real-valued expressions is established. A special emphasis is put on the treatment of second-order internal resonances and the specific case of a 1:2 resonance is made explicit. Finally, applications to large-scale models of Micro-Electro-Mechanical structures featuring 1:2 and 1:3 resonances are reported, along with considerations on computational efficiency.


Introduction
The rate of development of computational power has enabled the simulation of complex systems. However, computational costs for the analysis of mechanical components are still not competitive enough on an industrial scale, where first guesses on the performance of mechanical components must be obtained rapidly. This limitation led to the development of techniques aimed at reducing the computational cost of numerical models while retaining sufficient reliability, i.e. Model Order Reduction (MOR) techniques [1,2,3]. Restricting ourselves to the case of vibratory systems including large-amplitude displacements and thus exciting geometric nonlinearities, the first idea used since decades has been to project the nonlinear equations of motion onto a selected subset of the linear modes basis, see e.g. [4,5] to cite only two examples. Unfortunately, the loss of invariance of linear eigenspaces, expressed through important nonlinear coupling terms that may happen with high-frequency modes [6], makes this approach not efficient in terms of accuracy [7,8,9,2,10]. The Proper Orthogonal Decomposition (POD) method offers a gain since being able to modify slightly the orientation of the subspaces to better fit the curvatures of the nonlinear data, but it is still restricted to the use of linear orthogonal subspaces [11,12,13,14].
On the other hand, nonlinear reduction methods start by defining a nonlinear relationship between the original coordinates and those of the reduced dynamics, hence providing a more accurate treatment of the nonlinear trajectories and thus faster convergence with fewer master modes. First attempts can be traced back to the work by Rosenberg [15] who introduced the term Nonlinear Normal Mode (NNM) [16,17], which has been further developed in the works by Shaw and Pierre who first recognised that the concept of invariant manifold is key to compute accurate ROMs for nonlinear vibratory systems [7,8,18,19]. The normal form theory together with a truncation to a small subset of master modes have been proposed in [20,9,21,10]. On the mathematical point of view, invariant manifolds are investigated for a long time, see e.g. [22,23], and the methods have been written in a unified and abstract general formalism thanks to the parametrisation methods [24,25,26]. The book by Haro et al. [27] gives a clear presentation of these developments with a presentation oriented toward applications. In particular, it shows that the computations led by previous authors in vibration theory followed the guidelines of the two main parametrisation styles of invariant manifolds. While Shaw and Pierre used the graph style, the normal form style was used in [20,9,21], with a real formulation in order to better fit the classical oscillatory framework.
A remarkable advancement within the correct and formal definition and properties of nonlinear normal modes has been given by Haller et al. in a series of papers [28,29,30,31]. In the conservative framework, existence and uniqueness of the searched invariant structures are given by the Lyapunov center theorem [32], stating that under non-resonance conditions, a two-dimensional manifold densely filled with periodic orbits exists for each couple of imaginary eigenvalues. As remarked in [28], the picture is completely different for dissipative systems since the whole phase space is foliated by invariant manifolds tangent at the origin to linear subspace (see also [33,34]). In order to give a correct, accurate and unequivocal definition of the searched reduction subspaces, Haller et al. introduced the notion of spectral submanifolds (SSM) as the smoothest nonlinear continuation of a spectral subspace of the linearised system. They proved that existence and uniqueness of SSM are related to non-resonance conditions and a spectral quotient computed from the real part of the spectrum of the linearised system. In the conservative case, the manifolds are named as Lyapunov Subcenter Manifold (LSM) [35,30]. Retrospectively, the earlier computations of NNMs turned out to be computations of LSM when damping was not taken into account. On the computational point of view, Haller et al. used normal form style [27,29], and proposed an automated computational framework allowing to go to any order of asymptotic development while taking into account damping and various nonlinear terms.
Applications of methods inherited from dynamical systems theory to finite element (FE) problems has remained scarce until recently. A major motivation was that all the methods described in the last paragraph have as common starting point the mechanical system written in modal space. Unfortunately, FE methods, where meshes with millions of degrees-of-freedom (dofs) are routinely used, make this step out of reach. The direct computation of NNMs from large Finite Element Models [16,17,36] has been addressed with harmonic balance approaches or shooting procedures. Thanks to their purely computational nature, these techniques are very general and can be applied to a wide range of problems, but they require a certain amount of computational effort to be applied to realistic structures and fail at serving the aim of generating agile reduced order models (ROMs). On the other hand, for FE structures including geometric nonlinearities, a number of methods have been proposed in the last 20 years in order to formulate ROMs that can be easily computed, in the best case in a non-intrusive manner, i.e. without the need to enter new calculations at the elementary level in the code, such that any research-oriented or commercial FE code could be used as starting point. The stiffness evaluation procedure (STEP) uses the linear modes as projection basis [37,38], but can be simply used as a technique to compute non-intrusively linear and nonlinear characteristics. The implicit condensation and expansion (ICE) method relies on a set of applied forces to keep trace of non-resonant couplings among modes [39,40,41,42]. Finally modal derivatives (MD) have been introduced with the aim of taking the amplitude dependence of modes into account, thus pursuing the same goal as the NNMs [43,44,45,46]. ICE and MD assume the manifold to be velocity independent, assumption which is the more fulfilled, the larger the slow/fast separation between the slave and master coordinates ( [47,48,49]). In [48], it is estimated that a ratio between eigenfrequencies of the slave modes and those of the master modes of at least 3 ensures the correctness of MDs in the prediction of the hardening/softening behaviour ( [50,51,49]). The inclusion of velocity dependence in the MD approach have been proposed in [52] to overcome this limitation, and leads to similar formulations than those reported in [53].
Consequently there is still a real need of direct applications of invariant-based reduction methods, using the general theorems from dynamical systems theory, and that could be applied non-intrusively on the dofs of FEM problems in physical space. Pursuing previous developments using normal form theory [9,21], explicit and non-intrusive real formulas have been proposed in [53]. All the coefficients of second and third-order nonlinear mappings have been derived by rewriting the normal form approach proposed in [9,21] directly in the physical space. A real formalism was used so that all along the calculation, the results are expressed in the form of oscillator-like equations. Damping and forcing were taken into account and illustrative examples on blades and beams were reported in [53,50]. At the same time, the computation based on SSM led by the group of Haller also proposed a direct computation that has been written in the open code SSM tool 2.0 [54].
The present article aims at completely rewriting the direct normal form approach proposed in [53] in order to explain more clearly its settings and possible further developments. Secondly, the special case of second-order internal resonance, that had not been treated in [9,21,53], is here detailed. Finally, applications to large-scale Micro-Electro-Mechanical systems (MEMS) are derived in order to show the potentiality of the method to deal with millions of dofs yet providing fast and accurate ROMs. MEMS structures are generally actuated at resonance, they are subjected to geometric nonlinearities due to large transformations and have very small damping values as they operate in near-vacuum packages, hence showing highly nonlinear dynamical features that are rarely observed at the macro scale [55,56,57,58,59,60]. Furthermore, nonlinear dynamic properties of MEMS can be tailored to yield performance that would not be accessible through operation in the linear regime [61,62,63,64]. A remarkable example of successful application of nonlinear mode interaction for the development of highly efficient mechanical filters is for instance reported in [65]. This implies that the proposed method would play a major role in this field since it would ensure a fast and efficient estimation of the frequency response functions of structures within times that are compatible with industrial design requirements.
The paper is organised as follows. Section 2 details the reduction method based on direct normal form (DNF). A first-order (state-space) formulation along with a complex-valued formalism is introduced in Sections 2.1 and 2.2, with the general idea of rewriting the DNF in a more complete and symmetric formalism, a necessary step for further developments to include more easily new forces and/or going to higher-orders. Section 2.3 details the second and third-order homological equations, written from the physical space, such that the method can be easily understood and adopted by a wider group of researchers. These steps are important novelties as compared to the calculations presented in [53] where only the translations of the generic method presented in [9,21] for direct application in physical space was given, in a real formalism. Once the complete normal form is computed, the reduction method is explained in Section 3. More insights in the computations are given for the simple case of single master coordinate in Section 3.2. The case of second-order internal resonance is then tackled in Section 3.3. Finally, Section 4 explains how one can pass from the complex formalism detailed here to the real nonlinear mappings provided in [53], hence bridging earlier works.
Section 5 shows the numerical results obtained for three different MEMS structures. First, a MEMS micromirror undergoing large rotations is considered. In this case a single master mode is needed to correctly retrieve the hardening behaviour up to very large angles of rotations. This example is important since other reduction methods had been tested before, all of them giving incorrect results, as illustrated with application of the ICE method. Then a MEMS beam resonator that features a 1:3 internal resonance [66], and a MEMS arch resonator showing 1:2 internal resonance are selected in order to show the ability of the method to deal with systems that feature internally resonant modes. All numerical simulations in the present work are compared with full-order harmonic balance finite element simulations of the systems [67,68,69].

Equations of motion and direct normal form approach
The aim of this Section is to derive a general formalism to write direct expressions for the computation of the normal transform on the physical dofs. In particular, homological equations adapted to the framework of mechanical systems including geometric nonlinearities are given.

Equations of motion and linear decomposition
Let us consider the equation of motion associated with an undamped mechanical system subjected to geometric nonlinearities: where M denotes the mass matrix, U the vector of nodal displacements of dimension N,( ·) the time derivative operator, K the stiffness matrix, G(U, U) and H(U, U, U) the quadratic and cubic force terms, respectively. This starting point is common to any finite element discretisation of the linear momentum equation for mechanical systems. In particular in the framework of three-dimensional linear elasticity with large transformations, Eq. (1) is exact. For the sake of completeness, Appendix A gives some details on the adopted notations, and Appendix B recalls how the quadratic and cubic terms are obtained from the decomposition of the internal force vector. The solution of a generalised eigenvalue problem formulated on the linear part of Eq. (1) yields a finite set of real-valued mass-normalised eigenmodes φ s which are collected column-wise in the eigenvector matrix Φ. The following equalities hold: where Ω is a diagonal matrix that stores the eigenfrequencies ω s . Equation (1) can be mapped to modal coordinates u through the linear transformation U = Φu, leading to: with: Modal decomposition of Eq. (3) is computationally infeasible for large systems and must be avoided. In order to derive a general formalism for computing the normal transform of Eq. (1), let us first express Eq. (1) in state-space formalism, i.e. as a first-order dynamical system by introducing the velocity V as the time derivative of U: Note that the previous developments on the DNF approach proposed in [53] aimed at avoiding this first-order formulation in order to keep real expressions and oscillator-like equations, with second-order derivatives in time, throughout the calculations. This could be realised at the cost of an important loss in the symmetries of the problem. Moreover, generalisations of the second-order formalism to other linear and nonlinear forces, e.g. linear viscous damping, electrostatic and piezoelectric forces, become more difficult to handle compared to a first-order approach. Consequently, we propose in this article a complete rewriting of the DNF approach, using state-space formulation as starting point, but keeping in all calculations the link to the usual mechanical representation as given in Eq. (1), so that further generalisations will be easier to derive. For the sake of conciseness, the diagonalisation of the linear part of Eq. (5) is reported in Appendix C and leads to: with p the generalised coordinates, Λ a diagonal matrix, and Γ (p, p), ∆(p, p, p) nonlinear operators in generalised coordinates. In index notation, it reads: with λ s the diagonal entries of Λ, and Γ skl , ∆ sklm the scalar coupling coefficients. Note that due to the statespace formulation, k, l, m are summed from 1 to 2N in Eq. (7), and the system is now of size 2N. Further details on the derivation of Eq. (6) are reported in Appendix C. In particular, the structure of the diagonal Λ is such that: and the λ s coefficients read: with i the imaginary unit. Moreover displacements U, velocities V, modal displacements u = [u 1 , ..., u N ] T , and modal velocities v = [v 1 , ..., v N ] T are related to the generalised coordinates via the following relationships: As highlighted by Eq. (6), the diagonalisation of the linear part of Eq. (5) through a linear change of coordinates does not guarantee a decoupling of the nonlinear terms. This implies that a given couple p s and p s+N does not define an invariant subspace of the system. The next steps of the developments aim at deriving a nonlinear mapping that could express the dynamics in an invariant-based span of the phase space, such that once the nonlinear mapping is computed with correct truncation, reduced-order models could be easily obtained by keeping only a few master coordinates. This will be realised thanks to the normal form computation, directly applied to Eq. (5).

Direct Normal Form Setting
The normal form approach has been first introduced by Poincaré, leading to well-known theorems [70,71]. In the context of vibratory systems, it has been introduced for model order reduction purpose in [9,21,10], with the additional idea of truncating to a few subset of master coordinates once the full nonlinear mapping is computed, hence retrieving the parametrisation method of invariant manifold with normal form style [27,28], as will be further discussed in Sec. 3. In this contribution, we follow the guidelines of the normal transform by first deriving the complete mapping for all coordinates, and then truncating to obtain a ROM. Consequently a nonlinear relationship between the original (U, V) variables (displacement and velocity N-dimensional vectors) and the normal z coordinate, a 2N-dimensional vector describing the dynamics in an invariant-based span, is introduced as: with Ψ (z), Υ (z) nonlinear polynomial mappings. Mappings are expanded in terms of their components and written in either tensorial form: where the in-parenthesis upperscript refers to the order of the polynomial, or in indicial form: where a bold capital letter with one or more indices denotes a vector. Since modal coordinates define invariant subspaces in the absence of nonlinearities, the first-order maps should correspond to the usual linear decomposition. This is also in line with the general idea of finding a continuation of the linear mode subspace where the higher-order terms will account for the amplitude-dependence of modal quantities, given by the curvatures of the invariant manifolds [72,10,73]. Note that this is a common idea with the quadratic manifold approach including MD as proposed in [45,46], the only difference being that in the present derivation, velocities are fully taken into account. The linear terms of the mapping are thus simply expressed as: This choice makes first-order terms identical to the linear transformation that maps U and V to the generalised coordinates p (see Appendix C). In order to derive complete expressions for the homological equations [74,75,76], one needs to first express the dynamics of the normal variable z, i.e. the normal form of the initial problem, as:ż where f (z) is expressed as a polynomial function of z with unknown coefficients: or in index form: This first-order form for the dynamics of the normal coordinates is in full accordance with [27,28,29]. Since Ψ (1) s = φ s , then f (1) is taken equal to Λ, or equivalently f (1) s = λ s . This is the logical consequence of using an identity-tangent change of coordinates that leaves the linear diagonalised part of the dynamics unchanged. The homological equations are obtained by discarding the dependence on time and equating to zero the collections of terms with same powers in the asymptotic developments, order by order. To that purpose, the time derivatives of the mappings are obtained from Eqs. (11) and (15) as: where some properties defined in Appendix A have been used. Alternatively, using indicial notation: The same operation can be applied for system nonlinearities in order to correctly distinguish each order appearing in the expansions: As before the same equation in full indicial notation is also provided for clarity: where expansions of nonlinear terms up to fourth-order is done to better highlight the nested structure of the derived equations. It is worth stressing that the coefficients of the nonlinear mapping Ψ (z), Υ (z) are unknown at this stage, together with the expression of the normal form f (z). The aim is to find the simplest possible vector field f upon nonlinear transforms and this is achieved step-by-step by writing for each order the homological equations that collect the same powers of the normal variable z.

Homological Equations
Computation of Ψ (z), Υ (z), and f (z) is performed by substituting Eqs. (11) and (18) in Eq. (5) and by collecting terms of equal power in z in hierarchical order, hence providing for each order the corresponding homological equation [74,27]. The first-order homological equation reads: The 2N equations have 4N unknowns, that are the 2N coefficients of f (1) and N vector pairs (Υ (1) , Ψ (1) ). Therefore, the solution is not unique and any triplet (Υ (1) , Ψ (1) , f (1) ) that satisfies Eq. (22) is a valid choice. The second part of Eqs. (22) allows writing showing that a direct relationship exists between the two parts of the mapping (displacement/velocity). Using the identity-tangent solution as justified before leads to the choice f (1) = Λ such that in indicial notation one can simply write Υ which makes appear the usual linear eigenproblem in vibration theory, thus justifying again the adopted choices Collecting second-order terms yields the second-order homological equation written in a manner that explicitly encompasses the usual mechanical context as: Three remarks on this second-order homological equation can be made. First, one can observe that Ψ (2) and Υ (2) are expressed in terms of the lower order terms, a feature that is typical of asymptotic developments. Second, the quadratic part of the nonlinear internal force vector G now directly appears at the right-hand side of Eq. (25), underlining that the formalism can be extended in order to include different nonlinear forces, with more complex features like velocity-dependence or even more involved physics with new variables, e.g. electrostatic or piezo-electric couplings. Even though this is not addressed in the present contribution, we believe that the framework developed in this work is sufficiently general so that inclusions of new forces can be now simply incorporated by modifying the right-hand side. The last comment regards the dependence of the velocity mapping Υ on the displacement mapping Ψ . As in Eq. (22), the second part of (25) can be extracted to show that a linear relation between Υ (2) and Ψ (n) , f (n) with n ≤ 2, exists. This implies that in the present context of conservative vibratory systems, the velocity mapping Υ is not independent of the displacement mapping Ψ , thus allowing us to keep only Ψ at all orders, and then derive Υ thanks to the known relationships, only in case this mapping is needed. At second-order, the relationship between velocity and displacement mappings read: On the other hand the first lines of Eq. (25) can be written in terms of the displacement mapping Ψ only, as: As for Eq. (22), the system is underdetermined, i.e. the number of unknowns (Ψ (2) , Υ (2) , f (2) ) is larger that the number of equations. The aim of the normal form is to simplify as much as possible the resulting dynamical system, i.e. to choose the solution such that f has the smallest number of terms, and is, in the ideal case, simply linear. Retaining this choice to solve Eq. (27) leads to selecting f (2) skl = 0, ∀ s, k, l = 1, ..., 2N whenever it is possible, and thus to express Ψ (2) kl as: When Eq. (28) is solvable, the normal form of the system will be free of quadratic coupling terms. The condition for which (λ k + λ l ) 2 M + K is singular is related to the usual second-order resonances. As already remarked in previous developments with normal transforms [9], since the spectrum is composed of purely imaginary conjugate pairs, these resonances occur if ω 2 j = (±ω k ± ω l ) 2 . For the rest of the development presented herein, we now assume that no second-order resonance between the eigenfrequencies are present. The treatment of second-order resonances is further investigated in Sec. 3.3 to enlarge the scope of the present developments and explicitly show how to take them into account.
Since the quadratic terms can be cancelled under this assumption, and since mechanical systems are known to show nonlinear response at finite amplitudes (e.g. hardening/softening behaviour), it is important to derive at least the third-order dynamics. The third-order homological equation reads: which is reported in an extended form to highlight its hierarchical structure. Again, the second part of Eq. (29) provides a direct relationship between the velocity mapping Υ and its counterpart for displacements Ψ . This result is in fact true for any homological equation in direct form of a given order: every velocity mapping Υ (n) is a linear function of displacement mappings Ψ (m) with m ≤ n. This is in accordance with the results provided by the real-valued mapping introduced in [53], where equivalent findings were reported yet not systematically proven.
Eq. (29) can be simplified under the assumption of no second-order resonances. In that case, one has f (2) = 0, such that a compact and simplified third-order homological equation can be expressed as: One can notice in particular that Eq. (30b) gives the explicit relationship between velocity and displacement mappings, while Eq. (30a) has been rewritten as a function of displacement mapping only, i.e. by using Eq. (30b) to eliminate Υ terms. Let us define Ξ (3) klm as minus the right hand side of Eq. (30a) to improve readability of following equations: From Eq. (30), third order nonlinearities can be set to zero (f can be computed as: Contrary to Eq. (28), there are always (k, l, m) combinations such that the resulting system is singular regardless of the associated eigenfrequency, because of the purely imaginary complex conjugate spectrum. This is for instance observed for (k, l, m) combinations of the type (s, s, mod(s + N, 2N)) with mod(·, ·) modulo operation. Resonance conditions that do not depend on the eigenfrequency value are called trivial resonances. Their presence and consequences for the development of reduced models is detailed in Sec. 3, see also [10] for general considerations.
In the present work, application to large-scale models are reported by truncating the mappings at secondorder, i.e. (Ψ (n) , Υ (n) ) = 0, ∀ n ≥ 3. Dynamical coefficients are truncated at third-order, i.e. f (n) = 0, ∀ n ≥ 4. These assumptions lead to the second-order DNF method as introduced in [53], which has already proven to give accurate results. Extensions to higher orders are possible, see e.g. [29]; however for the present study we restrict ourselves to this second-order DNF in Section 5 where numerical examples are presented.
The complete normal form of the problem can be written by expressing the third-order coefficients f (3) . Due to the complexity of the calculations, it is easier to present them by using the modal basis. Under the hypothesis of no second-order resonances, i.e. f (2) = 0, the general third-order dynamics can be written by projecting Eq. (30) onto the modal basis, hence yielding: from which the following relationships are retrieved, thus providing explicit expressions of f sklm : ∀ k, l, m = 1, ..., 2N, ∀ s = 1, ..., N, One can note in particular that the previous equations can be used to write the complete normal form of the problem in the 2N-dimensional setting, and thus they involve the full matrix of eigenvectors Φ. However these operations are never computed in a reduction perspective, as it will be explained in Section 3. Indeed, reduction will be performed by selecting a few master coordinates so that Φ will be restricted to the selected master eigenvectors only. Examples will be provided with direct reduction to single-mode and multi-mode dynamics.
Following the procedure developed within the present section, the method of normal transform can be expanded up to arbitrary order, thus allowing algorithmic implementations, such as those reported in [29]. A brief summary of higher order expansions, highlighting the hierarchical structure of the procedure, is reported in Appendix D even if, as recalled, in the present work, only the second-order DNF is used as it proves highly efficient for the applications addressed in Section 5.
Finally, it is worth stressing that the generic structure of the method can be summarised in three steps: for each order n of the expansion, the homological equation of order n is written, thus providing mappings Ψ (n) , Υ (n) and dynamics coefficients f (n) that cannot be set to zero. Then, f (n+1) needs to be computed through projection of the homological equation of order n+1, since the order n terms will create new monomials at order n + 1 that need to be tracked properly. The procedure then starts again by identifying mappings and non zero terms at the next order. Implementation details of the method up to second order are reported in Appendix E, and the case of second-order internal resonance will be specifically detailed in Section 3.3, underlining how the different steps of the method may create higher-order terms in the asymptotics.

Model Order Reduction
In this section, the reduction method is explained in its general settings, then two particular cases are investigated for the sake of clarity. First, reduction to a single master mode is detailed in Section 3.2, while the case of second-order internal resonance, needing extra calculations and at least two master coordinates, is developed in Section 3.3.

Reduction and normal form style parametrisation
The normal form procedure as exposed in Section 2 is a complete nonlinear change of coordinates where the number of input and output variables is the same. As underlined in [9], reduction can then be applied by selecting only a subset of few master normal variables, and cancelling all the others, such that the size of the z vector falls from 2N to 2m with m N the number of master coordinates. By doing so, one modifies the nonlinear one-to-one diffeomorphism to a parametrisation procedure of the underlying invariant manifold (or Lyapunov subcenter manifold (LSM) in a conservative context), the existence of which is guaranteed by the Lyapunov center theorem [32,35]. The derived expressions are then identical to the one found by first assuming a parametrisation, then solving the resulting tangent and normal homological problems using the normal form style, as stated in [27].
In the framework of mechanical systems, the parametrisation of invariant manifolds has been extensively treated in [28,29]; the application of the normal form approach to model order reduction is instead adopted in the developments led in [9,21,77].
A geometrical interpretation is reported in Fig. 1(a) where a two-dimensional representation of the fourdimensional subspace associated to two modes φ s and φ m is reported. With modal coordinates, the dynamics is expressed within an orthogonal basis formed by the eigenmodes. Due to the presence of invariant-breaking terms [10], linear eigensubspaces lose the invariance property. Invariant manifolds (LSM) emerging from modal spectral subspaces are highlighted in red and green respectively, their curvature expressing the non-resonant coupling between modes and thus the amplitude-dependence of modal quantities in a nonlinear framework. The normal form allows expressing the nonlinear dynamics in an invariant-based span of the phase space. Once that obtained, reduction to an invariant subspace is simply given by cancelling slave normal variables. Fig. 1(b) illustrates the application of the nonlinear mappings. Thanks to normal form transformation, non-resonant coupling terms are removed, hence the invariant manifold (LSM) tangent at the origin to a given mode lies on the plane identified by the (z m , z m+N ) pair with associated Ψ Therefore, the normal transform allows expressing the dynamics with normal coordinates that are nonlinearly related to the original ones through the nonlinear mappings Ψ (z), Υ (z). In a more general framework, reduction is achieved by selecting a set Φ m of master modes and identifying the set Z of indices such that Ψ (1) s ∈ Φ m . All other coordinates are set to zero: which yields the following mapping expansions in indicial notation: and reduced dynamics: Note in particular that under the assumption of no internal resonance then the quadratic term vanishes, f skl = 0, and explicit expressions for f (3) sklm are given in Eq. (34). Also, the reduction procedure implies that only quantities associated to the master modes are required during computation, so that the complete eigenfunctions matrix Φ is not required. In the remainder of the paper, the case of the second-order DNF is selected. With this choice only secondorder internal resonances have to be tracked. On the other hand, the dynamics up to order three contains all the monomials, be they resonant or non-resonant. Applications of the third-order mapping would lead to cancellation of the non-resonant third-order monomials, but it would create higher-order terms. Keeping the development at second-order has the advantage that all higher-order internal resonance, starting from thirdorder, do not need a special attention since they are all preserved in the reduced dynamics. The next Section details the case of a single master mode while Sec. 3.3 is focused on the treatment of second-order internal resonance.

Single Mode Reduction
In this Section, the simplest case of a reduction to a single master coordinate is detailed, which is performed under the assumption of no internal resonance. Let us denote as φ m the master mode of interest, the master normal coordinates are thus selected as (z m , z m+N ). The reduction procedure consists in setting to zero all other coordinates: ∀p = m, p = 1, ..., N : z p = 0, z p+N = 0.
Since only z α , z β are different from zero, then only the projection on φ m is required: which yields the following relation to compute the reduced dynamics coefficients: The associated third-order reduced dynamics reads: Since G is symmetric with respect to permutations of its arguments, the reduced dynamics can also be written as: This last equation represents the nonlinear dynamics, up to the third-order, of the system along the corresponding m th invariant manifold, with the key feature that all coefficients can be computed directly with operations from the FE dofs in physical space. It also underlines that the derived formulas are explicit enough so that a non-intrusive implementation of all calculations can be targeted. Indeed, evaluation of G and H terms does not require an explicit decomposition of the internal force in their linear, quadratic and cubic parts, but they can be evaluated for instance through non-intrusive methods, retrieving the direct algorithm proposed in [53] with a real formalism. This approach is particularly appealing for application of the present method with commercial finite element software since there is no special need to implement finite element routines to obtain the ROMs given by the DNF approach. More details about the potential non-intrusive implementation of the method are provided in Sec. 4.

Second-order Internal Resonance and Multiple Master Modes
In this Section, more detailed explanations on handling the case of a second-order internal resonance are given. This case has not been treated in earlier developments of the normal form approach for reduced-order modelling in [9,21,53]. A difficulty relies in the fact that resonant monomials will stay in the second-order normal form, since one cannot simply cancel all quadratic terms by stating f (2) skl = 0, ∀ s, k, l = 1, ..., 2N as in Sec. 2.3. In turn, second-order components of the mappings are not solvable anymore, see Eq. (28), since a second-order internal resonance makes the matrix to invert singular. The consequence of these two facts is that new cubic terms will arise in the third-order reduced dynamics and they must be properly computed.
As anticipated in Sec. 2.3, since the spectrum is composed by purely imaginary conjugate pairs, a secondorder internal resonance stems from any combination of three frequencies of the type: Or by using the entries of Λ: If this condition is met, then Eq. (28) cannot be solved since the operation [(λ k + λ l ) 2 M + K] filters ω 2 r from the eigenspectrum, and the resulting matrix is singular. Furthermore, the number of linearly dependent columns of the system is equal to the number of modes with eigenfrequency ω r . This is better highlighted by projecting to the modal basis via: showing clearly that any diagonal entry of Ω 2 equal to ω 2 r is cancelled. Hence, all (k, l) such that Eq. (45) is fulfilled must be treated with care.
The treatment of internal resonance with a direct approach operating in physical space is more difficult as compared to normal form computation operating in the modal basis. Indeed, the tracking of monomials associated to the resonance is immediate in the modal basis, and the known remedy consisting in cancelling the term in the mapping, and thus keeping the resonant monomial in the reduced dynamics, can be transparently applied [9,27,29]. With a direct approach in physical space, this is no longer possible since a single Ψ (2) kl vector embeds terms required to cancel all kl monomials. This implies that if one sets crudely Ψ (2) kl = 0, then all z coordinates would remain coupled by the resulting f (2) skl monomials, whereas the correct strategy consists in tracking only the terms that are involved in the resonance. A proper approach is to force the solution of the system to be orthogonal to the kernel of the matrix to invert in Eq. (28). This requires identifying all triplets (r, k, l) involved in the resonance relationship (45), and enforcing Ψ (2) kl and Υ (2) kl to be mass-orthogonal to the set of involved modes. Let us define as Φ R the set of eigenmodes such that a given pair of indices (k, l) satisfies a resonance condition, so that the enforced orthogonality condition reads: Since the mapping is imposed to be orthogonal to Φ R , resonant monomials that remain in the reduced dynamics are obtained by projecting Eqs. (27) and (26) on each φ r ∈ Φ R . In the context of model order reduction, only the subset of modes defined by Z (1/2) , with Z (1/2) subset of indices of Z lower or equal to N, is considered, hence coefficients are estimated as: (r+N)kl = 0.
Implementation details of the orthogonality condition are reported in Appendix E.
The presence of f (2) in Eq. (15) modifies the third order homological equation given by Eq. (30), hence one needs to rely on its full form. In the present work, second-order DNF is used, hence Ψ (3) = 0 and the third-order homological equation is used only to compute monomials at order n + 1, i.e. at third-order. Equation (30) in presence of second-order resonances modifies to: where new terms MΥ ks f (2) slm in Eq. (50a) and Ψ (2) sm f (2) skl + Ψ ks f (2) slm in Eq. (50b) appear. These terms highlight that low order monomials affect the estimation of higher order mapping and reduced dynamics coefficients. The latter are estimated by projecting Eq. (50) onto the modal basis: Equation (51) shows that in presence of second-order resonances f (s+N)klm . This has consequences when mapping the reduced dynamics to real-valued quantities. Further details are reported in Appendix F. Overall, the concepts introduced up to this point allow deriving the second-order DNF in the most general case. Second-order mappings Ψ (2) are obtained from Eq. (28) with potential application of the constraint given by Eq. (48) if resonance conditions are met. Non-trivial resonance conditions in a FE framework are identified as detailed in Appendix G. Third-order reduced dynamics coefficients are then computed from Eq. (51), which reduces to Eq. (34) if f (2) = 0. In order to be more specific, let us give more detail on the case of a 1:2 internal resonance. Let us assume that φ I and φ II satisfy the relation ω II = 2ω I , such that a special treatment is needed for these two modes as compared to all others. As a result, combinations (λ I + λ I ), (λ I+N + λ I+N ) filter ω II from the matrix required to compute Ψ (2) , and combinations (λ I + λ II+N ), (λ II + λ I+N ) filter ω I , hence f (2) = 0. Due to the 1:2 resonance, pairs (z I , z I+N ) and (z II , z II+N ) do not individually define invariant subspaces: a strong coupling exists between the two and a reduced-order model cannot involve only one of these two modes, they are intimately related and create a four-dimensional invariant manifold. A ROM with these two master coordinates is simply built by selecting ∀ p = {I, II}, p = 1, ..., N : z p = 0, z p+N = 0.
The development of the reduced model is then the application of the same procedure detailed for single master mode reduction in Sec. 3.2 but with two modes and with proper treatment of resonance conditions as detailed in the present section. By defining α = I, β = II, γ = I + N, and δ = II + N, the reduced dynamics is then expressed as: which highlights the presence of second order resonant monomials. This last equation can be mapped to real quantities to provide an oscillator-like equation. This is developed in Appendix F.

Real-Valued Mappings and Reduced Dynamics
In Secs. 2-3 the system is parametrised using a complex-valued mapping in z and the resulting reduced dynamics is complex-valued as well, as a result of the initial choices with a starting point at first-order (state-space formulation). A different point of view had been developed in [9,21,53], where the choice of real-valued mappings and real-valued reduced dynamics was enforced all along the calculations in order to fit the more standard oscillator equations for vibratory systems. This section is devoted to clarify the link between the two approaches by giving formulas allowing to pass from one representation to another. Mapping the reduced dynamics to real-valued quantities is efficient since it presents results that are easier to interpret from a physical point of view. Furthermore, solution of a real-valued system is often preferred.
Real-valued parametrisation is based on the introduction of two real normal coordinates, namely the displacement r and velocity s, which are tangent at the origin to modal displacements and velocities. This choice is not unique but provides an immediate link with quantities that are more familiar to engineers, i.e. modal quantities. Following the direct mapping proposed in [53], one can write: withΨ (r, s),Υ (r, s) polynomial mappings in r, s. Truncation at second-order yields: â kl r k r l +b kl s k s l +ĉ kl r k s l , withâ kl ,b kl ,ĉ kl ,α kl ,β kl ,γ kl unknown maps. Equation (55) is equivalent to Eq. (11) even if the derivation of direct expressions for mappings is somehow more involved due to the anti-diagonal structure of the linear part of the associated dynamics: A link between z and the couple (r, s) can be established in analogy with the linear transform that relates p to modal displacements u and velocities v, as reported in Appendix C. By letting y = {s, r}, one has: which yields: Focusing on displacement mappings, the link between the complex and real-valued expressions at second-order is retrieved by introducing the following split: Let us define vectors Ψ From Eq. (58) it is then possible to establish the relationship between complex-valued and real-valued mappings as:â which corresponds to the results provided in [53]. This procedure can be further extended for any expansion order.

Real-Valued Single Mode Reduction
Further insights into the relationship between complex and real-valued formulations can be derived by focusing on the case of a single mode reduction. In this case the polynomial mapping truncated at second-order simply reads: where, from Eqs. (61), the relationships between real and complex mappings simplify to: Recalling Eqs. (60), the evaluation of Ψ (N) mm and Ψ (P) mm is performed through: The explicit expressions for Ψ (N) mm and Ψ (P) mm allows for a non-intrusive implementation of the method in any FE software. Focusing on Eqs. (64), they both consists in the solution of a linear system with unknown displacement-like vector Ψ mm and imposed force vector −G (φ m , φ m ). In particular, to solve Eq. (64b), a linear static analysis has to be performed, whereas to solve Eq. (64a) a linear combination of stiffness and mass matrix has to be created in the software before solving the linear system. Regarding the computation of the force vector G (φ m , φ m ), it can also be non-intrusively calculated by means of the STEP method [37,38]. For example, by imposing a displacement along the mode φ m , firstly with a positive, then with a negative modal amplitude, it is possible to extract the vectors G (φ m , φ m ) and H(φ m , φ m , φ m ) respectively.
The reduced dynamics reported in Eq. (44) in the case of single master mode m is obtained by selecting the master coordinates as r m = z α + z β and s m = iω m (z α − z β ). Application of Eq. (57) to the reduced dynamics in Eq. (44) then yields: where Ψ (1) m = φ m . Finally, using Eqs. (63a)-(63b) to makeâ mm andb mm appear, the third-order single-mode reduced dynamics finally reads:r where the introduced coefficients are given by: The explicit form of these equations makes the computation of the reduced dynamics possible without the need of extrapolating the full tensors G and H. In fact, the nonlinear force vectors H(φ m , φ m , φ m ), G(â mm , φ m ), and G(b mm , φ m ) can be either computed in an intrusive manner by integration over each element in the FE code, or non-intrusively by means of the STEP method. This result is also in full accordance with previous derivations reported in [53]. As detailed in [53], the STEP method allows extrapolating the required nonlinear force vectors in a non-intrusive manner from any finite element code by imposing a series of displacements on the structure and subsequently extracting the resulting forces. If mappings and reduced dynamics coefficients are not truncated, then higher order terms appear in the reduced dynamics and the whole process can be pushed further. It is possible to show that also in the case of higher order calculations, the method can be implemented non-intrusively even though the number of calculations required in the STEP method would grow consistently with the order.

Applications
The DNF method is here validated on three different MEMS structures with complex geometries, with and without internal resonance. MEMS are known to be operated in near-vacuum packages and are thus subjected to very small damping values. They are also generally actuated at resonance to fulfil technological requirements such as high sensitivity or large drive displacements. This in turn makes MEMS systems subjected to large displacements and geometric nonlinearities, hence making the DNF approach ideal for developing a predictive ROM strategy for design purposes. MEMS actuation is performed using either electrostatic, magnetic, or piezoelectric actuation. These types of actuation introduce nonlinearities that are not taken into account in the present work, their treatment via the DNF approach being postponed to further studies.
Let us denote as φ I the actuated mode, hereafter referred to as drive mode. We assume that the excitation is provided thanks to a body force proportional to φ I with a driving frequency in the vicinity of ω I , the eigenfrequency of the drive mode. The losses are modelled with the assumption of mass-proportional damping, with coefficient equal to ω I /Q, where Q is the quality factor of the drive mode. The resulting model for MEMS structures in the physical space of the dofs of the FE mesh then reads: with κ a scalar load multiplier representing the intensity of the forcing. The damping model is pragmatic and corresponds to an established practice in microsystems. Indeed simplified models are often utilised to identify a unique damping parameter, i.e. a single Q value, for the whole structure. This assumption is of current use in MEMS community, and it has shown to be reliable for a number of test cases [78,79]. The ROMS are derived thanks to the second-order DNF method, where the nonlinear mapping is truncated at second-order while the reduced dynamics is truncated at third-order. The reduced model is then obtained by mapping Eq. (15) to real variables, following the guidelines given in Section 4. Modal forcing is then simply added at the right-hand side of the reduced dynamics. This assumption has already been used in [9,21,80,53], and has been justified in [9,21] by the fact that the variations of the time-dependent manifolds are of second-order as compared to other effects. A mathematical justification of this assumption has also been given in [30]. Since a mass-proportional damping is selected and dissipation values are small, the simplest solution consists in adding directly the modal damping factors to the selected master modes. This simplification can be overcome by adding a more complex formulation of the damping, that takes into account the losses of all the slave modes to better represent the damping on the invariant manifold, following the general formula presented in [21,53], and leading to a nonlinear form of the damping in the reduced dynamics. As shown in [53], this more complete model is particularly meaningful when using stiffness-proportional damping. In case of massproportional Rayleigh damping, the decay rates of the slave modes are rapidly negligible so that the assumption of modal master damping is sufficient.
The proposed model is applied for the analysis of three MEMS structures showing different types of nonlinearity: a MEMS micromirror that undergoes large rotations, a beam resonator featuring 1:3 internal resonance, and a shallow arch resonator showing 1:2 internal resonance.

MEMS Micromirror: Single Mode Reduction
The first device addressed is a MEMS micromirror developed by STMicroelectronics ® . The device in shown in Fig. 2(a-b). The structure is composed of a circular reflective surface with a diameter of 3000 µm and is connected to the substrate through two torsional springs. Actuation of the device is performed with two pairs of lead-zirconate titanate PZT patches deposited on top of four trapezoidal beams. Patches are highlighted in orange in Fig. 2(a). Silicon is anisotropic and a general discussion of its mechanical properties can be found e.g. in [81,67].
The device is operated at resonance in the vicinity of its first mode with eigenfrequency f 0 = 2266 Hz, corresponding to a rotation of the main mirror around the axis that passes through the two torsional springs. The displacement field of the mode is represented in Fig. 2(c). The first six eigenfrequencies of the structure are reported in Table 1, highlighting that no second-order resonance condition occurs between the first and higher modes. The ROM is obtained by imposing Ψ = φ m with m = 1 (master mode), and by computing the mappings in Eq. (39). Noticeably, only two linear systems with symmetric matrices must be solved to obtain the second-order mapping and third-order reduced dynamics by exploiting the symmetries of the Λ operator. Furthermore, only the eigenvalue and the eigenfunction of the master mode is required. The operator β ) needs to be computed only for the master mode and this step has moderate impact on the total running time of the analysis. Overall the technique is suitable for heavy parallelisation and vectorisation on modern processors. In the present contribution, all calculations have been realised thanks to a custom FE code developed by the authors without resorting to any STEP-like computation.
Reduced dynamics coefficients are obtained as detailed in Sec. 3. The results of the proposed model are compared with full-order simulations of the device, which are performed with the Harmonic Balance Finite Element Method (HBFEM) with pseudo-arc length continuation. The solution of the HBFEM model is obtained with Fourier series expansion up to 7 to ensure model convergence. The geometry is discretised with quadratic (15 nodes) wedge elements and the HBFEM model has 15,341 nodes, corresponding to 690,345 Fourier coefficients. Also the ROM is solved with the Harmonic Balance (HB) method with an order 9 Fourier expansion.
The comparison between HBFEM and single-mode ROM solutions, obtained for four κ values: 0.1, 0.2, 0.3, and 0.4 µm/µs 2 , is reported in Fig. 3(a). The chart presents the maximum rotation angle reached by the mirror for a single steady-state oscillation cycle and for a given frequency value. The results show an excellent agreement between the HBFEM and ROM solution at any rotation angle, up to the curve with the highest κ value. Furthermore, the invariant manifold defined by the expansion coincides with the trajectories of the HBFEM solutions, as highlighted by the phase-space representation of the system in Fig. 3(b), where the velocity-dependence of the manifold is clearly observable. A slight departure between the ROM and the full-order solution appears at very large rotations, as a a consequence of the second-order DNF method used. Using higher-order developments would likely improve the match even at higher amplitudes.
The reduction to a single master mode offers the main benefit of embedding within the equation of motion of the master mode all contributions of non-resonant slave modes, hence accounting for the curvature of the manifold. This in turn yields impressive computational performance compared to the HBFEM simulations. The time required to integrate the ROM was less than 1 minute using a custom HB solver. The time required to compute the four HBFEM solutions reported in Fig. 3 was equal to 1 week, hence highlighting how the present method offers outstanding results with little computational resources. Stability is not reported in the present analysis since the HBFEM solution has not this feature implemented. Considerations on the stability are postponed to Section 5.2.3. In order to better highlight the quality of the results obtained with the DNF method, we compare with another reduction strategy using implicit condensation (IC) as proposed in [39,41]. As shown for example in [50,42] with the case of a cantilever beam, IC method fails at reproducing correctly the nonlinear dynamics of structures when inertia nonlinearities are importantly excited. In particular, the IC method condenses statically the contribution of slave modes on the master mode trajectory, hence neglecting velocity dependent terms [49]. The comparison between the DNF and IC methods is reported in Fig. 4, highlighting the benefits of adopting the proposed method for modelling systems subjected to large rotations and, more in general, to large transformations. The IC reduction method overestimates the hardening behaviour, and cannot be used in such context of large rotations, whereas the DNF approach is giving uniformly valid solutions without the need of any extra assumption thanks to the invariance property. Further insight can be gained by looking more closely at the individual values of the coefficients of the reduced dynamics with single mode, Eq. (66). The MEMS micromirror under investigation does not show mode interactions even at large rotation amplitudes, yet methods based on static condensation show difficulties in modelling the nonlinear response of the device. For this device, the physical cubic coefficient h m in Eq. (66) is equal to 6.479·10 −10 µN/µm 3 . On the other hand, the first correction term computed by the normal transform A m equals -6.479·10 −10 µN/µm 3 , hence their sum cancels out on the leading cubic term r 3 m in Eq. (66). This is a remarkable result from a physical standpoint since it can be stated that elastic nonlinearities associated to the excited mode are not the main source of nonlinear response of the structure. On the other hand, B m is equal to 4.53·10 −12 µNµs 2 /µm 3 . This last term is a velocity dependent nonlinearity which accounts for the change in inertia of the system during motion. Therefore, the nonlinear response of the structure is conveyed by the change in configuration of the structure rather than elastic nonlinearities. Interestingly, similar observations on the particular values of the coefficients have been reported in the case of a cantilever beam.
An example of the performance of the proposed method is reported in Fig. 5, where the CPU times required to obtain the mappings and reduced dynamics coefficients of the MEMS micromirror for different mesh refinements are reported. The computation is performed on a desktop workstation with an AMD Ryzen™ 5 1600 Six-Core Processor 3.20 GHz and 64 GB RAM. The computational times highlight how the proposed technique is scalable to large models even with moderate computational resources. The integration of the ROM with the HB is not reported since it is has a negligible cost.

Degrees of Freedom [-]
Time [min]

Internally resonant MEMS resonators
In this Section, two different internal resonance scenarios are investigated on MEMS resonators having a beamlike and an arch-like structure, respectively. In the first case, a 1:3 internal resonance is investigated, where the second-order DNF method does not need extra assumptions. Indeed the third-order monomials are not affected by the second order mappings and they all remain in the reduced dynamics, making the excitation of the 1:3 resonance possible. This advantage of the second-order DNF as compared to higher orders has been already underlined in [53] and will be further commented herein. The second case is a 1:2 resonance which needs extra developments of the DNF approach as discussed in Section 3.3, with the reduced dynamics given in Appendix F. In both cases, reduction to the two internally resonant modes is able to catch the complex nonlinear phenomena and retrieve the frequency-response functions of the structures. The dynamics is reduced to a four-dimensional invariant manifold.

MEMS beam resonator: 1:3 internal resonance
In this section the analysis of a 1:3 internally resonant beam is reported. The device is shown in Fig. 6(a). The structure is made of two beams with a length of 1000 µm and a thickness of 5 µm. The two parts are connected in three points. The structure is not symmetric with respect to the center since the external connection points are located at different distances (see the values of P 1 and P 2 reported in the caption of Fig. 6). The material is polycrystalline silicon which is modelled as isotropic with a Young's modulus of 167 GPa and a Poisson ratio of 0.22. The first six eigenfrequencies of the structure are reported in Table 2, which highlights an almost perfect 1:3 ratio between mode 1 and mode 4, whose displacement field is reported in Fig. 6(b-c). The ROM is built by taking the pairs (z p , z p+N ) and (z q , z q+N ) such that Ψ (1) All remaining terms of z are set to zero. This implies a moderate increase of the computational burden as compared to the previous case of single-mode reduction, since only mappings and reduced dynamics coefficients that multiply a non-zero entry of z need to be computed.
As for the MEMS micromirror, the ROM is validated with full order HBFEM solutions with a Fourier expansion coefficient of order 7. The geometry is discretised with quadratic (15 nodes) wedge elements and the resulting HBFEM model is made of 12,906 dofs, corresponding to 193,590 nodal unknowns of the HBFEM problem. The quality factor of the model is set to 3000. Analyses are performed for κ values equal to 0.00025, 0.0005, 0.001, and 0.0015 µm/µs 2 . The ROM is solved again with the HB method with a Fourier expansion order equal to 9. Stability is not reported.
The comparison between HBFEM and ROM solutions is presented in Fig. 7(a). The accuracy of the reduced model is remarkable for any κ value. The strong nonlinear interaction is put in evidence by the loop appearing in the frequency response function, even at small κ values. A further comparison is shown in Fig. 7(b), where three trajectories are represented in the space (u 1 , v 1 , u 4 ). They have been obtained for the largest κ value selected, and have been chosen in the vicinity of the 1:3 resonance loop, as marked by the coloured points in Fig. 7(a). While the orange point is before the 1:3 resonance, the black point is exactly at resonance, and the green point after. Comparison between full and reduced solutions shows a very good match for the trajectories. Interestingly, the trajectories show a cubic shape with either positive or negative linear term, and they can be related to the two families of periodic orbits (backbone curves) arising in the 1:3 internal resonance. As a matter of fact, the solutions keep these shapes all along the FRF, before and after the 1:3 resonance. At exact resonance when the forcing frequency is equal to ω 4 /3, a transition form is obtained and reported in black in Fig. 7(b-c). The gain in computational time is again impressive, with a factor of 3000 between the ROM and the full-order HBFEM solution since a single FRF of the device requires approximately two days, while the solution of the reduced model requires less than a minute.

MEMS arch resonator with 1:2 internal resonance
The structure under investigation is the MEMS arch resonator depicted in Fig. 8(a). The arch has a length of 530 µm and it is made by two arched beams with a thickness of 5 µm connected at their midpoint. The structure is made in polycrystalline which is modelled as isotropic with a Young's modulus of 167 GPa and a Poisson's ratio of 0.22.
As reported in Table 3, the structure shows an almost perfect 1:2 internal resonance between mode 1 and mode 6. The displacement field associated to the two modes is reported in Fig. 8(b-c). Due to the 1:2 internal resonance, the reduced model is built by taking the (z p , z p+N ) and (z q , z q+N ) such that Ψ    and Ψ (1) q+N = φ 6 , hence yielding a two oscillators model as the one reported in Sec. 3, while the reduced dynamical equations are given in Appendix F. Damping and forcing is added following the general guidelines given at the beginning of Section 5.
The results provided by the ROM are compared with the full-order HBFEM solution of the MEMS device. HBFEM Fourier expansion order is taken up to order 9 to ensure convergence of the method. The geometry is discretised with quadratic (15 nodes) wedge elements and and the resulting HBFEM model is made of 5,913 dofs, hence the total number of nodal unknowns is equal to 112,347. The ROM is solved with the HB method with order 9. Curves are computed assuming a quality factor Q equal to 500 and analyses are performed for four κ values: 0.05, 0.1, 0.15, and 0.2 µm/µs 2 . The comparison between HBFEM and ROM solution is reported in Fig. 9a. The chart presents the maximum displacement reached by the device during an oscillatory cycle for each frequency value. The data highlight a perfect agreement between HBFEM and ROM solutions for each κ value. Only at the highest amplitude a small discrepancy is observed.
A set of trajectories of the solution is reported in phase space in Fig. 9(b-c). The trajectories are selected from the points in the FRF shown by bullets in Fig. 9(a), obtained for the largest forcing amplitude κ=0.2 µm/µs 2 . Comparison between full and reduced-order models on the trajectories shows again an excellent agreement. Interestingly, the shape of the trajectories in phase space can be related to the properties of the underlying backbones of the 1:2 system, as investigated in [82]. Indeed, two families of periodic orbits exist in that case, and are named as parabolic p + and p − modes depending on the sign of the curvature. The left-hand part of the FRF follows the p − backbone such that the shape of the trajectories reproduce the negative parabola in the correct axis, as shown in Fig. 9(b-c). The right-hand part corresponds to p + solutions. When the forcing frequency equals ω 6 /2, a transition form is obtained.

Stability Analysis
A complete characterisation of the nonlinear dynamic response associated to the structures reported in this Section requires evaluating the stability of the response and the detection of bifurcation points. However, the direct implementation of stability in large scale numerical procedures such as the HBFEM is challenging and asks for special attention. In particular, in our implementation of HBFEM, the stability is not computed yet. On the other hand, a number of available open-source continuation codes can be applied to small-scale systems like the ROMs previously discussed. Herein, stability analyses are performed using the ManLab continuation package [83] on the reduced dynamics, in order to give insight to the results presented in previous sections. The 1:3 internally resonant beam is addressed in Fig. 10(a) for a load multiplier κ value equal to 0.0015 µm/µs 2 . Starting from lower frequency values and moving towards larger values, the data highlight the presence of an unstable region enclosed between two saddle-node (SN) bifurcations. Afterwards, the response becomes stable again, until the system encounters another unstable region enclosed between two Neimark-Sacker (NS) bifurcations, hence suggesting the onset of quasi-periodic response of the system. The system becomes then stable again until a new unstable region enclosed between another set of SN bifurcations is found. Finally, the system becomes stable again and it retains the usual quasi-static response at higher frequencies. Overall, stability and bifurcation analysis of the system highlights important features that are compulsory for correct estimation of the structure response. Furthermore, this result was obtained with a negligible computational cost.
The stability analysis of the 1:2 internally resonant arch is reported in Fig. 10(b) for κ equal to 0.0015 µm/µs 2 . Starting from low frequencies, the chart shows the same sequence of stable and unstable regions observed for the 1:3 internally resonant beam. Indeed, starting from lower frequencies and moving upward, the system starts as stable. Then, an unstable region enclosed between two SN bifurcations is found. Afterwards, the system becomes stable again. An unstable region enclosed between two NS bifurcations is then observed. Stability is then recovered, only to be broken once again by another unstable region enclosed between two other pairs of SN bifurcations. Therefore, both 1:2 and 1:3 internal resonance examples show two pairs of SN bifurcations and one pair of NS bifurcations, all enclosing an unstable region of the system.  Overall, the proposed reduction technique does not only allow for an efficient estimation of the frequency response function of the device, but it also enables refined stability analysis that would be computationally expensive for large scale models.

Conclusions
In this paper, a direct normal form approach for model order reduction of the discretised equation of linear momentum for mechanical systems subjected to geometric nonlinearities is proposed. The direct computation has been rewritten with a state-space formulation as starting point, resulting in a symmetric and complex formulation of all the main quantities (mappings and reduced dynamics). The major outcome is the generalisation of the method by a formulation of the homological equations that keeps trace of the mechanical setting, thus facilitating further developments of the method. The complete rewriting also allows to formulating a comparison with the real-valued approach of the DNF originally proposed in [53], underlining the equivalence of the different methodologies. A special emphasis has also been put on the treatment of second-order internal resonance, a feature that had not been developed in all previous studies using this approach.
For illustration purposes, only second-order DNF has been used for numerical examples in this contribution, where the nonlinear mapping is truncated at second-order while third-order dynamics is computed. The method proves efficient for developing ROMs of large-scale finite element models, showing excellent computational performance and accuracy. Furthermore, the method does not rely on the slow-fast assumption between master and slave modes, since the velocity is directly accounted by the parametrisation procedure. This is a major achievement since it offers a uniformly valid and simulation-free method that could be blindly applied to any structure without the need of extra assumptions [48,49,50,51], for the same computational cost of other methods. Normal form theory makes the distinction between resonant and non-resonant couplings. Non-resonant couplings are automatically embedded in the master mode thanks to the curvature of the invariant manifold, so that there is no need of computing extra vectors (like MDs or dual modes) to achieve convergence of the ROM. Resonant couplings create strong energy exchange between the internally resonant modes and drastically modify the nonlinear dynamics of the system, such that modes in internal resonance must be appended to the ROM dimension.
A further advantage of the method is that it does not need to rely on the knowledge of all nonlinearities of the full order model, since it needs to evaluate system nonlinearities only for the mapping vectors associated to the master modes, hence further enhancing the scalability of the method. Furthermore, it has been underlined that the explicit results of the mappings presented in this contribution open the door to a non-intrusive coding of the method, hence making the DNF easily integrable within finite element commercial software.
The efficiency of the method is proved by studying nonlinear structures with strong reduction to 1 or 2 master modes. Excellent accuracy of the results have been demonstrated on three different examples of increasing dynamical complexity. Higher order developments may be required for studying more complex structures, for instance structures that show an initial softening behaviour followed by hardening response, a topic that is left for future works, together with the treatment of damping and forcing within the procedure.

A Notation
Compact expressions are introduced throughout the paper in order to improve readability. For polynomial representation of nonlinear terms, functional and indicial formulations are linked with the following relationships for Ts, T kl , T klm vectors of coefficients: T(x, x, x) = k,l,m The distributive property can be expressed as: Moreover, since Ts, T kl , T klm are constant vectors of coefficients used for polynomial representation, the following definitions hold for time derivatives:

B Explicit expressions for geometric nonlinear terms
We develop analytical expressions for quadratic and cubic nonlinearities in mechanical systems subjected to large transformations (geometric nonlinearities). The weak form of the linear momentum conservation equation mapped to reference configuration is expressed as: with Ω 0 the domain in reference configuration, ρ 0 the reference density, u the displacement field,ũ the test field defined over the space C(0) of functions that vanish on the portion of boundary where Dirichlet boundary conditions are prescribed. S is the second Piola-Kirchhof stress tensor, F is the deformation gradient, i.e. I + ∇u with ∇(·) denoting material gradient. The first integral in Eq. (72) is the kinetic power while the second term represents the power of internal forces. We introduce the Green-Lagrange strain tensor as e = (F T F − I)/2. Using a Saint-Venant Kirchhoff constitutive model, i.e. S = A : e with A fourth order elasticity tensor, Eq. (72) is expressed as: We notice that since e is nonlinear with respect to the displacement, Eq. (73) is nonlinear. We express e and its first variation in terms of the displacement gradient: δe(u,ũ) = 1 2 ∇ũ + ∇ Tũ + ∇ Tũ ∇u + ∇ T u ∇ũ .
The first terms is linear, the second and the third are quadratic and the last term is cubic. The power of internal forces term can then be rewritten as: δe : A : e dΩ 0 = k(u,ũ) + g(u, u,ũ) + h(u, u, u,ũ).
By taking a triplet of displacement fields u k , u l , and um the three terms k, g, and h are expressed as: h(u k , u l , um,ũ) = 1 3 Ω 0 e ns (ũ, u k ) : A : e ns (u l , um) + e ns (ũ, u l ) : A : e ns (um, u k ) + e ns (ũ, um) : where Eq. (77a) is linear with respect to the displacement, Eq. (77b) is quadratic, and Eq. (77c) is cubic. Upon finite element discretisation of Eqs. (77), one obtains: whereŨ denotes the vector of nodal values for the test function and subscript (·) h identifies approximated quantities. In the present work, the terms in Eq. (77) are computed directly as proposed in [84] and in [85] for MITC shell elements. Non-intrusive methods can be used as well, as for instance the STEP [37,86].

C Physical System diagonalisation
In this appendix we recall, for the sake of completeness, some well-known results regarding the diagonalisation of the state-space formulation used in Eq. (5). Let the matrix Φ collect column-wise the real-valued eigenvectors φ s of the linear part of Eq. (1). By setting U = Φu and pre-multiplying Eq. (1) by Φ T , one gets: where Ω is the diagonal matrix that stores the eigenfrequencies of the system, and g(u, u), h(u, u, u) are the modal quadratic and cubic nonlinearities: Introducing the modal velocity v =u, Eq. (79) is written as a system of first-order differential equations. By letting x = {v, u}: with: Let us now introduce the matrix R defined as: with I the N × N identity matrix, and set x = Rp, with p the 2N-dimensional vector of generalised coordinates. In particular one has: us = ps + p s+N ∀ s = 1, ..., N, The linear operator a is diagonalised by R. Indeed, pre-multiplying Eq. (81) by R −1 , the following equation is obtained: with: In this framework, Λ is a diagonal matrix which reads

D Higher Order Expansions
In Secs. 2 and 3 the normal form approach has been introduced with focus on low order mapping and reduced dynamics. This choice is motivated by the simplicity of low order formulations, which in turn highlights the efficiency of the presented method for the development of reduced-order models. However, the method can be implemented in an algorithmic way for arbitrary order expansions. Indeed the general structure of homological equations are: with G (n) defined as the sum of all G(Ψ (a) , Ψ (b) ) such that a + b = n, and H (n) equal to the sum of all H(Ψ (a) , Ψ (b) , Ψ (c) ) such that a + b + c = n. Let Ξ (n) = G (n) + H (n) . For the sake of completeness, we detail hereafter the equations needed up to order five.
First order: Second order: Third order: Fourth order: Fifth order: (λ k + λ l + λm + λn + λp) MΥ This set of equations highlights the general structure of the problems to be solved and gives the framework for computation of higher-order automated solutions; the implementation of the method for a generic order is out of the scope of this work and will be the subject of future contributions by the authors.

E Implementation Algorithm
We provide herein some details of the implementation procedure, as schematically reported in Algorithm 1.
The input parameters of the model are a discretisation of the system domain (MESH), together with boundary conditions (BCS), and material properties (PROP). Furthermore, the user should provide the list of master modes to be selected in the reduction, collected in the subset Z (1/2) following the notation used in the main text.
The algorithm starts by computing the mass M and stiffness K matrices (LinInt) as well as the eigenmodes Φm and the eigenfrequencies Ωm (GenEig) of the selected master modes only. Since our basic assumption is that of a conservative vibratory system, all maps Ψ (2) are real-valued and the reduced dynamics coefficients f (2) , f (3) are purely imaginary. This implies that the formulation allows working on real double precision arithmetic even for the coefficients of the reduced dynamics by storing only the imaginary parts. First-order reduced dynamics coefficients f (1) are obtained as ±iΩ. First-order displacement mappings are the linear eigenmodes. Second order reduced dynamics coefficients f (2) and mappings Ψ (2) are then computed by iterating over all kl permutations. For n master modes, then each index spans from 1 to 2n. For each permutation, first we check that (λ k + λ l ) does not yield a resonance condition (CheckResonance2). Then, the matrix D required to compute Ψ (2) is assembled (AssemblyMapMatrix): where the two blocks MΦ R are used to impose mass-orthogonality between the mapping and the eigenmodes of the resonant monomials. The right hand side of Eq. (28), i.e. G(Ψ l ), is then assembled in a vector L of dimensions compatible with D (IntegrateFNL2). By solving the resulting linear system: with Λ R and Λ R+N being the diagonal matrices with entries λr and λ r+N such that Ψ (1) Rkl is the vector of monomials f (2) rkl that cannot be cancelled in the reduced dynamics. From Eq. 95 both mappings and reduced dynamics coefficients are obtained. For second order DNF, a second iteration is then performed to compute f (3) . Therefore all klm permutations are spanned. For each permutation, L is assembled as (IntegrateFNL3): Finally, f sklm terms are computed by projecting L onto the master modes and by exploiting Eq. (51b) (step "Project" in the algorithm).
Algorithm 1: Implementation algorithm. The linear mapping relating complex to real-valued quantities reported in Eq. (57) can be applied to reduced-order models with an arbitrary number of master modes. The case of a 1:2 internally resonant system is here considered as an illustration, since this case is explicitly investigated in the example of Section 5.2.2. Let us assume that the two modes satisfying the 1:2 relation are mode 1 and mode 2 for the sake of simplicity. Therefore, model order reduction with the present method requires selecting as master coordinates (z 1 , z 1+N ) and (z 2 , z 2+N ). The reduced dynamics given in Eq. (15) when mapped to real-valued quantities using Eq. (57) yields the following equations forṙ: r 1 + g 112 (4b 111 r 1 r 2 s 1 + 2b 122 r 2 1 s 2 ) = s 1 , (97a) r 2 + g 211 (4b 212 r 1 r 2 s 2 + 2b 212 r 2 1 s 1 + 2b 222 r 2 1 s 2 ) = s 2 .
An important remark compared to previous developments is that the normal velocity s is not equal to the time derivative of the normal displacementṙ. This is a noticeable result that cannot be observed for real-valued reduced dynamics truncated at third order if no second order resonances between master modes are observed. By taking the derivatives of Eq. (97) with substitution truncated at third-order: r 1 + g 112 (4b 111 r 2 s 2 1 + (4b 111 + 4b 122 )r 1 s 1 s 2 − (4b 111 ω 2 1 + 2b 122 ω 2 2 )r 2 1 r 2 ) =ṡ 1 , (99a) r 2 + g 211 (−2b 212 ω 2 1 r 3 1 − 2b 222 ω 2 2 r 2 1 r 2 − 4b 212 ω 2 2 r 1 r 2 2 + 4b 222 r 1 s 1 s 2 + 4b 212 r 2 s 1 s 2 + 4b 212 r 1 s 2 2 + 4b 212 r 1 s 2 1 ) =ṡ 2 . This last set of equations represent the exact normal form of the system in real-valued formalism if truncated to second-order, with the presence of the only two resonant monomials. Third order monomials are all kept in the reduced dynamics in accordance with the assumption of truncating the nonlinear mapping at second-order and the reduced dynamics at third. Writing the exact normal form up to the third-order would need application of third-order mappings that would cancel all non-resonnt third-order monomials in Eqs. (100). Second-order physical coefficients are obtained as usual with: which also applies to the cubic physical modal coupling coefficients h pklm , reading From G and mappingsâ,b, one obtains the values for A pklm and B pklm as These equations follows the general computation guidelines for all the coefficients used in the case where no internal resonance is present, i.e. by direct application of the formula given in [53,9]. The only difference being that the coefficientsâ klm andb klm corresponding to the 1:2 resonant monomials, are now equal to zero, which does not impact the general formula. Finally, new terms that arise due to the presence of non-zero second order terms in a 1:2 internally resonant system read: withâ mkl ,b mkl obtained by the pre-multiplying the real-valued mappingsâ kl ,b kl by M and subsequently projecting on φ m . These new terms are of particular importance and are completely related to the existence of a 1:2 resonance. They are the consequence of the remaining second-order monomials in the normal form, that in turn creates new cubic terms that have been here computed and made explicit.

G On the identification of Non Trivial Resonance Conditions for Discrete Models
The identification of non trivial internal resonance conditions in finite element systems needs to be decided by giving a tolerance since perfect integer ratios are not possible due to round-off errors. As detailed in Sec. 2, the generic mapping Ψ (n) can be estimated through solution of a linear system of the type: with E a symmetric matrix which is singular in presence of resonance conditions. For trivial resonances, then det(E) is always zero. For non trivial resonances in floating point arithmetic the determinant of E is not zero, yet E is ill-conditioned. For each mapping term Ψ (n) , one has to invert a matrix of the form: with σ summation of Λ entries. Resonance condition implies that σ 2 is equal to minus any of the entries of Ω 2 . Therefore, resonance conditions of a given order can be detected by taking the norm of the difference between σ 2 and the square of any eigenfrequency of the master modes ωr. If the value is below a given tolerance ε: then a resonance condition is assumed.