# Un schéma d'intégration temporelle multi-échelles appliqué à un problème poro-mécanique non-linéaire

Vincent Chiaruttini, Christian Rey

► To cite this version:

Vincent Chiaruttini, Christian Rey. Un schéma d'intégration temporelle multi-échelles appliqué à un problème poro-mécanique non-linéaire. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01812982

HAL Id: hal-01812982

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

Submitted on 12 Jun 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
# Un schéma d'intégration temporelle multi-échelles appliqué à un problème poro-mécanique non-linéaire

V. Chiaruttini*, Ch. Rey**

* Laboratoire de Modélisation Matériaux et Structures, UPMC/UCP
Case 161, 4 place Jussieu, F-75252 PARIS CEDEX 05
e-mail: chiarutt@ccr.jussieu.fr

** LMT Cachan, CNRS/ENS Cachan/UPMC
61 av. du Président Wilson, F-94235 CACHAN CEDEX
e-mail: rey@lmt.ens-cachan.fr

RÉSUMÉ. Une méthode efficace pour la résolution des problèmes mécaniques d'évolution non-linéaire est présentée. Basée sur un schéma d'intégration temporelle parallélisé, elle utilise une démarche multi-échelles dans un processus itératif pour atteindre des solutions fines aux problèmes de mécanique en évolution nécessitant des comportements non-linéaires complexes, avec un coût de calcul fortement réduit par rapport aux méthodes séquentielles. Les performances de cette approche seront illustrées par rapport aux méthodes classiques à travers plusieurs expérimentations numériques.

ABSTRACT. An efficient method to solve non-linear evolution mechanical problems is presented. It combines a parallelized temporal integration method with a multiscale strategy inside an iterative process in order to reach high-quality solutions to complex non-linear mechanical problems, with a highly reduced computational cost compared with sequential methods. The performance results of this approach will be highlighted and compared with traditional methods throughout several numerical assessements led on massive parallel computers.

MOTS-CLÉS : Parallélisme, intégration temporelle, multi-échelles, multi-physique, milieux poreux, mécanique non-linéaire des matériaux, plasticité.

KEYWORDS: Parallelism, temporal integration, multiscale, multiphysics, porous media, non-linear material mechanics, plasticity.
# 1. Introduction

La résolution efficace des systèmes d'équations différentielles non-linéaires représente un enjeu important en mécanique des structures. Pour répondre aux surcoûts créés par la complexité des comportements et la finesse des géométries actuelles, tous les aspects des processus de résolutions numériques tendent à être accélérés. Dans ce contexte, de nouvelles stratégies d'intégration temporelle (adaptées aux méthodes de sous-structuration [COM 02] ou aux phénomènes physiques [DUR 03, CHI 05]) ont été développées. Une autre démarche [LIO 01] propose une résolution parallélisée en temps avec des performances prometteuses mais un niveau d'accélération limité [FAR 03] eu égard aux méthodes de décomposition de domaine spatial.

L'objectif de notre étude est de proposer un algorithme de résolution parallélisé en temps permettant d'obtenir des gains significatifs en terme de performance numérique et des solutions de qualité a priori "illimitée", par la mise en œuvre d'une approche temporelle multi-échelles au sein d'une stratégie adaptée à la résolution des problèmes de mécanique non-linéaire à variables internes. L'ensemble de ces travaux est illustré par la résolution numérique de problèmes couplés en milieu poreux. Ces simulations présentent une dépendance temporelle forte due au transport de fluide ou de chaleur au sein de la structure et des comportements plastiques fortement non-linéaires, qui mettent à rude épreuve le processus de résolution numérique.

La section 2 présente la formulation d'un problème de poro-plasticité non-linéaire, la section 3 explicite l'algorithme de résolution parallélisée en temps, la section 4 illustre sa performance pour deux problèmes de poromécanique et la section 5 donne quelques conclusions et perspectives.

# 2. Formulation du problème de poro-plasticité non-linéaire

On considère l'évolution isotherme dans l'intervalle de temps $[0, t_f]$ d'une structure $\Omega$ formée d'un milieu poreux saturé, en petites déformations autour d'une position de référence [BIO 41]. Les hypothèses générales du modèle relèvent du comportement des Argilites de l'Est [CHA 04] : écoulement darcéen compressible, coefficient de Biot constant, loi de comportement élastoplastique avec un modèle de Drucker-Prager en contraintes effectives. On note $\vec{u}$ le déplacement et $\varepsilon$ le tenseur des déformations du squelette, $\sigma$ le tenseur des contraintes, $\sigma_d$ son déviateur, $p$ la pression du fluide, $\rho^{fl}$ la masse volumique du fluide, $\phi$ la porosité et $\vec{M}$ le flux de masse de fluide, et les notations usuelles correspondantes pour les conditions aux limites.

Deux équations de conservation gouvernent un milieu poreux en évolution isotherme, l'équilibre interne et la conservation de la masse de fluide (1) :

$$ \overrightarrow{\text{div}} \sigma + \rho \vec{g} = \vec{0} \qquad \frac{\partial \rho^{fl}}{\partial t} = -\text{div} \vec{M} \qquad (1) $$

où $\vec{g}$ représente l'accélération de la pesenteur et $\rho$ la masse volumique de la roche. L'écoulement est darcéen et compressible, ce qui conduit aux relations :

$$ \vec{M} = -k \left( \overrightarrow{\text{grad}} p - \rho^{fl} \vec{g} \right) \qquad \frac{d\rho^{fl}}{\rho^{fl}} = \frac{dp}{K_e} \qquad d\phi = (b - \phi) \left( \mathbf{I} : d\varepsilon + \frac{dp}{K_s} \right) \qquad (2) $$
où $k$ est la perméabilité du milieu, $K_e$ le module d'incompressibilité de l'eau, $b$ le module de Biot et $K_s$ le module d'incompressibilité des grains. Enfin, pour la loi de comportement mécanique, on fait l'hypothèse des contraintes effectives (en notant $\boldsymbol{\sigma}'$ le tenseur des contraintes effectives) et on introduit $\boldsymbol{\varepsilon}^p$ le tenseur des déformations plastiques :

$$ \boldsymbol{\sigma} = \boldsymbol{\sigma}' - b \boldsymbol{p} \mathbf{I} \qquad \boldsymbol{\sigma}' = \lambda \operatorname{Tr}(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^p) \mathbf{I} + 2 \mu_0 (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^p) \quad (3) $$

avec $\lambda$ et $\mu_0$ les coefficients de Lamé drainés. Le potentiel de dissipation plastique est régi par une fonction seuil $F$ du type Drücker-Pracker. Le modèle de plasticité est associé et les lois d'évolution dérivent normalement du potentiel.

$$ F = \sqrt{\frac{3}{2} \boldsymbol{\sigma}'_d : \boldsymbol{\sigma}'_d + \frac{2 \sin \Phi}{3 - \sin \Phi} \operatorname{Tr}(\boldsymbol{\sigma}') - \frac{6c \cos \Phi}{3 - \sin \Phi}} \quad (4) $$

où $\Phi$ représente l'angle de frottement, $c$ la cohésion de l'argilite.

Pour la formulation faible (5), les variables $\vec{u}$ et $p$ sont privilégiées. Les autres grandes sont, soit issues directement des relations de comportement ($y$ rassemble $\boldsymbol{\sigma}$, $\boldsymbol{\sigma}'$, $\vec{M}, \phi$), ou gérées comme des variables internes ($\nu$ contient $\rho^{fl}, \xi, \boldsymbol{\varepsilon}^p$) et vérifient localement la fonctionnelle $g$. Soit, $\forall t \in [0, t_f]$ :

$$ \int_{\Omega} \boldsymbol{\sigma} : \boldsymbol{\varepsilon}^* d\Omega = \int_{\Omega} \rho \vec{g} \cdot \vec{u}^* d\Omega + \int_{\partial_F \Omega} \vec{F}_d \cdot \vec{u}^* dS, \forall \vec{u}^* \in \mathcal{U}_0^{ad} $$

$$ \int_{\Omega} \vec{M} \cdot \text{grad} p^* d\Omega - \int_{\Omega} \frac{\partial \rho^{fl} \phi}{\partial t} p^* d\Omega = \int_{\partial_M \Omega} M_d p^* dS, \forall p^* \in \mathcal{P}_0^{ad} \quad (5) $$

$$ g(\nu, y, x, t, P) = 0, \forall P \in \Omega $$

$$ \vec{u}(P, t = 0) = \vec{u}_0(P), \quad p(P, t = 0) = p_0(P), \quad \nu(P, t = 0) = \nu_0(P), \quad \forall P \in \Omega $$

### 3. Mise en œuvre d'un schéma d'intégration temporelle multi-échelles pour les problèmes de mécanique à variables internes

À l'instar de l'approche décrite dans [LIO 01], la résolution du système différentiel non-linéaire (5) est menée par un algorithme d'intégration temporelle parallélisé à plusieurs échelles. L'intervalle temporel $[t_0, T_f]$ est décomposé en $N$ domaines : $t_0 = T^0 < T^1 < \dots < T^i < \dots < T^N = t_f$. Sur chaque intervalle $[T^i, T^{i+1}]$, le problème se réécrit $\forall i \in \{0, 1, \dots, N-1\}$ :

$$ \left\{ \begin{array}{l} f\left(\frac{dx^i}{dt}(t), x^i(t), t\right) = 0, \forall t \in [T_i, T_{i+1}] \\ g(\nu^i, y^i, x^i, t, P) = 0, \forall P \in \Omega, \forall t \in [T_i, T_{i+1}] \\ x^i(T^i) = h^i, \nu^i(T^i) = \mu^i \end{array} \right. \quad (6) $$

$x$ contient les inconnues principales ($\vec{u}$ et $p$), $\nu$ les variables internes, $h^i$ et $\mu^i$ les conditions initiales du sous-problème. Cette écriture est équivalente au problème original si et seulement si la connexion est assurée exactement entre chaque domaine :

$$ h^i = x^{i-1}(T^i), \quad \mu^i = \nu^{i-1}(T^i), \quad h^0 = x_0, \quad \mu^0 = \nu_0 \quad (7) $$
L'algorithme parallèle en temps repose sur la construction d'un processus itératif

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Grilles des intégrateurs temporels sur le domaine *i*.</figcaption></figure>

permettant de corriger la valeur initiale $(h^i, \mu^i)$ de chaque intervalle et de mener des résolutions en parallèle sur chaque domaine. La stratégie multi-échelles retenue utilise plusieurs intégrateurs $\mathcal{F}_j^i$, avec des discrétisations spécifiques (figure 1), tel que $(x, v)(t) = \mathcal{F}_j^i(h^i, \mu^i, t)$ représente une solution du problème (6) obtenue par un calcul sur une grille temporelle de pas de temps $\delta t_j^i$. Pour minimiser le temps de calcul et atteindre des solutions de plus en plus précises, à chaque nouvelle itération un intégrateur disposant d'une grille temporelle plus fine est utilisé.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Intégration des variables principales et du comportement sur le domaine *i*.</figcaption></figure>

La première étape consiste à initialiser les conditions initiales sur chaque intervalle en menant une résolution avec l'intégrateur $\mathcal{F}_0^i$ :

$$ egin{cases} (h_0^0, \mu_0^0) = (X_0^0, V_0^0) = (x_0, v_0) \\ (h_0^{i+1}, \mu_0^{i+1}) = (X_0^{i+1}, V_0^{i+1}) = \mathcal{F}_0^i(h_0^i, \mu_0^i, T^{i+1}), \forall i \in \{1, \dots, N-2\} \end{cases} \quad (8) $$

Puis on effectue un processus itératif à 4 étapes sur l'indice *j* :

1) Sur chaque domaine, le problème est résolu *en parallèle* avec l'intégrateur $\mathcal{F}_j^i$ spécifique à l'itération *j* : $(x_j^{i+1}, v_j^{i+1}) = \mathcal{F}_j^i(h_j^i, \mu_j^i, T^{i+1})$.

2) Pour améliorer la convergence, une correction est apportée à l'aide de l'intégrateur initial $\mathcal{F}_0^i$. Sur le premier domaine $i = 0$, on pose $(h_{j+1}^1, \mu_{j+1}^1) = (x_j^1, v_j^1)$, pour les domaines suivants, on résout puis corrige les conditions initiales des variables principales :

$$ (X_{j+1}^{i+1}, V_{j+1}^{i+1}) = \mathcal{F}_0^i(h_{j+1}^i, \mu_{j+1}^i, T^{i+1}) \qquad h_{j+1}^{i+1} = x_j^{i+1} + (X_{j+1}^{i+1} - X_j^{i+1}) \quad (9) $$

L'admissibilité des variables internes vis-à-vis des variables principales corrigées est obtenue par la résolution de la fonctionnelle *g* seule (figure 2) :

$$ \mu_{j+1}^{i+1} = \nu \text{ vérifiant : } \begin{cases} g(\nu, y, h_{j+1}^{i+1}, T^{i+1}, P) = 0, & \forall P \in \Omega \\ \nu(T^i) = \mu_{j+1}^i \end{cases} \quad (10) $$
3) La convergence de la méthode s'évalue à partir de la distance entre deux solutions successives sur l'interface entre domaines :

$$ \eta_j^x = \max_i ||h_{j+1}^i - h_j^i||, \quad \eta_j^\nu = \max_i ||\mu_{j+1}^i - \nu_j^i|| \quad (11) $$

4) La finesse de l'intégrateur pour l'itération suivante est déterminée en fonction du gain en précision entre 2 itérations successives : $\delta t_{j+1}^i = \delta t_j^i \eta_j^x / \eta_{j-1}^x$.

Quand la convergence est atteinte, un post-processing, mené en parallèle sur chaque domaine, effectue une correction linéaire de la solution afin d'assurer la continuité et la connexion avec les dernières conditions initiales corrigées. Cette étape améliore la qualité du résultat de façon comparable au calcul d'une nouvelle solution fine. La mise en œuvre numérique est obtenue avec des schémas du premier ordre (θ-méthode) pour les intégrations temporelles et une méthode de type Newton-Raphson pour les résolutions non-linéaires.

En termes de performances, le processus multi-échelles permet de diminuer considérablement le coût des premières itérations et rend la précision à convergence "illimitée". Si chaque itération apporte un gain de précision constant de facteur $k$ sur une grille à $k^j$ pas de temps, l'accélération théorique en temps de calcul par rapport à la résolution séquentielle de précision équivalente pour $J_0$ itérations s'écrit :

$$ S = \frac{N k^{J_0}}{N + J_0 + k \frac{k^{J_0-1}}{k-1}} \le N \frac{k-1}{k}, \text{ avec } N \text{ processeurs et } k \ge 2 $$

## 4. Mise en œuvre numérique

Pour évaluer la pertinence de la méthode, nous avons analysé la résolution de deux problèmes tests sur un calculateur parallèle (SMP NEC TX-7 à 32 processeurs) : une simulation de consolidation 3D en poro-élasticité, puis la simulation 2D d'une excava- tion en poroplasticité. Pour ces applications nous avons utilisé les paramètres matériau fournis dans le Benchmark Hydromécanique du GDR MOMAS [CHA 04].

En poroélasticité linéaire, l'algorithme présente des propriétés remarquables : la qua-

<table><caption>**Tableau 1.** Performances des méthodes d'intégration numérique en poroélasticité.</caption><thead><tr><th rowspan="2">Méthode</th><th rowspan="2">Itérations/<br>Intégrations</th><th colspan="2">Précision temporelle</th><th rowspan="2">Temps de<br>calcul</th></tr><tr><th>U</th><th>P</th></tr></thead><tbody><tr><td>Séquentielle</td><td>1372</td><td>1.35e-04</td><td>7.52e-04</td><td>T<sub>1</sub> = 5746s</td></tr><tr><td>Pararéelle 16 CPU</td><td>2/3456</td><td>1.64e-04</td><td>1.53e-03</td><td>T<sub>2</sub> = 2062s</td></tr><tr><td>Pararéelle 16 CPU</td><td>3/3456</td><td>5.59e-05</td><td>4.33e-04</td><td>T<sub>3</sub> = 3055s</td></tr><tr><td>Séquentielle</td><td>3456</td><td>3.47e-05</td><td>2.20e-04</td><td>T<sub>4</sub> = 14557s</td></tr><tr><td>Multi-échelles 16 CPU</td><td>3/3456</td><td>3.57e-05</td><td>2.33e-04</td><td>T<sub>5</sub> = 1281s</td></tr></tbody></table>

<table><caption>**Tableau 2.** Performances des méthodes d'intégration numérique en poro-plasticité.</caption><thead><tr><th rowspan="2">Méthode</th><th rowspan="2">Itérations/<br>Intégrations</th><th colspan="2">Précision temporelle</th><th rowspan="2">Temps de<br>calcul</th></tr><tr><th>U</th><th>P</th></tr></thead><tbody><tr><td>Séquentielle</td><td>512</td><td>4.03e-04</td><td>2.51e-03</td><td>T<sub>6</sub> = 2734s</td></tr><tr><td>Multi-échelles 8 CPU</td><td>3/512</td><td>4.80e-04</td><td>2.72e-03</td><td>T<sub>7</sub> = 617s</td></tr></tbody></table>
lité des solutions est équivalente aux résolutions séquentielles avec la grille de la der-
nière itération. Alors que la méthode "Pararéelle" à 16 domaines fournit des accé-
lérations limitées (entre $2,8 = T_1/T_2$ et $4,7 = T_4/T_3$) par rapport à la résolution
séquentielle de qualité équivalente, notre approche avec 16 processeurs fournit une
accélération de facteur $11,4 = T_4/T_5$.

Pour le problème fortement non-linéaire, notre approche constitue une première ten-
tative de parallélisation en temps d'un problème de mécanique complexe, le niveau
d'accélération obtenu est de $4,4 = T_6/T_7$ pour 8 processeurs.

**5. Conclusion**

Nous avons présenté une stratégie de résolution basée sur méthode intégration pa-
rallélisée en temps utilisant une approche itérative multi-échelles pour les problèmes
de mécanique non-linéaire. Les résultats obtenus sur des problèmes linéaires montrent
une efficacité indéniable de la méthode pour atteindre des solutions de grande préci-
sion à moindre coût. Pour le problème fortement non-linéaire à variables internes, la
démarche présentée procure un gain en temps de calcul et des évolutions sont à l'étude
en associant des accélérateurs de type multi-résolutions [GOS 02].

**Remerciements.** Nous remercions chaleureusement F. Feyel pour l'aide apportée à la mise en œuvre de nos approches au sein du code ZeBuLoN [BES 97]. Nous remercions également le pôle MESO Ile-de-France-Sud pour les moyens de calcul mis à notre disposition.

**6. Bibliographie**

[BES 97] BESSON J., FOERCH R., « Large Scale Object-Oriented Finite Element Code Design », *Comp. Meth. in Appl. Mech. and Eng.*, vol. 142, 1997, p. 165-187.

[BIO 41] BIOT M. A., « General theory of three dimensional consolidation », *Journal of Applied Physics*, vol. 12, n° 2, 1941, p. 188-164.

[CHA 04] CHAVANT C., MILLARD A., « Simulation d'excavation en comportement hydro mécanique fragile », *Benchmark proposé dans le cadre du GDR MOMAS*, vol. 4, 2004.

[CHI 05] CHIARUTTINI V., REY C., « A self-adaptive time-multiscale algorithm for coupled evolution problems », *Coupled Problems*, 2005.

[COM 02] COMBESCURE A., GRAVOUIL A., « A numerical scheme to couple subdomains with different time-steps for predominantly linear transient analysis », *Computer Methods in Applied Mechanics and Engineering*, vol. 191, 2002, p. 1129-1157.

[DUR 03] DUREISSEX D., LADEVÈZE P., NERON D., SCHREFLER B., « A multi-time-scale strategy for multiphysics problems : application to poroelasticity », *Int. Jal for Multiscale Computational Engineering*, vol. 1, n° 4, 2003, p. 387-400.

[FAR 03] FARHAT C., CHANDESRIS M., « Time-Decomposed Parallel Time-Integrators : Theory and Feasibility Studies for Fluid, Structure, and Fluid-Structure Applications », *International Journal for Numerical Methods in Engineering*, vol. 58, 2003, p. 1397-1434.

[GOS 02] GOSSELET P., REY C., « On a selective reuse of Krylov subspaces in Newton-Krylov approaches for nonlinear elasticity », *DD 14*, 2002, p. 419-426.

[LIO 01] LIONS J.-L., MADAY Y., TURINICI G., « Résolution d'EDP par un schéma en temps pararéel », *C. R. Acad. Sci. Paris*, vol. 333, n° 1, 2001, p. 1-6.