L Meyrand
email: meyrand@lma.cnrs-mrs.fr

E Sarrouy

B Cochelin
email: bruno.cochelin@centrale-marseille.fr

G Ricciardi

Nonlinear Normal Mode continuation through a Proper Generalized Decomposition approach with modal enrichment

Keywords: Continuation, Model reduction, Proper Generalized Decomposition (PGD), Nonlinear dynamics, Nonlinear Normal Mode, Harmonic Balance Method (HBM)

Nonlinear normal modes (NNMs) of a mechanical structure provide a mean to understand nonlinear vibrational phenomena measured experimentally. As they are an interesting extension of the linear normal modes (LNMs), numerical methods are needed to estimate these oscillatory motions, generally including continuation aspects with respect to the circular frequency. In this paper, a continuation method combined to a model reduction based on Proper Generalized Decomposition (PGD) technique is described. This formulation mixes PGD, harmonic balance method (HBM) and basic continuation techniques in order to reach a highly reduced description of the NNMs. This PGD/HBM-based continuation algorithm includes a modal enrichment as the NNM energy grows.

The method is applied to two conservative problems including a cubic spring and a unilateral contact.

Introduction

In the frame of nonlinear dynamics, the Nonlinear Normal Modes (NNMs) have been developed for several decades, in the aim to use them as their wellknown linear counterparts (LNMs for Linear Normal Modes). Their properties are investigated in order to enable faster computations using modal reduction and to ease the forced dynamics prediction and interpretation [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF][START_REF] Laxalde | Dynamics of a linear oscillator connected to a small strongly non-linear hysteretic absorber[END_REF]. Basically, NNMs are the sets of periodic solutions of an unforced and -generally -undamped mechanical system [START_REF] Rosenberg | On nonlinear vibrations of systems with many degrees of freedom[END_REF][START_REF] Shaw | Normal modes for non-linear vibratory systems[END_REF][START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF]. The last decade witnessed the development of numerical methods dedicated to NNMs but this research area has not yet reached maturity [START_REF] Peeters | Nonlinear normal modes, part ii: Toward a practical computation using numerical continuation techniques[END_REF][START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF]. Although algorithms can now be effectively applied to real engineering structures, the discrete systems to solve are often high-dimensional problems and their nonlinear contribution causes convergence issues: NNMs require large sets of variables to be properly described which make them hard to compute in short times and hard to reuse even if some works still take advantage of their computation [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steadystate vibrations in non-conservative systems[END_REF].

Model reduction techniques have taken an important place in the solving process of mechanical systems, as their numerical complexity is increasing and the models are becoming more and more detailed and realistic. The continuation of NNMs using reduced-order models has already been addressed, for instance in [START_REF] Kuether | A numerical approach to directly compute nonlinear normal modes of geometrically nonlinear finite element models[END_REF][START_REF] Sombroek | Numerical computation of nonlinear normal modes in a modal derivative subspace[END_REF]. Both these works use a shooting method combined to a pseudo-arclength continuation. Model reduction comes from the space where the q 0 point from which time integration starts is looked for. Kuether and Allen [START_REF] Kuether | A numerical approach to directly compute nonlinear normal modes of geometrically nonlinear finite element models[END_REF] build this reduced space by using either a set of LNMs directly or the static displacement resulting from a force based on a set of LNMs while Sombroek et al. [START_REF] Sombroek | Numerical computation of nonlinear normal modes in a modal derivative subspace[END_REF] use some LNMs and their modal derivatives. This paper proposes to compute NNMs using a Proper Generalized Decomposition (PGD) approach combined to a continuation method. The PGD methods are a class of reduction techniques based on the separation of the unknowns which allows dealing with smaller systems to reduce the algorithmic complexity and the computation times. A review of the PGD framework has been detailed by Chinesta in [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] and the literature shows that a wide range of problems can be managed with this approach [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | Pgd-based computational vademecum for efficient design, optimization and control[END_REF][START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF]. PGD was already used by Grolet and Thouverez [START_REF] Grolet | On the use of the proper generalised decomposition for solving nonlinear vibration problems[END_REF] to compute free and forced response of nonlinear systems. The interest of the paper is first to couple this approach to a continuation method with a progressive enrichment of the PGD modes family and second to improve some of the PGD steps to get a robust algorithm. Contrarily to methods proposed in [START_REF] Kuether | A numerical approach to directly compute nonlinear normal modes of geometrically nonlinear finite element models[END_REF][START_REF] Sombroek | Numerical computation of nonlinear normal modes in a modal derivative subspace[END_REF], the reduced space is not built using a priori quantities but on the fly. The resulting method is applied to two simple systems which allows for an accurate analysis of its performance. The detailed study of the physical meaning of NNMs and their internal resonances or their use for the computation of forced responses is beyond the scope of this article.

The reader is referred to [START_REF] Hill | Identifying the significance of nonlinear normal modes[END_REF][START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF][START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steadystate vibrations in non-conservative systems[END_REF] for detailed works on theses aspects.

The paper is organized as follows: the mechanical framework is described in Sec. 2, including a short review of the NNM definition. In Sec. 3 the PGD and the Harmonic Balance Method (HBM) frameworks are exposed. The differences with Grolet and Thouverez [START_REF] Grolet | On the use of the proper generalised decomposition for solving nonlinear vibration problems[END_REF] work are highlighted. Then the PGD/HBMbased continuation is described in Sec. [START_REF] Shaw | Normal modes for non-linear vibratory systems[END_REF]. This part is a proposition of a full algorithm to compute NNMs with on the fly enrichment of the reduced basis in conjunction with continuation techniques. The last section, Sec. 5 is dedicated to two illustrative cases.

Mechanical framework

This section first describes the general equations of a nonlinear dynamical problem and then defines the associated NNMs.

Nonlinear dynamical system

After a spatial discretisation, using for example the Finite Elements method, the differential equations governing the motion of a nonlinear dynamical system can usually take the following form [START_REF] Renson | Nonlinear Normal Modes of Nonconservative Systems[END_REF][START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steadystate vibrations in non-conservative systems[END_REF], where the nonlinear efforts are separated from the linear ones:

Mü + C u + Ku + f nl (u, u) = f (t) (1)
Eq. ( 1) is a set of N second order nonlinear differential equations, the vector of the unknowns u -or vector of the degrees of freedom (dof) -is a vector of R N . M, C and K are respectively the linear mass, damping and stiffness matrices, with {M, C, K} ∈ M N (R) 3 . f (t) is the vector of excitation forces and lies in R N too. Finally f nl (u, u) is the vector of nonlinear forces and belongs to R N .

The general system described by Eq. ( 1) represents a very wide class of mechanical problems. Indeed, the f nl term can be any nonlinear expression (e.g. cosine function, polynomial terms like u i u 3 j , non smooth function, etc.).

Short review of nonlinear normal modes

This paper focuses on the computation of nonlinear normal modes. NNMs can be defined as the periodic solutions of Eq. ( 1) in the particular case when no forcing neither damping is considered or as an extension of the definition of linear normal modes which takes into account the nonlinear vector f nl . Multiple definitions of NNMs can be found in the literature: NNMs were introduced by Rosenberg in the 1960's [START_REF] Rosenberg | On nonlinear vibrations of systems with many degrees of freedom[END_REF] for conservative systems as a periodic vibration in which all the degrees of freedom are synchronous (same origin, same extrema, same period) and can be expressed from only one of them. Shaw and Pierre [START_REF] Shaw | Non-linear normal modes and invariant manifolds[END_REF][START_REF] Shaw | Normal modes for non-linear vibratory systems[END_REF] extended this definition in the 1990's to all the autonomous systems, including the dissipative ones. A NNM is then defined as a two-dimensional invariant manifold in phase space. The reader can find more information about NNMs definitions in the works of Kerschen and Renson [19,[START_REF] Renson | Nonlinear Normal Modes of Nonconservative Systems[END_REF][START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF] for instance.

More recently, Haller and Ponsioen [START_REF] Haller | Nonlinear normal modes and spectral submanifolds: existence, uniqueness and use in model reduction[END_REF] introduced a new definition for NNMs and the notion of spectral submanifolds.

Assuming the Rosenberg framework -that is undamped modes -a NNM is a set of limit cycles -or periodic solutions -of the following system:

Mü + Ku + f nl (u) = 0 (2)
To shorten further expressions, a signal u(t) satisfying Eq. ( 2) will be equivalently denoted:

R(u(t)) = 0 (3)
Although they do not have orthogonality properties like LNMs and do not decouple the equations, NNMs can be useful to investigate the modal interaction between widely spaced modes or modal bifurcations [START_REF] Renson | Nonlinear Normal Modes of Nonconservative Systems[END_REF][START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF]. One of the major features of NNMs is that frequencies depends on the response amplitude on the contrary to LNMs.

According to the size and the complexity of the dynamical system studied, computing branches of NNMs can be expensive in terms of data storage and computation times [START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF]. This observation motivates the developments described in the rest of the paper.

PGD and HBM combination to describe periodic solutions

A HBM framework is interesting as NNMs are periodic solutions. Despite a truncation order dependence, it allows dealing with algebraic equations in the frequency domain instead of (nonlinear) ordinary differential equations in the time domain. However a simple HBM becomes costly as the number of dof increases and more harmonics are needed to accurately describe the behaviour of a system. A combination with a PGD reduction technique allows to compute smaller systems without losing HBM advantages.

In this section the HBM framework used to compute periodic solutions is first given. Then PGD principles are introduced: PGD temporal and spatial subproblems are defined in combination with the HBM. Finally the full PGD/HBM algorithm is described.

Harmonic Balance Method

Unlike shooting methods [START_REF] Sundararajan | Dynamics of forced nonlinear systems using shooting/arc-length continuation method-application to rotor systems[END_REF][START_REF] Peeters | Nonlinear normal modes, part ii: Toward a practical computation using numerical continuation techniques[END_REF], the HBM provides a mean to compute limit cycles without the need of a time integration scheme. Only main steps are briefly described here, introducing some notations which are later used. The reader is referred to the numerous works on this subject for detailed explanations [START_REF] Urabe | Galerkin's procedure for nonlinear periodic systems[END_REF][START_REF] Urabe | Numerical computation of nonlinear forced oscillations by galerkin's procedure[END_REF][START_REF] Nayfeh | Nonlinear Oscillations[END_REF][START_REF] Sarrouy | Non-linear periodic and quasi-periodic vibrations in mechanical systems -on the use of the harmonic balance methods[END_REF].

Basic principle

HBM is based on a Galerkin method: looking for a periodic signal u(t) with period T = 2π/ω, a truncated Fourier series is used as a parametrization:

u(t) = a 0 √ 2 + H k=1 (a k cos(kωt) + b k sin(kωt)) (4)
where a k and b k are respectively the cosine and sine coefficients of the Fourier series. This decomposition also provides expressions for any u derivative with respect to time. Injecting this decomposition in the equations governing u and its derivatives (i.e. Eq. (1) or Eq. ( 2)) and following the Galerkin procedure, the equations are projected onto each element of the Fourier basis [1/ √ 2, cos(ωt), sin(ωt), . . . , cos(Hωt), sin(Hωt)] using the inner product:

< f, g > T = 2 T T 0 f (t)g(t)dt (5)
This leads to an algebraic set of N (2H + 1) equations where N is the size of the time problem; the unknowns are a k and b k vectors and sometimes ω when there is no external forcing (e.g. Eq. ( 2)). These equations are denoted:

H(a 0 , a 1 , b 1 , . . . , ω) = 0 (6)
When the system has an external periodic forcing, that is when it is non autonomous, ω in the Fourier series is taken equal to the forcing circular frequency ω * ; the resulting system, with ω = ω * , is then square and has, locally, one solution only. In the case of NNMs computation, the system is autonomous and for a given ω = ω * value, although the system is still square, there is, locally, not only one solution but a family of solutions. An additional equation usually called phase condition [START_REF] Seydel | From Equilibrium to Chaos, Practical Bifurcation and Stability Analysis[END_REF] is then used to choose one solution among this family. Possible phase conditions and their consequences on the problem size are mentioned later in Sec. 3.1.3.

AFT method to get nonlinear forces Fourier decomposition

The expression of the projection of the linear part of Eq. ( 1) or (2) can be established analytically [START_REF] Sarrouy | Non-linear periodic and quasi-periodic vibrations in mechanical systems -on the use of the harmonic balance methods[END_REF]. When considering Eq. ( 2), it is H l (ω) u H , where:

         H l (ω) = diag(Λ 0 , Λ 1 , . . . , Λ H ) Λ 0 = K, Λ k = diag(K -(kω) 2 M, K -(kω) 2 M) u H = {a T 0 , . . . , a T k , b T k , . . . } T (7)
The difficult and time consuming step when solving this system is the eval- The algebraic nonlinear system to solve writes:

H l (ω) u H + H nl (u H , ω) = 0 (8)

Phase condition when circular frequency is unknown

When considering an autonomous system such as Eq. (2), no "initial time" exists and using u(t), t ∈ [0, T ] or u(t + t), t ∈ R to describe the limit cycle are equivalent. Hence, the Fourier coefficients are not uniquely defined and it is preferable to add an equation, frequently called phase condition so as to define a unique solution line of Eq. [START_REF] Kuether | A numerical approach to directly compute nonlinear normal modes of geometrically nonlinear finite element models[END_REF]. Different equations can be added: the amplitude or velocity of a given dof can be imposed when t = 0 in Eq. ( 4) [START_REF] Seydel | From Equilibrium to Chaos, Practical Bifurcation and Stability Analysis[END_REF]. This effectively adds an equation linking Fourier coefficients of this dof together. This may imply that solutions violating this condition will not be found. Another solution is to nullify one of the Fourier coefficients components;

this is a priori of a lower cost as it removes one unknown [START_REF] Sarrouy | Stochastic study of a non-linear selfexcited system with friction[END_REF]. Moreover, considering a conservative system with no symmetry, choosing to nullify one sine coefficient implies that all other sine coefficients are null, hence, problem size is divided by almost 2, having only to solve for (H + 1) a k coefficients only based on the N (H + 1) equations provided by projection on cosine functions only. Other equations are of course possible such as the CHBM developped by Coudeyras et al. [START_REF] Coudeyras | A new treatment for predicting the self-excited vibrations of nonlinear systems with frictional interfaces: The constrained harmonic balance method, with application to disc brake squeal[END_REF] which involves the eigenvalues analysis of the tangent system or the integral condition used in AUTO [START_REF] Doedel | Auto: A program for the automatic bifurcation analysis of autonomous systems[END_REF].

When adding an equation, the resulting system (6) for an arbitrarily chosen ω * value is not square anymore: the N (2H +1) Fourier coefficients are described by the N (2H + 1) + 1 equations (HBM equations plus phase condition). There are different ways to circumvent this problem as stated in [START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF].

The rest of the paper focuses on undamped NNMs which can be described by a collection of cosine coefficients only. Systems with out-of-unison behaviour as defined by Hill et al. [START_REF] Hill | Out-of-unison resonance in weakly nonlinear coupled oscillators[END_REF] are excluded. This is equivalent to a phase condition imposing a null velocity at initial time for one dof for systems without symmetries. Using this assumption, all the equations derived from the projection onto sine functions sin(kωt) are useless and then removed. Equation ( 8) still holds but with the following definitions replacing definitions given in Eq. ( 7):

         H l (ω) = diag(Λ 0 , Λ 1 , . . . , Λ H ) Λ 0 = K, Λ k = K -(kω) 2 M ∈ M N (R) u H = {a T 0 , . . . , a T H } T (9)
For a given ω * value, the system is then square, having N (H + 1) unknowns (a 0 , a 1 , . . . , a H ) driven by N (H + 1) equations (projections of Eq. ( 2) onto [1/ √ 2, cos(ωt), . . . , cos(Hωt)] functions).

PGD combined to HBM

According to Chinesta et al. [START_REF] Chinesta | Pgd-based computational vademecum for efficient design, optimization and control[END_REF] and Grolet and Thouverez [START_REF] Grolet | On the use of the proper generalised decomposition for solving nonlinear vibration problems[END_REF], a PGD process can simply be reduced to three fundamental steps: separating variables, obtaining as many subproblems as there are variables, and run an alternated directions algorithm which solves each subproblem with the other variables which are fixed, in a same loop iteration. In this section the temporal and the spa-tial problems are obtained and coupled with HBM features to produce a global PGD/HBM algorithm.

Separable variables assumption

As a first PGD requirement, we choose to look for a specific form of solutions including a space-time separation. u(t), solution of Eq. ( 2), is sought as:

u(t) ≈ m j=1 p j q j (t) ⇔ u(t) ≈ Pq(t) with P = [p 1 , . . . , p m ] ( 10
)
where m is a positive integer such as m << N and denotes the number of PGD modes used for u decomposition, P ∈ M N,m (R) is the matrix of the m PGD mode shapes p j , and q(t) ∈ R m is the vector containing the time dependence of each PGD mode. As a product of two data is involved, the spatial modes p j are normalized to 1 arbitrarily; amplitude information is then reflected by q j (t) only.

One of the main advantages of the PGD approach is to split the problem into two smaller problems: one in space returning a value of P and one in time to get q. As we are here looking for periodic solutions, q will be decomposed using HBM. Following Sec. 3.1, q(t) is approximated by its truncated Fourier series:

q(t) = a 0 √ 2 + H k=1 a k cos(kωt) (11)
A compact way to write q(t) is then:

q(t) = Q H h H (t, ω), where    Q H = [a 0 , a 1 , . . . ] h H = [1/ √ 2, cos(ωt), . . . ] T (12)
Computing q(t) is then equivalent to computing both Q H and ω. For examples of problems solved by PGD in the time domain, the reader can refer to Nouy [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF] for instance.

Temporal problem T m

In this part, the mode shapes P are assumed to be known. The objective is to compute q(t) -i.e. Q H and ω -with P as a parameter. To this end, a weak formulation is used: the test function u (t) is written u (t) = P q (t) and q(t) must satisfy Eq. ( 13) for an oscillation period

I t = [0, 2π/ω]. ∀t ∈ I t ∀q (t) It q T (t)P T R(Pq(t))dt = 0 (13)
A few matricial calculus are needed to get the temporal problem ( 14) from the weak formulation ( 13):

∀t ∈ I t M r q(t) + K r q(t) + f nlr (Pq(t)) = 0 (14)
where

{M r = P T MP, K r = P T KP} ∈ M m (R) 2 and f nlr = P T f nl ∈ R m .
This subproblem is a system of m second order ODEs with a nonlinear contribution solved using HBM as described in Sec. 3.1. A new algebraic system with m (H + 1) equations and m (H + 1) + 1 unknowns -Q H and ω -is then obtained. Indeed ω is unknown for a NNM computation and a constraint equation c(q H , ω) is needed to make the system square. Using the notations introduced in Eq. ( 8) and ( 9), the temporal problem T m can be defined as follows:

T m (Q H , ω|P) = 0 ⇔    H l (ω)q H + H nl (q H , ω) = 0 c(q H , ω) = 0 ( 15
)
where q H = {a T 0 , a T 1 , . . . } T contains the columns of Q H , H l is the linear con-185 tribution, H nl is the vector of Fourier coefficients of the nonlinear contribution of Eq. ( 14), and c(q H , ω) ∈ R is an arbitrary constraint equation which allows to get a square system. c(q H , ω) will be specified in Sec. 4.1 as a continuation criterion.

Spatial problem S m 190

In this part, q(t) -that is Q H and ω -are supposed to be known and only

P is modified. The test function is then u = m k=1
p k q k which leads to a set of m systems with N equations:

∀k ∈ [[1; m]], m j=1 It q k qj dt M + It q k q j dt K p j + It q k f nl (Pq(t)) dt = 0 (16)
This can be written in a more compact way as a unique N × m sized system defining S m :

S m (P|Q H , ω) = 0 ⇔ S l p + S nl (p) = 0 (17)
where

• p = {p T 1 , . . . , p T m } T ∈ R N ×m contains the columns of P;
• S l can be defined analytically:

S l = I 2 ⊗ M + I 0 ⊗ K (18)
with

I 0 = π w Q H Q T H , I 2 = -πω Q H D 2 Q T H and D = diag(0, 1, . . . , H); • S nl (p) = [ It q k f nl (Pq(t)
) dt] 1≤k≤m is computed using the AFT approach described in Sec. 3.1.

The size of this system depends on N so this is more expensive to solve than the temporal problem as m << N , especially in an industrial context.

PGD/HBM core global algorithm

Both temporal and spatial problem being properly defined, they are embedded in the global PGD/HBM algorithm which also addresses the convergence criterion issue and additional modes initialization. It is described in Algorithm 1.

First, an error measure ε is defined (lines 1 and 11). It compares the norm of the residue vector R(u(t)) where u(t) = PQ H h H (t, ω) to the norm of the linear stiffness contribution over a period I t . This lets us introduce a convergence criterion for the alternated directions algorithm loop (line 4) while Grolet

and Thouverez [START_REF] Grolet | On the use of the proper generalised decomposition for solving nonlinear vibration problems[END_REF] imposed to loop over k three times without checking the convergence. The loop is then broken as soon as the number of iterations k becomes larger than a given value k max or the error measure ε becomes smaller than a given limit ε max or when the current iteration did not provide enough error reduction, that is, when ε > r ad ε prev , with r ad < 1 an arbitrary positive factor, and ε prev being the error measured at the end of the previous iteration. Data: Initial values for m, P, Q H and ω.

1 ε ← max t∈It (||R(u(t))||/||Ku(t)||); 2 while m ≤ m max and ε > ε max do 3 k ← 0, ε prev ← 2 ε; 4 while k ≤ k max and ε ≥ ε max and ε ≤ r ad ε prev do /* AD loop */ 5 k ← k + 1, ε prev ← ε; 6 Q H , ω ← solutions of T m (Q H , ω|P) = 0; 7 P ← solution of S m (P|Q H , ω) = 0; 8 for j ← 1 to m do /*
(K -ω 2 k M)φ k = 0 (19)
can be easily computed and are known to form a basis. Using this spatial information then seems a better choice than using random values. When trying to compute the n 0 -th NNM, it seems relevant to look for an additional spatial mode close to the next LNM. Starting from this point, the algorithm then adjusts this spatial mode by solving S m subproblem.

Branches of solutions using continuation with PGD/HBM modes enrichment

In this section a full algorithm to compute NNMs is proposed. It uses PGD combined to HBM with on the fly enrichment in conjunction with a continuation technique. A highly reduced and economical description of NNMs is obtained: one PGD mode is required near the n 0 -th LNM and new PGD modes are added only when it is necessary along the path. This dimension varying continuation scheme keeps the number of variables used to describe the branch as small as possible while energy grows and the system behaviour becomes more complex.

The main interest of this work does not consist in the continuation technique in itself but in the handling of a variable size of the solution space as PGD modes are added along the path. In what follows, the main ingredients of a classical prediction-correction continuation scheme are briefly defined before introducing the choices made to tackle the variable size.

To simplify the notations, we will denote y n the n-th point on the current path. This point gathers the PGD modes as well as the circular frequency information, P, Q H and ω and hence has size N m + (H + 1) m + 1:

y n ≡ {P n , Q H n , ω n } (20)
Moreover, the distance d(y n-1 , y n ) between two consecutive points y n-1 and y n is defined using separated weights for each variable:

d(y n-1 , y n ) = α 2 P ||P n -P n-1 || 2 + α 2 Q ||Q H n -Q H n-1 || 2 + α 2 ω |ω n -ω n-1 | 2 (21)
where α P , α Q and α ω are strictly positive constants which provide a mean to manage the different orders of magnitudes of the three variables. Each of them can be set to 1 but they can also be adjusted so as to give more or less importance to each of these variables by compensating for the difference in order of magnitudes.

Continuation process

As stated above, a very simple continuation procedure is used. More complex procedures can be found in [START_REF] Allgower | Introduction to Numerical Continuation Methods[END_REF]. Let us assume that n points y 1 , . . . , y n have already been computed. The next point y is sought using a prediction step which returns y pred ; this point is used as an initialization point for the correction step which, in our case, relies on the coupled PGD/HBM method described in Sec. 3.2 with a specific constraint equation c(q H , ω) = 0 linked to the chosen correction method. The global algorithm is depicted in Algorithm 2.

The prediction is based on a secant method:

y pred = y n + ∆s y n -y n-1 d(y n-1 , y n ) ( 22
)
∆s defines "how far" the point is predicted. It is usually adapted along the path:

for areas where the solution does not vary a lot, a large value should be used to save computation time; on the contrary, if the solution does vary "rapidly" small values provide better predictions and fewer correction steps. Here, we use a rough adaptation scheme depending on the convergence of the correction step and the error measure value ε (see Algorithm This secant prediction implies that at least 2 points have already been computed. The first point of the path is taken equal to the LNM associated to the studied NNM with null amplitude. Using notations introduced in Eq. ( 19) and ( 20): y 1 ≡ {φ n0 , 0 1,H+1 , ω n0 }. To find the second point, we choose to impose the amplitude of the j-th dof at t = 0: u j (t = 0) = u j0 with u j0 = 0. This prevents the algorithm from returning the null (and trivial) solution. This condition is denoted c 2 (q H , ω) in Algorithm 2 and defines the constraint equation c(q H , ω) = 0 of Algorithm 1. The coupled PGD/HBM method can then be applied to find y 2 starting from y 1 .

The correction step consists in applying the PGD/HBM method described in Algorithm 1 with y pred as a startpoint and using a specific constraint equation

c(q H , ω) involved in the time problem T m .
This equation is denoted c p (q H , ω) in Algorithm 2. It defines the correction method and can be freely chosen among the classical possibilities [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. An arclength correction is used for applications studied in Sec. 5:

c p (q H , ω ) = d(y n , y ) 2 -∆s 2 (23)
This condition also involves P which is the last value of spatial modes returned after the spatial problem S m when processing the current temporal problem T m .

As said above, these steps are rather classical for a prediction/correction based continuation procedure. The important difference lies in the possibility of increasing the problem dimension when stepping from y n to y n+1 as the -k max , ε max and r ad ;

-LNMs: Φ = [φ n0 , . . . , φ n ];

-Constraint equations: c 2 (q H , ω) and c p (q H , ω).

Data: Circular frequency of n 0 -th mode ω n0 .

1 y 1 ← {φ n0 , 0 1,H+1 , ω n0 }; // First point 2 c(q H , ω) ← c 2 (q H , ω),

Management of a variable problem size during the continuation

Having to deal with different problem dimensions along the path is not a problem in itself but it implies to decide for rules when y n points with different sizes have to be considered at the same time. This may happen during the prediction step when y n has a greater size than y n-1 and during the correction procedure when trying to enforce the arclength condition Eq. [START_REF] Urabe | Numerical computation of nonlinear forced oscillations by galerkin's procedure[END_REF]. Each time, it is possible either to pad the point with the smallest size with zeros for missing components or to delete the supernumerary components of the point with the bigger size. The choices for each of the two steps may not be identical.

Prediction step

If y n-1 and y n have different sizes m the secant prediction needs a slight adaptation because the distance given by Eq. ( 21) need summations of matrices that do not have the same sizes. The choice made is to pad y n-1 missing components with zeros instead of truncating y n . Hence, all computed data are taken into account to reach a better prediction.

Correction step

The same dimension issue happens when it comes to evaluate c p (q H , ω).

When y is described with more PGD modes than y n a truncation is applied for the current point. The use of an arclenght correction requires the distance between y n and y so it is evaluated in a space which has the size of y n . This does not impose any conditions on the new component introduced in y .

Application to two beam cases

In this section the main branch of the NNM1 (nonlinear normal mode tangent to the first linear mode for small energy levels) of two mechanical systems is built using the PGD/HBM based continuation method described in previous Both studied cases are illustrated in Fig. 1. In the first case a polynomial nonlinearity is chosen by introducing a cubic stiffness k nl to the spring. The second case examines a unilateral contact problem: the spring has a linear stiffness but there is a gap g between the free end of the beam and the spring.

The parameters required by Algorithm 2 are given in Table 2 for both studied (a) (b) cases.

∆m

Bending beam with a cubic spring at its free end

The underlying linear system is made of the mass and stiffness matrices of the beam plus the linear stiffness of the spring -i.e. K 11 = K 11 beam + k. Its linear modal properties are given at the last row in Table 1. The nonlinear

contribution is f nl (u) = [k nl u 3 1 (t) 0 N -1 ]
T with k nl = 9.2 10 5 N.m -3 . On an exploratory basis, H = 50 harmonics are considered. It is further showed that no more than ten harmonics are really needed for this problem.

Fig. 2 provides an overview of the results in a Frequency-Energy plot (FEP), which is a convenient way to display the NNM branches [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF][START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF]. On this figure and all the following ones N pt is the index of converged solution points through the continuation. We will focus on the 35 first computed points of the branch (N pt ≤ 35, solid line on Fig. 2) for the sake of clarity in the analysis.

Only 4 PGD modes are needed to describe the studied branch (Fig. 3 and solid line on Fig. 2) and only 6 PGD modes are required to build the whole plotted branch (including dotted line on Fig. 2). The NNM shape at null velocity gives interesting data as the studied NNM framework is the Rosenberg's one [START_REF] Rosenberg | On nonlinear vibrations of systems with many degrees of freedom[END_REF]: all dofs reach their maximum at the same time. On Fig. 2, the more energy grows the more the nonlinear spring takes importance against the bending inertia of the beam, generating localized curvature inversions near the free end. Kerschen qualitatively gives similar results for another cantilever beam with cubic spring in [START_REF] Kerschen | Computation of Nonlinear Normal Modes through Shooting and Pseudo-Arclength Computation[END_REF]. The error measure ε along the path is depicted along with the number of PGD modes on Fig. 3. Its evolution is consistent with the evolution of m: PGD modes are added when ε reaches the limit value ε max which, in turns, decreases the error. Moreover Fig. 2 shows the PGD modes are gradually added along the branch as expected with a smooth cubic nonlinearity. Note that the error for the first point N pt = 1 is not displayed as it is null and cannot appear in a logarithmic scale.

Fig. 4 provides an analysis of PGD spatial modes p j . PGD modes globally keep their shape through continuation. An interesting feature proper to the presented method is the presence of PGD modes which involve the participation of several LNMs. Such a "combined mode" is here the PGD mode 3, which gathers contributions of the first 6 LNMs. With such an approach, high order LMNs can be taken into account in a reduced range of PGD mode shapes: here 4 PGD modes are sufficient to describe a system behaviour which involves at Finally the main objectives of the method are met: to describe a NNM branch with a few variables, to keep the number of PGD modes as small as possible by incrementing m only when the convergence criterion is not fulfilled, and obviously to retrieve all physical data belonging to the case under study.

In the next section the method is tested on a non polynomial nonlinearity.

Bending beam with a unilateral contact

In this example the underlying linear system is the Euler-Bernoulli cantilever beam. Its modal properties are given in the first row of Table 1. The unilateral contact is modelled as follows:

f nl (u) =    [k(u 1 (t) -g) 0 N -1 ] T , if u 1 (t) > g 0 N , otherwise
with g = 2 cm. This law is kept non smooth and the contact is quite hard as k/ 3EI L 3 = 94.56% with 3EI L 3 the spring rate of the cantilever beam taken alone. A similar contact problem from the stiffness point of view is a contact between this beam and another identical one. Looking for periodical solutions for this class of problem usually requires a wide range of harmonics so the arbitrary value H = 70 is chosen. The consistency of this value is shown a posteriori in Fig. 9. NNM1 branch displayed in Fig. 6 is built with only 6 PGD modes and 199 solution points for the set of parameters given in Table 2. On the FEP given on Fig. 6, from N pt = 1 to N pt = 21, NNM1 is exactly equal to LNM1 -a pure fundamental cosine -as the beam does not touch the spring yet. Then a transition zone occurs for 22 ≤ N pt ≤ 85 where 5 PGD modes are added on 5 close points. Finally the rest of the branch is computed. Note that two internal resonances extracted from two others simulations are drawn with dotted lines. These data are obtained by manually modifying the parameters of Table 2 at the beginning or through the simulation. Analogous FEPs of related contact systems are studied by Moussi and Cochelin [START_REF] Cochelin | Numerical computation of nonlinear normal modes using HBM and ANM[END_REF][START_REF] Moussi | Nonlinear normal modes of a two degrees-of-freedom piecewise linear system[END_REF] for example. In this example again the error measure ε follows the evolution of m as shown on Fig. 7. Here 5 PGD modes are required to significantly decrease the error.

Other simulations with different continuation parameters add these 5 modes in one single point when it is allowed by using ∆m = 5. Note that by contrast with the cubic spring case all the PGD modes needed to build the branch are added in the same narrow area in the FEP. This shows that describing sharp or gradual changes of the system behaviour by the PGD/HBM continuation As shown on Fig. 8, for this problem PGD mode shapes do not significantly vary except the two last modes. The fifth PGD mode shape strongly changes between N pt = 82 and N pt = 199, mostly because of LNM4 participation.

However its final shape is about to be reached at N pt = 85, when PGD mode 6 420 is introduced, and PGD mode 5 does not evolve a lot from this point. The last PGD mode has remarkable coupling with high order LNMs. It condenses in a single mode a lot of modal information which illustrates the economical benefits of the presented method. The harmonic decomposition associated to each PGD mode is given on Fig. 9. The temporal part is complex as the contact issue generates a lot of important high harmonics but global trends can be considered. First a predominant cosine in PGD mode 1 drives the linear behaviour at the beginning of the range, then the amplitude of the whole harmonic set increases as energy grows. The free flight of the beam during a period can help to explain why low non negligible static coefficients a 0 appear because of the static term introduced in the nonlinear vector f nl (u).

Numerical considerations

The main objective of reduced model approaches is reached as few descriptors 435 are processed and stored with the presented approach, as shown in Table 3.

A classical HBM continuation requires N × (H + 1) descriptors by solution point whereas the implemented PGD/HBM continuation only needs m × (N + H + 1) ones. Obviously calculating a solution point via PGD/HBM is relevant if and only if m × (N + H + 1) < N × (H + 1). An easy proof gives the necessary 440 condition for the PGD/HBM method to be more interesting than a classical HBM approach: m < H + 1. In other words, the PGD/HBM method is useful for problems which require a large number of harmonics. In both examples the saving of variables with respect to a classical HBM can be assessed with the ratio 1 -m×(N +H+1) N ×(H+1) . The PGD reduction is very efficient, with 88% and 82% less variables for the cubic spring and the contact cases respectively.

We now focus on the underlying numerical solver and its iteration data when the subproblems are solved. The trust-region algorithm used requires function and abscissa tolerances -set on 10 -9 here -and allows to set a maximum number of iterations N ite max . The value N ite max = 100 is set but this limit is never reached (max. 54 iterations for the spatial problem of the cubic spring case, when PGD mode 4 is about to be added). In As expected the spatial problem is heavier and the given maxima are reached when a PGD mode is about to be added. N ite max could be an interesting control parameter for bigger problems as one can keep it very small and save computation time by "converting" convergence issues into new PGD modes or a reduced arclength ∆s.

Another performance indicator is the CPU time. The total computational time spent to build NNM1 of the cubic spring case was approximately 380s

with the PGD/HBM approach, against 10400s for a classic HBM, that is about 27 times less. Computations used the same continuation scheme, 36 solutions points were found in both cases and branches were computed up to E max ≈ 4.1 10 4 J; they were run on a desktop computer with MATLAB R2017a, 32

GB of RAM and an Intel Core i7 CPU. In addition to these points, the main advantage of the presented non-random method against the random variant is the reproducibility of its results.

Conclusion

The computation of NNMs using a PGD/HBM continuation technique with on the fly enrichment derived in this paper has been presented and successfully applied on two classical nonlinear examples with either smooth or non smooth nonlinearities. The only limitation of the method is to have access to nonlinear forces value f nl (u) at some specific points u. A significant decrease of the number of needed variables has been noticed when compared with a classical HBM approach; this demonstrates that NNM branches can be described in a highly reduced manner. The PGD modes obtained provide relevant physical data as spatial initializations are the LNMs of the system. The numerical aspects developed here are highly adaptable as many parameters are let free (cf. Table 2), giving a certain flexibility to the algorithm.

The presented approach could be used to process bigger structures with different kinds of nonlinearities. It could be interesting for systems which present many contact issues for example. Variants of the algorithm could then be considered in order to alleviate the cost of the spatial subproblem by using a progressive PGD method which only addresses a spatial problem when a PGD mode is added. On the fly impoverishment of the PGD modes set could also be considered so as to keep it as small as possible.

Finally, future prospects could include the computation of damped NNM or forced responses of a structure.

Algorithm 1 :

 1 PGD/HBM global algorithm Parameters: Values for m max , k max , ε max and r ad ; Definition of c(q H , ω) used in T m ; LNMs: Φ = [φ n0 , . . . , φ n ].

280Algorithm 2 :

 2 correction step managed by the PGD/HBM method described in Algorithm 1 can add up to ∆m modes to y n+1 (see line 7 in Algorithm 2 which defines m max used in Algorithm 1). This leads to continuation points with possibly different sizes which has to be dealt with. The termination criterion in Algorithm 2 (line 6) involves the maximum 285 number of PGD mode m max tot . This criterion could of course be advantageously replaced by a more meaningful criterion involving the frequency or the amplitude of a given dof for example. Continuation with on the fly PGD enrichment algorithm Parameters: Values for ∆m, m max tot , r ε , r mult s

  sections. Both systems are Euler-Bernoulli bending beams with one end clamped and a spring at its free end. Two different classes of conservative nonlinearities are then added and detailed. These models are generated using MATLAB directly: the need for nonlinear forces f nl expression makes it difficult to use the method in conjunction with a commercial FEM software. Let us define the values of the beam properties: Young modulus E = 210 GPa, density ρ = 7800 kg.m -3 , area of the square section S = 9 10 -4 m 2 , second area moment I = 6.71 10 -8 m 4 for an in-plane bending and length L = 1 m. The linear stiffness of the spring is k = 4 10 4 N.m -1 . The beam is discretized with N e = 20 elements by FEM with two dof per node: a transverse displacement u and a rotation θ. This leads to N = 40 dof once the clamping is taken into account. The first eigenvalues of the beam alone and beam plus the linear spring are given

Figure 1 :

 1 Figure 1: Euler-Bernouilli cantilever beam (a) with a transverse cubic spring at dof 1 (b) with a transverse gap between dof 1 and the linear spring.

Figure 2 :

 2 Figure 2: Beam with cubic spring -(a) Main branch of NNM1 in a FEP and (b) its mode shape at null velocity points for three continuation points. Squares: solution points where a PGD mode is added.

Figure 3 :

 3 Figure 3: Beam with cubic spring -(a) Number of PGD modes m and (b) error measure ε against point index Npt.

Figure 4 :

 4 Figure 4: Beam with cubic spring -Analysis of PGD spatial modes. First row: PGD mode shapes. Normalization is made at the maximum deflection point. Dotted line: shape when the mode is first introduced; Solid line: shape at the end of the branch (Npt = 35). Second row: Participation factors of LNMs in each PGD mode. Crosses: factors when the mode is introduced; Circles: factors at the end of the branch (Npt = 35).

Figure 5 :

 5 Figure 5: Beam with cubic spring -Amplitudes of (cosine) Fourier coefficients of PGD modes. Crosses: coefficients when the mode is first introduced; Circles: coefficients at the end of the branch (Npt = 35).

Figure 6 :

 6 Figure 6: Beam with unilateral contact -NNM1 main branch (solid) and internal resonances (dotted) in a FEP. Squares: points where PGD modes are added.

Figure 7 :

 7 Figure 7: Beam with unilateral contact -(a) m and (b) residue ε against Npt. Dotted line: convergence threshold.

Figure 8 :

 8 Figure 8: Beam with unilateral contact -Analysis of PGD spatial modes. Left column: PGD mode shapes. Normalization is made at the maximum deflection point. Dotted line: shape when the mode is first introduced; Solid line: shape at the end of the branch (Npt = 199). Right column: Participation factors of LMNs in each PGD mode. Crosses: factors when the mode is introduced; Circles: factors at the end of the branch (Npt = 199).

Figure 9 :

 9 Figure 9: Beam with unilateral contact -Amplitudes of (cosine) Fourier coefficients of PGD modes. Crosses: coefficients when the mode is first introduced; Circles: coefficients at the end of the branch (Npt = 199).

Table 3 :

 3 Comparaison between the number of descriptors needed by PGD/HBM -m(N + H + 1) + 1 by point -and classical HBM -N (H + 1) + 1 by point -for both examples.

Finally

  we compare the numerical efficiency of Algorithm 1 when using the proposed algorithm and when processing the spatial problem S m first. As explained in Sec. 3.2.4, when S m is solved before T m , the lack of temporal data to initialize the spatial problem is filled by using random values for Q H . Hence, uniformly distributed random numbers between 0 and 1 are used to initialize Q H . The cubic spring case is computed with exactly the same parameters except the random new PGD modes. In the non-random case, 35 attempts are made to compute the 35 solution points: line 17 in Alg. 2 is never called. To compare the result, among 25 draws of the random case we focus on the number of converged points computed with 35 attempts, the maximum energy reached E max and the number of iterations data of all the draws. The iteration number of the failed attempts are counted into the iteration number of the next converged point. The results are given in

30 Unlike

 30 the non-random computation, there are always failed attempts. The number of computed iterations is far bigger. In particular the temporal problem is heavily penalized by the random coefficients. Otherwise no point reaches the final energy of the non-random branch. Some draws stay at a very low level of energy and 24 out of 25 draws have solutions points which go back into the NNM branch. This shows initializing new modes with random values leads to non robust predictions.

  Inside the alternated directions algorithm loop, temporal and spatial problems are solved (lines 6 and 7 respectively). Once a new P matrix is computed, each spatial mode p j is normalized as stated in Sec. 3.2.1.Once the loop is broken, if the convergence criterion ε ≤ ε max is not fulfilled, a new PGD mode is added (from line 13): m is increased and a new column is added to P matrix while a new row is added to Q H matrix. The choice of these additional values is important as it will greatly influence the convergence of the next resolution loop. Grolet and Thouverez[START_REF] Grolet | On the use of the proper generalised decomposition for solving nonlinear vibration problems[END_REF] chose to add random values for Q H , and processed the spatial problem before the temporal one. We propose here to process the temporal problem first based on the new spatial mode introduced φ n0+m , with n 0 being the index of the n 0 -th NNM. Indeed, we do not know much of the relevant temporal information to inject: having to initialize Q H is a real difficulty. However, the LNMs φ k solutions of:

	Q H 0 1,H+1	;

Mode normalization */ 9 p j ← p j /||p j ||; 10 end 11 ε ← max t∈It (||R(u(t))||/||Ku(t)||); 12 end 13 if ε ≥ ε max then /* Next mode initialization */ 14 P ← [P, φ n0+m ], Q H ← 15 m ← m + 1; 16 end 17 end 18 Return m, P, Q H , ω and ε;

  1 in Sec. 3.2.4) once a new point was computed: if the correction step did not provide a satisfying new point, ∆s is decreased by multiplying it by r div

s < 1 (see line 17) leading further to a new prediction point y pred closer to y n ; if the correction step provided a satisfying point for which the error measure is less than r ε ε max with r ε < 1, ∆s is multiplied by r mult s > 1 (see lines 13 and 14); otherwise, ∆s is kept constant for the next prediction step.

Table 1 :

 1 in Table1. First eigenvalues of the clamped-free beam alone and with the linear spring attached to the free end.

		i	1	2	3	4	5	6
	ω i (rad.s -1 )	Beam alone With spring 217 999 2769 5420 8960 13386 157 988 2765 5419 8959 13385

Table 2 :

 2 Algorithm 2 required parameters.

			mmax tot	rε	r mult s	r div s	∆s
	Cubic spring case	1	10	0.1	1	0.5	1.5
	Contact case	5	10	0.5	1.07	0.5	1
		kmax	εmax	r ad	y 2 : u 1 (t = 0)	{α P , αq, αω}
	Cubic spring case	10	10 -3	0.8	10 -2 m		{1, 10 -1 , 10 -4 }
	Contact case	10	10 -3	0.8	5 10 -4 m	{0, 10 -6 , 1}

Table 4 :

 4 Table. 4 are given mean and maximum numbers of solver iterations for temporal and spatial problems through the continuation, for both examples. Mean, standard deviation and maximum of iteration numbers for Tm and Sm problems, through the continuation.

		Mean	Standard	Maximum
		iteration number deviation iteration number
	Tm	6.23	4.89	22
	Cubic spring case
	Sm	16	10.92	54
	Tm	2.18	1.53	16
	Contact case
	Sm	9.77	8.68	49

Table 5 .

 5

				Standard
		Min.	Mean		Max.
				deviation
	Converged points	5	21.2	9.31	31
	Emax (J)	3.74 10 -10	3.71 10 3	4.4 10 3	1.99 10 4
	Iteration number for Tm	3	209.93	389.59	3000
	Iteration number for Sm	3	110.62	192.03	1463

Table 5 :

 5 Mean, standard deviation and extrema of simulation features for 25 random draws.

scheme with on the fly enrichment is feasible.

frequency harmonics are always significant through the continuation. Note that

Acknowledgements & Declaration of interest

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Declarations of interest: none.