# Le schéma de Lax-Wendroff pour les solides élastoplastiques

Thomas Heuzé

► To cite this version:

Thomas Heuzé. Le schéma de Lax-Wendroff pour les solides élastoplastiques. 14ème Colloque National en Calcul de Structures (CSMA 2019), May 2019, Giens, France. http://hal-02279631

HAL Id: hal-02279631

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

Submitted on 5 Sep 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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Le schéma de Lax-Wendroff pour les solides élastoplastiques

T. Heuzé<fnref n="1" />

**Résumé** — On présente dans ce travail une extension du schéma de Lax-Wendroff pour le calcul en dynamique rapide de la réponse de solides élastoplastiques. Dans la formulation développée ici, la ré-solution est effectuée sur les composantes de la déformation et de la vitesse. Le schéma est ici présenté dans sa version à deux étapes (dite de Richtmyer), permettant de projeter deux fois par pas de temps une solution prédite élastiquement sur la surface de charge. Ceci permet de diminuer le nombre de cel-lules nécessaires pour l'approximation des ondes plastiques discontinues, par rapport à une projection simple. Des simulations sur un cas test d'onde plane présentant une solution analytique sont conduites et montrent des résultats satisfaisants.

**Mots clés** — Lax-Wendroff, Elastoplasticité

## 1 Introduction

Les problèmes d'impact sur des structures sont modélisés par des équations différentielles hyperbo- liques. Les solutions de ces dernières peuvent être discontinues, même à partir de champs initialement continus. En mécanique, les fronts d'onde générés par des impacts peuvent propager des phénomènes irréversibles (déformations plastiques) et des déformations qui peuvent être importantes. La capacité à capturer numériquement à la fois les solutions continues et discontinues de ces équations est essentielle si on souhaite comprendre le séquencement et la physique des événements. Par ailleurs, un calcul précis des déformations plastiques permettra une bonne évaluation de l'état résiduel de la structure impactée.

Le calcul en dynamique rapide est encore principalement réalisé avec une discrétisation spatiale par éléments finis couplée à un schéma de différence finies de type Newmark ou différences centrées en temps [1]. Bien que cette approche permette le calcul de solides dissipatifs via des intégrateurs numé- riques appropriés [2], les intégrateurs temporels classiques introduisent du bruit de haute fréquence au voisinage de discontinuités du champ qu'il est difficile de gommer au seul moyen de viscosité artificielle sans détruire la précision de la solution. Par ailleurs, une formulation classique en déplacement conduit à ce que le champ (le déplacement) et ses gradients convergent à des taux différents.

D'autres schémas basés sur les volumes finis [3] sont cependant de plus en plus utilisés pour le calcul de la réponse des solides dissipatifs en dynamique rapide [4, 5, 6]. Ces schémas sont en général construits de sorte à inclure la structure caractéristique des équations hyperboliques, soit de façon explicite via un solveur de Riemann définis aux interfaces entre cellules pour les schémas 'cell-centered', soit de façon implicite pour les schémas centrés. Ces schémas présentent certains intérêts, comme l'obtention d'un même taux de convergence pour les champs de vitesse et de contrainte avec une formulation en lois de conservation [5], ou bien, en fonction de la technologie retenue, la construction de solutions non- oscillantes tout en préservant un haut ordre de précision dans les régions plus régulières [3, 4].

Etant donné qu'aucun solveur de Riemann élastoplastique n'existe (même approximé) pour un état de contrainte multidimensionnel, le calcul de la réponse des solides élastoplastiques avec ces schémas est généralement effectué soit en projetant *a posteriori* sur la surface de charge les contraintes élastique- ment prédites [4], soit en effectuant une étape de relaxation au moyen d'un schéma étagé, correspondant plutôt à de la viscoplasticité [6, 7]. Cette dernière solution converge vers la solution élastoplastique en se donnant un paramètre de relaxation très faible, pour peu que le schéma soit capable de la capturer de façon asymptotique. Ces approches nécessitent généralement pas mal de cellules pour capturer les ondes plastiques quand elles sont discontinues, comme c'est par exemple le cas avec une cinématique d'onde plane et un écrouissage linéaire.

On s'intéresse dans ce travail au schéma de Lax-Wendroff [8] dans sa version de Richtmyer [9] à deux

<footnote n="1">¹ GeM, Ecole Centrale Nantes, thomas.heuze@ec-nantes.fr</footnote>
étapes de calcul au cours d'un pas de temps. D'une part, c'est un schéma centré [10] qui ne nécessite donc pas de solveur de Riemann. D'autre part, on peut tirer profit de cette version du schéma pour projeter deux fois au cours d'un pas de temps une solution élastiquement prédite sur le critère de plasticité. L'objectif est de diminuer le nombre de cellules nécessaires pour l'approximation des ondes plastiques discontinues, par rapport à des schémas n'effectuant qu'une seule projection par pas de temps. Des simulations sur un cas test d'onde plane présentant une solution analytique sont conduites et montrent des résultats satisfaisants. On notera que la version proposée ici diffère de celle introduite dans [11].

## 2 Problème aux valeurs limites et initiales élastoplastique

On considère un milieu continu régi, dans le cadre de l'hypothèse de transformation infinitésimale isotherme, par le système de lois de conservation suivant :

$$ \frac{\partial \mathcal{Q}}{\partial t} + \nabla \cdot \mathcal{F} = 0 \qquad (1) $$

où $\mathcal{Q}$ et $\mathcal{F}$ représente les vecteurs des quantités conservées, et les flux associés, définis respectivement comme :

$$ \mathcal{Q} = \begin{Bmatrix} \mathbf{v} \\ \boldsymbol{\varepsilon} \end{Bmatrix}; \quad \mathcal{F} = \begin{Bmatrix} -\frac{\boldsymbol{\sigma}}{\rho} \\ -(\mathbf{v} \otimes \mathbf{1} + \mathbf{1} \otimes \mathbf{v}) \end{Bmatrix} \qquad (2) $$

où $\boldsymbol{\varepsilon}, \boldsymbol{\sigma}$ sont les tenseur des déformations linéarisées et des contraintes de Cauchy, $\mathbf{v}$ est la vitesse et l'opérateur $\otimes$ est défini tel que $\mathbf{1} \otimes \mathbf{v} = \delta_{ik}v_j e_i \otimes e_j \otimes e_k$, $\delta_{ik}$ étant le symbol de Kronecker. La première équation du système (1) est le bilan de la quantité de mouvement, la seconde est la compatibilité géométrique. Les flux sont fonction d'un vecteur auxiliaire

$$ \mathcal{W} = \begin{Bmatrix} \mathbf{v} \\ \boldsymbol{\sigma} \end{Bmatrix} \qquad (3) $$

contenant les contraintes. La fermeture du problème est réalisée au moyen d'équations constitutives, qu'on particularise ici à un matériau élastoplastique à écrouissage isotrope linéaire, obéissant au critère de Mises :

$$ \dot{\boldsymbol{\varepsilon}} = \dot{\boldsymbol{\varepsilon}}^e + \dot{\boldsymbol{\varepsilon}}^p \qquad (4) $$

$$ \dot{\boldsymbol{\sigma}} = \mathbf{C} : \dot{\boldsymbol{\varepsilon}}^e \qquad (5) $$

$$ f = \boldsymbol{\sigma}_{\text{eq}} - C p - \boldsymbol{\sigma}_{\text{y}}; \quad \boldsymbol{\sigma}_{\text{eq}} = \sqrt{\frac{3}{2}} \mathbf{s} : \mathbf{s}; \quad \mathcal{C} = \{(\boldsymbol{\sigma}, \dot{\boldsymbol{\varepsilon}}^p) | f \le 0\} \qquad (6) $$

$$ \dot{\boldsymbol{\varepsilon}}^p = \frac{3}{2} \dot{p} \frac{\mathbf{s}}{\boldsymbol{\sigma}_{\text{eq}}} \qquad (7) $$

$$ \dot{p} \ge 0; f \le 0; \dot{p} f = 0 \qquad (8) $$

où les équations (4), (5), (6), (7) et (8) réfèrent respectivement à la partition additive du taux de déformation $(\dot{\boldsymbol{\varepsilon}})$ en parties élastique $(\dot{\boldsymbol{\varepsilon}}^e)$ et plastique $(\dot{\boldsymbol{\varepsilon}}^p)$, la loi d'élasticité, le critère $f$ associé à la norme de Mises et le convexe d'élasticité $\mathcal{C}$, la loi d'écoulement plastique et les conditions de complémentarité de Kühn-Tucker. On note $\mathbf{s}$ le déviateur des contraintes de Cauchy, $\dot{p}$ le multiplicateur plastique, $\mathbf{C}$ le tenseur des rigidités élastiques, et $\mathbf{C}$ le paramètre d'écrouissage. Le système (1) associé aux équations constitutives (4)-(8) sont de plus complétées par des conditions aux limites et initiales appropriées.

## 3 Le schéma de Lax-Wendroff

On considère un maillage constitué de de $N_c$ cellules polygonales, de noeuds sommets $p = 1, \dots, N_n$. Suivant les notations introduites dans [4, 11], on définit $\mathcal{P}(c)$ l'ensemble des noeuds connectés à la cellule $c$, $\mathcal{C}(p)$ celui des cellules connectées au noeud $p$, $\mathcal{E}(c)$ les arêtes connectées à une cellule $c$ et $\mathcal{E}^p(p)$ les arêtes connectées à un noeud $p$. Le maillage dual est lui constitué de polygones fermé autour d'un noeud $p$, dont les bords sont formés par des séparateurs $spc-$ et $spc+$, de longueurs $l_{spc-}$, $l_{spc+}$ et de normale unitaire $\mathbf{n}_{spc-}$, $\mathbf{n}_{spc+}$ respectivement, reliant le barycentre des cellules $\mathcal{C}(p)$ aux points milieux des arêtes $\mathcal{E}^p(p)$ (voir figure 1). Le chevauchement entre une cellule $c$ et une cellule duale centrée sur un noeud $p$ forme la sous-cellule $\mathbf{pc}$ de volume (ou plutôt d'aire ici) $V_{\mathbf{pc}}$.
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Cellule duale connectée au noeud p.</figcaption></figure>

## 3.1 Première étape de Lax-Wendroff

L'écriture des équations de conservation (1) intégrée sur une cellule duale de volume $V_p$ donne :

$$V_p \mathbf{U}_p^{n+\frac{1}{2}} = \sum_{c \in \mathcal{C}(p)} \left( V_{pc} \mathbf{U}_c^n - \frac{\Delta t}{2} \mathbf{F}_c^n \cdot (l_{spc-} \mathbf{n}_{spc-} + l_{spc+} \mathbf{n}_{spc+}) \right) \quad (9)$$

où $\mathbf{U}_p$ et $\mathbf{U}_c$ sont les moyennes intégrales du vecteur $\mathcal{Q}$ sur une cellule duale et primale respectivement. Les flux $\mathbf{F}_c^n$ sont évalués tels que $\mathbf{F}_c^n = \mathcal{F}(\mathbf{W}_c^n)$, $\mathbf{W}_c^n$ étant la moyenne intégrale du vecteur auxiliaire $\mathcal{W}$ (3) dans la cellule primale $c$ au temps $t_n$. Les déformations calculées au noeud $p$ sont utilisées pour actualiser les contraintes et les variables internes au milieu du pas de temps. Ces dernières étant connues dans les cellules $\mathcal{C}(p)$ connectées au noeud $p$ au temps $t_n$, on approxime leurs valeurs à ce noeud par une loi des mélanges :

$$(\varepsilon)^n_p = \frac{\sum_{c \in \mathcal{C}(p)} V_{pc} \varepsilon_c^n}{V_p}; \quad \sigma_p^n = \frac{\sum_{c \in \mathcal{C}(p)} V_{pc} \sigma_c^n}{V_p}; \quad (\varepsilon^p)^n_p = \frac{\sum_{c \in \mathcal{C}(p)} V_{pc} (\varepsilon^p)^n_c}{V_p} \quad (10)$$

Ensuite, les équations constitutives peuvent être intégrées avec un intégrateur numérique disponible [2], de sorte que la procédure d'actualisation s'écrive formellement :

$$\left\{ \sigma_p^{n+\frac{1}{2}}, (\varepsilon^p)^{n+\frac{1}{2}} \right\} = \mathcal{F} \left( (\varepsilon)^{n+\frac{1}{2}}_p, \sigma_p^n, (\varepsilon^p)^n_p \right) \quad (11)$$

où $\mathcal{F}$ représente l'intégrateur numérique.

## 3.2 Seconde étape de Lax-Wendroff

L'écriture des équations de conservation (1) intégrée sur une cellule $c$ du maillage de volume $V_c$ donne :

$$\mathbf{U}_c^{n+1} = \mathbf{U}_c^n - \frac{\Delta t}{V_c} \sum_{e \in \mathcal{E}(c)} (l_e \mathbf{F}_e^{n+1/2} \cdot \mathbf{n}_e^{(c)}) \quad (12)$$

où $l_e$ est la longueur de l'arête $e \in \mathcal{E}(c)$. Les flux $\mathbf{F}_e^{n+1/2}$ définis à ces arêtes sont calculés par moyennation de ceux définis aux noeuds sommets de ces arêtes $p$ et $p^+$, tels que :

$$\mathbf{F}_e^{n+1/2} = \frac{\mathbf{F}_p^{n+1/2} + \mathbf{F}_{p^+}^{n+1/2}}{2} \quad (13)$$

Les flux $\mathbf{F}_p^{n+1/2}$, $\mathbf{F}_{p^+}^{n+1/2}$ sont calculés en fonction des vitesses et des contraintes définis aux noeuds à l'instant $t_n + \Delta t/2$. Ensuite, l'actualisation des contraintes et des variables internes à l'instant $t_{n+1}$ se fait au moyen de l'intégrateur numérique de la loi de comportement élastoplastique, soit formellement :

$$\{ \sigma_c^{n+1}, (\varepsilon^p)^{n+1} \} = \mathcal{F} (\varepsilon_c^{n+1}, \sigma_c^n, (\varepsilon^p)^n_c) \quad (14)$$
## 4 Onde plane en élastoplasticité avec conditions initiales de type Riemann

On considère un domaine rectangulaire $(x, y) \in ]0, L[ \times ]0, H[, L \gg H$, dont le comportement satisfait aux équations (4)-(8). Les extrémités gauche et droite sont supposées libres et des conditions de symétrie sont imposées sur les bords haut et bas, de sorte qu'avec une hypothèse de déformations planes des ondes planes puissent être simulées. Des conditions initiales de type Riemann sont définies : la composante longitudinale de la vitesse est imposée à $-\bar{v}$ sur la première moitié du domaine $x \in ]0, L/2[$, et à $\bar{v}$ sur la seconde moitié. Les contraintes et déformations sont supposées initialement nulles partout. La vitesse imposée est telle qu'un écoulement plastique soit généré :

$$ \bar{v} = 2 \frac{Y_H}{\rho c_L} \qquad (15) $$

où $Y_H = (\lambda + 2\mu)\sigma_y/2\mu$ est la limité d'élasticité d'Hugoniot, $c_L$ est la célérité des ondes élastiques de pression $c_L = \sqrt{(\lambda + 2\mu)/\rho}$, et $\lambda, \mu$ et $\rho$ réfèrent aux paramètres de Lamé et à la masse volumique respectivement. Une solution analytique à ce problème a été fournie dans [12]. Ce cas test permet de générer une charge élastoplastique de compression suivant une première charge en tension, après réflexion et interactions des ondes élastiques et plastiques. Les valeurs numériques des paramètres utilisés pour les simulations sont résumées dans la table 1. Les figures 2 et 3 montrent une comparaison portant sur les

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Comparaison à l'instant $t = 4.37 \times 10^{-4}$ sec entre les solutions numériques et analytique.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Comparaison à l'instant $t = 4.37 \times 10^{-4}$ sec entre les solutions numériques et analytique.</figcaption></figure>

champs de contrainte et de déformation plastique à deux instants entre plusieurs solutions numériques, obtenues avec un nombre de Courant de 0.9 et un maillage de 200 × 3 éléments quadrangulaires. La solution obtenue avec le schéma de Lax-Wendroff (LW) est comparée avec celle obtenue par éléments finis (avec Cast3M [13]) avec des éléments Q1 et un intégrateur temporel explicite, une solution LW élastique
<table><caption>TABLE 1 – Valeurs numériques des paramètres</caption><thead><tr><th>Matériau</th><th></th><th>Géométrie</th></tr></thead><tbody><tr><td><i>E</i> = 200GPa</td><td><i>C</i> = 10GPa</td><td><i>L</i> = 6m</td></tr><tr><td><i>ν</i> = 0.3</td><td><i>σ</i><sub>y</sub> = 400MPa</td><td><i>H</i> = 0.08m</td></tr><tr><td><i>ρ</i> = 7800kg.m<sup>-3</sup></td><td></td><td></td></tr></tbody></table>

suivie d'une projection sur le critère de plasticité, et avec la solution analytique de ce problème. On observe que le schéma de Lax-Wendroff avec une double projection sur le critère améliore la résolution de l'onde plastique discontinue par rapport à une simple projection, et obtient des résultats similaires à ceux de la méthode des éléments finis. En revanche, les éléments finis et la double projection surestiment un peu le plateau de déformation plastique lors de la charge en tension, qui est un peu mieux capturé par la simple projection.

## 5 Conclusion

Dans ce papier on a présenté une extension du schéma de Lax-Wendroff à l'élastoplasticité, permettant dans son implémentation de Richtmyer de projeter une solution prédite élastiquement sur le critère de plasticité deux fois par pas de temps. Ceci permet de diminuer le nombre de cellules nécessaires pour l'approximation des ondes plastiques discontinues, même si la capture des plateaux de déformations plastiques n'est pas améliorée (sur cet exemple). Plus de détails sur le schéma proposé peuvent être trouvés dans [14].

## Références

[1] T. Belytschko, W.K. Liu, and B. Moran. *Nonlinear finite elements for continua and structures*. Wiley, 2000.

[2] J.C. Simo and T.J.R. Hughes. *Computational inelasticity*. Springer, 1997.

[3] R.J. Leveque. *Finite volume methods for hyperbolic problems*. Cambridge University Press, 2002.

[4] P.H. Maire, R. Abgrall, J. Breil, R. Loubère, and B. Rebourcet. A nominally second-order cell-centered lagrangian scheme for simulating elastic-plastic flows on two dimensional unstructured grids. *Journal of Computational Physics*, 235:626-665, 2013.

[5] J. Haider, C.H. Lee, A.J. Gil, and J. Bonet. A first order hyperbolic framework for large strain computational solid dynamics: An upwind cell centered total lagrangian scheme. *International Journal for Numerical Methods in Engineering*, 109(3):407-456, 2017.

[6] S. Ndanou, N. Favrie, and S. Gavrilyuk. Multi-solid and multi-fluid diffuse interface model: Applications to dynamic fracture and fragmentation. *Journal of Computational Physics*, 295:523-555, 2015.

[7] T. Heuzé. Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes. *Advanced Modeling and Simulation in Engineering Sciences*, 5(1):9, 2018.

[8] P. Lax and B. Wendroff. Systems of conservation laws. *Communications on Pure and Applied Mathematics*, 13(2):217-237, 1960.

[9] R.D. Richtmyer. A survey of difference methods for non-steady fluid dynamics. Technical report, National Center for Atmospheric Research, 1962. http://dx.doi.org/10.5065/D67P8WCQ.

[10] E.F. Toro and S.J. Billett. Centred TVD schemes for hyperbolic conservation laws. *IMA Journal of Numerical Analysis*, 20:47-79, 2000.

[11] D. Fridrich, R. Liska, and B. Wendroff. Cell-centered lagrangian Lax-Wendroff HLL hybrid method for elastic-plastic flows. *Computers and Fluids*, 157:164-174, 2017.

[12] T. Heuzé. Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids. *Journal of Computational Physics*, 346:369-388, 2017.

[13] Cast3M, 2018. User's manual. www.cast3M-cea.fr.

[14] T. Heuzé. Lax-Wendroff schemes for elastic-plastic solids. *Journal of Computational Physics*, 2019. submitted.