# Numerical simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method

Thomas Heuzé

► To cite this version:

Thomas Heuzé. Numerical simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method. 6th ECCOMAS Thematic Conference on Computational Methods in Structural Dynamics and Earthquake Engineering (COMPDYN 2017), Jun 2017, Rhodes island, Greece. hal-02279799

HAL Id: hal-02279799

https://hal.science/hal-02279799

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 Attribution 4.0 International License
NUMERICAL SIMULATION OF IMPACTS ON
ELASTIC-VISCOPLASTIC SOLIDS WITH THE FLUX-DIFFERENCE
SPLITTING FINITE VOLUME METHOD

T. Heuzé¹

¹Research Institute in Civil and Mechanical Engineering (GeM, UMR 6183 CNRS)
École Centrale de Nantes, 1 rue de la Noë, F-44321 Nantes, France
e-mail: thomas.heuze@ec-nantes.fr

**Keywords:** Flux difference splitting finite volume method, Elastic-viscoplastic solids, Dynamic impact, Curvilinear structured meshes.

**Abstract.** A finite volume method is used and extended for the numerical simulation of impacts on elastic-viscoplastic solids with bidimensional curvilinear structured meshes. The formulation is based on the flux-difference splitting method [1], has second order accuracy through flux limiters, embeds the corner transport upwind method, and uses the second order accurate Strang splitting method to handle the right hand side of the system of balance laws. The approach is here derived with a Chaboche-type [2] elastic-viscoplastic solid within the small strain framework, and is illustrated on a problem of impact on a heterogeneous volume containing an inclusion. A comparison is performed with a finite element solution obtained with the finite element code Cast3M [3].
# 1 INTRODUCTION

The numerical simulation of hyperbolic initial boundary value problems including extreme loading conditions such as impacts requires the ability to accurately capture and track the wave front of shock waves induced in the medium. Indeed, this permits to correctly follow the path of waves and hence understand the mechanical phenomena occuring within that medium. For solid-type media, it allows also for an accurate assessment of the plastic strain field and hence that of residual stresses and distortions within the structure. High speed forming processes like electromagnetic material forming are some application examples of severe loading conditions in which the track of wave fronts is important both for understanding the development of plastic strain in the workpiece and optimizing its final shape. Hence, these problems require numerical schemes able to meet high orders of accuracy and a high resolution of discontinuity without any spurious oscillations.

The numerical simulation of impacts on dissipative solids has been and is again mainly performed with the classical finite element method coupled with centered differences or Newmark finite difference schemes in time [4]. Though this approach has well-known advantages, classical time integrators introduce high frequency noise in the vicinity of discontinuities which is hard to remove with artificial viscosity without destroying the accuracy of the numerical solution. The finite volume method, initially developed for the simulation of gas dynamics [1], has gained recently more and more interest for problems involving impacts on solid media [5, 6, 7, 8]. This method shows some advantages to achieve an accurate tracking of wavefronts; among others (i) continuity of fields is not enforced on the mesh in its cell-centered version, that allows for capturing discontinuous solutions, (ii) the characteristic structure of hyperbolic equations can be introduced within the numerical solution, either through the explicit solution of a Riemann problem at cell interfaces, or in an implicit way through the construction of the numerical scheme, (iii) the amount of numerical viscosity introduced can be controlled locally as a function of the local regularity of the solution, so that to permit the elimination of spurious numerical oscillations while preserving a high order of accuracy in more regular zones.

We are here interested in elastic-viscoplastic systems, leading to a non-homogeneous system of partial differential equations, generating a system of weakly discontinuous waves beyond the viscoplastic yield, following a discontinuous (elastic) wave due to the transition between elastic and elastic-viscoplastic ranges. In this work, we use and extend a finite volume method applied to these dissipative media with bidimensional curvilinear structured meshes of quadrangles. The formulation is based on the flux-difference splitting method [1], has second order spatial accuracy through flux limiters, embeds the corner transport upwind method, and uses the second order accurate Strang splitting method to handle the right hand side of the system of balance laws. The approach is here derived with a Chaboche-type [2] elastic-viscoplastic solid within the small strain framework, and is illustrated on a problem of impact on a heterogeneous volume containing an inclusion. A comparison is performed with a finite element solution obtained with the finite element code Cast3M [3].

# 2 ELASTIC-VISCOPLASTIC INITIAL BOUNDARY VALUE PROBLEM

Let's consider an elastic-viscoplastic solid with an Odqvist creep law [2] (or Chaboche-type model) and a linear kinematic hardening, the associated constitutive equations written within
the small strain framework are summarized below:

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

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

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

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

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

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

where equations (1), (2), (3), (4), (5) and (6) respectively refer to the additive partition of the strain rate ($\dot{\boldsymbol{\epsilon}}$) into elastic ($\dot{\boldsymbol{\epsilon}}^e$) and viscoplastic ($\dot{\boldsymbol{\epsilon}}^p$) parts, the elastic law, the yield criterion $f$ (or viscous stress) associated to Mises' norm and the elastic convex $\mathcal{C}$, the viscoplastic flow rule, the evolution law of the kinematic hardening and the creep law. One defines $\boldsymbol{\xi} = \boldsymbol{s} - \boldsymbol{X}$, where $\boldsymbol{s}$ is the deviatoric part of Cauchy stresses $\boldsymbol{\sigma}$, $\boldsymbol{X}$ the variable defining the center of the elastic convex, $\langle x \rangle = (x + |x|)/2$ the positive part of $x$, $\boldsymbol{C}$ the elastic stiffness tensor, $\dot{p}$ the rate of cumulated viscoplastic strain, and $\boldsymbol{D}$ the hardening parameter.

To these constitutive equations, are added the equations of motion and the relations between the linearized strain rate tensor $\dot{\boldsymbol{\epsilon}}$ and the velocity field. One can show that we obtain the following non-homogeneous system of partial differential equations, written here for the bidimensional plane strain case:

$$ \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 t) \in (\Omega \times ]0, t_{\text{end}}]) \qquad (7) $$

defined within a domain $\Omega$ and a time interval of duration $t_{\text{end}}$, and has to be supplemented with appropriate initial and boundary conditions. The vector of unknowns $\mathbf{u}$, fluxes $\mathbf{f}_x$ and $\mathbf{f}_y$, and the right hand side $\mathbf{S}(\mathbf{u})$ are defined as follows:

$$ \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_{\text{eq}}} \left\langle \frac{f}{K} \right\rangle^n \begin{bmatrix} \xi_{xx} \\ \xi_{yy} \\ \xi_{xy} \\ 0 \\ 0 \end{bmatrix} \qquad (8) $$

where $\lambda$ and $\mu$ denote the Lamé's parameters, $u$ and $v$ the velocity components, and $\rho$ the mass density. The eigen-analysis of the jacobian matrix associated to the homogeneous part of (7)

$$ \mathcal{J} = \frac{\partial(\mathcal{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) $$

gives the characteristic speeds $\lambda_1 = -c_L$, $\lambda_2 = c_L$, $\lambda_3 = -c_s$, $\lambda_4 = c_s$, $\lambda_5 = 0$, where $c_L = \sqrt{(\lambda + 2\mu)/\rho}$ and $c_s = \sqrt{\mu/\rho}$ are the celerities associated to tensile/compression and shear waves respectively; and associated characteristic directions (eigenvectors) $\mathbf{K}^{(p)}$, $1 \le p \le 5$.
# 3 FINITE VOLUME METHOD

The finite volume method is based on subdividing the computational domain in elementary cells within which, for the cell-centered version, an approximation $U_i$ of the field $u$ is defined in the cell $i$ by integral averaging. This approximation is updated by means of balance laws (7) written in integral form on each grid cell, coupled to an explicit time integration scheme.

## 3.1 Strang splitting for the right hand side

The solution procedure to account for the right hand side of the differential system (7) can be carried out by means of a splitting algorithm [1]. For first order accuracy (Godunov splitting), a first step consists in solving the homogeneous part of (7) on the time step associated to the initial condition known at time $t_n$ ($n$ referring to the time step number). The obtained solution $\bar{u}^{n+1}$ is then used as an initial condition for the following initial value problem:

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

This decomposition is exact for a linear advection-reaction equation. For the nonlinear system considered here, it is an approximation. Naming $C$ and $S$ operators associated to the solution of homogeneous and non-homogeneous parts, the update of the unknown field reads

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

In order to improve the accuracy of the computed viscoplastic strains, the second order accurate Strang splitting algorithm is used here:

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

## 3.2 Flux-difference splitting formulation

We use here a finite volume formulation based on flux-difference splitting [1], applied to curvilinear structured meshes of quadrangles. This version falls within the class of total variation non-increasing methods, that allow to meet both a high order of accuracy in zones where the solution field is regular and a high resolution of discontinuities without spurious numerical oscillations. Their strength relies on their ability to introduce a controlled amount of numerical viscosity locally, so that to adapt to the local regularity of the solution. This formulation is applied to the homogeneous part of the differential system (7) so that to be in accordance with a non-increasing total variation.

Let's consider the grid cell $i$ shown in figure 1, of area $|A_i|$, which edge $s$ ($1 \le s \le 4$) joining points $P_s$ and $P_{s+1}$ is of length $L_s$, has an outward unit normal $\mathbf{n}_s$, and has left ($L$) and right ($R$) states known at time $t_n$. The update of the cell state at time $t_{n+1}$ reads:

$$ \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) \quad (13) $$

The above formula consists of two groups of terms which are detailed below.
<figure><img src="image_2.png" /><figcaption>Figure 1: Quadrangular finite volume</figcaption></figure>

### 3.2.1 Godunov flux

The first group of terms in equation (13) shows fluctuations, denoted by operators $A_k^{\pm} \Delta U_k$:

$$A_k^{\pm} \Delta U_k = \sum_{p=1}^{M_w} \lambda_p^{\pm} W_k^{(p)} = \sum_{p=1}^{M_w} \lambda_p^{\pm} \alpha_k^{(p)} K_k^{(p)} \quad (14)$$

quantifying the effect of all waves travelling rightward or leftward respectively in the local frame of edge $k$. These fluctuations are defined with positive and negative parts of characteristic speeds $\lambda_p^{\pm}$, and associated characteristic directions $K_k^{(p)} \equiv K^{(p)}(n_k)$. Each wave is weighted with a coefficient $\alpha_k^{(p)}$ determined by the projection of the jump of the unknown vector $\Delta U_k = (U_R - U_L)_k$ across the edge $k$ onto the characteristic basis $K_k^{(p)}$ ($1 \le p \le M_w$, $M_w = 5$ here):

$$\Delta U_k = \sum_{p=1}^{M_w} W_k^{(p)} = \sum_{p=1}^{M_w} \alpha_k^{(p)} K_k^{(p)} = K \alpha_k \quad (15)$$

The wavestrength is defined by $W_k^{(p)} = \alpha_k^{(p)} K_k^{(p)}$. This projection amounts to solve a (linear) Riemann problem associated to edge $k$. The fluctuations associated to each edge are summed to compute the contribution of first order terms to the update of grid cell $i$; this summation is performed on negative fluctuations for edges having an outward normal, and on positive fluctuations for edges having an inward normal, so that $P + Q = N$, $N$ being the number of edges of grid cell $i$. The formulation with fluctuations is actually equivalent [1] to applying the Godunov's method which interface fluxes are computed with the steady part ($x/t = 0$) of the Riemann problem solution.

### 3.2.2 Additional fluxes

The second group of terms in equation (13) involve additional fluxes defined at edges, either inward (in) or outward (out) depending on the edge normal, that consist of two types of contributions

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

#### High order fluxes

The first contribution allows to reach a higher order (HO) of accuracy (order two here),
defined with wavestrength $W_l^{(p)}$ that have been limited, hence denoted $\tilde{W}_l^{(p)} = \tilde{\alpha}_l^{(p)} 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{W}_l^{(p)} \quad (17) $$

where $\Delta s_l$ refers to the distance between barycenters of grid cells sharing edge $l$. Waves are limited based on an upwind ratio $\theta_l^{(p)}$ defined for wave $p$ on edge $l$ as:

$$ \theta_l^{(p)} = \frac{W_J^{(p)} \cdot W_l^{(p)}}{\|W_l^{(p)}\|^2} \quad (18) $$

where $J$ denotes the upwind edge, that is the opposed edge to $l$ belonging to grid cell $L$ located to the left in the local frame of edge $l$ (see figure 1) if $\lambda_l^{(p)} > 0$, or the opposed edge belonging to grid cell $R$ located to the right if $\lambda_l^{(p)} < 0$. However, for noncartesian quadrangles upwind and downwind edges do not necessarily have the same normal. Thus, the computation of the upwind ratio (18) has to be performed with wavestrengths computed in a same local reference frame, that of edge $l$. Wavestrengths $W_J^{(p)}$ are hence computed with weighting coefficients $\alpha_J^{(p)}$ recomputed in the local frame of edge $l$:

$$ \alpha_J^{(p)} = \mathbf{K}^{-1}(\mathbf{n}_l) \cdot \Delta \mathbf{U}_J \quad (19) $$

where $\Delta \mathbf{U}_J$ denotes the jump across edge $J$ of the unknown vector. The upwind ratio (18) can be understood as a certain measure of the local regularity of the solution, the wavestrength of wave $p$ associated to edge $l$ is then limited by means of a limiting function $\phi(\theta_l^{(p)})$:

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

Many limiting functions exist and permit to obtain different known finite volume schemes [9]. Some of them enable the numerical scheme to satisfy a non-increasing total variation, so that the appearance of spurious numerical oscillations can be avoided in the vicinity of discontinuities. The limiter Superbee defined by $\phi(\theta) = \max(0, \min(1, 2\theta), \min(2, \theta))$ falls in this family, and is used in this work. More generally, the limitation of wavestrength amounts to add locally some numerical viscosity, and to locally lower the order of accuracy to properly capture discontinuities. In zones where the solution field is more regular, the limitation is not active and an accuracy of order two can be reached.

## Transverse fluxes

The second contribution to additional fluxes (16) enables to improve the stability of the numerical scheme, so that the Courant number can be set at one. These fluxes allows to account for informations travelling in bias with respect to the considered grid cell; this is the contribution to the grid cell to be updated of a cell only sharing a node (but not an edge) with it. This method is called the Corner Transport Upwind (CTU) method [1]. These fluxes are of great importance to ensure numerical stability for elastic media; indeed elasticity couples strain components through Poisson's effect, so that a transverse information to the considered grid cell should be introduced in the numerical scheme.

Let's consider the patch of grid cells shown in figure 2. One focuses on the edge denoted (i) which local frame $(\mathbf{n}_i, \mathbf{t}_i)$ is shown. This edge gives rise to the computation of normal
<figure><img src="image_3.png" /><figcaption>Figure 2: Normal and transverse fluctuations defined from edge i.</figcaption></figure>

fluctuations $A_i^+ \Delta U$ and $A_i^- \Delta U$ contributing to grid cells $R$ et $L$ respectively. These normal fluctuations leads to the computation of transverse fluctuations giving contribution to neighboring cells across edges (j) and (k) for cell $L$, and across edges (m) and (l) for cell $R$. These transverse fluctuations are computed by projecting normal fluctuations on the characteristic basis associated to the Riemann problem defined on the adjacent edge; it appears as a transverse Riemann solver. For instance the negative normal fluctuation is decomposed on the characteristic basis associated to edge (j) as

$$A_i^- \Delta U = \sum_{p=1}^{M_w} \beta_p K_j^{(p)} = K_j \beta \quad (21)$$

where $K_j$ accounts here for the normal $n_j$ of edge (j), but also of different material properties between grid cells $L$ and $T$. Coefficients $\beta_p$ are determined analytically for bidimensional plane strain elasticity equations. For outward normals shown in figure 2, the transverse fluctuations are computed with the positive operator $B^+$, that is only waves with positive characteristic speeds will contribute to this transverse fluctuation

$$B_j^+ A_i^- \Delta U = c_L W^{(2)} + c_s W^{(4)} \quad (22)$$

An additional numerical flux defined at edges is hence built from these transverse fluctuations:

$$\tilde{\mathbf{F}}_j^{\text{tran}} = \frac{\Delta t}{2\Delta s_j} B_j^+ A_i^- \Delta U_i \quad (23)$$

which contributes to grid cell $T$. The flux $\tilde{\mathbf{F}}_l^{\text{tran}}$ associated to edge $l$ appearing in equation (16) thus denotes transverse contributions from adjacent grid cells to cell $i$. More generally, this finite volume scheme is linked to a Taylor expansion of the solution in the vicinity of a grid cell. In this framework, Godunov fluxes are first order terms, high order fluxes are second order terms, and transverse fluxes correspond to cross-derivatives.

## 4 SUDDEN VELOCITY LOADING AND UNLOADING OF A HETEROGENEOUS VOLUME

Let's consider a square heterogeneous volume, of side length $2a$, with an inclusion of circular cross-section of radius $R$ centered within the volume, which mesh is shown in figure 3.
This volume in an initial natural state is suddenly loaded on its left side at time $t = 0$ with a constant first component of velocity $\bar{u}$. After time $\tau$, the applied velocity is set to zero. Normal velocity components at the top and bottom sides of the volume are set to zero while transmittive boundary conditions have been set on the right side. We assume an arbitrary heterogeneous material, that consists of an inclusion in aluminium embedded in a matrix made of steel. The analysis is carried out in plane strain. Numerical data used for computations are summarized in table 1.

<figure><img src="image_4.png" /><figcaption>Figure 3: Mesh</figcaption></figure>

<table><caption>Table 1: Numerical data used for computations.</caption><thead><tr><th>Geometry</th><th>Matrix</th><th>Inclusion</th></tr></thead><tbody><tr><td>$a = 10^{-3}$ m</td><td>$E_M = 200$ GPa</td><td>$E_I = 700$ GPa</td></tr><tr><td>$R = 5 \times 10^{-3}$ m</td><td>$\nu_M = 0.3$</td><td>$\nu_I = 0.34$</td></tr><tr><td>Loading</td><td>$\rho_M = 7800$ kg.m<sup>-3</sup></td><td>$\rho_I = 2700$ kg.m<sup>-3</sup></td></tr><tr><td>$\bar{u} = 40$ m.s<sup>-1</sup></td><td></td><td>$\sigma_{yI} = 350$ MPa</td></tr><tr><td>$\tau = 7 \times 10^{-8}$ s</td><td></td><td>$D = 2$ GPa</td></tr></tbody></table>

The matrix is assumed to remain elastic, only the inclusion undergoes viscoplastic strains. For comparison purpose, a finite element solution is obtained with the finite element code Cast3M [3], computed with an implicit time discretization and bilinear quadrangular finite elements. This numerical solution is compared to the finite volume one in figures 4, 5, and 6 at times $1 \times 10^{-7}$, $1.9 \times 10^{-7}$ et $3 \times 10^{-7}$ seconds respectively. The Courant number has been set at 0.9, the same time step being used for both numerical solutions. The comparison is performed on

<figure><img src="image_5.png" /><figcaption>Figure 4: Time $t = 1 \times 10^{-7}$ s</figcaption></figure>

the component $\sigma_{11}$ (in Pa) of the stress field and on the cumulated viscoplastic strain $p$. These two quantities are plotted along the horizontal median line, and superposed for the finite volume solution (DataSet0) and the finite element solution (DataSet1). Isovalues on the top left of figures 4 to 6 are associated to the finite volume solution while these at the top right are associated
<figure><img src="image_6.png" /><figcaption>Figure 5: Time $t = 1.9 \times 10^{-7}$s</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Time $t = 3 \times 10^{-7}$s</figcaption></figure>

to the finite element solution. One can observe in figure 4 the slot of negative normal stress formed by the prescribed horizontal velocity on the left side. Its falling front has interacted with the front interface between the matrix and the inclusion, and has been both reflected and transmitted, generating an intermediate state of higher stress (lower in absolute value) and a higher first component of velocity since the impedance of the matrix is higher than that of the inclusion. The finite element stress field shows spurious numerical oscillations in the vicinity of discontinuities, especially close to the top and bottom sides where tensile spurious stress states appear. Cumulated viscoplastic strain are almost identical for both numerical solutions.
In figure 5, one can see that the reflected wave has been transformed into tensile wave after the interaction with the increasing stress front associated to velocity unloading; then has been reflected on the left side of the volume. The cumulated viscoplastic strains propagation is delayed with respect to the propagation of the negative stress slot, which is due to viscoplasticity. At last, one can see in figure 6 that this slot has interacted with the back side of the inclusion interface, generating intermediate states of normal compression stress, while cumulated viscoplastic strains grow more importantly close to this area. Generally speaking, the finite volume solution

<figure><img src="image_8.png" /><figcaption>Figure 7: Shear stress component $\sigma_{12}$ at time $t = 3 \times 10^{-7}$s</figcaption></figure>

allows to obtain the same viscoplastic strains than these of the finite element solution without the spurious numerical oscillations on stresses obtained with the finite element solution. Figure 7 shows the comparison of the shear component of stress $\sigma_{12}$ at time $t = 3 \times 10^{-7}$s, on which these oscillations are visible.

## 5 CONCLUSIONS

It has been presented in this work the use and extension of a finite volume method for elastic-viscoplastic solids with curvilinear structured meshes of quadrangles. The formulation is based on the flux-difference splitting method, and uses the second order accurate Strang splitting approach for the computation of the right hand side of the differential system of equations. Comparison is performed with a finite element solution obtained with the finite element code Cast3M [3] on a suddenly loaded then unloaded heterogeneous volume. The comparison shows identical cumulated viscoplastic strains computed. However, the finite volume solution enables to overcome spurious numerical oscillations on the stress and velocity fields thanks to flux limiters, applied on the homogeneous part of (7), satisfying a non-increasing total variation.

## REFERENCES

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

[2] J. Lemaitre, J.L. Chaboche. *Mechanics of solid materials*. Cambridge university press, 1994.

[3] Cast3M. User's manual. www.cast3M-cea.fr, 2016.

[4] T. Belytschko, W.K. Liu, B. Moran. *Nonlinear finite elements for continua and structures*. Wiley, 2000.
[5] 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.

[6] 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.

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

[8] 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.

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