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

HAL
open science

# Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes

Thomas Heuzé

► To cite this version:

Thomas 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, 2018, 5 (1), 10.1186/s40323-018-0101-z. doi:10.1186/s40323-018-0101-z

HAL Id: hal-01790480

https://hal.science/hal-01790480

Submitted on 12 May 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 Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

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

# Simulation of impacts on elastic–viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes

Thomas Heuzé*

## Abstract

The flux-difference splitting finite volume method (Leveque in J Comput Phys 131:327–353, 1997); (Leveque in Finite volume methods for hyperbolic problems. Cambridge: Cambridge University Press, 2002) is here employed to perform numerical simulation of impacts on elastic–viscoplastic solids on bidimensional non-uniform quadrilateral meshes. The formulation is second order accurate in space through flux limiters, embeds the corner transport upwind method, and uses a fractional-step method to compute the relaxation operator. Elastic–viscoplastic constitutive models falling within the framework of generalized standard materials (Halphen and Nguyen in J Mech 14:667–688, 1975) in small strains are considered. Many test cases are proposed and two particular viscoplastic constitutive models are studied, on which comparisons with finite element solutions show a very good accuracy of the finite volume solutions, both on stresses and viscoplastic strains.

**Keywords:** Elastic–viscoplastic solids, Finite volume method, Flux-difference splitting, Non-uniform quadrilateral meshes, Generalized standard materials

## 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 fronts of shock waves induced in the medium. Indeed, this permits to correctly follow the path of waves and hence understand the mechanical phenomena occurring within that medium. For solid-type media, it also allows for an accurate assessment of the propagation of irreversible strains and hence of residual stresses and distortions within the structure. High speed forming processes like electromagnetic material forming [1–3] are some application examples of severe loading conditions in which the track of wave fronts is important both for understanding the development of irreversible strains in the workpiece and optimizing its final shape. Hence, these problems require numerical schemes capable to rewrite the film of history of loading undergone by any material point with sufficient precision to permit the understanding of particular physical phenomena of interest, while

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

<footnote>*Correspondence:
thomas.heuze@ec-nantes.fr
GeM, UMR 6183 CNRS, École
Centrale de Nantes, 1 rue de la
Noë, 44321 Nantes, France</footnote>
freeing oneself from any numerical disturbance that might impair that understanding. In particular, numerical schemes able to represent regular as well as discontinuous solutions are of interest; more precisely, they should meet both high orders of accuracy in regions where the solution is smooth and a high resolution of discontinuities when they occur without any numerical spurious oscillations appearing in their vicinity.

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,5], which is implemented in many industrial codes. Indeed, the finite element method is still popular in the solid mechanics community for, among others, its easy implementation of nonlinearities of partial differential equations, that is for solid-type media it enables to account easily for history-dependent constitutive equations through appropriate integration algorithms [6] and storage of internal variables at integration points in each element. However, on the one hand the amount of artificial viscosity added to numerical time integrators required to reduce the high frequency noise in the vicinity of shocks is hard to assess properly in order to remove the sole spurious oscillations, without destroying the accuracy of the numerical solution. On the other hand, finite elements do not use any feature of the characteristic structure of the set of hyperbolic equations, and is hence not the best suited method to accurately capture discontinuous solutions.

The finite volume method, initially developed for the simulation of gas dynamics [7,8], has gained recently more and more interest for problems involving impacts on solid media (see e.g. [9–17]). This family of methods show some advantages to achieve an accurate tracking of wavefronts; among others (i) the 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 same order of convergence is achieved for both the velocity and stress fields [12], and (iv) 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.

Since the early work of Wilkins [18] and Trangenstein et al. [19], several authors have proposed many ways to simulate impacts on dissipative solid media, such as elastic–plastic and elastic–viscoplastic solids, with this class of methods. These can be merely classified into Eulerian approaches, generally based on a fractional-step method to treat the irreversible processes [9, 10, 13, 17, 20] and used for extremely high strain, strain rate and pressure problems, and lagrangian approaches [14, 16, 21] that allow to follow the path of material particles and hence account for refined history-dependent constitutive equations though limited by mesh entanglement, both being coupled with an approximate Riemann or WENO solver. Eulerian approaches are written in a conservative form with a relaxation operator containing inelastic terms [22], and are often based on the so-called Maxwell-type relaxation approach [9–11, 13, 17, 23] which actually refers to an adapted version of Perzyna’s elastic–viscoplastic solids [24], in fact in its perfect viscoplasticity version since no hardening rules is generally accounted for, obtained by means of a relaxation process of stresses [11, 17]. Lagrangian approaches have been less investigated and have been so far more treated in elastoplasticity [14,16] using classical integration of constitutive
equations [6] coupled with acoustic Riemann solvers, or using simplified elastic-plastic Riemann solvers [21].

We are interested in this work in elastic-viscoplastic systems, whose study is here focused on the isothermal and linearized geometrical framework, leading to a nonhomogeneous 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. This work intends to apply the flux-difference splitting finite volume method, whose formalism has been made popular by Leveque [7,25], for the simulation of impacts on elastic-viscoplastic solid media on bidimensional non-uniform quadrilateral meshes. Its derivation for these unstructured meshes follows classical ones [26,27] for first order terms, but the process of limitation of waves required to achieve high resolution methods here accounts for different orientations between the current and upwind edges. Moreover, the approach is here derived using the class of generalized standard materials [28] (GSM) that describes a convenient framework to define thermodynamically consistent viscoplastic constitutive models, which can embed refined viscoplastic models with respect to these already used with such approach [9,17], some particular creep and hardening rules being considered in this work. The viscoplastic relaxation system is solved by means of a fractional-step method, whose convection part is solved with the flux-difference splitting formalism.

The paper is organized as follows. First, the elastic-viscoplastic constitutive model, the governing balance laws and the characteristic analysis are presented in “Elastic-viscoplastic Initial Boundary Value Problem” section. Next, the flux-difference splitting finite volume method is presented for bidimensional non-uniform quadrilateral meshes in “The flux-difference splitting finite volume method” section. “Computation of the viscoplastic part” section discusses the asymptotic limit of the elastic-viscoplastic system, and the fractional-step method used to compute the viscoplastic part of the behaviour. At last, several test cases are presented in “Applications” section, mainly conducted within the two-dimensional plane strain assumption, on which comparisons with finite element solutions allow to show the good accuracy of finite volume solutions, both on stresses and viscoplastic strains. The viscoplastic flow computed with two viscoplastic constitutive models is compared on the last example. Especially, a very simple Chaboche-type [29] viscoplastic model coupled with Prager’s linear kinematic [30] hardening, and a more refined Chaboche-Nouailhas’ [31] one coupled with the Armstrong-Frederick’s [32,33] kinematic nonlinear hardening law are considered.

## Elastic-viscoplastic Initial Boundary Value Problem

### Elastic-viscoplastic constitutive model

Following the local accompanying state approach [29,34], the thermodynamical state of the material is described by a set of state variables that consists of the linearized strain tensor $\boldsymbol{\varepsilon}$ and the temperature $T$, plus some internal state variables $Z$ which describe the evolution of internal microstructure and stored energy due to plastic deformation and other irreversible processes. Under small strain assumption, the total strain $\boldsymbol{\varepsilon}$ is additively decomposed into an elastic strain ($\boldsymbol{\varepsilon}^e$) and a plastic strain ($\boldsymbol{\varepsilon}^p$):

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

The set of internal state variables $Z$ consists here of the plastic strain $\boldsymbol{\varepsilon}^p$ plus some additional variables $\alpha_I$, $1 \le I \le N$. We assume the existence of a Helmholtz free-energy
density potential $\psi(\boldsymbol{e}, T, \boldsymbol{Z}) \equiv \psi(\boldsymbol{e} - \boldsymbol{e}^p, T, \alpha_I, 1 \le I \le N)$, concave with respect to temperature and convex with respect to other variables, from which the state laws can be derived:

$$ \boldsymbol{\sigma} = \rho \frac{\partial \psi}{\partial \boldsymbol{e}^p}; \quad s = -\frac{\partial \psi}{\partial T}; \quad A_I = \rho \frac{\partial \psi}{\partial \alpha_I}, \quad 1 \le I \le N \qquad (2) $$

where $\boldsymbol{\sigma}$, $s$ and $A_I$ refer to the Cauchy stress tensor, the entropy density and the force conjugated to variable $\alpha_I$ respectively. Thus, the mechanical dissipation reads:

$$ \mathcal{D} = \boldsymbol{\sigma} : \dot{\boldsymbol{e}}^p - A_I \cdot \dot{\alpha}_I = \boldsymbol{Y} \cdot \dot{\boldsymbol{Z}} \ge 0 \qquad (3) $$

with

$$ \boldsymbol{Y} = \{\boldsymbol{\sigma}, A_I\}; \quad \boldsymbol{Z} = \{\boldsymbol{e}^p, -\alpha_I\} \qquad (4) $$

where the vector $\boldsymbol{Y}$ denotes the thermodynamic forces conjugated to internal variables $\boldsymbol{Z}$, and the dot ($\dot{\cdot}$) applied to quantity $\square$ stands for a time rate. Following the framework of GSM [28], we assume the existence of a dissipation pseudo-potential $\Phi(\boldsymbol{Y})$, convex with respect to flux variables $\boldsymbol{Y}$, and containing the origin, such that:

$$ \dot{\boldsymbol{Z}} = \frac{\partial \Phi}{\partial \boldsymbol{Y}} \qquad (5) $$

which ensures the mechanical dissipation

$$ \mathcal{D} = \boldsymbol{Y} \cdot \frac{\partial \Phi}{\partial \boldsymbol{Y}} \ge 0 \qquad (6) $$

to be non-negative. The dissipation pseudo-potential $\Phi(\boldsymbol{Y})$ defines a family of equipotential surfaces on which any point yields the same dissipation and effective viscoplastic strain rate, and may depend on flux variables $\boldsymbol{Y}$ through a yield function $f(\boldsymbol{Y})$, so that the surface of zero potential delimits the elastic convex

$$ \mathcal{C} = \{\boldsymbol{Y}|f(\boldsymbol{Y}) \le 0\} \qquad (7) $$

The flow rule thus reads:

$$ \dot{\boldsymbol{Z}} = \frac{\partial \Phi}{\partial f} \frac{\partial f}{\partial \boldsymbol{Y}} = \dot{p} \frac{\partial f}{\partial \boldsymbol{Y}} \qquad (8) $$

where

$$ \dot{p} = \frac{\partial \Phi}{\partial f} \qquad (9) $$

denotes the viscoplastic flow intensity or effective viscoplastic strain rate, and $\partial f / \partial \boldsymbol{Y}$ is the flow direction, normal to the yield function $f$.

### Balance laws and quasi-linear form

Let's consider a continuum body $\Omega$, of boundary $\partial\Omega$, and current coordinates $\boldsymbol{x}$. The initial boundary value problem driving the motion of the solid in the small strains framework must satisfy the conservation of linear momentum, written here with neglected body forces, the geometric relationships between the strain rate $\dot{\boldsymbol{e}}$ and the velocity field $\boldsymbol{v}$, plus the constitutive equations that consist of the strain rate partition (1) combined with the elastic law (2)$_1$, the viscoplastic flow rule, and evolution equations of other internal parameters:

$$ \rho \dot{\boldsymbol{v}} = \nabla \cdot \boldsymbol{\sigma} \qquad (10) $$

$$ \dot{\boldsymbol{e}} = \frac{1}{2}(\nabla \boldsymbol{v} + (\nabla \boldsymbol{v})^T) \qquad (11) $$
$$ \dot{\sigma} = C : (\dot{\epsilon} - \dot{\epsilon}^p) \qquad (12) $$

$$ \dot{\epsilon}^p = \frac{\partial \Phi}{\partial \sigma} \qquad (13) $$

$$ -\dot{\alpha}_I = \frac{\partial \Phi}{\partial A_I}, \quad 1 \le I \le N \qquad (14) $$

which hold $\forall (\mathbf{x}, t) \in \Omega \times ]0, T]$, supplemented with appropriate initial and boundary conditions, and $\mathbf{C}$ denotes the fourth order elastic stiffness tensor in Eq. (12). In the above system, Eqs. (10) to (13) can be combined in order to write a system of balance equations of the form:

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

where $\mathcal{Q}$, $\mathcal{F}$ and $\mathcal{S}$ denote the arrays of balanced quantities, associated fluxes, and the source term respectively, defined as follows assuming a homogeneous medium:

$$ \mathcal{Q} = \begin{Bmatrix} \sigma \\ \mathbf{v} \end{Bmatrix}; \quad \mathcal{F} = \begin{Bmatrix} -\mathbf{C} : \mathbf{v} \otimes \mathbf{1} \\ -\frac{\sigma}{\rho} \end{Bmatrix}; \quad \mathcal{S} = \begin{Bmatrix} -\mathbf{C} : \frac{\partial \Phi}{\partial \sigma} \\ \mathbf{0} \end{Bmatrix} \qquad (16) $$

where $\mathbf{1}$ is the second order identity tensor. More precisely, the system (15) can be derived in full matrix form with cartesian coordinates as:

$$ \frac{\partial \mathcal{Q}}{\partial t} + \sum_{k=1}^{3} \frac{\partial \mathcal{F}_k}{\partial x_k} = \mathcal{S}(\mathcal{Q}) \qquad (17) $$

with vector components such that

$$ \begin{aligned} \mathcal{Q} &= \begin{Bmatrix} \{\boldsymbol{\sigma} \cdot \mathbf{e}_i\}_{1 \le i \le 3} \\ \mathbf{v} \end{Bmatrix}; \quad \mathcal{F}_k = \begin{Bmatrix} \{-\mathbf{C} : \mathbf{v} \otimes \mathbf{e}_k \cdot \mathbf{e}_i\}_{1 \le i \le 3} \\ -\frac{\boldsymbol{\sigma} \cdot \mathbf{e}_k}{\rho} \end{Bmatrix} = \mathcal{F} \cdot \mathbf{e}_k \\ \mathcal{S} &= \begin{Bmatrix} \{-\mathbf{C} : \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} \cdot \mathbf{e}_i\}_{1 \le i \le 3} \\ \mathbf{0} \end{Bmatrix} \end{aligned} \qquad (18) $$

where $\mathbf{e}_i$ ($1 \le i \le 3$) refer to the cartesian basis vectors, and redundant equations due to symmetry of the stress tensor are of course not considered. The term $\mathbf{C} : \mathbf{v} \otimes \mathbf{e}_k \cdot \mathbf{e}_i$ in the $k$th component of the fluxes simplifies into $C_{ijkl} \nu_j \mathbf{e}_p$, and if elastic isotropy is assumed (leading to the following expression $C_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu(\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk})$, $\lambda$ and $\mu$ denoting the Lamé's parameters), it simplifies as follows:

$$ \mathbf{C} : \mathbf{v} \otimes \mathbf{e}_k \cdot \mathbf{e}_i = \lambda \nu_k \mathbf{e}_i + \mu (\mathbf{v} \delta_{ik} + \nu_i \mathbf{e}_k) \qquad (19) $$

where $\delta_{ik}$ denotes the Kronecker delta symbol. From the system of balance equation (15), a quasi-linear form is conveniently written as follows:

$$ \frac{\partial \mathcal{Q}}{\partial t} + \sum_{k=1}^{3} \mathcal{J}_k \frac{\partial \mathcal{Q}}{\partial x_k} = \mathcal{S} \qquad (20) $$

where

$$ \mathcal{J}_k = \frac{\partial \mathcal{F}_k}{\partial \mathcal{Q}} = - \begin{bmatrix} \mathcal{O} & \mathbf{C} : \mathbf{e}_k \otimes \mathbf{1} \\ \mathbf{e}_k \cdot \mathcal{I} & \mathbf{0} \end{bmatrix} \qquad (21) $$

where $\mathcal{O}$ and $\mathcal{I}$ denote fourth order zero and identity tensors respectively, the same in full matrix form can be derived from Eq. (17):

$$ \mathcal{J}_k = \frac{\partial \mathcal{F}_k}{\partial \mathcal{Q}} = - \begin{bmatrix} \mathbf{0}_{9 \times 9} & [C_{p1qk}]_{1 \le p,q \le 3} \\ \frac{\delta_{1k} \mathbf{1}_{3 \times 3}}{\rho} & \frac{\delta_{2k} \mathbf{1}_{3 \times 3}}{\rho} & \frac{\delta_{3k} \mathbf{1}_{3 \times 3}}{\rho} \\ \frac{\delta_{1k} \mathbf{1}_{3 \times 3}}{\rho} & \frac{\delta_{2k} \mathbf{1}_{3 \times 3}}{\rho} & \frac{\delta_{3k} \mathbf{1}_{3 \times 3}}{\rho} \end{bmatrix} \qquad (22) $$
## Characteristic analysis

Since irreversible viscoplastic effects only occur in the source term (16) by means of the viscoplastic flow rule (13), the homogeneous part of this system is governed by the sole elastic part of the elastic-viscoplastic behavior, meaning that any information and in particular irreversible processes propagate along elastic characteristic curves. Let's consider the coordinate $X = x \cdot n$ and $n$ refers to any direction of propagation, along the characteristic curves $dX = \lambda_p dt$, $1 \le p \le 9$, the eigensystem

$$ \mathcal{J}_N \mathbf{K}^{(p)} = \lambda_p \mathbf{K}^{(p)}; \quad \mathbf{K}^{(p)} = \begin{bmatrix} \boldsymbol{\sigma}^{(p)} \\ \mathbf{v}^{(p)} \end{bmatrix} \qquad (23) $$

formed with the Jacobian matrix associated to the quasi-linear form (20)

$$ \mathcal{J}_N = \frac{\partial(\mathcal{F} \cdot \mathbf{n})}{\partial \mathcal{Q}} = \mathcal{J}_k n_k = - \begin{bmatrix} \mathcal{O} & \mathbf{C} : \mathbf{n} \otimes \mathbf{1} \\ \frac{\mathbf{n} \cdot \mathbf{I}}{\rho} & \mathbf{0} \end{bmatrix}, \qquad (24) $$

$\mathbf{K}^{(p)}$ denoting the right eigenvector associated to the $p$th eigenvalue (or characteristic speed) $\lambda_p$ of the Jacobian matrix $\mathcal{J}_N$, should be diagonalizable with real eigenvalues for the system (20) to be hyperbolic. The system (23) can be decoupled into the two following systems:

$$ - \mathbf{C} : \mathbf{n} \otimes \mathbf{v}^{(p)} = \lambda_p \boldsymbol{\sigma}^{(p)} \qquad (25) $$

$$ - \frac{\boldsymbol{\sigma}^{(p)} \cdot \mathbf{n}}{\rho} = \lambda_p \mathbf{v}^{(p)} \qquad (26) $$

Eliminating $\boldsymbol{\sigma}^{(p)}$ by inserting (26) into (25) yields the symmetric eigenvalue problem

$$ \mathbf{C}_{NN} \cdot \mathbf{Y}_R^{(K)} = \omega^{(K)} \mathbf{Y}_R^{(K)} \quad 1 \le K \le 3 \qquad (27) $$

where $[\mathbf{C}_{NN}]_{ij} = C_{ijkl}n_j n_l$ denotes the acoustic elastic tensor, which admits three nonzero eigenvalues and associated right eigenvectors $\{\omega^{(K)}; \mathbf{Y}_R^{(K)}\}_{1 \le K \le 3}$ thanks to the coercivity property of the elastic stiffness tensor. This leads to six distinct characteristic speeds $\lambda_p = \pm\sqrt{\omega_K/\rho}$, $1 \le p \le 6$, and right eigenvectors $\mathbf{v}^{(2K)} = \mathbf{v}^{(2K-1)} = \mathbf{Y}_R^{(K)}$. The characteristic fields can be rewritten for the non zero eigenvalues as:

$$ \left\{ \pm \sqrt{\frac{\omega^{(K)}}{\rho}}; \begin{bmatrix} -\mathbf{C} : \mathbf{n} \otimes \mathbf{Y}_R^{(K)} \\ \pm \sqrt{\frac{\omega^{(K)}}{\rho}} \mathbf{Y}_R^{(K)} \end{bmatrix} \right\}_{1 \le K \le 3} \qquad (28) $$

These are ordered in three pairs corresponding to pressure waves of speed $U_p = \sqrt{(\lambda + 2\mu)/\rho}$ and shear waves of speed $U_s = \sqrt{\mu/\rho}$ for an isotropic medium, such that $\lambda_{1,2} = \pm U_p$ and $\lambda_{3,4} = \lambda_{5,6} = \pm U_s$. The last step consists in completing the characteristic basis by finding three independent vectors associated to the kernel of the Jacobian matrix, solving equation (26) for a null eigenvalue.

Combining (23) and (20) yields the characteristic equations satisfied along elastic characteristics lines:

$$ \mathbf{L}^{(p)} \cdot \left( \frac{d\mathcal{Q}}{dt} - \mathcal{J}(\mathcal{Q}) \right) = 0 \qquad (29) $$

where $\mathbf{L}^{(p)}$ refers to the $p$th left eigenvector of the Jacobian matrix $\mathcal{J}_N$. The characteristic equations are particularized for an elastic-viscoplastic medium as

$$ \mathbf{Y}_L^{(K)} \otimes \mathbf{n} : \left( d\boldsymbol{\sigma} + \mathbf{C} : \frac{\partial\boldsymbol{\phi}}{\partial\boldsymbol{\sigma}} dt \right) \pm \rho \lambda_p \mathbf{Y}_L^{(K)} \cdot d\mathbf{v} = 0, \quad 1 \le p \le 9 \qquad (30) $$
where $Y_L^{(K)}$ refers to the Kth left eigenvector of the elastic acoustic tensor. Observe that time explicitly occurs in Eq. (30), which is a direct consequence of viscoplasticity, hence in general they do not admit first integrals and thus Riemann invariants.

## The flux-difference splitting 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 vector of balanced quantities $\mathcal{Q}$ is defined in the cell $i$ by integral averaging. Let's consider the quadrangular grid cell $i$ shown in Fig. 1, of area $|A_i|$, whose edge $s$ ($1 \le s \le 4$) joining points $P_s$ and $P_{s+1}$ is of length $L_s$, and has an outward unit normal $n_s$. Integrating a system of conservation laws, i.e. the homogeneous part of (15), on grid cell $i$ yields the following system of ordinary differential equations:

$$ \left(\frac{dU}{dt}\right)_i = -\frac{1}{|A_i|} \sum_{s=1}^{N} L_s F_s \qquad (31) $$

where $F_s$, $1 \le s \le N$, denote numerical fluxes defined at cell interfaces. The order of accuracy, the physical content and the computation cost of the finite volume method essentially rely on the definition of these numerical fluxes.

Commonly, the approach consists in exploiting the solution of Riemann problems defined at cell interfaces to determine these fluxes. In particular the stationary solution ($x/t = 0$) of the Riemann problem yields the well known Godunov's method [35]. The flux-difference splitting formulation, whose formalism has been made popular by Leveléque [7,25], amounts to rephrase the Godunov's method by splitting numerical fluxes defined at cell interfaces in terms of waves contributions, known as *fluctuations*, denoted by operators $A_k^{\pm} \Delta U_k$:

$$ \sum_{s=1}^{N} L_s F_s = \sum_{k=1}^{P} L_k A_k^{+} \Delta U_k + \sum_{l=1}^{Q} L_l A_l^{-} \Delta U_l \qquad (32) $$

where $P + Q = N$, $N$ being the number of edges of grid cell $i$. Assuming edge $k$ has left (L) and right (R) states known (Fig. 1), fluctuations read:

$$ 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)} \qquad (33) $$

and quantify the effect of all $M_w$ 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)$, computed at each edge. Each wave is weighted with a coefficient $\alpha_k^{(p)}$ determined by the

<figure><img src="image_6.png" /><figcaption>Fig. 1 Quadrangular finite volume</figcaption></figure>
projection of the jump of the state 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$ being the number of waves, equal to five for bidimensional plane strain case for example):

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

The wave strength 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 the state of grid cell $i$; this summation is performed on negative fluctuations for the Q edges having an outward normal, and on positive fluctuations for the P edges having an inward normal as shown in Fig. 2 for a non-cartesian quadrangle.

However, the above first order scheme can be improved using the class of total variation non-increasing methods [7,25], 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 when they occur. The strength of these methods relies on their ability to introduce a controlled amount of numerical viscosity locally, so that to adapt to the local regularity of the solution. One way among many others to implement them amounts to add additional fluxes to these defined from fluctuations in Eq. (33), which are limited so that a non-increasing total variation of the numerical solution be satisfied at each time step, thus ensuring that no new extrema appear that would not have already exist previously. These additional fluxes are either inward (in) or outward (out) ones, depending on the edge normal, and consist of two types of contributions

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

## High order fluxes

The first contribution allows to reach a higher order (HO) of accuracy (order two here), defined with wave strength $W_l^{(p)}$ that has been limited, hence denoted $\tilde{W}_l^{(p)} = \tilde{\alpha}_l^{(p)} K_l^{(p)}$:

$$ \tilde{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 (36) $$

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

<figure><img src="image_7.png" /><figcaption>Fig. 2 Fluctuations</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 3 Upwind edge for wave comparison</figcaption></figure>

$$ \theta_l^{(p)} = \frac{W_J^{(p)}(\mathbf{n}_l) \cdot W_l^{(p)}}{\|W_l^{(p)}\|^2} \qquad (37) $$

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 Fig. 3) if $\lambda_l^{(p)} > 0$, or the opposed edge belonging to grid cell $R$ located to the right if $\lambda_l^{(p)} < 0$. The upwind ratio (37) can be understood as a certain measure of the local regularity of the solution. However, for non-cartesian quadrangles, upwind and downwind edges do not necessarily have the same normal. Thus, the computation of the upwind ratio (37) is here performed with wave strengths computed in a same local reference frame, that of edge $l$. Wave strengths $W_J^{(p)}$ are hence computed from weighting coefficients $\alpha_J^{(p)}$ recomputed in the local frame of edge $l$:

$$ \alpha_J(\mathbf{n}_l) = \mathbf{K}^{-1}(\mathbf{n}_l) \cdot \Delta U_J \qquad (38) $$

where $\Delta U_J$ denotes the jump across edge $J$ of the state vector. Wave strengths associated to edge $J$ expressed in the frame of edge $l$ are then corrected as:

$$ \mathbf{W}_J(\mathbf{n}_l) = \mathrm{diag}(\alpha_J(\mathbf{n}_l)) \cdot \mathbf{K}(\mathbf{n}_l) = [\mathrm{diag}([\mathbf{K}(\mathbf{n}_l)]^{-1} \cdot \Delta U_J)] \cdot \mathbf{K}(\mathbf{n}_l) \qquad (39) $$

where $\mathbf{W}_J(\mathbf{n}_l)$ is the matrix made of wave strength vectors $W^{(p)}$, $1 \le p \le M_w$. The upwind ratio (37) can thus be correctly computed from (39). The wave strength 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)} \qquad (40) $$

Many limiting functions exist and permit to obtain different known finite volume schemes [36]. 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 Superbee limiter defined by $\phi(\theta) = \max(0, \min(1, 2\theta), \min(2, \theta))$ falls in this family, and is used in the following of this work. More generally, the limitation of wave strength 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 (35) enables to improve the stability of the numerical scheme, so that the Courant number can be set at one. These fluxes allow to
account for information 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 known from [37] as the Corner Transport Upwind (CTU) method. 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 Fig. 4. One focuses on the edge denoted (i) whose local frame ($\mathbf{n}_i, \mathbf{t}_i$) is shown. This edge gives rise to the computation of normal fluctuations $A_i^+ \Delta U$ and $A_i^- \Delta U$ contributing to grid cells R and L respectively. These normal fluctuations lead 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 (41)$$

where $K_j$ accounts here for the normal $\mathbf{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 Fig. 4, 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 = \sum_{p=1}^{M_w} \lambda_p^+ \beta_p K_j^{(p)} \quad (42)$$

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

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

<figure><img src="image_9.png" /><figcaption>Fig. 4 Normal and transverse fluctuations defined from edge i</figcaption></figure>
which contributes to grid cell T. The flux $\tilde{F}_i^{\text{tran}}$ associated to edge l appearing in Eq. (35) thus denotes transverse contributions from adjacent grid cells to cell i.

Gathering first order fluctuations and additional fluxes, and considering an explicit Euler time integration, the state of grid cell i is updated at time $t_{n+1}$ through the following formula:

$$ \mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{|A_i|} \left( \sum_{k=1}^{P} L_k A_k^+ \Delta \mathbf{U}_k + \sum_{l=1}^{Q} L_l 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 (44) $$

More generally, this finite volume scheme is linked, in a cartesian case, 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.

## Computation of the viscoplastic part

### Asymptotic limit of the elastic-viscoplastic relaxation system

The elastic-viscoplastic system of Eqs. (10)–(14) appears to be a system of balance equations due to the viscoplastic flow rule (13) and evolution equation (14), leading to the appearance of the source term $\mathcal{S}(\mathcal{Q})$ in the system (15). This source term can be identified to a relaxation operator of the form

$$ \mathcal{S}(\mathcal{Q}) = \frac{\mathcal{R}(\mathcal{Q})}{\tau} \quad (45) $$

where $\mathcal{R}(\mathcal{Q})$ is the relaxation term, and $\tau$ the relaxation time or stiffness parameter. The relaxation term $\mathcal{R}(\mathcal{Q}) \equiv \mathcal{R}(f(\mathcal{Q}))$ depends on $\mathcal{Q}$ through the yield function $f(\mathbf{Y})$, and is such that

$$ \frac{\partial \mathcal{R}(\mathcal{Q})}{\partial f} < 0 \quad (46) $$

leading to the relaxation of the state of the system towards an equilibrium one associated to the yield condition $f = 0$. The system (10)–(14) can thus be identified to a relaxation system [38,39], that is a system of hyperbolic conservation laws with relaxation. The latter consists of conservation laws coupled with rate equations, here the set of viscoplastic constitutive equations (12)–(14). The relaxation time $\tau$ governs the time evolution of the viscoplastic strain and internal parameters, and determines how quickly these nonconserved quantities approach their respective equilibrium values. Moreover, in the context of geometric linearization, this system of conservation laws is linear.

In the asymptotic limit $\tau \to 0$, the system (15) tends to the equilibrium system

$$ \frac{\partial \mathcal{Q}}{\partial t} + \nabla \cdot \mathcal{G}(\mathcal{Q}) = \mathbf{0} \quad (47) $$

consistent with the enforced yield (“equilibrium”) condition $f = 0$. More precisely, it is well known that viscoplasticity tends to rate-independent plasticity for the limit case of vanishing viscosity, or equivalently here for a vanishing relaxation time. In particular, Haupt [40] gave the solution of the evolution of the overstress for thermomechanical processes, and studied its limit for both slow processes and vanishing viscosity. In these cases, the viscoplastic flow rule (13) reduces to the rate-independent plastic one, within
which the expression of the plastic multiplier follows from the consistency condition $\dot{f} = 0$. Hence, in the case of plastic loading, the fluxes $\mathcal{G}$ of the equilibrium system (47) are defined from these of $\mathcal{F}$ (15) by replacing the elastic stiffness tensor by the tangent moduli $\mathbf{H}$ of the rate-independent plastic constitutive model obtained for the limit case of vanishing viscosity. Another way to show this result is to perform an expansion analogous to that of Chapman–Enskog for relaxation systems [38,39,41–43], here considered in the vicinity of the equilibrium condition $f = 0$, as shown in the one-dimensional case in [44] with the elastic–viscoplastic model of [45]. The leading term of this expansion yields (47), while the $O(\tau)$ correction gives an additional dissipative contribution.

Note that the equilibrium system (47) satisfies the subcharacteristic condition [38], meaning that the characteristic speeds of the equilibrium system (47) should be interlaced between these of the relaxation system (15). If isotropic elastic and plastic behaviour is assumed, this means that both plastic shear $c_s^P$ and pressure $c_L^P$ wave velocities should be bounded by their elastic counterparts ($c_s$ and $c_L$):

$$|c_s^P| \le |c_s| \quad \text{and} \quad |c_L^P| \le |c_L| \qquad (48)$$

This arises from the definition of tangent moduli $\mathbf{H}$, whose absolute values of eigenvalues are by construction lower than these of the elastic stiffness tensor. The subcharacteristic condition can be understood as a stability condition, and the consequences of its violation have been studied for the linear case in [41].

## Properties of numerical schemes

A system of hyperbolic conservation laws with relaxation is said to be stiff when the relaxation time $\tau$ is small compared to the time scale determined by the characteristic speeds of the system and some appropriate length scales [42], or put in another way, the wave-propagation behavior of interest occurs on a much slower time scale than the fastest time scales of the ordinary differential equation (ODE) arising from the source term. This means that if the solution is perturbed away from its equilibrium condition, then it rapidly relaxes back towards the equilibrium.

Solving stiff hyperbolic equations with relaxation can be even more challenging than solving stiff ODEs. Indeed in a stiff hyperbolic equation, the fastest reactions are often not in equilibrium everywhere. The stiffness of the relaxation operator can cause many difficulties to numerical schemes, particularly on coarse, underresolved ($\Delta t \gg \tau$) grids. Stiff source terms and underresolved numerical methods, though stable, may yield spurious nonphysical or poor numerical solutions [46]. Accordingly, numerical schemes should be designed so that to satisfy some particular properties to ensure asymptotic convergence, accuracy, and stability.

In particular, numerical schemes should (i) use coarse grids [42] that do not resolve the small relaxation time $\tau$, and still remain bounded by the Courant–Friedrichs–Lewy (CFL) stability constraint, governed by the sole convection part of the system, (ii) be asymptotic preserving [42,47] meaning that the numerical scheme should be consistent with the asymptotic limit $\tau \to 0$ for fixed $\Delta x, \Delta t$, that is the limiting scheme is a good discretization of the equilibrium system (47) even if the source term is underresolved. The numerical schemes should also (iii) be asymptotic accurate [47], that is preserve the order of accuracy in the stiff limit, (iv) strong stability preserving [47], strong stability is maintained at discrete level, and (v) well-balanced [43], preserving steady state numerically.
## Fractional-step or splitting methods

Fractional-step methods are the most commonly used approaches to solve systems of balance equations, and consist in solving alternatively a system of conservation laws with no source term, and a system of ordinary differential equations. The main idea is to take advantage of the numerical methods and mathematical backgrounds already developed both for systems of conservation laws and for stiff ODEs. The simplest (Godunov) fractional-step method takes the form

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

where $C^{(\Delta t)}$ denotes the numerical solution operator over a time step $\Delta t$ for the homogeneous part of (15), and $S^{(\Delta t)}$ that for the ODE system

$$ \frac{d\mathcal{Q}}{dt} = \mathcal{S}(\mathcal{Q}) \quad (50) $$

However, the Godunov splitting is only first order accurate, and the well-known Strang splitting [48]:

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

provides second order accuracy if each step is at least second order accurate. However, it reduces to first order for very stiff problems [42]. The implicit integration of (50) eliminates any influence of the relaxation time on the CFL condition, which thus depends solely on the convection part of (15). However, any implicit ODE solver may not yield a correct solution. For instance, although the trapezoidal method is second-order accurate and A-stable, it is only marginally stable in the stiff case [7]. Rather, the L-stability property is required [49], for example the simplest L-stable scheme is the Euler implicit method.

In this work, the Godunov splitting (49) coupled with an implicit backward Euler ODE solver will be used for stiff problems (e.g. see [9, 17]), while the Strang splitting and a backward differentiation formula at order two (BDF2) will be used for less stiff problems. However, more complex and higher order time integrators for stiff relaxation terms exist, but are not considered in this work. These are in general not based on splitting methods, and are still the purpose of current researches. Among others, the family of implicit-explicit IMEX Runge-Kutta schemes [47] allows to define high order schemes using an explicit time discretization for numerical flux and an implicit (DIRK [49]) one for the relaxation operator. Other approaches like ADER-WENO schemes [43] are also available.

## Applications

### Plane waves in a one-dimensional finite medium with Riemann-type initial conditions

### Chaboche's viscoplastic constitutive model

As a simple, guideline, viscoplastic constitutive model, a Chaboche-type one [29] is considered, derived from Perzyna's work [24], itself generalizing to three dimensions the one-dimensional formulation initially proposed by Sokolowski [50] and Malvern [45], and is based on the following expression of the pseudo-potential of dissipation:

$$ \Phi = \frac{K}{N+1} \left\langle \frac{f}{K} \right\rangle^{N+1} \quad (52) $$

where $K$ and $N$ denote viscosity and sensitivity material parameters, $\langle x \rangle = (x + |x|)/2$ the positive part of $x$, and $f$ is the yield criterion (or viscous stress) defined from the elastic convex $\mathcal{C}$ and is here associated to Mises' norm:

$$ f = \sigma_{\text{eq}}(\xi) - \sigma_y; \quad \sigma_{\text{eq}} = \sqrt{\frac{3}{2} \xi : \xi}; \quad \xi = \mathbf{s} - \mathbf{X}; \quad \mathcal{C} = \{(\sigma, \mathbf{X}) | f \le 0\} \quad (53) $$
where s is the deviatoric part of Cauchy stresses σ, X the variable defining the center of the elastic convex, and σ_y the tensile yield stress. From (9) and the potential (52), the effective viscoplastic strain rate reads:

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

This is coupled with a Prager's [30] linear kinematic hardening:

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

where D is the hardening parameter. Defining

$$ \tau = \left( \frac{K}{\sigma_y} \right)^N \qquad (56) $$

the source term of the viscoplastic relaxation system can be written as a relaxation operator (45) with the following relaxation term:

$$ \mathcal{R}(\mathcal{D}) = \left\{ \begin{array}{l} -3\mu \left( \frac{\sigma_{eq}}{\sigma_y} - 1 \right)^N \frac{\xi}{\sigma_{eq}} \\ \mathbf{0} \end{array} \right\} \qquad (57) $$

where elastic isotropy has been considered and μ is the elastic shear modulus.

### Numerical elastic-plastic asymptotic limit

Let's consider a one-dimensional finite medium of length L = 6 m with free boundaries at its two ends, made of the above elastic-viscoplastic material with linear kinematic hardening. Riemann-type initial conditions are prescribed, the velocity is prescribed to $-\dot{v}$ in the first half of the medium $x \in [0, L/2[$, and to $\dot{v}$ in the second half $x \in ]L/2, L]$, while the stress is considered to be zero everywhere initially. The prescribed velocity is set so that viscoplastic flow occurs:

$$ \dot{v} = 2 \frac{Y_H}{\rho c_L} \qquad (58) $$

where $Y_H = (\lambda + 2\mu)\sigma_y/2\mu$ denotes the Hugoniot elastic limit and $c_L$ is the elastic pressure wave velocity. The analytical solution of that problem for an elastic-plastic material with linear (kinematic and/or isotropic) hardening has been introduced in [21]. The solution first consists of two elastic and plastic waves travelling from the middle in opposite directions leading to tensile stress states. These waves are elastically reflected at both free ends, and then interact at the middle of the medium leading to a compressive reloading, first elastically, then plastically.

As $\tau$ tends to zero, the computed elastic-viscoplastic solution should tend to the elastic-plastic one. On this test case, numerical solutions computed with the Strang and Godunov splitting are compared to that computed with classical P1-finite elements. The ODE system (50) computed for the finite volume numerical solution is solved by means of an implicit backward Euler scheme for the Godunov splitting and a backward differentiation formula at order two (BDF2) for the Strang splitting. Then the viscoplastic strain is updated explicitly (with a forward Euler scheme) at the end of the time step with the viscoplastic flow rule [(combining (13), (52) and (53)]. The finite element solution is coupled with a central difference explicit time integrator, uses a lumped mass matrix [4], and the constitutive equations are integrated with a radial return algorithm [6].
Figure 5 shows a comparison between elastic-viscoplastic finite element (FEM), finite volume (FV) with Strang and Godunov splitting solutions and the elastic-plastic analytical solution. This comparison is performed on axial stress and viscoplastic strain fields at two different times. The first time corresponds to elastic and plastic waves travelling away from the center of the one-dimensional medium, while the second pertains to the compressive reloading with the same wave path. Numerical solutions in Fig. 5 are computed on a mesh that consists of 200 grid cells/elements. The time step is set so that the Courant-Friedrich-Lewy (CFL) number be equal to one. Comparison is performed using the values extracted at integration points for the finite element solution, consistently with centroid values of cells for finite volume solutions. In addition, solutions shown in Fig. 5a-c are computed with a decreasing value of the relaxation parameter $\tau$, while keeping $N$ fixed. Parameters of Table 1 correspond to results shown in Fig. 5b.

In Fig. 5a, a moderately low relaxation time is considered, viscous effects are observable since the numerical elastic-viscoplastic solutions are smoother than the elastic-plastic one. Moreover, an increased apparent tensile yield stress is observable. FEM and Godunov splitting solutions are superposed while the Strang one appears slightly in advance. In Fig.

<figure><img src="image_10.png" /></figure>
<table><caption>**Table 1 Material parameters**</caption><tr><td>E = 200 GPa</td><td>σ<sub>y</sub> = 400 MPa</td></tr><tr><td>ν = 0.3</td><td>D = 10 GPa</td></tr><tr><td>ρ = 7800 kg/m<sup>3</sup></td><td>K = 24.3 MPa s</td></tr><tr><td></td><td>N = 4.37</td></tr></table>

5b, a lower relaxation time $τ$ is set; viscous effects are now less apparent, but the same global behaviour than previously is observed. For a very low relaxation time (see Fig. 5c), the numerical solutions now conform as well as their respective possibility to the elastic-plastic analytical solution. A local overshoot occurs in the viscoplastic strains computed by the Godunov splitting scheme, and small oscillations appear: they do not subsist if the viscoplastic flow rule is solved implicitly together with the source term. The Strang solution shows an error in computing the sound speed of plastic waves, actually it shows a stiff behaviour. Indeed, this very low relaxation time leads to compute a numerical solution on an underresolved grid, since it becomes far smaller than the time step dicted by the convection part of the system. It is well-known that though Strang splitting can yield second-order accuracy for smooth solutions and non-stiff problems, it may fail to correctly compute wave speeds on underresolved grids, even with L-stable ODE solvers (see e.g. [46]).

**Energy balance**

Energy balance can also be performed in order to compare these schemes. Indeed, the central difference explicit time integrator coupled to FEM is known to conserve the total energy, while some energy is expected to be numerically dissipated for finite volume upwind schemes. The total energy is computed as the sum of kinetic $W_{kin}$ and strain $W_{int}$ energies, the latter is integrated in time with a trapezoidal rule, as discussed in [4]. These are plotted in Fig. 6 for the three schemes, computed with a CFL number set at unity, and for three values of the relaxation parameter $τ$. Figure 6a–c show three stages in their time evolution, respectively associated to the first tensile waves go, the unload return due to reflected waves at the boundaries, and the compressive reloading go from the middle. First, it can be observed that the FEM solution almost conserves the total energy as expected, while this total energy is not constant for finite volume solutions especially during the first and third stages for which viscoplastic flow occur. It is constant during the second stage because no viscoplastic flow occur, and since the CFL number is set at one, flux limiters are not activated. Therefore, the energy numerically dissipated shown in Fig. 6 arises from the fractional-step algorithms. The total energy of the FEM solution slightly changes during this stage because small oscillations appear with elastic unloading disturbances (see Fig. 5), they do not subsist for a lower CFL number (see Fig. 7).

As the relaxation parameter $τ$ decreases, the energy numerically dissipated by the Godunov splitting slightly increases, while the Strang time evolution of the total energy decreases markedly: it is associated to the appearance of its stiff behaviour (Fig. 5c).

Time evolutions of the total energy are also shown in the stiff case for lower values of the CFL number in Fig. 7a, b. First, a lower CFL number decreases the time step, therefore the Strang splitting do not show a stiff behaviour anymore for that value of relaxation time since this fast time scale can be solved. Second, the energy numerically dissipated by the Godunov splitting is reduced because the numerical error in integrating the ODE system
<figure><img src="image_11.png" /><figcaption>**Fig. 6** Time evolution of the total energy (kinetic plus strain energies) shown for different values of the relaxation time $\tau$, computed with a CFL number of one, with the elastic–viscoplastic finite element (FEM), finite volume (FV) with Strang and Godunov splitting numerical solutions. **a** $\tau = 10^{-4}$ s. **b** $\tau = 4.8 \times 10^{-6}$ s. **c** $\tau = 10^{-7}$ s</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 7** Time evolution of the total energy shown for different values of the CFL number, computed with a relaxation parameter $\tau = 1 \times 10^{-7}$ s, with the elastic–viscoplastic finite element (FEM), finite volume (FV) with Strang and Godunov splitting numerical solutions. **a** CFL = 0.4. **b** CFL = 0.04</figcaption></figure>

reduces with the time step. Third, since no viscoplastic evolution occur during the second stage, the observed numerical dissipation in that stage only arises from flux limiters which are activated due to the lower values set for the CFL number. However, this numerical dissipation appears negligible with respect to that developed during the propagation of viscoplastic flow, generated by the fractional-step algorithms.

## Convergence analysis

To complete the comparison between these schemes, convergence analyses have been performed for both discontinuous and smooth solutions, and for different values of the relaxation parameter $\tau$. These comparisons are done at the fixed time $t = 4.34 \times 10^{-4}$ s (first go of tensile waves), provided a CFL number set at one. Smooth solutions are obtained by smoothing the velocity profile of the initial condition with a portion of sinus on a length $L/4$, centered with respect to the middle of the one-dimensional medium.
Figure 8a–c show the convergence curves associated to discontinuous solutions for the three values of the relaxation parameter $\tau$. A very close convergence rate of about 0.5 for the stress and 0.8 for the velocity is observed for the three schemes in the case of initial discontinuous profile of the velocity. Indeed, because of the discontinuous solution, the source term is expected to be active only over thin regions where there are fast transients that cannot be resolved with a high accuracy [7]. Convergence curves are almost superposed for the highest value of the relaxation parameter. However, the constant in the Strang splitting convergence curves (Fig. 9b, c) becomes larger than these of the FEM and the Godunov splitting for lower values of $\tau$, yielding a bigger error. Indeed, as $\tau$ decreases, the Strang splitting shows a too stiff behaviour (Fig. 5c) on an underresolved grid (for example $\Delta x = 3 \times 10^{-2}$ m, 200 grid cells). As expected, this behaviour is improved for a lower grid size associated to a lower time step of the order of the time scale of the source term (for example $\Delta x = 3 \times 10^{-3}$ m $\rightarrow \Delta t \approx 5 \times 10^{-7}$ s $\sim \tau$). But in general we do not want to use such a fine grid, because larger time steps are preferred to save computational cost.

Figure 9a–c show the convergence curves associated to smooth initial profile of the velocity for the three values of the relaxation parameter $\tau$. Now, convergence rates of about 1.1 and 1.5 for stress and the velocity are observed, and are higher than these provided for discontinuous solutions. The FEM error appears globally to be the smallest, though that of the Godunov splitting solution is close to it in each case. However, the Strang splitting solution appears less efficient than the two others, although its error has decreased with respect to discontinuous solutions. But, as $\tau$ decreases, it shows again a too stiff behaviour.

<figure><img src="image_13.png" /><figcaption>**Fig. 8** Relative $L^2$ errors of the stress and velocity computed at time $4.34 \times 10^{-4}$ s with a CFL number set at one for a discontinuous profile of the initial velocity field. Comparison is performed for three values of the relaxation parameter $\tau$, with the FEM, the Strang and Godunov splitting finite volume solutions. **a** $\tau = 10^{-4}$ s. $\mathbf{b}$ $\tau = 4.8 \times 10^{-6}$ s. $\mathbf{c}$ $\tau = 10^{-7}$ s.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 9** Relative $L^2$ errors of the stress and velocity computed at time $4.34 \times 10^{-4}$ s with a CFL number set at one for a smooth profile of the initial velocity field. Comparison is performed for three values of the relaxation parameter $\tau$, with the FEM, the Strang and Godunov splitting finite volume solutions. **a** $\tau = 10^{-4}$ s. **b** $\tau = 4.8 \times 10^{-6}$ s. **c** $\tau = 10^{-7}$ s</figcaption></figure>

As a first draw of conclusion, for very stiff problems, the Godunov splitting should be preferred to the Strang one. For non-stiff ones, both can be used. However, it should be noticed that the present case of a plane wave is a hard test because of its one-dimensional strain state. Two-dimensional cases below yield less stiff solutions due to their multi-dimensional strain state.

## Partial impact on a plane

Let's now consider the bidimensional square domain shown in Fig. 10, submitted to an impact on a part of its top face, by means of a step function of a pressure $p$. A symmetry condition is considered on the left part, free boundaries are set at the bottom face and on the remaining part of the top face, and a perfect transmission condition is set on its right face. This problem is treated in the bidimensional plane strain framework, for which the expressions of characteristic speeds $\lambda_p$ and directions $K^{(p)}$ (28), $1 \le p \le 5$, are detailed in [7]. For this test case, a comparison between a finite volume numerical solution with Strang splitting and Q1-finite elements is considered. The finite element solution is obtained with the code Cast3M [51], computed with an implicit time stepping, and an absorbing boundary to compute the perfect transmission condition set on the right face. The previously introduced elastic-viscoplastic constitutive model with linear kinematic hardening is also considered with parameters listed in Table 1, to which are added these associated to geometry and loading of the partial impact test case, summarized in Table 2.

First, the non-uniform quadrilateral mesh shown in Fig. 11 is considered, generated with the free finite element mesh generator Gmsh [52], defined so that to achieve refined elements close to the impact area and coarse ones far from this area. Both numerical solutions are computed with this mesh, and the CFL number is set at one for both solutions.
<figure><img src="image_15.png" /><figcaption>Fig. 10 Partial impact on a plane</figcaption></figure>

<table><caption>**Table 2 Parameters of the partial impact test case**</caption><thead><tr><td>Geometry (m)</td><td>Loading</td></tr></thead><tbody><tr><td>a = 0.5</td><td rowspan="2">p = 1.5γ<sub>H</sub> = 1.05 GPa</td></tr><tr><td>L = 2</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>Fig. 11 Non-uniform quadrilateral mesh</figcaption></figure>

Figures 12 and 13 show the comparison between finite element and finite volume numerical solutions at two different times of computation. These figures consist on the one hand of numerical isovalues of the normal stress $\sigma_{22}$ and cumulated viscoplastic strains $\varepsilon_{eq} \equiv p$, for which the left-half of the domain shows the finite element solution while the right one shows the finite volume solution, and on the other hand of superposed plots of the normal stress $\sigma_{22}$ and cumulated viscoplastic strains $\varepsilon_{eq} \equiv p$ along the symmetry line. At time $2.58 \times 10^{-4}$ s (see Fig. 12), pressure and shear waves have been generated by the step pressure prescribed on a part of the top face, and have propagated downward. The curved wave front of the pressure wave is observable on stress isovalues. This front has generated and propagated viscoplastic strains. Finite element and finite volume numerical solutions look close, up to a slight undershoot of the finite element normal stress at the wave front and a slightly delayed progression of viscoplastic strains for the finite volume solution. At time $9.25 \times 10^{-4}$ s (see Fig. 13), the pressure wave has travelled one round trip from the top face, and still propagates forth. The solution is now smoother, but viscoplastic strains have developed. In particular, close to the connection between loaded and free parts of the top face of the domain, shear has been undergone leading to the spread of viscoplastic strains
<figure><img src="image_17.png" /><figcaption>**Fig. 12** Normal stress $\sigma_{22}$ and cumulated viscoplastic strains $\varepsilon_{eq} \equiv p$ at time $2.58 \times 10^{-4}$ s, computed with finite elements (left-half domain) and finite volumes (right-half domain). Superposed plots are made along the symmetry line (Additional file 1)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 13** Normal stress $\sigma_{22}$ and cumulated viscoplastic strains $\varepsilon_{eq} \equiv p$ at time $9.25 \times 10^{-4}$ s, computed with finite elements (left-half domain) and finite volumes (right-half domain). Superposed plots are made along the symmetry line</figcaption></figure>

from that area. Second, the free boundary defined at the bottom face has led to convert the initial pressure wave into a tensile one after reflection, which has generated viscoplastic strains close to that bottom face. Again, finite element and finite volume solutions look close up to few details here and there on stress and cumulated viscoplastic strains fields.
<figure><img src="image_19.png" /><figcaption>**Fig. 14** Normal stress $\sigma_{22}$ and cumulated viscoplastic strains $\varepsilon_{eq}$ at time $4.89 \times 10^{-4}$ s, computed on a 100 × 100 refined cartesian mesh with finite elements (left-half domain) and finite volumes (right-half domain). Superposed plots are made along the symmetry line</figcaption></figure>

But a clear improvement provided by finite volumes over the finite element solution can be observed by using a refined 100 × 100 quadrangle cartesian mesh, which provided the CFL number set at one, enables to properly capture the elastic discontinuity. Figure 14 shows the tensile wave reflected from the bottom side at time $9.25 \times 10^{-4}$ s, captured in one cell with the finite volume solution, while finite elements exhibit spurious numerical oscillations around this wavefront. Finite volumes benefit here from flux limiters designed to achieve a nonincreasing total variation, applied to the convection part of the system (15).

Double-notched specimen with tensile initial velocity

We consider now a double-notched specimen, whose quarter of geometry is shown in Fig. 15. This double-notched specimen is submitted to a nonzero initial tensile velocity in its extremal parts, within the shaded area (see Fig. 15). This test case is actually extracted from [17], but is here treated with the previously defined elastic–viscoplastic material with linear kinematic hardening. Due to the two planes of symmetry of the problem, the sole quarter of the double-notched specimen is meshed. Its geometrical and loading parameters are listed in Table 3. The finite volume solution is still computed with the Strang splitting, and both finite element and finite volume solutions are computed with a CFL number set at 0.9. Figures 16 and 17 show their comparison at two different times of computation. These figures are organized in the same way than for the partial impact test case, except the normal stress $\sigma_{11}$ is now considered, and the finite element solution is shown on the half-bottom part of the domain and the finite volume one is shown on the half-top part. Superposed plots are performed along the horizontal symmetry line.

The nonzero initial velocity generates a tensile wave, which is first reflected on the right free boundary, as shown in Fig. 16. The finite element solution shows spurious oscillations upstream of the left wavefront, essentially due to the implicit time stepping,
<figure><img src="image_20.png" /><figcaption>Fig. 15 One quarter of a double notched specimen, submitted to an initial tensile velocity on its right part</figcaption></figure>

<table><caption>Table 3 Parameters of the double-notched specimen test case</caption><thead><tr><td>Geometry</td><td>Loading</td></tr></thead><tbody><tr><td>a = 1.5 × 10<sup>-2</sup> m</td><td>v̄ = 40 m/s</td></tr><tr><td>L = 3.73 × 10<sup>-2</sup> m</td><td></td></tr><tr><td>h = 1 × 10<sup>-2</sup> m</td><td></td></tr><tr><td>α = π/6</td><td></td></tr><tr><td>l = 3 × 10<sup>-3</sup> m</td><td></td></tr></tbody></table>

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

<figure><img src="image_22.png" /><figcaption>Fig. 16 Normal stress σ₁₁ and cumulated viscoplastic strains ε<sub>eq</sub> at time 3.45 × 10<sup>-6</sup> s, computed with finite elements (bottom-half domain) and finite volumes (top-half domain). Superposed plots are made along the horizontal symmetry line</figcaption></figure>

while the viscoplastic strains are close for both numerical solutions and show a conical spread pattern due to the plane of symmetry. Then, the left front of the tensile wave is reflected both on the notch and left plane of symmetry, which leads the normal stress σ₁₁ to double on the symmetry line and to concentrate at the notch corner. After few waves interactions, Fig. 17 shows viscoplastic strains which have much increased at the notch corners, so does for the normal stress which shows shearing pattern due to multiple wave
<figure><img src="image_23.png" /><figcaption>**Fig. 17** Normal stress $\sigma_{11}$ and cumulated viscoplastic strains $\varepsilon_{eq}$ at time $9.54 \times 10^{-6}$ s, computed with finite elements (bottom-half domain) and finite volumes (top-half domain). Superposed plots are made along the horizontal symmetry line</figcaption></figure>

reflexions. Essentially, the two numerical solutions fit well, though the finite volume one shows less numerical spurious oscillations.

## Sudden velocity loading and unloading of a heterogeneous volume

### Chaboche's viscoplastic constitutive model

The last test case considered in this work is a square heterogeneous volume, of side length $2a$, with an inclusion of circular cross-section of radius $R$ centered within the volume. 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 $\tilde{v}$. After time $t_u$, the applied velocity is set to zero. Symmetry conditions are set at the top and bottom sides of the volume, while transmissive boundary conditions have been set on the right side. Due the symmetry of the problem, only one half of the domain is meshed as shown in Fig. 18. We assume an arbitrary heterogeneous material, that consists of an inclusion made of aluminium embedded in a matrix made of steel. The matrix is here assumed to remain elastic, only the inclusion undergoes viscoplastic strains, driven by the previously introduced Chaboche's viscoplastic constitutive model. The analysis is also carried out in plane strain, with a CFL number set at one. Numerical data used for computations are summarized in Table 4. Figures 19 and 20 show the comparison between finite element and finite volume numerical solutions at two different times of computation. The normal stress $\sigma_{11}$ and the cumulated viscoplastic strains are shown, and plotted along symmetry line.

A slot of compressive normal stress is first formed by the horizontal component of the velocity prescribed on the left side, and travels rightward in the matrix. Then, its first front interacts with the front interface between the matrix and the inclusion, and generates an
<figure><img src="image_24.png" /><figcaption>**Fig. 18** Half of a heterogeneous volume that consists of a circular inclusion embedded in a matrix, suddenly loaded then unloaded after a time $t_u$ on its left side by a prescribed velocity $\tilde{v}$</figcaption></figure>

<table><caption>**Table 4** Parameters values for the heterogeneous volume test case</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td><strong>Geometry (m)</strong></td><td><strong>Loading</strong></td></tr><tr><td>a = 10<sup>-3</sup></td><td>v̄ = 40 m/s</td></tr><tr><td>R = 5 × 10<sup>-4</sup></td><td>t<sub>u</sub> = 7 × 10<sup>-8</sup> s</td></tr><tr><td><strong>Matrix</strong></td><td><strong>Inclusion</strong></td></tr><tr><td>E<sub>M</sub> = 200 GPa</td><td>E<sub>I</sub> = 70 GPa</td></tr><tr><td>v<sub>M</sub> = 0.3</td><td>v<sub>I</sub> = 0.34</td></tr><tr><td>ρ<sub>M</sub> = 7800 kg/m<sup>3</sup></td><td>ρ<sub>I</sub> = 2700 kg/m<sup>3</sup></td></tr><tr><td></td><td>σ<sub>yI</sub> = 350 MPa</td></tr><tr><td></td><td>D<sub>I</sub> = 10 GPa</td></tr><tr><td></td><td>K<sub>I</sub> = 24.3 MPa s</td></tr><tr><td></td><td>N<sub>I</sub> = 4.37</td></tr></tbody></table>

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

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

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

<figure><img src="image_28.png" /><figcaption>**Fig. 19** Normal stress $\sigma_{11}$ and cumulated viscoplastic strains $\varepsilon_{eq}$ at time $1.52 \times 10^{-7}$ s, computed with finite elements (bottom-half domain) and finite volumes (top-half domain). Superposed plots are made along the symmetry line</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>**Fig. 20** Normal stress $\sigma_{11}$ and cumulated viscoplastic strains $\varepsilon_{eq}$ at time $3.12 \times 10^{-7}$ s, computed with finite elements (bottom-half domain) and finite volumes (top-half domain). Superposed plots are made along the symmetry line</figcaption></figure>

intermediate state of stress and velocity due to the mismatch of elastic impedances of the matrix and the inclusion. In Fig. 19, the second front of the former stress slot interacts with the generated intermediate state and yields a tensile stress wave, while the first compressive loading keeps on travelling within the inclusion, and propagates viscoplastic strains. One can observe that the finite element stress field shows spurious numerical oscillations in the vicinity of discontinuities, especially close to lateral boundaries where tensile spurious stress states appear. Cumulated viscoplastic strains are almost identical for both numerical solutions, except close to the matrix/inclusion interface. Once the tensile wave has reflected on the left side of the volume, it propagates rightward, following the initial compression slot, as shown in Fig. 20. The former compression slot has interacted with the back side of the inclusion interface, generating an important growth of cumulated viscoplastic strains close to this area. Note also that the front of the tensile slot has been curved after reflexion first on the circular matrix/inclusion interface, second on the boundary of the volume. Generally speaking, the finite volume solution 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.

### Chaboche-Nouailhas' viscoplastic constitutive model

The previous version of Chaboche's viscoplastic constitutive model can be enriched to show the generality and robustness of the method. A viscoplastic model extracted from that of Chaboche and Nouailhas [29,31] introduces an exponential term in the expression of the pseudo-potential of dissipation:

$$ \Phi = \frac{K}{\alpha(N+1)} \exp \left( \alpha \left\langle \frac{f}{K} \right\rangle^{N+1} \right) \qquad (59) $$
which allows, through a parameter $\alpha$, to saturate the overstress as the effective viscoplastic strain rate increases, which now reads:

$$ \dot{p} = \left\langle \frac{f}{K} \right\rangle^N \exp \left( \alpha \left\langle \frac{f}{K} \right\rangle^{N+1} \right) \qquad (60) $$

The saturation of the overstress driven by expression (60) is shown in Fig. 21, and superposed with the previous one (54). It amounts to consider a power law with a varying apparent exponent:

$$ N^* = \frac{d \ln \dot{p}}{d \ln f} = N + \alpha(N + 1) \left(\frac{f}{K}\right)^{N+1} \qquad (61) $$

The variation of this exponent is also shown in Fig. 21. This creep law is coupled with the nonlinear kinematic hardening law of Armstrong–Frederick [32,33]:

$$ \dot{X} = \frac{2}{3} D \dot{\varepsilon}^p - \gamma X \dot{p} \qquad (62) $$

where $\gamma$ denotes an additional hardening parameter. A tensile/compression cycle is depicted in Fig. 22 for a strain rate $\dot{\varepsilon} = 1/s$, and $\gamma = 200$. The ODE system (50) computed for finite volume numerical solutions is now solved together with the viscoplastic flow rule [combining (13), (59) and (53)] and the hardening rule (62), provided the expression (60) of the effective viscoplastic strain rate, and is still discretized with an implicit backward Euler scheme. In the bidimensional plane strain case, the nonlinear system of equations to be solved consists of eight equations associated to components 11, 22, 12, 33 of $\sigma$ and $X$. Moreover, the finite volume numerical solution is computed with a Godunov splitting, numerical values used in plots of Figs. 21 and 22, and the same CFL number set at one.

Figure 23 shows the cumulated viscoplastic strain field in the inclusion after the compression slot is passed, computed by means of finite volume method with Chaboche and Chaboche–Nouailhas’ viscoplastic constitutive models. The latter leads to much more important cumulated viscoplastic strain than the former. Indeed, for an overstress close to $10^2$ MPa, the effective cumulated strain rate has largely increased for that model (see Fig. 21), which yields a more important viscoplastic flow. As the compression slot passes within the inclusion, viscoplastic flow concentrates in the inclusion close to the rear part

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

<figure><img src="image_31.png" /><figcaption>Fig. 21 Overstress evolutions obtained with the creep laws (54) and (60), and numerical values of Table 1, and $\alpha = 10^{-3}$</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>**Fig. 22** Tensile/compression stress-strain curve at $\dot{\epsilon} = 1 \text{ s}^{-1}$</figcaption></figure>

Time: 3.53506e-7 sec.

<figure><img src="image_33.png" /><figcaption>**Fig. 23** Cumulated viscoplastic strain $\epsilon_{eq} \equiv p$ in the inclusion at time $3.35 \times 10^{-7}$ s, computed with Chaboche (right) and Chaboche–Nouailhas’ (left) viscoplastic constitutive models. Isovalues have been warped according to the magnitude of $\epsilon_{eq}$ with the same scale factor 0.2</figcaption></figure>

of the interface with the matrix (see Fig. 23). This is provided by the circular geometry of the inclusion that defines a concentrating profile as the compression wave travels rightward. In Figs. 24 and 25, isovalues of the effective viscoplastic strain rate are shown at two successive instants, as the compression slot passes. The observed profiles appear quite different in terms of the chosen viscoplastic constitutive model. In particular, the saturation of the overstress predicted by expression (60) combined with the particular rear geometry of the interface inclusion/matrix yield an effective strain rate that reaches very high numerical values on a narrow band. This small example illustrates the importance of the chosen viscoplastic constitutive model on the propagated viscoplastic strains in a dynamic process.
<figure><img src="image_34.png" /><figcaption>**Fig. 24** Effective viscoplastic strain rate $\dot{\mu}$ in the inclusion at time $1.89 \times 10^{-7}$ s, computed with Chaboche and Chaboche–Nouailhas’ viscoplastic constitutive models. Both unwarpped (left) and warpped (right) isovalues are presented, the latter being scaled with the same scale factor $10^{-8}$</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>**Fig. 25** Effective viscoplastic strain rate $\dot{\mu}$ in the inclusion at time $2.25 \times 10^{-7}$ s, computed with Chaboche and Chaboche–Nouailhas’ viscoplastic constitutive models. Both unwarpped (left) and warpped (right) isovalues are presented, the latter being scaled with the same scale factor $10^{-9}$</figcaption></figure>

Conclusion

In this work, the flux-difference splitting finite volume method [7,25] has been employed to perform numerical simulation of impacts on elastic-viscoplastic solids on bidimensional non-uniform quadrilateral meshes. The elastic-viscoplastic system of equations
identifies itself with a relaxation system with threshold, whose asymptotic limit yields an elastic–plastic system. The linearized geometrical framework considered here leads to a linear hyperbolic system with a nonlinear source term, driven by the viscoplastic part of the behaviour. This relaxation system is solved by means of a fractional-step method (Strang or Godunov splitting), whose convection part is solved with the flux-difference splitting formalism applied here to bidimensional non-uniform quadrilateral meshes. Several test cases have been proposed, and show the good accuracy of the computed finite volume solutions in terms of both stresses and viscoplastic strains with respect to finite element ones. The flux limiters used to compute the convection part enable to remove spurious numerical oscillations shown in the finite element solution close to the elastic discontinuity. In addition, two viscoplastic constitutive models have been tested to illustrate the genericity of the approach, and their influence on the viscoplastic flow has been shown on the heterogeneous volume test case. Notice that this work is straightforward extendable to three-dimensional meshes following [53] for example.

**Additional file**

**Additional file 1.** Animation of the wave propagation.

**Authors' contributions**

TH developed the idea, conducted numerical experiments and wrote the paper. The author read and approved the final manuscript.

**Acknowledgements**

None.

**Competing interests**

The author declares that he has no competing interests.

**Availability of data and materials**

Not applicable.

**Consent for publication**

Not applicable.

**Ethics approval and consent to participate**

Not applicable.

**Funding**

None.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Received: 9 November 2017 Accepted: 3 April 2018

Published online: 02 May 2018

**References**

1. Psyk V, Rich D, Kinsey BL, Tekkaya AE, Kleiner M. Electromagnetic metal forming—a review. J Mater Process Technol. 2011;211:787–829.

2. Thomas JD, Triantafyllidis N. On electromagnetic forming processes in finitely strained solids: theory and examples. J Mech Phys Solids. 2009;57:1391–416.

3. Heuzé T, Leygue A, Racineux G. Parametric modeling of an electromagnetic compression device with the proper generalized decomposition. Int J Mater Form. 2015;9:101–13.

4. Belytschko T, Liu WK, Moran B. Nonlinear finite elements for continua and structures. New York: Wiley; 2000.

5. Benson DJ. Computational methods in Lagrangian and Eulerian hydrocodes. Comput Methods Appl Mech Eng. 1992;99:235–394.

6. Simo JC, Hughes TJR. Computational inelasticity. Berlin: Springer; 1997.
7. Leveque RJ. Finite volume methods for hyperbolic problems. Cambridge: Cambridge University Press; 2002.

8. Toro EF. Riemann solvers and numerical methods for fluid dynamics. Berlin: Springer; 2013.

9. Barton PT, Drikakis D, Romenskii EI. A high-order Eulerian Godunov method for elastic-plastic flow in solids. Int J Numer Methods Eng. 2010;81:453–84.

10. Hill DJ, Pullin DI, Ortiz M, Meiron DI. An Eulerian hybrid WENO centered-difference solver for elastic-plastic solids. J Comput Phys. 2010;229:9053–72.

11. Favrie N, Gavrilyuk S. Mathematical and numerical model for nonlinear viscoplasticity. Philos Trans R Soc A. 2011;369:2864–80.

12. Lee CH, Gil AJ, Bonet J. Development of a cell centred upwind finite volume algorithm for a new conservation law formulation in structural dynamics. Comput Struct. 2013;118:13–38.

13. Ortega AL, Lombardini M, Pullin DI, Meiron DI. Numerical simulation of elastic-plastic solid mechanics using an Eulerian stretch tensor approach and HLLD Riemann solver. J Comput Phys. 2014;257:414–41.

14. Maire PH, Abgrall R, Breil J, Loubère R, Rebourcet B. A nominally second-order cell-centered Lagrangian scheme for simulating elastic-plastic flows on two dimensional unstructured grids. J Comput Phys. 2013;235:626–65.

15. Aguirre M, Gil AJ, Bonet J, Carreno AA. A vertex centred finite volume Jameson-Schmidt-Turkel (JST) algorithm for a mixed conservation formulation in solid dynamics. J Comput Phys. 2014;259:672–99.

16. Aguirre M, Gil AJ, Bonet J, Lee CH. An upwind vertex centred finite volume solver for Lagrangian solid dynamics. J Comput Phys. 2015;300:387–422.

17. Ndanou S, Favrie N, Gavrilyuk S. Multi-solid and multi-fluid diffuse interface model: applications to dynamic fracture and fragmentation. J Comput Phys. 2015;295:523–55.

18. Wilkins ML. Calculation of elastic–plastic flow. Methods in computational physics. New York: Academic Press; 1964. p. 211–63.

19. Trangenstein JA, Collela P. A higher-order Godunov method for modeling finite deformation in elastic-plastic solids. Commun Pure Appl Math. 1991;47:41–100.

20. Miller GH, Collela P. A high-order Eulerian Godunov method for elastic-plastic flow in solids. J Comput Phys. 2001;167:131–76.

21. Heuzé T. Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids. J Comput Phys. 2017;346:369–88.

22. Plohr BJ, Sharp DH. A conservative formulation for plasticity. Adv Appl Math. 1992;13:462–93.

23. Godunov SK, Romenskii EI. Elements of continuum mechanics and conservation laws. New York: Kluwer Academic Plenum Publishers; 2003.

24. Perzyna P. Fundamental problems in viscoplasticity. Adv Appl Mech. 1966;9:243–77.

25. Leveque RJ. Wave propagation algorithms for multidimensional hyperbolic systems. J Comput Phy. 1997;131:327–53.

26. Leveque RJ. High resolution finite volume methods on arbitrary grids via wave propagation. J Comput Phys. 1988;78:36–63.

27. Barth TJ, Jerpersen DC. The design and application of upwind schemes on unstructured meshes. In: AIAA 89-0366. 1989. p. 36–63.

28. Halphen B, Nguyen QS. Sur les matériaux standards généralisés. J Mech. 1975;14(1):667–88 (in French).

29. Lemaitre J, Chaboche JL. Mechanics of solid materials. Cambridge: Cambridge University Press; 1994.

30. Prager W. Recent developments in the mathematical theory of plasticity. J Appl Phys. 1949;20:235–41.

31. Nouailhas D. Unified modelling of cyclic viscoplasticity: application to austenitic stainless steels. Int J Plast. 1989;5:501–20.

32. Frederick CO, Armstrong PJ. A mathematical representation of the multiaxial Bauschinger effect. Berkeley: Central Electricity Generating Board; 1966.

33. Frederick CO, Armstrong PJ. A mathematical representation of the multiaxial Bauschinger effect. Mater High Temp. 2007;24:1–26.

34. Germain P, Nguyen QS, Suquet P. Continuum thermodynamics. J Appl Mech. 1983;50:1010–20.

35. Godunov SK. Finite difference method for the numerical computation of discontinuous solutions of the equations of fluid dynamics. Matematicheskii Sbornik. 1959;47:271–306.

36. Sweby PK. High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM J Numer Anal. 1984;21:995–1011.

37. Collela P. Multidimensional upwind methods for hyperbolic conservation laws. J Comput Phys. 1990;87:171–200.

38. Liu TP. Hyperbolic conservation laws with relaxation. Commun Math Phys. 1987;108:153–75.

39. Chen GQ, Levermore CD, Liu TP. Hyperbolic conservation laws with stiff relaxation terms and entropy. Commun Pure Appl Math. 1994;47(6):787–830.

40. Haupt P, Kamlah M, Tsakmakis C. On the thermodynamics of rate-independent plasticity as an asymptotic limit of viscoplasticity for slow processes. In: Finite inelastic deformations—theory and applications. Berlin: Springer; 1992. p. 107–116.

41. Leveque RJ, Jinghua W. Nonlinear hyperbolic problems: theoretical, applied, and computational aspects. In: Donato A, Oliveri F, editors. A linear hyperbolic system with stiff source terms. Berlin: Springer; 1992. p. 401–8.

42. Jin S. Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms. J Comput Phys. 1995;122:51–67.

43. Dumbser M, Eaux C, Toro EF. Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws. J Comput Phys. 2008;227:3971–4001.

44. Tzavaras A. Viscosity and relaxation approximation for hyperbolic systems of conservation laws. In: Lecture Notes in Computational Science and Engineering, vol. 5. 1999. p. 73–122.

45. Malvern LE. The propagation of longitudinal waves of plastic deformation in a bar of material exhibiting a strain-rate effect. J Appl Mech Trans ASME. 1951;18(2):203–8.

46. Leveque RJ, Yee HC. A study of numerical methods for hyperbolic conservation laws with stiff source terms. J Comput Phys. 1990;86:187–210.
47. Pareschi L, Russo G. Implicit-explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation. J Sci Comput. 2005;25:129–55.

48. Strang G. On the construction and comparison of difference schemes. SIAM J Numer Anal. 1968;5:506–17.

49. Asher UM, Petzold LR. Computer methods for ordinary differential equations and differential-algebraic equations. Philadelphia: SIAM; 1998.

50. Sokolowskii VV. Propagation of elasto-visco-plastic waves in bars. Prik Mat Mekh. 1948;12(2):261 (in Russian).

51. Cast3M. Finite element code developed by the CEA, (French atomic energy commission). 2018. http://www-cast3m.cea.fr. Accessed 21 Apr 2018.

52. Remacle JF, Lambrechts J, Seny B, Marchandise E, Johnen A, Geuzaine C. Blossom-Quad: a non-uniform quadrilateral mesh generator using a minimum-cost perfect-matching algorithm. Int J Numer Methods Eng. 2011;89:1102–19.

53. Langseth JO, Leveque RJ. A wave propagation method for three-dimensional hyperbolic conservation laws. J Comput Phys. 2000;165:126–66.

Submit your manuscript to a SpringerOpen® journal and benefit from:
▶ Convenient online submission
▶ Rigorous peer review
▶ Open access: articles freely available online
▶ High visibility within the field
▶ Retaining the copyright to your article

Submit your next manuscript at ► springeropen.com