Arthur Givois
email: arthur.givois@ensam.eu

Aurélien Grolet

Olivier Thomas

Jean-François Deü

J.-F Deü

On the frequency response computation of geometrically nonlinear flat structures using reduced-order finite element models

Keywords: Geometric nonlinearities, Non-intrusive stiffness evaluation procedure, Reduced-order finite element model, Continuation method

This paper presents a general methodology to compute nonlinear frequency responses of flat structures subjected to large amplitude transverse vibrations, within a finite element context. A reduced-order model (ROM) is obtained by an expansion onto the eigenmode basis of the associated linearized problem, including transverse and in-plane modes. The coefficients of the nonlinear terms of the ROM are computed thanks to a non-intrusive method, using any existing nonlinear finite element code. The direct comparison to analytical models of beams and plates proves that a lot of coefficients can be neglected and that the in-plane motion can be condensed to the transverse motion, thus giving generic rules to simplify the ROM. Then, a continuation technique, based on an asymptotic numerical method and the harmonic balance method, is used to compute the frequency response in free (nonlinear mode computation) or harmonically forced vibrations. The whole procedure is tested on a straight beam, a clamped circular plate and a free perforated plate for which some nonlinear modes are computed, including internal resonances. The convergence with harmonic numbers and

Introduction

Recent advances in non-intrusive reduced-order finite element modeling of geometrically nonlinear structures offer new perspectives for nonlinear prediction in structural computation. These geometrical nonlinearities are due to large transverse displacements of thin structures which induce a non-negligible stretching in the plane of the structure known as membranebending coupling phenomena. They are involved in a wide range of applications, among others: microelectromechanical systems (MEMS) developments [START_REF] Chen | Direct observation of coherent energy transfer in nonlinear micromechanical oscillators[END_REF][START_REF] Thomas | Efficient parametric amplification in MEMS with integrated piezoelectric actuation and sensing capabilities[END_REF][START_REF] Dezest | Wafer-scale fabrication of self-actuated piezoelectric nanoelectromechanical resonators based on lead zirconate titanate (PZT)[END_REF], energy harvesting [START_REF] Quinn | Energy harvesting from impulsive loads using intentional essential nonlinearities[END_REF], vibroacoustics of music instruments [START_REF] Ducceschi | Modal approach for nonlinear vibrations of damped impacted plates: application to sound synthesis of gongs and cymbals[END_REF][START_REF] Monteil | Identification of mode couplings in nonlinear vibrations of the steelpan[END_REF], dynamics of blades and rotors [START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: application to a simplified model[END_REF], and aerospace engineering [START_REF] Renson | Complex dynamics of a nonlinear aerospace structure: numerical continuation and normal modes[END_REF]. In a finite element context, since geometrical nonlinearities couple all the displacement degrees of freedom, the computation of models with a large number of variables is found to be a tough challenge, for which strategies for model reduction have been the focus of sustained efforts since the 1980s [START_REF] Noor | Recent advances in reduction methods for nonlinear problems[END_REF][START_REF] Slaats | Model reduction tools for nonlinear structural dynamics[END_REF]. In vibration analyses, the expansion of displacements on a reduced basis of vibration eigenmodes is an efficient and well-known approach to study structural dynamics at low and middle frequencies. This article addresses the use of this kind of modal reduced-order models (ROMs) to compute the frequency response of geometrically nonlinear thin structures.

Modal ROMs are the most used technique to solve analytical models of thin structures like beams, plates and shells. Nonlinear models based on the Von Kármán assumptions [START_REF] Von Kármán | Festigkeitsprobleme im Maschinenbau[END_REF] are well known and able to capture the membrane-bending coupling. The interested reader can refer to [START_REF] Nayfeh | Nonlinear Oscillations[END_REF] for an overview of the most popular cases, [START_REF] Chu | Influence of large amplitudes on free flexural vibrations of rectangular elastic plates[END_REF][START_REF] Thomas | Geometrically non-linear flexural vibrations of plates: in-plane boundary conditions and some symmetry properties[END_REF] for plate models, [START_REF] Woinowski-Krieger | The effect of axial force on the vibration of hinged bars[END_REF] for beam models and [START_REF] Eisley | Nonlinear vibration of beams and rectangular plates[END_REF] that shows the similarities between beam and plate models. Since the coefficients of the ROM depend on the structure's geometries as well as its boundary conditions, a huge literature is available for a large set of particular cases. Among others, the reader can refer to [START_REF] Ho | Non-planar, non-linear oscillations of a beam. part II: Free motions[END_REF] for beams with doubly hinged ends, to [START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: application to a simplified model[END_REF] for an assembly of clamped/free boundary conditions with cyclic symmetry, to [START_REF] Sridhar | Non-linear resonances in the forced responses of plates, part II: asymmetric responses of circular plates[END_REF] for clamped circular plates, to [START_REF] Touzé | Asymmetric non-linear forced vibrations of free-edge circular plates, part I: theory[END_REF] for free circular plates and to [START_REF] Ducceschi | Nonlinear dynamics of rectangular plates: investigation of modal interaction in free and forced vibrations[END_REF] for rectangular plates with several boundary conditions. In those analytical models, the in-plane inertia is usually neglected, which enables to use a mixed formulation, in which a transverse motion (described with a displacement field) is coupled to the in-plane motion described by a force unknown: an Airy stress function in the case of plates [START_REF] Thomas | Geometrically non-linear flexural vibrations of plates: in-plane boundary conditions and some symmetry properties[END_REF] and the axial force for beams [START_REF] Woinowski-Krieger | The effect of axial force on the vibration of hinged bars[END_REF].

To address complex geometries out of the scope of analytical models, modal ROMs can be obtained from numerical models. It can be first achieved by explicit computations, based on finite element discretization [START_REF] Capiez-Lernout | Computational stochastic statics of an uncertain curved structure with geometrical nonlinearity in three-dimensional elasticity[END_REF][START_REF] Touzé | Direct finite element computation of non-linear modal coupling coefficients for reduced-order shell models[END_REF]. p-type finite elements (high-order elements of large size) have also been applied to geometrically nonlinear vibrations of beams and plates in a family of papers [START_REF] Ribeiro | Non-linear vibration of beams with internal resonance by the hierarchical finite element method[END_REF][START_REF] Stoykov | Periodic geometrically nonlinear free vibrations of circular plates[END_REF]. However, with the increasing spread of commercial finite element codes, nonintrusive methods-i.e., using the finite element code as a "black box" without modifying it-have been proposed. Two main families of methods have been introduced, both relying on static finite element computations, since the geometrical nonlinearities are concentrated in the stiffness part of the model (the internal force vector). On the one hand, one can impose static loads as linear combination of eigenmodes [START_REF] Mcewan | A combined modal/finite element analysis technique for the dynamic response of a non-linear beam to harmonic excitation[END_REF], and on the other hand, static displacement fields, also as linear combination of eigenmodes [START_REF] Muravyov | Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures[END_REF], can be prescribed. These methods-and particularly the second one called "STiffness Evaluation Procedure" (STEP) on which this paper focuses-have retained a large attention in the past two decades (see, e.g., the review [START_REF] Mignolet | A review of indirect/non-intrusive reduced-order modeling of nonlinear geometric structures[END_REF]), motivated by the perspective to compute nonlinear responses of structures with arbitrary geometries. More recently, a variant of the STEP has been proposed in [START_REF] Perez | Nonintrusive structural dynamic reduced-order modeling for large deformations: enhancements for complex structures[END_REF], which enables to reduce the number of static computations by using the tangent stiffness matrix. In the past few years, those reduced-order models have been validated by combining computational and experimental investigations [START_REF] Murthy | Uncertainty-based experimental validation of nonlinear reduced-order models[END_REF][START_REF] Claeys | Multi-harmonic measurements and numerical simulations of nonlinear vibrations of a beam with non-ideal boundary conditions[END_REF]. They have been applied to problems of increasing difficulty, such as fatigue life estimation of thin gauges [START_REF] O'hara | Modeling vibratory damage with reduced-order models and the generalized finite element method[END_REF], calibration of a perforated plate [START_REF] Ehrhardt | Finite element model calibration of a nonlinear perforated plate[END_REF] and dynamical simulations of composite plates [START_REF] Mignolet | A review of indirect/non-intrusive reduced-order modeling of nonlinear geometric structures[END_REF][START_REF] Kim | Nonlinear reduced-order modeling of isotropic and functionally graded plates[END_REF], of a nine-bay panel structure [START_REF] Perez | Nonintrusive structural dynamic reduced-order modeling for large deformations: enhancements for complex structures[END_REF] or of an electromechanical nanobeam [START_REF] Lazarus | Finite elements reduced order models for nonlinear vibrations of piezoelectric layered beams with applications to NEMS[END_REF].

In the above-cited methods, one major issue lies in the selection of the linear modal basis to perform nonlinear structural calculations. In-plane modes must be included to account for bending-membrane coupling, whereas only transverse modes constitute most of the times the modal basis in linear dynamic analysis since in-plane modes are located at higher frequencies than transverse ones. To capture these effects, several methods have been proposed. The most natural one consists in keeping in the ROM basis several in-plane modes in addition to transverse one. Since numerous in-plane modes are necessary, some authors propose to reduce the ROM size by condensing the inplane modes into the transverse oscillators by neglecting the membrane inertia [START_REF] Mignolet | A review of indirect/non-intrusive reduced-order modeling of nonlinear geometric structures[END_REF][START_REF] Hollkamp | Reduced-order models for nonlinear response prediction: implicit condensation and expansion[END_REF][START_REF] Kim | Nonlinear reduced order modeling of flat cantilevered structures[END_REF]. In the case of flat/transversely symmetric structures, it is also possible to reduce the number of computed coefficients of the ROM by observing that a lot of them are zero. This sparsity property has been mentioned in [START_REF] Mignolet | A review of indirect/non-intrusive reduced-order modeling of nonlinear geometric structures[END_REF][START_REF] Nash | Nonlinear structural dynamics by finite element modal synthesis[END_REF], but no precise proof is available, to the knowledge of the authors. It is also possible to select the most relevant in-plane modes by using proper/smooth orthogonal decomposition (P/SOD) [START_REF] Rizzi | System identification-guided basis selection for reduced-order nonlinear response analysis[END_REF][START_REF] Przekop | Alternative modal basis selection procedures for reduced-order nonlinear random response simulation[END_REF]. These strategies are efficient but rely on the drawback of requiring the preliminary time integration of the full-scale model. Some authors propose strategies to naturally compute relevant in-plane displacement field with prescribed modal loads, with so-called implicit condensation with expansion approach [START_REF] Hollkamp | Reduced-order models for nonlinear response prediction: implicit condensation and expansion[END_REF] or dual modes strategies [START_REF] Kim | Nonlinear reduced-order modeling of isotropic and functionally graded plates[END_REF]. These methods have the drawback that the ROM coef-ficients are sensitive to the size of the transverse modal basis and the amplitudes of the prescribed loads [START_REF] Kuether | Evaluation of geometrically nonlinear reduced-order models with nonlinear normal modes[END_REF]. A last family of method is based on the addition of a modal derivative subspace to the transverse linear basis [START_REF] Idelsohn | A reduction method for nonlinear structural dynamic analysis[END_REF]. It enables to capture efficiently the relevant in-plane displacement associated with a transverse linear mode [START_REF] Sombroek | Numerical computation of nonlinear normal modes in a modal derivative subspace[END_REF]. Since the size of the subspace grows quadratically with the number of investigated eigenmodes, the concept of quadratic manifold has been proposed to minimize the number of unknowns [START_REF] Jain | A quadratic manifold for model order reduction of nonlinear structural dynamics[END_REF]. The computation of the modal derivatives with a commercial finite elements code is based on a finite difference scheme that requires the extraction of the tangent stiffness matrix, which renders the strategy more intrusive [START_REF] Rutzmoser | Generalization of quadratic manifolds for reduced order modeling of nonlinear structural dynamics[END_REF]. An asymptotic expansion of a full finite element discretization proposed in [START_REF] Boumediene | Nonlinear forced vibration of damped plates by an asymptotic numerical method[END_REF] can lead to build a different reduction basis than the modal one which gives satisfactory results [START_REF] Boumediene | Nonlinear forced vibration of damped plates coupling asymptotic numerical method and reduction models[END_REF].

When the reduced-order model has been built, the dynamical nonlinear responses of the structures at low and middle frequencies can be investigated. Most of the early studies about modal ROMs investigate transient responses and rely on direct time integration [START_REF] Mignolet | A review of indirect/non-intrusive reduced-order modeling of nonlinear geometric structures[END_REF][START_REF] Kim | Nonlinear reduced-order modeling of isotropic and functionally graded plates[END_REF][START_REF] Hollkamp | Reduced-order models for nonlinear response prediction: implicit condensation and expansion[END_REF][START_REF] Przekop | Alternative modal basis selection procedures for reduced-order nonlinear random response simulation[END_REF]. Conversely, it is possible to investigate the frequency response by computing the nonlinear modes of the structure, because they constitute the skeleton of the dynamics, since the harmonically forced responses lies around [START_REF] Kerschen | Nonlinear normal modes, part I: a useful framework for the structural dynamicist[END_REF]. The so-called nonlinear normal modes-(NNMs) can be represented as frequency energy plots (or amplitude-frequency plots, also called "backbone curves" [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]) and computed by various analytical [START_REF] Lamarque | An upper bound for validity limits of asymptotic analytical approaches based on normal form theory[END_REF] or numerical [START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF] methods. Among these last, a method combining shooting and pseudo-arclength continuation proposed in [START_REF] Peeters | Nonlinear normal modes, part II: toward a practical computation using numerical continuation techniques[END_REF] enables to compute NNMs from reduced-order finite element models: with a homemade finite element beam formulation [START_REF] Sombroek | Numerical computation of nonlinear normal modes in a modal derivative subspace[END_REF][START_REF] Sombroek | Bridging the gap between nonlinear normal modes and modal derivatives[END_REF] or some finite element commercial codes [START_REF] Kuether | A numerical approach to directly compute nonlinear normal modes of geometrically nonlinear finite element models[END_REF], the energy is computed by injecting the result of the ROM in the full FE model. Among the huge literature described above, the present article addresses the computation of the nonlinear frequency response of thin structures using a continuation method and finite element modal reduced-order models, obtained by the STEP. Within this subject, we focus on two main points: (i) a precise validation of the STEP using reference analytical models of flat structures (beams and plates) and (ii) the modal convergence of the ROMs, especially when internal resonance is involved. This study is restricted to flat structures since simple analytical models are available and in order to provide reference benchmarks that can be used to test algorithms, potentially adapted to more complex structures such as shells or thin structures with stiffeners. Moreover, the flat geometry gives a natural separation of in-plane and transverse modes in terms of the ROM convergence, which is especially addressed. Since inplane modes are condensed into the transverse motion, we first analyze the convergence of the ROM in terms of in-plane modes at the coefficient level and its convergence in terms of transverse modes is studied in a second step, with frequency response computations. For the STEP validation, several topics are proposed: (a) a systematic validation of the ROM coefficients value with reference analytical models; (b) the proof of the sparsity structure of the coefficients matrices of the ROM and its consequences on the computation time saving; (c) the proof of the robustness of the STEP in terms of static prescribed displacement amplitude; (d) an heuristic rule to filter significant ROM coefficients to reduce the ROM size. We show that the STEP is perfectly validated on two reference cases (a clampedclamped beam and a clamped circular plate) for which analytical models are known. Finally, the method is applied to a perforated plate, for which no analytical solutions are known. Convergence results are obtained, showing that very reduced ROM, composed of a few transverse oscillators coupled by cubic terms only, is efficient to compute nonlinear modes even in the case of internal resonances and vibration localization.

Reference analytical models of flat structures

The purpose of this section is to recall classical analytical models for flat structures, firstly to serve as a reference for the validations of the STEP proposed in the upcoming sections and secondly to systematically prove that a lot of coefficients of the modal reducedorder models are null. Here, only models relying on the von Kármán theory have been selected since they are valid for most of practical cases of 2D (plate) structures and 1D (beam) structures with immovable ends in the axial direction. In those cases, large rotation "geometrically exact" models give the same results in the same range of amplitude of vibration and are thus not necessary. Von Kármán models can be viewed as an intermediate and consistent model between linear models and full large rotation models [START_REF] Ciarlet | A justification of the von-Kármán equations[END_REF][START_REF] Millet | A classification of thin plate models by asymptotic expansion of non-linear three-dimensional equilibrium equations[END_REF]. On the contrary, large rotation models have to be considered in the other cases (beams with ends free to move such as cantilever beams, see [START_REF] Lacarbonara | Refined models of elastic beams undergoing large in-plane motions: theory and experiment[END_REF][START_REF] Thomas | Hardening/softening behaviour and reduced order modelling of nonlinear vibrations of rotating cantilever beams[END_REF]), out of the scope of the article.

Straight beam model

In this section, we consider a non-prestressed beam of length L of uniform rectangular cross section of area S and second moment of area I , made in an homogeneous and isotropic material of Young's modulus E and density ρ. Only the bending motion is under concern, and w(x, t) and u(x, t) denote, respectively, the transverse and axial displacements at time t and abscissa x.

Following the von Kármán assumptions and with a linear elastic constitutive law, the equations of motions can be written:

ρ S ẅ + E I w ,x x x x -E S u ,x w ,x + 1 2 w 3 ,x ,x = p(x, t), (1a) ρ S ü -E S u ,x x + w ,x w ,x x = n(x, t). ( 1b
)
where

• ,x = ∂ • /∂ x.
All details about this formulation can be found in "Appendix A".

To obtain a generic formulation as well as to scale the different terms in the equations, the following dimensionless parameters are introduced:

w = w h , x = x L ū = L h 2 u, t = 1 L 2 E I ρ S t, p = L 4 E I h p, (2)
n = L 3 E Sh 2 n, ε = Sh 2 I , μ = I SL 2 . ( 3
)
where the scaling of the transverse displacement has been chosen as the thickness h of the beam. Then, Eqs. (1a, b) are rewritten as:

ẅ + w ,x x x x -ε u ,x w ,x + 1 2 w 3 ,x ,x = p(x, t), (4a) μ ü -(u ,x x + w ,x w ,x x ) = n(x, t), ( 4b
)
where the overbars have been dropped for simplification.

The bending and axial displacements are expanded onto a linear basis of, respectively, N w transverse modes Φ k and N u axial modes Ψ p :

w(x, t) = N w k=1 Φ k (x)q k (t) and u(x, t) = N u p=1 Ψ p (x)η p (t) (5)
with time dependencies governed by unknown modal coordinates q k and η p . The modes are chosen to satisfy the following two eigenproblems associated with the linear parts of Eqs. (4a, b):

Φ k,x x x x -ω 2 k Φ k = 0 ( 6 a ) Ψ p,x x -γ 2 p Ψ p = 0 (6b)
where ω k and γ p are the transverse and axial dimensionless eigenfrequencies. Finally, by multiplying (4a, b), respectively, by Φ k and Ψ p , integrating on the length of the beam and using the orthogonality properties, the following ROM is obtained, for all k = 1, . . .

N w , p = 1, . . . N u : qk + ω 2 k q k + ε N u p=1 N w i=1 C k pi η p q i + ε N w i=1 N w j=1 N w l=1 D k i jl q i q j q l = Q k (7a) μ ηp + γ 2 p η p - N w i=1 N w j=1 G p i j q i q j = E p (7b)
where the forcing terms Q k (t), E p (t) and the coefficients of the nonlinear terms C k pi , D k i jl and G p i j are defined in "Appendix A". Notice that the value of those coefficients depends on the mode normalization.

The above dimensionless form of the equations shows that because of the slender geometry of the beam (h/L 1), the axial displacement u is one order of magnitude smaller than the transverse displacement w, which leads to neglect the axial inertia because μ 1. In this case, if no axial external forcing is prescribed (n(x, t) = 0 ⇒ E p = 0), Eq. (7b) leads to write the axial modal coordinates modes as:

η p = 1 γ 2 p N w i=1 N w j=1 G p i j q i q j , ( 8
)
and to condense them in the transverse equation of motion (7a) to obtain, for all k = 1, . . . N w :

qk + ω 2 k q k + ε N w i=1 N w j=1 N w l=1 Γ k i jl q i q j q l = Q k (9)
with

Γ k i jl = D k i jl + N u p=1 G p i j C k pi γ 2 p . ( 10
)
Initial problem (7a, b) has been replaced by a reduced one composed of only N w transverse equations with cubic nonlinear terms only.

Plate model

In the same way as the beam, we consider here a thin plate of thickness h with a boundary of arbitrary shape, made in an homogeneous and isotropic material of Young's modulus E, Poisson's ratio ν and density ρ.

The transverse and in-plane displacements are written w(r, t) and u(r, t) at time t and position r. The vectors u and r are two-dimensional vectors of the mid-plane of the plate. The equations of motions can be written [START_REF] Bilbao | Conservative numerical methods for the full von Kármán plate equations[END_REF]:

ρh ẅ + D w -div (N∇w) = p(r, t), (11a) ρh ü -div N = n(r, t). ( 11b
)
where N is the membrane forces tensor, function of u and w, div and div are the scalar/vector divergence of a vector/tensor field, is the Laplacian of a scalar field,

D = Eh 3 /[12(1 -ν 2 )
] is the bending stiffness and p, n are transverse and in-plane external forces by unit area. By eliminating N between Eqs. ( 50) and (11a, b), one obtains a (u, w)-formulation written in terms of displacement unknowns only. Then, by expanding w and u on a linear basis of, respectively, N w transverse modes Φ k and N u axial modes Ψ p , one obtains exactly the same ROM than in the case of beams, namely Eqs. (7a, b), which has the form of a set of N w transverse and N u in-plane oscillators, coupled by quadratic and cubic terms of coefficients C k pi , D k i jl and G p i j . Finally, it can also be proven that for the same reasons, the in-plane unknowns can be condensed in the transverse equations to obtain exactly the same ROM [START_REF] Noor | Recent advances in reduction methods for nonlinear problems[END_REF], composed of only N w transverse equations with cubic nonlinear terms. All details are provided in "Appendix B".

Synthesis on analytical models for flat structures

The previous two sections have shown two main results, valid for any flat structure (beam or flat plate).

-First, in the case of a ROM written in term displacements unknowns only (the (u, w)-formulation for beams and the (u, w)-formulation for plates) expanded on the linear eigenmodes basis, the form of the ROM is the one of Eqs. (7a, b), with only a few nonzero nonlinear terms:

-quadratic terms in the transverse-q k oscillators that couple the transverse and in-plane coordinates: C k pi η p q k ; -quadratic terms in the in-plane-η p oscillators that involve only transverse coordinates: G p i j q i q j ; -cubic terms in the transverse-q k oscillators that couple only transverse coordinates: D k i jl q i q j q l . -Second, by neglecting the in-plane inertia, it is possible to condense the in-plane coordinates into the transverse oscillators, which leads to set (9), a ROM composed of N w transverse oscillators only, without additional oscillators for the in-plane motion, coupled by cubic terms only.

Another interesting result comes from the dimensionless form of the equations. For beams, if the axial inertia is neglected (μ 1), there is only one free parameter in the equations of motion (4a, b), which is ε = h 2 /i 2 with i = √ I /S the radius of gyration of the beam's cross section. If the transverse motion would have been scaled by i [ w = w/i in place of Eq. ( 2)], ε would be equal to 1, which means that the nonlinear dynamics of any straight slender beam with a uniform cross section depends only on its boundary conditions, whatever be the geometry of the cross section and its length. For plates, an analogous result holds if w is scaled by h/ 12(1 -ν 2 ). In this case, ε = 1 in Eq. ( 59) and the nonlinear dynamics of any flat plate of uniform thickness depends only on the geometry of the mid-plane (circular, rectangular…) and the boundary conditions. That is why, the numerical values of the nonlinear coefficients of Tables 9 and10 can serve as references.

3 Finite element reduced-order model

Finite element model and expansion on a linear modal basis

The general expression of a conservative nonlinear dynamical system with smooth (geometrical) nonlinearities, discretized with a finite element method, can be written in the following form:

M ẍ + f i (x) = f e , ( 12
)
where x, M, f i and f e denote, respectively, the vector of unknown displacements, the mass matrix, the internal and external force vectors. This model is obtained after any discretization of the displacement field in a finite element context. In this article, shell finite elements have been used, but the method is valid in the case of any displacement-based finite elements (beams, 3D…). It is convenient to split the internal force vector into a linear part and a purely nonlinear part. Assuming that the equilibrium point is at x = 0, the tangent stiffness matrix K writes: [START_REF] Chu | Influence of large amplitudes on free flexural vibrations of rectangular elastic plates[END_REF] so that the nonlinear internal force vector is defined by:

K = ∂ f i ∂ x x=0
f nl (x) = f i (x) -K x, ( 14
)
and the equations of motion are rewritten:

M ẍ + K x + f nl (x) = f e , ( 15
)
where the geometrically nonlinear part of the problem is concentrated in the nonlinear internal force vector f nl (x).

The displacements are expanded on a basis of N eigenmodes (ω r , Φ r ) of the linearized counterpart of (15):

x(t) = N r =1 Φ r q r (t). ( 16
)
where q r (t) is the r th unknown modal coordinate. The eigenmodes are solutions of:

(K -ω 2 r M)Φ r = 0. ( 17
)
By introducing Eq. ( 16) into Eq. ( 15), multiplying by Φ T r and using the orthogonality properties, one obtains the following set of differential equations verified by the modal coordinates q r (t):

qr + ω 2 r q r + N α=1 N β=α a r αβ q α q β + N α=1 N β=α N γ =β b r αβγ q α q β q γ = Q r , ( 18
)
where Q r = Φ T r f e /m r is the modal force and m r = Φ T r MΦ r is the r th modal mass, usually equal to 1 with a suitable normalization of the mode shapes Φ r . In set [START_REF] Sridhar | Non-linear resonances in the forced responses of plates, part II: asymmetric responses of circular plates[END_REF], we have written the nonlinear force vector expansion Φ T r f nl in the form of quadratic and cubic coupling terms. It is the result of an implicit Taylor expansion of f nl (x) around x = 0, truncated to the third order in q r , r = 1, . . . N . Depending on the formulation of the geometrical nonlinearities in the finite element discretization, this expansion can be exact (it is the case for 3D elements or for shell/plate/beam elements with a von Kármán strain/displacement nonlinear law [START_REF] Touzé | Direct finite element computation of non-linear modal coupling coefficients for reduced-order shell models[END_REF][START_REF] Lazarus | Finite elements reduced order models for nonlinear vibrations of piezoelectric layered beams with applications to NEMS[END_REF]) or truncated, for instance in the case of large rotation elements (such as geometrically exact beam elements for which the internal force vector include sine and cosine functions of the rotation degrees of freedom of the cross section [START_REF] Thomas | Hardening/softening behaviour and reduced order modelling of nonlinear vibrations of rotating cantilever beams[END_REF][START_REF] Cottanceau | A finite element/quaternion/asymptotic numerical method for the 3D simulation of flexible cables[END_REF]). Notice that Eq. ( 18) is written in an upper triangular form by organizing the coupling terms such that α ≤ β ≤ γ .

Computation of the nonlinear coefficients

In ROM [START_REF] Eisley | Nonlinear vibration of beams and rectangular plates[END_REF][START_REF] Sridhar | Non-linear resonances in the forced responses of plates, part II: asymmetric responses of circular plates[END_REF], parameters ω r , Φ r and Q r are obtained by the linear modal analysis of Eq. ( 17), available in any finite element code. The main issue is thus to compute the coefficients a r αβ and b r αβγ of the nonlinear terms. As stated in the Introduction section, several methods are available for this and we propose here to use the STiffness Evaluation Procedure (STEP) introduced in [START_REF] Muravyov | Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures[END_REF]. It consists in two steps.

-First, one has to prescribe a set of displacement vectors x which are linear combination of selected linear mode shapes Φ r and computing the corresponding nonlinear force vector f nl (x). This operation is available in any finite element code and is a simple evaluation of f nl (x) by knowing x. It thus does not require to solve any nonlinear system, for instance with a Newton-Raphson technique. -Second, it is shown that the unknown coefficients a r αβ and b r αβγ are solutions of linear systems for which the second members are known functions of the f nl (x) computed at the previous step.

To illustrate the method, we first show the computation of coefficients a r αα and b r ααα . We prescribe the following time-independent displacements to the structure:

x 1 = ±λΦ α ⇒ q α = λ q i = 0 ∀i = α ( 19
)
where λ refers to an amplitude coefficient of the mode Φ α whose value will be addressed hereafter. The second part of the above equations comes from Eq. ( 16), which states that since the modes are orthogonal, imposing x on the αth is equivalent to setting to zero all modal coordinates but the αth. Since x 1 is time independent, introducing Eq. ( 19) into Eqs. [START_REF] Woinowski-Krieger | The effect of axial force on the vibration of hinged bars[END_REF] and [START_REF] Sridhar | Non-linear resonances in the forced responses of plates, part II: asymmetric responses of circular plates[END_REF] leads to, for all r = 1, . . . N :

λ 2 a r αα + λ 3 b r ααα = Φ T r f nl (λΦ α )/m r , ( 20a
)
λ 2 a r αα -λ 3 b r ααα = Φ T r f nl (-λΦ α )/m r . ( 20b
)
As announced, the problem thus reduces to solve the above linear system which unknowns are the nonlinear coefficients a r αα and b r ααα . On a same way, the expansion on two or three modes (α, β, γ ) of the prescribed displacement

x 2 = ±λ α Φ α ± λ β Φ β (21a) x 3 = ±λ α Φ α ± λ β Φ β ± λ γ Φ γ , ( 21b
)
leads to other linear systems enabling to compute the a r αβ and b r αβγ with α = β = γ (see [START_REF] Muravyov | Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures[END_REF] for details). In the above developments, the values of the λ coefficients can be arbitrarily chosen, which means that for a given number of nonlinear coefficients to estimate, an infinite number of equations can be took into consideration. As Eqs. ( 19) and [START_REF] Capiez-Lernout | Computational stochastic statics of an uncertain curved structure with geometrical nonlinearity in three-dimensional elasticity[END_REF] show, the overdetermination comes from the amplitude as well as from the sign of the modal amplitude coefficients λ. Thus, the linear system whose unknowns are the nonlinear coefficients could be overdeterminated as well as minimal. In the first case, the overdetermined system is solved in a least square sense, as it has been done in [START_REF] Chang | Reduced order modelling for the nonlinear geometric response of some curved structures[END_REF][START_REF] Grolet | Dynamique non-linéaire des structures mécaniques: application aux systèmes à symétrie cyclique[END_REF]. The second case needs to dispose the minimal set of prescribed loads and nonlinear static computations. It is the approach chosen in this work.

For a linear mode basis of size N , the number N c of coefficients to compute is equal to:

N c = N N 2 + N 3 = N 4 6 + N 3 + 5N 2 6 , ( 22
)
where the number of combinations of size k with repetition in a given set of size n is the multiset number:

n k = n + k -1 k = (n + k -1)! k! (n -1)! ,
where the binomial coefficient is

n k = n! k! (n -k)!
and •! is the factorial of integer •. Equation [START_REF] Touzé | Direct finite element computation of non-linear modal coupling coefficients for reduced-order shell models[END_REF] shows that the computational burden increases with the power four of N . Figure 1 gives the evolution of N c as a function of N . One has to notice that the number of static computations necessary to obtain the N c coefficients is much less than N c : In Eq. ( 20), computing solely the two vectors f nl (±λΦ α ) leads to obtain the 2N coefficients a r αα and b r ααα with r = 1, . . . N .

Significant coefficients for flat structures

As stated in the previous section, the number of coefficients a r αβ and b r αβγ greatly increases with the number N of retained modes in the reduction basis. A question arises naturally: Does all the coefficients have the same impact on the nonlinear dynamics of the structure, or, to begin with, is there some of the coefficients that are simply zero ? To precise this idea, we compare the 22); " ," " ," " ," " ": N c from Eq. ( 24), with N = N u + N w and with, respectively,

N u = N w , N u = 2N w , N u = 5N w and N u = 10N w
present finite element formulations with the analytical models of flat structures described in Sect. 2. Since the structures under consideration are flat, the bending and in-plane motion are uncoupled at the linear stage and their modes Φ r [solutions of eigenproblem [START_REF] Ho | Non-planar, non-linear oscillations of a beam. part II: Free motions[END_REF]] are either pure bending modes Φ (w) r or pure in-plane modes Φ (u) r . It is then possible to separate them in the modal basis, which is now written

{Φ (w) 1 , . . . Φ (w) N w , Φ (u) 1 , . . . Φ (u)
N u } with N w and N u the number of bending and in-plane mode, respectively. In practice and in the case of plate/shell/beam finite elements, these modes can be calculated separately in most finite element codes, by separating the degrees of freedom involved in bending and in in-plane. This approach considerably reduces the computation time of in-plane modes which are located at higher frequency than bending modes.

Then, since the results obtained in Sect. 2 for analytical models of flat plates are valid for any flat structure with an arbitrary edge geometry, reduced-order model ( 7) can be identified to the one obtained with present finite element formulation [START_REF] Sridhar | Non-linear resonances in the forced responses of plates, part II: asymmetric responses of circular plates[END_REF], showing that a large number of coefficients a r αβ and b r αβγ are in fact zero for flat structures. More precisely, it depends on the values of α, β and γ that correspond to a bending mode Φ (w) r or a in-plane mode Φ (u) r . Table 1 gathers those results.

As a consequence, we rewrite the finite element ROM of Eqs. ( 18) by separating the modal coordinates into bending ones and in-plane ones and keeping only nonzero nonlinear coefficients. For this, indices i, j, k, l ∈ {1, . . . N w } will refer to bending coordinates (BC), whereas p ∈ {N w + 1, . . . N w + N u } will refer to in-plane coordinates (MC). One obtains:

BC: ∀k 1,...Nw qk + ω 2 k q k + Nw+Nu p=Nw +1 Nw i=1 a k pi q i q p + Nw i=1 Nw j=i Nw l= j b k i jl q i q j q l = Q k , ( 23a
) MC: ∀ p Nw+1,...Nw+Nu qp + ω 2 p q p + Nw i=1 Nw j=i a p i j q i q j = Q p , ( 23b
)
which is a set of oscillators of the same form than the one obtained for analytical models (Eqs. ( 7)).

The number of nonzero coefficients to be computed with this approach is:

N c = N w N w N u + N w 3 + N u N w 2 , = N 4 w 6 + N 3 w 2 + 3N 2 w N u 2 + N 2 w 3 + N u N w 2 . ( 24
)
In practice-it will be discussed in Sects. 4, 5 and 6one often includes a larger number of in-plane modes than transverse ones in the modal basis, which means N w N u . To illustrate this point, Fig. 1 shows N c as a function of the number of modes in the ROM expansion basis N = N u + N w , for several values of N u as a function of N w (N u = i N w with i ∈ {1, 2, 5, 10}). To illustrate the gain in the number of coefficients to be computed, N c as a function of N , from Eq. ( 22), is also plotted in Fig. 1. One can observe that the number N c of coefficients to be computed in the present approach can be largely reduced. Moreover, for a given modal basis size N = N u + N w , the number of coefficients N c to be computed largely decreases with the ratio i = N u /N w . In addition, as shown in "Appendix E", the required number of static computations N sc is one order of power smaller than N c :

N sc = N 3 w 6 + N 2 w 2 + N w 3 + N u N w . ( 25
)
In practice, N c can be reduced to a factor of about 1000 if N u = 10N w . For instance, if 4 transverse modes and 40 in-plane modes are retained (N w = 4, N u = 40) and the total number of modes is N = 44, the total number of coefficients is N c = 711,480, which becomes N c = 1120 by considering only the nonzero coefficients, with a requirement of only N sc = 180 static computations. or in-plane

Φ (u) r ) modes involved Type of coupled LM Type of oscillator q r (α, β, γ ) Bending r = k In-plane r = p Cubic coeff. b r αβγ 3 Bending LMs (i, j, l) b k i jl = 0 0 1 In-plane-2 bending LMs (m, i, j) 0 0 2 In-plane-1 bending LMs (m, p, i) 0 0 3 In-plane LMs (m, p, s) 0 0 Quadratic coeff. a r αβ 2 Bending LMs (i, j) 0 a p i j = 0 1 In-plane-1 bending LMs ( p, i) a k pi = 0 0 2 In-plane LMs (m, p) 0 0

Condensation of the in-plane coordinates

A further simplification of the finite element ROM can be obtained. In the same way than what was done for analytical models (Sect. 2), if the in-plane inertia is neglected and if no in-plane forcing is applied, the inplane coordinates are found to be quadratically linked to the transverse oscillators by Eq. ( 23b), for all p = 1, . . . N u :

q p = - 1 ω 2 p N w j=1 N w l= j a p jl q j q l . ( 26
)
The above in-plane coordinates can thus be condensed in transverse oscillators (23a), to obtain the following set of equations:

qk + ω 2 k q k - N w i=1 N w j=1 N w l= j ⎛ ⎝ N u +N w p=N w +1 a k pi a p jl ω 2 p ⎞ ⎠ q i q j q l + N w i=1 N w j=i N w l= j b k i jl q i q j q l = Q k ( 27
)
which leads to:

qk + ω 2 k q k + N w i=1 N w j=i N w l= j Γ k i jl q i q j q l = Q k , ( 28
)
where the coefficients of the cubic terms are defined in accordance with the upper triangular form chosen for the above equation, with l > j > i:

Γ k iii = b k iii - Nu +Nw p=Nw +1 a k pi a p ii ω 2 p , ∀k, i (29a) Γ k ii j = b k ii j - Nu +Nw p=Nw +1 a k pi a p i j +a k pj a p ii ω 2 p ∀k, i, j > i (29b) Γ k i j j = b k i j j - Nu +Nw p=Nw +1 a k pi a p j j +a k pj a p i j ω 2 p ∀k, i, j > i (29c) Γ k i jl = b k i jl - Nu +Nw p=Nw +1 a k pi a p jl +a k pj a p il +a k pl a p i j ω 2 p , ∀k, i, j > i, l > j (29d)
ROM [START_REF] Perez | Nonintrusive structural dynamic reduced-order modeling for large deformations: enhancements for complex structures[END_REF], valid if the in-plane inertia is neglected, is the most reduced since it involves only transverse coordinates, with all the in-plane motion embedded in the cubic coupling terms. It is also the most natural since it is very similar to classical linear modal models in which one considers only transverse modes. The issue of truncation of this model remains, but has been simplified since the choice of the in-plane modes to keep in the basis has been solved by the condensation.

Validation on beam and plate test cases

Definition of the two test cases

In order to test the finite element reduced-order modeling procedure of Sect. 3, two test cases have been selected, for which analytical solutions are at hand with the models of Sect. 2. Those two test cases are described in Table 2: a clamped-clamped straight beam of length The reduction basis consists of N u in-plane modes and N w bending modes, with a total of N = N w + N u linear modes (LMs) kept in the ROM of Eq. [START_REF] Sridhar | Non-linear resonances in the forced responses of plates, part II: asymmetric responses of circular plates[END_REF]. The first transverse and in-plane LMs of the two test cases, solutions of eigenproblem [START_REF] Ho | Non-planar, non-linear oscillations of a beam. part II: Free motions[END_REF], are shown in Tables 3 and4 . Since one of the purpose of the present study is to evaluate the effect on the ROM quality of the inplanes modes kept in the reduced basis, a more larger number of in-plane than transverse LMs is took into account, as shown in Table 5.

Geometry L × b × h = 1 × 0.05 × 0.001 m R × h = 0.3 × 0.001 m Material E = 210 GPa, ρ = 7800 kg/m 3 ν = 0 ν = 1/3 FE nodes 5 × 1 4 1 4 5 0 5 =

Numerical validation of significant coefficients

The purpose of this section is to numerically validate the results of Sect. 3.3 in which structures with a flat 1 In fact, the in-plane frequencies normalized values depend on parameter μ defined by Eqs. ( 2), [START_REF] Ciarlet | A justification of the von-Kármán equations[END_REF] but have a priori a negligible influence on the transverse motion. geometry are considered. The results presented in this paragraph were computed on the clamped beam, with a discretization consisting of 100 element in the length and 4 element in the width of a beam (see the mesh in Table 2), and maximum of modal amplitudes for the static nonlinear STEP computations fixed at

λ r max i |[Φ r ] i | = h/20, (30)
with r = {α, β, γ } (see Sect. 3.2) and [Φ r ] i the ith displacement2 component of vector Φ r In order to validate the numerical values of coefficients of nonlinear terms obtained with the STEP, they are compared to those obtained by the analytical model of Sect. 2.1 (the numerical values of some of them are given in Table 9). Since the coefficients value depend on the mode shapes normalization, the same normalization must be applied to both FE and analytical mode shapes. For this, the following dimensionless parameters are defined and applied to the FE model:

x = x x r , t = t t 0 Φr = υ r Φ r , ( 31
)
where x r , t 0 are defined in the same way than for analytical models (Eq. ( 2)): x r = h (resp. x r = h 2 /L) for the components of x related to transverse (resp. in-plane) motion; t 0 = L 2 √ ρ S/E I . The mode shape scaling υ r is chosen to obtain the same maximal amplitude for the analytical and FE mode shapes. Introducing those parameters into Eq. [START_REF] Eisley | Nonlinear vibration of beams and rectangular plates[END_REF][START_REF] Sridhar | Non-linear resonances in the forced responses of plates, part II: asymmetric responses of circular plates[END_REF] leads to replace them by

qr = q r υ r x r , ār αβ = t 2 0 x α x β x r υ α υ β υ r a r αβ , br αβγ = t 2 0 x α x β x γ x r υ α υ β υ γ υ r b r αβγ . ( 32
)
The two formulations are compared by analyzing the coupling coefficients of an expansion basis composed of three modes of the beam: the first transverse mode and the first two axial modes (in Table 3: mode 1 of the first column and modes 1,2 of the second column). Table 6 gathers the results. It shows that the estimation of the coefficients from the two formulations (FE with the STEP and analytic (u, w) -formulation, see "Appendix A" for details) is equivalent. In particular, the coefficients computed by the STEP that should be nonzero according to the theory of Sect. 3.3 are accurately obtained with at least 3 identical digits. For those that should be zero, values less than 10 -10 are obtained. One can notice that the analytical model shows that some coefficients G p i j and C k pi are zero: This is due to the chosen modes and the symmetry properties of the structure.

Sensitivity of the coefficients to the finite element discretization

In this section, we analyze the convergence of the numerical values of the coefficients computed by the STEP as a function of the FE mesh refinement. Figure 2 shows the convergence of (a) the eigenfrequencies ω k and (b) the cubic coefficient D k kkk ≡ b k kkk , for k = 1 . . . 5, computed with the FE code, as a function of the number of finite elements in the length of the beam. While the first eigenfrequencies are fastest converging toward the analytic ones, the convergence of the numerical estimations of cubic coefficients D k kkk is slower. A non-monotonic convergence of the coefficients is also noticed. It could be explained by the dependency of the nonlinear coefficients to the nodal force evolution along the axial direction, governed in the case of the beam by Eq. (44b): The zeros of the internal axial force would match the number of nodes of the mesh to be correctly approximated. Those tendencies are found identical to the one noticed in [START_REF] Lazarus | Finite elements reduced order models for nonlinear vibrations of piezoelectric layered beams with applications to NEMS[END_REF] in the case of a beam with a Euler-Bernoulli kinematics.

Sensitivity of the coefficients to the modal displacement amplitude

As explained in Sect. 3.2, the computation of the coefficients a k i j and b k i jl depends on free parameters λ r , and 19), [START_REF] Capiez-Lernout | Computational stochastic statics of an uncertain curved structure with geometrical nonlinearity in three-dimensional elasticity[END_REF], which have to be chosen with an amplitude able to activate the geometrical nonlinearities of the stiffness part of the model. In particular, if λ r is too small, the nonlinear part f nl (x) of Eq. ( 15) would be negligible with respect to the linear one K x and linear system (20) would be badly conditioned, the values of a k i j and b k i jl depending on divisions of 0 by 0.

Φ 3 = Φ (u) 2 Quadratic STEP a 1
To verify this point, Fig. 3 shows the evolution of some typical quadratic and cubic coefficients as a function of λ r . More precisely, measures of the amplitude of the prescribed displacement are defined as w max , u max = max(λ r Φ r ) with Φ r being, respectively, a bending mode 3 or an in-plane mode, scaled by the thickness h of the structure: w max / h, u max / h. According to Fig. 3 (left), the coefficients D i iii ≡ b i iii are merely constant on a large range of displacement. For very small transverse amplitudes (w max / h < 10 -4 ), the estimated coefficients present a random variation, probably due to the fact that nonlinearities are not activated enough. On the contrary, if displacements are too large (w max / h > 100), estimated values of nonlinear coefficients are underestimated, in a smooth tendency. This can be explained by considering that expansion [START_REF] Sridhar | Non-linear resonances in the forced responses of plates, part II: asymmetric responses of circular plates[END_REF], which includes quadratic and cubic terms only, is exact for moderate rotation models of the von Kármán type, as shown in Sect. 2. Those von Kármán models are known to be valid for moderate transverse amplitudes, of the order of the thickness h of the structure. On the contrary, the FE used for the STEP includes a large rotation formulation (the GROT_GDEP option is used in Code_Aster; see Sect. 4.1), valid for any displacement amplitude, that would activate geometrical nonlinearities of a higher order than 3, not included in the model of Eq. ( 18).

In the case of coefficients C i pi ≡ b i pi [Fig. 3(right)], which involves both axial and bending modes, there is also an accuracy plateau that cease to exist for u max / h 10. Above this critical amplitude, an erratic behavior is observed, which can probably be explained by observing that this critical amplitude is very large (ten times the thickness h of the beam), since the axial motion is physically one order of magnitude in h/L smaller than the transverse motion.

Those numerical experiments tend to show that the thickness h of the structure is the average order of magnitude of proper values of w max and u max . Consistently to these results, modal displacements amplitudes will thereafter be chosen as 1/20 of the thickness of the investigated structures (Eq. ( 30)).

Convergence of the cubic coefficients and in-plane motion condensation

As explained in Sect. 3.4, the in-plane/bending coupling inherent to the correct simulation of the geometrical nonlinearities in thin structures can be included in the cubic Γ k i jl coefficients, which depend on the quadratic coefficients a k pi and a p i j associated with the in-plane modes of the expansion basis of the ROM, as defined by Eqs. (29a-d). In this section, we analyze the convergence of the numerical values of some typical Γ k i jl with the number of in-plane modes retained in the summations of Eqs. (29a-d). The two test structures of Table 2 are considered here with the reference numerical values of the Γ k i jl coefficients, from the analytical models of Sect. 2, given in Table 10.

Fig. 5 Evolution of nonlinear plate cubic coefficients Γ i

iii with regard to number of in-plane modes retained in the summations of Eqs. (29a), for the six transverse modes of the 1st column of Table 4, i = 1, 2, 4, 6, 15, 30 Figure 4, in the case of the clamped-clamped beam, and Fig. 5, in the case of the clamped circular plate, highlight that a set of in-plane LMs are needed to converge to a nearly constant value of coefficients Γ i iii . It indicates also which of these oscillators participate in the nonlinear response of the structure. Thus, in the case of the beam, it is possible to associate a given transverse mode to a set of axial LMs, thus giving an indication of which axial modes are nonlinearly linked to a particular bending mode at large amplitude. Indeed, Fig. 4 shows that the first bending mode is mainly linked to the fourth axial mode, the second bending mode is notably associated with the second, fourth and sixth axial ones, whereas the third bending mode is linked to the second, sixth and eight axial modes. For higher orders, the convergence could be nearly reached by retaining only two axial modes in the condensation. Figure 4 also shows that the second axial mode has an effects on all cubic coefficient Γ i iii , although for some of them it could be considered as negligible. In the case of the plate, analogous results can be drawn by analyzing Fig. 5.

Since Figs. 4 and5 show coefficients Γ i iii computed with the STEP normalized by their analytical value (obtained by Eq. ( 10) and given in "Appendix D"), one can notice that the STEP has a tendency to underesti-mate the values of the coefficients, since the asymptotic value in Fig. 4 is slightly smaller than 1, the more and more as the order of the involved modes is increasing. It can be due to the mesh refinement, which is correct for low-order modes with a large wavelength and that brings errors for higher-order modes with a smaller wavelength.

Notice that this condensation has always the effect of decreasing the absolute value of the cubic coefficients by taking into account the contribution of inplane modes. In the case of the beam, most of the time, more than one axial mode has a significant effect on the transverse motion, except in the case of the first bending mode. This holds true as well for the circular plate-see Fig. 5-whose modal convergence needs more in-plane modes because of the strong density of asymmetric LMs.

Dynamical responses of structures

Solving of the nonlinear system

A continuation method is used to compute the dynamical response of reduced-order model [START_REF] Perez | Nonintrusive structural dynamic reduced-order modeling for large deformations: enhancements for complex structures[END_REF]. The con-tinuation method used in this text is the asymptotic numerical method (ANM) coupled to the harmonic balance method (HBM) [START_REF] Cochelin | A high order purely frequencybased harmonic balance formulation for continuation of periodic solutions[END_REF], implemented in the software Manlab [START_REF] Cochelin | Méthode asymptotique numérique[END_REF]. For this, the system is first rewritten in the state space. Secondly, a quadratic form is used, which considerably improve the performance of the continuation procedure. By introducing the state variable v k = qk and the product of two modal coordinates S i j = q i q j , j ≥ i, set ( 28) is replaced by:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ qk = v k vk = Q k -ω 2 k q k -2ξ k ω k v k -ε N w i=1 N w j=i N w l= j
Γ k i jl S i j q l 0 = S i jq i q j [START_REF] Kim | Nonlinear reduced-order modeling of isotropic and functionally graded plates[END_REF] with all the unknowns gathered in the following unknown vector:

u = q 1 q 2 . . . q N w v 1 v 2 . . . v N w S 11 S 12 . . . S N w N w (34)
A linear modal damping terms of factor ξ k has been added in the model. Periodic solutions of this system are computed. Following the HBM, u is expanded in Fourier series:

u(t) = u 0 + H h=1 u c h cos(hΩt) + u s h sin(hΩt) ( 35
)
where Ω is the angular frequency of the periodic motion. Introducing the above equation into [START_REF] Kim | Nonlinear reduced-order modeling of isotropic and functionally graded plates[END_REF], one obtains an algebraic quadratic system:

R(U, Ω) = 0, (36)
where

U = u 0 u c 1 u s 1 . . . u c H u s H , ( 37
)
is the vector of unknowns that gathers all the Fourier components of the unknowns.

In the software Manlab [START_REF] Cochelin | Méthode asymptotique numérique[END_REF] used for the computations, system (33) is given as input. Then, the application of the HBM to write algebraic system (36) is automatic and final algebraic system [START_REF] Kim | Nonlinear reduced order modeling of flat cantilevered structures[END_REF] is solved by the ANM, which computes U and Ω as Taylor series of a pseudo-arclength parameter a [START_REF] Cochelin | A high order purely frequencybased harmonic balance formulation for continuation of periodic solutions[END_REF]:

U = U 0 + aU 1 + • • • + a n U n , Ω = Ω 0 + aΩ 1 + • • • + a n Ω n ( 38
)
where n is the order of the series, chosen by the user (n = 20 in the present study). Then, (U 0 , Ω 0 ) are solutions of a nonlinear algebraic system solved by a Newton-Raphson algorithm, whereas U i , i = 1, . . . , n are solutions of a cascade of linear systems.

In this text, two kinds of computations are done:

-the first one is the free conservative response of the investigated system, i.e., the backbone curve of a given nonlinear mode (NNM), which means the system of Eqs. [START_REF] Kim | Nonlinear reduced-order modeling of isotropic and functionally graded plates[END_REF] with Q k = 0 and ξ k = 0 ∀k ∈ {1, ..., N w }. In this case, as proposed in [START_REF] Karkar | A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities[END_REF], a perturbation term -λv k is added to the second line of Eq. ( 33) at each oscillator k, where λ is an explicit continuation parameter, which enables to add a phase condition. According to [START_REF] Munoz-Almaraz | Continuation of periodic orbits in conservative and Hamiltonian systems[END_REF], it ensures that the computed solutions are those of the unperturbed and conservative system of Eqs. ( 28) and [START_REF] Kim | Nonlinear reduced-order modeling of isotropic and functionally graded plates[END_REF]. The phase conditions used in this work consist in setting to zero the sine component of the first harmonic of q i in (35) if the ith NNM is computed. Moreover, the initial Newton-Raphson algorithm (to compute U 0 in Eq. ( 38)) needs to be initialized in the vicinity of a given NNM. In this work, for the ith NNM, the unknowns are initialized on the corresponding linear mode: q i (t) = q 0 cos ω i t, v i (t) = -ω i q 0 sin ω i t, S ii = 1 2 (1+cos 2ω i t), Ω = ω i and all the other unknowns are set to zero. The amplitude q 0 is chosen small with respect to the plate thickness h.

-The second one is the harmonically forced response of the system, considering consistent values of Q k and ξ k .

The size of algebraic system (36) increases fast with the number of harmonics and transverse oscillators retained in the basis. In practice, there is generally no simplification allowing to neglect the coupling between two oscillators, which infers that all additional variables S i j must be included in u. Then, the size N s of the final unknown vector U of a system with N w oscillators is given by:

N s = (2H + 1) 2N w + N w (N w + 1) 2 . ( 39
)

Modal bases of the ROMs

Since condensed model ( 9) is used for the computations, the expansion basis has to be selected in two steps:

-firstly, one has to select a given set of N w transverse modes to be retained in ROM [START_REF] Noor | Recent advances in reduction methods for nonlinear problems[END_REF]; -secondly, one has to select a larger set of N u N w modes used to compute the quadratic coefficients a k pi and a p ii necessary for the estimations of the cubic coefficients Γ k i jl (see Eqs. ( 28)) with the condensation of the in-plane motion into the bending motion.

In the results presented in the next paragraphs and figures, we assume that all cubic coefficients Γ k i jl have attained convergence with a suitable number N u of in-plane modes. Those computations are made offline before the computation of the dynamical response with Manlab.

For the beam examples, the basis composed of the first 10 bending modes has been used. For the circular plate simulations, one has to consider that asymmetric modes appear in pair in the basis (so-called companion modes). Three different bases are used, in regard to possible modal interactions between axisymmetric and companion asymmetric modes. The first modal basis is presented in Fig. 6a and is composed of the first 5 axisymmetric modes; the second one includes the two first pairs of asymmetric modes (4 modes) also shown in Fig. 6a; the third one is composed of the first 12 LMs of the circular plate, but includes only one of the companion modes associated with each asymmetric mode shape, as shown in Fig. 6b.

Nonlinear normal modes and forced responses

A useful property of NNM is that all nonlinear forced responses of the system are distributed around the computed free undamped response. This is illustrated in Fig. 7 which presents the forced responses of the clamped beam and the clamped circular plate around their first NNM.

More precisely, a point harmonic excitation at location ζ of amplitude F and frequency Ω has been selected, for which the value of the modal forcing in Eqs. ( 28), [START_REF] Kim | Nonlinear reduced-order modeling of isotropic and functionally graded plates[END_REF] is

Q k = Φ k (ζ )P cos(Ωt). ( 40
)
where Φ k (ζ ) is the transverse modal amplitude at location ζ . The forced response of the beam is computed with a damping value of each LMs fixed at ξ k = 2% and the excitation point chosen at x = L/π , such that all LMs are involved in the displacement response. For the plate, the damping values are ξ k = 5% and the excitation point is located at r = 0.5461. A tongue is noticeable on the free response of the two structures. It corresponds to 1:5 and 1:7 internal resonances, respectively, for the beam and the plate. This phenomena are not visible on the forced response, and it would appear if the dampings of the interacting LMs were reduced.

Figure 7 enables to draw generic results. As announced, the free response (also known as the NNM or the backbone curve) represents the skeleton of the dynamics since it connects all the forced responses at their resonance point. Since it is obtained in a single continuation, it appears to be an efficient way of characterizing the dynamics of the structure. In the following, the free responses of the beam and the plate are considered to analyze the behavior and the efficiency of the numerical method.

Nonlinear modes of the plane structures

In this section, we present an overview of the dynamical response of the clamped beam and the circular plate, through their first four NNMs, shown in Figs. 8 and9 as frequency energy plots (FEP). In fact, NNMs are the free conservative oscillations of the structure (with Q k = 0 and ξ k = 0 ∀k ∈ {1, ..., N w } in Eqs. ( 33)), for which the mechanical energy H is an invariant at any given amplitude. Here, the energy can be written:

H = 1 2 N w k=1 q2 k + ω 2 k q 2 k + V nl ( 41
)
where V nl denotes the nonlinear part of the energy defined in "Appendix F". In practice, the backbone curves shown in Figs. 8 and9 are computed with the ANM/HBM method (see Sect. 5.1) and H is computed in post-processing in the time domain, with the Fourier coefficients of q k and qk , by taking the mean value over one period of Eq. ( 41). H = 10 harmonic components are used for the computations. The 10 first LMs compose the reduction basis for the computation of the beam FEP (Fig. 8). For the plate FEP (Fig. 9), the first basis (Fig. 6a) is used for the axisymmetric NNMs (0,1) and (0,2), whereas the second one is used for the asymmetric ones (1,1) and (2,1). All the NNMs show the same qualitative behavior. First, since both structures are plane, all the NNMs have a hardening behavior, which globally bends the backbone curves toward the high frequencies (indeed, the Γ i iii coefficients are all positive; see "Appendix D"). Secondly, some internal resonances modify the topology of the backbone curve, which is described more precisely in the next section.

A focus on internal resonances

On the FEP of Figs. 8 and9 , only simple internal resonances of the form 1:n with n ∈ N are observed. They are characterized by a strong coupling between two NNMs (the ith and the jth) when the following frequency relation is fulfilled:

ω jnl nω inl , ( 42
)
where ω inl , ω jnl are the oscillation frequencies of the two involved NNMs, which are read on the x-axis of the FEP. This kind of internal resonance appears here in two distinct ways, which can be gathered in two families.

The first family gathers 1:n internal resonance that naturally appears at high energy levels because of the bending of the backbone curves. At low energy, ω inl and ω jnl are close to the natural frequencies ω i and ω j , which are far from fulfilling ω j = nω i . Then, because of the bending of the backbones, the oscillations frequencies ω inl , ω jnl change and Eq. ( 42) is fulfilled at high energy. This kind of internal resonances always shows the same features, illustrated by Figs. 10 and11: -if an IR tongue emerges from the main backbone curve in the FEP, it is characterized in terms of frequency content by a transfer of energy between the first harmonics q i1 (H1) of mode i and the nth harmonics q jn (Hn) of mode j. -there is always a particular point on the IR tongue, close to its tip, for which Eq. ( 42) is an equality (ω jnl = nω inl ). At this point, denoted B in Figs. 10 and 11 , the tongue of the ith NNM is connected to the backbone curve of the jth NNM. Precisely, it is the curve H = f (ω inl ) that is connected to the curve H = f (ω jnl /n), shown in gray dashed line in Figs. 8, 9, 10 and 11 . -when following the tongue from its emergence from the ith backbone (point A) to point B, a transfer of energy is observed from the first harmonics of mode i to the nth harmonics of mode j. Indeed, the latter increases while the former decreases, giving a continuous change of the dynamics between mode i oscillating at Ω and mode j oscillating at nΩ. At point B, mode i is zero and the dynamics only involves mode j oscillating at nΩ. -since the observed internal resonances tongues are with n odd, they naturally emerge from the main backbone without any bifurcation point.

The second family gathers the 1:1 internal resonances observed in the backbone curves of the asymmetric modes of the circular plate (modes (1,1) and (2,1) of Figs. 9 and 12), for which the natural frequencies of the companion modes are naturally tuned and fulfill Eq. [START_REF] Sombroek | Numerical computation of nonlinear normal modes in a modal derivative subspace[END_REF]. In this case, well documented in [START_REF] Touzé | Asymmetric non-linear forced vibrations of free-edge circular plates, part I: theory[END_REF][START_REF] Manevitch | Free oscillations in conservative and dissipative symmetric cubic two-degree-of-freedom systems with closed natural frequencies[END_REF], the backbone curve of one of the companion modes becomes unstable after a pitchfork bifurcation (denoted by PF in Fig. 12), appearing at very low energy. It gives rise to a new branch, for which both companion modes have nonzero energy and oscillate with a quadrature phase shift, giving rise to a traveling wave.

On the FEPs of Figs. 8 and9 , the internal resonances are labeled as "1:n IR," with the mode in interaction specified by "NNM j Ω/n." For instance, the first NNM of the beam shows a 1:5 internal resonance with the third NNM.

Convergence with the number of harmonic components

Figure 13 presents the free response of the beam around its first NNM with a varying number of harmonic components retained in the Fourier series expansion of Eq. ( 35), with a converged (see next section) number of transverse modes [N w = 10 in Eqs. ( 28), [START_REF] Kim | Nonlinear reduced-order modeling of isotropic and functionally graded plates[END_REF] Frequency dim] 0 0.5

Amplitude w/h

H = 1 H = 2 H = 3 H = 10
Fig. [START_REF] Chu | Influence of large amplitudes on free flexural vibrations of rectangular elastic plates[END_REF] Convergence with the number of harmonic components of the free responses of the beam around its first NNM. Maximum amplitude of the transverse displacement of the periodic response computed at x = L/π with N w = 10 transverse LMs. The x-axis is scaled by the first eigenfrequency ω 1 of the beam sufficient to predict correctly the backbone curve anywhere out of the internal resonance tongue. One would have expect that H = 1 harmonics would have been sufficient with the HBM, whereas H = 3 is necessary here. This can be explained by the quadratic recast imposed by the HBM (see Eq. ( 33)), which imposes the introduction of the quadratic variables S i j which need at least two harmonics to accurately embed the nonlinear response. When more harmonics are retained, the free response changes only with the apparition of the internal resonance. In this case, since it is a 1:5 internal resonance, at least H = 5 harmonics are necessary to correctly compute it.

Modal convergence of NNMs

We investigate now the convergence of the computation of the free response of the first NNMs of the beam and the plate as a function of the number of modes retained in the ROM basis. For the first NNM of the clamped beam, Fig. 14a shows that a correct computation is obtained with only the first oscillator-the first LM of the clamped beam-for a range of transverse displacements w until half the thickness h. The contribution of the second LM is negligible, whereas the effect of the third LM is significant on a larger range of displacement. The convergence of the backbone curve is thus achieved with the first and third linear modes only. The internal resonances need to retain the LMs interacting with the first oscillator: Indeed, the second internal resonance involves here a LM of higher order than the third and then cannot be computed with three oscillators. The second NNM shows analogous results. The third NNM needs to take more LMs into account to limit the error in comparison with the reference solution, as shown in Fig. 14c. Adding an oscillator to the system has either a softening or an hardening effect on the response: It is associated with the sign of the nonlinear coefficients Γ which involve the retained LMs. In particular, the LMs of lower order than the computed NNM must be retained in the basis: The softening effect of the first oscillator on the third NNM cannot be neglected.

In the case of the plate, the convergence is shown in Fig. 15 for the first two axisymmetric NNMs (mode (0,1) and mode (0,2)) and the first asymmetric one (mode (1,1)), by comparing several modal basis and H = 3 and H = 10 harmonic components. As for the beam case, the reduction to only one oscillator is correct on a moderate range of displacements, the responses needing to include more LMs in the basis after a certain value of the displacements. One can notice that an example of modal interaction between axisymmetric and asymmetric modes is exhibited in the case of NNM (0,2) (Fig. 15b): The internal resonance which is noticeable here needs asymmetric modes in the basis to be obtained. In this section, we test the proposed method on a plane structure with an original geometry, for which we want to compute the nonlinear dynamics. It has the form of a perforated rectangular plate with three rectangular holes that delineate two beams at the center, as shown in Fig. 16. The main parameters of the model are gathered in Table 7. The dimensions of this structure have been adjusted in order to obtain two linear modes (LM) which present mode shapes similar to those of the first LM of a clamped-clamped beam, respectively, with an in-phase and an out-of-phase motion between the two beams. This is verified by analyzing the modal basis of the structure, shown in Table 8, for which the 7th and 8th modes have the required shape. They are similar to the two first modes of a structure composed of two clamped-clamped beams connected by a linear spring of small stiffness. This particular structure has been selected to analyze and compute the nonlinear dynamics of those two beam-like modes and their bifurcation to a localized motion on one of the beams [START_REF] King | Mode localization in a system of coupled flexible beams with geometric nonlinearities[END_REF][START_REF] Grolet | Computing multiple periodic solutions of nonlinear vibration problems using the harmonic balance method and Groebner bases[END_REF].

Computation of the nonlinear coefficients of the ROM

Since the studied structure is a plate with a particular boundary, all the theory described in the present article can be applied. We then compute its nonlinear dynamics with the reduced-order model of Eq. ( 28) in which the Γ k i jl includes the in-plane motion effect, 1) with the STEP.

Frequency
Figure 17 shows the convergence of the condensation of coefficients Γ k i jl as a function of the number N u 123 Fig. 17 Evolution of the physical nonlinear coefficients Γ i iii of the perforated plate in regard to the number of in-plane modes of in-plane modes retained in Eq. ( 29). This number has been chosen voluntarily high-1000 in-planes LMsin order to discuss the convergence of cubic coefficients throughout the condensation. Although numerous in-plane modes are included in the modal basis, the evolution of the nonlinear cubic condensed coeffi-cients presented in Fig. 17 seems to need more in-plane modes to converge toward a stable value. This could be explained by the high density of the platelike LMs. As for the beam and circular plate, the value of the coefficients decreases with the number of condensed in-plane LMs, which then have a softening effect on the non- 43), for the clamped-clamped beam, as a function of k. (left) ROM with 4 bending modes i, j, l, k ∈ {1, . . . 4}; (left) ROM with 16 bending modes i, j, l, k ∈ {1, . . . 16} linear system. Notice that the ratio between the initial and final value of the cubic coefficients is greater than those of the beam and the plate. The low-order in-plane modes have a strong influence on the cubic coefficients related to platelike transverse modes. For the coefficients Γ 7 777 and Γ 8 888 , there is a large range of in-plane LMs which does not influence the condensation.

Sorting and selection of the ROM coefficients

After those operations, the computation time of the nonlinear dynamics by the ANM/HBM is mainly sensitive to the number of nonlinear terms associated with coefficients Γ k i jl in Eqs. [START_REF] Kim | Nonlinear reduced-order modeling of isotropic and functionally graded plates[END_REF]. Among this set of cubic coefficients, many of them are still negligible. Consequently, an heuristic filtering procedure is applied. First, coefficients Γ k i jl are scaled with Eqs. ( 31), [START_REF] Ehrhardt | Finite element model calibration of a nonlinear perforated plate[END_REF], with x r = h, t 0 = 1/ω 1 and v r = 1. Then, observing Eq. ( 32) naturally shows that the nonlinear coefficients are scaled by t 2 0 , equivalent to the inverse of the square of a frequency and because Γ k i jl depends on modes i, j, l, k, we propose to remove from system [START_REF] Perez | Nonintrusive structural dynamic reduced-order modeling for large deformations: enhancements for complex structures[END_REF] all coefficients such that:

Γ k i jl ωi ω j ωl ωk < , with ωi = ω i ω 1 , Γ k i jl = h 2 ω 2 1 Γ k i jl . ( 43
)
with the filtering criteria chosen as = 10 -3 .

To verify this scaling, Fig. 18 shows the application of this filtering process on the cubic coefficients of the clamped-clamped beam of Sects. 2 and 4 (see Table 10). First of all, one can observe that the scaling of Eq. ( 43) almost counterbalances the increasing the cubic coefficient with their order k: Coefficients Γ k kkk are of the order of magnitude of 1 independently of k. Then, the other scaled coefficients decrease with the natural frequency distance between the involved modes i, j, l and k. With = 10 -3 , a 1000 factor is kept between the higher and the smaller cubic coefficient. In this case, all beam coefficients for i, j, l, k ∈ {1, . . . 4} are kept, whereas the smallest of them are neglected for higher ranges of variation of i, j, l, k.

A two LMs reduced-order model

At first, a reduced-order model composed of the 7th and 8th oscillators only, associated with mode 7 and 8 of interest, has been selected. The two backbone curves of the corresponding NNMs are shown in Fig. 19 for an increasing number of harmonics. This plot shows the classical picture of a nonlinear localization of the motion. NNM7, of lower frequency, has a standard hardening backbone curve. On the contrary, NNM8 shows a pitchfork bifurcation that give birth to two branches on which the motion is a coupling between mode 7 and 8, which are associated with a localization of the motion on beam 1 or beam 2 (depending on the considered branch). A minimal number of three 123 Fig. [START_REF] Ducceschi | Nonlinear dynamics of rectangular plates: investigation of modal interaction in free and forced vibrations[END_REF] Modal convergence with the number of transverse LMs retained in the basis, for NNMs 7 and 8. Displacements at the center of one of the beams harmonic components in the harmonic balance formulation is necessary to reach convergence.

Convergence in terms of LMs

We now consider an enrichment of the 2 degrees of freedom basis of the previous section by including the first 16 transverse modes in the ROM basis. The Fig. 20 shows that this enrichment does not significantly change the geometry of the backbone curves but renders the computations tedious since more degrees of freedom are involved. It thus validates, in this case, the two LM model of the previous section.

Conclusion

In this paper, a complete procedure to compute the geometrically nonlinear frequency response of flat structures by using reduced finite element models has been proposed.

A study at the coefficient level has confirmed that a large part of the ROM coefficients are equal to zero, which enables to significantly reduce the computational burden of the stiffness evaluation procedure, in terms of number of coefficients to compute (and stock) and number of nonlinear static problems. The resulting reduced-order model can be further simplified by condensing the in-plane motion to the transverse one, such that only transverse modes coupled by cubic nonlinearities appear in the system. A convergence of the cubic coefficients is reached when a sufficient number of in-plane modes is retained. At this stage of the procedure, all coefficients of the ROMs of a straight beam and circular plate have been validated by a direct comparison with analytical models. It is shown that the ROM coefficient value is independent of the prescribed displacement amplitude on a large range, attesting the robustness of the method and providing some guidelines for the STEP. Moreover, an heuristic rule, to filter the significant ROM coefficients, is also provided. It is validated for flat structures and could be applied with no restriction to more complex geometries, extending the relevance of the present study.

An original coupling of the obtained nonlinear systems with a ANM-HBM following-path method allows to equally compute the forced as well as the free response of the plane structures. In accordance with previous studies, the distribution of the forced response around the free response and the backbone curve was attested with this method.

Finally, the computation of the free responses-also known as nonlinear modes-highlights many typical nonlinear phenomena as hardening behaviors and internal resonances. These latter can complexify the resolution procedure when too many oscillators and harmonics are retained in the reduced-order model, because of the computation of all modal interactions which occur between the retained oscillators. At each resonance, a very reduced model, where only one transverse oscillator is kept and with H = 3 harmonics, has been found sufficient to compute the backbone curves with a negligible error on an acceptable range of motion amplitude. In the case of a 1:1 internal resonance, in axisymmetric structures or linked to vibration localization, a model with only the two transverse modes involved in the internal resonance seems sufficient in a large range of amplitude.

= 0 C 1 11 = 0 C 1 21 = 9.663 C 1 31 = 0 C 1 41 = -102.7 C 1 51 = 0 ω 2 = 61.67 D 2 222 = 1773.7 G 1 22 = 0 G 2 22 = 55.02 G 3 22 = 0 G 4 22 = -130.8 G 5 22 = 0 C 2 12 = 0 C 2 22 = -110.0 C 2 32 = 0 C 2 42 = 261.6 C 2 52 = 0 ω 3 = 120.9 D 3 333 = 7970.2 G 1 33 = 0 G 2 33 = 91.97 G 3 33 = 0 G 4 33 = 55.24 G 5 33 = 0 C 3 13 = 0 C 3 23 = -183.9 C 3 33 = 0 C 3 43 = -110.5 C 3 53 = 0
axial mode shapes are solutions of Eq. (6b) and can be written:

Ψ p (x) = √ 2 sin γ p x, γ p = pπ. ( 49
)
They naturally verify Eq. ( 47). The values of some of these coefficients are given in Table 9.

Appendix B: On the plate (u, w)-formulation ROM

Some details about the analytical plate model of Sect. 2.2 are given here. The in-plane forces in the plate are represented by a two-dimensional tensor N, which, with the von Kármán assumptions, can be written [START_REF] Thomas | Geometrically non-linear flexural vibrations of plates: in-plane boundary conditions and some symmetry properties[END_REF]:

N = A [(1 -ν) + ν tr 1] , = 1 2 ∇u + ∇ T u + ∇w ⊗ ∇w (50)
where is the in-plane strain tensor, A = Eh/(1ν 2 ) is the in-plane stiffness of the plate, ∇ denotes the vector/tensor gradients of scalar/vector fields, tr is the trace of a tensor ⊗ is the tensor product of two vectors and 1 is the identity tensor. By eliminating N between Eqs. ( 50) and (11a, b), one obtains the following (u, w)-formulation:

ρh ẅ + D w-A (1-ν) ∇u : ∇∇w + u • ∇w + ν div u w + ∇div u • ∇w + ∇w • ∇∇w • ∇w + 1 2 (∇w) 2 w = 0. (51a) ρh ü - A 2 (1 -ν) u + (1 + ν)∇div u + (1 -ν) w∇w + (1 + ν)∇∇w • ∇w = 0. ( 51b
)
where : denotes the doubly contracted product of two tensors, the vector Laplacian of a vector field and • the dot product.

With the following dimensionless variables:

w = w h , r = r L F = F Eh 3 , ε = 12(1 -ν 2 ), t = 1 L 2 D ρh t, p = L 4 Dh p, ( 52
) ε = Ah 2 D = ε 1 -ν 2 = 12, μ = D AL 2 , ū = L h 2 u, N = L 2 Eh 3 N. ( 53
)
where L is a characteristic dimension of the middle plane of the plate (its diameter for instance). The (u, w)-formulation is then rewritten: (

ẅ + w -ε (1 -ν) ∇u : ∇∇w + u • ∇w + ν div u w + ∇div u • ∇w + ∇w • ∇∇w • ∇w + 1 2 (
) 55
The modes are chosen to satisfy the following two eigenproblems associated with the linear parts of Eqs. (54a, b):

Φ k -ω 2 k Φ k = 0 (56a) (1 -ν) Ψ p + (1 + ν)∇div Ψ p -γ 2 p Ψ p = 0 (56b)
where ω k and γ p are the transverse and in-plane dimensionless eigenfrequencies (the frequency equation of the in-plane vibrations of a clamped circular plate has been derived in [START_REF] Park | Frequency equation for the in-plane vibration of a clamped circular plate[END_REF]). Then, by multiplying (54a, b), respectively, by Φ k and Ψ p , integrating on the midplane surface S and using the orthogonality properties, the ROM of Eqs. (7a, b) is obtained, with:

C k pi = -S (1 -ν) ∇Ψ p : ∇∇Φ i + Ψ p • ∇Φ i + ν div Ψ p Φ i + ∇div Ψ p • ∇Φ i Φ k dS S Φ 2 k dS , ( 57a
)
G p i j = S (1 -ν) Φ i ∇Φ j + (1 + ν)∇∇Φ i • ∇Φ j • Ψ p dS S Ψ 2 p dS , ( 57b
)
D k i jl = -S ∇Φ i • ∇∇Φ j • ∇Φ l + 1 2 ∇Φ i • ∇Φ j Φ l Φ k dS S Φ 2 k dS . ( 57c
)
The above developments are not analytically investigated in this study: They have been reported here for a sake of completeness and to precisely demonstrate the results of Sect. 2.3.

Appendix C: On the plate (w, F)-formulation ROM

In the literature dealing with analytical models of plates, the (u, w)-formulation of the previous section is never used since it would lead to tedious-maybe impossible-computations in an analytical context (see Eq. ( 51)) and also because there are three unknown displacement fields (w and the two components of u). In most of the cases, the in-plane inertia is neglected and Eq. (11b) reduces to div N = 0, which enables the introduction of an Airy stress function F(r, t) defined by:

F1 + ∇∇F = N. ( 58
)
Then, by adding a compatibility condition, one obtains the following (w, F)-formulation (all details can be found in [START_REF] Thomas | Geometrically non-linear flexural vibrations of plates: in-plane boundary conditions and some symmetry properties[END_REF]):

ẅ + w = εφ(w, F) + p(r, t) (59a)

F = - 1 2 φ(w, w) (59b)
where the overbars have been dropped and where φ(•, •) is a bilinear operator [START_REF] Thomas | Geometrically non-linear flexural vibrations of plates: in-plane boundary conditions and some symmetry properties[END_REF]. Similarly to the beam case, the transverse displacements and the force function are expanded on a linear basis of vibration modes Φ k and Υ j :

w(r, t) = N w k=1
Φ k (r)q k (t) and

F(r, t) = N F j=1
Υ j (r)η j (t) [START_REF] Grolet | Dynamique non-linéaire des structures mécaniques: application aux systèmes à symétrie cyclique[END_REF] where the q k and η j are time-dependent modal coordinates, respectively, associated with bending and inplane eigenmodes, which verify:

Φ k -ω 2 k Φ k = 0, (61a)
Υ j + ζ 4 j Υ j = 0. ( 61b
)
After multiplying equation (59a) by a mode Φ l , equation (59b) by a mode Υ m and integrating on the surface of the plate, the orthogonality property of the eigenmodes yields:

qk + ω 2 k q k -ε N w p=1 N F j=1 E k r j q r η j = Q k (t) (62a) ζ 4 j η j = - 1 2
N w p=1 N w q=1 H j pq q p q q , (62b)

where nonlinear quadratic coefficients (E k r j , H j pq ) and modal load Q k are defined in "Appendix C". Finally, eliminating η j (t) between Eqs (62a, b) leads to the same set of N w coupled transverse oscillators with cubic nonlinearities than in the case of beams, Eq. ( 9), with the following values of the condensed cubic coefficients:

Γ k pqr = N F j=1 E k r j H j pq 2ζ 4 j . ( 63
)
Again, the above-defined nonlinear coefficients depend on the mode normalization, defined by Eq. ( 65).

The analytical expressions of the coefficients of nonlinear terms in Eq. ( 63) are:

E k r j = S Φ k L(Φ r , Υ j ) dS Φ 2 dS , H j pq = S Υ j L(Φ p , Φ q ) dS S Υ 2 j dS , Q k (t) = S Φ k p(r, t) dS S Φ 2 k dS . ( 64
)
The modes are here normalized as:

S Φ 2 k dS = S Υ 2 p dS = 1. ( 65
)
Step 1 Imposing a displacement on a single BM, x 1 = λΦ α , and expanding the result on either a BM Φ k or a MM Ψ p leads to:

λ 3 D k ααα = Φ T k f nl (λΦ α ), ( 66
)
λ 2 G p αα = Ψ T p f nl (λΦ α ). ( 67
)
Consequently, only N w static computations are necessary to obtain coefficients D k ααα and G p αα .

Step 2 Imposing displacements on two BM, x 2 = λ α Φ α ± λ β Φ β , with β > α, and expanding the result on either a BM Φ k or a MM Ψ p leads to: Step 4 Imposing a displacement on a three BM, x 4 = λ α Φ α + λ β Φ β + λ γ Φ γ , and expanding the result on a BM Φ k leads to:

λ 2 α λ β D k ααβ + λ α λ 2 β D k αββ =Φ T k f nl (λ α Φ α + λ β Φ β ) -λ 3 α D k ααα -λ 3 β D k βββ , (68)
λ α λ β λ γ D k αβγ = Φ T k f nl (λ α Φ α + λ β Φ β + λ γ Φ γ ) - i, j,l∈{α,β,γ } l≥ j≥i i jl =αβγ λ i λ j λ l D k i jl . ( 72
)
This computation of the D k αβγ coefficients thus requires N 3 w /6 -N 2 w /2 + N w /3 static computations.

Appendix F: Computation of the energy for the FEPs

The nonlinear part of the potential energy in Eq. ( 41) writes:

V nl = N w i=1 Γ i iii q 4 i 4 + N w i=1 N w j>i Γ i
ii j q 3 i q j 3 + Γ j ii j q 2 i q 2 j + Γ j i j j q 3 j q i 3

+ N w i=1 N j>i N w k> j Γ i i jk
2 q 2 i q j q k + Γ j i jk 2 q i q 2 j q k + Γ k i jk

2 q i q j q 2 k + N w i=1 N w j>i N w k> j N w l=1
Γ l i jk q i q j q k q l . (73)

To obtain this equation, the starting point is the general expression of a dynamical system with cubic nonlinearities:

qk + ω 2 k q k + N w i=1 N w j=i N w l= j
Γ k i jl q i q j q l = Q k , (74) which can also be written q + f (q) = 0, where f denotes the internal forces vector. It derives from a potential and can be integrated over time if and only if the following condition on the crossed derivatives is verified:

∂ f i ∂q j = ∂ f j ∂q i ∀ i = j. ( 75
)
According to the developments of [START_REF] Thomas | Geometrically non-linear flexural vibrations of plates: in-plane boundary conditions and some symmetry properties[END_REF], it is established for some classical boundary conditions of platesincluding the clamped one-that the symmetry properties of the bilinear operator L lead to the following equality between the nonlinear coefficients of the (w, F)-formulation ROM presented in "Appendix B":

H j pq = E q pj , H j pq = H j qp . ( 76
)
It results after some developments to some equalities between the non-upper triangular form coefficients Γ k i jl , i, j, l, k ∈ [1, N w ], which are written:

Γ k i jl = Γ k jil = Γ i k jl = Γ i kl j = Γ j lki = Γ j kli = Γ l jki = Γ l k ji , (77)
thus, the upper triangular form yields

Γ i ii j = 3Γ j
iii , Γ i i j j = Γ j ii j , Γ j i j j = 3Γ i j j j ∀ j > i Γ j i jl = 2Γ i j jl = 2Γ l i j j , Γ i i jl = 2Γ l ii j = 2Γ j iil , Γ l i jl = 2Γ

j ill = 2Γ i jll ∀ l > j > i Γ k i jl = Γ l i jk = Γ j ilk = Γ i jlk ∀ k > l > j > i, ( 78
)
which allows to verify the condition of Eq. (75). The general expression of the energy in Eq. ( 73) is then obtained by multiplying the equation of the kth oscillator of Eq. (74) by the modal velocity qk and integrating over time.

Fig. 1

 1 Fig.[START_REF] Chen | Direct observation of coherent energy transfer in nonlinear micromechanical oscillators[END_REF] Order of coefficients of the ROM as a function of the number N of modes in the expansion basis. '•': N c from Eq. (22); " ," " ," " ," " ": N c from Eq. (24), with N = N u + N w and with, respectively, N u = N w , N u = 2N w , N u = 5N w and N u = 10N w

Fig. 2 Fig. 3

 23 Fig. 2 Convergence of the first eigenfrequencies ω k (left) and of the cubic coefficients D k kkk ≡ b k kkk (right) as a function of the number of FE in the length of the beam, computed by the FE

Fig. 4

 4 Fig. 4 Evolution of nonlinear beam cubic coefficients Γ iiii with regard to the number of axial modes retained in the summations of Eqs. (29a), for the first six transverse modes of the beam (1st column of Table3, i = 1, . . .[START_REF] Monteil | Identification of mode couplings in nonlinear vibrations of the steelpan[END_REF]

Fig. 6

 6 Fig. 6 Organization of the two bases of the ROMs for the circular plate computations: eigenfrequencies of the modes as a function of there number. a First basis and b second basis

60 Fig. 7 Fig. 8

 6078 Fig. 7 Maximum amplitude of the transverse displacement of the periodic response of the beam (left) and circular plate (right) around their first NNM. The x-axis is scaled by the first eigenfrequency ω 1 of the structure. Free response in black and forced

Fig. 10 Fig. 11

 1011 Fig. 10 Topology of the 1:5 internal resonance of NNM1 with NNM3 of the clamped-clamped beam. a Zoom of the FEP of Fig. 8 and b harmonics amplitude of NNM1 (H1 of LM1 in blue

Fig. 14

 14 Fig. 14 Modal convergence with the number of transverse modes of the free responses of the beam around the a first, b second and c third NNMs of the beam. Maximum amplitude of

Fig. 15 Fig. 16

 1516 Fig.[START_REF] Woinowski-Krieger | The effect of axial force on the vibration of hinged bars[END_REF] Modal convergence with the number of transverse modes of the free responses of the circular plate around some NNMs: a first axisymmetric (0,1) (NNM 1), b second axisymmetric (0,2) (NNM 6) and c first asymmetric (1,0) (NNM 2).

Fig. 18

 18 Fig. 18 Scaled cubic coefficients Γ k i jl of Eq. (43), for the clamped-clamped beam, as a function of k. (left) ROM with 4 bending modes i, j, l, k ∈ {1, . . . 4}; (left) ROM with 16 bending modes i, j, l, k ∈ {1, . . . 16}

Fig. 19

 19 Fig.[START_REF] Touzé | Asymmetric non-linear forced vibrations of free-edge circular plates, part I: theory[END_REF] Response of the two LMs reduced-order model. a Convergence as a function of the number of harmonics components of the modal coordinates q 7 and q 8 on the backbones of NNMs

  ∇w) 2 w = 0. (54a)μ ü -(1 -ν) u + (1 + ν)∇div u + (1 -ν) w∇w + (1 + ν)∇∇w • ∇w = 0. (54b)The bending and in-plane displacements are expanded onto a linear basis of, respectively, N w transverse modes Φ k and N u axial modes Ψ p :w(r, t) = N w k=1 Φ k (r)q k (t) and u(r, t) = N u +N w p=N w +1 Ψ p (x)η p (t).

-λ 2

 2 α λ β D k ααβ + λ α λ 2 β D k αββ =Φ T k f nl (λ α Φ α -λ β Φ β ) p f nl (λ α Φ α + λ β Φ β ) -λ 2 α G p αα + λ 2 β G p ββ (70)where the last part of the second members is known from the first step. At this step, since β > α, N w (N w -1) static computations are required.Step 3 Imposing displacements on one BM and one MM, x 3 = λ α Φ α + λ β Ψ β , and expanding the result on a BM Φ k leads to:λ α λ β C k αβ = Φ T k f nl (λ α Φ α + λ β Φ β ) -λ 3 α D k ααα (71)Since α = 1, . . . N w and β = 1, . . . N u , the computation of the coefficients C k αβ thus requires N w N u static computations.

Table 1

 1 Nonzero coefficients a r αβ and b r αβγ of nonlinear terms for flat structures as a function of the type (bending Φ

	(w) r

Table 2

 2

	Mesh, geometrical
	and mechanical
	characteristics of the finite
	element beam and plate test
	cases
	The FE dof. row gives the
	number of degrees of
	freedom of the finite
	elements model

Table 3

 3 Order of appearance in the basis (#), eigenfrequencies and shapes of the first transverse and axial modes for the clamped-clamped beam. The axial mode shapes are shown with arrows showing the in-plane vector field corresponding to Φ

	(u) p ,

Table 4 Label

 4 The labels of the transverse modes are a doublet (k, n) with k the number of nodal diameters and n the number of nodal circles

			Bending modes			In-plane modes
		#	ω k /(2π) ω/ω 1	Shape	ω k /(2π)	ω/ω 1	Shape
			[Hz]			[Hz]
	(0,1)	1	28.69	1		5673.9	197.74
	(1,1)	2	59.70	2.08		6463.2	222.25
	(2,1)	4	97.90	3.41		8808.7	306.99
	(0,2)	6	111.68	3.89		9077.1	316.43
	(0,3) 15	143.15	4.99		11,189.3	389.96
	(0,4) 30	170.71	5.95		11,413.1	397.76

, number of appearance in the basis (#), eigenfrequencies and shapes of the first transverse and in-plane modes for the clamped circular plate. The in-plane mode shapes are shown with arrows showing the in-plane vector field corresponding to Φ (u) p , with colors related to their length.

Table 5

 5 Modal basis of the two test cases for which the coefficients of the ROMs are computed

	Test case	Type of LM	Number of LMs
			in the ROM
	Beam	Transverse	10
		Axial	31
	Plate	Axisymmetric transverse	5
		Asymmetric transverse	2 × 23
		Total transverse	51
		In-plane	750
	identical ones that now depend on scaled modal coor-
	dinate qr and coefficients ār αβ and br αβγ , defined by:

Table 6

 6 Comparison of nonlinear stiffness coefficient estimations with one bending mode Φ 1 = Φ

	(w) 1 and two axial modes Φ 2 = Φ	(u) 1

Table 3 ,

 3 i = 1, . . . 6)

		1.45											1.6												1.6
		1.35 1.4		P	l	t a	e	d o m	e	, 0 (	1	)	1.5		P	l	t a	e	m	o	e d	1 (	0 ,	)	1.5		P	l	t a	e	d o m	e	2 (	0 ,	)
	Γ 111 1,FE / Γ 111 1,AN	1.1 1.3 1.15 1.2 1.25										Γ 222 2,FE / Γ 222 2,AN	1.1 1.2 1.3 1.4											Γ 444 4,FE / Γ 444 4,AN	1.2 1.3 1.4
		1 1.05											1												1.1
		0.95											0.9												1
			0	200	400				600				800	0	200	400				600				800	0	200	400				600	800
				In-plane mode number					In-plane mode number					In-plane mode number
			1.6											1.6												1.5
			1.5		Plate mode (0,2)	1.5		Plate mode (0,3)	1.4		Plate mode (0,4)
	Γ 666 6,FE / Γ 666 6,AN	1.1 1.2 1.4 1.3										Γ 151515 15,FE / Γ 151515 15,AN	1.1 1.2 1.3 1.4											Γ 303030 30,FE / Γ 303030 30,AN	1.1 1.2 1.3
			1											1												1
			0.9											0.9												0.9
			0	200	400				600				800	0	200	400				600				800	0	200	400				600	800
				In-plane mode number					In-plane mode number					In-plane mode number

Table 7

 7 Parameters of the model

	Parameters	Values
	Thickness	1 mm
	Young modulus	210 GPa
	Poisson's ratio	0.38
	Number of transverse LMs	16
	Number of in-plane LMs	1000
	which has been condensed thanks to Eqs. (29a-d). Only
	the nonzero coefficients a k i j and b k i jl are computed (see
	Table

Table 8

 8 Natural frequencies and shapes of the first twelve transverseLMs of the free perforated plate

				Order of appearance #	1	2	3	4
				Natural frequency [Hz]	184.8	189.1	405.3	480.1
				ω/ω 1	1.00	1.023	2.193	2.598
				Shape
				Order of appearance #	5	6	7	8
				Natural frequency [Hz]	508.3	561.7	736.7	867.3
				ω/ω 1	2.750	3.039	3.987	4.693
				Shape
	The 7th and 8th LMs, investigated in this study,	Order of appearance # Natural frequency [Hz]	9 918.6	10 1017.0	11 1156.0	12 1309.1
	are shown in bold		ω/ω 1	4.970	5.503	6.255	7.083
				Shape
		8
		6
	Γ 111 1	4
		2
		0
		0	500	1000
		In-plane mode number

Table 9

 9 Numerical values, computed with a trapezoidal rule, of some nonlinear coefficients of a the clamped-clamped beam model (Eqs. (46a, b)). Scaling of Eqs. (2) and mode normalization of Eqs.[START_REF] Kerschen | Nonlinear normal modes, part I: a useful framework for the structural dynamicist[END_REF]

		γ 1 = 3.141	γ 2 = 6.283	γ 3 = 9.425	γ 4 = 12.566	γ 5 = 15.708
	ω 1 = 22.37	D 1 111 = 111.2
		G 1 11 = 0	G 2 11 = -4.832	G 3 11 = 0	G 4 11 = 51.35	G 5 11

We consider here only the transverse displacement degrees of freedom and not the eventual rotation ones, encountered in Timoshenko (shear deformation) 1D finite elements.

See footnote 2.

Acknowledgements

The French Ministry of Research is warmly thanked for the financial support of this study, through the Ph.D. Grant of the first author.

Compliance with ethical standards

Conflict of interest

The authors declare that they have no conflict of interest

Appendix A: The beam (u,w) model

The beam model of Eq. ( 1) is here discussed. Using Hamilton principle, one obtained the classical equations of motion for the axial and transverse motion of the beam, which are:

where p(x, t) and n(x, t) are distributed transverse and in-plane forces per unit length and where the axial force in the beam is:

Equations ( 45) and (44a, b) constitute the classical von Kármán model describing the transverse motion of a beam with geometrical nonlinearities [START_REF] Nayfeh | Nonlinear Oscillations[END_REF][START_REF] Woinowski-Krieger | The effect of axial force on the vibration of hinged bars[END_REF][START_REF] Eisley | Nonlinear vibration of beams and rectangular plates[END_REF]. Then, eliminating N between Eqs. ( 45) and (44a, b) leads to the (u, w) model of Eq. ( 1).

After the modal expansion of Sect. 2.1, the analytical expressions of the coefficients of nonlinear terms in Eq. (7a, b) are:

The modes are here normalized as:

In the case of a clamped-clamped boundary conditions, the transverse mode shapes are solution of Eq. (6a) and can be written:

k for all k > 1 and κ k is numerically computed to verify Eq. [START_REF] Kerschen | Nonlinear normal modes, part I: a useful framework for the structural dynamicist[END_REF]. The

Appendix D: Numerical values of beam and plate cubic ROM coefficients

The numerical values of the Γ k i jl coefficients of Eq. ( 9), in the case of the clamped-clamped beam and the circular plate, are gathered in Table 10. They have been computed by the analytical models (Eq. ( 10) and ( 63)) and are identical to those computed by the STEP (Sect. 3.2).

Table 10 Numerical values of the cubic coefficients Γ k

i jl of Eq. ( 9), for the clamped-clamped beam [(u, w)-formulation, Eq. ( 10), scaling of Eq. ( 2) and mode normalization of Eq. ( 47)] and the clamped circular plate [(w, F)-formulation of section C, Eq. ( 63), scaling of Eqs. [START_REF] Kuether | A numerical approach to directly compute nonlinear normal modes of geometrically nonlinear finite element models[END_REF] and mode normalization of Eqs. ( 65

Appendix E: Details of STEP

We develop here the details to compute only the nonzero coefficients of ( 23) with the method of [START_REF] Muravyov | Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures[END_REF]. We refer to the notations C k pi , D k i jl and G p i j of Eq. ( 7). Four distinct steps are necessary, by separating bending modes (BM) and in-plane modes (MM). We also consider unitary modal masses (m r = 1 ∀r ).

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.