# Extended-PGD Model Reduction for Nonlinear Solid Mechanics Problems Involving Many Parameters

Pierre Ladevèze, Charles Paillet, David Néron

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

## To cite this version:

Pierre Ladevèze, Charles Paillet, David Néron. Extended-PGD Model Reduction for Nonlinear Solid Mechanics Problems Involving Many Parameters. Advances in Computational Plasticity, pp.201-220, 2017, 10.1007/978-3-319-60885-3_10. hal-01926425

HAL Id: hal-01926425

https://hal.science/hal-01926425

Submitted on 19 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Extended-PGD Model Reduction for Nonlinear Solid Mechanics Problems Involving Many Parameters

P. Ladevèze, Ch. Paillet and D. Néron

**Abstract** Reduced models and especially those based on Proper Generalized Decomposition (PGD) are decision-making tools which are about to revolutionize many domains. Unfortunately, their calculation remains problematic for problems involving many parameters, for which one can invoke the “curse of dimensionality”. The paper starts with the state-of-the-art for nonlinear problems involving stochastic parameters. Then, an answer to the challenge of many parameters is given in solid mechanics with the so-called “parameter-multiscale PGD”, which is based on the Saint-Venant principle.

## 1 Introduction

Numerical simulation has made a forceful entry into design and analysis offices. This revolution, which is anything but complete, has entered in a new stage, called simulation-driven “robust” design. It leads to a major scientific challenge: simulations should be performed in quasi real-time. The key is a new generation of reduced-order methods which comprises essentially the Proper Orthogonal Decomposition (POD), the Proper Generalized Decomposition (PGD) and the Reduced Basis Method (RB), the basis and recent developments of which are given in [6]. These are problems in which uncertainties or variations in parameters are to be taken into account, or problems with very high number of degrees of freedom, with multiple scales or interactions between several physics. These methods, together with the notions of offline and online calculations, also open the way to new approaches where simulation and analysis of structures can be carried out almost in real time.

The object of this work is the PGD, that was introduced in [13, 14] for the treatment of nonlinear time-dependent problems of solid mechanics. Many developments have been made over the last thirty years: multi-scale, multiphysics, stochastic or non-stochastic parameters, acoustics, large displacements and deformations ...

P. Ladevèze (✉) · Ch. Paillet · D. Néron

LMT (ENS Paris-Saclay, CNRS, Université Paris-Saclay), 61 Avenue du Président

Wilson, 94235 Cachan Cedex, France

e-mail: ladeveze@lmt.ens-cachan.fr
In [18], one can find a synthesis of most of the developments in Cachan, where the method LATIN plays a central role. Currently, a number of tools have become mature and have been applied to industrial cases and then are competitors of classical computational methods (see book [7]). This is the case for the calculation of (visco)-plastic structures with less than ten parameters. The PGD not only makes it possible to construct reduced models that can be used in real time, but it also makes it possible to reduce drastically the calculation time in many situations. However, a major limitation is still the number of the parameters that can be involved (no more than 10). The paper starts with the state-of-the-art for nonlinear problems involving parameters which could be stochastic. It appears that the key to extend the PGD to a large number of parameters and then to cancel the “curse of dimensionality” is the same for linear and nonlinear problems.

Several attempts have been introduced for problems with a large number of parameters: Tucker tensors, Tensor Train tensors, Hierarchical Tucker tensors or more general tree-based Hierarchical Tucker tensors [3, 8, 9, 12, 24]. As they are generic a priori approximations, we believe that their improvement which is real, is limited. Here, in this paper we are following quite a different way. The idea is to notice that the operator underlying the problem treated is a binder between the parameters. Thus, at the level of the solution, the “curse of dimensionality” is no longer relevant. More, we go further with the parameter-multiscale PGD which takes its source in a fine analysis of the solution according to parameters, analysis built on the Principle of Saint-Venant. The basic idea is to introduce a two-level description of each parameter, the “macro” scale and the “micro” scale as one considers two scales for the space or the time. Also, to implement this vision, we have been led to use completely discontinuous approximations [18]. To carry out these idea, we use the Weak-Trefftz Discontinuous Method introduced in [16] and applied in [20] for the calculation of “medium frequency” phenomena. This paper is limited to fundamental aspects and the first numerical experiments for linear problems. More details can be founded in [25].

## 2 Nonlinear Problems with Parameters: The State-of-the-art

We have been working on ROM computation for 30 years with the so-called LATIN-PGD and what we are doing at the present time is the result of many works. PGD means Proper Generalized Decomposition and LATIN denotes the computational method which is nonincremental. The LATIN-PGD method was introduced in [13] for viscoplastic materials whose constitutive relations are described using a functional approach. Its extension to modern material descriptions involving internal variables, still for viscoplastic materials, was proposed in [14]. A number of mathematical properties regarding convergence and error indicators were proved in the book [15]. Overview could be found in [15, 18]. Originally, PGD was called radial
loading approximation, which, to us, meant a “mechanics” approximation in solid mechanics. LATIN-PGD leads today to a general and robust PGD computation technique which is based on an “abstract” reformulation of parametric nonlinear solid mechanics problems defined over a time-space-parameter domain. This work is the follow-up to [11, 19, 22, 27]. Today, there are few other works except the works of Ryckelynck and its group done with POD [28]. However, things are changing and today there are more and more POD-approaches developed in relation with the homogenization technique FE2 [10, 21, 26]. Additional reduction or interpolation are then introduced to reduce their computation cost. This is done offline with the PGD computation. Our answer is the so-called Reference Point Method (RPM) [4].

## 2.1 The Reference Problem and Its Reformulation

### 2.1.1 Notations

With the assumption of small perturbations, let us consider the quasi-static and isothermal evolution of a structure defined over the time-space domain $[0, T] \times \Omega$ and depending of the parameters $\mu \in \Sigma_{\mu}$. This structure is subjected to prescribed body forces $f_{\underline{d}}$, traction forces $F_{\underline{d}}$ over a part $\partial_2 \Omega$ of the boundary, and displacements $u_{\underline{d}}$ over the complementary part $\partial_1 \Omega$ (see Fig. 1). All the data and the material characteristics depend on the parameters $\mu \in \Sigma_{\mu}$.

The state of the structure is defined by the set of the fields $s = (\dot{\epsilon}_p, \dot{X}, \sigma, Y)$ (where the dot notation $\dot{\phantom{X}}$ denotes the time derivative), in which:

* $\epsilon_p$ designates the inelastic part of the strain field $\epsilon$ which corresponds to the displacement field $\underline{u}$, uncoupled into an elastic part $\epsilon_e$ such that $\epsilon_p = \epsilon - \epsilon_e$; $X$ designates the remaining internal variables;

* $\sigma$ designates the Cauchy stress field and $Y$ the set of variables conjugate of $X$ ($Y$ and $X$ have the same dimension). $X$ could be hardening variables, damage variables, chemical variables, ...

<figure><img src="image_3.png" /><figcaption>**Fig. 1** The reference problem</figcaption></figure>
• generalized quantities over the parameter-time-space domain.

All these quantities are defined over the parameter-time-space domain $\Sigma_{\mu} \times [0, T] \times \Omega$ and assumed to be sufficiently regular. For the sake of simplicity, the displacement $\underline{u}$ alone is assumed to have a nonzero initial value, denoted $\underline{u}_0$. Introducing the following notations for the primal fields:

$$ \mathbf{e}_p = \begin{bmatrix} \epsilon_p \\ -\mathbf{X} \end{bmatrix}, \quad \mathbf{e} = \begin{bmatrix} \epsilon \\ 0 \end{bmatrix} \quad \text{and} \quad \mathbf{e}_c = \begin{bmatrix} \epsilon_c \\ \mathbf{X} \end{bmatrix} \quad \text{so that} \quad \mathbf{e}_p = \mathbf{e} - \mathbf{e}_c \quad (1) $$

and for the dual fields:

$$ \mathbf{f} = \begin{bmatrix} \sigma \\ Y \end{bmatrix} \quad (2) $$

The mechanical dissipation rate for the entire structure $\Omega$ is:

$$ \int_{\Omega} (\dot{\mathbf{e}}_p : \boldsymbol{\sigma} - \dot{\mathbf{X}} \cdot \mathbf{Y}) d\Omega = \int_{\Omega} (\dot{\mathbf{e}}_p \circ \mathbf{f}) d\Omega \quad (3) $$

where $\cdot$ denotes the contraction adapted to the tensorial nature of $X$ and $Y$. Notation $\circ$ denotes the contraction operator for generalized quantities. Let us introduce the following fundamental bilinear “dissipation” form:

$$ \langle \mathbf{s}, \mathbf{s}' \rangle = \int_{\Sigma_{\mu} \times [0, T] \times \Omega} \left(1 - \frac{t}{T}\right) (\dot{\mathbf{e}}_p \circ \mathbf{f}' + \dot{\mathbf{e}}_p' \circ \mathbf{f}) d\Omega dt d\mu \quad (4) $$

along with $\mathbf{E}$ and $\mathbf{F}$, the spaces of the fields $\dot{\mathbf{e}}_p$ and $\mathbf{f}$ which are compatible with (4). These spaces enable us to define $\mathbf{S} = \mathbf{E} \times \mathbf{F}$, the space in which the state $\mathbf{s} = (\dot{\mathbf{e}}_p, \mathbf{f})$ of the structure is being sought.

### 2.1.2 The State Equations over $\Sigma_{\mu} \times [0, T] \times \Omega$

Following [15], a normal formulation with internal state variables is used to represent the behavior of the material. If $\rho$ denotes the mass density of the material, from the free energy $\rho\Psi(\epsilon_c, X)$ with the usual uncoupling assumptions, the state law yields:

$$ \sigma = \rho \frac{\partial \psi}{\partial \epsilon_c} = \mathbf{K} \epsilon_c \quad (5) $$

$$ Y = \rho \frac{\partial \psi}{\partial X} = \Lambda X $$

where the Hooke's tensor $\mathbf{K}$ and the constant, symmetric and positive definite tensor $\Lambda$ are material characteristics.
### 2.1.3 The State Evolution Laws over $\Sigma_{\mu} \times [0, T] \times \Omega$

The state evolution laws can be written:

$$ \dot{\mathbf{e}}_{\mathbf{p}} = \mathbf{B}(\mathbf{f}) \quad \text{with} \quad \mathbf{e}_{\mathbf{p}}|_{t=0} = 0 \qquad (6) $$

where $\mathbf{B}$ is a positive operator which is also for most viscoplastic models maximal monotone [15]. Let us introduce now the space $\mathcal{W}_{\mu,ad}^{[0,T]}$ of admissible displacement fields $\underline{u}$ defined over $\Sigma_{\mu} \times [0, T] \times \Omega$ and $\mathcal{W}_{\mu,0}^{[0,T]}$ the associated vectorial space. The compatibility equation can be written as:

$$ \begin{aligned} \text{Find } \underline{u} \in \mathcal{W}_{\mu,ad}^{[0,T]} \text{ such that } & \forall \underline{\mu} \in \Sigma_{\mu}, \forall \underline{u}^{*} \in \mathcal{W}_{\mu,0}^{[0,T]}, \\ & \int_{[0,T] \times \Omega} \text{Tr}[\underline{\epsilon}(\underline{u})\mathbf{K}\underline{\epsilon}(\underline{u}^{*})]d\Omega dt = \int_{[0,T] \times \Omega} \text{Tr}[\underline{\epsilon}_{\mathbf{p}}\mathbf{K}\underline{\epsilon}(\underline{u}^{*})]d\Omega dt + \\ & \int_{[0,T] \times \Omega} f_{-d} \cdot \underline{u}^{*} d\Omega dt + \int_{[0,T] \times \partial_{2}\Omega} F_{-d} \cdot \underline{u}^{*} dS dt \quad (7) \end{aligned} $$

It follows that the stress $\sigma = \mathbf{K}(\underline{\epsilon}(\underline{u}) - \underline{\epsilon}_{\mathbf{p}})$ can be written:

$$ \sigma = \Omega \underline{\epsilon}_{\mathbf{p}} + \mathbf{r}_{d} \qquad (8) $$

where $\Omega$ is a linear given operator and $\mathbf{r}_d$ is a prestress depending on the data. Introducing the generalized stress, the admissibility conditions can be written as:

$$ \mathbf{f} = \mathbf{Q} \mathbf{e}_{\mathbf{p}} + \mathbf{r}_{d} \qquad (9) $$

with

$$ \mathbf{Q} = \begin{bmatrix} \mathbf{\Omega} & \mathbf{0} \\ \mathbf{0} & \mathbf{\Lambda} \end{bmatrix} \quad \text{and} \quad \mathbf{r}_d = \begin{bmatrix} \mathbf{r}_d \\ \mathbf{0} \end{bmatrix} \qquad (10) $$

where $\mathbf{Q}$ is a linear symmetric positive operator. Finally, the problem to solve, which is defined over $\Sigma_{\mu} \times [0, T] \times \Omega$, is:

$$ \begin{gathered} \text{Find } \mathbf{s} = (\dot{\mathbf{e}}_{\mathbf{p}}, \mathbf{f}) \in \mathbf{S} \text{ such that:} \\ \mathbf{f} = \mathbf{Q}\mathbf{e}_{\mathbf{p}} + \mathbf{r}_{d} \quad \text{and} \quad \dot{\mathbf{e}}_{\mathbf{p}} = \mathbf{B}(\mathbf{f}) \quad \text{with} \quad \mathbf{e}_{\mathbf{p}}|_{t=0} = 0 \end{gathered} \qquad (11) $$

Consequently, one has to solve a first order differential equation with an initial condition; The operators $\mathbf{Q}$ and $\mathbf{B}$ as well as the right-hand-side member $\mathbf{r}_d$ depends on the parameter $\underline{\mu}$ belonging to the parameter set $\Sigma_{\mu}$.
## 2.2 The LATIN Solver for PGD Computation

### 2.2.1 The LATIN Solver

The LATIN method is an iterative strategy which differs from classical incremental or step-by-step techniques in that, at each iteration, it produces an approximation of the complete structural response over the whole loading history being considered. A review of the state-of-the-art and more recent extensions could be found in [18]. The LATIN method is designed as a mechanics-based computational strategy whose aim is to achieve the best possible performance level for solid mechanics problems. Consequently, this alternative approach is rooted in some remarkable properties which are verified by most of the models encountered in structural mechanics.

The LATIN method operates here over the parameter-time-space domain $\Sigma_{\mu} \times [0, T] \times \Omega$, and its first principle (P1) consists in separating the difficulties. Thus, the equations are divided into:

* a set of linear equations which can be global in the space variables: the equilibrium and compatibility equations, the state equations;

* a set of equations which are local in the space variables but can be nonlinear: the state evolution laws.

The reformulation (11) of the reference problem enters into this framework because:

* Q is a linear operator;

* B is at least local in space variables.

In the geometric representation given Fig. 2, $A_d$ and $\Gamma$ represent the solutions of the first and second set respectively, $A_d$ being associated to the free energy and $\Gamma$ with the dissipation. The LATIN second principle (P2) is also very natural. It consists in solving the two sets of equations alternatively until practical convergence. In order to do that, one uses search directions given as parameters of the LATIN method. One possible choice (Newton search direction) consists of the tangent direction and its conjugate direction (see Fig. 3).

**Local stage principle at iteration n + 1**—Find $\hat{s}_{n+1/2} = (\hat{e}_{p,n+1/2}, \hat{f}_{n+1/2}) \in S$ such that:

<figure><img src="image_4.png" /><figcaption>**Fig. 2** The geometric representation associated to the reformulation of the reference problem</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Iteration n + 1 of the LATIN method over $[0, T] \times \Omega$</figcaption></figure>

$$ 
\begin{aligned}
\hat{\mathbf{e}}_{p,n+1/2} &= \mathbf{B}(\hat{\mathbf{f}}_{n+1/2}) \quad \text{with} \quad \hat{\mathbf{e}}_{p,n+1/2} = 0 \text{ at } t=0 \\
\hat{\mathbf{e}}_{p,n+1/2} - \dot{\mathbf{e}}_{p,n} + \mathbf{H}^{+}(\hat{\mathbf{f}}_{n+1/2} - \mathbf{f}_{n}) &= 0
\end{aligned}
\qquad (12) $$

The search direction $\mathbf{H}^+$ is a LATIN parameter. Practically, one takes a linear positive operator which is local both on time and space variables. This local stage is very suitable for parallel computing.

**Linear stage principle at iteration n + 1**—Find $\mathbf{s}_{n+1} = (\dot{\mathbf{e}}_{p,n+1}, \mathbf{f}_{n+1}) \in \mathbf{S}$ such that:

$$ 
\begin{aligned}
\mathbf{f}_{n+1} &= \mathbf{Q}\mathbf{e}_{p,n+1} + \mathbf{r}_d \\
\dot{\mathbf{e}}_{p,n+1} - \dot{\mathbf{e}}_{p,n+1/2} - \mathbf{H}^{-}(\mathbf{f}_{n+1} - \hat{\mathbf{f}}_{n+1/2}) &= 0 \quad \text{with} \quad \mathbf{e}_{p,n+1} = 0 \text{ at } t=0
\end{aligned}
\qquad (13) $$

The search direction $\mathbf{H}^-$ is a LATIN parameter. This is a linear positive operator which is local both on time and space variables. It is associated to the material operator **B**. One has to solve a first order linear differential equation with an initial condition, the operator **Q** being non-explicit.

**Remark** In practice, $\mathbf{H}^-$ is chosen close to the tangent to the manifold $\mathbf{\Gamma}$ at the point $\hat{\mathbf{s}}_{n+1/2} = (\hat{\mathbf{e}}_{p,n+1/2}, \hat{\mathbf{f}}_{n+1/2})$. For $\mathbf{H}^+$, one takes $\mathbf{0}$ or $\mathbf{H}^-$. The convergence of the iterative process has been proved in the case of non-softening materials and contacts without friction [15]. The distance between two successive approximations gives a good and easily computed error indicator. Let us also note that one often uses an additional relaxation with a coefficient equal to 0.8.

### 2.2.2 PGD Computation

Let us recall that the problem is defined over the parameter-time-space domain $\Sigma_{\mu} \times [0, T] \times \Omega$.

**Linear stage at iteration n + 1**—Let us introduce corrections:

$$ 
\begin{aligned}
\Delta \dot{\mathbf{e}}_{\mathrm{p}} &= \dot{\mathbf{e}}_{\mathrm{p},n+1} - \dot{\mathbf{e}}_{\mathrm{p},n} \\
\Delta \mathbf{f} &= \mathbf{f}_{n+1} - \mathbf{f}_{n}
\end{aligned}
\qquad (14) $$
where $\mathbf{s}_{n+1} = (\dot{\mathbf{e}}_{p,n}, \mathbf{f}_n)$ has been computed at iteration $n$. The problem to solve over $\Sigma_{\mu} \times [0, T] \times \Omega$ at iteration $n+1$ is then:

Find $\Delta\mathbf{s} = (\Delta\dot{\mathbf{e}}_p, \Delta\mathbf{f}) \in \mathbf{S}$ such that

$$ \Delta\mathbf{f} = \mathbf{Q}\Delta\dot{\mathbf{e}}_p \quad (15) $$

$$ \Delta\dot{\mathbf{e}}_p - \mathbf{H}^{-} \Delta\mathbf{f} = \mathbf{R}_d \quad \text{with} \quad \Delta\mathbf{e}_p = 0 \text{ at } t=0 $$

The main idea is to interpret the search direction as a linear constitutive relation, the operator $\mathbf{H}^{-}$ being local both on time and space variables and positive definite as the Hooke tensor. Consequently, one introduces the associated constitutive relation error which will be minimized:

$$ r(\Delta\mathbf{s}, t) = \frac{1}{2} \int_{\Omega} [\Delta\dot{\mathbf{e}}_p - \mathbf{H}^{-} \Delta\mathbf{f} - \mathbf{R}_d](\mathbf{H}^{-})^{-1} [\Delta\dot{\mathbf{e}}_p - \mathbf{H}^{-} \Delta\mathbf{f} - \mathbf{R}_d] d\Omega \quad (16) $$

and

$$ R(\Delta\mathbf{s}) = \int_{\Sigma_{\mu} \times [0, T]} \left(1 - \frac{t}{T}\right) r(\Delta\mathbf{s}, t) dt d\mu \quad (17) $$

with $\Delta\mathbf{s} = (\Delta\dot{\mathbf{e}}_p, \Delta\mathbf{f}) \in \mathbf{S}$. The problem (15) becomes:

Find $\Delta\mathbf{s} \in \mathbf{S}$ minimizing

$$ \Delta\mathbf{s} \in \mathbf{S} \mapsto R(\mathbf{s}) \in \mathbb{R} \quad (18) $$

with the constraints $\Delta\mathbf{f} = \mathbf{Q}\Delta\dot{\mathbf{e}}_p$ and $\Delta\mathbf{e}_p = 0$ at $t=0$

One only prescribes that:

$$ \Delta\mathbf{e}_p = \sum_{i=1}^{m} \lambda_i(t) \gamma_i(\underline{\mu}) \mathbf{g}_i(\underline{x}) \quad (19) $$

with $\lambda_i(0) = 0$ (initial condition), $\mathbf{g}_i \in L^2(\Omega)$, $\lambda_i(t) \in L^2[0, T]$ and $\gamma_i(\underline{\mu}) \in L^2(\Sigma_{\mu})$. It follows, using admissibility conditions, that:

$$ \Delta\mathbf{f} = \sum_{i=1}^{m} \lambda_i(t) \gamma_i(\underline{\mu}) \mathbf{Q} \mathbf{g}_i(\underline{x}) \quad (20) $$

where $\mathbf{Q}\mathbf{g}_i(\underline{x})$ are computed solving several elasticity problems. For the sake of simplicity, let us consider now that $\mathbf{Q}$ does not depend on $t$ and $\underline{\mu}$ belonging to $\Omega \times [0, T] \times \Sigma_{\mu}$. The general case does not involve serious difficulties, the $\mathbf{Q}$-constraint being satisfied in a mean sense.

A greedy algorithm with updating is used to solve the minimization problem. More details can be founded in [17, 18]. For few parameters, it could be advantageous to describe point-by-point the parameter space using the remarkable property
of the LATIN method: the initialization of the iterative process can be any function
defined over $[0, T] \times \Sigma_{\mu}$ [27].

**Local stage at iteration n + 1**—One has the solve a problem which is local over the parameter-time-space domain for which “hyper-reduction” techniques are also welcome. A very popular technique is the Empirical Interpolation Method (EIM) [2] and its discrete version named DEIM [5]. The Hyperreduction method [28] makes the most of a restricted subdomain of the space domain. Here we use the Reference Point Method (RPM) developed in [4]. Let us note that such method can be also used for solving the linear stage.

## 2.3 An Engineering Illustration

To illustrate the use of the technique to deal with parametrized problems, we consider an example issued from [22] and which is freely inspired from a blade of the Vulcain engine of the Ariane 5 launcher. The geometry, boundary conditions and mesh are presented on Fig. 4. A four-sinusoidal-cycles displacement with is prescribed on the lower part. The total number of DOFs is 141,500 and the time interval is discretized using 120 time steps. The material coefficients used for the Marquis-Chaboche elastic-viscoplastic material are typical of a Titanium TA6V material at 500° K. The parametric study is defined by Table 1. It concerns the influence of the loading amplitude, of the limit stress and the power coefficient in the evolution law. The range of variation of each parameter was discretized into 10 values, leading to 1000 different problems. Figure 5 shows the total agreement of the results obtained with the LATIN-PGD compared to the one obtained with ABAQUS.

<figure><img src="image_6.png" /><figcaption>Fig. 4 Geometry, boundary conditions and mesh of the blade test-case</figcaption></figure>
<table><caption>**Table 1** Range of variation of loading amplitude, R₀ and γ</caption><thead><tr><td>Parameter</td><td>Min. value</td><td>Max. value</td><td>Step</td><td>Range of variation (%)</td></tr></thead><tbody><tr><td>±loading amplitude</td><td>0.1 mm</td><td>0.19 mm</td><td>0.01 mm</td><td>±31</td></tr><tr><td>Limit stress R₀</td><td>20 MPa</td><td>29 MPa</td><td>1 MPa</td><td>±18</td></tr><tr><td>Power γ</td><td>285</td><td>330</td><td>5</td><td>±7</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 5** Stress versus strain curves at the most loaded gauss point</figcaption></figure>

To illustrate the performances of LATIN-PGD, we compare the CPU times:

* about 50 days (estimated time) are necessary to complete the 1000 resolution with ABAQUS;

* less than 17 h are necessary to complete the 1000 resolution with the multiple runs algorithm.

The gain is about 70 using the multiple runs algorithm, but can achieve more than 700 when using also the RPM strategy. The important point is that, once this parametric study has been performed, a reduced-model of the nonlinear problem is built. For example, some stochastic studies can be performed very easily. Assuming a probabilistic distribution of the three parameters, a probabilistic distribution of the quantity of interest can be computed in quasi-real time. Let us consider now that the first parameter is fuzzy and the two others follow a normal law. Then, it is easy to compute the probability law related to the maximum Mises stress over [0, T] × Σμ (see Fig. 6).
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Probability of the maximum Mises stress</figcaption></figure>

# 3 The Parameter-Multiscale PGD

## 3.1 Model Problem

The basic principles of parameters-multiscale-PGD are given in [17]. One considers an elastic media which occupies the domain $\Omega$ divided into sub-domains or elements $\Omega_E, E \in E$. The parameters $\mu_E \in \mathbb{R}^{q_E}$ are associated to the rigidity of the volume $\Omega_E$. For the sake of simplicity, we consider that $\mu_E$ is a scalar; it belongs to $[-1/2; 1/2]$. Let us introduce $\mu = \{\mu_E \mid E \in E\}$, the corresponding space being $\Sigma_\mu$. The problem to solve can be written as:

Find $\underline{X}(\mu) \in U^\mu$ where $U = \mathbb{R}^N$ such that:

$$ \forall \underline{\mu} \in \underline{\Sigma}_{\mu} \qquad A \underline{X}(\underline{\mu}) = \underline{F}_{d} $$

where $A$ is a linear positive definite operator depending on $\underline{\mu}$. $\underline{F}_d$ is a given loading which could also depend on $\underline{\mu}$.

**Numerical illustration with the standard PGD**—The standard parameter-PGD has been introduced by [1]. Stochastic framework has been considered by [23]. An overview is given in [6]. Figure 7 shows an illustration of the model problem. This is a cube submitted to a uniaxial traction displacement, the opposite face being clamped. The parameters $\underline{\mu}$ could be interpreted as damage intensity.

Convergence curves for the classical PGD are given Fig. 8 for 8, 27 and 64 parameters; they show that convergence cannot be obtained when the number of parameters is more than 30.

## 3.2 The Key Idea

One has to solve:

$$ A(\underline{\mu})\underline{X}(\underline{\mu}) = \underline{F}_d $$
<figure><img src="image_9.png" /><figcaption>Fig. 7 Model problem with 27 parameters: particular solution and highlighting of a subdomain</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8 Convergence curve of the classical PGD</figcaption></figure>

Let us introduce $\mathbf{A}^0 = \mathbf{A}(\underline{0})$ and let us suppose for the sake of simplicity that $\underline{F}_d$ does not depend on $\underline{\mu}$. One has:

$$ \mathbf{A} = \mathbf{A}^0 \left[ 1 - \underbrace{(1 - \mathbf{A}^{0-1} \mathbf{A})}_{\Delta} \right] $$

Let be $\underline{X}_0 = \mathbf{A}^{0-1} \underline{F}_d$; the solution can be written as:

$$ \underline{X}(\underline{\mu}) : \underline{X}_0 + \underbrace{\Delta \underline{X}_0}_{\underline{X}_1(\underline{\mu})} + \underbrace{\Delta^2 \underline{X}_0}_{\underline{X}_2(\underline{\mu})} + \dots \qquad (21) $$
where $\underline{X}_1$ is $\mu$-linear and $\underline{X}_2$ $\mu$-quadratic. It follows that with only few terms, one gets a good approximation and then, the curse of dimensionality disappears; the operator $\mathbf{A}$ links the parameters. Let us go further: on each element, one has:

$$ \mathbf{A}_E = \mathbf{A}_E^0 + \mu_E \bar{\mathbf{A}}_E $$

and then

$$ \mathbf{A} = \bigwedge_{E \in E} \mathbf{A}_E $$

where $\mathbf{A}_{E \in E}$ is the finite element assembly operator. Introducing the operator $\mathbf{I}_E$, giving the restriction $\underline{V}_E$ of a spatial vector $\underline{V}$ on the subdomain $\Omega_E$ through the relation $\underline{V} = \mathbf{I}_E \underline{V}_E$, one can write: $\mathbf{A} = \mathbf{A}_{E \in E} \bar{\mathbf{A}}_E = \sum_{E \in E} \mathbf{I}_E \mathbf{A}_E \mathbf{I}_E^T$. It follows:

$$ \underline{X}_1(\underline{\mu}) = -\mathbf{A}^{-1} \left[ \sum_{E \in E} \mu_E \mathbf{I}_E \bar{\mathbf{A}}_E \mathbf{I}_E \underline{X}_0 \right] = \sum_{E \in E} -\mu_E \mathbf{A}^{-1} \mathbf{I}_E \underbrace{\left[ \bar{\mathbf{A}}_E \underline{X}_{0,E} \right]}_{\underline{Z}_E} $$

The term $\mathbf{I}_E \underline{Z}_E$ is associated to a self-equilibrated loading. Therefore, from the Principle of Saint-Venant, the solution is localized in the neighborhood of the element $E$, essentially over the elements sharing a common point with $E$ noted $C_E$. Let be:

$$ \underline{Z}_{1,E} = \mathbf{A}^{-1} \mathbf{I}_E \underline{Z}_E $$

which can be seen as negligible over the complement of $C_E$. One has:

$$ \underline{Z}_{1,E} \equiv \mathbf{A}^{-1} \mathbf{I}_E \underline{Z}_E \approx \mathbf{I}_{E'} \underline{Z}_{1,E'E} $$

and $\underline{X}_1(\underline{\mu}) = - \sum_{E \in E} \mu_E \underline{Z}_{1,E}$; $\underline{X}_1$ is then $\mu$-linear and $\underline{Z}_{1,E}$ is localized over the neighborhood $C_E$ of $E$. A similar property holds for $\underline{X}_2$:

$$ \begin{align*} \underline{X}_2(\underline{\mu}) &= \sum_{E' \in E} \sum_{E \in E} \mu_{E'} \mu_E \mathbf{A}^{-1} \mathbf{I}_{E'} \bar{\mathbf{A}}_{E'} \mathbf{I}_{E'} \underline{Z}_{1,E} \\ &\approx \sum_{E' \in E} \sum_{E \in E} \mu_{E'} \mu_E \mathbf{A}^{-1} \mathbf{I}_{E'} \bar{\mathbf{A}}_{E'} \mathbf{I}_{E'} \mathbf{I}_{E''} \underline{Z}_{1,E''E} \\ &\approx \sum_{E' \in C_E} \sum_{E \in E} \mu_{E'} \mu_E \mathbf{A}^{-1} \mathbf{I}_{E'} \bar{\mathbf{A}}_{E'} \underline{Z}_{1,E'E} \end{align*} $$
$$ \approx \sum_{E \in \mathcal{E}} \mu_E^2 \mathbf{A}^{0-1} \mathbf{I}_E \bar{\mathbf{A}}_E Z_{1,EE} + \sum_{E \in \mathcal{E}} \sum_{E' \in \bar{\mathcal{C}}_E} \mu_E \mu_{E'} \mathbf{A}^{0-1} \mathbf{I}_{E'} \bar{\mathbf{A}}_{E'} Z_{1,E'E} $$

$$ \approx \sum_{E \in \mathcal{E}} \mu_E^2 Z_{2,E} + \sum_{E \in \mathcal{E}} \sum_{E' \neq E} \mu_E \mu_{E'} Z_{2,EE'} $$

It follows from the Saint-Venant principle that the quadratic term of $X_2(\mu)$ is such that $Z_{2,EE'}$ is in practice localized over a close neighborhood of $C_E$. The second term is linear with respect to each parameter; more, the $\mu_E$-contribution concerns essentially the neighborhood of $C_E$.

The parameter-multiscale PGD that we propose is based on these remarks. Thus we introduce two scales, micro and macro, to describe the parameter space $\Sigma_\mu$. In this way, we propose as approximation:

$$ \underline{X}_E(\underline{\mu}) = \sum_{i=1}^{m} \tilde{X}_E^{(i)} \prod_{E'' \notin \bar{\mathcal{C}}_E} f_{E''}^{M(i)}(\mu_{E''}) \prod_{E' \in \bar{\mathcal{C}}_E} g_{EE'}^{m(i)}(\mu_{E'}) \quad (22) $$

where $f^M$ and $g^m$ are respectively “macro” and “micro” functions. $\bar{C}_E$ denotes a chosen neighborhood of $E$ defining the “micro” impact in space of the parameter $\mu_E$. One can take the elements having a common point with $E$. Here, we will choose a linear discretization of the “macro” functions $f_{E''}^M$ and thus consider only their value on two points, $\mu_{E''} = \{\pm 1/2\}$.

However, there is a difficulty: $X(\mu)$ is discontinuous from an element to another. It is removed using the so-called *Weak-Trefftz Discontinuous Galerkin method* (WTDG) proposed in [19] and extended to the quasistatic loadings in [18]. This approach is introduced in the next paragraph.

## 3.3 The WTDG Method

The domain is still split into elements or subdomains $\Omega_E, E \in \mathcal{E}$ on which the Hooke tensor $K_E$ is constant but depends on $\mu_E$, belonging to $[-1/2, 1/2]$. The classical WTDG is modified here by adding a regularization term which assure the positivity of the operator, even if the rigidity vanishes over one or several subdomains.

The admissible space associated to $E \in \mathcal{E}$ is denoted: $\mathcal{W}_{E,ad}^{h,\mu}$ and the associated vectorial space: $\mathcal{W}_{E,0}^{h,\mu}$.

The problem to solve is then:
Find $\underline{U}_E \in \mathcal{W}_{E,ad}^{h,\mu}$, $E \in \mathbf{E}$ such that $\forall \underline{U}_E^* \in \mathcal{W}_{E,0}^{h,\mu}$ one has:

$$ \begin{aligned} \forall \underline{\mu} \in \Sigma_{\underline{\mu}} & \sum_{E \in \mathcal{E}} \left[ \sum_{E' \in \bar{C}_E} \left[ \int_{\Gamma_{EE'}} dS \left[ (\underline{\sigma n} + \underline{\sigma'} n') \cdot \frac{(\underline{U}^* + \underline{U}^{*'})}{2} + (\underline{U} - \underline{U}') \cdot \frac{(\underline{\sigma}^* n - \underline{\sigma'}^* n')}{2} \right. \right. \\ & \qquad \left. \left. + k(\underline{U} - \underline{U}') \cdot (\underline{U}^* - \underline{U}^{*'}) \right] \right] - \frac{1}{2} \int_{\Omega_E} \left[ (\underline{div}(\underline{\sigma}) + \underline{f}_d) \cdot \underline{U}^* + \underline{div}(\underline{\sigma}^*) \cdot \underline{U} \right] d\Omega \\ & \qquad + \sum_{\Gamma_{EE} \subset \partial_1 \Omega} \int_{\Gamma_{EE}} \left[ (\underline{U} - \underline{U}_d) \cdot \underline{\sigma}^* n + k(\underline{U} - \underline{U}_d) \cdot \underline{U}^* \right] dS \\ & \qquad + \sum_{\Gamma_{EE} \subset \partial_2 \Omega} \int_{\Gamma_{EE}} \left[ (\underline{\sigma n} - \underline{F}_d) \cdot \underline{U}^* \right] dS \right] = 0 \end{aligned} \quad (23) $$

with $\underline{\sigma} = \underline{K\epsilon}(\underline{U})$. $\Gamma_{EE'}$ is the common boundary of two adjacent subdomain $E$ and $E'$; $\Gamma_{EE}$ is the common boundary of $\Omega_E$ and $\partial\Omega$. One consider here approximations such that the interior equation

$$ \underline{div}(\underline{\sigma}) + \underline{f}_d = 0 $$

is satisfied in an average sense, i.e. in resultant and moment over $\Omega_E$. Two elements WP1 and WP2 can be easily associated to the classical elements P1 and P2. With these elements, one gets a coercivity property leading to the unicity of the solution [17].

On each element $E$, the WTDG leads to a contribution:

$$ \sum_{E' \in \bar{C}^E} \underline{A}_{EE'} \underline{X}^{E'} $$

The symmetric part is $\underline{A}_E$ and its value for $\underline{\mu} = \underline{0}$ is $\underline{A}_E^0$. The generalized force given by the WTDG is then:

$$ \underline{A}_{E \in \mathcal{E}} \sum_{E' \in \bar{C}^E} \underline{A}_{EE'} \underline{X}^{E'} $$

## 3.4 Computational Method

### 3.4.1 Basic Operators

* Computation of $\underline{I}_E(\underline{R}^E) \equiv \underline{\tilde{Z}}^E \in \mathbb{R}^n$

$\underline{R}$ is a known residual, the contribution to the subdomains $E \in \mathbf{E}$ being $\underline{R}^E$. One computes here a search direction in space using as conditioner the symmetric operator $\underline{A}_E^0$.

Let be $(\gamma^E(\underline{\mu}), \underline{\tilde{Z}}^E) \in \mathbb{R} \times \mathbb{R}^n$; the corresponding space is $\Gamma \times \mathbb{R}^n$. One defines:
$$ (\gamma^E(\underline{\mu}), \tilde{\underline{Z}}^E) = \arg \left[ \min_{\gamma' \in \Gamma, \tilde{\underline{Z}}' \in \mathbb{R}^n} \left\langle (\underline{R}^E - \gamma' \mathbf{A}_E^0 \tilde{\underline{Z}}') \mathbf{A}_E^{0-1} (\underline{R}^E - \gamma' \mathbf{A}_E^0 \tilde{\underline{Z}}') \right\rangle \right] $$

$$ \text{where } \langle \cdot, \cdot \rangle = \int_{\Sigma_{\mu}} \cdot d\underline{\mu} $$

It follows:

$$ \tilde{\underline{Z}}^E = \arg \max_{\tilde{\underline{Z}}' \in \mathbb{R}^n, \tilde{\underline{Z}}' \neq \underline{Z}} \underline{Z}'^T \langle \underline{R}^E \underline{R}^E \rangle \underline{Z}' $$

The $\underline{\mu}$-integration can be done using the macro description.

$$ \bullet \text{ Computation of the extension } \tilde{\underline{X}} = \bigwedge_{E \in \mathbb{E}} a_E \tilde{\underline{Z}}^E \text{ with } \underline{a} = \mathbf{J}(\underline{R}, \tilde{\underline{Z}}) $$

From $\tilde{\underline{Z}}^E$, $E \in \mathbb{E}$, one defines a space search direction over $\Omega$. One has:

$$ (a_E, E \in \mathbb{E}) = \arg \left[ \min_{\gamma' \in \Gamma, a'_E \in \mathbb{R}} \sum_{E \in \mathbb{E}} \left\langle (\underline{R}^E - \gamma' a'_E \mathbf{A}_E^0 \tilde{\underline{Z}}^E) \mathbf{A}_E^{0-1} (\underline{R}^E - \gamma' a'_E \mathbf{A}_E^0 \tilde{\underline{Z}}^E) \right\rangle \right] $$

which is equivalent to find the maximum of the Rayleigh quotient:

$$ \frac{\underline{a}^T \mathbf{M} \underline{a}}{\underline{a}^T \underline{a}} $$

which is relatively easy to compute.

$$ \bullet \text{ Computation of the "macro" functions } f^M = \mathbf{K}(\underline{R}, \tilde{\underline{X}}) $$

One minimize the residual:

$$ f^M = \arg \left[ \min_{f' \in \Gamma^M} \sum_{E \in \mathbb{E}} \left\langle (\underline{R}^E - [\mathbf{A}(f'\tilde{\underline{X}})]_{|E})^T \mathbf{A}_E^{0-1} (\underline{R}^E - [\mathbf{A}(f'\tilde{\underline{X}})]_{|E}) \right\rangle \right] $$

and then

$$ f^M = \arg \left[ \max_{f' \in \Gamma} \mathcal{R}(f') \right] $$

where

$$ \mathcal{R}(f) = \frac{\left\langle f \sum_{E \in \mathbb{E}} \underline{R}^E \underline{A}_E^{0-1} \tilde{\underline{X}}^E \right\rangle}{\left\langle f^2 [\mathbf{A}(\tilde{\underline{X}})]_{|E}^T \mathbf{A}_E^{0-1} \right\rangle [\mathbf{A}(\tilde{\underline{X}})]_{|E}} $$

The classical PGD technique is used with for example two complete iterations over the parameter space.
• Computation of the “micro” functions {$g^{EE'}, E' \in \bar{C}_E$} = $L_E(\underline{R}, \tilde{\underline{X}}, f^M)$
Let us consider that the initial residual is:

$$ \underline{R}^0 = \underline{R} - A(f^M \tilde{\underline{X}}) $$

One introduces a new approximation which is equal to $f^M \tilde{\underline{X}}$ over the complementary part of $\bar{C}_E$:

$$ \gamma^E(\underline{\mu}) = f^M(\underline{\mu}) \prod_{E' \in C_E} \frac{g^{EE'}(\mu_{E'})}{f^M_{E'}(\mu_{E'})} $$

This residual minimization problem which is here a small problem is solved classically.

### 3.4.2 Computational Strategy

The computational strategy uses error indicators as:

• Global error: $\epsilon_n = \frac{|||(\underline{F}_d - \underline{A}\tilde{\underline{X}}^n)|||$}{|||\underline{F}_d|||}$

• Norm of the residual $\underline{R}^{E,n}: r_{E,n} = \frac{|||\underline{R}^{E,n}|||_E}{\sup_{E \in \bar{E}} |||\underline{F}_d|||_E}$

where $||| \cdot |||_E^2 = \langle \cdot^T \underline{A}_E^0 \cdot^{-1} \rangle$

The optimal computational strategy is under study. More details can be found in [17, 25].

## 3.5 A First Illustration

The parameter-multiscale PGD has been implemented on a monodimensional example. A displacement is imposed at the extremity of a cantilever beam and each element has an independent Young modulus. A solution for a random set of parameters can be seen Fig. 9. This one-dimension example respects the Saint-Venant principle in stress, but not in displacement, requiring the addition of a rigid-body displacement associated to each micro function. This specificity will disappear in a two or three dimensional space.

First, we have computed the errors associated to the different approximations following the development (21) in the case where the variations of $E$ are $\pm 50\%$. Let us introduce the following norms and errors which are $E$-independent for the studied problem:
<figure><img src="image_11.png" /><figcaption>Fig. 9 Displacement in a beam, each element has an independent and random Young Modulus</figcaption></figure>

$$ \cdot \varepsilon^2 = \frac{\sup_{E \in \mathbb{E}} \int_{\Sigma_u} ||\sigma_{ex} - \sigma||_E^2 d\mu}{\sup_{E \in \mathbb{E}} \int_{\Sigma_\mu} ||\sigma_{ex}||_E^2 d\mu} $$

$$ \cdot \frac{\bar{\varepsilon}^2}{\varepsilon} = \frac{\sup_{\mu \in \Sigma_\mu, E \in \mathbb{E}} ||\sigma_{ex} - \sigma||_E^2}{\sup_{\mu \in \Sigma_\mu, E \in \mathbb{E}} ||\sigma_{ex}||_E^2} $$

$$ \text{with } ||\sigma||_E^2 = \int_E \frac{1}{E_0} \sigma(U)^2 dx $$

One gets:

* 0-order approximation: $\bar{\varepsilon} = 0.5 \quad \varepsilon = 0.29$

* 1-order approximation: $\bar{\varepsilon} = 0.25 \quad \varepsilon = 0.10$

* 2-order approximation: $\bar{\varepsilon} = 0.11 \quad \varepsilon = 0.042$

Then, the parameter-multiscale PGD has been implemented with the following strategy on a problem involving 64 parameters. For each element, one starts the computation with a micro-function corrected by a macro-one. The error at this stage is already small:

$$ \cdot \bar{\varepsilon} = 0.12 \quad \varepsilon = 0.07 $$

Let us note that such an approximation gives the same error for any number of parameters. After introducing macro-functions, the error continues to decrease but slowly.
# 4 Conclusion

The parameter-multiscale PGD seems to be a very promising way to build reduced order model for problems involving a large number of parameters. Further work will be devoted to the derivation of verification tools and to the extension to nonlinear problems such as viscoplastic ones.

# References

1. A. Ammar, B. Mokdad, F. Chinesta, R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids. Part II: Transient simulation using space-time separated representations. J. Non-Newton. Fluid Mech. **144**(2–3), 98–121 (2007)

2. M. Barrault, Y. Maday, N.C. Nguyen, A.T. Patera, An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations. C. R. Math. **339**(9), 667–672 (2004)

3. M. Billaud-Friess, A. Nouy, O. Zahm, A tensor approximation method based on ideal minimal residual formulations for the solution of high-dimensional problems. ESAIM. Math. Mode. Numer. Anal. **48**(6), 1777–1806 (2014)

4. M. Capaldo, P.-A. Guidault, D. Néron, P. Ladevèze, The Reference Point Method, a "hyper-reduction" technique: application to PGD-based nonlinear model reduction. Comput. Method. Appl. Mech. Eng. **322**, 483–514 (2017)

5. S. Chaturantabut, D.C. Sorensen, Nonlinear model reduction via discrete empirical interpolation. SIAM J. Sci. Comput. **32**(5), 2737–2764 (2010)

6. F. Chinesta, P. Ladevèze, (eds.) *Separated Representations and PGD-Based Model Reduction: Fundamentals and Applications*, vol. CISM 554 (Springer, 2014)

7. E. de Souza Neto, D. Peric, D.R.J. Owen, *Computational Methods for Plasticity*, vol. 55 (2008)

8. A. Falco, W. Hackbusch, A. Nouy, *Geometric Structures in Tensor Representations* (2015) pp. 1–50 (Work document)

9. W. Hackbusch, *Tensor Spaces and Numerical Tensor Calculus* (Springer, 2012)

10. J.A. Hernandez, J. Oliver, A.E. Huespe, M.A. Caicedo, J.C. Cante, High-performance model reduction techniques in computational multiscale homogenization. Comput. Methods Appl. Mech. Eng. **276**, 149–189 (2014)

11. C. Heyberger, P.A. Boucard, D. Néron, A rational strategy for the resolution of parametrized problems in the PGD framework. Comput. Methods Appl. Mech. Eng. **259**, 40–49 (2013)

12. S. Holtz, T. Rohwedder, R. Schneider, On manifolds of tensors of fixed TT-rank. Numer. Math. **120**(4), 701–731 (2012)

13. P. Ladevèze, On Algorithm Family in Structural Mechanics. Comptes rendus des séances de l'Académie des sciences. Série 2, **300**(2) (1985)

14. P. Ladevèze, The large time increment method for the analyse of structures with nonlinear constitutive relation described by internal variables. Comptes rendus des séances de l'Académie des sciences. Série 2, **309**(2), 1095–1099 (1989) (in french)

15. P. Ladevèze, *Nonlinear Computational Structural Mechanics: New Approaches and Non-incremental Methods of Calculation* (Springer, New York, 1999)

16. P. Ladevèze, New variational formulations for discontinuous approximations. Technical Report, LMT Cachan, 2011, (in french)

17. P. Ladevèze, A new method for the ROM computation: the parameter-multiscale PGD, Technical report, LMT Cachan, 2016a, (in french)

18. P. Ladevèze, On reduced models in nonlinear solid mechanics. Eur. J. Mech. A/Solids **60**, 227–237 (2016b)
19. P. Ladevèze, J.C. Passieux, D. Néron, The LATIN multiscale computational method and the proper generalized decomposition. Comput. Methods Appl. Mech. Eng. **199**(21–22), 1287–1296 (2010)

20. P. Ladevèze, H. Riou, On Trefftz and weak Trefftz discontinuous Galerkin approaches for medium-frequency acoustics. Comput. Methods Appl. Mech. Eng. **278**, 729–743 (2014)

21. E. Monteiro, J. Yvonnet, Q.C. He, Computational homogenization for nonlinear conduction in heterogeneous materials using model reduction. Comput. Mater. Sci. **42**(4), 704–712 (2008)

22. D. Néron, P.A. Boucard, N. Relun, Time-space PGD for the rapid solution of 3D nonlinear parametrized problems in the many-query context. Int. J. Numer. Methods Eng. **103**(4), 275–292 (2015)

23. A. Nouy, A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations. Comput. Methods Appl. Mech. Eng. **199**(23–24), 1603–1626 (2010)

24. I.V. Oseledets, Tensor-train decomposition. SIAM J. Sci. Comput. **33**(5), 2295–2317 (2011)

25. C. Paillet, P. Ladevèze, D. Néron, A Parametric-multiscale PGD for problems with a large number of parameters (2017) (In preparation)

26. A. Radermacher, S. Reese, POD-based model reduction with empirical interpolation applied to nonlinear elasticity. Int. J. Numer. Methods Eng. **107**(6), 477–495 (2016)

27. N. Relun, D. Néron, P.A. Boucard, A model reduction technique based on the PGD for elastic-viscoplastic computational analysis. Comput. Mech. **51**(1), 83–92 (2013)

28. D. Ryckelynck, Hyper-reduction of mechanical models involving internal variables. Int. J. Numer. Methods Eng. **77**(1), 75–89 (2009)