# Méthode PGD pour non linéarités matériaux appliquée à la fonderie

Pierre Despret, Jean-Luc Dulong, Frédéric Druesne, Pierre Villon

► To cite this version:

Pierre Despret, Jean-Luc Dulong, Frédéric Druesne, Pierre Villon. Méthode PGD pour non linéarités matériaux appliquée à la fonderie. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722053

HAL Id: hal-01722053

https://hal.science/hal-01722053v1

Submitted on 2 Mar 2018

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 CC0 1.0 - Universal - International License
# Méthode PGD pour non linéarités matériaux appliquée à la fonderie.

Pierre DESPRET <fnref n="1" />,<fnref n="2" />, Jean-Luc DULONG * <fnref n="1" />, Frédéric DRUESNE <fnref n="1" />, Pierre VILLON <fnref n="1" />

**Résumé** — La méthode de réduction de modèles PGD (Proper Generalized Decomposition) avec non linéarités matériaux présente des difficultés. Nous proposons un complément de formulation par une discrétisation spatio-temporelle des propriétés matériaux pour prendre en compte ces cas de non linéarités en fonction de la température. Cette formulation de la PGD sera appliquée à l'équation de la chaleur pour l'étude de la phase de solidification d'une éprouvette en aluminium utilisée en fonderie.

**Mots clés** — Réduction de modèles, séparation de variables, équation de la chaleur, solidification de l'aluminium.

## 1. Introduction

La simulation de la coulée en fonderie repose sur deux étapes : la phase de remplissage puis la phase de solidification lors du refroidissement. Lors de cette seconde phase, l'objectif est de détecter l'apparition des défauts de fonderie. Ces défauts sont localisés au moyen de l'évolution temporelle des champs de fraction liquide/solide du métal. Sur des pièces industrielles comme des culasses automobiles, de telles simulations sont réalisées avec des maillages très fins, de l'ordre de $10^6 – 10^7$ degrés de libertés (ddl) ainsi que des pas de temps réduits. Ceci nécessite le stockage d'une très grande quantité de données et, par conséquent, un temps de calcul important.

On se propose de réduire ce volume de stockage et le temps de calcul au moyen de la méthode de réduction de modèles PGD introduite par Ammar [1]. Cette méthode a déjà été employée pour la résolution de problèmes thermiques non linéaires mais théoriques.

Les travaux de Pruliere [8] traitent un cas où la conductivité du matériau varie linéairement en fonction de la température. Ammar [2] offre une proposition de résolution avec une non linéarité appliquée au second membre sous forme parabolique, tout comme Joyot [5]. Ghnatios [4] propose de résoudre l'équation de chaleur avec la conductivité thermique considérée comme un paramètre recherché supplémentaire. Nguyen [5] quant à lui, propose un couplage de la méthode PGD et de la Méthode Asymptotique Numérique pour la résolution d'un problème transitoire, où la non linéaire intervient au second membre, dans le terme source. Il applique également la PGD a un problème de thermodiffusion. Enfin, Nouy [7] étudie les cas $k(x)$ ou $k(t)$ mais pas $k(T(x,t))$.

Pour mieux prendre en compte les non linéarités matériaux dans le contexte industriel, nous nous proposons de mettre en place un complément de formulation à l'algorithmes PGD, notamment à travers une discrétisation élémentsofinis (EF) spatio-temporelle des propriétés matériaux afin de mettre en place la séparation de variables sur ceux-ci.

Nous présenterons nos résultats sur un cas de fonderie aluminium avec des propriétés matériaux non linéaires sans prendre en compte, dans un premier temps, l'effet de la chaleur latente. Ce travail est effectué dans le cadre du projet régional Picardie-ProFond avec deux partenaires industriels (Montupet SA et E.J.) leaders dans le domaine de la fonderie.

<footnote n="1">1 Laboratoire Roberval, Département GSM UTC, {prénom nom}(at) utc fr. rue Personne de Roberval, 60200 Compiègne</footnote>

<footnote n="2">2 Montupet SA, Service R&D. 3 rue de Nogent, 60290 Laigneville</footnote>

<footnote>* Auteur correspondant</footnote>
## 2. Méthodes

Les modélisations numériques des procédés de fonderie font appel à l'équation de la chaleur (Eq 1) pour simuler le changement de phase lors du refroidissement et de la solidification. Cette modélisation s'effectue sur l'intervalle temporel $\Omega t$ et sur un domaine spatial $\Omega e$, de frontière $\Gamma$, avec $\Gamma_D$ la frontière à température imposée $T_D(\mathbf{x})$ et $\Gamma_N$ la frontière à flux convectif imposé avec coefficient d'échange $h$ constant, de normale sortante $\mathbf{n}$ et à température extérieure $T_{ext}$ constante.

$$ 
\begin{align} 
\rho(T(\mathbf{x}, t))c_p(T(\mathbf{x}, t))\frac{\partial T(\mathbf{x}, t)}{\partial t} &= \nabla \cdot [k(T(\mathbf{x}, t))\nabla T(\mathbf{x}, t)] + Q \tag{1} \\ 
k(T(\mathbf{x}, t))\nabla T(\mathbf{x}, t) \cdot \mathbf{n} &= h(T(\mathbf{x}, t) - T_{ext}) \quad \text{sur } \Gamma_N \times \Omega t \nonumber \\ 
T(\mathbf{x}, t) &= T_D(\mathbf{x}, t) \quad \text{sur } \Gamma_D \times \Omega t \nonumber \\ 
T(\mathbf{x}, 0) &= T_o(\mathbf{x}) \quad \text{sur } \Omega e \nonumber 
\end{align} 
$$

Cette équation de la chaleur repose sur les propriétés matériaux de masse volumique $\rho$, de conductivité thermique $k$ et de capacité calorifique spécifique $c_p$, qui dépendent de la température $T(\mathbf{x}, t)$. Leurs valeurs sont d'origine expérimentale. Pour des raisons de propriété industrielle nous les avons normées (Fig 1). Les valeurs intermédiaires sont interpolées linéairement.

Afin de modéliser l'effet thermique de la chaleur latente de fusion $L$ lors du changement de phase, le terme source $Q$ est supposé proportionnel à la dérivée temporelle de la fraction solide $f_s$ et s'exprime sous la forme $\rho(T)L\frac{\partial(f_s(T))}{\partial t}$. D'après Dantzig [3], ce terme source peut être placé dans la capacité calorifique spécifique dont l'expression devient alors $c_p^{app}(T) = c_p(T) - L\frac{\partial(f_s(T))}{\partial T}$. Ainsi, nous nous ramenons à l'équation de la chaleur standard corrigée au moyen du $c_p^{app}$ :

$$ \rho(T(\mathbf{x}, t))c_p^{app}(T(\mathbf{x}, t))\frac{\partial T(\mathbf{x}, t)}{\partial t} = \nabla \cdot [k(T(\mathbf{x}, t))\nabla T(\mathbf{x}, t)]. \quad (2) $$

<figure><img src="image_2.png" /><figcaption>Fig 1: Évolution des propriétés matériaux en fonction de la température : conductivité thermique ($k$), masse volumique ($\rho$), fraction solide ($f_s$) et la capacité calorifique spécifique ($c_p$).</figcaption></figure>
La méthode PGD repose sur l'hypothèse d'une séparation des variables spatiales $\mathbf{x}$ et temporelles $t$.

Ainsi, nous pouvons approximer $T(\mathbf{x}, t) = \sum_{i} \varphi_{i}(\mathbf{x}) a_{i}(t)$ où $\varphi_{i}(\mathbf{x})$ est un mode spatial et $a_{i}(t)$ son amplitude au cours du temps. A l'issue de l'enrichissement $m$, nous connaissons $\varphi_{i}(\mathbf{x})$ et $a_{i}(t) \forall i \le m$. Nous recherchons :

$$R(\mathbf{x}), S(t) \text{ tels que } T^{m+1}(\mathbf{x}, t) = T^m + R(\mathbf{x})S(t) = \sum_{i=1}^{m} \varphi_i(\mathbf{x}) a_i(t) + R(\mathbf{x})S(t). \quad (3)$$

La formulation variationnelle de l'équation de la chaleur (Eq 2) s'écrit, pour tout champ $T^*(\mathbf{x}, t)$, comme suit :

$$\int_{\Omega_t} \int_{\Omega_e} \rho c_p^{app}(T^m) \frac{\partial T^{m+1}}{\partial t} T^* dV dt + \int_{\Omega_t} \int_{\Omega_e} k(T^m) \nabla T^{m+1} \nabla T^* dV dt = - \int_{\Omega_t} \int_{\Gamma_D} k(T^m) \nabla T_D \nabla T^* dV dt + \int_{\Omega_t} \int_{\Gamma_N} T^* h(T - T_{ext}) dS dt. \quad (4)$$

On notera que la non linéarité matériau est gérée de manière incrémentale [1], c'est à dire que les coefficients matériaux sont calculés à partir du résultat de l'itération précédente $m$.

Dans le cadre de la résolution du couple $R(\mathbf{x})$, $S(t)$ par la méthode du point fixe, nous obtenons deux équations. La première pour $T^*(\mathbf{x}, t) = R^*(\mathbf{x})S(t)$, permet de trouver $R(\mathbf{x})$ :

$$\int_{\Omega_t} \int_{\Omega_e} \rho c_p^{app}(T^m) R \frac{\partial S}{\partial t} R^* S dV dt + \int_{\Omega_t} \int_{\Omega_e} k(T^m) S^2 \nabla R \nabla R^* dV dt + \int_{\Omega_t} \int_{\Gamma_N} R^* R S^2 h dS dt = - \int_{\Omega_t} \int_{\Omega_e} k(T^m) S \nabla T_D \nabla R^* dV dt - \int_{\Omega_t} \int_{\Omega_e} \rho c_p^{app}(T^m) T^m R^* S dV dt - \int_{\Omega_t} \int_{\Omega_e} k(T^m) S T^m \nabla R^* dV dt + \int_{\Omega_t} \int_{\Gamma_N} h T^m R^* S dS dt + \int_{\Omega_t} \int_{\Gamma_N} h T_{ext} R^* S dS dt. \quad (5)$$

La seconde est pour $T^*(\mathbf{x}, t) = R(\mathbf{x})S^*(t)$ et permet de trouver $T^*(\mathbf{x}, t) = R(\mathbf{x})S^*(t)$ :

$$\int_{\Omega_t} \int_{\Omega_e} \rho c_p^{app}(T^m) R^2 \frac{\partial S}{\partial t} S^* dV dt + \int_{\Omega_t} \int_{\Omega_e} k(T^m) S S^* \nabla R \nabla R dV dt + \int_{\Omega_t} \int_{\Gamma_N} R^2 S S^* h dS dt = - \int_{\Omega_t} \int_{\Omega_e} k(T^m) S^* \nabla T_D \nabla R dV dt - \int_{\Omega_t} \int_{\Omega_e} \rho c_p^{app}(T^m) T^m R S^* dV dt - \int_{\Omega_t} \int_{\Omega_e} k(T^m) S^* T^m \nabla R dV dt + \int_{\Omega_t} \int_{\Gamma_N} h S^* T^m R dS dt + \int_{\Omega_t} \int_{\Gamma_N} h T_{ext} R S^* dS dt. \quad (6)$$

Ces deux équations sont résolues par une discrétisation EF spatio-temporelle dont les fonctions de formes sont respectivement $N_i(\mathbf{x})$ au nœud $i$ et $P_j(t)$ au piquet de temps $j$ ; par exemple $R(\mathbf{x}) = R^T N(\mathbf{x})$.

Dès lors, les coefficients matériaux s'écrivent $\rho c_p^{app}(T(\mathbf{x}, t)) = N(\mathbf{x})^T [C] P(t)$ et $k(T(\mathbf{x}, t)) = N(\mathbf{x})^T [K] P(t)$. Avec $[C]$ et $[K]$ des matrices de taille (nombre de degrés de liberté $x$ nombre de piquets de temps) qui représentent respectivement les propriétés matériaux de capacité calorifique ($\rho c_p$) et de conductivité matériaux ($k$). Cette écriture évite l'utilisation de la SVD telle que proposée par Joyot [5].

**Rq 1:** Les matrices $[C]$ et $[K]$ ne sont jamais assemblées totalement mais sont calculées uniquement au niveau élémentaire en chaque point de Gauss.

Ainsi, on peut dissocier les intégrales en temps et en espace de l'équation 5 et réaliser une séparation de variables sur $T^m$ (Eq 3) :
$$
\begin{align}
& R^{*T} \int_{\Omega_e} N N^T R N^T C \left[ \int_{\Omega_t} P P^T S S^T \frac{\partial P}{\partial t} dt \right] dV + R^{*T} \int_{\Omega_e} \nabla N \nabla N^T R N^T K \left[ \int_{\Omega_t} P S^T P P^T S dt \right] dV \nonumber \\
& + R^{*T} \int_{\Gamma_N} N N^T R h \left[ \int_{\Omega_t} P^T S S^T P dt \right] dS = - R^{*T} \int_{\Omega_e} \nabla N \nabla N^T T_D N K \left[ \int_{\Omega_t} P P^T S dt \right] dV \nonumber \\
& - R^{*T} \sum_i \int_{\Omega_e} N N^T \varphi_i N^T C \left[ \int_{\Omega_t} P P^T S a_i^T \frac{\partial P}{\partial t} dt \right] dV + R^{*T} \sum_i \int_{\Omega_e} \nabla N \nabla N^T \varphi_i N^T K \left[ \int_{\Omega_t} P a_i^T P P^T S dt \right] dV \tag{7} \\
& \qquad - R^{*T} \sum_i \int_{\Omega_e} N N^T \varphi_i h \left[ \int_{\Omega_t} S^T P P^T a_i dt \right] dS - R^{*T} \int_{\Gamma_N} N h T_{ext} \left[ \int_{\Omega_t} P^T S dt \right] dS. \nonumber
\end{align}
$$

Les intégrales temporelles (entre crochets) sont calculées tout d'abord à chaque piquet de temps.
L'assemblage spatial est réalisé ensuite.

Nous ferons de même pour l'équation 6. On obtient finalement un système linéaire d'une part en $R(\mathbf{x})$ :
$A_1 R = b_1$ et d'autre part en $S(t) : A_2 S = b_2$. La succession de ces deux systèmes jusqu'à convergence du
couple $R(\mathbf{x})$, $S(t)$ permet de calculer les nouveaux modes $\varphi_i^{m+1}(\mathbf{x}) = R(\mathbf{x})$ et $a_i^{m+1}(t) = S(t)$.

**Rq 2** : La convergence de cet algorithme, bien que théoriquement non acquise, ne pose en pratique pas de difficulté dans notre cas.

**Rq 3** : En raison de la non linéarité des matériaux, toutes les intégrales doivent être recalculées systématiquement à chaque itération du point fixe, contrairement à l’approche habituelle utilisée sur des problèmes linéaires [1].

3. Résultats

Une éprouvette à retassures 2D a été modélisée. Celle-ci est utilisée en métallurgie afin de provoquer des
retassures, c'est à dire des lacunes de matière, qui proviennent de la contraction du métal. La géométrie de
l'éprouvette présente une partie utile au dessous de la masselotte avec échancrure. Cette zone échancrée
solidifiera rapidement, tout comme les bords de l'éprouvette.

Deux maillages sont utilisés : le premier -1- (Fig 2 -a) présente 803 ddl spatiaux et le second -2- (Fig 2 -b) 2961 ddl. Ces maillages 2D sont constitués uniquement d'éléments Q4 linéaires. La simulation consiste à imposer une température initiale uniforme $T_o$ de 720°C dans toute l'éprouvette. Sur la totalité des bords sont imposés des flux convectifs avec une température extérieure $T_{ext}$ fixée à 100°C. Sur la surface supérieure (Fig 3 -b rouge gras), en contact avec l'air, est appliqué le coefficient d'échange convectif $h_1$ est égal à 200 W.m⁻².K⁻¹ (l'air est isolant). Sur les autres faces en contact avec le moule sable, le coefficient d'échange convectif $h_2$ est égal à 1000 W.m⁻².K⁻¹. Dans un premier temps, nous négligerons l'effet de la chaleur latente $L$, ceci revient à utiliser la capacité calorifique spécifique $c_p$.

<figure><img src="image_3.png" /><figcaption>**Fig 2** : Deux maillages 2D sont utilisés : le premier -1- (Fig 2 -a) présente 803 ddl spatiaux et le second -2- (Fig 2 -b) 2961 ddl.</figcaption></figure>
Afin d'obtenir une solution de référence, notée REF, ces modèles sont d'abord résolus au moyen d'un schéma de résolution standard (schéma implicite en temps, linéarisé par la méthode de Newton-Raphson et discrétisé spatialement en EF). Ensuite, la méthode PGD décrite ci-dessus est utilisée. Toutes les simulations sont réalisées sous MATLAB. La différence entre la solution de référence $T_{ref}(\mathbf{x}, t)$ et la solution PGD $T_{PGD}(\mathbf{x}, t)$ est caractérisée par l'erreur relative maximale sur tous les piquets de temps :

$$e = \max_t \frac{\| T_{ref}(\mathbf{x}, t) - T_{PGD}(\mathbf{x}, t) \|_{L2}}{\| T_{ref}(\mathbf{x}, t) \|_{L2}} \cdot 100. \qquad (8)$$

La figure 3 représente l'erreur *e* calculée pour chaque enrichissement *m*. Nous observons pour chacun des maillages (-1- et -2-)une stabilisation de l'erreur : à plus de 3 % après 5 modes pour le maillage -1- et à moins de 1 % après 11 modes pour le maillage -2-. Notons également que 3 enrichissements suffisent pour réduire sensiblement l'erreur.

<figure><img src="image_4.png" /><figcaption>Fig 3 : Erreur relative (%) en fonction du nombre d'enrichissements.</figcaption></figure>

Intéressons nous maintenant au temps de calcul. Nous pouvons observer dans le tableau 1 ci-dessous un très fort gain de temps, de l'ordre de 200 par rapport à la référence, pour une erreur relative acceptable à moins de 1 % pour le maillage -2-. Dans le cas d'une erreur relative acceptable de 3 %, le rapport de gain en temps est alors de plus de 2 000.

<table><caption>Tab 1 : Comparaison du temps de calcul des différents modèles.</caption><thead><tr><th></th><th>Maillage -1-<br>pour 5 modes</th><th>Maillage -2-<br>pour 3 modes</th><th>Maillage -2-<br>pour 11 modes</th></tr></thead><tbody><tr><td>Nb de ddl<br>spatiaux</td><td>803</td><td>2961</td><td>2961</td></tr><tr><td>REF</td><td>946 s</td><td>17 h</td><td>17 h</td></tr><tr><td>PGD</td><td>56 s</td><td>31 s</td><td>287 s</td></tr><tr><td>Rapport des<br>temps de calcul</td><td>15</td><td>2 000</td><td>200</td></tr><tr><td> Erreur relative<br>(%)</td><td>3,2 %</td><td>2,8 %</td><td>&lt; 1 %</td></tr></tbody></table>
Les figures 4 a-c représentent la fraction solide de l'éprouvette échancrée avec le maillage -2- et un seuil de coupure $f_s=0,7$. Pour $f_s \ge 0,7$, nous considérons la phase déjà solide et la représentons en blanc.

<figure><img src="image_5.png" /><figcaption>Fig 4 : Les figures 4 a – c représentent la fraction solide de l'éprouvette échancrée (à t = 14s) avec le maillage -2- (2961 ddl) et un seuil de coupure $f_s=0,7$. La couleur blanche représente une $f_s \ge 0,7$. La figure 4 -a représente notre solution de référence, la figure 4 -b, la solution PGD décrite avec 3 modes et la figure 4 -c, la solution PGD décrite avec 11 modes.</figcaption></figure>

Nous pouvons observer que la solution à 3 modes (Fig 4 -b) présente des différences sensibles avec la référence (Fig 4 -a). En effet, les deux poches encore liquides sont plus étroites. Malgré tout, la localisation des défauts paraît globalement très correcte. La solution à 11 modes (Fig 4 -c), quant à elle, ne présente pas de différence clairement identifiable par rapport à notre référence.

Pour affiner encore l'analyse de nos résultats, observons le comportement thermique de l'éprouvette au cours du temps par le biais de l'évolution temporelle d'un nœud particulier (Fig 2 -b). Nous observons (Fig 5) que l'évolution globale des trois courbes est très similaire sauf aux tout premiers instants. Le changement de phase localement prend très peu de temps. Aussi, dans notre cas, ce sont les 20 premières secondes qui déterminent la qualité de la solidification en ce point.

<figure><img src="image_6.png" /><figcaption>Fig 5 : Évolution de la température au cours du temps en un point de l'éprouvette.</figcaption></figure>
Un agrandissement des 25 premières secondes de la solidification (Fig 6) nous permet d'observer plus en détail les pentes à l'origine. Nous constatons à t = 0s que la solution à 3 modes présente un gradient thermique plus fort que celle de référence ou celle à 11 modes.

<figure><img src="image_7.png" /><figcaption>Fig 6 : Évolution de la température au cours du temps en un point de l'éprouvette.</figcaption></figure>

De plus, nous voyons plus clairement se dessiner l'erreur commise pour 3 ou 11 enrichissements. En terme de maximum d'erreur absolue, après 3 enrichissements : $\Delta T_{max} = \max_t (T_{REF} - T_{PGD}) = -18^\circ C$ à t = 5s alors qu'après 11 enrichissements : $\Delta T_{max} = -4^\circ C$ à t = 3s (soit 0,6%) !

## 4. Conclusion

La méthode PGD décrite ci-dessus fait appel à une formulation discrétisée par EF spatio-temporelle des propriétés matériaux. Dans le cadre de notre problème, la résolution obtenue est correcte. En effet, la méthode de réduction de modèles PGD décrite fournit des résultats de qualité (erreur relative < 1 %, erreur maximale absolue < 0,6 %) pour un faible nombre d'enrichissements (11 modes). Dès 3 enrichissements, la localisation des dernières poches liquides devient possible. De plus, le gain de temps est de l'ordre de 200 par rapport au temps de calcul standard. Ces travaux montrent l'efficacité de la PGD dans le cas de non linéarités des coefficients matériaux.

D'un point de vue industriel, il est souhaitable d'obtenir un pré-résultat permettant de décider s'il est pertinent de continuer une simulation qui prendra 24, 48 heures ou plus. En effet, s'il est possible de connaître rapidement l'évolution globale du champ de température au cours du temps, moyennant une marge d'erreur acceptable, et d'en déduire la localisation possible des défauts alors il devient possible d'optimiser le nombre de simulations nécessaire au cours du développement d'un nouveau produit. Dès lors, le gain de temps devient sensible et se matérialise sous forme de compétitivité.

Fort de ces résultats positifs, l'étape suivante consistera à prendre en compte l'effet thermique de la chaleur latente de fusion $L$. La capacité calorifique se trouvera nettement affectée par cette chaleur latente et sera d'autant plus non linéaire. Il sera certainement nécessaire d'adapter une nouvelle fois notre formulation afin que soit bien prise en compte cette énergie libérée lors de la solidification.
**5. Remerciements**

Ce travail est réalisé dans le cadre d'une thèse CIFRE Montupet SA – UTC pour le Projet 100 % Picard : ProFond. Ceci avec le soutien actif de la Picardie. Montupet SA (Laigneville) et E.J. (Méru) fournissent les cas tests industriels et les propriétés matériaux.

**6. Références bibliographiques**

[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 modeling of complex fluids, Journal of Non-Newtonian Fluid Mechanics, 139, 153-176, 2006.

[2] A. Ammar, M. Normandin, F. Chinesta, P. Villon. Non-incremental strategies based on separated representation : application in computational rheology, Commun. Math. Sci., Vol 8 N° 3, 671-695, 2010.

[3] J. A. Dantzig, M. Rappaz. Solidification, EPFL Press, 2009.

[4] Ch. Ghnatios, F. Chinesta, E. Cueto, A. Leygue, A. Poitou, P. Breitkopf, P. Villon. Methodological approach to efficient modeling and optimization of thermal processes taking place in a die: Application to pultrusion, Composites, 2011.

[5] P. Joyot, G. Bonithon, F. Chinesta, P. Villon. La méthode PGD-BEM appliquée à l'équation de la chaleur non linéaire, CSMA, 2011.

[6] Ngyen. La Décomposition Propre Généralisée pour la résolution de problèmes multiphysiques transitoires couplés dédiés à la mécanique des matériaux -Maillage adaptatif et couplage avec la MAN, ENSMA, 2012.

[7] Nouy. A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations, Comput. Methods Appl. Mech. Engrg., 2010.

[8] E. Pruliere, F. Chinesta, A. Ammar. On the deterministic solution of multidimensional parametric models using the PGD, Math. Comput. Simulat., 81-4, 791-810, 2010.