Direct computation of nonlinear mapping via normal form for reduced-order models of finite element nonlinear structures
Alessandra Vizzaccaro, Yichang Shen, Loïc Salles, Jiří Blahoš, Cyril Touzé
aa r X i v : . [ c s . C E ] S e p Direct computation of nonlinear mapping via normal form for reduced-ordermodels of finite element nonlinear structures
Alessandra Vizzaccaro · Yichang Shen · Lo¨ıc Salles · Jiˇr´ıBlahoˇs · Cyril Touz´eAbstract
The direct computation of the third-order normal form for a geometrically nonlinear structure dis-cretised with the finite element (FE) method, is detailed. The procedure allows to define a nonlinear mappingin order to derive accurate reduced-order models (ROM) relying on invariant manifold theory. The proposedreduction strategy is direct and simulation free, in the sense that it allows to pass from physical coordinates(FE nodes) to normal coordinates, describing the dynamics in an invariant-based span of the phase space. Thenumber of master modes for the ROM is not a priori limited since a complete change of coordinate is proposed.The underlying theory ensures the quality of the predictions thanks to the invariance property of the reducedsubspace, together with their curvatures in phase space that accounts for the nonresonant nonlinear couplings.The method is applied to a beam discretised with 3D elements and shows its ability in recovering internalresonance at high energy. Then a fan blade model is investigated and the correct prediction given by the ROMsare assessed and discussed. A method is proposed to approximate an aggregate value for the damping, thattakes into account the damping coefficients of all the slave modes, and also using the Rayleigh damping modelas input. Frequency-response curves for the beam and the blades are then exhibited, showing the accuracy ofthe proposed method.
Keywords
Reduced Order modelling · Normal Form · Geometric nonlinearities · Nonlinear mapping
Model order reduction methods for geometrically nonlinear structures is an active field of research and numeroustechniques have been proposed in the past [1, 2, 3, 4, 5]. Among them, nonlinear mappings and reduction to
A. VizzaccaroImperial College London, Exhibition Road, SW7 2AZ London, UKE-mail: [email protected]. ShenIMSIA, ENSTA Paris, CNRS, EDF, CEA, Institut Polytechnique de Paris, 828 Boulevard des Mar´echaux, 91762 PalaiseauCedex, FranceL. SallesImperial College London, Exhibition Road, SW7 2AZ London, UKJ. BlahoˇsImperial College London, Exhibition Road, SW7 2AZ London, UKC. Touz´eIMSIA, ENSTA Paris, CNRS, EDF, CEA, Institut Polytechnique de Paris, 828 Boulevard des Mar´echaux, 91762 PalaiseauCedex, France Alessandra Vizzaccaro et al. invariant manifolds attracts a special attention since they are able by nature to provide better results than anylinear method [6, 7, 8]. Indeed, the invariance property is key to ensure that trajectories from the ROM alsoexist for the full model, and the curvature of the reduction spaces allows the use of less coordinates to describethe dynamics than any linear decomposition.Nonlinear normal modes (NNMs) offers a framework for such reduction methods. Since their first definitionby Rosenberg in the 1960s [9], they have witnessed numerous developments and definitions: family of periodicorbits [10, 11], invariant manifold tangent at the origin to the linear eigenspaces [6, 12], as well as numerouscomputational methods, from asymptotic developments [6, 13, 14, 15, 16] to numerical solution [17, 18, 19, 20, 21]including continuation of periodic orbits [22, 23, 24, 25]. Recently, Haller and coworkers revisited the mathemat-ical definition of NNMs in order to settle down a unified framework [5]. They defined a spectral submanifold(SSM) as the smoothest member of an invariant manifold family, tangent to a linear eigenstructure. In thecontext of conservative systems, this definition reduces to the Lyapunov subcenter manifolds that are filled withperiodic orbits, thus allowing unification of the different definitions given in the past, from Rosenberg to Shawand Pierre approach.Based on the introduction of invariant manifolds, NNMs have been used in order to derive reduced-ordermodels for nonlinear dynamical systems and geometrically nonlinear structures, see e.g. [18, 26]. A real approachto normal form has also been derived in order to make clear the connection with the invariant manifold approach,and generalise the derivation of reduced-order models thanks to a complete nonlinear change of coordinates [15,27, 28]. Successful applications to shells have been reported with this method [29, 30]. SSMs have been developedwith efficient parametrisation allowing to deal with large amplitude and produce efficient predictive models [8,31].Application of these methods to finite-element (FE) structures with geometric distributed nonlinearity,remains however scarce [26, 32, 33]. On the other hand, other methods have been proposed in the FE communityin order to tackle geometric nonlinearity, with a particular emphasis to non-intrusive or indirect methods [3], thatcan be derived from any commercial FE software by using only standard operations, and thus without the needto implement new calculations at the elementary level. The STEP (Stiffness Evaluation Procedure) allows for anon-intrusive computations of the modal nonlinear coupling coefficients [34, 35], but as such is not a reductionmethod since it allows only projecting the equations of motion onto the linear modes basis. Improvements of themethod have been proposed, see e.g. the dual modes proposed in [36], the combination of STEP with a POD(Proper Orthogonal Decomposition) [37], or a modified STEP where displacements are imposed on selecteddegrees of freedom only [38, 39]. Implicit condensation and expansion (ICE) have been proposed in order toeliminate the high-frequency components [40, 41, 42] thanks to a two-step procedure where a stress manifold isfitted after a series of static loading computations. Finally, modal derivatives, first introduced by Idelsohn andCardona [43] have also been used with this objective of finding out a complementary projection subspace withfaster convergence than linear modes [44, 45]. Recent extension proposes to use modal derivatives in a nonlinearmapping perspective to compute a quadratic manifold as reduction subspace [46, 47].Recent studies compare all these approaches mainly developed in the computational mechanics communityto invariant and spectral submanifolds in order to better understand their mathematical foundations, advantagesand drawbacks. Haller and Ponsioen derived general theorems showing that condensation and modal derivativesneed a slow/fast assumption between master and slave coordinates in order to ensure convergence to a correctsolution [48]. Elaborating on this idea, the quadratic manifold (QM) has been compared to the normal formapproach, showing again how the QM tends to invariant manifold only when slow/fast assumption holds [49].A numerical criterion was established, stating that this slow/fast separation could be claimed as soon as theeigenfrequencies of the slave coordinates are 4 times larger than those of the master. A particular emphasis wasalso placed on both cases of static (SMD) and full modal derivatives (MD), showing how quadratic nonlinearitiescan be treated incorrectly with the SMD approach, even if the slow/fast assumption is fulfilled. The ICE methodhas also been compared to invariant manifolds in [50], again underlining the requirement of slow/fast separation.Importantly, ICE and QM methods produce reduction subspaces that are not velocity-dependent and neitherinvariant, and these simplifications become an important drawback in the prediction given by the ROMs [49,50]. irect normal form for reduced-order models of finite element nonlinear structures 3
These results advocate better application of reduction methods based on invariant manifold for FE structuresas well, that would take properly into account the specificity of FE discretisation and also open the doors topossible non-intrusive procedures that could inherit the intrinsic properties of reduced subspace ( e.g. invariance,velocity dependence). Also, a main advantage of invariant-based techniques is that they are by nature simulationfree , i.e. without the need of resorting to a priori computations or outputs from the full-order model, as is thecase for example for the POD or the PGD (Proper Generalized Decomposition) methods [51, 52, 53, 54]. Suchan effort has been recently tackled for methods based on SSM [55]. By using a single master coordinate, theauthors gave general formula, up to the third-order, allowing for a direct computation of a ROM from thephysical coordinates. The proposed method is however limited to a single master mode in the proposed version.Even though the extension to multiple modes is theoretically possible, following e.g. [56], the procedure needsto be revised and developed further for direct computation.The objective of this contribution is to extend the results obtained thanks to the real normal form ap-proach [15, 27] to FE structures. A main limitation of the method as it was presented in [15, 27] is to rely on theequations of motion written in the modal basis. Since this projection step is out of reach for a number of FEproblems involving very large number of degrees of freedom, the method has been completely rewritten in orderto compute the needed quantities directly from the FE unknowns. In that manner, a full nonlinear mappingis introduced, allowing one to compute directly the reduced dynamics in an invariant-based span of the phasespace, up to the third-order.The paper is organised as follows. Section 2 details the method, from the theoretical foundations to practicalconsiderations on implementation. More particularly, Section 2.1 recalls the general framework on geometricnonlinearities used in this contribution, while Section 2.2 briefly recalls the nonlinear mapping defined from thenormal form approach, using modal coordinates as starting point. Then section 2.3 describes the calculationsneeded to perform the direct computation of the second-order normal form and gives its associated reduceddynamics. It also offers a short comparison to the quadratic mapping proposed in the framework of modalderivatives in order to understand the gain brought by the normal form approach. Section 2.4 extends themethod to the third-order. All these results being obtained in the case of no internal resonance between theeigenfrequencies, section 2.5 explains how the method needs to be adapted to account for possible internalresonance. Since most of the engineering application are for forced and damped systems, the addition of externalforce and damping is detailed in section 2.6. For the damping, a simplified expression is provided, that allowsto take into account the case of lightly damped systems, and is specified to the case of Rayleigh damping,which is the most common form of losses customary introduced in FE codes. Section 2.7 closes the theoreticalpart by giving more details on the implementation and comments on the non intrusiveness characteristics ofthe technique, thanks to an algorithmic presentation. Section 3 shows simulation results obtained on differentstructures. The conservative backbone curves (frequency-amplitude relationships) for a clamped-clamped beamand a fan blade are shown, illustrating the versatility of the method in taking into account internal resonancesbetween nonlinear frequencies by adding more master coordinates in the ROM. Finally, frequency-responsefunctions (FRFs) including damping and harmonic forcing are shown for both cases, illustrating how the methodcan handle these additional terms and produce efficient and reliable ROMs. X gathers all the degrees of freedom (dofs) ofthe model (displacements/rotations at each nodes) and is N -dimensional. The equation of motion reads: M ¨ X + KX + G ( X , X ) + H ( X , X , X ) = , (1)where M is the mass matrix, and K the tangent stiffness matrix. Geometric nonlinearity involves a nonlinearstrain/displacement relationship and a linear elastic behaviour law. In the framework of three-dimensional Alessandra Vizzaccaro et al. elasticity, it implies that only quadratic and cubic polynomial terms have to be taken into account [3, 57, 32, 33],they are expressed thanks to the terms G ( X , X ) and H ( X , X , X ), using a functional notation for the quadraticand cubic terms with coefficients gathered in third-order tensor G and fourth-order tensor H . The explicitindicial expression reads: G ( X , X ) = N X r =1 N X s =1 G rs X r X s , (2) H ( X , X , X ) = N X r =1 N X s =1 N X t =1 H rst X r X s X t , (3)where G rs the N-dimensional vector of coefficients G prs , for p = 1 , ..., N . In the first part of the paper, themain results are given for a conservative system, without taking into account a damping term in the equationof motion. However, section 2.6 will be devoted to damping and forcing, and an approximate method will bederived in order to show how light damping and forcing can also be taken into account in the reduced-ordermodeling strategy.For the following development, the eigenproblem needs to be defined. Let ( ω i , φ i ) be the couple eigenfrequency-eigenvector verifying the problem: ( K − ω i M ) φ i = . (4)Using normalisation with respect to mass, then one has the following equations fulfilled: V T MV = I , and V T KV = Ω , (5)with V the matrix of all eigenvectors, V = ( φ , ..., φ N ), I the identity matrix, and Ω a diagonal matrixcomposed of the square of the eigenpulsations, Ω = diag( ω i ). Since most of the presented results will useprevious calculations of normal form from the modal basis as a starting point, the equations of motion in themodal coordinates are used in the rest of the text in order to draw out parallels between the methods. Using thelinear change of coordinates X = Vx , with x the N -dimensional vector of modal displacements, the dynamicsreads: ¨ x + Ω x + g ( x , x ) + h ( x , x , x ) = , (6)where the third- and fourth-order tensors g and h expresses the nonlinear modal coupling coefficients. They arelinked to their equivalent G and H in the physical basis via: g ij = V T G ( φ i , φ j ) , (7a) h ijk = V T H ( φ i , φ j , φ k ) . (7b)In this contribution, the internal force vector is assumed to derive from a potential, and symmetry relationshipsexist between the nonlinear quadratic and cubic coefficients. For the quadratic coefficients, one has in particularthat G ( φ i , φ j ) = G ( φ j , φ i ), from which follows g ij = g ji , as well as g ijk = g jki = g kij , while for the cubic terms: h ijkl = h ijlk = h iklj = h ilkj .The main difficulties for deriving reduced-order models for geometric nonlinearity are linked, in a FE context,with the number of dofs N which can be prohibitively large. The number of coefficients involved in G and H tensors scales respectively as N and N , and due to the distributed nature of the nonlinearity, all the oscillatorsare nonlinearly coupled. On the other hand, the nonlinear dynamics exhibited by such structures are often simpleso that they can be captured by elementary systems of coupled nonlinear oscillators. However, in order to deriveefficient ROMs, the projection subspace giving rise to such low-order dynamics has generally a complex, curvedshape in the phase space. Obtaining such a projection in a straightforward and simulation-free context is theaim of a nonlinear mapping that could lead directly from the FE discretisation to variables describing thedynamics in a curved manifold. In the next sections, such nonlinear mappings will be given and detailed. Theyare based on the normal form approach previously developed from the modal coordinates, see [15, 27, 28], andbriefly recalled in section 2.2. irect normal form for reduced-order models of finite element nonlinear structures 5 e.g. [59] and [29] respectively for beams and shells, [7] for a comparisonwith the POD method, [30] for the correct prediction of the hardening/softening behaviour of shallow sphericalshells, and [28] for an overview.The nonlinear mapping is defined up to the third order and is based on an asymptotic expansion thatcan be pushed further if needed. The theory of normal form; Poincar´e and Poincar´e-Dulac’s theorems, lay thefoundation for its derivation. It is also defined for both displacements and velocities. Previous derivations showedthat the results are equivalent to the center manifold approach used by Shaw and Pierre [1]. Hence the nonlinearmapping allows one to express the dynamics in an invariant-based span of the phase space. This property iscrucial in order to defined accurate ROMs. Let us denote as y = ˙ x the velocity vector in physical space, R i and S j = ˙ R j the normal coordinates (respectively displacement and velocity). Let us also assume that the ROMis composed of n master modes. Even if it is generally assumed that n ≪ N , the method is given for a free n and can also be used with n = N , hence defining a complete change of coordinates between modal and normalcoordinates. The nonlinear mapping reads: x = n X i e i R i + n X i =1 n X j =1 a ij R i R j + n X i =1 n X j =1 b ij S i S j + n X i =1 n X j =1 n X k =1 r ijk R i R j R k + n X i =1 n X j =1 n X k =1 u ijk R i S j S k , (8a) y = n X i e i S i + n X i =1 n X j =1 γ ij R i S j + n X i =1 n X j =1 n X k =1 µ ijk S i S j S k + n X i =1 n X j =1 n X k =1 ν ijk S i R j R k . (8b)In these expressions, e i represents the unit eigenvector of the modal basis, with length N , composed of zerosexcept for the i -th entry which is equal to 1. The full expressions of all the reconstruction vectors a , b , γ , r , u , µ , and ν are respectively given in Appendix A and B, recalling the detailed expressions demonstrated in [15].The sole difference between these expression and those provided in [15], lies in the treatment of the symmetricterms in the summations. In these expressions, full summations are used with lower index covering from 1 to N . The other choice, often selected in such case, is to use upper-diagonalised forms for the tensors with orderedsummations ( s ≥ r and t ≥ s ), using the fact that the usual product is commutative. This choice is not retainedhere since it has been found easier to handle the expressions with full summations.The following important properties of the nonlinear mapping from Eqs. (8) are recalled: (i) the change ofcoordinate is identity-tangent in the sense that the first order is colinear to a given eigenmode. The correctingterms allows to take into account the curvature of the NNM (invariant manifold). (ii) the number of mastermodes n can be selected freely and in the most exhaustive case one can have n = N . The main advantage ofthese formula is that there is no need to recompute all the quantities when adding a new master mode. Thisis in contrast with methods based either on invariant manifolds, see e.g. [60], or the recently proposed methodbased on SSM [55].The reduced-order dynamics is written for the general case where no internal resonance exists between theeigenfrequencies of the system. When an internal resonance is present, some terms are vanishing in Eqs. (8),leading to extra terms staying in the normal form of the system. The theory is detailed in [15, 27, 28], and furthercomments will be provided in the next sections. In any case the dynamics onto the third-order 2 n -dimensionalinvariant manifold corresponding to the n selected master modes writes, in this general case without internalresonance, ∀ r = 1 ...n :¨ R r + ω r R r + ( A rrrr + h rrrr ) R r + ( B rrrr ) R r ˙ R r + R r n X j = r ( A rjjr + A rjrj + A rrjj + 3 h rrjj ) R j + R r n X j = r ( B rrjj ) ˙ R j + ˙ R r n X j = r ( B rjjr + B rjrj ) R j ˙ R j = 0 . (9) Alessandra Vizzaccaro et al.
As clearly emphasised in Eq. (9) where R r and ˙ R r terms have been factorised, the reduced dynamics doesnot contain invariant-breaking terms anymore. An important remark is also that since no internal resonancehave been assumed between the eigenfrequencies, no quadratic terms are present in Eq. (9), only cubic termscorresponding to trivial resonance stay in the normal form. Finally, one can note the presence of velocity-dependent terms in the reduced dynamics, reflecting the velocity-dependence of the invariant manifolds. Eventhough displacement and velocities are used as independent variables in the process, the method can simplyreduce to displacements only, one has just to replace S r = ˙ R r everywhere to obtain full expressions dependingon displacements only. In the same line, Eq. (8b) allows for a better reconstruction of the velocities but is notmandatorily needed and can be deduced from (8a).New fourth-order tensors A and B appear in these equations, their expressions from the modal basis can befound in [15] and are here recalled: A ijk = N X s =1 ( g is + g si ) a sjk = N X s =1 g is a sjk , (10a) B ijk = N X s =1 ( g is + g si ) b sjk = N X s =1 g is b sjk , (10b)where the last simplification stems from the symmetry of the quadratic tensor g .In a FE context, the main problem is that applicability of this reduction procedure needs as a first step thefull equations in modal basis, Eq. (6), with all the coefficients g and h known and computed. Even though theo-retically feasible by resorting to the Stiffness Evaluation Procedure (STEP) as described in [34], this computationis generally out of reach for much of the FE structures since the associated computational cost is prohibitivelylarge. Also, an a priori selection of some of the important coupled eigenmodes, which could indicate a way ofsolving that dimensionality issue, is clearly out of reach, as underlined for example for 3D elements in [39],where non-negligible couplings with very high frequency thickness modes have been exhibited. Hence there is aclear need for a direct computation of the normal form from the FE discretisation, since it will allow one to takedirectly into account all the coupled modes without any a priori or assumptions to formulate, and will define asimulation-free method applicable for efficient reduced-order modeling with geometric nonlinearity.2.3 Second-order direct normal formIn order to introduce progressively the cubic nonlinear mapping, we begin with a first step where only thesecond-order terms of the normal form are computed directly from the FE discretisation. This first step allowsfor a better understanding of the physical meaning of the involved coefficients, as well as some short comparisonswith the quadratic mapping introduced from the modal derivatives in [46, 47]. The reduced dynamics in thiscase will also be specified and the meanings of this assumption further detailed. In section 3 where numericalresults will be provided, this assumption will be again discussed and highlighted on given test cases.The quadratic mapping, deduced from Eq. (8), from the physical X coordinates of the FE model, reads: X = n X i =1 φ i R i + n X i =1 n X j =1 ¯ a ij R i R j + n X i =1 n X j =1 ¯ b ij S i S j , (11a) Y = n X i =1 φ i S i + n X i =1 n X j =1 ¯ γ ij R i S j . (11b)In these expressions, n stands for the number of master modes retained for building the ROM, and φ i is the corresponding eigenvector. The normal coordinates ( R i , S j ) have the same meaning as in Section 2.2,since they describe the dynamics in the same invariant-based span of the phase space. The expressions of thereconstruction vectors ¯ a , ¯ b and ¯ γ , are the equivalent to those obtained in modal basis, the only difference beingthat the new ones can now be computed directly from the FE nodes. As already remarked in section 2.2, the irect normal form for reduced-order models of finite element nonlinear structures 7 second equation (11b) is not mandatorily needed for the reduction process. Indeed, (11b) can be deduced from(11a), by adding the first-order assumption ˙ S p = − ω r R r .The detailed expressions of the second-order tensors ¯ a , ¯ b , ¯ γ , are now given. They have been derived fromthe coefficients obtained in [15] where the starting point was the modal coordinates. Appendix A explains howthis operation is done with a particular emphasis on the ¯ a tensor, also giving important computational detailswith regard to the direct computation in a FE context. Let us first define the two vectors ¯Zs ij and ¯Zd ij as: ¯Zs ij = ((+ ω i + ω j ) M − K ) − G ( φ i , φ j ) , (12a) ¯Zd ij = (( − ω i + ω j ) M − K ) − G ( φ i , φ j ) , (12b)which are needed to arrive at a compact expression for ¯ a , ¯ b and ¯ γ . These vectors encompass all the possiblesecond-order internal resonance thanks to the appearance of the sum ω i + ω j and difference − ω i + ω j of twoeigenfrequencies (thus giving the names of these ¯Z vectors with s (summation) and d (difference) as subscripts).Indeed, the matrix to invert will become singular in case of existence of a second-order internal resonance. Theseterms come from the denominators of the modal normal form and vanish in case of internal resonance, see [15,28] and section 2.5 for more details. They have been written with these expressions because they are moreconvenient on the computational viewpoint: ¯Zs ij and ¯Zd ij can be computed easily by solving a linear system,thus avoiding the matrix inversion. This will be further commented in Section 2.7 where specific algorithmicdetails will be highlighted.The expressions for the three vectors ¯ a , ¯ b and ¯ γ deduce from the two ¯Zs ij and ¯Zd ij only, underlining againthe fact that the second equation on Y is not independent from the first one, so that in a simplified version ofthe method it can be neglected. They read: ¯ a ij = 12 ( ¯Zd ij + ¯Zs ij ) , (13a) ¯ b ij = 12 ω i ω j ( ¯Zd ij − ¯Zs ij ) , (13b) ¯ γ ij = ω j − ω i ω j ¯Zd ij + ω j + ω i ω j ¯Zs ij . (13c)Vectors ¯ a ij , ¯ b ij and ¯ γ ij can be seen as correction vectors needed to take into account the quadratic curvatureof the invariant manifold. They thus share a common interpretation with the concepts of modal derivatives,as it has been introduced e.g. in [43, 44] . Indeed modal derivatives aimed at taking into account the nonlineardependence of an eigenvector with respect to perturbations in other modal directions. It has been recognisedthat this can be linked to the Hessian of the internal force vector and represents a manifold curvature in phasespace [46, 47, 49]. In order to make the connection with static modal derivative (SMD) more clear, let us assumethat a single master coordinate, say i , is retained for building the ROM. Restricting to a single master modemotion, then Eqs. (12) simplifies to: ¯Zs ij = ((2 ω i ) M − K ) − G ( φ i , φ i ) , (14a) ¯Zd ij = − K − G ( φ i , φ i ) . (14b)One can observe that Eq. (14b) is fully equivalent to the usual definition of static modal derivative θ ii onecan found e.g. in [44, 46, 49], to the multiplicative factor 2, i.e. ¯Zd ii = 2 θ ii . This means that in this simplifiedcase, θ ii exactly recovers the correct direction in phase space that express the quadratic couplings. However,this is not true anymore for cross-coupled SMD θ ij , i = j , since SMD is not able to make appear the differencebetween eigenfrequencies, nor the sum.A quite similar comment can also be addressed about ¯Zs ii , with the additional assumption of a slow/fastseparation between master and slave coordinates. Indeed, Eq. (14a) makes appear implicitly, in the terms of thematrix to inverse, the differences (2 ω i ) − ω s between the master coordinate i and all the slave modes s . Theslow/fast assumption requires that the slave modes are much more stiff than the master, implying ∀ s, ω s ≫ ω i .Hence if this assumption is well fulfilled, then ¯Zs ii will also tend to the direction pointed by ¯Zd ii = 2 θ ii . Thequadratic manifold produced by SMD will then tend to the one provided by direct normal form, since ¯ a will Alessandra Vizzaccaro et al. tend to twice the SMD and ¯ b to zero. This simple comparison underlines the common point between quadraticmanifold from SMD and the one proposed from direct normal form. In the general case, one understands thatthe formulas given by Eqs. (13) are more general, and only degenerate to the formulas proposed in [46, 47]when specific assumptions are met. As a conclusion, one can state that ¯ a and ¯ b vectors are in fact the correctcorrections that would have needed to be defined as the modal derivatives, since it gives the proper curvaturesof the invariant manifolds (NNMs) which are defined as the continuation of the underlying eigenmodes, tangentto the linear subspaces at origin. The interested reader is also referred to [49] where a complete comparisonbetween the two methods is provided.In the case considered in this section where a second-order normal form is studied, it is also interesting toderive the reduced-order dynamics arising from this choice, which will express the dynamics onto second-orderapproximations of the invariant manifolds. The dynamics is nonetheless expressed up to cubic order and reads, ∀ r ∈ [1 , n ]: ¨ R r + ω r R r + n X i =1 n X j =1 n X k =1 ( A rijk + h rijk ) R i R j R k + B rijk R i ˙ R j ˙ R k = 0 . (15)As compared to Eq. (9), this reduced dynamics does not contain quadratic terms also, as a consequence ofthe fact that no quadratic terms are trivially resonant. As long as no second-order internal resonance exists,these terms can be cancelled out. However all the cubic terms are present, without any simplification. This is themain contrast with Eq. (9), where all the non-trivial resonant monomial terms have been cancelled thanks tothe third-order term in the nonlinear mapping. Consequently Eq. (9) contains much less cubic terms than (15),and in particular all the invariant coupling terms can be cancelled. This point will be further commented insections 2.4 and 3. One can already note that reduced-order dynamics can be simulated with Eq. (15), meaningthat: (i) the invariant manifolds are approximated to the second-order only, (ii) all the resonant monomial terms,including non-trivial ones, are present. This might appear interesting in certain cases where strong higher orderinternal resonance are present, especially when they involve resonance between the nonlinear frequencies, thatone cannot easily foresee from a linear analysis. This point will be highlighted in section 3.1.In the particular case where a single master mode motion is selected, it is important to notice that Eqs. (9)and (15) reduce to the same, so that the reduced-order dynamics on a single invariant manifold is equivalent withsecond and third-order normal form. Assuming that the master coordinate has label r , then ∀ k = r, R k = S k = 0,and the dynamics of R r writes:¨ R r + ω r R r + ( A rrrr + h rrrr ) R r + B rrrr R r ˙ R r = 0 . (16)Indeed, in this simple case, the only cubic monomial terms staying in (16) are trivially resonant, and cannotbe cancelled by the cubic order terms from the nonlinear mapping. This means that if one is interested to thereduction to a single invariant manifold, then there is no need to compute further the cubic terms from thenormal form.The last point to address is the direct computation of the A and B , expressed in Eqs. (10), and appearingin the reduced dynamics. Since their definitions make appear the g tensor of quadratic coefficients from themodal basis, a direct computation is needed to circumvent the problem of computing all of these coefficients.For that purpose, let us define ¯A and ¯B the equivalent tensors in physical coordinates of A and B , which canbe simply obtained by premultiplying A and B by V − T . On the physical point of view, they are homogeneousto forces and not to displacements. Indeed, the linear change of coordinates for displacements reads X = Vx ,whereas for a force one has f = V T F , where f is the modal force and F the nodal force. Since, by definitionof the quadratic tensors in modal and physical coordinates, one has V − T g is = G ( φ i , φ s ), then: ¯A ijk = N X s =1 G ( φ i , φ s ) a sjk , (17a) ¯B ijk = N X s =1 G ( φ i , φ s ) b sjk . (17b) irect normal form for reduced-order models of finite element nonlinear structures 9 To arrive at a finalised expression involving only terms from the physical basis, one has to notice that P Ns =1 φ s a sjk = ¯ a jk (and similarly for ¯ b ), by again using the simple projection rule. The coefficients a sjk and b sjk can then bebrought into the brackets of G ( φ i , φ s ) leading to the final expressions: ¯A ijk = 2 G ( φ i , ¯ a jk ) , (18a) ¯B ijk = 2 G ( φ i , ¯ b jk ) . (18b)Eqs. (18) show how ¯A and ¯B can be computed directly by simple manipulations of the FE code (more algorithmicdetails will be given in Section 2.7). Once they have been obtained, their equivalent in modal basis simply read: A ijk = Φ T ¯A ijk , (19a) B ijk = Φ T ¯B ijk , (19b)where Φ is the reduced matrix of the master eigenvectors, which should not be confused with V : Φ containsonly the n master modes extracted from V , whereas V contains all the N eigenvectors.To conclude this section, we have derived general second-order formula that give rise to a quadratic mappingbased on the second-order normal form theory. The expressions have been compared to previous works usingstatic modal derivatives and the difference between the two methods have been underlined. In particular, it hasbeen shown how the proposed direct normal form (DNF) computation generalise the quadratic manifold SMDapproach and allows to take properly into account the internal resonance, the curvature of the phase space andthe velocity dependence. In the next section, we push the method further and give the detailed formulas neededto obtain the cubic terms of the nonlinear mapping.2.4 Third Order direct normal formThe computation of the third-order tensors is more difficult for two different reasons. First the expressions arelonger and more tedious. Second and most importantly, at the cubic order, one has to distinguish betweentrivial and non-trivial internal resonance. The reader is referred to [15, 27, 28] for detailed discussions, whilesections 2.4.2 and 2.5 recalls important definitions. Roughly speaking, if at second-order all the terms canbe cancelled as long as no order-two internal resonance exists, this is not the case anymore at the third-order,otherwise this would mean that the normal form (the reduced-order dynamics) is linear. As explained e.g. in [28],cubic monomial term X p on oscillator p is trivially resonant, and on the physical point of view this is meaningfulsince it will bend the frequency-response curves so as to produce hardening or softening behaviour [15].The nonlinear mapping up to order three reads: X = n X i =1 φ i R i + n X i =1 n X j =1 ¯ a ij R i R j + n X i =1 n X j =1 ¯ b ij S i S j , + n X i =1 n X j =1 n X k =1 ¯ r ijk R i R j R k + n X i =1 n X j =1 n X k =1 ¯ u ijk R i S j S k , (20a) Y = n X i =1 φ i S i + n X i =1 n X j =1 ¯ γ ij R i S j + n X i =1 n X j =1 n X k =1 ¯ µ ijk S i S j S k + n X i =1 n X j =1 n X k =1 ¯ ν ijk S i R j R k . (20b)Even though the second equation on Y is not mandatorily needed in the reduction technique, it is givenfor the following reasons. First, the original method has been developed including this second equation. Eventhough it can be deduced from the first in order to have expressions involving only displacements, the neededapproximation ˙ S r = − ω r R r is just the leading-order and might encounter limitations when the third order isincluded. Second, giving Eq. (20b) allows for direct reconstruction of the velocities, without deducing it fromthe first with an assumption, and for a light added computational effort. Finally, full dependency on the secondequation from the first could be lost in more complex cases, as e.g. including the damping (see [27] for theexpressions of the mapping where all the monomials are involved) or other forces in the starting equation. All these points need further analytical investigations that are beyond the scope of the present study. For all thesereasons, we keep the change of coordinates with both displacements and velocities.The expressions of the newly introduced coefficients ¯ r ijk , ¯ u ijk , ¯ µ ijk , ¯ ν ijk are given in the next subsection.The presentation will first consider the simplest case of the terms corresponding to non-trivial internal resonance.These terms can always be cancelled so that no specific treatment is needed as compared to the calculationspresented in the previous section. Then the trivially resonant monomials are considered. These terms shouldnot appear in the mapping so that the corresponding monomials will stay in the reduced dynamics (normalform). This specific treatment is detailed in section 2.4.2.The reduced dynamics after the cubic mapping is given by Eq. (9) in case of no internal resonance. It differsfrom (15) (reduced dynamics after second-order normal form) only by the fact that all the non resonant cubicterms have been cancelled. In some sense the order-three reduced dynamics appears simpler after the cubictreatment, and contains less possible nonlinear interactions. As a matter of fact, since the reduced dynamics istruncated at order three, the third-order treatment shown here is somehow incomplete. Indeed, non-resonantcubic monomial terms have been discarded, but the higher-order terms stemming from the nonlinear change ofcoordinates (quartic and quintic terms), have not been taken into account. We will show in section 3 that thismight have some important consequences when dealing with higher-order internal resonance. One must keepin mind that processing the cubic terms creates quartic terms that have not been considered in the presentstudy, thus reducing the quality of the third-order reduced dynamics. A full treatment would require to go tohigher orders, for example for a correct treatment of 5:1 resonance. Hence the third-order must be viewed asthe needed first step toward this achievement, which is however beyond the scope of the present study. In this subsection, the non-resonant terms are first considered. They correspond to the cases where none of thepossible combinations of the frequencies ± ω i ± ω j ± ω k is equal to another eigenfrequency ω s of the system. Thegeneral expressions for the fourth-order tensors ¯ r , ¯ u , ¯ µ , ¯ ν appearing in Eq. (20) have been obtained following asimilar analysis as the one shown in Appendix A for the quadratic terms. This analysis is detailed in Appendix B.They are here expressed from four independent vectors ¯Z0 ijk , ¯Z1 ijk , ¯Z2 ijk and ¯Z3 ijk given in Eqs. (22), in asimilar fashion as the quadratic tensors have been computed from ¯Zs ij and ¯Zd ij in the previous section.To simplify the presentation, it is convenient to define the following vectors obtained as a combination ofcubic forces: ¯P0 ijk = ¯A ijk + ¯A jki + ¯A kij + 3 H ( φ i , φ j , φ k ) − ω j ω k ¯B ijk − ω k ω i ¯B jki − ω i ω j ¯B kij , (21a) ¯P1 ijk = ¯A ijk + ¯A jki + ¯A kij + 3 H ( φ i , φ j , φ k ) − ω j ω k ¯B ijk + ω k ω i ¯B jki + ω i ω j ¯B kij , (21b) ¯P2 ijk = ¯A ijk + ¯A jki + ¯A kij + 3 H ( φ i , φ j , φ k ) + ω j ω k ¯B ijk − ω k ω i ¯B jki + ω i ω j ¯B kij , (21c) ¯P3 ijk = ¯A ijk + ¯A jki + ¯A kij + 3 H ( φ i , φ j , φ k ) + ω j ω k ¯B ijk + ω k ω i ¯B jki − ω i ω j ¯B kij . (21d)From them, one can obtain four independent cubic vectors of displacement that will be used to build the tensors: ¯Z0 ijk = ((+ ω i + ω j + ω k ) M − K ) − ¯P0 ijk , (22a) ¯Z1 ijk = (( − ω i + ω j + ω k ) M − K ) − ¯P1 ijk , (22b) ¯Z2 ijk = ((+ ω i − ω j + ω k ) M − K ) − ¯P2 ijk , (22c) ¯Z3 ijk = ((+ ω i + ω j − ω k ) M − K ) − ¯P3 ijk . (22d)Note that in the computational algorithm, each of these vectors is obtained by solving a linear system ofequations, and not by actually inverting the matrices, which is much more meaningful on the computationalpoint of view. The matrices of these systems are linear combinations of K and M and, similarly to the ¯Zs and ¯Zd vectors of the second order, cover all the possible combinations of ω i , ω j , and ω k . Four of them are neededin order to obtain all possible cases, corresponding to the splitting of the denominators appearing in the modalnormal form approach. The case where one of these combinations is equal to another eigenfrequency ω s of the irect normal form for reduced-order models of finite element nonlinear structures 11 system, thus making one of the matrices singular, will be specifically treated in Sec. 2.4.2. Here we assumedthat it is not the case, therefore the vectors ¯Z can be computed.Finally the expressions of the cubic tensors read: ¯ r ijk = 112 (cid:0) ¯Z0 ijk + ¯Z1 ijk + ¯Z2 ijk + ¯Z3 ijk (cid:1) , (23a) ¯ u ijk = 14 ω j ω k (cid:0) − ¯Z0 ijk − ¯Z1 ijk + ¯Z2 ijk + ¯Z3 ijk (cid:1) , (23b) ¯ µ ijk = 112 ω i ω j ω k (cid:0) − (+ ω i + ω j + ω k ) ¯Z0 ijk + ( − ω i + ω j + ω k ) ¯Z1 ijk + (+ ω i − ω j + ω k ) ¯Z2 ijk + (+ ω i + ω j − ω k ) ¯Z3 ijk (cid:1) , (23c) ¯ ν ijk = 14 ω i (cid:0) +(+ ω i + ω j + ω k ) ¯Z0 ijk − ( − ω i + ω j + ω k ) ¯Z1 ijk + (+ ω i − ω j + ω k ) ¯Z2 ijk + (+ ω i + ω j − ω k ) ¯Z3 ijk (cid:1) . (23d)Before generalizing these computations to the case of trivial internal resonance, one can note that the ¯Z vectors introduced in Eqs.(22) (where ¯Z is a shortcut notation for all ¯Z0 ijk , ¯Z1 ijk , ¯Z2 ijk and ¯Z3 ijk ) have somesymmetry relationships: ¯Z0 ijk = ¯Z0 ikj = ¯Z0 jki = ¯Z0 jik = ¯Z0 kij = ¯Z0 kji , (24a) ¯Z1 ijk = ¯Z1 ikj = ¯Z3 jki = ¯Z2 jik = ¯Z2 kij = ¯Z3 kji , (24b) ¯Z2 ijk = ¯Z3 ikj = ¯Z1 jki = ¯Z1 jik = ¯Z3 kij = ¯Z2 kji , (24c) ¯Z3 ijk = ¯Z2 ikj = ¯Z2 jki = ¯Z3 jik = ¯Z1 kij = ¯Z1 kji . (24d)Like the second order tensors, the third order tensors needed to build the velocity term Y are not inde-pendent from those to build X . More specifically, the expression for ¯ ν ijk , ¯ µ ijk , and all the ¯ ν and ¯ µ obtainedby permutations of the indexes ijk , can be expressed as a linear combination of the four linearly independentvectors ¯ r ijk , ¯ u ijk , ¯ u jki , ¯ u kij . In turn, these vectors can be unequivocally expressed as a linear combination of thefour linearly independent vectors ¯Z0 ijk , ¯Z1 ijk , ¯Z2 ijk , ¯Z3 ijk . This particularly means that only four independentvectors ¯Z0 ijk , ¯Z1 ijk , ¯Z2 ijk , ¯Z3 ijk are needed to compute 8 different vectors: ¯ r ijk , ¯ u ijk , ¯ u jki , ¯ u kij (equation on X ) and ¯ µ ijk , ¯ ν ijk , ¯ ν jki , ¯ ν kij (equation on Y ), leading to the conclusion that equation on Y is not independentfrom the one on X . Trivial resonances are the consequence of the purely imaginary eigenspectrum ± iω p of any conservative vibratorysystem. Then for two different indices ( i, j ) one is always able to create a third-order trivial internal resonancerelationship since iω j = + iω i − iω i + iω j . The mathematical consequence is that the treatment of the thirdorder terms is complicated by taking them into account. The physical consequence is that in vibration theory,backbone curves are bent to create either hardening or softening behaviour [15, 28].To better present the workaround that is needed to compute the third order tensors in case of a trivialresonance in a direct way, it is far more understandable to show how to compute them in modal basis first. Thenwe will extend the treatment to the direct computation. Even though trivial resonance can appear only whentwo indices are present instead of four (for a general cubic term from a fourth-order tensor), the presentationbelow is given for four different ( s, i, j, k ). This allows to maintain generality in the derivation, as well as moresimple coding since these calculations are typically run in nested loops.Let us start by writing Eqs. (22) in modal basis and by component for a non-resonant row s : Z sijk = ((+ ω i + ω j + ω k ) − ω s ) − P sijk , (25a) Z sijk = (( − ω i + ω j + ω k ) − ω s ) − P sijk , (25b) Z sijk = ((+ ω i − ω j + ω k ) − ω s ) − P sijk , (25c) Z sijk = ((+ ω i + ω j − ω k ) − ω s ) − P sijk . (25d) In modal basis, mass and stiffness matrices reduce to simple scalars thanks to their diagonalisation, and thebar over the variables are simply removed. Let us now suppose that the s -th row is resonant, meaning that ω s isequal to one of the combinations of ω i , ω j , and ω k . Then, one of the terms at the denominator in Eqs.(25) willbe zero and the inversion will no longer be possible. In such case, the terms Z sijk (where again the simplifiednotation Z sijk refers to Z sijk to Z sijk ) that aim at cancelling A sijk + h sijk and B sijk must be set equal to zero,thus avoiding the singularity (problem known as small denominator), and the above mentioned A sijk + h sijk and B sijk , that are no longer cancelled, must remain in the reduced dynamics.In order to write this set of operations in compact form, one has to first define the four denominator matrices: D = ( ± ω i ± ω j ± ω k ) I − Ω , (26)one for each combination of ± . We leave it generic because this operation will apply to each D . One of thesewill be singular because its s -th row will be zero. In order to be able to invert all other nonzero rows whileenforcing the constraint that the s -th row of Z must be zero, one has to fulfil the following system: ( DZ = P − e s P s , e T s Z = 0 . (27)In this way, the s -th row of the first system is the identity 0 = 0, the others satisfy Eqs (25), and the lastequality ensures Z s = 0. To express this system in compact form, the matrices must be augmented by the e s vector so that the system becomes: (cid:20) Z P s (cid:21) = (cid:20) D e s e T s (cid:21) − (cid:20) P (cid:21) (28)For a non-resonant row labelled p , the system is coherent with Eqs.(25); and for the resonant row labelled s ,the system enforces Z s = 0. One must notice that not only the matrix D has been augmented by a row and acolumn, but also the RHS vector is made of the force vector P plus a zero on its last row, as well as the solutionvector will be one row longer than the sought vector Z , with the last row being the reduced dynamics termthat cannot be cancelled. The last step is then to express the same operation in physical basis. The equivalentof D in physical basis is ¯D = ( ± ω i ± ω j ± ω k ) M − K and the equivalent of the vector e s in physical basis isthe vector M φ s . In fact, the constraint on ¯Z is expressed in physical basis by the mass orthogonality between ¯Z and the s -th mode φ s . The expression for the generic ¯Z in direct form finally reads: (cid:20) ¯Z P s (cid:21) = (cid:20) ¯D M φ s ( M φ s ) T (cid:21) − (cid:20) ¯P . (cid:21) (29)This augmentation procedure thus allows one to have a direct computation of all the needed trivially resonantthird-order terms. Computationally speaking, since the evaluation of the Z tensors is performed inside nestedloops spanning over three indexes i, j, k , the augmentation procedure has to be actuated every time two of thethree indexes are equal, say i = j , with the augmenting vector φ s being the non equal remaining index s = k .2.5 Internal resonanceWhen deriving the theory of normal form, as already stated in the previous sections and fully detailed in [15, 27,28], one has to take care of the occurrence of internal resonance. The general relationship for nonlinear resonanceis found from Poincar´e and Poincar´e-Dulac’s theorems. In the context of nonlinear vibratory systems composedof a purely imaginary eigenspectrum, it is customary to distinguish trivial internal resonance from non-trivialone. Third-order trivial resonances always exist. Consequently the normal form contains cubic-order terms thatcannot be cancelled, see Eq. (9). Each monomial term staying in the normal form is linked to a specific trivialinternal resonance, and the previous section details how to handle these terms in the processing of the cubicterms.Non-trivial internal resonances have not been taken into account yet since the main assumption holding fromthe beginning is that of a system free of these relationships. In case a third-order non-trivial internal resonance irect normal form for reduced-order models of finite element nonlinear structures 13 exist ( e.g. a relationship of the form ω p = ± ω k ± ω i ± ω j ), then exactly the same procedure as the one explainedin the previous section for trivial resonance has to be applied. This will let the corresponding cubic term of thenonlinear mapping to zero and make appear the corresponding monomial in the normal form. Consequentlyadapting the general reduced dynamics given by Eq. (9) to a case with e.g. e.g. ω p = 2 ω j or a combinationresonance ω p = ω j ± ω i ) then a special care has to be taken. Indeed Eq. (9) cannot be applied directly since thechange in the second-order tensors will create new and important changes in the cubic terms. At present themethod needs further refinement in order to compute fully the third-order tensor, so in this case it is advisedto use only the second-order DNF to create a ROM.2.6 Damping and forcingIn most of the cases one is interested in computing frequency-response curves of nonlinear structures usinga dedicated ROM, and not only to predict the backbone curve of the unforced and undamped system. Thissection is devoted to explain how external forcing and viscous damping can be taken into account in the proposedreduction strategy using direct normal form, based on previous results already shown in [15, 27, 28].For the external forcing, the strategy had already been proposed in [27, 28]. Since the reduced-order modelused normal coordinates linked to the invariant manifolds that are tangent to their linear counterpart at origin,one can simply add the modal force at the right-hand side. This approximation does not lead to an exactsolution since the correct formalism should take into account time-dependent manifolds. However, numericalresults reported in [61, 18, 62, 63] clearly shows that the time-dependent variation of the manifold are negligibleso that the results from unforced problem is meaningful. More recently, a mathematical proof that at theleading-order, this assumption of using a phenomenological forcing aligned with a curvilinear coordinates, iscorrect, has been given in [55], hence justifying again that this method can be used safely for including theforcing.For the damping, general formulas for computing the normal form with modal damping have already beenderived in [27]. However, the coefficients are more complex to handle and a number of new terms appear also inthe equations, see [27] for a more detailed comparison between conservative and damped normal form. Hencederiving the general expressions for a damped system overshoot the mark of the present study and is postponedto further developments. However, in order to take into account in the ROM a meaningful damping term thataggregates the losses of all the slave modes, a simplified formula is given here and will be used in the simulationresults presented in section 3.3. The general idea of the present development is to restrict to the case of lightlydamped systems, which is meaningful since geometric nonlinearities develop more easily in this context, whereasincreasing damping generally favours linear behaviours. Also, since the goal of the present work is to give aROM strategy applicable for FE structures, a special emphasis is put on the case of Rayleigh damping whichis currently used in FE codes, and general formula for a direct computation of new coefficients that take intoaccount Rayleigh damping is derived.The proposed strategy fully relies on the general results given in [27] that are simplified and modified inorder to tackle the problem at hand. More specifically, the proposed formula are given here only for the case ofthe second-order normal form. Indeed, computing the third-order terms in the nonlinear mapping with dampingincluded needs a further development that needs to be accomplished in a further work. The main assumptionis that of light damping. This means that, in an asymptotic development relative to damping factors, only theleading order term is considered, thus simplifying the general formulas given in [27].The dynamics of the structure in FE nodes now reads: M ¨ X + C ˙ X + KX + G ( X , X ) + H ( X , X , X ) = , (30)where the only additional term as compared to Eq. (1) is C ˙ X , and where the damping matrix C reads: C = ζ M M + ζ K K , (31) following the definition of Rayleigh damping, with ζ M and ζ K two free parameters. If expressed in the modalbasis, the equations of motion for a generic mode s reads:¨ x s + ζ s ˙ x s + ω s x s + N X i =1 N X j =1 g sij x i x j + N X i =1 N X j =1 N X k =1 h sijk x i x j xk = 0 . (32)The parameter ζ s can also be rewritten as ζ s = 2 ξ s ω s in order to make appear the ξ s modal damping ratio. If C = ζ M M + ζ K K , then ζ s = ζ M + ζ K ω s and ξ s = ( ζ M /ω s + ζ K ω s ) /
2. Note that in [27], analytical asymptoticexpansions for small damping are conducted by assuming the modal damping ratio to be small, ξ s ≪ ¯ α , ¯ β , ¯ c have to be added in the nonlinearmapping up to the second order, which now reads: X = n X i =1 φ i R i + n X i =1 n X j =1 ¯ a ij R i R j + n X i =1 n X j =1 ¯ b ij S i S j + n X i =1 n X j =1 ¯ c ij R i S j , (33a) Y = n X i φ i S i + n X i =1 n X j =1 ¯ α ij R i R j + n X i =1 n X j =1 ¯ β ij S i S j + n X i =1 n X j =1 ¯ γ ij R i S j . (33b)In the general case of arbitrary damping values, all the tensors are modified, meaning that ¯ a , ¯ b , ¯ γ needsalso to be recomputed in order to take properly into account the damping. However, as shown in [27], if lightdamping is considered, ξ s ≪
1, then all the expressions can be simplified, and one finds that the leading-orderterm for the ¯ a , ¯ b , ¯ γ is the same as the one in the conservative case (the first adjustment due to damping beingat order ξ s ). On the other hand, ¯ α , ¯ β , ¯ c involve only even powers of the damping in the asymptotics, so thatthe leading order is at order ξ s only: this term is thus now taken into account.The expressions up to first order for the tensor ¯ c with the assumption of Rayleigh damping reads: ¯ c ij = ( ζ M + 3 ω i ζ K ) ¯ b ij − (2 ζ K ) ¯ a ij + ( − ζ M + 2 ω i ζ K )( ¯Zss + ¯Zdd ) + ( − ζ M + 2 ω j ζ K )( ω i /ω j )( ¯Zss − ¯Zdd ) (34)with ¯ a ij and ¯ b ij given by Eq. (13) (conservative values), and the new terms as: ¯Zss ij = ((+ ω i + ω j ) M − K ) − M ¯Zs ij , (35a) ¯Zdd ij = (( − ω i + ω j ) M − K ) − M ¯Zd ij . (35b)As one can observe, the structure of the equations are close to the undamped case, the difference being that ¯Zss ij and ¯Zdd ij are more involved, but implying the same kind of operations as ¯Zs ij and ¯Zd ij . The full expression for ¯ c ij clearly underlines that it is a first-order term on the damping coefficients. Interestingly, only the dampingcoefficients of the master modes i, j appear in Eq. (34), meaning in fact that the assumption of small dampingneeds only to be retained for the master modes. This is particularly significant since ROMs are generally built forlow-frequency modes that have the smallest damping ratios. It underlines that the assumption of light dampingshould not be too restrictive since slave modes can have larger damping ratios.From Eq. (34) defining ¯ c ij , the computation of the two other second-order tensors simply reads: ¯ α ij = − ω i ¯ c ij , (36a) ¯ β ij = ¯ c ij − ( ζ i + ζ j ) ¯ b ij . (36b)To conclude on the nonlinear mapping, in case of damping, one simply has to keep the ¯ a , ¯ b , ¯ γ as in the conser-vative case, and compute the three new tensors ¯ α , ¯ β , ¯ c from (34)-(36).The reduced-order dynamics is also modified by the damping terms. Following [27], and restricting to thecase of the second-order normal form, the reduced dynamics now reads:¨ R r + (cid:16) ζ M + ζ K ω r (cid:17) ˙ R r + ω r R r + n X i =1 n X j =1 n X k =1 (cid:2) ( A rijk + h rijk ) R i R j R k + B rijk R i ˙ R j ˙ R k + C rijk R i R j ˙ R k (cid:3) = 0 . (37) irect normal form for reduced-order models of finite element nonlinear structures 15 Apart from the linear damping term, the proposed strategy to take damping into accounts leads to the ap-pearance of the term C rijk R i R j ˙ R k . The coefficient C rijk is a summation on all the modes including the slaveones. This added term can thus be understood as an aggregate nonlinear damping that takes into account thedamping of all the slave modes in order to better approximate the losses on the invariant manifold. From [27],the new coefficient reads: C rijk = N X s =1 g ris c sjk , (38)The last needed formula is the equivalent of Eq. (18) within the direct formulation in order to compute C rijk ,which is analogous to the fourth-order tensors A and B , in order to avoid the full computation of all the g rrs appearing in (38). After the same kind of manipulations as those reported in section 2.3, explicit expression for C rijk finally reads: C rijk = 2 φ T r G ( φ i , ¯ c jk ) , (39)hence allowing to express all the needed quantities so that they can be computed in a direct way from the FEdiscretisation.Note that all other terms appearing in the reduced dynamics are left unchanged as compared to (15). Inthe case of a single master coordinates labelled p , then only one extra term has to be taken into account in thereduced dynamics, C pppp R p ˙ R p . In case of more than one master coordinate, a simplified approach could havebeen to consider only the C pppp terms in the reduced dynamics ( e.g. in case of two master modes 1 and 2, only C , C could have been added). However the cross couplings terms have been found to be important in thereduction process, this will be illustrated in Section 3.3.2.7 Practical implementation and algorithmIn this section, a simplified overview of the computational algorithm is given, together with further considerationson the implementation and the non-intrusiveness of the method. In order to simplify the presentation, only thetreatment of the conservative terms is made explicit. Adding the damping terms can be done easily followingthe guidelines given in the previous section.In the present framework of distributed geometric nonlinearity where all the oscillator equations are fullycoupled with quadratic and cubic terms, direct extractions from the finite element model of the full tensors G and H is not feasible. The proposed method is direct in the sense that only the nonlinear force tensors G ( φ i , φ j )and H ( φ i , φ j , φ k ) are required with i, j, k ∈ [1 , n ] only spanning the master modes.These can be easily obtained thanks to the known STEP method, that computes the quadratic and cubicforces generated in the structure when a displacement along a combination of eigenvectors is imposed, see e.g. [34, 33]. The number of STEP method calculations is then equal to n ( n + 1)( n + 2) / n ) for the H tensor and n ( n + 1) / n ) for the G tensor if the symmetry property of the tensors are exploited.Also one of the main advantage of the proposed method is to rely on a direct computation solely involvinga selection of a priori master coordinates from usual physical considerations (on the nature of the forcing or thedynamics at hand). Since the ROM is directly expressed in the invariant-based span of the phase space, thecurvatures induced by the nonlinear mapping allows to directly compute the important, non-resonant couplings.There is no extra need to seek the slave modes coupled to the master ones by a convergence study: they aredirectly embedded in the nonlinear mapping and thus taken into account. The consequence is that there is noneed to compute the full eigenvector matrix V , only the master modes are needed and stored in the matrix Φ .A simplified overview of the algorithm is shown in Alg. 1. It relies on four main separate functions: Eig , StEP , Calc O2 and
Calc O3 . Eig is the usual eigensolver and is not explicated further.
StEP implements theStEP method following the general and known guidelines, see for example [34, 35, 33]. The starting point ofthe algorithm is given by the user, with the selection of the master modes of interest, from which the matrix Φ , together with their companion eigenfrequencies, are formed. The STEP method is then applied in order toderive the quadratic and cubic terms of the master modes of interest. Algorithm 1:
Direct Normal Form
Input : K , MFunctions :
Eig , StEP , Calc O2 , Calc O3
Output : Φ , ω , ¯ a , ¯ b , ¯ γ , A , B , hOptionalOutput: ¯ r , ¯ u , ¯ µ , ¯ ν Compute eigenproblem only for the master modes φ r , ω r ← Eig ( K , M ) Compute StEP method on the master modes Φ G ( φ i , φ j ) , H ( φ i , φ j , φ k ) ← StEP ( φ i , φ j , φ k ) Compute second order mapping tensors ¯ a , ¯ b , ¯ γ ¯ a ij , ¯ b ij , ¯ γ ij ← Calc O2 ( K , M , ω i , ω j , G ( φ i , φ j )) Compute StEP method on Φ , ¯ a and , Φ , ¯ b G ( φ i , ¯ a jk ) ← StEP ( φ i , ¯ a jk ) G ( φ i , ¯ a jk ) ← StEP ( φ i , ¯ b jk ) Obtain third order forces tensors ¯A and ¯B ¯A ijk ← G ( φ i , ¯ a jk ) ¯B ijk ← G ( φ i , ¯ b jk ) if DNF up to Second Order then Compute full third order reduced dynamics tensors h rijk ← φ T r H ( φ i , φ j , φ k ) A rijk ← φ T r ¯A ijk B rijk ← φ T r ¯B ijk else Compute third order mapping tensors ¯ r ijk , ¯ u ijk , ¯ µ ijk , ¯ ν ijk ← Calc O3 ( K , M , ω i , ω j , ω k , ¯A ijk , ¯B ijk , H ( φ i , φ j , φ k )) Fill third order reduced dynamics tensors solely with trivially resonant terms h rijk ← φ T r H ( φ i , φ j , φ k ) A rijk ← φ T r ¯A ijk B rijk ← φ T r ¯B ijk end The operations contained in the functions
Calc O2 and
Calc O3 have been already detailed in Sec. 2.3 andSec. 2.4, but their practical implementation is yet to be discussed. For the sake of simplicity, the discussion ishere restricted to the case of a single master mode, but of course the same considerations also apply to the moregeneral case of multiple master modes.The outputs of
Calc O2 are the ¯ a , ¯ b , ¯ γ second-order tensors. As shown in Section 2.3, they are derived fromthe computation of internal variables ¯Zs and ¯Zd . Let us detail this computation and comment its effectivenesswith regard to a non-intrusive method, in the case of a single master mode i . Eqs. (14) gives the expressions of ¯Zs and ¯Zd , which are written in their explicit form for a direct solution, i.e. with the unknown vectors aloneon the left hand side. In the actual implementation, there is however no need to perform a matrix inversionbecause, in the FE solver, a linear algebraic system is solved instead:((2 ω i ) M − K ) ¯Zs = G ( φ i , φ i ) , (40a) K ¯Zd = − G ( φ i , φ i ) . (40b)In Eq. (40b), the linear system can be solved directly by performing a linear static analysis that computes theunknown displacement vector ¯Zd when the structure is subjected to the force vector − G ( φ i , φ i ). This type ofanalysis is a standard operation in every FE software and its computational cost is quite small even for verylarge models. The same kind of operation has already been used in the context of SMD and is considered as a irect normal form for reduced-order models of finite element nonlinear structures 17 non-intrusive calculation. The linear system in Eq. (40a) is a bit different and involves a linear combination of K and M , but without an increase in the size of the system to solve. Despite this operation is not a standardone in a FE software, its computational cost is again quite low. To be able to perform this operation, the FEsoftware must allow the user to script and to do matrix operations online. In this way, the K and M matricesare never exported from the software but only the resulting vector is.The next steps of the algorithm is then to compute the ¯A and ¯B tensors following Eqs. (18). This is performedin lines 8 and 9 of the algorithm, by using the STEP function. Indeed, one can note that STEP is a calculationmethod, that can be applied with any input vectors needed. Even though, in its first derivation given in [34],it was only thought for a non-intrusive computation of the modal nonlinear coupling coefficients, so that theentries of the procedure were prescribed displacements along given eigenmodes, the procedure can also be usedmore generally, as done here with entries composed of one eigenvector and a vector ¯ a ij or ¯ b ij . Once ¯A and ¯B computed, their counterpart A and B are found by using a projection involving Φ , which is needed for thereduced-order dynamics.The last part of the algorithm distinguishes if the user only needs the second-order normal form or thethird-order. In case the third-order is needed, then the function Calc O3 is called, producing the computationsexplained in Section 2.4.For all our computations, the open software Code Aster [64] has been used. All the calculations have beensimply implemented using external scripts driving the master code, without the need of entering intrusively inthe code. e.g. [65, 11, 25].This salient feature shares common points with internal resonance between the eigenfrequencies of the system,already discussed in section 2.5, for which special care needs to be taken since corresponding to a resonantmonomial term in the normal form. In the present case, the resonance occurs between the nonlinear frequencies,at large amplitudes. Such a beam example has also been studied recently in [66], where these resonance loopshave been reported in Frequency-Energy Plots (FEP). Consequently the results shown in [66] will serve here asa guideline in order to detect if the same behaviour can be retrieved with the DNF.The model in [66] used simple 1D beam elements thus resorting to a simplified kinematics, and 30 elementsgave the space discretisation. In order to cope with a more realistic mesh, closer to general assumptions ofelasticity, 3D block elements with three displacements per nodes, have been used for meshing the clamped-clamped beam. The FE model consists of a mesh of 80 hexahedrons (20 along the axis, 2x2 in the section)with 20 nodes each, for a total number of 621 nodes, and 1863 degrees-of-freedom (dofs). The dimensions ofthe beam are L = 1 m along the z axis, b = h = 0 .
01 m. The material properties are the following: E = 210GPa, ρ = 8750 kg/m , ν = 0 .
3. The mesh is shown in Fig. 1 together with the eigenmodeshapes of modes 1 to4 and 6. In the rest of the paper, the displacement along x will be denoted as u , while the displacement alongthe in-plane direction z , is denoted as w . No motion is allowed along y .The six lowest eigenfrequencies of the beam, corresponding to bending modes, are reported in the first row ofTable 1. In the second row, the frequency ratio with respect to mode 1 is given, showing that the eigenfrequencyof mode 3 is a bit larger than 5 times the first. At the linear level, a 5:1 internal resonance does not exist.However, since the beam has a hardening behaviour, increasing the amplitude will make the fulfilment of 5:1ratio possible, so one could expect for the first NNM, a strong interaction in 5:1 ratio with mode 3, a result thathas indeed already been reported in [66]. The last row of Table 1 shows the frequency ratios with respect tomode 2. One can see that a 3:1 internal resonance may be excited once the hardening behaviour has sufficientlyincreased the nonlinear frequencies. So a 3:1 internal resonance loop with an interaction with mode 4 can be expected for NNM 2. As reported in [66], this will be indeed the case, as well as, at larger amplitudes, a 5:1resonance with mode 6. XY Z F
Fig. 1 – Mesh selection and geometry of the straight clamped-clamped beam, as well as eigenmode shapes correspondingto bending modes number 1 to 4 and mode 6. The bending vibration direction is along x while the in-plane direction is z .Mode Number 1 2 3 4 5 6Frequency (Hz) 50.900 140.74 277.09 460.64 692.93 975.85Frequency Ratio with mode 1 1 5.44Frequency Ratio with mode 2 1 3.27 6.93 Table 1 – Linear modes frequencies and frequency ratios with mode 1 and 2.
Fig. 2 shows some of the components of the quadratic and cubic tensors used in the nonlinear mapping,Eq. (20). The first line shows quadratic ¯ a ij vectors for some specific combinations: ¯ a , ¯ a , and ¯ a , and givealso a direct comparison with the corresponding static modal derivative (SMD) θ , θ and θ . As explained insection 2.3, in a simplified case where a slow/fast assumption can be assumed, then ¯ a ij vectors will point to thesame direction as static modal derivatives, and ¯ b ij vectors should be negligible. This property is indeed verifiedin this example: all three selected ¯ a ij examples are exactly the same as their SMD counterparts. Interestingly, ¯ a ij vectors are able to recover the most important quadratic couplings between eigenmodes. As analysed e.g. in [39], for flat structures such as the straight beam considered, quadratic couplings arise only in the oscillatorequations governing the in-plane modes dynamics, and involve two bending modes. With this respect, oneunderstands why the ¯ a ij vectors only involve in-plane motions, such that their main component is along w ( z direction) in Fig. 2. Also, slow/fast assumption in order to use safely SMD has been estimated in [49] as soonas the ratio between slave and master mode is larger than 4. Here the ratio between the first bending modeand the fourth axial mode is approximately 192 .
4. Second, ¯ a allows to directly recover the fact that for thefirst bending mode of a beam, the most important coupling is with the 4th axial mode, see e.g. [33, 39, 50] fornumerical examples highlighting this result. It is thus fully logical to observe that ¯ a have the shape of thisfourth axial. The direct consequence is that there is no need to analyse beforehand the mode couplings, sincethey are automatically embedded in the added quadratic tensors. Finally, as shown in [33], bending mode 3has a strong quadratic coupling with in-plane modes number 2, 6 and 8. Consequently, ¯ a shows an axialdeformation being a combination of these three mode shapes.The second line of Fig. 2 shows the same plot as the first line but now for the ¯ b ij vectors: ¯ b , ¯ b and ¯ b . These additional vectors are not taken into account if one uses the quadratic manifold from SMD. In thisspecific case of a straight beam where the slow/fast assumption is very well fulfilled, one can observe that the ¯ b ij vectors have, as anticipated, very small amplitudes, and can thus be safely neglected. These first results,highlighting how the nonlinear couplings are embedded into the additional vectors of the second-order termin the mapping, explains why in this simple case of a straight beam, the results provided by SMD are very irect normal form for reduced-order models of finite element nonlinear structures 19 z -0.500.5 w z -505 w z -505 w z -505 w x10 -5 z -505 w x10 -3 z -0.0200.02 w Fig. 2 – Illustration of the physical content of some of the quadratic and cubic tensors used in DNF approach. First line: ¯ a , ¯ a (black lines) compared to static modal derivatives θ , θ and θ (dashed green line). Second line: quadraticvectors ¯ b , ¯ b and ¯ b . Third line: cubic vectors: ¯ r , ¯ µ and ¯ ν . good and allows to retrieve correct predictions, as reported in [66]. However, adding curvature to the structuresand thus more complex couplings between modes and disappearance of the slow/fast assumption, it may beanticipated that SMD will not produce correct results anymore, as shown e.g. in [49]. The reflection of that factshould then appear in the behaviours of ¯ b ij vectors. This specific illustration is reported to a future work butpreliminary comparisons are already reported in [49].Finally, the third line of Fig. 2 shows three of the third-order vectors: ¯ r , ¯ µ and ¯ ν . Again, for aflat structure, cubic couplings involve only bending modes. Consequently, these correction vectors are in the x direction and involve mostly bending displacement u . Their shape also underlines the fact that, due to symmetryreasons, a clear separation exist between odd and even modes, the two families being coupled together with nocross-coupling between them. Consequently, ¯ r , ¯ µ and ¯ ν shows only combinations of odd modes andunderlines why mode 1 couples only with modes 3, 5, and so on, and the same for mode 2 with mode 4, 6 andso on.We now turn to the analysis of the backbone curves for the first and the second mode. In each case, areference, full-order simulation, is performed thanks to numerical continuation on all the degrees of freedom ofthe structure. A parallel implementation using harmonic balance method and pseudo arc-length continuation,as reported in [67], is used. In this beam case, due to some limitations of this algorithm, the beginning of thebackbone is numerically found by using a small amount of harmonic forcing. Then, as soon as the first pointsare found, forcing is removed and continuation restarted.Fig. 3(a) shows the reference result for the first mode obtained with the full-order model. As anticipated, aresonance loop appears along the backbone due to the excitation of the 5:1 internal resonance between mode1 and mode 3. Interestingly, two results from the computation of the full-order model are shown in Fig. 3(a).The first one has been obtained by using 9 harmonics in the numerical solution for continuation. As observed,this solution is not converged since the complete resonance loop is partially missed, in contrast to the convergedsolution obtained with 15 harmonics. This result underlines that in the resonance loops, complex high-order
320 330
340 350 360 370 380 390 4000123456 10 -3
320 330 340 350 360 370 380 390 4000123456 10 -3 Fig. 3 – Amplitude-frequency backbone curve for the first NNM of the clamped-clamped beam. (a) reference full-ordersolution obtained with 9 harmonics (black line) and with 15 harmonics (purple). Stability is not reported. (b) comparisonof the reference solution to different ROMs computed with the direct normal form (DNF). Light blue curve: direct normalform with a single master coordinate corresponding to mode 1. Dashed blue line: second-order direct normal form withtwo master coordinates 1 and 3, green line: third-order direct normal form with two master modes (1 and 3) and inclusionof the resonant monomials corresponding to 3:1 internal resonance. interactions are activated between the harmonics of the solution and a too crude truncation oversimplify theanalysis and misses important detail of the nonlinear dynamics.In Fig. 3(b), the reference solution is compared to different ROMs computed with the direct normal form(DNF). The first ROM is obtained by retaining a single master coordinate in the reduction strategy, corre-sponding to mode 1 (DNF single master mode, light blue curve). As it could be expected, the single-modesolution offers a drastic reduction, is able to recover the correct hardening behaviour of the beam, but is notable to retrieve the internal resonance, since a strong interaction with mode 3 is activated. The minimal reducedmodel should consists of at least two master coordinates. Consequently different ROMs including two mastercoordinates corresponding to mode 1 and 3 are tested.The second ROM is computed thanks to the second-order normal form presented in section 2.3, and twomaster coordinates (modes 1 and 3). The main advantage of this strategy is that all possible combinations ofhigher-order resonance (from order three) have not been dealt with the nonlinear mapping. Consequently theyare still present in the cubic terms and can express themselves if needed. This strategy is particularly meaningfulin the case considered since a 5:1 internal resonance is at hand, which would formally need to push all the normalform development up to order five, which is beyond the scope of the present study. The ROM obtained withtwo master modes and second order normal form is able to properly retrieve the 5:1 internal resonance, withoutany extra effort from the analysis. This is explained by the fact that the second-order resonances have beencorrectly treated, but then stopping at this order let all other possibilities free.On the other hand, using the third-order normal form for getting back this 5:1 relationship opens the doorsto new questions. In the best setting, one should go to order five and add only the resonant monomial termscorresponding to the 5:1 relationship. Since this calculation is far more difficult, different options are considered.The first, already commented option, is to stop at second order, leaving all higher-order terms there. This optiongives very good result in this case. Its main drawbacks are twofold. First an order of accuracy in the nonlinearmapping has been lost as compared to what can be easily computed, meaning that the approximation of theunderlying invariant manifold is less accurate. Second, the number of cubic terms in the ROM is more important,which is not a big issue since only two modes are considered. A second option would be to select two mastercoordinates and use the normal form up to the third order, using the solution obtained without consideringinternal resonance, see Eq. (9) for the reduced dynamics. In this case the same solution as with a single mastercoordinate will be retrieved, because no invariant-breaking terms are present in Eq. (9). irect normal form for reduced-order models of finite element nonlinear structures 21
A third option consists in adding invariant-breaking terms in order to excite the coupling and retrieve the 5:1internal resonance loop with only cubic terms. Since internal resonance appears due to nonlinear interactionsbetween harmonics of the solution, it is meaningful to add invariant-breaking terms related to lower-orderinternal resonance, in this case 3:1 and 1:1. In the first case (3:1 resonance), two terms are added to thedynamics due to 3:1 resonance: a R R term for the first oscillator (and all their ancillary terms involvingvelocities), and a R on the second master coordinate (corresponding to mode 3, again with all its companionterms with velocities). The dynamics of these ROMs are fully explicated in Appendix C so that the reader canget a better understanding of the equations used.The resulting backbone is shown in Fig. 3(b), and is named DNF third order + 3:1. Interestingly, this ROMcan recover the internal resonance loop of the 5:1, with only a slight departure from the reference solution.Finally if one also considers the terms due to 1:1, it means adding the R term on R equation, and R R monomial on R equation (see Appendix C for details). In this case the model is simply the same as the oneobtained with second-order normal form, since all the cubic terms are now present. So this case is equivalent tothe result already obtained with the second-order normal form, and gives a perfect match.
900 1000 1100 1200 1300 1400 15000
900 1000 1100 1200 1300 1400 15000 (cid:28)(cid:29)(cid:30)(cid:31) !"
Fig. 4 – Amplitude-frequency backbone curve for the second NNM of the clamped-clamped beam. (a) Reference solutionobtained with 15 harmonics (purple) compared to ROM considering only master coordinate of mode 2 (DNF singlemaster mode, light blue), and a ROM with the third-order normal form with modes 2 and 4 and the resonant monomialscorresponding to 3:1 internal resonance (DNF third-order + 3:1, dashed red curve). (b) Reference solution is comparedto a ROM obtained with the second-order normal form with 3 master coordinates corresponding to modes 2, 4 and 6(dashed blue), and to a ROM obtained with the third-order normal form, 3 master coordinates, and resonant monomialscorresponding to 3:1 between mode 2 and 4 and 2 and 6.
Fig. 4 shows the results obtained for the second NNM. As expected from the linear analysis and from theresults already reported in [66], the full-order solution shows two successive loops, corresponding to the fulfilmentof 3:1 internal resonance between the nonlinear frequencies of mode 2 and mode 4, and the 5:1 resonance betweenmode 2 and mode 6. Again, different reduced-order models are compared in order to better understand theircapabilities. A first ROM is built using a single master coordinate corresponding to mode 2, and compared tofull-order solution in Fig. 4(a). This ROM allows to predict the correct hardening behaviour, but cannot catchthe resonance loops since the resonant dynamics is intrinsically of higher dimension. A second ROM with twomaster coordinates corresponding to mode 2 and 4 is then built, using the third-order normal form. Since a 3:1internal resonance is expected, the resonant monomials corresponding to 3:1 resonance between modes 2 and 4are taken into account in this reduced dynamics. The result is shown in Fig. 4(a). It underlines that this ROMis capable of reproducing the loop of 3:1 internal resonance with a fair accuracy, but of course is not able tocatch the second loop since a nonlinear interaction with mode 6 comes into play.
Fig. 4(b) compares the results provided by two ROMs composed of three master coordinates correspondingto modes 2, 4 and 6. The first one is computed with the second-order normal form. As in the case of thefirst NNM, it gives a perfect match to the reference solution and is able to recover the two loops of internalresonances corresponding to 3:1 resonance with mode 4 and 5:1 resonance with mode 6. Note that in some sensethis ROM is minimal since only three master coordinates have been used. To our point of view, it is not possibleto reproduce such backbones with less than three master coordinates since nonlinear interactions are betweenmode 2 and 4 and then between mode 2 and 6. On the other hand, the same behaviour has been successfullyreported in [66], but their model were composed of 9 modes (3 linear modes plus 6 static modal derivatives).The last ROM tested still contains three master coordinates and applies the third-order normal form.Resonant monomial corresponding to 3:1 resonance between mode 2 and 4 have been added to reproduce thefirst resonance loop. For the second loop implying a 5:1 resonance, the same problem appears as with the case ofNNM 1: a perfect recovering of this resonance with the normal form theory would need to push the asymptoticdevelopment up to order five. Consequently, in the line of the previous study on NNM 1, only the resonantmonomials corresponding to a 3:1 resonance between mode 2 and 6 have been added. This ROM contains farless nonlinear terms in the reduced dynamics than the second-order normal form. It is able to well reproducethe two resonance loops, however with a small difference as compared to the full-order solution.These examples clearly underline that selecting the second-order normal form is generally a very goodoption when odd resonances are present. Indeed, comparing in this case the reduced dynamics given by secondand third-order DNF shows that the only difference is that the non-resonant terms have been removed in thedynamics for the third-order. On the other hand, the higher-order corrective terms that should come into playsince the computation of normal form is intrinsically an asymptotics, have not been taken into account, basedon the fact that the dynamics is truncated up to order three. This means that in the present version, thesecond-order normal form contains more information and is capable of reproducing higher-order resonance asshown here with the 5:1 case. The main drawback is that in the reconstruction process, the nonlinear mappingis less accurate since truncated to order two. But the gain in the reduced dynamics, which contains much moreterms, seems to be more important.3.2 FE model of a fan bladeIn this section, a FE model of a fan blade is selected in order to test the method with an application ofhigh interest in engineering, with a more complex geometry and large-amplitude vibrations exciting geometricnonlinearities. The backbone curve of the first mode of the structure will be computed with the DNF approachand compared to a reference solution obtained from direct continuation of all the dofs of the structure.The FE model of the fan blade is shown in Fig. 5a, and consists of a mesh of 2041 tetrahedrons with 10nodes each, for a total number of 3895 nodes (11685 dofs). The blade is clamped at its root as shown in Fig. 5a.The dimension of the bounding box containing the blade and the portion of disk at its root is: ∆X = 0 .
41 m, ∆Y = 1 .
14 m, ∆Z = 0 .
35 m. The span of the blade is approximately 1 .
02 m, the chord measured from trailingto leading edge is 0 .
48 m and the thickness in the centre of the profile is 0 .
01 m. The material properties arethe following: E = 104 GPa, ρ = 4400 kg/m , ν = 0 .
3. The frequency of the first five modes are reported inTable 2. The mode shapes are reported in Fig. 5b.
Mode Number 1 2 3 4 5Frequency (Hz) 14.592 40.571 49.055 170.36 197.03Frequency Ratio 1 2.78 3.36 11.67 13.50
Table 2 – First five eigenfrequencies of the fan blade and frequency ratios with the first mode.
This example is challenging for a number of reasons. First the nonlinear behaviour of the blade is close toa cantilever beam, for which inertia nonlinearity plays an important role [68]. In this case, reduction methods irect normal form for reduced-order models of finite element nonlinear structures 23 (a)
Blade mesh isometric and top view.
XYZ (b)
Blade modes
Fig. 5 – Blade mesh and first five modes. In Fig. 5a, the blade mesh is shown from two different views. Clamped faceshighlighted in yellow. In Fig. 5b, the first five modes of the blade are reported. based on static condensation are not reliable anymore since neglecting the most important phenomenon. Second,as compared to the beam, the shape of the fan is twisted and curved, making the problem closer to a shell thanto a flat structure. Consequently quadratic nonlinearities are much more pronounced and slow/fast separation isless fulfilled. For all these reasons, computing the correct type of nonlinearity (hardening/softening behaviour)of this kind of structure, still remains difficult for a ROM. As a matter of fact, no clear numerical demonstrationthat a reduced-order model is able to accurately compute the non-linear vibration of a Fan Blade has beenproposed so far. The existing methods published in literature for vibration of large 3D models cannot achievea better accuracy that 10% of error, in the best case, as compared to the full model [69].In Fig. 6, the second order tensors ¯ a , ¯ b relative to the first mode of the blade are reported together with thereduced dynamics parameters of the ROM for a single mode motion. The shape of the ¯ a vector shows that ashortening motion (proportional to R ) takes place as the vibration amplitude of the first mode (proportionalto R ) increases. This is a fairly expectable behaviour occurring in every cantilever-like structure. The correctionprovided by ¯ a thus serves the purpose of relaxing the stiffening effect that one would observe when this shorteningmotion is erroneously locked, see e.g. [47, 70]. The relaxation provided by ¯ a is what generates the A termin the reduced dynamics of the first mode. Importantly, the absolute values of A and h are nearly thesame and their sum is almost vanishing, meaning that the R term is small in the reduced dynamics. Unlikeoverconstrained structures, cantilever-like structures are isostatic therefore the slave modes quadratically coupledto the master are free to move, generating this particular behaviour that we have also observed in other cases.Since the correction brought by A makes the term A + h very small, one understands that muchlarger amplitudes needs to be reached in order to observe a remarkable curvature in the backbone curve.Secondly, the B term in the reduced dynamics becomes more important and convey a significant informationin order to retrieve the correct hardening/softening behaviour. These remarks underline why reduction methodssuch as static and implicit condensation cannot catch properly the correct behaviour in this particular case ofcantilever structure [38]. Indeed, since these methods neglect velocities and inertia effects in their development,referring to a static manifold instead of the correct invariant manifold, they can retrieve the relaxation on cubicterm but will neglect the inertial effect appearing in B , linked to a R ˙ R monomial term. As reported in [47],static modal derivatives also encounters convergence issues with cantilever structures when taking into accounttwo master modes in the ROM.The computation of the full-order solution also presents some challenges. As known from nonlinear continu-ation analysis on such complex structures, increasing the number of harmonics can make the solution more andmore difficult to achieve since more and more complex internal resonance with numerous loops appears in thebackbone. This has also been the case with this fan blade. A full model computed using a harmonic balancecontinuation method with 5 harmonics has been found unable to achieve moderate amplitudes of vibrations,since the computation was not able to get out of numerous loops of high-order internal resonance occurring atsmall amplitudes. For such complex structures with large modal density, it is known that an impressive number of high-orderinternal resonance are possible creating an intricate web of loops. Internal resonance has already been observedand studied on a reduced order model of a full 3D fan blade [71]. Although the dynamics observed in this paperis interesting, the selected reduced order model is too simple to validate the conclusions. The authors usedindeed a reduced-order model based on two linear modes with computation of the quadratic and cubic termusing the STEP method.Also, a number of these loops occur on very short parameter range and are often not robust to e.g. addition ofdamping. It results that they can be interpreted as obstacle for full order simulation with continuation, that arein general not meaningful for the global observable dynamics in forced-damped case. A computation with only3 harmonics has been able to get rid of this limitation and attains large amplitude by filtering the higher orderinteractions. Hence this solution will be taken as reference, and is shown in Fig. 7, for a vibration amplitude upto 0.18 m. Here the reported amplitude corresponds to the displacement in the x direction, at the point shownin Fig. 5a. Given the slenderness of the blade, 0 . m vibration amplitude at this point corresponds to 0 . m atthe leading edge which is indeed a large magnitude with regard to usual applications. In this range, a softeningbehaviour is reported, in the line of usual results reported in different studies for fan blades, and showing theimportance of the quadratic nonlinearity and shell-like behaviour of this structure. XYZ
Fig. 6 – Second-order correction vector shapes ¯ a and ¯ b used in the reduction procedure for the first mode of the fanblade φ . The vector ¯ b presents a strong component of the mode φ , highlighting the important quadratic self-couplingterm g due to the curvature. In the table are listed the values of the parameters needed to construct the reduced-orderdynamics for the fundamental mode. Three different ROMs are built in order to test their ability to recover the backbone of the reference solution.The first ROM is composed of a single master mode. As shown in Fig. 7, it recovers directly the correct softeningbehaviour, as it could have been expected from previous theoretical derivations [15, 30]. This good predictionis an excellent result since many reduction method encounters difficulties in retrieving such a behaviour [47,69]. However, one can see that, when vibration amplitude at the forcing point reach values at about 0.1 m, thesingle NNM solution starts to depart from the reference solution, underlining that a new stiffening effect comesinto play and needs to be taken into account on order to obtain a ROM prediction closer to the reference.As shown in Table 2, modes 2 and 3 of the blade have close eigenfrequencies as compared to the first. Also,they are both not far from showing a 3:1 resonance relationship with the first eigenfrequency. Indeed, the ratioare respectively 2 .
78 and 3 .
36. Depending on the hardening/softening behaviour of each mode, it is likely thatthe exact fulfilment of 3:1 resonance might appear at larger amplitudes. Consequently, two other ROMs havebeen built, including these two additional modes. The first one is thus composed of master modes 1, 2 and 3,and uses second-order normal form, while the second one selects the same master modes, but uses third-ordernormal form with additional resonant terms corresponding to 3:1 resonance relationships. These two ROMspredicts backbone curves that are closer to the reference, showing that the addition of these supplementarymodes is meaningful for achieving a correct solution on a larger range of amplitudes. They nevertheless show irect normal form for reduced-order models of finite element nonlinear structures 25
85 86 87 88 89 90 91 92 93 94 ?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]
Fig. 7 – Backbone curve of the fundamental mode of the fan blade. Amplitude of vibrations (in m ) computed at the nodewhere the force vector is shown, see Fig. 5a, and in the same direction x . a slight departure from reference at large amplitude, underlining the fact that the full-order solution is againmore complex, with possibly other higher-order resonances appearing.3.3 Frequency-response curve with damping and forcingIn this last section, frequency-response functions (FRFs) are computed including damping and forcing. The twostructures of the previous sections, namely the straight beam and the fan blade, are used for illustrative purposes.The aim is to show how on can easily build ROM including damping and forcing from the DNF technique, usingthe guidelines provided in section 2.6. First the fan blade is used to highlight the efficiency of the dampingmodel proposed, by comparing the results when using only either mass-proportional or stiffness-proportionaldamping. Then the beam is considered in order to show the behaviour of the FRF when approaching the 3:1internal resonance, in the case of the second NNM. The frequency response curve of the fan blade investigated in section 3.2 are computed and comparison betweena full-order simulation and reduced-order models composed of a single master mode using DNF with dampingand forcing, is produced. The forcing is harmonic with forcing frequency Ω and magnitude F . It is alignedwith the direction x and located at the point shown in Fig. 5a, a few centimetres below the tip. To build theROM, the guidelines provided in section 2.6 are applied. A single master coordinate corresponding to mode 1is selected. A modal force corresponding to modal projection of the external forcing is added at the right-handside of the dynamics. A linear damping term, directly deduced from the Rayleigh damping coefficients, is takeninto account, as well as the additional nonlinear damping term involving the coefficient C , and aggregatingthe contributions of the damping factors of the slave modes. In the remainder, comparisons will be drawn outby considering this additional term or not.In the first case under study, only mass-proportional damping is considered. The amplitude of the forcing isset at F = 30 N , and the Rayleigh coefficients are such that ζ M = 3 E − ζ K = 0. In the case of only mass-proportional damping, the modal damping factors reads ξ p = ζ M / ω p , where ξ p is the nondimensional factorappearing when expressing the linear damping term as 2 ξ p ω p ˙ X p . Consequently this damping factor decreaseswith the frequency, which is not a correct representation of real structures where losses generally increases withfrequency. In this specific case, the slave modes are thus less and less damped and their influence on the masterlow-frequency mode that is directly excited is awaited to be negligible.
85 90 95 100 ^_‘abcdefghijklmnopqrstuvwx Fig. 8 – Frequency response functions of the fan blade subjected to harmonic forcing in the vicinity of the first bendingeigenfrequency. (a) Only mass-proportional Rayleigh damping is taken into account, with ζ M = 3 E − ζ K = 0.Amplitude of the forcing F = 30 N . Comparison between full-order model (purple), and ROM composed of a single mastermode computed with direct normal form, either with only linear damping term (DNF linear damping, blue curve), orby taking into account the nonlinear aggregated damping term (DNF nonlinear damping, yellow dashed curve). (b) Onlystiffness-proportional damping with ζ K = 2 E − ζ M = 0. Amplitude of forcing F = 2 N . (c) Only stiffness-proportionaldamping with ζ K = 1 E − ζ M = 0, amplitude of forcing F = 10 N . Fig. 8(a) shows the FRFs obtained and compares the full-order solution to two different ROMs, one whereonly the linear damping term ζ M ˙ R is taken into account in the ROM (DNF linear damping, blue curve), andanother one where the aggregated nonlinear damping term is added (DNF nonlinear damping, dashed yellowcurve). In this case the two ROMs reproduce very well the solution of the full-order model, and the differencebetween the two is negligible, meaning that the aggregated C R ˙ R is very small and can be easily neglected.This is the direct consequence of the choice of only mass-proportional damping and the fact that modal dampingratios are decreasing with frequencies.Fig. 8(b) and (c) shows two different cases where only stiffness-proportional damping is selected. This caseis much more interesting since it corresponds to a more physically relevant situation. Also, since the modaldamping factors are linearly increasing in this case, this means that the damping factors of the slave modesare more and more important, so that neglecting their effect in a ROM would lead to discrepancies. Fig. 8(b)shows a case of an important stiffness-proportional damping with ζ K = 2 E − ζ M = 0. The forcing is setto F = 2 N . One can now observe a very important difference between the FRFs computed by the two ROMs.Neglecting the aggregated nonlinear damping term leads to important overestimation of the amplitude of thesolution in the resonant region. On the other hand, the ROM with the nonlinear damping exactly reproducesthe amplitude of the full model, and enforces the structure to vibrate in an almost linear regime as it shouldbe. Fig. 8(c) shows a case with smaller value of damping, ζ K = 1 E − ζ M = 0, and larger amplitude of forcing: F = 10 N . Even in this case the behaviour of the full-order model is almost linear with a small bump closeto linear resonance. Interestingly, the prediction given by taking into account only the linear damping term iscompletely wrong and gives amplitudes that are huge as compared to the reference. This clearly shows thatneglecting the damping of the slave modes can lead to erroneous predictions. On the other hand, the nonlineardamping term exactly recovers the correct amplitude and is able to perfectly match the reference solution. Theseexamples shows that even in a system with complex geometry and 3D elements, the proposed ROMs are ableto compute with great accuracy the FRF with minimal models. irect normal form for reduced-order models of finite element nonlinear structures 27 In this last example, FRFs of the straight clamped-clamped beam considered in section 3.1 are considered.The idea is to show how the method behaves in case of internal resonance where there is the need to take intoaccount more than one master coordinate in the ROM. For that purpose, a harmonic forcing in the vicinity ofthe second eigenfrequency is considered. Two different values of forcing are selected. The first case is that of asmall amplitude of forcing, such that the FRF does not enter in the loop of the 3:1 resonance. In this case, onlyone master mode is awaited to produce a correct prediction. A second case of larger amplitude is then selectedso that the FRF enters the resonance loop where 3:1 with mode 4 is excited. In this case, ROMs built from twomaster coordinates will be selected. In each case the forcing is located at z = 0 .
275 from the clamp as in Fig. 1,and excites the beam in the x direction. Finally, only the most interesting case of only stiffness-proportionaldamping is considered, since it is closer to real situations and give rise to more important differences betweenmodels. Fig. 9 – (a) Frequency response functions of the clamped-clamped beam subjected to harmonic forcing in the vicinity ofthe second eigenfrequency. Only stiffness-proportional damping is considered with ζ K = 13 E − ζ M = 0. Amplitudeof the forcing F = 10 N . Full-order solution is compared to ROMs obtained with a single master coordinate (mode 2), withonly linear damping term (blue curve) and with additional nonlinear damping term (dashed yellow curve). (b) close-upview near the maximum amplitude. Fig. 9 shows the results obtained when selecting ζ K = 13 E − F = 10 N , such that the FRF is below theinternal resonance loop. The reference solution is obtained by taking into account 6 harmonics in the expansion.Two different ROMs are contrasted, both containing a single master coordinate corresponding to mode 2, andone with only the linear damping term while the second one takes into account the aggregated nonlinear dampingterm. As in the previous case, the difference between the two ROMs are important and taking into account thenonlinear damping term is very important in order to correctly predict the maximum amplitude of the FRF.Fig. 10 now considers the case of a smaller damping coefficient ζ K = 3 E − F = 30 N . The reference solution now enters the loop so that the 3:1 internal resonance with mode 4 is activated.Consequently, the reference solution has been computed with a larger number of harmonics (15), in order toachieve convergence. Three different ROMs are compared, all of them being built with two master coordinatescorresponding to modes 2 and 4, and second-order normal form. The first ROM contains only the linear viscousdamping term. Since the underlying conservative system with second-order normal form is able to perfectly wellreproduce the resonance loop of the 3:1 resonance, this solution also shows this typical feature, highlightingthat the dynamics is more complex and involves more than one mode. However, the damping is not importantenough, so that a slight departure is observed from the reference solution. A second ROM is built by addingnonlinear damping terms, following Eq. (37). However a simplification is introduced. In the summation terms Fig. 10 – (a) Frequency response functions of the clamped-clamped beam subjected to harmonic forcing in the vicinity ofthe second eigenfrequency. Stiffness-proportional damping is considered with ζ K = 3 E − ζ M = 0. Amplitude of theforcing F = 30 N . Full-order solution is compared to different ROMs obtained with a two master coordinates correspondingto modes 2 and 4. Reference solution obtained with 15 harmonics (purple). All the ROMs are obtained with two modesand second-order normal form. Blue curve: only linear damping term considered. Yellow dashed curve: full considerationof nonlinear damping including cross-coupling terms. Brown curve: nonlinear damping neglecting the cross-couplings. (b)close-up view in the 3:1 resonance loop. involving C pijk terms, only the self-nonlinear damping terms C pppp are considered. This means that only two termshave been added to the reduced dynamics: C R ˙ R for the equation of mode 2 (first master coordinate) and C R ˙ R for the equation of mode 4 (second master coordinate). The cross-coupling nonlinear damping termsinvolving C pijk with different indexes have been discarded. In this case, one can observe a better behaviour ascompared to the linear case, meaning that the two added terms convey important information on the damping ofthe slave modes, that cannot be neglected. However the result is still a bit different from the reference solution.Finally a ROM with two modes and considering all the nonlinear damping terms including the cross-couplingsone, is reproduced in Fig. 10, and shows an almost perfect comparison with the reference. This example showsthat these cross-coupling terms also convey an important information, needed in order to reproduce with greataccuracy the FRF. A nonlinear mapping for the derivation of reduced-order models for geometrically nonlinear structures discretisedby the FE method, has been presented. It relies on the normal form theory and presents a firm theoreticalbackground in order to cope with nonlinear dynamics phenomena such as internal resonance. The method allowsfor a direct computation of the normal form, from a FE discretisation. Computational details are provided anda full algorithm has been presented with related discussion on the non-intrusiveness of the method and codingadvice for practical implementation. In essence, the technique allows one to go from the physical space (nodes ofthe FE structure given by the original mesh) to an invariant-based span of the phase space, with a third-orderapproximation of the invariant manifolds (nonlinear normal modes) of the structure. It is a simulation-freemethod that directly computes the third-order approximation of the reduced dynamics, with no restriction onthe number of selected master modes.The main developments have been proposed for a correct estimation of the nonlinear force of the ROM, i.e. in a conservative framework. Two declinations of the method have been proposed: a second-order normal formand a third-order mapping, with, in each case, the associated reduced dynamics. The second-order normal formproposes a quadratic mapping approximating the underlying invariant manifold up to the second-order only. irect normal form for reduced-order models of finite element nonlinear structures 29
Consequently, the reduced dynamics (up to order three) contains all possible third- and higher-order internalresonances embedded, and has no quadratic nonlinearity if no second-order internal resonance are present. Thethird-order normal form distinguishes the trivial and non-trivial resonant terms. The main improvement consistsin offering a better approximation of the invariant manifolds, and a more simplified reduced dynamics with onlytrivially resonant monomial cubic terms included. On the other hand, the reduced dynamics contains less termssince the higher-orders have been neglected. Consequently the dynamical solutions of the second-order DNFcan produced more accurate results as highlighted in the examples.A methodology has been proposed in order to take into account external forcing and damping in the ROMso as to produce frequency-response functions (FRF). For the damping, the proposed technique relies on afirst-order approximation of the general result derived in [27]. It has been rewritten for Rayleigh damping law,the most commonly retained assumption to model losses in FE context. The main interest relies in proposing anaggregated damping force on the reduced dynamics, taking into account the losses of all the slave (neglected)modes. In the present version, the general formulas for taking properly into account the damping have beengiven only for the second-order DNF.The method has been applied to a clamped-clamped beam and a FE model of a fan blade. Backbone curveshave been extracted and discussed. In particular, as known from previous studies on normal form for modelorder reduction, the method predicts directly the correct hardening/softening behaviour. Also, in the beamcase, it has been shown that the ROM can recover internal resonance loops due to the fulfilment of an exactcommensurability of nonlinear frequencies. FRF have been computed, underlining the importance of the added,aggregated damping, as well as the very good match between full and reduced models.Future work will apply the methodology to different structures having also internal resonance betweeneigenfrequencies of the system, in order to test the predictions of the ROM with an increasing number of mastermodes. Improvements of the method can also be investigated. For example higher-order reduced dynamicscan be computed by pushing the normal form up to order five for better accuracy. Taking into account thefull damped version of the normal form as developed in [27] could also be helpful. Finally, different internalforces could also be taken into account by enlarging the application field to e.g. gyroscopic forces, electrostaticforces (for MEMs applications), non conservative follower forces (for fluid-structure interaction problems) orpiezoelectric couplings.
Conflict of interest
The authors declare that they have no conflict of interest.
Codes availability statement
The codes written to run most of the simulations presented in this paper can be available upon simple requestto the authors. Note that for all the simulations the open-source software code aster has been used. Futuredevelopments of this code should include the reduction method proposed in this paper as a meta-command.
Acknowledgements
The author A. Vizzaccaro is thankful to Rolls-Royce plc for the financial support. The author Y. Shen wishes tothank China Scholarship Council (No.201806230253). The author L. Salles is thankful to Rolls-Royce plc and theEPSRC for the support under the Prosperity Partnership Grant Cornerstone: Mechanical Engineering Science toEnable Aero Propulsion Futures, Grant Ref: EP/R004951/1. The author J. Blahoˇs thank the European UnionsHorizon 2020 Framework Programme research and innovation programme under the Marie Sklodowska-Curieagreement No 721865.
References
1. S. W. Shaw and C. Pierre. Non-linear normal modes and invariant manifolds.
Journal of Sound and Vibration , 150(1):170–173, 1991.2. A. Steindl and H. Troger. Methods for dimension reduction and their applications in nonlinear dynamics.
InternationalJournal of Solids and Structures , 38:2131–2147, 2001.3. M. P. Mignolet, A. Przekop, S. A. Rizzi, and S. M. Spottswood. A review of indirect/non-intrusive reduced order modelingof nonlinear geometric structures.
Journal of Sound and Vibration , 332:2437–2460, 2013.4. A. Roberts.
Model emergent dynamics in complex systems . SIAM, Philadelphia, 2014.5. G. Haller and S. Ponsioen. Nonlinear normal modes and spectral submanifolds: existence, uniqueness and use in modelreduction.
Nonlinear Dynamics , 86(3):1493–1534, 2016.6. S. W. Shaw and C. Pierre. Normal modes for non-linear vibratory systems.
Journal of Sound and Vibration , 164(1):85–124,1993.7. M. Amabili and C. Touz´e. Reduced-order models for non-linear vibrations of fluid-filled circular cylindrical shells: comparisonof pod and asymptotic non-linear normal modes methods.
Journal of Fluids and Structures , 23(6):885–903, 2007.8. S. Ponsioen, T. Pedergnana, and G. Haller. Automated computation of autonomous spectral submanifolds for nonlinearmodal analysis.
Journal of Sound and Vibration , 420:269 – 295, 2018.9. R. M. Rosenberg. The normal modes of nonlinear n-degree-of-freedom systems.
Journal of Applied Mechanics , 29:7–14,1962.10. A. F. Vakakis, L. I. Manevitch, Y. V. Mikhlin, V. N. Philipchuck, and A. A. Zevin.
Normal modes and localization innon-linear systems . Wiley, New-York, 1996.11. G. Kerschen, M. Peeters, J.C. Golinval, and A.F. Vakakis. Non-linear normal modes, part I: a useful framework for thestructural dynamicist.
Mechanical Systems and Signal Processing , 23(1):170–194, 2009.12. S. W. Shaw and C. Pierre. Normal modes of vibration for non-linear continuous systems.
Journal of Sound and Vibration ,169(3):85–124, 1994.13. A. H. Nayfeh and S. A. Nayfeh. On nonlinear normal modes of continuous systems.
Trans. ASME/Journal of Vibrationand Acoustics , 116:129–136, 1994.14. G. Rega, W. Lacarbonara, and A. H. Nayfeh. Reduction methods for nonlinear vibrations of spatially continuous systemswith initial curvature.
Solid Mechanics and its applications , 77:235–246, 2000.15. C. Touz´e, O. Thomas, and A. Chaigne. Hardening/softening behaviour in non-linear oscillations of structural systems usingnon-linear normal modes.
Journal of Sound and Vibration , 273(1-2):77–101, 2004.16. X. Liu and D. J. Wagg. Simultaneous normal form transformation and model-order reduction for systems of coupled nonlinearoscillators.
Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences , 475(2228):20190042, 2019.17. J.C. Slater. A numerical method for determining nonlinear normal modes.
Nonlinear Dynamics , 10(1):19–30, 1996.18. E. Pesheck, C. Pierre, and S. Shaw. A new Galerkin-based approach for accurate non-linear normal modes through invariantmanifolds.
Journal of Sound and Vibration , 249(5):971–993, 2002.19. D. Noreland, S. Bellizzi, C. Vergez, and R. Bouc. Nonlinear modes of clarinet-like musical instruments.
Journal of Soundand Vibration , 324(3):983 – 1002, 2009.20. F. Blanc, C. Touz´e, J.-F. Mercier, K. Ege, and A.-S. Bonnet Ben-Dhia. On the numerical computation of nonlinear normalmodes for reduced-order modelling of conservative vibratory systems.
Mechanical Systems and Signal Processing , 36(2):520– 539, 2013.21. L. Renson, G. Kerschen, and B. Cochelin. Numerical computation of nonlinear normal modes in mechanical engineering.
Journal of Sound and Vibration , 364:177 – 206, 2016.22. R. Lewandowski. Computational formulation for periodic vibration of geometrically nonlinear structures, part II: numericalstrategy and examples.
International Journal of Solids and Structures , 34:1949–1964, 1997.23. R. Arquier, S. Bellizzi, R. Bouc, and B. Cochelin. Two methods for the computation of nonlinear modes of vibrating systemsat large amplitudes.
Computers & Structures , 84(24):1565 – 1576, 2006.24. B. Cochelin and C. Vergez. A high order purely frequency-based harmonic balance formulation for continuation of periodicsolutions.
Journal of Sound and Vibration , 324(1):243 – 262, 2009.25. M. Peeters, R. Vigui´e, G. S´erandour, G. Kerschen, and J.C. Golinval. Non-linear normal modes, part II: toward a practicalcomputation using numerical continuation techniques.
Mechanical Systems and Signal Processing , 23(1):195–216, 2009.26. P. Apiwattanalunggarn, C. Pierre, and D. Jiang. Finite-element-based nonlinear modal reduction of a rotating beam withlarge-amplitude motion.
Journal of Vibration and Control , 9:235–263, 2003.27. C. Touz´e and M. Amabili. Non-linear normal modes for damped geometrically non-linear systems: application to reduced-order modeling of harmonically forced structures.
Journal of Sound and Vibration , 298(4-5):958–981, 2006.28. C. Touz´e. Normal form theory and nonlinear normal modes: theoretical settings and applications. In G. Kerschen, editor,
Modal Analysis of nonlinear Mechanical Systems , pages 75–160, New York, NY, 2014. Springer Series CISM courses andlectures, vol. 555.29. C. Touz´e, M. Amabili, and O. Thomas. Reduced-order models for large-amplitude vibrations of shells including in-planeinertia.
Computer Methods in Applied Mechanics and Engineering , 197(21-24):2030–2045, 2008.30. C. Touz´e and O. Thomas. Non-linear behaviour of free-edge shallow spherical shells: effect of the geometry.
InternationalJournal of Non-linear Mechanics , 41(5):678–692, 2006.31. T. Breunung and G. Haller. Explicit backbone curves from spectral submanifolds of forced-damped nonlinear mechanicalsystems.
Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences , 474(2213):20180083, 2018.irect normal form for reduced-order models of finite element nonlinear structures 3132. C. Touz´e, M. Vidrascu, and D. Chapelle. Direct finite element computation of non-linear modal coupling coefficients forreduced-order shell models.
Computational Mechanics , 54(2):567–580, 2014.33. A. Givois, A. Grolet, O. Thomas, and J.-F. De¨u. On the frequency response computation of geometrically nonlinear flatstructures using reduced-order finite element models.
Nonlinear Dynamics , 97(2):1747–1781, 2019.34. A. A. Muravyov and S. A. Rizzi. Determination of nonlinear stiffness with application to random vibration of geometricallynonlinear structures.
Computers & Structures , 81(15):1513–1523, 2003.35. R. Perez, X. Q. Wang, and M. P. Mignolet. Nonintrusive Structural Dynamic Reduced Order Modeling for Large Deforma-tions: Enhancements for Complex Structures.
Journal of Computational and Nonlinear Dynamics , 9(3), 2014.36. K. Kim, A. G. Radu, X.Q. Wang, and M. P. Mignolet. Nonlinear reduced order modeling of isotropic and functionallygraded plates.
International Journal of Non-Linear Mechanics , 49:100 – 110, 2013.37. M. Balmaseda, G. Jacquet-Richardet, A. Placzek, and D.-M. Tran. Reduced order models for nonlinear dynamic analysiswith application to a fan blade. In
Proceedings of ASME Turbo Expo 2019 , Paper No. GT2019-90813, June 17–21, Phoenix,Arizona, 2019.38. K. Kim, V. Khanna, X.Q. Wang, and M.P. Mignolet. Nonlinear reduced order modeling of flat cantilevered structures. In
Proceedings of the 50th Structures, Structural Dynamics, and Materials Conference , AIAA Paper AIAA-2009-2492., May4–7, Palm Springs, California, 2009.39. A. Vizzaccaro, A. Givois, P. Longobardi, Y. Shen, J.-F. De¨u, L. Salles, C. Touz´e, and O. Thomas. Non-intrusive reduced ordermodelling for the dynamics of geometrically nonlinear flat structures using three-dimensional finite elements.
ComputationalMechanics , 2020.40. J. J. Hollkamp and R. W. Gordon. Reduced-order models for non-linear response prediction: Implicit condensation andexpansion.
Journal of Sound and Vibration , 318:1139–1153, 2008.41. R. J. Kuether, B. J. Deaner, J. J. Hollkamp, and M. S. Allen. Evaluation of geometrically nonlinear reduced-order modelswith nonlinear normal modes.
AIAA Journal , 53(11):3273–3285, 2015.42. A. Frangi and G. Gobat. Reduced order modelling of the non-linear stiffness in MEMS resonators.
International Journalof Non-Linear Mechanics , 116:211 – 218, 2019.43. S. R. Idelsohn and A. Cardona. A reduction method for nonlinear structural dynamic analysis.
Computer Methods inApplied Mechanics and Engineering , 49(3):253 – 279, 1985.44. O. Weeger, U. Wever, and B. Simeon. On the use of modal derivatives for nonlinear model order reduction.
InternationalJournal for Numerical Methods in Engineering , 108(13):1579–1602, 2016.45. L. Wu and P. Tiso. Nonlinear model order reduction for flexible multibody dynamics: a modal derivatives approach.
Multibody Syst Dyn , 36:405–425, 2016.46. S. Jain, P. Tiso, J. B. Rutzmoser, and D. J. Rixen. A quadratic manifold for model order reduction of nonlinear structuraldynamics.
Computers & Structures , 188:80–94, 2017.47. J. B. Rutzmoser, D. J. Rixen, P. Tiso, and S. Jain. Generalization of quadratic manifolds for reduced order modeling ofnonlinear structural dynamics.
Computers & Structures , 192:196–209, 2017.48. G. Haller and S. Ponsioen. Exact model reduction by a slowfast decomposition of nonlinear mechanical systems.
NonlinearDynamics , 90:617–647, 2017.49. A. Vizzaccaro, L. Salles, and C. Touz´e. Comparison of nonlinear mappings for reduced-order modelling of vibrating struc-tures: normal form theory and quadratic manifold method with modal derivatives.
Nonlinear Dynamics , 2020.50. Y. Shen, N. B´ereux, A. Frangi, and C. Touz´e. Reduced order models for geometrically nonlinear structures: assessmentof implicit condensation in comparison with invariant manifold approach.
European Journal of Mechanics, A/Solids ,submitted, 2020.51. P. Krysl, S. Lall, and J. E. Marsden. Dimensional model reduction in non-linear finite element dynamics of solids andstructures.
International Journal for Numerical Methods in Engineering , 51(4):479–504, 2001.52. M. Amabili, A. Sarkar, and M. P. Pa¨ıdoussis. Reduced-order models for nonlinear vibrations of cylindrical shells via theproper orthogonal decomposition method.
Journal of Fluids and Structures , 18(2):227–250, 2003.53. F. Chinesta, P. Ladeveze, and E. Cueto. A short review on model order reduction based on proper generalized decomposition.
Archives of Computational Methods in Engineering , 18(4):395, 2011.54. L. Meyrand, E. Sarrouy, B. Cochelin, and G. Ricciardi. Nonlinear normal mode continuation through a proper generalizeddecomposition approach with modal enrichment.
Journal of Sound and Vibration , 443:444 – 459, 2019.55. Z. Veraszto, S. Ponsioen, and G. Haller. Explicit third-order model reduction formulas for general nonlinear mechanicalsystems.
Journal of Sound and Vibration , 468:115039, 2020.56. E. Pesheck, N. Boivin, C. Pierre, and S. Shaw. Nonlinear modal analysis of structural systems using multi-mode invariantmanifolds.
Nonlinear Dynamics , 25:183–205, 2001.57. A. Lazarus, O. Thomas, and J.-F. De¨u. Finite element reduced order models for nonlinear vibrations of piezoelectric layeredbeams with applications to NEMS.
Finite Elements in Analysis and Design , 49:35–51, 2012.58. C. Touz´e. A normal form approach for non-linear normal modes. Technical report, Publications du LMA, num´ero 156,(ISSN: 1159-0947, ISBN: 2-909669-20-3), 2003.59. C. Touz´e, O. Thomas, and A. Huberdeau. Asymptotic non-linear normal modes for large amplitude vibrations of continuousstructures.
Computers and Structures , 82(31-32):2671–2682, 2004.60. N. Boivin, C. Pierre, and S. Shaw. Non-linear normal modes, invariance, and modal dynamics approximations of non-linearsystems.
Nonlinear Dynamics , 8:315–346, 1995.61. E. Pesheck.
Reduced-order modeling of nonlinear structural systems using nonlinear normal modes and invariant manifolds .PhD thesis, University of Michigan, 2000.2 Alessandra Vizzaccaro et al.62. D. Jiang.
Nonlinear modal analysis based on invariant manifolds. Application to rotating blade systems . PhD thesis,University of Michigan, 2004.63. D. Jiang, C. Pierre, and S. Shaw. Nonlinear normal modes for vibratory systems under harmonic excitation.
Journal ofSound and Vibration , 288(4-5):791–812, 2005.64. Electricit´e de France. Finite element code aster
International Journal of Solids and Structures , 34:1925–1947, 1997.66. C. S. M. Sombroek, P. Tiso, L. Renson, and G. Kerschen. Numerical computation of nonlinear normal modes in a modalderivative subspace.
Computers & Structures , 195:34 – 46, 2018.67. J. Blahoˇs, A. Vizzaccaro, F. El Haddad, and L. Salles. Parallel harmonic balance method for analysis of nonlinear dynamicalsystems. In
Proceedings of ASME Turbo Expo 2020 , Paper No. GT2020-15392, Sep 21–25, London, UK, 2020.68. H. Farokhi and M. H. Ghayesh. Geometrically exact extreme vibrations of cantilevers.
International Journal of MechanicalSciences , 168:105051, 2020.69. M. Balmaseda, G. Jacquet-Richardet, A. Placzek, and D.-M. Tran. Reduced order models for nonlinear dynamic analysiswith application to a fan blade.
Journal of Engineering for Gas Turbines and Power , 142(4), 2020.70. J. B. Rutzmoser.
Model order reduction for nonlinear structural dynamics: simulation-free approaches . PhD thesis, Tech-nischen Universit¨at M¨unchen (TUM), 2018.71. N. Di Palma, A. Martin, F. Thouverez, and V. Courtier. Nonlinear harmonic analysis of a blade model subjected to largegeometrical deflection and internal resonance. In
Proceedings of ASME Turbo Expo 2019 , Paper No. GT2019-91213, June17–21, Phoenix, Arizona, 2019.72. A. Vizzaccaro, Y. Shen, L. Salles, and C. Touz´e. Model order reduction methods based on normal form for geometricallynonlinear structures: a direct approach. In proc. of Euromech Non-linear Dynamics Conference, ENOC 2020 , Lyon, July2020.
A Derivation of second-order tensor
In this appendix, we start by recalling the general formulas needed to compute the second-order tensors from the normal formin modal basis, namely a , b , γ . These have been first derived in [58,15], they are here expressed in a different setting allowingfor better comparison with the tensors obtained with the direct method.The derivation of the second order tensors in modal basis, is obtained by substituting the nonlinear mapping equationsfor x , y into the equations of motion for ˙ x , ˙ y , by replacing each ˙ R r (and ˙ S r ) term with the first order assumption ˙ R r = S r (and ˙ S r = − ω r R r ), and finally by balancing the second order terms. This process will lead to four equations for each couple ofindexes i, j , one for the terms multiplying R i R j , one for S i S j , one for R i S j , and one for S i R j . The equations obtained permitto derive a ij , b ij , γ ij , γ ji ; the system reads: a ij − ω i b ij − γ ij = , a ij − ω j b ij − γ ji = , − Ω a ij + ω j γ ij + ω i γ ji = g ij , Ω b ij + γ ij + γ ji = . (41)The solution of this system of equations in the unknown tensors can be written as: D ij a ij =((+ ω i + ω j ) I − Ω ) g ij , (42a) D ij b ij =2 g ij , (42b) D ij γ ij =(( − ω i + ω j ) I − Ω )2 g ij , (42c) D ij γ ji =((+ ω i − ω j ) I − Ω )2 g ij , (42d)where the left-hand side term D ij is the denominator of the coefficients. It makes appear the second order internal resonancerelationship (small denominator problem), and reads: D ij = ((+ ω i + ω j ) I − Ω )(( − ω i + ω j ) I − Ω ) . (43)If no second order internal resonance occurs, in Eqs. (42) the denominator matrix can be inverted and the expression for thetensors made explicit. With these expressions, one is then able to compute all the coefficients needed for the modal normal form.irect normal form for reduced-order models of finite element nonlinear structures 33In order to shed more light on the close relationship between the coefficients of the modal and the direct normal form, letus show on a specific term a sij of the second order tensor a ij how one can go from one form to another. The explicit expressionfor a sij can be deduced from the previous equation or found in [15]. It reads: a sij = ((+ ω i + ω j ) − ω s )((+ ω i + ω j ) − ω s )(( − ω i + ω j ) − ω s ) g sij , (44)Let us first rewrite this formula so that the denominator, which is of order four in ω , is split into two factors of ordertwo, reminding that this denominator has a physical meaning. Indeed, the two factors represents the two possible second-orderinternal resonances between the modes of index i, j, s . One obtains: a sij = (cid:18) ω i + ω j ) − ω s + 1( − ω i + ω j ) − ω s (cid:19) g sij . (45)This can be rewritten in a more compact form for the vector a ij for all s : a ij = 12 ((+ ω i + ω j ) I − Ω ) − g ij + 12 (( − ω i + ω j ) I − Ω ) − g ij . (46)It is now possible to express a in terms of K and M by using Eqs. (5) stating V T MV = I , and V T KV = Ω , so that: a ij = 12 ((+ ω i + ω j ) V T MV − V T KV ) − g ij + 12 (( − ω i + ω j ) V T MV − V T KV ) − g ij . (47)Each factor in brackets can be rewritten as: (cid:16) V T ( σ ij M − K ) V (cid:17) − = V − ( σ ij M − K ) − V − T (48)where for the sake of readability the term σ ij replaced the summation of eigenvalues ( ± ω i + ω j ). One can see that the fulleigenvector matrix appears twice in each factor. Consequently a simplification arises by premultiplying Eq. (47) by V , whichalso leads to make appear ¯ a ij since ¯ a ij = Va ij . Also, by noticing that V − T g ij = G ( φ i , φ j ), the post-multiplied term can besimplified, so that a final expression involving only terms from the physical basis is obtained: ¯ a ij = 12 ((+ ω i + ω j ) M − K ) − G ( φ i , φ j ) + 12 (( − ω i + ω j ) M − K ) − G ( φ i , φ j ) . (49)Lastly, by defining the two vectors ¯Zd ij , ¯Zs ij the expressions given in text can be finally retrieved. Moreover, by applying asimilar procedure to the other tensors b , γ , one can show that the four tensors can all be expressed as linear combinations of ¯Zd ij , ¯Zs ij , hence only two operations are needed to derive four tensors.Another possible procedure to express the tensors in physical coordinates, is now described in order to show why the splittingoperation done in Eq. (45) has permitted to drastically reduce the computational cost. In fact, to derive the expression of thetensors in physical basis, the splitting operation is not strictly needed and another expression for them can be found startingfrom Eq. (44). From Eq. (44), the expression for a ij in compact form reads: a ij = (( − ω i + ω j ) I − Ω ) − ((+ ω i + ω j ) I − Ω ) − ((+ ω i + ω j ) I − Ω ) g ij . (50)Again, its equivalent in physical coordinates is the tensor ¯ a ij = Va ij , and the equivalent of the forcing term in physicalcoordinates is G ( φ i , φ j ) = V − T g ij . By introducing these two equivalences into Eq. (50) one obtains: ¯ a ij = V (( − ω i + ω j ) I − Ω ) − ((+ ω i + ω j ) I − Ω ) − ((+ ω i + ω j ) I − Ω ) V T G ( φ i , φ j ) . (51)This equation still presents two matrices in modal coordinates: the full eigenvectors matrix V and the full eigenvalues matrix Ω . To eliminate these two matrices and make K and M appear, one last operation is needed. Introducing the matrix: O = M − K , (52)one has that the relationship between O and Ω is: Ω = V − O V . (53)Finally, by substituting into Eq. (51), the last expression for Ω and simply V − V for I , the direct equation for a ij can beobtained: ¯ a ij = (( − ω i + ω j ) I − O ) − ((+ ω i + ω j ) I − O ) − ((+ ω i + ω j ) I − O ) M − G ( φ i , φ j ) . (54)This expression of a has been given in [72] for the case of equal indexes ii . Interestingly, it bears a strong resemblance with theexpression given in [55] where the ROM is built using a direct approach based on spectral submanifold (SSM).4 Alessandra Vizzaccaro et al.Comparing this expression for a with the one obtained after the splitting operation (Eq. (49)), and bearing in mind thatthey are completely equivalent, one sees that from a computational point of view, Eq. (54) is much more expensive. In Eq. (49),two linear systems have to be solved to find ¯Zd and ¯Zs and once they are found, the four tensors a , b , γ can be obtained. InEq. (54), three linear systems have to be solved, one for each − M needed to build the matrix O .This cost of this operation is equivalent to the cost of solving a linear system for each degree of freedom of the structure, makingprocedure infeasible for any engineering application. For this reason, the splitting operation that brought to Eq. (49), althoughnot strictly necessary from a theoretical point of view, is instead crucial from a computational one. B Derivation of third-order tensor
In this appendix, a similar approach to that of Appendix A is used to generate the expressions of the third order tensors r , u , µ , ν .We start by recalling the general expressions from the modal approach already developed in [58,15] .In absence of second order internal resonances, no second order terms are present in the reduced dynamics equations andthe first order assumptions ( ˙ R r = ˙ S r and ˙ S r = − ω r R r ), can still be used in the derivation of the third order tensors. Thesystem resulting from the balance of each third order term for a given triplet of indexes i, j, k reads: µ ijk − u ijk − u jki − u kij = , r ijk − ν kij − ω i u jki − ω j u ijk = , r ijk − ν jki − ω i u kij − ω k u ijk = , r ijk − ν ijk − ω j u kij − ω k u jki = ,ω i ν ijk + ω j ν jki + ω k ν kij − Ω r ijk = 3 h ijk + A ijk + A jki + A kij , − ν jki − ν kij + 3 ω i µ ijk − Ω u ijk = B ijk , − ν ijk − ν kij + 3 ω j µ ijk − Ω u jki = B jki , − ν ijk − ν jki + 3 ω k µ ijk − Ω u kij = B kij . (55)The expressions resulting from the solution of this system are much more lengthy than those of second order tensors but theypossess a similar structure. For instance, the expression for the tensor r is in the following form: D (8) ijk r ijk = Q (6) ijk ω j ω k B ijk + Q (6) jki ω k ω i B jki + Q (6) kij ω i ω j B kij + P (6) ijk ( A ijk + A jki + A kij + 3 h ijk ) . (56)The numerator matrices Q (6) and P (6) that multiplies the forcing terms A and B are in the form: P (6) ijk = p (6) ijk I + p (4) ijk Ω + p (2) ijk ( Ω ) + ( Ω ) , (57) Q (6) ijk = q (6) ijk I + q (4) ijk Ω + q (2) ijk ( Ω ) + ( Ω ) , (58)with p ( × ) ijk and q ( × ) ijk , both × -order polynomials in the sole parameters ω i , ω j , ω k . The numerator matrices varies for the differentthird order tensors. Conversely, the denominator matrix is the same for each tensor and reads: D (8) ijk = ((+ ω i + ω j + ω k ) I − Ω )(( − ω i + ω j + ω k ) I − Ω )((+ ω i − ω j + ω k ) I − Ω )((+ ω i + ω j − ω k ) I − Ω ) . (59)Once again, the denominator matrix carries a physical meaning in that every possible third order internal resonance betweenthe frequency appears in it.Following a similar approach than that on second order tensors, it is convenient to rewrite the expressions for the third ordertensors in such a way that the denominator is split. Taking a combination of indexes ijk with no internal resonance betweentheir eigenvalues, we can retrieve the basic vectors that represents the equivalent of ¯Zs and ¯Zd for the third order tensors inmodal basis. If one defines: Z0 ijk = ((+ ω i + ω j + ω k ) I − Ω ) − (cid:0) A ijk + A jki + A kij + 3 h ijk − ω j ω k B ijk − ω k ω i B jki − ω i ω j B kij (cid:1) (60a) Z1 ijk = (( − ω i + ω j + ω k ) I − Ω ) − (cid:0) A ijk + A jki + A kij + 3 h ijk − ω j ω k B ijk + ω k ω i B jki + ω i ω j B kij (cid:1) (60b) Z2 ijk = ((+ ω i − ω j + ω k ) I − Ω ) − (cid:0) A ijk + A jki + A kij + 3 h ijk + ω j ω k B ijk − ω k ω i B jki + ω i ω j B kij (cid:1) (60c) Z3 ijk = ((+ ω i + ω j − ω k ) I − Ω ) − (cid:0) A ijk + A jki + A kij + 3 h ijk + ω j ω k B ijk + ω k ω i B jki − ω i ω j B kij (cid:1) (60d)irect normal form for reduced-order models of finite element nonlinear structures 35it is possible to demonstrate that the third order tensors of modal normal form are equivalently expressed as a linear combinationof the above defined vectors and their expressions read: r ijk = Z0 ijk + Z1 ijk + Z2 ijk + Z3 ijk
12 (61a) u ijk = − Z0 ijk − Z1 ijk + Z2 ijk + Z3 ijk ω j ω k (61b) µ ijk = − (+ ω i + ω j + ω k ) Z0 ijk + ( − ω i + ω j + ω k ) Z1 ijk + (+ ω i − ω j + ω k ) Z2 ijk + (+ ω i + ω j − ω k ) Z3 ijk ω i ω j ω k (61c) ν ijk = +(+ ω i + ω j + ω k ) Z0 ijk − ( − ω i + ω j + ω k ) Z1 ijk + (+ ω i − ω j + ω k ) Z2 ijk + (+ ω i + ω j − ω k ) Z3 ijk ω i (61d)It is then straightforward to express Eqs. (60) in physical basis by using the same strategy that we applied to the secondorder tensors. Also in this case, there is only one linear system to solve to find the four independent tensors relative to the triplet i, j, k and from them all the eight third order tensors ¯ r , ¯ u , ¯ µ , ¯ ν can be derived. As a final remark, we want to point out thatan expression in physical coordinates of Eq. (56) could have been also derived, by replacing the matrix Ω with its equivalent V − O V but in that case the computational cost of building all the tensors ¯ r , ¯ u , ¯ µ , ¯ ν would have been even higher than thatof the second order one. C Reduced-order models for the beam case
This appendix gives the detailed equations governing the reduced-order dynamics in the beam case, in order to better understandwhich terms are present in the different ROMs, and which one are added when taking internal resonances into account. Thediscussion is restricted to the case of the the first NNM, where the main interesting feature is the presence of a 5:1 internalresonance (involving the nonlinear frequencies) between mode 1 and mode 3. Consequently ROMS have been built includingthese two modes as master coordinates. In the remainder, R and R will thus refer to the normal coordinates linked respectivelyto mode 1 and mode 3.The simplest case is that of the ROM composed of a single master coordinate R . The reduced dynamics reads, followingEq. (16): ¨ R + ω R + ( A + h ) R + B R ˙ R = 0 . (62)The second case discussed in section 3.1 is that of a third-order normal form composed of modes 1 and 3, without takinginto account any other terms due to the presence of an internal resonance. In this case the reduced dynamics deduces fromEq. (9) and directly reads:¨ R + ω R + ( A + h ) R + B R ˙ R + (cid:0) A + A + A + 3 h (cid:1) R R + B R ˙ R + (cid:0) B + B (cid:1) ˙ R R ˙ R = 0 , (63a)¨ R + ω R + ( A + h ) R + B R ˙ R + (cid:0) A + A + A + 3 h (cid:1) R R + B R ˙ R + (cid:0) B + B (cid:1) ˙ R R ˙ R = 0 . (63b)These equations does not contain any invariant-breaking terms. Consequently they are unable to create a coupling betweenthe invariant manifolds dynamics of mode 1 and 3. Since the ROM wants to follow the backbone curve of mode 1, adding thesecond equation on mode 3 without any invariant-breaking term is in this case meaningless, and this ROM will produce exactlythe same backbone as Eq. (62). In order to counteract this property, the resonant monomials corresponding to the 1:3 internalresonance can be added to these equations, leading to the ROM with two modes and inclusion of new resonant terms. Thereduced dynamics of this ROM now reads:¨ R + ω R + ( A + h ) R + B R ˙ R + (cid:0) A + A + A + 3 h (cid:1) R R + B R ˙ R + (cid:0) B + B (cid:1) ˙ R R ˙ R + (cid:0) A + A + A + 3 h (cid:1) R R + B R ˙ R + (cid:0) B + B (cid:1) R ˙ R ˙ R = 0 , (64a)¨ R + ω R + ( A + h ) R + B R ˙ R + (cid:0) A + A + A + 3 h (cid:1) R R + B R ˙ R + (cid:0) B + B (cid:1) R ˙ R ˙ R + ( A + h ) R + B R ˙ R = 0 . (64b)As explained in the main text, the only extra terms added to this ROM are the third line in each oscillator equation. ForEq. (64a), the cubic resonant monomial corresponding to 1:3 resonance is R R . This term needs to be accompanied with allhis ancillary ones conveying the same resonance but with the velocity terms ˙ R and ˙ R , resulting in the addition of the other6 Alessandra Vizzaccaro et al.resonant monomials R ˙ R and R ˙ R ˙ R . For the second equation (64b), the resonant monomial implying displacements onlyis R , which must also be complemented with all its companion terms sharing the same frequency resoannce and involcingvelocities, thus the R ˙ R term. One can note that R on the second equation is an invariant-breaking term. Consequently thissystem is prone to activate the resonance and break the invariance property of the single-mode manifold.The last system investigated includes as extra terms those corresponding to the 1:1 resonance. This final two-dofs reduceddynamics finally reads:¨ R + ω R + ( A + h ) R + B R ˙ R + (cid:0) A + A + A + 3 h (cid:1) R R + B R ˙ R + (cid:0) B + B (cid:1) ˙ R R ˙ R + (cid:0) A + A + A + 3 h (cid:1) R R + B R ˙ R + (cid:0) B + B (cid:1) R ˙ R ˙ R + +( A + h ) R + B R ˙ R = 0 , (65a)¨ R + ω R + ( A + h ) R + B R ˙ R + (cid:0) A + A + A + 3 h (cid:1) R R + B R ˙ R + (cid:0) B + B (cid:1) R ˙ R ˙ R + (cid:0) A + A + A + 3 h (cid:1) R R + B R ˙ R + (cid:0) B + B (cid:1) ˙ R R ˙ R + ( A + h ) R + B R ˙ R = 0 ..