# POD et Méthode LATIN multiéchelle en temps et en espace

Jean-Charles Passieux, Pierre Ladevèze, David Néron

► To cite this version:

Jean-Charles Passieux, Pierre Ladevèze, David Néron. POD et Méthode LATIN multiéchelle en temps et en espace. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01422171

HAL Id: hal-01422171

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

Submitted on 23 Dec 2016

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 CC BY 4.0 - Attribution - International License
# POD et Méthode LATIN multiéchelle en temps et en espace

J.C. Passieux¹, P. Ladevèze¹,², D. Néron¹

¹ LMT (ENS de Cachan/CNRS/UPMC/PRES UniverSud Paris)
61, avenue du président Wilson, 94235 CACHAN CEDEX
{passieux,ladeveze,neron}@lmt.ens-cachan.fr

² Chaire de la fondation EADS « Advanced Computational Structural Mechanics »

**Résumé** — Ce travail concerne une technique d'approximation adaptative de type POD appliquée à la résolution d'une famille de problèmes par la méthode LATIN multiéchelle avec homogénéisation en temps et en espace. L'effet de synergie de ces deux méthodes procure à l'approche une grande efficacité dans le cadre des problèmes non linéaires d'évolution multiéchelles.

**Mots clés** — approximation radiale, séparation de variables, décomposition de domaine.

## Introduction

La simulation du comportement des matériaux complexes à des échelles fines conduisent à des problèmes à très grand nombre de degrés de liberté, dont le coût de calcul est généralement prohibitif avec les moyens de calculs actuels.

Dans [1], une nouvelle stratégie de calcul multiéchelle basée sur la méthode LATIN [2] a été proposée pour l'analyse de structures dont le comportement est décrit à une échelle fine, tant en temps qu'en espace. Cette stratégie de décomposition de domaine conduit à la résolution d'un grand nombre de problèmes d'évolutions indépendants par sous-structures à chaque itérations.

Si elle est menée avec une méthode incrémentale classique, la résolution de ces problèmes peut représenter un coût prohibitif. Dans ce cadre, une technique de résolution réduite, appelée approximation radiale, a été proposée dans [2, 3] et développée dans [4]. Il s'agit d'une méthode de construction *a priori* d'une décomposition de type POD [5] de la solution d'un problème d'évolution. Dans ce travail, nous proposons une reformulation de cette approximation dans le cadre de modèles de matériaux à variables internes [6]. Par ailleurs, au cours des itérations une partie de l'information générée est stockée dans une base réduite de fonctions. Sa réutilisation dans une phase préliminaire à chaque itération, sous forme d'un modèle réduit, est souvent suffisante à la construction d'une approximation de très bonne qualité. Dans le cas contraire, l'enrichissement de la base réduite est gérée par des indicateurs basés sur l'erreur en relation de comportement qui assurent à la stratégie de construire le minimum de fonctions nécessaires à la convergence de la méthode.

## 1 La méthode LATIN multiéchelle en temps et en espace

Dans cette partie, les grandes lignes de la méthode LATIN multiéchelle en temps et en espace sont présentées. Pour plus de détails, le lecteur pourra se référer à [2, 6]
# 1.1 Problème de décomposition de domaine

## 1.1.1 Problème de départ

Nous considérons l'évolution d'une structure visco-élastique définie par le domaine $[0, T] \times \Omega$. Dans le cas d'autre comportement matériaux, on pourra se reporter à [6]. La structure est soumise à des efforts extérieurs surfaciques $F_d$ sur son bord $\partial_2\Omega$. Sur la partie complémentaire du bord $\partial_1\Omega$, le déplacement $U_d$ est imposé (cf. Figure 1). Le tenseur des contraintes de Green-Lagrange

<figure><img src="image_2.png" /><figcaption>Figure 1 – problème avec décomposition de domaine</figcaption></figure>

est noté $\boldsymbol{\sigma}$. La déformation, notée $\boldsymbol{\varepsilon}$, est le gradient symétrique du déplacement $\underline{U}$. Elle se compose d'une partie élastique $\boldsymbol{\varepsilon}_e$ vérifiant la loi d'état $\boldsymbol{\sigma} = \mathbf{K}\boldsymbol{\varepsilon}_e$, et d'une partie anélastique $\boldsymbol{\varepsilon}_p$ suivant la relation de comportement $\dot{\boldsymbol{\varepsilon}}_p = \mathbf{B}(\boldsymbol{\sigma})$ éventuellement non linéaire. Afin de simplifier l'exposé, les conditions initiales sont supposées nulles, excepté le déplacement $\underline{U}_E|_{t=0} = \underline{U}_E^0$.

## 1.1.2 Décomposition du domaine

Le premier point de la méthode consiste d'une part à considérer la structure $\Omega$ comme un ensemble de sous-domaines $\Omega_E$ séparés par des interfaces $\Gamma_{EE'}$ (cf. Figure 1) et d'autre part, à diviser le domaine d'étude $[0, T]$ en un ensemble de sous-intervalles grossiers $I_i^G = [t_i^G, t_{i+1}^G]$. L'état de chaque sous-structure est défini par $\dot{\boldsymbol{\varepsilon}}_{pE}$ et $\boldsymbol{\sigma}_E$ (où $\square_E$ indique la restriction de $\square$ à $\Omega_E$). Chaque interface est caractérisée par $\underline{W}_{EE'}$, la restriction du déplacement $\underline{U}$ à $\Gamma_{EE'}$, et $\underline{F}_{EE'}$ la contrainte normale sur $\Gamma_{EE'}$. Une relation entre ces quantités est introduite pour traduire le comportement $\mathbf{b}_{EE'}$ de l'interface :

$$ \mathbf{b}_{EE'} (\dot{\underline{W}}_{EE'}, \dot{\underline{W}}_{E'E}, \underline{F}_{EE'}, \underline{F}_{E'E}) = 0 \quad (1) $$

Pour une interface parfaite, on aura (1) $\Leftrightarrow \dot{\underline{W}}_{EE'} = \dot{\underline{W}}_{E'E}$ et $\underline{F}_{EE'} + \underline{F}_{E'E} = 0$. On notera enfin $\underline{W}_E$ (resp. $\underline{F}_E$) l'ensemble des champs $\underline{W}_{EE'}$ (resp. $\underline{F}_{EE'}$) d'une même sous-structure $\Omega_E$.

## 1.1.3 Description à deux échelles

L'approche multiéchelle consiste à introduire une description à deux échelles des quantités d'interface $\underline{W}_E$ et $\underline{F}_E$. Une quantité $\underline{W}_E$ (resp. $\underline{F}_E$) est la somme de sa partie macro $\underline{W}_E^M \in \mathcal{W}_E^M$ (resp. $\underline{F}_E^M \in \mathcal{F}_E^M$) et de son complément micro $\underline{W}_E^m$ (resp. $\underline{F}_E^m$).

$$ \underline{F}_E = \underline{F}_E^M + \underline{F}_E^m \quad \text{avec} \quad \forall \underline{W}_E^* \in \mathcal{W}_E^M, \int_{I_i^G \times \Gamma} (\underline{F}_E - \underline{F}_E^M) \cdot \underline{W}_E^* ds = 0 \quad \text{et} \quad \underline{F}_E^m = \underline{F}_E - \underline{F}_E^M $$

Les espace $\mathcal{W}_E^M$ et $\mathcal{F}_E^M$ peuvent être choisis arbitrairement, mais sont pris égaux à l'espace produit des espaces vectoriels engendrés par une base de fonctions linéaires en espace $\mathcal{E}_i(x)$ et
quadratique en temps $\underline{f}_j(t)$ (cf. [1]). Ainsi :

$$ \forall \underline{F}_E^M \in W_E^M, \quad \underline{F}_E^M = \sum_i \sum_j \alpha_{ij} \underline{\varepsilon}_i(x) \underline{f}_i(t) $$

Le choix de l'espace d'admissibilité des efforts macro est un point majeur de la méthode et lui confère son caractère multiéchelle. En effet, on choisit d'imposer aux parties macro des efforts d'interface $\underline{F}^M = \{\underline{F}_E^M\}_E$ de vérifier systématiquement les conditions de transmissions aux interfaces. Le sous-espace de $\mathcal{F}_E^M$ correspondant est noté $\mathcal{F}_{E,ad}^M$. On défini $W_{E,ad}^M$ de manière analogue.

## 1.2 La stratégie de calcul multiéchelle avec homogénéisation

On définit $\mathbf{A}_d$, l'espace des solutions des équations linéaires, éventuellement globales et $\Gamma$ l'espace des solutions des équations locales, éventuellement non linéaires. La solution du problème sous-structuré est naturellement définie par $\mathbf{s}_{\text{exact}} = \mathbf{A}_d \cap \Gamma$. Ce qui revient à faire la séparation suivante :

$$ \mathbf{A}_d \quad \left| \begin{array}{ll} \text{- admissibilité cinématique} & \text{- loi d'évolution} \\ \text{- admissibilité statique} & \text{- comportement des interfaces} \\ \text{- loi d'état} & \text{- conditions initiales} \\ \text{- conditions initiales} & \text{- conditions initiales} \end{array} \right. \quad \text{et} \quad \Gamma $$

Dans ces conditions, on peut appliquer l'algorithmes LATIN [2], qui est un solveur itératif pour problèmes non linéaires d'évolution et qui a la particularité d'opérer sur l'espace-temps complet à chaque itération. Une itération de la méthode consiste à déterminer alternativement une solution $\hat{\mathbf{s}}_{n+1/2}$ de $\Gamma$ puis une solution $\mathbf{s}_{n+1}$ de $\mathbf{A}_d$ en suivant des directions de recherches $\mathbf{E}^+$ et $\mathbf{E}^-$. Les deux étapes correspondantes sont appelées *étape locale* et *étape linéaire*.

### 1.2.1 Étape locale

Cette étape consiste à déterminer $\hat{\mathbf{s}}_{n+1/2} \in \Gamma$, connaissant $\mathbf{s}_n \in \mathbf{A}_d$ en suivant la direction de recherche $\mathbf{E}^+$ définie par la relation suivante, $\forall \Omega_E$ :

$$ \mathbf{E}^+ \quad \left\{ \begin{array}{l} (\hat{\mathbf{e}}_{pE} - \hat{\mathbf{e}}_{pE}) + \mathbf{H}(\hat{\mathbf{\sigma}}_E - \mathbf{\sigma}_E) = 0 \\ (\hat{\mathbf{W}}_E - \hat{\mathbf{W}}_E) + \mathbf{h}(\hat{\mathbf{F}}_E - \mathbf{F}_E) = 0 \end{array} \right. \quad (2) $$

où $\mathbf{H}$ et $\mathbf{h}$ sont des opérateurs symétriques définis positifs qui sont des paramètres de la méthode. A cette étape, l'ensemble des équations étant locales en espace et en temps (si tant est que $\mathbf{B}$ et $\mathbf{b}$ le soient aussi), on doit résoudre un certain nombre d'équations non linéaires, mais locales en chaque point d'interpolation.

### 1.2.2 Étape linéaire

Ici, le problème consiste à trouver $\mathbf{s}_{n+1} \in \mathbf{A}_d$, connaissant $\hat{\mathbf{s}}_{n+1/2} \in \Gamma$ et en suivant la direction de recherche $\mathbf{E}^-$. À cause de la condition d'admissibilité des efforts macro la direction de recherche est globale en espace et en temps. Un multiplicateur de Lagrange $\dot{\underline{\mathbf{W}}}_E^M = (\dot{\underline{\mathbf{W}}}_E^M \in W_{E,ad}^{M*})_{\Omega_E \subset \Omega}$ est alors introduit pour garantir l'admissibilité des efforts macro au sens faible (voir [1]). La direction $\mathbf{E}^-$ s'écrit donc pour chaque $\Omega_E$ de la manière suivante :

$$ \mathbf{E}^{-} \quad \left\{ \begin{array}{l} (\hat{\mathbf{e}}_{pE} - \hat{\mathbf{e}}_{pE}) - \mathbf{H}(\hat{\mathbf{\sigma}}_E - \hat{\mathbf{\sigma}}_E) = 0 \\ (\hat{\mathbf{W}}_E - \hat{\mathbf{W}}_E) + \mathbf{h}(\underline{\mathbf{F}}_E - \underline{\mathbf{F}}_E) = \dot{\underline{\mathbf{W}}}_E^M \end{array} \right. \quad (3) $$

L'étape linéaire se résume donc à résoudre deux types de problèmes : un problème macro homogénéisé défini sur l'ensemble des inconnues macroscopiques d'interface ; et un problème micro indépendant par sous-structure tous deux définis sur l'intervalle grossier $I_i^G$.
## 2 Approximation radiale

À chaque itération, un ensemble de problèmes micro définis sur les sous-structures d'espace-temps $\Omega_E \times I_i^G$ doit donc être résolu : les problèmes variant peu d'une itération à l'autre. Les techniques incrémentales classiques ne tirent pas parti de cette propriété, et peuvent conduire à un coût de calcul prohibitif pour l'étape linéaire. Dans ce contexte, une technique de résolution alternative basée sur le concept de fonctions radiales (technique de résolution efficace en visco-plasticité [3]) a été proposée [2]. Elle consiste à rechercher la solution d'un problème d'évolution sous la forme d'un produits de m fonctions découplées :

$$f(\underline{x}, t) \approx \sum_{i=1}^{m} \lambda_{i}(t) \Lambda_{i}(\underline{x}) \qquad (4)$$

Dans lequel les fonctions radiales $\lambda_i$ et $\Lambda_i$ sont inconnues a priori. Cette écriture générale ap-partient à la famille des méthodes de type POD (pour Proper Orthogonal Decomposition, voir [5]), dont font partie les méthodes de décomposition a posteriori KLE (Expansion de Karhunen-Loeve), SVD (Décomposition en Valeurs Singulières), PCA (Analyse en Composantes Principales). Ici, il s'agit d'une méthode de résolution a priori qui a été utilisé avec succès dans le cadre du recalage, dans le domaine stochastique [7] ainsi que dans le cas de problèmes multi-dimensionnels [8].

### 2.1 Approximation des problèmes micro de la stratégie

**Problème de minimisation.** À l'itération $n$, la solution est recherchée sous la forme $\mathbf{s}_n = \mathbf{s}_{n-1} + \Delta\mathbf{s}_n$ où seule la correction de solution $\Delta\mathbf{s}_n \in \mathbf{A}_d^*$ est recherchée sous la forme radiale sus-citée. La direction de recherche (3) pouvant être interprétée comme une relation de comportement linéaire, les problèmes micro consistent, pour chaque sous-domaine $\Omega_E \times I_i^G$, à déterminer $\Delta\mathbf{s}_{n,E}$ qui minimise sur $\mathbf{A}_d^*$ l'erreur en relation de comportement $e_{RdC}^2(\mathbf{s}_E)$ associée à la direction de recherche (3), définie par :

$$e_{RdC}^2(\mathbf{s}_E) = \|\dot{\mathbf{e}}_{pE} - \mathbf{H}\mathbf{\sigma}_E + \bar{\mathbf{\Delta}}\|_{\mathbf{H},E}^2 + \|\underline{\dot{\mathbf{W}}}_E + \mathbf{h}\underline{\mathbf{F}}_E - \bar{\mathbf{\delta}}\|_{\mathbf{h},E}^2 \quad (5)$$

où $\bar{\mathbf{\Delta}}$ et $\bar{\mathbf{\delta}}$ ne dépendent que de $\mathbf{s}_{n-1}$ et $\hat{\mathbf{s}}_{n-1/2}$ connus à cette étape ; et avec les normes :

$$\|\square\|_{\mathbf{H},E}^2 = \int_{I_i^G \times \Omega_E} \left(1 - \frac{t}{T}\right) \square \circ \mathbf{H}^{-1} \square d\Omega dt \qquad \|\square\|_{\mathbf{h},E}^2 = \int_{I_i^G \times \partial\Omega_E} \left(1 - \frac{t}{T}\right) \square \cdot \mathbf{h}^{-1} \square dS dt$$

**Représentation des inconnues.** Il s'agit alors de représenter les inconnues minimales du problème sous forme radiale. Le choix qui est fait ici est une amélioration de la version proposée dans [4] dans le cas de modèles de matériaux à variables internes. Le premier point consiste à représenter la déformation inélastique et le déplacement d'interface de cette façon :

$$\Delta\varepsilon_p(t, \underline{\mathbf{M}}) = \sum_{k=1}^{p} \alpha^k(t) \mathbf{E}_p^k(\underline{\mathbf{M}}) \quad \text{et} \quad \Delta\underline{\mathbf{W}}(t, \underline{\mathbf{M}}) = \sum_{k=1}^{p} \alpha^k(t) \underline{\mathbf{Z}}^k(\underline{\mathbf{M}})$$

En utilisant les équations d'admissibilité, les autres variables $\Delta\varepsilon, \Delta\sigma, \Delta\underline{\mathbf{F}}$ peuvent être déduites des premières facilement. Selon nos expérimentations dans le cas viscoélastique, cette technique conduit approximativement au calcul de deux fois moins de fonctions temporelles.

**Algorithme.** Pour déterminer l'approximation d'ordre $m = 1$, l'idée consiste à minimiser alternativement par rapport aux fonctions spatiales et temporelles, les autres étant fixées. Ce qui conduit à l'alternance de deux types de problèmes :

- **Problème spatial** : $\alpha(t)$ étant fixée, on cherche $(\mathbf{E}_p, \underline{\mathbf{Z}})$ minimisant $e_{RdC}^2(\Delta\mathbf{s}_E)$.

- **Normalisation des fonctions spatiales**.
– *Problème temporel* : ($E_p, Z$) étant fixées, on cherche $\alpha(t)$ minimisant $e_{RdC}^2(\Delta s_E)$.

Le problème spatial est un problème éléments finis standard indépendant du temps. Le problème temporel, quant à lui, se résume à une simple équation différentielle ordinaire scalaire. En pratique cet algorithme converge rapidement. Cette méthode peut être vue comme une méthode des puisance généralisée sur un opérateur dépendant du problème et de son chargement. L'approximation d'ordre *m*, on construit les fonctions les unes après les autres par le même algorithme.

## 2.2 Une stratégie adaptative

Les fonctions spatiales étant, de loin, les plus coûteuses à calculer, elles seront stockées dans une base réduite. L'idée consiste à réutiliser au maximum l'information générée aux cours des itérations précédentes. On commencera donc chaque itération par un étape dite préliminaire dont le but est de trouver l'ensemble des fonctions temporelles $\alpha_k(t)$ associée à l'ensemble de fonctions spatiales de la base. Cette étape se traduit par la résolution d'un système d'équations différentielles du second ordre avec conditions initiales et finales.

À l'itération *n*, la résolution d'un problème micro se fait donc en deux étapes :

- **prédiction**. La résolution d'une étape préliminaire permet de déterminer une solution approchée pour un faible coût.

- **amélioration de la base**. Si la qualité de la prédiction n'est pas suffisante, une correction peut y être apportée par le calcul d'une nouvelle approximation d'ordre 1. La fonction générée est orthogonalisée par rapport aux fonctions de la base, puis si la norme du résidu n'est pas trop faible, elle est stockée dans la base réduite ; sinon, on conserve la prédiction et l'on passe à l'itération suivante.

<figure><img src="image_3.png" /><figcaption>Figure 2 – Influence de l'indicateur de qualité sur le comportement de la méthode et sa convergence. (a) indicateur LATIN par itérations LATIN ; (b) indicateur LATIN vs taille de la base réduite *m* ; nombre de nouvelle fonctions générées (c) et rejetée par le critère d'orthogonalité (d) par itérations</figcaption></figure>

*Choix de l'indicateur de qualité*. S'agissant d'un problème de minimisation, un indicateur naturel consiste à prendre la valeur de l'erreur $e_{RdC}^2(\Delta s_E)$ normée. On examine l'exemple d'une structure carrée 2D encastrée à sa base et soumise à un effort à son sommet. La structure est décomposée en 25 sous-domaines. La Fig.2 représente l'évolution de l'algorithme si l'on utilise cet indicateur. Sur 2(d)(gris), on constate que l'algorithme conduit au calcul d'un grand nombre de fonctions rejetées par le critère d'orthogonalité. On propose de modifier ce critère en le normant par le plus grand dénominateur sur l'ensemble des $\Omega_E$, ce qui permet de réduire considérablement le nombre de fonctions générées inutilement (cf. Fig.2 (d) (noir)), tout en ne dégradant pas la convergence globale de la stratégie (cf. Fig.2 (a) (trait noir discontinu)).

On peut remarquer de plus que le critère est bien adapté à la stratégie, puisque l'effort n'est mené que sur les sous-structures les plus mal convergées ; c'est-à-dire, ici, celles proches des conditions aux limites et celles proches des singularités géométriques (cf. Fig.3).
<figure><img src="image_4.png" /><figcaption>Figure 3 – Nombre de nouvelles fonctions par itération et par sous-structure</figcaption></figure>

## 3 Conclusion

Une véritable stratégie de calcul multiéchelle en temps et en espace adaptative est proposée. Elle consiste à mettre en synergie la méthode LATIN multiéchelle et la POD. L'approche est basée sur la réutilisation des fonctions issues de résolutions aux itérations précédentes par l'approximation radiale. Un indicateur basé sur l'erreur en relation de comportement a été mis en place afin de piloter la stratégie adapative et enrichir, quand il le faut, l'espace d'approximation des problèmes micro. Dès lors, seuls les calculs strictement nécessaires à la convergence de la méthode sont menés.

## Références

[1] P. Ladevèze and A. Nouy. On a multiscale computational strategy with time and space homogenization for structural mechanics. *Computer Methods in Applied Mechanics and Engineering*, 192:3061-3087, 2003.

[2] P. Ladevèze. *Nonlinear Computationnal Structural Mechanics - New Approaches and Non-Incremental Methods of Calculation*. Springer Verlag, 1999.

[3] P. Ladevèze. Sur une famille d'algorithmes en mécanique des structures. *Compte rendu de l'académie des Sciences*, 300(2):41-44, 1985.

[4] A. Nouy and P. Ladevèze. On a multiscale computational strategy with time and space homogenization for structural mechanics. *Computer Methods in Applied Mechanics and Engineering*, 192:3061-3087, 2004.

[5] Anindya Chatterjee. An introduction to the proper orthogonal decomposition. *Current Science*, 78(7):808-817, 2000.

[6] P. Ladevèze, D. Néron, and J.-C. Passieux. On multiscale computational mechanics with time-space homogenization. In *Multiscale methods—Bridging the scales in Science and Engineering*, chapter Space Time Scale Bridging methods. Oxford University Press. Ed. J. Fish, to appear 2009.

[7] A. Nouy. A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations. *Computer Methods in Applied Mechanics and Engineering*, 196(45-48):4521-4537, 2007.

[8] F. Chinesta, A. Ammar, F. Lemarchand, P. Beauchene, and F. Boust. Alleviating mesh constraints: Model reduction, parallel time integration and high resolution homogenization. *Computer Methods in Applied Mechanics and Engineering*, 197:400-413, 2008.