# PGD-BEM Applied to the Nonlinear Heat Equation

Pierre Joyot, Nicolas Verdon, Gaël Bonithon, Francisco Chinesta, Pierre Villon

## ► To cite this version:

Pierre Joyot, Nicolas Verdon, Gaël Bonithon, Francisco Chinesta, Pierre Villon. PGD-BEM Applied to the Nonlinear Heat Equation. ASME 2012 11th Biennial Conference on Engineering Systems Design and Analysis, Jul 2012, Nantes, France. pp.205-213, 10.1115/ESDA2012-82407. [hal-01008663](https://hal.science/hal-01008663)

HAL Id: hal-01008663

https://hal.science/hal-01008663

Submitted on 17 May 2023

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
PGD-BEM APPLIED TO THE NONLINEAR HEAT EQUATION

Pierre Joyot*

Nicolas Verdon

ESTIA-Recherche
Technopole izarbel, 64210 Bidart
France

Email: {pjoyot,nverdon}@estia.fr

Gaël Bonithon

GeM CNRS-ECN
Ecole Centrale de Nantes
44321 Nantes Cedex 3

France

Email: gael.bonithon@ec-nantes.fr

Francisco Chinesta

EADS Corporate Foundation International Chair
GeM CNRS-ECN

44321 Nantes Cedex 3
France

Email: Francisco.Chinesta@ec-nantes.fr

Pierre Villon

Laboratoire Roberval UMR 6253
University of Technology of Compiègne
60200 Compiègne

France

Email: pierre.villon@utc.fr

ABSTRACT

The Boundary Element Method (BEM) allows efficient solution of partial differential equations whose kernel functions are known. The heat equation is one of these candidates when the thermal parameters are assumed constant (linear model). When the model involves large physical domains and time simulation intervals the amount of information that must be stored increases significantly. This drawback can be circumvented by using advanced strategies, as for example the multi-poles technique. We propose radically different approach that leads to a separated solution of the space and time problems within a non-incremental integration strategy. The technique is based on the use of a space-time separated representation of the unknown field that, introduced in the residual weighting formulation, allows to define a separated solution of the resulting weak form. The spatial step can be then treated by invoking the standard BEM for solving the resulting steady state problem defined in the physical space. Then, the time problem that results in an ordinary first order differential equation is solved using any standard approach.

backward finite differences). When considering the nonlinear heat equation, the BEM cannot be easily applied because its Green's kernel is generally not known but the use of the PGD presents the advantage of rewriting the problem in such a way that the kernel is now clearly known. Indeed, the system obtained by the PGD is composed of a Poisson equation in space coupled with an ODE in time so that the use of the BEM for solving the spatial part of the problem is efficient. During the solving, we must however separate the nonlinear term into a space-time representation that can limit the method in terms of CPU time and storage, that is why we introduce in the second part of the paper a new approach combining the PGD and the Asymptotic Numerical Method (ANM) in order to efficiently treat the nonlinearity.

INTRODUCTION

Solving the nonlinear heat equation with the Boundary Element Method is not straightforward because it is generally not possible to determine its Green's kernel. A direct resolution method based on the kernel of linear heat equation has been

<footnote>*Address all correspondence to this author.</footnote>
proposed in [1] and some other approaches have been used, we could just mention here the homotopy method used quite recently in [2]. The approach we use in this paper is radically different because it involves Reduced Order Modelling techniques, from which we can obtain a separated representation of the solution. One of the pioneer ROM methods, the Proper Orthogonal Decomposition, allows to obtain a space-time decomposition of the solution from a set of samples of the solution called snapshots. But this sampling step limits the performances of the method that is why the Proper Generalized Decomposition is nowadays preferred. Based on a *a priori* separated representations of the unknown, the PGD allows to treat high dimensional problems arising from multidimensional physics and its range of application is very large, see [3] for a recent review. For the application in the nonlinear heat equation, it will consist in determining a space-time separated representation with an iterative procedure: the algorithm computes at each iteration a couple of modes from the modes that have already been calculated. More details about the PGD can be found in [4–6]. The advantage of coupling the PGD and the BEM lies on the fact that is no more necessary to know the space-time kernel of the heat equation but only the Green’s kernel of the Poisson equation. The method, called in the later PGD-BEM, has been first introduced in [7] for the linear heat equation and we show in this paper that its algorithm is easily adapted to the nonlinear case. Even if the developed algorithm shows impressive results, it remains some limitations due to the cost of the space-time separation of the nonlinear term, that is why we focus our work in the second part of the paper on the Asymptotic Numerical Method, see for instance [8]. Originally designed for tracking bifurcation points, it has been lately considered as a good alternative to circumvent the problems due to high nonlinearity. Combined with the POD, Yvonnet et al. [9] and more recently Niroomandi et al. [10] applied with success a POD-ANM strategy for respectively post-buckling analysis and large deformations of hyperelastic materials. In this paper, we will present an approach coupling PGD and ANM that takes advantage of both methods. The performances of the different approaches will be outlined with a numerical example.

## THE PGD-BEM METHOD

In this section, we present the PGD-BEM method for solving the nonlinear heat equation. Here, the nonlinearity comes from the thermal conductivity function which depends on the temperature. Hence, the problem that will be treated throughout this paper is recalled here:

$$ \left\{ \begin{array}{ll} u_{,t} - (k(u)u_{,i})_{,i} = b & \text{in } \Omega \times \tau \\ u = 0 & \text{in } \partial\Omega \times \tau \\ u(t=0) = 0 & \text{in } \Omega \end{array} \right. \quad (1) $$

where $\Omega$ is the spatial domain, $\partial\Omega$ its boundary, $\tau$ the time interval, $u$ the temperature field, $b$ the source term and $k$ the thermal conductivity function.

The main idea of the PGD is to decompose the solution of the problem (1) as a sum of products of spatial functions $X^\alpha(\mathbf{x})$ and temporal functions $T^\alpha(t)$, namely we seek $u$ as follows:

$$ u(\mathbf{x},t) \approx \sum_{\alpha=1}^{N} X^{\alpha}(\mathbf{x}) T^{\alpha}(t) \quad (2) $$

In order to explain what we will call in the following the direct PGD, let us assume the $p-1$ first couples of $X^\alpha T^\alpha$ known. The solution $u^{(p)}$ at iteration $p$ is then expressed $u^{(p)} = u^{(p-1)} + XT$ and the functions $X$ and $T$ are obtained as solutions of the following nonlinear system:

$$ \int_{\tau} (u_{,t}^{(p)} T) d\tau - \int_{\tau} \left( \left( k(u^{(p)}) u_{,i}^{(p)} \right)_{,i} T \right) d\tau = \int_{\tau} b T d\tau \quad (3) $$

$$ \int_{\Omega} (u_{,t}^{(p)} X) dx - \int_{\Omega} \left( \left( k(u^{(p)}) u_{,i}^{(p)} \right)_{,i} X \right) dx = \int_{\Omega} b X dx \quad (4) $$

Equation (3) is a PDE in space that will be solved using the BEM whereas equation (4) is an ODE in time. This system of coupled equations is solved using a fixed-point method that will be discussed later.

## Spatial resolution with the BEM

As the spatial equation of the system is solved using the BEM, we need to rewrite (3) in a way that allows to compute Green’s functions. Hence, the system (3)-(4) is rewritten as follows:

$$ \begin{aligned} - \left[ \int_{\tau} (k(u^{(p)}) T T) d\tau \right]_{,ii} &= \int_{\tau} B T d\tau - \left[ \int_{\tau} T_{,t} T d\tau \right] X \\ &\quad - \sum_{\alpha=1}^{p-1} \left[ \int_{\tau} T_{,t}^{\alpha} T d\tau \right] X^{\alpha} \\ - \left[ \int_{\Omega} (k(u^{(p)}) X_{,ii} X) dx \right] T &= \int_{\Omega} B X dx - \left[ \int_{\Omega} X X dx \right] T_{,t} \\ &\quad - \sum_{\alpha=1}^{p-1} \left[ \int_{\Omega} X^{\alpha} X dx \right] T_{,t}^{\alpha} \end{aligned} \quad (5) $$

where $B = b + k'(u^{(p)})u_{,i}^{(p)}u_{,i}^{(p)} + k(u^{(p)})u_{,ii}^{(p-1)}$. For the case of homogeneous boundary conditions, we can write $\int_{\tau} T_{,t} T d\tau = \frac{1}{2} T_f^2$ where $T_f = T(\tau_f)$. The functions $k(u^{(p)})$ and $B$ depend on $x$ and $t$ that is why they must be decomposed into separated forms
for computing the integrals in system (5). This is practically achieved by using a SVD, so that if we note:

$$k(u^{(p)}) \approx \sum_{\alpha=1}^{p} K_x^{\alpha} K_t^{\alpha} \quad \text{and} \quad B \approx \sum_{\alpha=1}^{p} B_x^{\alpha} B_t^{\alpha} \qquad (6)$$

the system can finally be written as:

$$- \left[ \sum_{\alpha=1}^{p-1} K_x^{\alpha} \int_{\tau} K_t^{\alpha} T \, d\tau \right] X_{,ii} - \left[ \sum_{\alpha=1}^{p-1} B_x^{\alpha} \int_{\tau} B_t^{\alpha} T \, d\tau \right] + \left[ \frac{1}{2} T_f^2 \right] X + \sum_{\alpha=1}^{p-1} \left[ \int_{\tau} T_t^{\alpha} T \, d\tau \right] X^{\alpha} = 0 \qquad (7)$$

$$- \left[ \sum_{\alpha=1}^{p-1} K_t^{\alpha} \int_{\Omega} K_x^{\alpha} X_{,ii} X \, dx \right] T - \left[ \sum_{\alpha=1}^{p-1} B_t^{\alpha} \int_{\Omega} B_x^{\alpha} X \, dx \right] + \left[ \int_{\Omega} X X \, dx \right] T_t + \sum_{\alpha=1}^{p-1} \left[ \int_{\Omega} X^{\alpha} X \, dx \right] T_t^{\alpha} = 0 \qquad (8)$$

## The fixed-point algorithm

The system of equations (7)-(8) is then solved using a fixed-point method. Assuming that all fields are known at iteration $q$, the way of computing them at iteration $q+1$ is described hereafter.

We first determine $X^{q+1}$ by solving equation (7), note that it is a Poisson-like equation and hence is efficiently solved by the BEM. So, we have to solve:

$$-X_{,ii}^{q+1} = \frac{1}{\sum_{\alpha=1}^{p-1} K_x^{\alpha} \int_{\tau} K_t^{\alpha} T^q T^q d\tau} \left[ \sum_{\alpha=1}^{p-1} B_x^{\alpha} \int_{\tau} B_t^{\alpha} T^q d\tau - \left(\frac{1}{2}(T_f^q)^2\right) X - \sum_{\alpha=1}^{p-1} \left[ \int_{\tau} T_t^{\alpha} T^q d\tau \right] X^{\alpha} \right] \qquad (9)$$

The terms $B_x^\alpha$ and $B_t^\alpha$ are then evaluated by achieving a SVD over the following form of $B$:

$$B = b + k' \left( u^{(p-1),q,q} \right) u_{,i}^{(p-1),q,q} u_{,i}^{(p-1),q,q} + k \left( u^{(p-1),q,q} \right) u_{,ii}^{(p-1)} \qquad (10)$$

where $u^{(p-1),q,q} = u^{(p-1)} + X^q T^q$. In a similar manner, $K_x^\alpha$ and $K_t^\alpha$ are obtained from $k(u^{(p-1),q,q})$. In a second step, $T^{q+1}$

is obtained by solving (8), which gives at iteration $q+1$:

$$- \left[ \sum_{\alpha=1}^{p-1} K_t^{\alpha} \int_{\Omega} K_x^{\alpha} X_{,ii}^{q+1} X^{q+1} dx \right] T^{q+1} - \left[ \sum_{\alpha=1}^{p-1} B_t^{\alpha} \int_{\Omega} B_x^{\alpha} X^{q+1} dx \right] + \left[ \int_{\Omega} X^{q+1} X^{q+1} dx \right] T_{,t}^{q+1} + \sum_{\alpha=1}^{p-1} \left[ \int_{\Omega} X^{\alpha} X^{q+1} dx \right] T_{,t}^{\alpha} = 0 \qquad (11)$$

Once again, $B_x^\alpha$ and $B_t^\alpha$ are evaluated from the following form:

$$B = b + k' \left( u^{(p-1),q+1,q} \right) u_{,i}^{(p-1),q+1,q} u_{,ii}^{(p-1)} \qquad (12)$$

whereas $K_x^\alpha$ and $K_t^\alpha$ are obtained from $k(u^{(p-1),q+1,q})$, where $u^{(p-1),q+1,q} = u^{(p-1)} + X^{q+1}T^q$.

Contrary to the linear PGD-BEM presented in [7], it is here necessary to compute the separated forms of $B$ and $k(u)$ during the fixed-point loop. The way of evaluating the nonlinear term is very important for the efficiency of the algorithm but also affects the computational effort needed for the whole method. This aspect has been widely discussed in [11] where they test different linearization using the solution at previous iteration. Here, we linearize in a more implicit way because we separate the nonlinear term inside the fixed-point loop with a SVD. It involves then a number of SVD which is equal to the number of iterations until convergence of the fixed-point algorithm times the number of PGD modes we keep in the solution. The CPU time and storage can hence become non negligible so that we propose in the following another approach that avoids dealing with these separations.

## AN ALTERNATIVE: THE PGD-ANM APPROACH

In this section we present the methodology we have employed in order to circumvent the need of separating the conductivity function into a product of space and time functions.

## The Asymptotic Numerical Approach

The Asymptotic Numerical Method (ANM) has been introduced in the 90's to deal with strong nonlinear problems. It consists in transforming the nonlinear problem into a recursive sequence of well posed linear problems having the same tangent operator. In comparison with classical methods that are used for solving nonlinear problems, e.g. Newton-Raphson methods, the ANM has proved to be much more efficient in terms of computational costs. More details about the wide range of applications can be found in [8].

For generality's sake, let us briefly explain the principle of the
method by considering the following nonlinear problem written in quadratic form:

$$ \left\{ \begin{array}{ll} \mathcal{L}(u) + \mathcal{Q}(u, u) = \lambda b & \text{in } \Omega \\ u = 0 & \text{in } \partial\Omega \end{array} \right. \qquad (13) $$

where $\mathcal{L}(u)$ represents a linear operator and $\mathcal{Q}(u, u)$ a quadratic one, $b$ is the second member and $\lambda$ a load parameter. The main idea of the ANM is the asymptotic expansion around a known solution $(u^n, \lambda_0)$ – in our case $\lambda_0$ is always zero – at step $n$ into power series with respect to a path parameter $a$:

$$ \left\{ \begin{array}{l} u^{n+1}(a) = u^n + au_1 + a^2 u_2 + \dots \\ \lambda(a) = a\lambda_1 + a^2\lambda_2 + \dots \end{array} \right. \qquad (14) $$

When truncating the series at order $N$, the problem (13) gives:

$$ \mathcal{L} \left( u^n + \sum_{p=1}^{N} a^p u_p \right) + \mathcal{Q} \left( u^n + \sum_{p=1}^{N} a^p u_p, u^n + \sum_{p=1}^{N} a^p u_p \right) = \left( \sum_{p=1}^{N} a^p \lambda_p \right) b \qquad (15) $$

and by identifying the terms of same power $p$, we obtain a recursive sequence of linear problems:

$$ \mathcal{L}_p(u_p) = \lambda_p b + \mathcal{F}_p \qquad (16) $$

where $\mathcal{F}_p$ is a nonlinear second member term that depends on terms until order $p-1$. An additional condition is obtained when defining $a$ as a pseudo arc-length:

$$ a = (u^{n+1}(a) - u^n) \cdot u_1 + \lambda(a) \lambda_1 \qquad (17) $$

Relations (16) and (17) after a given discretization procedure allow to write a sequence of linear problems in the form:

$$ \text{order 1} \left\{ \begin{array}{l} \mathbf{K}_t \mathbf{u}_1 = \lambda_1 \mathbf{f} \\ \mathbf{u}_1^T, \mathbf{u}_1 + \lambda_1^2 = 1 \end{array} \right. \qquad (18) $$

$$ \text{order } p \left\{ \begin{array}{l} \mathbf{K}_t \mathbf{u}_p = \lambda_p \mathbf{f} + \mathbf{f}_p \\ \mathbf{u}_p^T \mathbf{u}_1 + \lambda_p \lambda_1 = 0 \end{array} \right. \qquad (19) $$

where $\mathbf{K}_t$ is the tangent stiffness matrix, $\mathbf{f}$ is the loading vector and $\mathbf{f}_p$ is the second term vector obtained from the discretisation of $\mathcal{F}_p$. These problems are solved using a prediction-correction approach, their solutions are obtained as follows:

$$ \text{order 1} \left\{ \begin{array}{l} \hat{\mathbf{u}}_1 = \mathbf{K}_t^{-1} \mathbf{f} \\ \lambda_1 = \frac{1}{\sqrt{(\hat{\mathbf{u}}_1, \hat{\mathbf{u}}_1)} + 1} \\ \mathbf{u}_1 = \lambda_1 \hat{\mathbf{u}}_1 \end{array} \right. \qquad (20) $$

$$ \text{order } p \left\{ \begin{array}{l} \hat{\mathbf{u}}_p = \mathbf{K}_t^{-1} \mathbf{f}_p \\ \lambda_p = -\lambda_1 (\hat{\mathbf{u}}_p, \mathbf{u}_1) \\ \mathbf{u}_p = \frac{\lambda_p}{\lambda_1} \mathbf{u}_1 + \hat{\mathbf{u}}_p \end{array} \right. \qquad (21) $$

## Coupling between PGD and ANM

There are obviously several ways of coupling the PGD and the ANM. A first strategy consists in applying the PGD on the suite of linear systems that would be obtained by the ANM on the nonlinear heat equation. The other approach takes advantage of the ANM as a nonlinear solver, so that the order of the steps is different: the ANM is applied inside the PGD algorithm. The first approach is a work currently in progress and we will here briefly present the most important features of the second approach we have implemented.

The first steps of the algorithm remain unchanged, so that we deal with the system of coupled equations (3)-(4) that are obtained thanks to the PGD procedure. Here, instead of separating the nonlinear term $k(u)$ into a couple of space-time functions, we propose to consider it as a new unknown that will be expressed as an asymptotic expansion. If we note $v = k(u)$, the expansion will be now written around the known solution $(u^n, v^n, \lambda^n)$ at step $n$:

$$ \begin{pmatrix} \mathbf{u}^{n+1} \\ \mathbf{v}^{n+1} \\ \lambda \end{pmatrix} = \begin{pmatrix} \mathbf{u}^n \\ \mathbf{v}^n \\ 0 \end{pmatrix} + \sum_{p=1}^{N} a^p \begin{pmatrix} \mathbf{u}_p \\ \mathbf{v}_p \\ \lambda_p \end{pmatrix} \qquad (22) $$

so that the PGD method will determine at the same time the separated representation of $u$ and $k(u)$. Typically, knowing $u$ and $v$ at iteration $p-1$, we seek their expression at iteration $p$ as:

$$ u^{(p)} = u^{(p-1)} + XT \qquad (23) $$

$$ v^{(p)} = v^{(p-1)} + RS \qquad (24) $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. EXACT AND COMPUTED TEMPERATURE.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>FIGURE 2. ABSOLUTE ERROR IN THE TEMPERATURE.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 3. FIRST THREE COUPLES OF MODES $(X^\alpha, T^\alpha)_{\alpha=1,2,3}$.</figcaption></figure>
and the ANM will be applied inside the PGD step over X and R so that $u^{(p)}$ and $v^{(p)}$ will be then expressed as:

$$u^{(p)} = u^{(p-1)} + \left( \sum_{\alpha=1}^{N} a_{\alpha} X^{\alpha} \right) T \quad (25)$$

$$v^{(p)} = v^{(p-1)} + \left( \sum_{\alpha=1}^{N} a_{\alpha} R^{\alpha} \right) S \quad (26)$$

Rewriting these expressions into the direct PGD algorithm leads finally to the PGD-ANM method.

## RESULTS

In the later, the nonlinear heat equation will be solved in the spatial domain $\Omega = [0, 1] \times [0, 1]$ and over the time domain $\tau = [0, 1]$. For comparison reasons, the analytic solution is imposed to be $u^{ex} = x(x-1)y(y-1)t$ and the second member is computed from this solution, depending on the value of the nonlinear term $k(u)$. For the PGD-BEM, we consider the nonlinear term of the following form $k(u) = (u^2 + 1)$. Each side of the domain $\Omega$ is discretized using 8 nodes and the time interval is decomposed into 129 time steps.

The figure 1 shows the results obtained with the PGD-BEM approach compared to the analytic solution at times $t = \{0s, 7.8 \times 10^{-3}s, 0.25s, 0.5s, 0.75s, 1s\}$ and figure 2 represents the corresponding absolute error computed from the analytic solution:

$$e(\mathbf{x}, t) = \left| \frac{u(\mathbf{x}, t) - u^{ex}(\mathbf{x}, t)}{u^{ex}(\mathbf{x}, t)} \right| \quad \forall (\mathbf{x}, t) \in \Omega \times \tau \quad (27)$$

At last, figure 3 shows the first three couples of modes that are obtained at the end of the algorithm, note here that the first mode is able to represent the structure of the solution whereas the additional modes are useful for catching accurately the details.

Now, we present the results obtained by applying the PGD-ANM method compared to the direct PGD, recalling here that this method avoids achieving a large number of SVD which is the case when using the direct PGD. We keep the same value of the analytic solution but now we consider the nonlinear term defined as $k(u) = (u + \varepsilon)^n$ where $n$ takes the values $n = 1$ or $2$. Figure 4 presents the evolution of the error in a $L^2_{\Omega \times \tau}$ sense with the number of PGD modes for the direct PGD compared to the PGD-ANM method at respectively order 1 and 2. We observe in this figure that the PGD-ANM at first order slightly improves the results in comparison with the direct PGD whereas they are significantly better at second order, even for the most nonlinear case. For this case, with 5 PGD modes, the error is divided by two when using the PGD-ANM instead of the direct PGD, which is already a quite interesting results.

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

<figure><img src="image_6.png" /><figcaption>FIGURE 4. COMPARISON BETWEEN PGD DIRECT AND PGD-ANM FOR DIFFERENT VALUES OF $k(u)$.</figcaption></figure>

## CONCLUSIONS

The work presented in this paper addresses some issues or at least interesting perspectives in the field of efficient solving of nonlinear equations. Indeed, the PGD-BEM algorithm has been successfully adapted for solving the nonlinear heat equation, without increasing the level of complexity in the methodology for the linear case. The main advantage of this approach lies on the fact that only the Green's kernel of Poisson equations is required to solve the nonlinear heat equation with the Boundary Element Method. Furthermore, coupling the PGD with the ANM allows to take advantage of both techniques, and the gains we can expect are really encouraging. Optimizations in the algorithm are currently in progress in order to make it even more accurate and CPU time efficient.
REFERENCES

[1] Goto, T., and Suzuki, M., 1996. "A boundary integral equation method for nonlinear heat conduction problems with temperature-dependent material properties". *International Journal of Heat and Mass Transfer*, **39**(4), Mar., pp. 823-830.

[2] Shijun, L., 2010. "An optimal homotopy-analysis approach for strongly nonlinear differential equations". *Communications in Nonlinear Science and Numerical Simulation*, **15**(8), Aug., pp. 2003-2016.

[3] Chinesta, F., Ladeveze, P., and Cueto, E., 2011. "A short review on model order reduction based on proper generalized decomposition". *Archives of Computational Methods in Engineering*, **18**(4), Nov., pp. 395-404.

[4] Ladevèze, P., Passieux, J. C., and Neron, D., 2010. "The latin multiscale computational method and the proper orthogonal decomposition". *Computer Methods in Applied Mechanics Engineering*, **199**(21-22), pp. 1287-1296.

[5] Chinesta, F., Ammar, A., Falco, A., and Laso, M., 2007. "On the reduction of stochastic kinetic theory models of complex fluids". *Modeling Simulation in Materials Science Engineering*, **15**, pp. 639-652.

[6] Nouy, A., 2010. "A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations". *Computer Methods in Applied Mechanics and Engineering*, **199**(23-24), pp. 1603-1626.

[7] Bonithon, G., Joyot, P., Chinesta, F., and Villon, P., 2011. "Non-incremental boundary element discretization of parabolic models based on the use of the proper generalized decompositions". *Engineering Analysis with Boundary Elements*, **35**(1), Jan., pp. 2-17.

[8] Abichou, H., Zahrouni, H., and Potier-Ferry, M., 2002. "Asymptotic numerical method for problems coupling several nonlinearities". *Computer Methods in Applied Mechanics and Engineering*, **191**(51-52), Dec., pp. 5795-5810.

[9] Yvonnet, J., Zahrouni, H., and Potier-Ferry, M., 2007. "A model reduction method for the post-buckling analysis of cellular microstructures". *Computer Methods in Applied Mechanics and Engineering*, **197**(1-4), Dec., pp. 265-280.

[10] Niroomandi, S., Alfaro, I., Cueto, E., and Chinesta, F., 2010. "Model order reduction for hyperelastic materials". *Int. J. Numer. Methods Eng.*, **81**(9), pp. 1180-1206.

[11] Pruliere, E., Chinesta, F., and Ammar, A., 2010. "On the deterministic solution of multidimensional parametric models using the proper generalized decomposition". *Mathematics and Computers in Simulation*, **81**(4), Dec., pp. 791-810.