On the Use of the Proper Generalised Decomposition for Solving Nonlinear Vibration Problems

This paper presents the use of the so called Proper Generalized Decomposition method (PGD) for solving nonlinear vibration problems. PGD is often presented as an a priori reduction technique meaning that the reduction basis for expressing the solution is computed during the computation of the solution itself. In this paper, the PGD is applied in addition with the Harmonic Balance Method (HBM) in order to ﬁnd periodic solutions of nonlinear dynamic systems. Several algorithms are presented in order to compute nonlinear normal modes and forced solutions. Application is carried out on systems containing geometrical nonlinearity and/or friction damping. We show that the PGD is able to compute a good approximation of the solutions event with a projection basis of small size. Results are compared with a Proper Orthogonal Decomposition (POD) method showing that the PGD can sometimes provide an optimal reduction basis relative to the number of basis components.


INTRODUCTION
The objective of this paper is to evaluate the efficiency of the proper generalized decomposition (PGD) for solving nonlinear vibration problems in the context of reduced order modeling.Most of the reduction method (e.g.krylov subspace [1,2], proper orthogonal decomposition [3], dual modes [4]) aim at reducing the model by first computing a set of vector which will constitute a reduced basis for solving the problem.This basis is computed relatively to some information already known about the system behaviour.For instance, the Proper Orthogonal Decomposition (POD) uses preliminary simulations (or mesurements) to extract a basis by using data analysis tools such as singular value decomposition.Thus those methods are often called a posteriori reduction methods.Contrary to those methods, the PGD aims at finding a reduced basis during the computation of the solution itself without any previous information on the system, thus making it an a priori reduction method.PGD was introduced under a different name (radial time-space approximation) by P.Ladeveze in the LATIN method [5].In the last decade, A.Nouy wrote a series of paper on the use of PGD for solving linear stochastic problems [6][7][8] and first order (in time) linear partial differential equations [7,9].The core of the method resides in separation of variable (also called radial decomposition by P.Ladeveze [5]), alternate Galerkin projections in time and space and fixed point computation.In this paper the PGD is used for finding steady state solutions of nonlinear second order differential equation arising for instance when modeling nonlinear structural dynamic problems with finite element methods.We will show how the harmonic balance method (HBM) can be used in conjunction with PGD for finding periodic solution of the problem.A comparison will be drawn between PGD and POD on a numerical example.The paper is organized as follows: the first section presents generically the kind of nonlinear problems studied in this paper along with the description of Galerkin projections and the particular case of the Harmonic Balance Method.In the second section the POD and the PGD are presented in detail along with the modifications used to find periodic solutions.The last section is dedicated to numerical examples and to the comparison between POD and PGD.

NONLINEAR VIBRATION PROBLEMS AND GALERKIN PROJECTIONS
In this section we presents the kind of nonlinear problems studied in the paper and the application of Galerkin projection on such problems.

Nonlinear vibrations
Problems considered in this paper are related to the study of nonlinear vibration of structural systems.After discretization of the problem (e.g. by finite elements methods) one have to solve a systems of n second order nonlinear differential equations, where n is the number of degree of freedom (dof) of the model.Generally speaking the system can be put in the following form: where M, C, K are respectively the mass, damping and stiffness matrices of size n × n, u u u is the vector of dof, f f f (t) is the vector of excitation forces which will be assumed to be harmonic with period T , and g g g is the vector of nonlinear forces acting on the system.In this paper we will consider a polynomial forces of maximum degree 3 corresponding to geometric nonlinearity, this force can be expressed as follows [4]: where k (1) , k (2) , k (3) are tensors of coefficient computed by the finite element model.

Galerkin projections
Here we briefly recall the results given by Galerkin projection applied to the system in Eq.( 1) as they will be used through the paper.Galerkin projection assumes a separation of variables such that a solution u u u can be expressed under the form u u u = Pq q q(t) = ∑ r i=1 p p p i q i (t) with P a matrix of size n × r and q q q(t) a vector of size r.A time (resp.space) projection basis of size r is then set and projections are computed in order to derive a purely spatial (resp.temporal) system.

Time projections
In the case of a projection in time, the basis q q q is set and we search for P such that the following holds for 1 ≤ i ≤ r: r ∑ j=1 [ I q i q j M + I q i q j C + I q i q j K]p p p j + I q i g g g(Pq q q, P q q q) = I q i f f f (3) which can be rewritten as: where p p p = [p p p T 1 , . . ., p p p T r ] T , g g g( p p p) = [ I q i g g g] 1≤i≤r and f f f = [ I q i f f f ] 1≤i≤r are vector of size nr and I k = [ I q i d k q j dt k ] are matrices of size r × r.The approximated solution u u u a (t) = Pq q q(t) is then computed by solving the nr nonlinear algebraic equations Eq.( 4).A typical use of Galerkin projection in time is the harmonic balance method described later in the paper.

Space Projections
In the case of a projection in space, the basis P is set and we search for q q q(t) such that the following holds: P T (MP q q q + CP q q q + KPq q q + g g g(Pq q q, P q q q)) = P T f f f (5) or M r q q q + C r q q q + K r q q q + g g g r (q q q, q q q) = f f f r (6) where A r = P T AP, A = M, C, K are reduced matrices of sizes r × r and g g g r , f f f r reduced vector of size r.The approximated solution is then computed by solving the r non linear differential equations Eq.( 6).

Harmonic Balance Method
The harmonic balance method is a particular case of Galerkin projection in time, and a widespread way to find periodic solutions of nonlinear systems such as Eq.(1) [10][11][12][13].
The basis used for this projection is a truncated Fourier basis up to H harmonics T T T (t) = [1, cos(ωt), sin(ωt), . . ., cos(Hωt), sin(Hωt)] and the time interval for integration is reduced to I = [0, 2π/ω].In this case the integrals I k Eq.( 4) are simplified into I 0 = I d (identity matrix), to the following algebraic system with n = 2H + 1 equations: )n, g g g and f f f vectors of size (2H + 1)n corresponding respectively to the nonlinear forces and the excitation forces in the frequency domain.When no analytical expressions can be derived for the nonlinear forces in the frequency domain g g g, an alternate frequency time (AFT) procedure is used to approximate g g g [14].
Let distinguish the two cases of forced and free solutions.In the case of forced solution, the frequency ω is taken equal to the frequency of the harmonic force f f f , that is ω = 2π/T .In the case of free solution, the excitation forces are null, and the frequency ω is considered as an unknown.In this case a phase equation is added (for instance null velocities at t = 0) along with an equation controlling the amplitude of the displacements in order to avoid obtaining the null solution.

MODEL REDUCTION METHODS: POD AND PGD Proper Orthogonal Decomposition
The proper orthogonal decomposition is a model reduction technique which has been employed in various domain such as fluids [15] or structural [3] dynamics.It consist in extracting a reduction basis from data acquired during numerical simulations or measurements.From a solution u u u = u u u(x x x,t) the reduction basis V(x) is computed such that it maximize the average in time of the inner product u u u(x,t), V .The result is a set of proper orthogonal modes (POM) V V V i (x x x) associated to a set of proper orthogonal values (POV) λ i .The participation of POM number i in the solution is given by the ratio λ i / ∑ j λ j .A convenient way to compute the reduced basis in the discrete case is the so called snapshot methods [3].From a set of data U = u u u(x x x i ,t j ), the POD is obtained by performing a singular value decomposition of the matrix U, i.e.U = VSW T where V is the matrix of proper orthogonal modes and S the matrix containing the singular values on its diagonal.When carrying the POD, one observe that the magnitude of the POV is rapidly decreasing by several order of magnitude.This fact allow for model reduction by building a reduction basis P composed of only the r most significant POMs.The reduction is then carried out by searching an approximated solution u u u a of Eq.( 1) under the form u u u a (t) = Pq q q(t), and performing a Galerkin projection in space as described in Eq.( 6).The resolution of the reduced r differential equations finally allow to find the approximated solution u u u a .
Although the POD can give good approximations, its majors drawback resides in the fact that data is needed to build the reduced model.If the data comes from numerical simulations, this mean that one has to solve the n dof system in Eq.( 1) (e.g. by numerical integration, or HBM) for various loading case f f f .Those simulations can time consuming and are only valid for an excitation force f f f lying in the neighborhoods of the load cases used to derived the reduced basis.

Proper generalized decomposition
The PGD is a numerical method allowing to reduce models during the computation of the solution itself.The method and the related algorithms have been described in numerous papers [6][7][8]16] and especially in [9].It is based on variable separation and on an iterative procedure including alternate Galerkin projection in space and in time.In this section we describe how we applied the PGD to find approximated solutions of system in Eq. (1).Two variants are presented, namely optimal Galerkin PGD (oPGD) and progressive PGD (pPGD), along with their corresponding algorithms.

Optimal
Galerkin PGD During the optimal Galerkin PGD, solutions are sought under the form u u u m (t) = ∑ m k=1 p p p k q k (t) = Pq q q(t).The matrix P and the time evolution vector q q q(t) are computed such that the following criterion hold simultaneously: B(u u u m , P * q q q(t)) = L(P * q q q(t)), ∀P * ∈ E oPGD S (8a) where E oPGD S is the set of real-valued matrix of size n × m, E oPGD T is the set of periodic functions vector of size m with period T = 2π/ω continuously derivable and B, L are two applications defined by the following relations: We define the two applications S S S m : q q q → P such that Eq.(8a) holds (the space problem), and T T T m : P → q q q such that Eq.(8b) holds (the time problem).The couple (P m , q q q m (t)) is given when stationarity has been reached, i.e when P m = S S S m (q q q m (t)) and q q q m (t) = T T T m (P m ).
When applied to problem in Eq.( 1), this procedure reduces to perform Galerkin projections in time and in space alternatively.The space problem corresponds to solve equation Eq.( 4) and the time problem corresponds to solve Eq.( 6).
Algorithm 1 gives the algorithm corresponding to oPGD.At each step k, the computation of the couple (P m , q q q m (t)) requires solving nm nonlinear algebraic equation (space problem in Eq.( 4)) and m nonlinear differential equations (time problem in Eq.( 6)).The time problem is quite small (only m equations) and can solved by HBM efficiently.The resolution of the space problem can become very costly when dealing with large number of

Algorithm 1 ALGORITHM FOR OPTIMAL PGD
for m = 1 to m max do Initialize q q q m (t) for k = 1 to k max do Compute P m = S S S m (q q q m m m ) by solving Eq.( 4) Orthonormalize P m [optional] Compute q q q m (t) = T T T m (P m ) by solving Eq.( 6) end for Set u u u m (t) = P m q q q m (t) and check convergence end for dof n.To improve performances, solutions of the space problem can be searched for in a smaller space, for instance by searching solutions as a linear combination of linear modes shapes [17].
Checking convergence is done by monitoring the relative change ε between two solutions given by ε = u u u m+1 − u u u m / u u u m .
Progressive PGD During the progressive PGD, solutions are computed iteratively.At step m, the previously computed solution u u u m−1 is completed with a new term p p p m q m (t) such that u u u m (t) = u u u m−1 (t) + p p p m q m (t) with u u u m−1 = ∑ m−1 k=1 p p p k q k (t) = P m−1 q q q m−1 (t).The vector p p p and the time evolution q(t) are computed such that the following orthogonality criterion hold simultaneously: B(u u u m , p p p * q(t)) = L(p p p * q(t)), ∀p p p * ∈ E pPGD S B(u u u m , p p pq * (t)) = L(p p pq * (t)), ∀q * (t) ∈ E pPGD T (10) where E pPGD S is the set of real-valued vector of size n, E pPGD T is the set of periodic function of period T = 2π/ω, B(u u u, v v v) and L(v v v) are given in Eq.( 9).We define the two applications S m−1 : q → p p p such that Eq.(10.1)holds (the space problem), and T m−1 : p p p → q such that Eq.(10.2) holds (the time The couple (p p p m , q m (t)) is given when stationarity has been reached, i.e when p p p m = S m (q m (t)) and q m (t) = T m (p p p m ).
This procedure again results in Galerkin projections, but the projection is carried out on only one vector of the basis (either p p p m or q m (t)), leading to the following equation for the time (Eq.(11)) and the space problem (Eq.( 12)): p T m (Mp p p m qm + Cp p p m qm + Kp p p m q m ) + g g g(Pq q q + p p p m q m ) = p T m f f f − p T m (MP q q q + CP q q q + KPq q q) (11) I q m qm M + I q m qm C + I q m q m K]p p p m + I q m g g g(Pq q q + p p p m q m ) = Algorithm 2 ALGORITHM FOR PROGRESSIVE PGD for m = 1 to m max do Initialize q m (t) for k = 1 to k max do Compute p p p m = S m (q m ) by solving Eq.( 11) Orthonormalize p p p m relative to p p p 1 , . . ., p p p m−1 [optional] Compute q m (t) = T m (p p p m ) by solving Eq.( 12) end for if update then Update time functions q q q m = T T T m (P m ) by solving Eq.( 6) Set u u u m (t) = P m q q q m (t) and check convergence else Set u u u m = u u u m−1 + p p p m q m (t) and check convergence end if end for The algorithm used to implement pPGD is depicted in Algorithm 2. At each step m the computation of the new term to be added consist in solving k max times an algebraic problem with n unknowns and a differential equation with one variable, where k max is the number of iteration needed to achieve stationarity.Since we are searching for periodic solutions, HBM is used to solve time problems in Eq. (12).The number of pPGD step m max can be determined by monitoring the relative change in norm between two solutions (as in oPGD), or can be set a priori by the user.The results of the pPGD can be improved by updating the whole set of time functions q q q after a new couple (p p p m , q m ) has been computed.This update is done by solving the time problem associated with the mapping T T T : P → q q q(P), i.e. by solving Eq. (6).pPGD can also be improved by computing several term at the same times.This method is based on the progressive PGD, but instead of computing one term (p p p m , q m ) at each iteration we choose to compute r terms simultaneously (with r remaining relatively small) by solving spacial and time problems derived from oPGD.

Remarks about the time and space problems
In all PGD methods described above, time problems (Eqs.(6) or ( 12)) are solved by HBM with H harmonics in order to find steady   6) or ( 12)) is transformed into a system of nonlinear algebraic equation having the form of Eq.( 7) Since time problems only consists in a small set of differential equations of size m (m ≤ n), a large number of harmonics H can be retained, leading to solve a system of algebraic equations with m(2H + 1) unknowns at each iteration of oPGD (or at each update of pPGD).
When building space problems (Eqs.(4) or ( 11)) integrals involving q(t), q(t) and q(t) are easily computed in the frequency domain.Indeed, the solution of time problems by HBM gives q i (t) = T(t) q q q i where T(t) contains the 2H + 1 vectors of Fourier basis.Integrals can then be approximated by: where

NUMERICAL APPLICATION Beam featuring geometric nonlinearities
The system considered here consists in an Euler-Bernoulli beam as depicted in Fig. 1.Since the beam is clamped/clamped, interactions can occur between axial and transverse displacements leading to geometric nonlinearities.The numerical values of physical and geometrical parameters are given in Table 1.The beam is discretized by mean of finite element method with three dof per node (axial displacement u, transverse displacement v, rotation θ ) and n e = 20 elements.

FIGURE 2. DISPLACEMENT OF THE CENTER OF THE BEAM AS A FUNCTION OF FREQUENCY FOR THE REFERENCE SOLU-TION
After assembling and adding a damping term C = 3M, the resulting vibration equation for the beam is in the form of Eq.( 1) with n = 60 dof and with the nonlinear term g g g arising from geometric nonlinearities given by Eq.( 2).

Forced Solution
In this section we consider that the beam is excited at its center by an purely transverse harmonic force f (t) = A cos(ωt) with A = 200N.The reference solution u u u re f is computed by HBM with H = 3 harmonics for ω 2π from 155 to 180 Hz (upward direction) and for ω 2π from 190 to 162 Hz(downward direction) see Fig. 2.
For POD, at each frequency step, a POD is computed from the reference solution, and an approximated solution u u u m POD is computed by projecting the system onto a basis with m elements (this ideal version of the POD have no practical interests but for comparison with other methods).
For optimal PGD the initialization of q q q ω 1 m (t) (Alg.1) is chosen randomly for the first frequency step, and for all subsequent steps it is initialized with the previously computed value q q q ω i−1 m .The same strategy is used for progressive PGD (Alg.2) along with the update step each time a new vector has been computed.Solutions computed by oPGD (resp.pPGD)with m terms are denoted u u u m oPGD (resp.uu u m pPGD ).Since no theoretical results are available for proving the convergence of our fixed point problem, we did some preliminary studies about the numerical convergence of the fixed point problem.It turns out that the convergence is actually very fast and can be achieve most of the time within 3 iterations.Consequently, we choose to set the number of maxi-mum iterations to k max = 3. Finally the threshold for the relative change is set to ε = 0.005.
Solution are compared relative to the reference solution by using the relative difference ε r = u u u re f (t) − u u u(t) / u u u re f .Figure 3 depicts the error relative to the reference solution for various number of retained mode m ≤ 4 for the three methods POD, oPGD and pPGD.For POD the error is clearly decreasing each time a new mode is retained.For oPGD we can see that increasing the number of maximum mode does not always decrease the error since the number of mode is actually controlled by the relative change in solution ε.The same remark holds for progressive PGD.We can see that oPGD gives very similar result than POD in term of error, and can sometimes provides better result than POD for a fixed number of modes, for example oPGD with m = 2 produces a smaller error than POD with the same number of mode (this is also happening for m = 4 see Fig. 3) (a similar fact can be observed for pPGD with m = 2 which produces a smaller error than POD with m = 2).We believe this phenomenon is due to the fact that oPGD introduces axial displacements sooner that POD (recall that nonlinearity arise from axial/transverse coupling), axial displacements are present in the second mode of oPGD and in the third mode of POD (see Figs. 4,6).
In term of performance, solving the spatial problem is very costly since one need to solve a nonlinear algebraic set of n × m with Newton Raphson method.One way to improve performance would be to consider an incremental or a Newton linearization of the space problem as proposed in [18].Another mean of improving performance is to seek the solution of the spatial problem in a smaller subspace as explained in [17].

Free Solutions
In this section we search for free solutions of Eq.( 1) i.e. nonlinear normal modes (NNM).We only concentrate on the first mode of vibration but this procedure can be applied for any mode.A reference solution is computed by HBM with H = 3 harmonics and with cosine term only, this last condition corresponds to a phase condition in which the velocities are set to zeros for t = 0.As in the forced case, a POD basis with m components is computed from the reference solution for each frequency step, and an approximated solution is obtained by projecting the system onto this basis.Finally, a PGD approximation of the NNM is computed by optimal PGD.For the first frequency step, the solution is initialized with the linear mode and then the amplitude of the time evolution of the first mode is controlled in order to avoid null solution.The parameter used for this computation with oPGD are k max = 3 and ε = 0.005.
Figure 7 shows the backbone curve of the first mode obtained for the various methods (HBM,POD,oPGD) and Fig. 8 show the relative error with respect to the reference solution.Once again we see that oPGD gives results very close to POD  when using only one mode in the decomposition.For m = 2 we can see that oPGD produces a smaller error than POD, and that this solution is very close to the reference one.

CONCLUSION
In this paper proper generalized decomposition was applied to solve nonlinear vibration problem.During PGD the solution of the space problems was found by using the harmonic balance method thus allowing to search for periodic solutions.Several algorithm were applied through a numerical example consisting in beam featuring geometric nonlinearity, modeled by finite element method.Comparison were drawn between solution computed by HBM, POD and PGD.Results show that optimal Galerkin gives very similar result than POD and can sometimes produce a smaller error than POD for the same number of mode in the solution.Although PGD gives results, work has to be done on performance and particularly on solving the space problem which can be very costly for large dof systems.

FIGURE 1 .
FIGURE 1. SCHEME OF THE CLAMPED/CLAMPED BEAM AND REFERENCE ELEMENT

FIGURE 4 .
FIGURE 4. DECOMPOSITION OF THE SOLUTION COMPUTED BY oPGD FOR m = 4 AT ω 2π = 180 (shape have been normed for the purpose of comparison)

FIGURE 7 .FIGURE 8 .
FIGURE 7. BACKBONE CURVE OF THE FIRST NONLINEAR MODE COMPUTED BY HBM, POD and oPGD

TABLE 1 .
NUMERICAL PARAMETERS USED IN THE BEAM MODEL