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

Aurélien Grolet, Fabrice Thouverez

► To cite this version:

Aurélien Grolet, Fabrice Thouverez. On the Use of the Proper Generalised Decomposition for Solving Nonlinear Vibration Problems. ASME 2012 International Mechanical Engineering Congress and Exposition, 2012, Houston, United States. 10.1115/IMECE2012-87538. hal-02117586

HAL Id: hal-02117586

https://hal.science/hal-02117586

Submitted on 14 Dec 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ON THE USE OF THE PROPER GENERALISED DECOMPOSITION FOR SOLVING NONLINEAR VIBRATION PROBLEMS

Aurelien Grolet *

Laboratoire de Tribologie et Dynamique des Systemes
Ecole Centrale de Lyon
Ecully, 69130
France
aurelien.grolet@ec-lyon.fr

Fabrice Thouverez

Laboratoire de Tribologie et Dynamique des Systemes
Ecole Centrale de Lyon
Ecully, 69130
France
fabrice.thouverez@ec-lyon.fr

## ABSTRACT

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 find 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 measurements) 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–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

<footnote>*Address all correspondence to this author.</footnote>
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 present 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:

$$ M\ddot{u} + C\dot{u} + Ku + g(u, \dot{u}) = f(t) \quad (1) $$

where $M, C, K$ are respectively the mass, damping and stiffness matrices of size $n \times n$, $u$ is the vector of dof, $f(t)$ is the vector of excitation forces which will be assumed to be harmonic with period $T$, and $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]:

$$ g_i(u) = \sum_j k_{ij}^{(1)} u_j + \sum_{j,k} k_{ijk}^{(2)} u_j u_k + \sum_{j,k,l} k_{ijkl}^{(3)} u_j u_k u_l, \quad 1 \le i \le n \quad (2) $$

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$ can be expressed under the form $u = Pq(t) = \sum_{i=1}^{r} p_i q_i(t)$ with $P$ a matrix of size $n \times r$ and $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$ is set and we search for $P$ such that the following holds

for $1 \le i \le r$:

$$ \sum_{j=1}^{r} \left[ \int_I q_i \ddot{q}_j \mathbf{M} + \int_I q_i \dot{q}_j \mathbf{C} + \int_I q_i q_j \mathbf{K} \right] \mathbf{p}_j + \int_I q_i \mathbf{g}(\mathbf{P} \mathbf{q}, \mathbf{P} \dot{\mathbf{q}}) = \int_I q_i \mathbf{f} \quad (3) $$

which can be rewritten as:

$$ (\mathbf{I}_2 \otimes \mathbf{M} + \mathbf{I}_1 \otimes \mathbf{C} + \mathbf{I}_0 \otimes \mathbf{K}) \tilde{\mathbf{p}} + \tilde{\mathbf{g}}(\tilde{\mathbf{p}}) = \tilde{\mathbf{f}} \quad (4) $$

where $\tilde{\mathbf{p}} = [\mathbf{p}_1^T, \dots, \mathbf{p}_r^T]^T$, $\tilde{\mathbf{g}}(\tilde{\mathbf{p}}) = [\int_I q_i \mathbf{g}]_{1 \le i \le r}$ and $\tilde{\mathbf{f}} = [\int_I q_i \mathbf{f}]_{1 \le i \le r}$ are vector of size $nr$ and $I_k = [\int_I q_i \frac{d^k q_j}{dt^k}]$ are matrices of size $r \times r$. The approximated solution $\mathbf{u}_a(t) = \mathbf{P} \mathbf{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 $\mathbf{P}$ is set and we search for $\mathbf{q}(t)$ such that the following holds:

$$ \mathbf{P}^T (\mathbf{M} \ddot{\mathbf{q}} + \mathbf{C} \dot{\mathbf{q}} + \mathbf{K} \mathbf{P} \mathbf{q} + \mathbf{g}(\mathbf{P} \mathbf{q}, \mathbf{P} \dot{\mathbf{q}})) = \mathbf{P}^T \mathbf{f} \quad (5) $$

or

$$ \mathbf{M}_r \ddot{\mathbf{q}} + \mathbf{C}_r \dot{\mathbf{q}} + \mathbf{K}_r \mathbf{q} + \mathbf{g}_r(\mathbf{q}, \dot{\mathbf{q}}) = \mathbf{f}_r \quad (6) $$

where $\mathbf{A}_r = \mathbf{P}^T \mathbf{A} \mathbf{P}, \mathbf{A} = \mathbf{M}, \mathbf{C}, \mathbf{K}$ are reduced matrices of sizes $r \times r$ and $\mathbf{g}_r, \mathbf{f}_r$ reduced vector of size $r$. The approximation 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–13]. The basis used for this projection is a truncated Fourier basis up to $H$ harmonics $T(t) = [1, \cos(\omega t), \sin(\omega t), \dots, \cos(H\omega t), \sin(H\omega t)]$ and the time interval for integration is reduced to $I = [0, 2\pi/\omega]$. In this case the integrals $I_k$ Eq.(4) are simplified into $I_0 = I_d$ (identity matrix), $\mathbf{I}_1 = \mathbf{D} = \omega \operatorname{diag}(0, \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix}, \dots, H \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix})$, $\mathbf{I}_2 = \mathbf{D}^2$ leading to the following algebraic system with $\tilde{n} = 2H + 1$ equations:

$$ \mathbf{L}(\omega) \tilde{\mathbf{p}} + \tilde{\mathbf{g}}(\tilde{\mathbf{p}}) = \tilde{\mathbf{f}} \quad (7) $$
with $L = D^2 \otimes M + D \otimes C + I_d \otimes K$ a matrix of size $(2H+1)n \times (2H+1)n$, $\tilde{g}$ and $\tilde{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 $\tilde{g}$, an alternate frequency time (AFT) procedure is used to approximate $\tilde{g}$ [14].

Let distinguish the two cases of forced and free solutions. In the case of forced solution, the frequency $\omega$ is taken equal to the frequency of the harmonic force $f$, that is $\omega = 2\pi/T$. In the case of free solution, the excitation forces are null, and the frequency $\omega$ 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 consists in extracting a reduction basis from data acquired during numerical simulations or measurements. From a solution $u = u(x,t)$ the reduction basis $V(x)$ is computed such that it maximizes the average in time of the inner product $\langle u(x,t), V \rangle$. The result is a set of proper orthogonal modes (POM) $V_i(x)$ associated to a set of proper orthogonal values (POV) $\lambda_i$. The participation of POM number $i$ in the solution is given by the ratio $\lambda_i/\sum_j \lambda_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(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 observes that the magnitude of the POV is rapidly decreasing by several order of magnitude. This fact allows 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_a$ of Eq.(1) under the form $u_a(t) = Pq(t)$, and performing a Galerkin projection in space as described in Eq.(6). The resolution of the reduced $r$ differential equations finally allows to find the approximated solution $u_a$.

Although the POD can give good approximations, its major drawback resides in the fact that data is needed to build the reduced model. If the data comes from numerical simulations, this means that one has to solve the $n$ dof system in Eq.(1) (e.g. by numerical integration, or HBM) for various loading case $f$. Those simulations can time consuming and are only valid for an excitation force $f$ lying in the neighborhoods of the load cases used to derive 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–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_m(t) = \sum_{k=1}^m p_k q_k(t) = Pq(t)$. The matrix $P$ and the time evolution vector $q(t)$ are computed such that the following criterion hold simultaneously:

$$B(\mathbf{u}_m, \mathbf{P}^* \mathbf{q}(t)) = L(\mathbf{P}^* \mathbf{q}(t)), \quad \forall \mathbf{P}^* \in \mathcal{E}_S^{\text{oPGD}} \quad (8a)$$

$$B(\mathbf{u}_m, \mathbf{P} \mathbf{q}^*(t)) = L(\mathbf{P} \mathbf{q}^*(t)), \quad \forall \mathbf{q}^*(t) \in \mathcal{E}_T^{\text{oPGD}} \quad (8b)$$

where $\mathcal{E}_S^{\text{oPGD}}$ is the set of real-valued matrix of size $n \times m$, $\mathcal{E}_T^{\text{oPGD}}$ is the set of periodic functions vector of size $m$ with period $T = 2\pi/\omega$ continuously derivable and $B, L$ are two applications defined by the following relations:

$$B(\mathbf{u}, \mathbf{v}) = \int_I \mathbf{v}^T (\mathbf{M}\ddot{\mathbf{u}} + \mathbf{C}\dot{\mathbf{u}} + \mathbf{K}\mathbf{u}) + \int_I \mathbf{v}^T \mathbf{g}(\mathbf{u}, \dot{\mathbf{u}}) \quad (9) \\ L(\mathbf{v}) = \int_I \mathbf{v}^T \mathbf{f}$$

We define the two applications $S_m : q \mapsto P$ such that Eq.(8a) holds (the space problem), and $T_m : P \mapsto q$ such that Eq.(8b) holds (the time problem). The couple $(P_m, q_m(t))$ is given when stationarity has been reached, i.e when $P_m = S_m(q_m(t))$ and $q_m(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_m(t))$ requires solving $nm$ nonlinear algebraic equations (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 be 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<sub>max</sub> **do**

Initialize **q**<sub>m</sub>(t)

**for** k = 1 **to** k<sub>max</sub> **do**

Compute **P**<sub>m</sub> = **S**<sub>m</sub>(**q**<sub>m</sub>) by solving Eq.(4)

Orthonormalize **P**<sub>m</sub> [optional]

Compute **q**<sub>m</sub>(t) = **T**<sub>m</sub>(**P**<sub>m</sub>) by solving Eq.(6)

**end for**

Set **u**<sub>m</sub>(t) = **P**<sub>m</sub>**q**<sub>m</sub>(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**<sub>m+1</sub> - **u**<sub>m</sub>||/||**u**<sub>m</sub>||.

**Progressive PGD** During the progressive PGD, solutions are computed iteratively. At step m, the previously computed solution **u**<sub>m-1</sub> is completed with a new term **p**<sub>m</sub>**q**<sub>m</sub>(t) such that **u**<sub>m</sub>(t) = **u**<sub>m-1</sub>(t) + **p**<sub>m</sub>**q**<sub>m</sub>(t) with **u**<sub>m-1</sub> = Σ<sub>k=1</sub><sup>m-1</sup> **p**<sub>k</sub>**q**<sub>k</sub>(t) = **P**<sub>m-1</sub>**q**<sub>m-1</sub>(t). The vector **p** and the time evolution **q**(t) are computed such that the following orthogonality criterion hold simultaneously:

$$
\begin{align}
B(\mathbf{u}_m, \mathbf{p}^* q(t)) &= L(\mathbf{p}^* q(t)), \quad \forall \mathbf{p}^* \in \mathcal{E}_S^{\text{PGD}} \nonumber \\
B(\mathbf{u}_m, \mathbf{p} q^*(t)) &= L(\mathbf{p} q^*(t)), \quad \forall q^*(t) \in \mathcal{E}_T^{\text{PGD}} \tag{10}
\end{align}
$$

where $\mathcal{E}_S^{\text{PGD}}$ is the set of real-valued vector of size n, $\mathcal{E}_T^{\text{PGD}}$ is the set of periodic function of period $T = 2\pi/\omega$, $B(\mathbf{u}, \mathbf{v})$ and $L(\mathbf{v})$ are given in Eq.(9). We define the two applications $\mathcal{I}_{m-1}: \mathbf{q} \mapsto \mathbf{p}$ such that Eq.(10.1) holds (the space problem), and $\mathcal{I}_{m-1}: \mathbf{p} \mapsto \mathbf{q}$ such that Eq.(10.2) holds (the time problem). The couple $(\mathbf{p}_m, \mathbf{q}_m(t))$ is given when stationarity has been reached, i.e when $\mathbf{p}_m = \mathcal{I}_m(\mathbf{q}_m(t))$ and $\mathbf{q}_m(t) = \mathcal{I}_m(\mathbf{p}_m)$.

This procedure again results in Galerkin projections, but the projection is carried out on only one vector of the basis (either **p**<sub>m</sub> or **q**<sub>m</sub>(t)), leading to the following equation for the time (Eq.(11)) and the space problem (Eq.(12)):

$$
\begin{align}
& \mathbf{p}_m^T (\mathbf{M} \mathbf{p}_m \ddot{\mathbf{q}}_m + \mathbf{C} \mathbf{p}_m \dot{\mathbf{q}}_m + \mathbf{K} \mathbf{p}_m \mathbf{q}_m) + \mathbf{g} (\mathbf{P} \mathbf{q} + \mathbf{p}_m \mathbf{q}_m) \nonumber \\
& = \mathbf{p}_m^T \mathbf{f} - \mathbf{p}_m^T (\mathbf{M} \mathbf{P} \ddot{\mathbf{q}} + \mathbf{C} \mathbf{P} \dot{\mathbf{q}} + \mathbf{K} \mathbf{P} \mathbf{q}) \tag{11}
\end{align}
$$

$$
\begin{align}
& \int_I q_m \ddot{\ddot{q}}_m \mathbf{M} + \int_I q_m \dot{\dot{q}}_m \mathbf{C} + \int_I q_m q_m \mathbf{K} \mathbf{p}_m + \int_I q_m \mathbf{g} (\mathbf{P} \mathbf{q} + \mathbf{p}_m \mathbf{q}_m) \nonumber \\
& = \int_I q_m \mathbf{f} - \sum_{k=1}^{m-1} \left[ \int_I q_m \ddot{q}_k \mathbf{M} + \int_I q_m \dot{q}_k \mathbf{C} + \int_I q_m q_k \mathbf{K} \mathbf{p}_k \right] \tag{12}
\end{align}
$$

**Algorithm 2 ALGORITHM FOR PROGRESSIVE PGD**

**for** m = 1 **to** m<sub>max</sub> **do**

Initialize **q**<sub>m</sub>(t)

**for** k = 1 **to** k<sub>max</sub> **do**

Compute **p**<sub>m</sub> = $\mathcal{I}_m(\mathbf{q}_m)$ by solving Eq.(11)

Orthonormalize **p**<sub>m</sub> relative to **p**<sub>1</sub>, ..., **p**<sub>m-1</sub> [optional]

Compute $\mathbf{q}_m(t) = \mathcal{I}_m(\mathbf{p}_m)$ by solving Eq.(12)

**end for**

**if** update **then**

Update time functions **q**<sub>m</sub> = **T**<sub>m</sub>(**P**<sub>m</sub>) by solving Eq.(6)

Set **u**<sub>m</sub>(t) = **P**<sub>m</sub>**q**<sub>m</sub>(t) and check convergence

**else**

Set **u**<sub>m</sub> = **u**<sub>m-1</sub> + **p**<sub>m</sub>**q**<sub>m</sub>(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<sub>max</sub> times an algebraic problem with n unknowns and a differential equation with one variable, where k<sub>max</sub> 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<sub>max</sub> 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** after a new couple (p<sub>m</sub>, q<sub>m</sub>) has been computed. This update is done by solving the time problem associated with the mapping **T**: **P** ↦ **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<sub>m</sub>, q<sub>m</sub>) 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
<figure><img src="image_2.png" /><figcaption>**FIGURE 1.** SCHEME OF THE CLAMPED/CLAMPED BEAM AND REFERENCE ELEMENT</figcaption></figure>

<table><caption>**TABLE 1.** NUMERICAL PARAMETERS USED IN THE BEAM MODEL</caption><thead><tr><th>L (m)</th><th>h (m)</th><th>ES (Pa.m<sup>2</sup>)</th><th>EI (Pa.m<sup>4</sup>)</th><th>ρS (kg.m<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1</td><td>3 · 10<sup>-2</sup></td><td>1.89 · 10<sup>8</sup></td><td>1.41 · 10<sup>4</sup></td><td>7.02</td></tr></tbody></table>

state solutions. This means that the set of differential equations defining the time problem (Eqs.(6) or (12)) is transformed into a system of nonlinear algebraic equations having the form of Eq.(7). Since time problems only consist in a small set of differential equations of size $m$ ($m \le 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)$, $\dot{q}(t)$ and $\ddot{q}(t)$ are easily computed in the frequency domain. Indeed, the solution of time problems by HBM gives $q_i(t) = T(t)\tilde{q}_i$ where $T(t)$ contains the $2H + 1$ vectors of Fourier basis. Integrals can then be approximated by:

$$I_k = \int_I q_i(t) \frac{dq_j^k}{dt^k}(t) dt \approx \tilde{q}_i^T H (D^k \tilde{q}_j) \quad (13)$$

where $H = \int_I T^T(t)T(t)dt$ is a $2H + 1 \times 2H + 1$ matrix.

## 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 $\theta$) and $n_e = 20$ elements.

<figure><img src="image_3.png" /><figcaption>**FIGURE 2.** DISPLACEMENT OF THE CENTER OF THE BEAM AS A FUNCTION OF FREQUENCY FOR THE REFERENCE SOLUTION</figcaption></figure>

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$ arising from geometric nonlinearities given by Eq.(2).

## Forced Solution

In this section we consider that the beam is excited at its center by a purely transverse harmonic force $f(t) = A \cos(\omega t)$ with $A = 200N$. The reference solution $u_{ref}$ is computed by HBM with $H = 3$ harmonics for $\frac{\omega}{2\pi}$ from 155 to 180 Hz (upward direction) and for $\frac{\omega}{2\pi}$ 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_{POD}^m$ 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_m^{\omega_1}(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_m^{\omega_{i-1}}$. 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_{oPGD}^m$ (resp. $u_{pPGD}^m$). 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 achieved 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 $\epsilon = 0.005$.

Solution are compared relative to the reference solution by
using the relative difference $\epsilon_r = ||\mathbf{u}_{ref}(t) - \mathbf{u}(t)||/||\mathbf{u}_{ref}||$.

Figure 3 depicts the error relative to the reference solution for various number of retained mode $m \le 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 $\epsilon$. 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 \times m$ equations 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 $\epsilon = 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

<figure><img src="image_4.png" /><figcaption>FIGURE 3. RELATIVE ERROR $\epsilon_r$ FROM THE REFERENCE SOLUTION FOR POD, oPGD AND pPGD FOR $m = 1, 2, 3, 4$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4. DECOMPOSITION OF THE SOLUTION COMPUTED BY oPGD FOR $m = 4$ AT $\frac{\omega}{2\pi} = 180$ (shape have been normed for the purpose of comparison)</figcaption></figure>

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
<figure><img src="image_6.png" /><figcaption>**FIGURE 5.** DECOMPOSITION OF THE SOLUTION COMPUTED WITH pPGD FOR $m = 4$ AT $\frac{\omega}{2\pi} = 180$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**FIGURE 6.** DECOMPOSITION OF THE SOLUTION COMPUTED WITH oPOD FOR $m = 4$ AT $\frac{\omega}{2\pi} = 180$</figcaption></figure>

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 acceptable 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><img src="image_8.png" /><figcaption>**FIGURE 7.** BACKBONE CURVE OF THE FIRST NONLINEAR MODE COMPUTED BY HBM, POD and oPGD</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**FIGURE 8.** RELATIVE ERROR $\epsilon_r$ FROM THE REFERENCE SOLUTION FOR POD AND oPGD FOR $m = 1, 2$ (MNL computation)</figcaption></figure>

## ACKNOWLEDGMENT

Thanks go to SNECMA for its technical and financial support.

## REFERENCES

[1] Bai, Z., and Yangfeng, S., 2005. “Soar: A second order arnoldi method for the solution of the quadratic eigenvalue problem”. *Matrix Anal. Appl.*, **26**, pp. 640–659.

[2] Bai, Z., and Li, R., 2005. “Strcuture preserving model re-
duction using a krylov subspace projection formulation”. *Communication in mathematical science*, 3, pp. 179–199.

[3] Kerschen, G., Golinval, J., Vakakis, A., and Bergman, L., 2005. “The method of proper orthogonal decomposition for dynamical characterisation and order reduction of mechanical systems: an overview”. *Nonlinear dynamics*, **41**, pp. 147–169.

[4] Chang, Y., Lernout, E. C., Mignolet, M., and Soize, C., 2011. “Reduced order modelling for the nonlinear geometric response of some curved structures”.

[5] Ladevèze, P., Passieux, J., and Neron, D., 2010. “The latin multiscale computational method and the proper generalized decomposition”. *Comput. Methods Appl. Mech. Engrg*, **199**, pp. 1287–1296.

[6] Nouy, A., 2008. “Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms”. *Comput. Methods Appl. Mech. Engrg.*, **197**, pp. 4718–4736.

[7] Nouy, A., 2007. “A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations”. *Comput. Methods Appl. Mech. Engrg.*, **196**, pp. 4521–4537.

[8] Nouy, A., and LeMaitre, O., 2009. “Generalized spectral decomposition for stochastic nonlinear problems”. *Journal of computational physics*, **228**, pp. 202–235.

[9] Nouy, A., 2010. “A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations”. *Comput. Methods Appl. Mech. Engrg.*, **199**, pp. 1603–1626.

[10] K.Y. Sanliturk, D. E., 1996. “Modelling two-dimensional friction contact and its application using harmonic balance method”. *Journal of sound and vibration*, **193**(2), pp. 511–523.

[11] Groll, G., and Ewins, D., 2001. “The harmonic balance method with arc-length continuation in rotor stator contact problems”. *Journal of sound and vibration*, **241** (2), pp. 223–233.

[12] Grolet, A., and Thouverez, F., 2012. “Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model”. *Journal of sound and vibration*, **331**, pp. 2911–2928.

[13] Cochelin, B., and Vergez, C., 2009. “A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions”. *Journal of sound and vibrations*, pp. 243–262.

[14] Cameron, T., Griffin, J., and Grant, P., 1989. “An alternating frequency time domain method for calculating the steady state response of nonlinear dynamic systems”. *Journal of applied mechanics*, **56**(1), pp. 149–154.

[15] Sterian, D., and Mihai, L., 2010. “Low pressure centrifugal compressor”. *Wseas transaction on fluids mechanics*, **5**(3), pp. 226–235.

[16] Chinesta, F., Ammar, A., and E., E. C., 2009. “Proper generalized decomposition of multiscale models”. *Int. J Numer. Mech. Engrg.*

[17] Chinesta, F., Prulière, E., Ammar, A., and E., E. C., 2009. “Non incermental strategies for simulating thermomechanical models with uncertainty”. *Int. J Mater. Form.*, **2**, pp. 563–566.

[18] Chinesta, F., Ammar, A., and E., E. C., 2011. “Recent advance and new challenges in the use of the proper generalized decomposition for solving multi dimensional models”. *Arch. Comput. Methods Engs*, **17**, pp. 327–350.