# Simulation d'impacts sur des solides élastique-viscoplastiques avec la méthode des volumes finis

Thomas Heuzé

► To cite this version:

Thomas Heuzé. Simulation d'impacts sur des solides élastique-viscoplastiques avec la méthode des volumes finis. 13e colloque national en calcul des structures, Université Paris-Saclay, May 2017, Giens, Var, France. hal-01899360

HAL Id: hal-01899360

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

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Simulation d'impacts sur des solides élastique-viscoplastiques avec la méthode des volumes finis

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

**Résumé** — Une méthode des volumes finis est utilisée et étendue dans ce travail pour la simulation d'impacts sur des solides élastique-viscoplastiques pour des maillages de quadrangles, bidimensionnels structurés et non cartésiens. La formulation est basée sur la décomposition de différence des flux aux interfaces entre cellules [1], est du second ordre au moyen de limiteurs de flux, intègre la contribution de flux transverses générés par les cellules de coin, et utilise un algorithme fractionné du second ordre pour le calcul du second membre. La méthode est ici illustrée sur un problème d'impact sur un volume élémentaire hétérogène contenant une inclusion. Une comparaison est effectuée avec une solution éléments finis obtenue avec le code Cast3M.

**Mots clés** — Volumes finis, élastoviscoplasticité, impact.

## 1 Introduction

La simulation numérique de problèmes aux valeurs limites et initiales hyperboliques incluant des conditions de chargement extrêmes comme des impacts nécessite la capacité de capturer précisément les fronts d'onde générés dans le milieu considéré. Cela permet de suivre les trajets d'onde et de comprendre de façon plus fine les phénomènes mécaniques ayant lieu. Pour les milieux solides, cela permet aussi un calcul précis des déformations inélastiques, et donc celui des contraintes et distorsions résiduelles. Les procédés de formage à grande vitesse comme le formage magnétique sont des exemples d'applications dans lesquels des conditions sévères de chargement conduisent à la génération de fronts d'onde qu'il est important de suivre pour comprendre la génération des déformations plastiques dans la pièce à former, et ainsi pouvoir à terme optimiser correctement sa forme finale. Ces objectifs nécessitent la construction de schémas numériques capable d'allier un haut ordre de précision et une haute résolution des discontinuités sans oscillations numériques parasites.

La simulation numérique d'impacts sur des solides dissipatifs est encore principalement réalisée 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. Bien que cette approche possède des atouts bien connus, 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. La méthode des volumes finis, initialement développée pour la dynamique des gaz [1], est de plus en plus utilisée pour la simulation d'impacts sur des milieux solides [2, 3, 4, 5]. Cette méthode présente un certain nombre d'intérêts pour assurer un suivi précis des fronts d'onde ; entre autres (i) la continuité des champs n'est pas imposée sur le maillage dans sa version cellule-centrée, ce qui permet la représentation de solutions discontinues, (ii) la structure caractéristique des équations hyperboliques est introduite dans la solution numérique, soit au travers de la résolution explicite de problèmes de Riemann aux interfaces entre cellules, soit de façon implicite de part la construction du schéma numérique, (iii) enfin la quantité de viscosité numérique introduite peut-être dosée localement en fonction de la régularité locale de la solution, de sorte à pouvoir éliminer les oscillations numériques parasites tout en préservant une précision de haut ordre dans les zones ou le champ est plus régulier.

On s'intéresse dans ce travail à des systèmes élastique-viscoplastiques conduisant à un système différentiel non-homogène, générant un système d'ondes faiblement discontinues après une onde discontinue associée à la transition élastique-élastoviscoplastique. On utilise et étend ici une méthode des volumes finis appliquée à ces milieux pour des maillages de quadrangles, bidimensionnels structurés et non cartésiens.

<footnote n="1">¹ GeM, Ecole Centrale Nantes, thomas.heuze@ec-nantes.fr</footnote>
siens. La formulation est basée sur la décomposition de différence des flux aux interfaces entre cellules [1], est du second ordre grâce à l'utilisation de limiteurs de flux, intègre la contribution de flux trans-verses générés par les cellules de coin, et utilise un algorithme fractionné du second ordre pour le calcul du second membre. L'approche est ici restreinte au cadre des petites perturbations, et est illustrée sur un problème d'impact sur un volume élémentaire hétérogène contenant une inclusion. Une comparaison est effectuée avec une solution éléments finis obtenue avec le code Cast3M.

## 2 Problème aux valeurs limites et initiales élastique-viscoplastique

On considère un solide élastique-viscoplastique avec une loi de fluage de type Odqvist [6] et un écrouissage cinématique linéaire, dont l'ensemble des équations de comportement s'écrit dans le cadre des petites perturbations :

$$ \dot{\epsilon} = \dot{\epsilon}^e + \dot{\epsilon}^p \qquad (1) $$

$$ \dot{\sigma} = C : \dot{\epsilon}^e \qquad (2) $$

$$ f = \sigma_{eq}(\xi) - \sigma_y; \quad \sigma_{eq} = \sqrt{\frac{3}{2}\xi : \xi}; \quad \xi = s - X; \quad \mathcal{C} = \{(\sigma, X, \dot{\epsilon}^p) | f \le 0\} \qquad (3) $$

$$ \dot{\epsilon}^p = \frac{3}{2} \dot{p} \frac{\xi}{\sigma_{eq}} \qquad (4) $$

$$ \dot{X} = \frac{2}{3} D \dot{\epsilon}^p \qquad (5) $$

$$ \dot{p} = \left\langle \frac{f}{K} \right\rangle^n \qquad (6) $$

où les équations (1), (2), (3), (4), (5) et (6) réfèrent respectivement à la partition additive du taux de dé-formation ($\dot{\epsilon}$) en parties élastique ($\dot{\epsilon}^e$) et viscoplastique ($\dot{\epsilon}^p$), la loi d'élasticité, le critère $f$ (ou contrainte visqueuse) associé à la norme de Mises et le convexe d'élasticité $\mathcal{C}$, la loi d'écoulement viscoplastique, la loi d'évolution de l'écrouissage cinématique et la loi de fluage. On note $\xi = s - X$, où $s$ est le déviateur des contraintes de Cauchy $\sigma$, $X$ la variable paramétrant la position du centre du convexe d'élasticité, $\langle x \rangle = (x + |x|)/2$ la partie positive de $x$, $C$ le tenseur des rigidités élastiques, $\dot{p}$ le taux de déformation viscoplastique cumulé, et $D$ le paramètre d'écrouissage.

À ces équations de comportement sont ajoutées les équations locales du mouvement et la définition du tenseur linéarisé des déformations ; on peut montrer qu'on obtient le système différentiel du premier ordre non-homogène suivant, écrit ici dans le cadre bidimensionnel des déformations planes :

$$ \frac{\partial \mathbf{u}}{\partial t} + \frac{\partial \mathbf{f}_x}{\partial x} + \frac{\partial \mathbf{f}_y}{\partial y} = \mathbf{S}(\mathbf{u}) \quad \forall (\mathbf{x} \times \mathbf{t}) \in (\Omega \times ]0, t_{\text{fin}}]) \qquad (7) $$

défini dans un domaine $\Omega$ et sur un intervalle de temps de durée $t_{fin}$, auquel il convient d'ajouter les conditions aux limites et initiales nécessaires. Les vecteurs inconnus $\mathbf{u}$, flux $\mathbf{f}_x$ et $\mathbf{f}_y$, et le second membre $\mathbf{S}(\mathbf{u})$ sont définis comme suit :

$$ \mathbf{u} = \begin{bmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \\ u \\ v \end{bmatrix} ; \quad \mathbf{f}_x = \begin{bmatrix} -(\lambda + 2\mu)u \\ -\lambda u \\ -\mu v \\ -\sigma_{xx}/\rho \\ -\sigma_{xy}/\rho \end{bmatrix} ; \quad \mathbf{f}_y = \begin{bmatrix} -\lambda v \\ -(\lambda + 2\mu)v \\ -\mu u \\ -\sigma_{xy}/\rho \\ -\sigma_{yy}/\rho \end{bmatrix} ; \quad \mathbf{S}(\mathbf{u}) = -3 \frac{\mu}{\sigma_{eq}} \left\langle \frac{f}{K} \right\rangle^n \begin{bmatrix} \xi_{xx} \\ \xi_{yy} \\ \xi_{xy} \\ 0 \\ 0 \end{bmatrix} \qquad (8) $$

où $\lambda$ et $\mu$ sont les paramètres de Lamé, $u$ et $v$ les composantes du champ de vitesse, et $\rho$ la masse volumique. L'analyse spectrale de la jacobienne de la partie homogène du système :

$$ g = \frac{\partial(F \cdot \mathbf{n})}{\partial \mathbf{u}} = \frac{\partial \mathbf{f}_x}{\partial \mathbf{u}} n_x + \frac{\partial \mathbf{f}_y}{\partial \mathbf{u}} n_y = \mathbf{A} n_x + \mathbf{B} n_y = - \begin{bmatrix} 0 & 0 & 0 & n_x(\lambda + 2\mu) & \lambda n_y \\ 0 & 0 & 0 & n_x\lambda & n_y(\lambda + 2\mu) \\ 0 & 0 & 0 & n_y\mu & n_x\mu \\ \frac{n_x}{\rho} & 0 & \frac{n_y}{\rho} & 0 & 0 \\ 0 & \frac{n_y}{\rho} & \frac{n_x}{\rho} & 0 & 0 \end{bmatrix} \qquad (9) $$
donne les vitesses caractéristiques $\lambda_1 = -c_L$, $\lambda_2 = c_L$, $\lambda_3 = -c_s$, $\lambda_4 = c_s$, $\lambda_5 = 0$, où $c_L = \sqrt{(\lambda + 2\mu)/\rho}$ et $c_s = \sqrt{\mu/\rho}$ désignent les célérités associées aux ondes de traction/compression et de cisaillement respectivement; et les directions caractéristiques (vecteurs propres) associés $K^{(p)}$, $1 \le p \le 5$.

# 3 La méthode des volumes finis

La méthode des volumes finis est basée sur la sous-division du domaine spatial en cellules au sein desquelles, pour sa version cellule-centrée, une approximation $U_i$ du champ $\mathbf{u}$ est définie dans la cellule $i$ par moyennation. Cette approximation est actualisée au moyen des lois de bilan (7) écrites sous forme intégrale dans la cellule, couplées à un schéma d'intégration explicite en temps.

## 3.1 Algorithme fractionné pour le second membre

La prise en compte du second membre du système différentiel (7) peut être opérée au cours de la résolution au moyen d'un algorithme fractionné, ou étagé [1]. Pour un algorithme d'ordre un, la première étape consiste à résoudre la partie homogène de (7) sur un pas de temps, associée à la condition initiale du champ connu à l'instant $t_n$ ($n$ référant au numéro du pas de temps). La solution obtenue $\bar{\mathbf{u}}^{n+1}$ est alors utilisée comme condition initiale du problème aux valeurs initiales suivant :

$$ 
\begin{aligned}
\frac{d\mathbf{u}}{dt} &= \mathbf{S}(\mathbf{u}) \\
\mathbf{u}(0) &= \bar{\mathbf{u}}^{n+1}
\end{aligned}
\qquad (10) $$

On montre par ailleurs que cette décomposition est exacte pour une équation scalaire linéaire d'advection-réaction. Pour le système non linéaire considéré, elle n'est en revanche qu'une approximation. En nommant $C$ et $S$ les opérateurs respectivement associés aux parties homogène et non-homogène, l'actualisation prend la forme

$$ \mathbf{u}^{n+1} = S^{(\Delta t)} C^{(\Delta t)} \mathbf{u}^n \qquad (11) $$

Mieux, un algorithme du second ordre

$$ \mathbf{u}^{n+1} = S^{(\Delta t/2)} C^{(\Delta t)} S^{(\Delta t/2)} \mathbf{u}^n \qquad (12) $$

permet un calcul plus précis des déformations viscoplastiques et sera ici utilisé.

## 3.2 Formulation en décomposition de différence de flux

On présente ici une formulation des volumes finis basée sur la décomposition de différence des flux aux interfaces entre cellules [1], appliquée à des maillages structurés de quadrangles, ayant a priori une répartition non cartésienne. Cette version intègre le cadre des méthodes de haut ordre à variation totale non-croissante, ayant l'intérêt d'allier à la fois un haut ordre de précision dans les zones où la solution est régulière, et une haute résolution des discontinuités du champ sans oscillations parasites. Leur intérêt réside dans leur capacité à introduire une quantité de viscosité numérique controllée de façon locale, de sorte à pouvoir s'adapter localement à la régularité de la solution. Cette formulation s'applique à la partie homogène de (7), de sorte qu'une variation totale non-croissante de la solution ait un sens.

Soit la cellule élémentaire représentée sur la figure 1, d'aire $|A_i|$, dont l'arête $s$ ($1 \le s \le 4$) délimitée par les points $P_s$ et $P_{s+1}$, de longueur $L_s$, de normale extérieure $\mathbf{n}_s$, et dont les états à gauche ($L$) et à droite ($R$) sont connus à l'instant $t_n$. L'actualisation de l'état de cette cellule à l'instant $t_{n+1}$ s'écrit :

$$ \mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{|A_i|} \left( \sum_{k=1}^{P} L_k \mathcal{A}_k^+ \Delta \mathbf{U}_k + \sum_{l=1}^{Q} L_l \mathcal{A}_l^- \Delta \mathbf{U}_l \right) - \frac{\Delta t}{|A_i|} \left( \sum_{k=1}^{P} L_k \tilde{\mathbf{F}}_k^{\text{out}} - \sum_{l=1}^{Q} L_l \tilde{\mathbf{F}}_l^{\text{in}} \right) \qquad (13) $$

Cette équation comprend deux groupes de termes qui sont détaillés ci-dessous.
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Volume fini bidimensionnel quadrangulaire</figcaption></figure>

### 3.2.1 Flux de Godunov

Le premier groupe de terme dans l'équation (13) fait apparaître les fluctuations, représentées par les opérateurs $\mathcal{A}_k^{\pm} \Delta \mathbf{U}_k$ :

$$ \mathcal{A}_k^{\pm} \Delta \mathbf{U}_k = \sum_{p=1}^{M_w} \lambda_p^{\pm} \mathcal{W}_k^{(p)} = \sum_{p=1}^{M_w} \lambda_p^{\pm} \alpha_k^{(p)} \mathbf{K}_k^{(p)} \quad (14) $$

quantifiant l'effet de toutes les ondes se propageant respectivement vers la droite ou vers la gauche dans le repère local de l'arête $k$. Ces fluctuations font apparaître les parties positives (resp. négatives) des vitesses caractéristiques $\lambda_p^{\pm}$, ainsi que les directions associées $\mathbf{K}_k^{(p)} \equiv \mathbf{K}^{(p)}(\mathbf{n}_k)$. Chaque onde est pondérée par un coefficient $\alpha_k^{(p)}$ déterminé par la projection du saut au travers de l'arête $k$ du vecteur inconnu $\Delta \mathbf{U}_k = (\mathbf{U}_R - \mathbf{U}_L)_k$ sur sa base caractéristique $\mathbf{K}_k^{(p)}$ ($1 \le p \le M_w$, $M_w = 5$ ici) :

$$ \Delta \mathbf{U}_k = \sum_{p=1}^{M_w} \mathcal{W}_k^{(p)} = \sum_{p=1}^{M_w} \alpha_k^{(p)} \mathbf{K}_k^{(p)} = \mathbf{K} \alpha_k \quad (15) $$

On définit ainsi la "force" d'une onde par $\mathcal{W}_k^{(p)} = \alpha_k^{(p)} \mathbf{K}_k^{(p)}$. Cette projection revient à résoudre un problème de Riemann (linéaire) associé à l'arête $k$. Les fluctuations associées à chaque arête sont sommées pour calculer la contribution des termes du premier ordre à l'actualisation de la cellule $i$; cette sommation est effectuée sur des fluctuations négatives pour les arêtes ayant une normale sortante, et respectivement positives pour les arêtes ayant une normale entrante, de sorte que $P + Q = N$, $N$ nombre d'arêtes de la cellule. On peut montrer [1] que ces termes (les fluctuations) sont équivalents à appliquer la méthode de Godunov dont les flux aux interfaces des cellules sont évalués en fonction de la partie stationnaire $(x/t = 0)$ de la solution du problème de Riemann.

### 3.2.2 Flux additionnels

Le second groupe de termes dans l'équation (13) représentent des flux additionnels définis aux arêtes, soit entrant (in) soit sortant (out) en fonction de la normale de l'arête, qui sont composés de deux types de contributions

$$ \tilde{\mathbf{F}}_l^{\text{in}} = \tilde{\mathbf{F}}_l^{\text{HO}} + \tilde{\mathbf{F}}_l^{\text{tran}}. \quad (16) $$

#### Flux de haut ordre

La première de ces contributions permet d'atteindre un plus haut ordre (HO) de précision (ordre deux ici), définie en fonction d'ondes dont la "force" $\mathcal{W}_l^{(p)}$ est limitée, alors notée $\tilde{\mathcal{W}}_l^{(p)} = \tilde{\alpha}_l^{(p)} \mathbf{K}_l^{(p)}$ :

$$ \tilde{\mathbf{F}}_l^{\text{HO}} = \frac{1}{2} \sum_{p=1}^{M_w} |\lambda_l^{(p)}| \left( 1 - \frac{\Delta t}{\Delta s_l} |\lambda_l^{(p)}| \right) \tilde{\mathcal{W}}_l^{(p)} \quad (17) $$

où $\Delta s_l$ désigne la distance entre les barycentres des deux cellules partageant l'arête $l$. Les ondes sont limitées en fonction d'un ratio de décentrement $\theta_l^{(p)}$ défini pour l'onde $p$ sur l'arête $l$ comme :

$$ \theta_l^{(p)} = \frac{\mathcal{W}_l^{(p)} \cdot \mathcal{W}_l^{(p)}}{\| \mathcal{W}_l^{(p)} \|^2} \quad (18) $$
où J désigne l'arête amont, c'est-à-dire l'arête opposée à l appartenant à la cellule L située à gauche dans le repère local de l'arête (voir figure 1) si $\lambda_l^{(p)} > 0$, ou bien l'arête opposée à l appartenant à la cellule R située à droite dans le repère local de l'arête si $\lambda_l^{(p)} < 0$. Cependant, pour des quadrangles non cartésiens les arêtes amont et aval de l'arête l ne possèdent pas nécessairement la même normale. Aussi, le calcul du ratio (18) doit être effectué avec des "forces" d'onde calculées dans un même repère local, celui de l'arête l. Les "forces" des ondes $W_j^{(p)}$ sont donc calculées à partir de coefficients pondérants $\alpha_j^{(p)}$ recalculés dans la base de l'arête l :

$$ \alpha_j^{(p)} = \mathbf{K}^{-1}(\mathbf{n}_l) \cdot \Delta U_j \quad (19) $$

où $\Delta U_j$ désigne le saut au travers de l'arête J du vecteur inconnu. Le ratio de décentrement (18) peut être compris comme une certaine mesure de la régularité locale de la solution, la "force" de l'onde p associée à l'arête l est alors limitée au travers d'une fonction limitante $\phi(\theta_l^{(p)})$ :

$$ \tilde{\alpha}_l^{(p)} = \phi(\theta_l^{(p)}) \alpha_l^{(p)} \quad (20) $$

De nombreuses fonctions limitantes existent et permettent de reconstruire différents schémas volumes finis existants [7]. Certaines permettent au schéma numérique de satisfaire à une variation totale non-croissante, c'est-à-dire que cela permet d'éviter l'apparition d'oscillations numériques parasites au voisinage de discontinuités du champ. On utilise dans ce travail le limiteur Superbee, défini par $\phi(\theta) = \max(0, \min(1, 2\theta), \min(2, \theta))$, qui intègre ce cadre. D'une façon générale la limitation locale de la "force" des ondes revient à ajouter localement un peu de viscosité numérique au schéma, et ainsi à abaisser l'ordre du schéma pour pouvoir capturer correctement les discontinuités. Dans les zones plus régulières, la limitation n'est pas effective et une précision d'ordre 2 peut être atteinte.

## Flux transverses

La seconde contribution aux flux additionnels (16) permet d'améliorer la stabilité du schéma numérique, de sorte qu'un nombre de Courant théorique égal à un peut être utilisé. Ces flux permettent de prendre en compte l'information se propageant en biais par rapport à la cellule considérée ; c'est la contribution à la cellule dont on souhaite actualiser l'état d'une cellule ne partageant qu'un noeud (mais pas d'arête) avec celle-ci. C'est la Corner Transport Upwind (CTU) method [1]. Ces flux sont d'une grande importance pour assurer la stabilité du calcul pour des système élastiques ; l'élasticité couple en effet les composantes de déformation au travers de l'effet Poisson, de sorte qu'une information transverse à la cellule considérée doit être intégrée au schéma numérique.

Considérons le patch de cellules de la figure 2. On considère l'arête dénommée (i) dont le repère local $(\mathbf{n}_i, \mathbf{t}_i)$ est représenté. Cette arête donne lieu au calcul de fluctuations normales $A_i^+ \Delta U$ et $A_i^- \Delta U$

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Fluctuations normales et tangentielles définies à partir de l'arête (i).</figcaption></figure>

contribuant respectivement aux cellules R et L. Ces fluctuations donnent lieu à l'existence de fluctuations
transverses contribuant aux cellules voisines au travers des arêtes (j) et (k) pour la cellule L, et au travers des arêtes (m) et (l) pour la cellule R. Ces fluctuations transverses sont calculées en décomposant les fluctuations normales sur la structure caractéristique liée au problème de Riemann défini sur l'arête adjacente. On parle de *solveur de Riemann transverse*. Par exemple la fluctuation normale négative se décompose sur la structure associée à l'arête (j) comme

$$ \mathcal{A}_i^- \Delta \mathbf{U} = \sum_{p=1}^{M_w} \beta_p \mathbf{K}_j^{(p)} = \mathbf{K}_j \beta \quad (21) $$

où $\mathbf{K}_j$ tient compte ici de la normale $\mathbf{n}_j$ de l'arête (j), mais aussi des propriétés matériaux différentes entre les cellules L et T. La détermination des coefficients $\beta_p$ s'obtient ici de façon analytique pour les équations de l'élasticité plane linéaire. Par convention, on considère dans ce travail la normale extérieure, aussi les fluctuations transverses sont calculées avec l'opérateur positif $\mathcal{B}^+$, c'est-à-dire que seules les ondes ayant une célérité positive contribueront à cette fluctuation transverse :

$$ \mathcal{B}_j^+ \mathcal{A}_i^- \Delta \mathbf{U} = c_L \mathcal{W}^{(2)} + c_s \mathcal{W}^{(4)} \quad (22) $$

Un flux numérique additionnel est donc construit à partir de ces fluctuations transverses :

$$ \tilde{\mathbf{F}}_j^{\text{tran}} = \frac{\Delta t}{2\Delta s_j} \mathcal{B}_j^+ \mathcal{A}_i^- \Delta \mathbf{U}_i \quad (23) $$

qui contribuera à la cellule T. Le flux $\tilde{\mathbf{F}}_j^{\text{tran}}$ associé à l'arête l apparaissant dans l'équation (16) représente donc les contributions transverses des cellules adjacentes. Plus généralement, ce schéma de volumes finis est lié à un développement de Taylor de la solution au voisinage d'une cellule. Dans ce cadre, les flux de Godunov correspondent aux termes du premier ordre, les flux dits de haut ordre à ceux du second ordre, et les flux transverses aux dérivées croisées.

## 4 Volume élémentaire hétérogène chargé par un créneau de vitesse

On considère un volume bidimensionnel élémentaire, hétérogène, carré de coté $2a$, avec une inclusion de section circulaire de rayon $R$ centrée dans le volume, dont le maillage est représenté sur la figure 3. Ce volume est chargé sur son côté gauche par un créneau positif de la composante horizontale de vitesse sur un temps $\tau$, la vitesse normale est imposée nulle sur ses faces inférieure et supérieure, et une frontière absorbante est prescrite sur la face de droite. On suppose un matériau hétérogène arbitraire, constitué d'une inclusion en aluminium et d'une matrice en acier. L'analyse est menée en déformations planes. Les données numériques utilisées pour cette étude numérique sont regroupées dans la table 1.

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Maillage</figcaption></figure>

<table><caption>TABLE 1 – Données numériques de l'analyse.</caption><thead><tr><th>Géométrie</th><th>Matrice</th><th>Inclusion</th></tr></thead><tbody><tr><td>a = 10<sup>-3</sup>m<br>R = 5 × 10<sup>-3</sup>m</td><td>E<sub>M</sub> = 200GPa<br>v<sub>M</sub> = 0.3</td><td>E<sub>I</sub> = 700GPa<br>v<sub>I</sub> = 0.34</td></tr><tr><td>Chargement</td><td>ρ<sub>M</sub> = 7800kg.m<sup>-3</sup></td><td>ρ<sub>I</sub> = 2700kg.m<sup>-3</sup></td></tr><tr><td>v<sub>d</sub> = 40m.s<sup>-1</sup><br>τ = 7 × 10<sup>-8</sup>s</td><td></td><td>σ<sub>yl</sub> = 350MPa<br>D = 2GPa</td></tr></tbody></table>

La matrice est supposée rester élastique, seule l'inclusion subit des déformations viscoplastiques. À titre de comparaison, une solution éléments finis est obtenue avec le code Cast3M, calculée avec une discrétisation temporelle implicite et des éléments bilinéaires. Cette solution est comparée à la solution volumes finis sur les figures 4, 5, et 6 aux instants respectifs $1 \times 10^{-7}$, $1.9 \times 10^{-7}$ et $3 \times 10^{-7}$ secondes. Le nombre de Courant est prescrit à 0.9 pour la solution volumes finis, et le même pas de temps est utilisé pour la solution éléments finis. Cette comparaison porte sur la composante $\sigma_{11}$ (en Pa) du champ de contrainte et sur la déformation viscoplastique cumulée $p$. Ces deux quantités sont tracées le long de la ligne médiane horizontale, et superposées pour les solutions volumes finis (DataSet0) et éléments finis (DataSet1). Les
<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Temps $t = 1 \times 10^{-7}$s</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5 – Temps $t = 1.9 \times 10^{-7}$s</figcaption></figure>

isovaleurs en haut à gauche représentent la solution volumes finis, celles en haut à droite la solution éléments finis. On observe sur la figure 4 le créneau de contrainte normale négative qui s'est formé par la condition de vitesse imposée sur la face gauche. Son front descendant s'est réfléchi à l'interface avec l'inclusion générant un état intermédiaire de contrainte supérieur (inférieur en valeur absolu) et un état de vitesse supérieur puisque l'impédance de la matrice est supérieure à celle de l'inclusion. Le champ de contrainte éléments finis montre des oscillations au voisinage des discontinuités, particulièrement au proche des faces supérieure et inférieure du volume où ces oscillations sont telles que des états de traction apparaissent. Les déformations viscoplastiques cumulées sont par contre quasi-confondues pour les deux solutions numériques. Sur la figure 5, on observe que l'onde réfléchie s'est transformée en onde de traction suite à l'interaction avec le front montant, puis s'est réfléchie sur la face gauche du volume. Les déformations viscoplastiques se propagent dans l'inclusion et progressent de façon retardée par rapport au créneau négatif de contrainte. Enfin sur la figure 6, ce même créneau se réfléchit sur la face arrière de l'inclusion générant des états de compression intermédiaire, alors que les déformations viscoplastiques finissent par croître sensiblement au proche de la face arrière de l'inclusion. D'une façon générale, la solution volumes finis permet d'obtenir les mêmes déformations viscoplastiques sans les oscillations sur les contraintes obtenues avec les éléments finis, la figure 7 illustre le champ de contrainte de cisaillement sur lequel on peut les observer.
<figure><img src="image_7.png" /><figcaption>FIGURE 6 – Temps $t = 3 \times 10^{-7}$s</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7 – Contrainte de cisaillement $\sigma_{12}$ à $t = 3 \times 10^{-7}$s</figcaption></figure>

## 5 Conclusion

On a présenté dans ce travail l'utilisation et l'extension d'une méthode de volumes finis pour des solides élastique-viscoplastiques pour des maillages de quadrangles structurés mais de géométrie non cartésienne. La technologie est basée sur la décomposition de différence des flux aux interfaces entre cellules et utilise un algorithme fractionné d'ordre deux pour le calcul du second membre. La comparaison avec une solution éléments finis obtenue avec le code Cast3M, sur un volume hétérogène contenant une inclusion chargé par un créneau de vitesse, montre que les déformations viscoplastiques calculées sont identiques. Cependant la solution volumes finis permet de s'affranchir des oscillations numériques parasites sur les champs de contrainte (et de vitesse) grâce à l'emploi de limiteurs de flux satisfaisant une variation totale non-croissante (de la partie homogène de (7)).

## Références

[1] R. Leveque, Finite volume methods for hyperbolic problems. Cambridge University PressMcGraw Hill, 2002.

[2] P.T. Barton, D. Drikakis, E.I. Romenskii. A high order eulerian godunov method for elastic-plastic flow in solids, International Journal for Numerical Methods in Engineering, 81, 453-484, 2010.

[3] A.L. Ortega, M. Lombardini, D. Pullin, D. Meiron, Numerical simulation of elastic-plastic solid mechanics using an eulerian stretch tensor approach and HLLD Riemann solver. Journal of Computational Physics, 257, 414-441, 2014.

[4] N. Favrie, S. Gavrilyuk. Mathematical and numerical model for nonlinear viscoplasticity, Phil. trans. R. Soc. A, 369, 2864-2880, 2011.

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

[6] J. Lemaitre, J.L. Chaboche, A. Benallal, R. Desmorat. Mécanique des matériaux solides, Dunod, 3ème édition, 2009.

[7] P.K. Sweby High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM Journal on Numerical Analysis, 21, 995-1011, 1984.