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

HAL
open science

# La Méthode des Points Matériels en Galerkin Discontinue pour la simulation d'impact sur des solides

Adrien Renaud, Thomas Heuzé

► To cite this version:

Adrien Renaud, Thomas Heuzé. La Méthode des Points Matériels en Galerkin Discontinue pour la simulation d'impact sur des solides. 13e colloque national en calcul des structures, Université Paris-Saclay, May 2017, Giens, Var, France. http://hal-01899330

HAL Id: hal-01899330

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

Submitted on 19 Oct 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# La Méthode des Points Matériels (MPM) en Galerkin Discontinue (DG) pour la simulation d'impact sur des solides.

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

**Résumé** — La Méthode des Points Matériels est étendue dans ce travail au carde de la méthode de Galerkin Discontinue pour la simulation d'impact sur des solides élastiques. La formulation est basée sur l'affaiblissement de la forme conservative des équations de la dynamique et sur le calcul des flux de Godunov aux interfaces entre les cellules d'une grille de calcul eulérienne. La méthode est ici illustrée sur deux cas unidimensionnels sur lesquels elle est comparée avec la MPM, et sur un cas d'impact bidimensionnel.

**Mots clés** — Méthode des Points Matériels, méthode de Galerkin Discontinue, impact, solide élastique.

## 1 Introduction

Les procédés de mise en forme à grande vitesse des matériaux métalliques comme le formage électromagnétique font intervenir des grandes transformations et de grandes vitesses de déformation (10³ – 10⁴ s⁻<fnref n="1" />). Le chargement rapide implique la présence de fronts d'onde qu'il est important de capter pour comprendre les phénomènes physiques présents, bien estimer les déformations et contraintes résiduelles, et optimiser le procédé. Cependant, les grandes transformations rendent complexe la construction d'un schéma numérique adapté. L'objectif de ce travail est de construire une méthode numérique permettant un suivi précis des fronts d'onde dans un milieu subissant de grandes transformations avec une description lagrangienne. La méthode des Éléments Finis Galerkin Discontinue (DGFEM) et la Méthodes des Points Matériels (MPM) remplissent chacune un de ces objectifs, qu'on propose ici de coupler.

La MPM [1], développée dans les années 90, est une adaptation des méthodes de type *Particle In Cell* (PIC) [2] à la mécanique du solide dans laquelle une collection de points matériels - ou particules - est mobile dans une grille de calcul eulérienne - ou maillage - fixe ou non. Les points matériels transportent avec eux l'intégralité des champs du problème et donc son histoire. La résolution du système numérique est faite sur la grille, ce qui implique des étapes de transfert entre les particules et les nœuds du maillage. Ainsi, la méthode permet de s'affranchir du calcul de termes de convection puisque l'information est transportée par les points matériels.

La DGFEM a été développée dans les années 70 pour le transport des neutrons [3]. Depuis, elle a été appliquée à la mécanique des fluides et des solides pour la résolution de problèmes elliptiques et hyperboliques. Contrairement à une approche *éléments finis* (FEM) classique, la DGFEM utilise des fonctions de forme discontinues à l'interface entre deux éléments [4]. L'interaction entre les éléments à travers les interfaces qu'ils partagent est prise en compte par l'intermédiaire de flux numériques qui peuvent être de plusieurs natures.

On développe en premier abord la DGMPM à partir d'une forme faible qui diffère de celles de la MPM et de la DGPIC [5], et qui fait apparaître entre les éléments d'un maillage, des flux que l'on calcule par la méthode de Godunov après résolution d'un problème de Riemann. On teste ensuite cette méthode en 1D et en 2D sur un solide élastique pour enfin discuter des perspectives de ce premier travail.

<footnote n="1">¹ GeM, Ecole Centrale Nantes, {adrien.renaud,thomas.heuze}@ec-nantes.fr</footnote>
## 2 Problème aux valeurs limites et initiales élastique

### 2.1 Système hyperbolique dans le cas de l'élasticité linéaire en petites perturbations

Soit le système d'équations suivant régissant le comportement dynamique d'un système élastique :

$$ \left. \begin{array}{l} \rho \frac{\partial \mathbf{v}}{\partial t} - \nabla \cdot \boldsymbol{\sigma} = \mathbf{0} \\ \frac{\partial \boldsymbol{\sigma}}{\partial t} - \mathbb{C} : \nabla_s \mathbf{v} = \mathbf{0} \end{array} \right\} \quad \forall \mathbf{x} \in \Omega, \forall t \in [0, t_{fin}] $$

où $\boldsymbol{\sigma}$ et $\mathbf{v}$ représentent le tenseur des contraintes de Cauchy et le vecteur vitesse, $\rho$ et $\mathbb{C}$ la masse volumique et le tenseur d'élasticité supposés constants par rapport au temps et aux coordonnées de l'espace. Le système d'équations précédent peut aussi s'écrire sous la forme conservative :

$$ \frac{\partial \mathbf{q}}{\partial t} + \sum_l \frac{\partial \mathbf{F}_l(\mathbf{q})}{\partial x_l} = \mathbf{0} \quad \forall \mathbf{x} \in \Omega, \forall t \in [0, t_{fin}] \quad (1) $$

où $\mathbf{q}$ est appelé le vecteur des quantités conservées $\mathbf{q}^T = \langle \boldsymbol{\sigma}, \rho \mathbf{v} \rangle$. La forme quasi-linéaire associée s'écrit :

$$ \frac{\partial \mathbf{q}}{\partial t} + \sum_l \mathbf{a}_l \cdot \frac{\partial \mathbf{q}}{\partial x_l} = \mathbf{0} \quad \forall \mathbf{x} \in \Omega, \forall t \in [0, t_{fin}] $$

avec $\mathbf{a}_l$ la jacobienne définie par :

$$ \mathbf{a}_l = \frac{\partial \mathbf{F}_l}{\partial \mathbf{q}} $$

Pour une barre unidimensionnelle, le vecteur flux prend la forme :

$$ \mathbf{F}(\mathbf{q}) = - \begin{bmatrix} E_v \\ \sigma \end{bmatrix} \quad \text{et} \quad \mathbf{q} = \begin{bmatrix} \boldsymbol{\sigma} \\ \rho v \end{bmatrix} $$

La généralisation au cas multidimensionnel est donnée dans [6]. On utilisera dans la suite des grandeurs spécifiques, c'est à dire par unité de masse :

$$ \mathbf{q}^s = \frac{\mathbf{q}}{\rho} \qquad \mathbf{F}_l^s = \frac{\mathbf{F}_l}{\rho} \qquad (2) $$

### 2.2 Méthode des Points Matériels [1] étendue au cadre Galerkin Discontinu

Soit une grille de calcul eulérienne composée de $N$ éléments, ou cellules, telle que $\Omega = \cup_{e=1}^N \Omega_e$. Le formalisme DG utilise une approximation discontinue aux interfaces entre les éléments qui conduit à une forme faible de l'équation (1) sur chacune des cellules :

$$ \int_{\Omega_e} \frac{\partial \mathbf{q}}{\partial t} \cdot \mathbf{w} \, d\Omega + \int_{\Omega_e} \sum_l \frac{\partial \mathbf{F}_l}{\partial x_l} \cdot \mathbf{w} \, d\Omega = 0 \quad \forall \mathbf{w} $$

et après intégration par partie, il vient :

$$ \int_{\Omega_e} \frac{\partial \mathbf{q}}{\partial t} \cdot \mathbf{w} \, d\Omega - \int_{\Omega_e} \sum_l \mathbf{F}_l \cdot \frac{\partial \mathbf{w}}{\partial x_l} \, d\Omega + \int_{\Gamma_e} \mathbf{F}_n \cdot \mathbf{w} \, d\Gamma = 0 \quad \forall \mathbf{w} \quad (3) $$

où $\mathbf{F}_n$ désigne le flux normal. Sur la grille, la masse volumique est approximée en fonction des masses portées par les $N_p$ points matériels :

$$ \rho(x) = \sum_{p=1}^{N_p} m_p \delta(x - x_p) \qquad (4) $$

où $\delta$ est la distribution dirac. La forme faible (3) écrite avec les grandeurs spécifiques (2) conduit à :

$$ \int_{\Omega} \rho \frac{\partial \mathbf{q}^s}{\partial t} \cdot \mathbf{w} \, d\Omega - \int_{\Omega} \rho \sum_l \mathbf{F}_l^s \cdot \frac{\partial \mathbf{w}}{\partial x_l} \, d\Omega + \int_{\Gamma_e} \mathbf{F}_n \cdot \mathbf{w} \, d\Gamma = 0 \quad \forall \mathbf{w} $$
et, en injectant la relation (4), les intégrales sur $\Omega_e$ sont transformées en des sommes finies sur les points matériels appartenant à $\Omega_e$ :

$$ \sum_{p=1}^{N_p} m_p \frac{\partial \mathbf{q}^s(x_p)}{\partial t} \cdot \mathbf{w}(x_p) - \sum_{p=1}^{N_p} m_p \sum_l \mathbf{F}_l(x_p) \cdot \frac{\partial \mathbf{w}(x_p)}{\partial x_l} + \int_{\Gamma_e} \mathbf{F}_n \cdot \mathbf{w} \, d\Gamma = 0 \quad \forall \mathbf{w} \quad (5) $$

## 2.3 Forme semi-discrète

La forme faible (5) fait appel à des champs évalués aux points matériels, qui sont exprimés par interpolation des quantités nodales grâce aux fonctions de forme $S(x)$. En injectant l'interpolation des champs dans (5) il vient la forme semi-discrète du problème :

$$ \begin{gathered} \sum_{p=1}^{N_p} S_{ip} m_p S_{jp} \frac{\partial \mathbf{q}_j^s}{\partial t} = \sum_l \sum_{p=1}^{N_p} \frac{\partial S_{ip}}{\partial x_l} m_p S_{jp} \mathbf{F}_{l,j}^s - \int_{\Gamma_e} S_i(\mathbf{x}) S_j(\mathbf{x}) \, d\Gamma \mathbf{F}_{n,j} \\ M_{ij} \frac{\partial \mathbf{q}_j^s}{\partial t} = \sum_l K_{ij}^l \mathbf{F}_{l,j}^s + \mathbf{F}_i^* \end{gathered} $$

où $\mathbf{F}_{l,j}^s$ est le $i$'ième flux spécifique évalué au nœud $j$, et $S_{ip}$ est la fonction de forme attachée au $i$'ième nœud évaluée au $p$'ième point matériel, $S_{ip} = S_i(x_p)$. La transformation des intégrales en sommes peut conduire à une sous-intégration de la matrice masse si un seul point matériel est contenu dans une cellule [7]. La matrice obtenue est alors singulière et on compense ce défaut en condensant cette dernière sur sa diagonale. L'intégration temporelle est effectuée grâce à un algorithme explicite de type Euler ou Runge-Kutta d'ordre deux sur l'intervalle de temps : $[0, T] = \cup_n [t^{n-1}, t^n]$, $n \in \mathbb{N}^*$.

## 2.4 Transfert entre les nœuds et les points matériels

La grille de calcul a pour objet de permettre le calcul des gradients. Les champs sont projetés aux nœuds pour la résolution des équations de la dynamique par interpolation sur chaque cellule des quantités conservées $\mathbf{q}$ :

$$ \mathbf{q}_i^n = \sum_{p=1}^{N_p} S_{ip} \mathbf{q}_p^n \quad \rightarrow \quad m_i \mathbf{q}_i^{s,n} = \sum_{p=1}^{N_p} S_{ip} m_p \mathbf{q}_p^{s,n} $$

Cette interpolation contient une pondération par les masses nodales qui stabilise la MPM. La projection inverse utilise l'interpolation des quantités spécifiques :

$$ \mathbf{q}_p^{s,n+1} = \sum_{i=1}^{N_n} S_{ip} \mathbf{q}_i^{s,n+1} $$

où $N_n$ est le nombre de nœuds de l'élément contenant le point matériel $p$.

# 3 Calcul des flux aux interfaces

L'approche DG relaxe la condition de continuité des fonctions de forme entre deux éléments, les champs peuvent donc être discontinus sur les bords des cellules et des flux aux interfaces apparaîssent. Pour obtenir un suivi des fronts d'onde efficace, il est intéressant d'utiliser la structure caractéristique des équations hyperboliques et de l'introduire dans le schéma numérique. Une manière de faire est de définir des problèmes de Riemann consistant à trouver $\mathbf{q}$ satisfaisant le problème aux valeurs initiales :

$$ \begin{gathered} \frac{\partial \mathbf{q}}{\partial t} + \frac{\partial \mathbf{F}}{\partial x} = \mathbf{0}, \quad -\infty < x < \infty, t > 0 \\ \mathbf{q}(x, 0) = \mathbf{q}_0(x) = \begin{cases} \mathbf{q}_L \text{ si } x < 0 \\ \mathbf{q}_R \text{ si } x > 0 \end{cases} \end{gathered} \quad (6) $$
où les indices *L* et *R* indiquent respectivement une solution présente à gauche (Left) ou à droite (Right) de la discontinuité. Les flux aux interfaces sont calculés par la méthode de Godunov [6, 8] correspondant à la solution stationnaire :

$$F^* = F(q^*(x/t = 0)) = F(q^*) \quad (7)$$

# 4 Cas tests unidimensionnels

On considère dans un premier temps un mouvement de solide rigide de translation pour observer l'influence du passage d'un point matériel d'une cellule à une autre sur la conservation des champs. Ensuite, l'impact de deux barres animées de vitesses initiales opposées est comparé avec la solution analytique [9] et aux résultats donnés par la MPM.

## 4.1 Mouvement de corps rigide

On considère une unique particule ayant une vitesse initiale $v_0 = 5$. Les positions successives de la particule au cours du temps ainsi que les champs de vitesse obtenus avec la DGMPM et la MPM sont tracés sur la figure (1).

<figure><img src="image_3.png" /><figcaption>FIGURE 1 – Résultats de la simulation du mouvement de solide rigide de la gauche vers la droite (MPM et DGMPM).</figcaption></figure>

On note en premier lieu que les deux méthodes testées offrent les mêmes résultats pour l'évolution de la particule : la vitesse initiale donnée à la particule est conservée tout au long de la translation rigide malgré plusieurs traversées de cellule. Ceci est un préalable nécessaire.

## 4.2 Impact de deux barres

### 4.2.1 Problème en petites perturbations

On considère à présent deux barres, constituées d'un même matériau élastique et possédant des vitesses initiales de directions opposées :

<figure><img src="image_4.png" /></figure>
où $v_0 = \frac{\sigma}{2\rho c}$, $\sigma < \sigma^y$, est une vitesse assurant de rester dans le domaine élastique, et $c = \sqrt{\frac{E}{\rho}}$ la célérité des ondes dans une barre. On montre sur la figure (2) le champ de contrainte analytique dans les deux barres à différents instants et les solutions numériques obtenues avec la DGMPM et la MPM.

<figure><img src="image_5.png" /><figcaption>FIGURE 2 – Solutions numériques et analytique du champ de contrainte pour la simulation d'un impact de deux barres : maillage de 150 points matériels et 150 éléments de longueur $\Delta x$, $E = 2.10^{11}$ Pa, $\rho = 7800$ kg·m$^{-3}$, $\sigma = 400.10^7$ Pa.</figcaption></figure>

La solution analytique est une discontinuité qu’un schéma numérique ne peut représenter correctement que si le nombre de Courant est égal à 1 ($\Delta t = \alpha \Delta x/c$, où $\alpha = 1$ est le nombre de Courant). On peut montrer par une analyse de von-Neumann que la DGMPM est stable à l’ordre 1 temporel pour $\alpha \le 1$ avec un seul point matériel par cellule (quelque soit leur position dans les éléments). Une conséquence de cette propriété est que la solution obtenue avec la DGMPM se superpose à la solution analytique.
La MPM en revanche ne peut pas représenter cette discontinuité car elle est instable pour $\alpha = 1$. Aussi, la formulation DGMPM réduit de manière significative les oscillations présentes dans la méthode de base, ce qui permet un suivi plus précis des fronts d'onde.

### 4.2.2 Introduction aux grandes transformations

On reprend le problème précédent avec cette fois une actualisation des positions des points matériels à chaque pas de temps. La mesure de déformation que l'on utilise n'est pas cohérente avec les grandes déformations, toutefois, cette simulation est un premier pas vers les grandes transformations. La figure (3) montre les positions de 10 points matériels situés dans la zone rouge et les champs de contrainte de Cauchy sur tout le domaine obtenus par la DGMPM et la MPM à des instants successifs ($\alpha = 0.5$ pour la MPM).

<figure><img src="image_6.png" /><figcaption>FIGURE 3 – Solutions numériques de la contrainte de Cauchy et positions successives d'un ensemble de points pour la simulation d'un impact de deux barres élastiques.</figcaption></figure>

Cette simulation montre l'apparition d'une instabilité pour les deux méthodes. Au début, la DGMPM décrit une discontinuité malgré le déplacement des points matériels. A la fin du calcul, certains points matériels partagent une cellule et on peut montrer que cette configuration entraîne une baisse du nombre de Courant pour assurer la stabilité. Cependant, l'analyse de von-Neumann permet d'obtenir une expression du nombre de Courant en fonction du nombre de points matériels présents dans une cellule et de leurs positions au sein de cette cellule pour les ordres d'intégration temporelle 1 et 2.
# 5 Impact 2D sur une plaque

Cette section a pour objet de sortir du cadre undimensionnel pour une première approche qualitative sur un problème élastique bidimensionnel en déformations planes. Soit une plaque rectangulaire de section $l \times h = 4 \times 3 \, m^2$ qui subit un impact sur une partie de son bord gauche ($a = 1 \, m$) :

<figure><img src="image_7.png" /><figcaption>FIGURE 4 – Plaque élastique impactée sur une partie de son bord gauche. Paramètres : $E = 2.10^{11}$ Pa, $\rho = 7800 \, kg \cdot m^{-3}$, $\sigma^d = 2.10^8$ Pa.</figcaption></figure>

La plaque est représentée par 33 points matériels dans la longueur et 25 dans la hauteur, avec un seul point par cellule. Un algorithme de Runge-Kutta d'ordre 2 est utilisé pour l'intégration temporelle.

<figure><img src="image_8.png" /><figcaption>FIGURE 5 – Champ de contrainte $\sigma_{xx}$ dans la plaque élastique (intégration temporelle RK2, $\alpha = 1$)</figcaption></figure>

La figure (5) montre la répartition de la contrainte $\sigma_{xx}$ dans la plaque avant et après réflexion de l'onde longitudinale sur le bord droit. Chaque point matériel est représenté par un carré sur lequel les champs sont constants. Ce problème implique la présence d'ondes longitudinales et transverses que la DGMPM peut capturer correctement avec une intégration temporelle d'ordre 2.

# 6 Intérêts de la méthode - Conclusion

On a présenté dans ce travail une nouvelle formulation de la Méthode des Points Matériels étendue au cadre Galerkin Discontinu qui diffère de l'originale sur plusieurs points. Premièrement, l'équation sur laquelle est basée la forme faible est la forme conservative (1), ce qui implique que le champ de contraintes est lui aussi défini aux nœuds. Cela constitue une méthode qui est plus proche de la *Fluid Implicit Particle* [10] que de la MPM. La base d'approximation continue par élément fait apparaître des flux aux interfaces entre les cellules qui sont ici calculés par la méthode de Godunov après résolution d'un problème de Riemann élastique. La DGPIC [5] n'utilise pas les grandeurs spécifiques dans la forme faible et conduit à une forme semi-discrète dans laquelle la matrice masse ne prend pas en compte les positions des points matériels, aussi, une diffusion numérique plus forte et un transfert des champs plus compliqué rendent moins efficace la méthode.
Cette nouvelle approche permet de capter correctement des fronts d'onde, ce qui n'est pas réalisable avec la MPM. Ceci est possible grâce à l'utilisation du solveur de Riemann qui permet d'introduire la solution analytique dans la solution numérique. La réduction du domaine d'influence des points matériels grâce aux fonctions de forme continues sur un élément joue également un rôle important dans le suivi plus précis des fronts d'onde. La capacité des particules à se déplacer d'une cellule à une autre est prometteuse pour une formulation en grandes transformations. En revanche, à l'instar de le FEM classique, la DGFEM est sensible aux distorsions de maillage. Toutefois, on peut tirer profit de cette dernière puisque la propriété des fonctions de forme continues par élément nous autorise à monter localement en ordre d'approximation, notamment dans le cas où une discontinuité est présente dans une cellule.

L'extension à des solides élastoviscoplastiques peut être réalisée à l'aide d'algorithmes fractionnés, comme montré dans [11]. La formulation en grandes transformations est en cours et les simulations en 2D sont en voie d'extension pour gérer la présence de plusieurs points matériels dans une cellule. L'analyse de stabilité est en finalisation pour les problèmes 1D et sera menée pour les problèmes 2D.

Références

[1] Deborah Sulsky, Shi-Jian Zhou, and Howard L Schreyer. Application of a particle-in-cell method to solid mechanics. *Computer Physics Communications*, 87(1-2):236-252, 1995.

[2] Francis H Harlow. The particle-in-cell method for numerical solution of problems in fluid dynamics. Technical report, Los Alamos Scientific Lab., N. Mex., 1962.

[3] Wm H Reed and TR Hill. Triangularmesh methodsfor the neutrontransportequation. *Los Alamos Report LA-UR-73-479*, 1973.

[4] Jan S Hesthaven and Tim Warburton. *Nodal discontinuous Galerkin methods : algorithms, analysis, and applications*. Springer Science & Business Media, 2007.

[5] GB Jacobs and Jan S Hesthaven. High-order nodal discontinuous galerkin particle-in-cell method on unstructured grids. *Journal of Computational Physics*, 214(1):96-121, 2006.

[6] Randall J Leveque. *Finite volume methods for hyperbolic problems*. Cambridge university press, 2002.

[7] E. Love and D.L. Sulsky. An unconditionally stable, energy-momentum consistent implementation of the material-point method. *Computer Methods in Applied Mechanics and Engineering*, 195(33-36):3903-3925, 2006.

[8] Eleuterio F Toro. *Riemann solvers and numerical methods for fluid dynamics : a practical introduction*. Springer Science & Business Media, 2013.

[9] Lili Wang. *Foundations of stress waves*. Elsevier, 2011.

[10] JU Brackbill and HM Ruppel. Flip : A method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions. *Journal of Computational Physics*, 65(2):314-343, 1986.

[11] Thomas Heuzé. Simulation d'impacts sur des solides élastique-viscoplastiques avec la méthode des volumes finis. In *Colloque National en Calcul des Structures*, Giens 2017.