# Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids

Thomas Heuzé

► To cite this version:

Thomas Heuzé. Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids. Journal of Computational Physics, 2017, 346, pp.369-388. 10.1016/j.jcp.2017.06.027. hal-01790248

HAL Id: hal-01790248

https://hal.science/hal-01790248

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

Distributed under a Creative Commons Attribution 4.0 International License
Lax-Wendroff and TVD finite volume methods for unidimensional
thermomechanical numerical simulations of impacts on
elastic-plastic solids

Thomas 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

Abstract

We present in this work two finite volume methods for the simulation of unidimensional impact problems, both for bars and plane waves, on elastic-plastic solid media within the small strain framework. First, an extension of Lax-Wendroff to elastic-plastic constitutive models with linear and nonlinear hardenings is presented. Second, a high order TVD method based on flux-difference splitting [1] and Superbee flux limiter [2] is coupled with an approximate elastic-plastic Riemann solver for nonlinear hardenings, and follows that of Fogarty [3] for linear ones. Thermomechanical coupling is accounted for through dissipation heating and thermal softening, and adiabatic conditions are assumed. This paper essentially focuses on one-dimensional problems since analytical solutions exist or can easily be developed. Accordingly, these two numerical methods are compared to analytical solutions and to the explicit finite element method on test cases involving discontinuous and continuous solutions. This allows to study in more details their respective performance during the loading, unloading and reloading stages. Particular emphasis is also paid to the accuracy of the computed plastic strains, some differences being found according to the numerical method used. Lax-Wendoff two-dimensional discretization of a one-dimensional problem is also appended at the end to demonstrate the extensibility of such numerical scheme to multidimensional problems.

Keywords: Thermo-elastic-plastic solids, Finite volume method, Impacts, Lax-Wendroff, High order TVD method

**I. 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
[4, 5, 6] 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 [7, 8], 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 [9] 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. 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 [1, 10], has gained recently more and more interest for problems involving impacts on solid media. The characteristic structure of the set of hyperbolic equations can be accounted for by the solution of a Riemann problem at interfaces between cells, and the same order of convergence is achieved for both the velocity and stress fields [11]. Since the early work of Wilkins [12] and Trangenstein et al. [13], several authors have proposed many ways to simulate impacts on dissipative solid media such as elastic-plastic solids, these can be merely classified into eulerian approaches, generally based on a fractional-step method to treat the plasticity [14, 15, 16, 17, 18] and used for extremely high strain and strain rate problems, and lagrangian approaches [19, 20] 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. It should be noted that the eulerian approaches are often based on the so-called Maxwell-type relaxation approach [21, 15, 16, 17, 18] that rather refer to elastic-viscoplastic solids than elastic-plastic ones [9][22] since the constitutive response actually depends on time and Kuhn-Tucker complementarity conditions are not enforced at each time step. However, some of these authors [14, 23] actually enforce these conditions in that framework as both rate-dependent and rate-independent inelastic media can be written in a conservative form in the eulerian setting with a right-hand-side containing inelastic terms [24]. For Lagrangian approaches, such a fractional-step method is only possible for elastic-viscoplastic media since they exhibit inelastic phenomena occurring within a right-hand-side [25], while elastic-plastic media do not and have been so far treated within the lagrangian setting [19, 20] using classical integration of constitutive equations [9] coupled with acoustic Riemann solvers. Such a distinction of solid behaviour can be
viewed of few importance for high pressure and extremely high strain rate applications
for which shear effects are of less importance than pressure ones, making the importance
of the equation of state preeminent. However for low and moderate pressure levels for
which shear effects are of great importance, these two sets of constitutive equations lead to
different inelastic strains since they have different characteristic structures [25], and only
converge in the limit when the iso-surfaces of the viscoplastic potential in the space of
thermodynamic flux-type variables become superimposed to the indicator function of the
elastic convex. One of these two constitutive models may indeed be more in accordance
with experimental data depending on the considered metallic alloy (see [26]).

This work intend to present two finite volumes schemes for the numerical simulation of impacts on thermo-elastic-plastic solids following the framework of Generalized Standard Materials [27], within the small strain framework, both for linear and nonlinear hardenings. This paper essentially focuses on one-dimensional problems (bars and plane waves) since analytical solutions exist or can easily be developed. It allows to study in more details the performance of numerical methods during the loading, unloading and reloading stages, and to compare the computed plastic strains to analytical ones, which is generally few carried out. We present first an extension of Lax-Wendroff to elastic-plastic constitutive models. Second, a high order Total Variation Diminishing (TVD) method based on flux-difference splitting [1] and flux limiters (in particular Superbee [2]) is coupled with an approximate elastic-plastic Riemann solver for nonlinear hardenings, or follows that of Fogarty [3] for linear ones. Since the material is deformed at high strain rate, adiabatic conditions are assumed to compute the rise of temperature through dissipation heating, and thermal softening is accounted for. The paper is organized as follows. First, the elastic-plastic constitutive model and governing conserved laws, both specified for bar and plane waves, are presented in section II. Next, the two finite volume schemes are introduced in section III. At last in section IV, these two methods are compared to analytical solutions and to the explicit P1-finite element method on three test cases involving either discontinuous or both discontinuous and continuous solutions. Particular emphasis is paid to the accuracy of the computed plastic strains, of particular importance for the computation of residual stresses and distortions. Indeed, some differences are found according the numerical method used. At the end of this section, two-dimensional results of the Lax-Wendoff elastic-plastic scheme obtained on a one-dimensional problem are shown to demonstrate the extensibility of such numerical scheme to multidimensional problems.

## II. Initial Boundary Value Problem for a thermo-elastic-plastic solid in one dimension

### I. *Elastic-plastic constitutive model*

Following the local accompanying state approach [28, 22], the thermodynamical state
of the material is described by a set of state variables consisting of the strain $\varepsilon$, the tem-
perature $T$ and a set of internal state variables $Z$ which describe the evolution of internal
microstructure and stored energy due to plastic deformation and other irreversible pro-
cesses. Following the framework of Generalized Standard Materials [27], we assume the
elastic-plastic material described by a Helmholtz free-energy density potential W(ε, T, Z) and the von Mises yield function f(σ, Y) ≤ 0, where Y denotes the variables energetically conjugated to Z. Under small strain assumption, the total strain ε is additively decomposed into an elastic strain (ε^e) and a plastic strain (ε^p):

$$ \epsilon = \epsilon^e + \epsilon^p \quad (1) $$

The set of internal state variables Z consists here of the plastic strain ε^p plus two additional variables pertaining to the kinematic and isotropic hardenings. The state laws, the plastic flow rule as well as the evolution equations of internal (hardening) parameters can be derived in a classical manner [9] [22] from the free-energy and the von Mises yield function.

In the case of a bar, a unidimensional stress and strain state is considered, the set of constitutive equations read:

$$ \dot{\epsilon} = \dot{\epsilon}^e + \dot{\epsilon}^p \quad (2) $$

$$ \dot{\sigma} = E\dot{\epsilon}^e \quad (3) $$

$$ \dot{\epsilon}^p = \dot{p} \operatorname{sign}(\sigma - X) \quad (4) $$

$$ f = |\sigma - X| - \left[ \sigma_0 \left( 1 - A \frac{(T - T_0)}{T_0} \right) + R \right] \le 0 \quad (5) $$

$$ \dot{p} \ge 0; f \le 0; f\dot{p} = 0 \quad (6) $$

where the equations (2), (3), (4), (5) and (6) refer to the additive partition of the total strain rate (ε̇) into elastic (ε^e) and plastic (ε^p) parts, the elastic law, the plastic flow rule, the (von Mises) criterion and the Kuhn-Tucker complementarity conditions. In these equations, the dot (·) stands for a time rate, but rather denotes here increments since elastic-plastic equations are rate-independent. In (3), E is the Young's modulus, and σ the Cauchy axial stress component. In (4) and (6), ṗ denotes the rate of cumulated plastic strain defined for a bar as ṗ = |ε^p|. The tensile yield stress consists of two parts, the first one associated to the initial tensile yield stress σ₀, submitted to thermal softening assumed here to decrease it linearly with temperature by means of a decrease coefficient A, T₀ being here a reference temperature; and the second is the isotropic hardening variable R, assumed temperature-independent here. Moreover, the criterion (5) involves the back stress X within the von Mises norm, associated to kinematic hardening to model the Bauschinger effect. The above set of equations should be supplemented with evolution laws for hardening parameters. Linear kinematic and isotropic hardening evolution laws may be considered in the sequel:

$$ X = D\epsilon^p \quad (7) $$

$$ R = Qp \quad (8) $$

where D and Q are hardening moduli leading to the constant tangent modulus

$$ \frac{d\sigma}{d\epsilon} = \frac{E(D + Q)}{E + D + Q'} \quad (9) $$
or a nonlinear isotropic hardening following a power law:

$$R = Bp^m \qquad (10)$$

where B is a hardening parameter and m < 1 a sensitivity exponent, so that the material tangent modulus decreases with the increase of strain ($d^2\sigma/d\epsilon^2 \le 0$).

For plane waves, the strain state remains unidimensional while the stress state is three-dimensional, denoting by x, r, θ the axial, radial and hoop coordinates respectively, these are expressed as:

$$\varepsilon = \varepsilon \mathbf{e}_x \otimes \mathbf{e}_x \qquad (11)$$

$$\sigma = \sigma_x \mathbf{e}_x \otimes \mathbf{e}_x + \sigma_r (\mathbf{e}_r \otimes \mathbf{e}_r + \mathbf{e}_\theta \otimes \mathbf{e}_\theta) \qquad (12)$$

According to the plastic strain incompressibility, the plastic strain and the back stress tensors read:

$$\varepsilon^p = \varepsilon^p \left( \mathbf{e}_x \otimes \mathbf{e}_x - \frac{1}{2} (\mathbf{e}_r \otimes \mathbf{e}_r + \mathbf{e}_\theta \otimes \mathbf{e}_\theta) \right) \qquad (13)$$

$$\mathbf{X} = X \left( \mathbf{e}_x \otimes \mathbf{e}_x - \frac{1}{2} (\mathbf{e}_r \otimes \mathbf{e}_r + \mathbf{e}_\theta \otimes \mathbf{e}_\theta) \right); \quad X = D\varepsilon^p \qquad (14)$$

Hence, the equations (2), (3), (4) and (5) are replaced by the following set of equations:

$$\dot{\sigma}_r = \lambda \dot{\varepsilon} + \mu \dot{\varepsilon}^p \qquad (15)$$

$$\dot{\sigma}_x = (\lambda + 2\mu)\dot{\varepsilon} - 2\mu\dot{\varepsilon}^p \qquad (16)$$

$$\dot{\varepsilon}^p = \dot{p} \operatorname{sign} \left( \sigma_x - \sigma_r - \frac{3}{2} X \right) \qquad (17)$$

$$f = \left| \sigma_x - \sigma_r - \frac{3}{2} X \right| - (\sigma_0 + R) \left( 1 - A \frac{T - T_0}{T_0} \right) \le 0 \qquad (18)$$

where $\sigma_x \equiv \sigma$ is the axial stress component and $\sigma_r (= \sigma_\theta)$ the radial one, $\lambda$ and $\mu$ are the Lamé's parameters. Besides, one can show that the tangent modulus is changed into

$$H = \frac{2\lambda\mu + (\lambda + 2\mu)(\mu + Q + \frac{3D}{2})}{3\mu + Q + \frac{3D}{2}} \qquad (19)$$

for a plane wave with isotropic and kinematic linear hardenings.

## II. Conservation laws

A set of conservation laws that consists of the axial momentum balance and the geometric strain compatibility

$$\frac{\partial \varepsilon}{\partial t} = \frac{\partial v}{\partial x'} \qquad (20)$$

combined with the incremental stress-strain relationship, can be written as the following system of first-order partial differential equations:

$$\frac{\partial \mathbf{u}}{\partial t} + \frac{\partial \mathbf{f}(\mathbf{u})}{\partial x} = \mathbf{0} \quad \forall (x, t) \in ]0, L[ \times ]0, t_{\text{end}}] \qquad (21)$$
within a domain of length L and a duration $t_{end}$, the unknown vector $\mathbf{u}$ and the flux $\mathbf{f}(\mathbf{u})$ are defined as

$$ \mathbf{u} = \begin{bmatrix} \sigma \\ v \end{bmatrix}; \quad \mathbf{f} = \begin{bmatrix} -Hv \\ -\sigma/\rho \end{bmatrix} \qquad (22) $$

where $\rho$ is the mass density and $\sigma$ and $v$ denote the axial stress and velocity components respectively. The conservation laws (21) should be supplemented with appropriate boundary and initial conditions. For a purely elastic material, the modulus $H$ is equal to $E$ (the Young's modulus) if a bar is considered (1D stress state), or to $(\lambda + 2\mu)(\lambda$ and $\mu$ being the two Lamé's parameters) if a plane wave is considered (1D strain state). For an elastic-plastic material, the modulus $H$ stands for the respective tangent moduli.

The temperature is solution of the heat equation, assuming adiabatic conditions, written with the sole mechanical dissipation on the right hand side:

$$ \rho C \dot{T} = \mathcal{D}_{\text{mech}} \qquad (23) $$

where $C$ is the specific thermal heat capacity, and the mechanical dissipation is given by $\mathcal{D}_{\text{mech}} = (\sigma - X)\dot{\epsilon}^p - R\dot{p}$ for a bar, $\mathcal{D}_{\text{mech}} = (\sigma_x - \sigma_r - \frac{3}{2}X)\dot{\epsilon}^p - R\dot{p}$ for a plane wave.

### III. Characteristic analysis

The jacobian matrix associated to the system (21), given by:

$$ \mathbf{A} = \frac{\partial \mathbf{f}(\mathbf{u})}{\partial \mathbf{u}} = \begin{bmatrix} 0 & -H \\ -1/\rho & 0 \end{bmatrix}, \qquad (24) $$

admits real eigenvalues function of the modulus $H$. If elastic loading ($f < 0$) or unloading ($f = 0$ and $\dot{f} < 0$) is considered, classical elastic characteristic line equations and compatibility equations on characteristic lines are solution:

$$ \left\{ \begin{array}{l} dx = \pm c_d dt \\ d\sigma = \pm \rho c_d dv \end{array} \right. ; \quad c = \sqrt{\frac{H}{\rho}} \qquad (25) $$

where $H$ is equal to $E$ or $(\lambda + 2\mu)$ depending on whether bar or plane wave is considered, $c$ is the elastic celerity. If elastic-plastic loading is considered ($f = \dot{f} = 0$), plastic characteristic line equations and plastic compatibility equations are solution [25]:

$$ \left\{ \begin{array}{l} dx = \pm c_p dt \\ dv = \pm d\varphi; \quad \varphi = \int_0^\sigma \frac{d\sigma'}{\rho c_p (\sigma')} \quad ; \quad c_p = \sqrt{\frac{H}{\rho}} \end{array} \right. \qquad (26) $$

where $H$ is now the tangent modulus $d\sigma/d\epsilon$. This tangent modulus depends explicitly on the cumulated plastic strain in the case of nonlinear hardening, and may be influenced indirectly in that case by the temperature that leads to a greater plastic flow through thermal softening, and thus to a decrease of the tangent modulus and that of the plastic celerity.
### III. Finite volume spatial discretization

The finite volume methods are based on subdividing the spatial domain into grid cells [1] of length $\Delta x$ for one-dimensional problems, and defining the approximation $\mathbf{U}_i$ of $\mathbf{u}$ within the $i$th grid cell by integral averaging. These are then updated using the conservation laws (21) written in integral form on one cell, generally coupled with an explicit time integration scheme. We present here the application of the Lax-Wendroff method and a high-order Total Variation Diminishing (TVD) method using the Superbee flux limiter and flux-difference splitting [2, 1] for a thermo-elastic-plastic solid.

#### I. Lax-Wendroff

The Lax-Wendroff method is a second order accurate method, presented here in its Richtmyer two-steps version. The first step amounts to approximate $\mathbf{u}$ at the midpoint in time $t_{n+1/2} = t_n + \Delta t/2$ ($n$ referring to the time step number) and at each interface by:

$$ \mathbf{U}_{i+\frac{1}{2}}^{n+\frac{1}{2}} = \frac{\mathbf{U}_i^n + \mathbf{U}_{i+1}^n}{2} + \frac{\Delta t}{2\Delta x}(\mathbf{F}_i^n - \mathbf{F}_{i+1}^n) \quad (27) $$

where $\mathbf{U}_i^n, \mathbf{U}_{i+1}^n$ and $\mathbf{F}_i^n, \mathbf{F}_{i+1}^n$ denote respectively the approximation of $\mathbf{u}$ and $\mathbf{f}(\mathbf{u})$ in grid cells $i$ and $i+1$ at time $t_n$. The second step requires to evaluate the flux at this point

$$ \mathbf{F}_{i+\frac{1}{2}}^{n} = \mathbf{f}\left(\mathbf{U}_{i+\frac{1}{2}}^{n+\frac{1}{2}}\right) \quad (28) $$

in order to update the unknowns at time $t_{n+1}$ with the following conservative scheme:

$$ \mathbf{U}_{i}^{n+1} = \mathbf{U}_{i}^{n} - \frac{\Delta t}{\Delta x} \left( \mathbf{F}_{i+\frac{1}{2}}^{n} - \mathbf{F}_{i-\frac{1}{2}}^{n} \right) \quad (29) $$

For implementing a thermo-elastic-plastic material, a dedicated elastic-plastic Riemann solver should be considered. At each step of the Richtmyer algorithm, a prediction-correction scheme should be set in order to handle the Kuhn-Tucker complementarity conditions (6).

##### I.1. Lax-Wendroff first step

An elastic constitutive behaviour is first assumed, leading to the following elastic trial solution:

$$ \left(\sigma_{i+\frac{1}{2}}^{n+\frac{1}{2}}\right)^{\text{trial}} = \frac{\sigma_i^n + \sigma_{i+1}^n}{2} + \frac{c^2 \Delta t}{2 \Delta x} (\rho v_{i+1}^n - \rho v_i^n) \quad (30) $$

$$ v_{i+\frac{1}{2}}^{n+\frac{1}{2}} = \frac{v_i^n + v_{i+1}^n}{2} + \frac{\Delta t}{2\Delta x} \left( \frac{\sigma_{i+1}^n - \sigma_i^n}{\rho} \right) \quad (31) $$

where $c$ stands for the elastic sound speed defined by equation (25), generally set at $\Delta t/\Delta x$ if the Courant number is set at one. Notice that the elastic trial state computed with the
first step of Lax-Wendroff is identical to the method of characteristics applied with an elastic material. Then, the criterion (5) (or (18)) is computed with that trial state within the grid cells i and i + 1, provided their respective kinematic and isotropic variables values $(X, R)_{i,i+1}^n$, respectively linked to the plastic strain $(\epsilon^p)_{i,i+1}^n$ and cumulated plastic strain values $p_{i,i+1}^n$ through the evolution laws (7) and (8) or (10), and temperature values $T_{i,i+1}^n$ at time $t_n$. For plane waves, the von Mises norm arising in the criterion (18) is combined with the non-differentiated form of (15) and (16), and with the evolution law (7) so that only axial components appear:

$$ \left| \sigma_x - \sigma_r - \frac{3}{2}X \right| = \left| \frac{2\mu}{\lambda + 2\mu} \sigma - \alpha \epsilon^p \right|; \quad \alpha = \frac{\mu(3\lambda + 2\mu)}{\lambda + 2\mu} + \frac{3D}{2} \quad (32) $$

Two cases thus arise for each grid cell: either the plastic criterion is satisfied, and hence the evolution is actually elastic $\sigma_{i+1/2}^{n+1/2} = (\sigma_{i+1/2}^{n+1/2})^{\text{trial}}$, or the plastic criterion is violated $f((\sigma_{i+1/2}^{n+1/2})^{\text{trial}}) > 0$, and a plastic correction needs to be carried out.

Let's restart from the geometric strain compatibility equation (20) written in integral form on the discrete volume $[x_i, x_{i+1}] \times [t_n, t_{n+1/2}]$:

$$ \int_{t_n}^{t_{n+1/2}} \int_{x_i}^{x_{i+1}} \left( \frac{\partial \epsilon}{\partial t} - \frac{\partial v}{\partial x} \right) dxdt = 0 \quad (33) $$

Introducing the incremental stress-strain relationship yields:

$$ \int_{t_n}^{t_{n+1/2}} \int_{x_i}^{x_{i+1}} \left( \frac{1}{H} \frac{\partial \sigma}{\partial t} - \frac{\partial v}{\partial x} \right) dxdt = 0 \quad (34) $$

where $H$ denotes either the elastic or elastic-plastic tangent modulus. Separating the space integral according to cell interface coordinate $x_{i+1/2}$, and considering the modulus $H$ homogeneous in each cell, one gets:

$$ \int_{t_n}^{t_{n+1/2}} \frac{\Delta x}{2} \left( \frac{\dot{\sigma}_i}{H_i} + \frac{\dot{\sigma}_{i+1}}{H_{i+1}} \right) dt = \frac{\Delta t}{2} (v_{i+1}^n - v_i^n) \quad (35) $$

where

$$ \sigma_i = \frac{1}{\frac{\Delta x}{2}} \int_{x_i}^{x_{i+1/2}} \sigma(x,t)dx; \quad v_i^n = \frac{1}{\frac{\Delta t}{2}} \int_{t_n}^{t_{n+1/2}} v(x,t)dt. \quad (36) $$

Defining

$$ \sigma_{i+1/2}^{n+1/2} = \frac{1}{\frac{\Delta x}{2}} \int_{x_i}^{x_{i+1/2}} \sigma(x, t_{n+1/2}) dx, \quad (37) $$

it comes the update formula of the first step of Lax-Wendroff for an elastic-plastic solid:

$$ \Delta x \left( \int_{\sigma_i^n}^{\sigma_{i+1/2}^{n+1/2}} \frac{d\sigma}{H_i} + \int_{\sigma_{i+1}^n}^{\sigma_{i+1/2}^{n+1/2}} \frac{d\sigma}{H_{i+1}} \right) = \Delta t (v_{i+1}^n - v_i^n) \quad (38) $$
It is worth noting that the above results is exactly what the method of characteristics gives. Indeed, the elastic-plastic Riemann invariants eqrefPCE written in a discrete manner defined on the two integral curves passing by the states of the grid cells i and i + 1 at time $t_n$ (see figure 1), read:

$$v_{i+\frac{1}{2}}^{n+\frac{1}{2}} = v_i^n + \frac{\Delta x}{\Delta t} \int_{\sigma_i^n}^{\sigma_{i+1}^{n+\frac{1}{2}}} \frac{d\sigma}{\rho c_i^2} \quad (39)$$

$$v_{i+\frac{1}{2}}^{n+\frac{1}{2}} = v_{i+1}^{n} - \frac{\Delta x}{\Delta t} \int_{\sigma_{i+1}^{n}}^{\sigma_{i+1}^{n+\frac{1}{2}}} \frac{d\sigma}{\rho c_{i+1}^2} \quad (40)$$

where $\rho c_{i,i+1}^2 = H_{i,i+1}$. The matching of (39) and (40) yields (38). Accounting for the plastic threshold, it gives:

$$\int_{\sigma_i^n}^{\sigma_i^*} \frac{d\sigma}{\rho c_i^2} + \int_{\sigma_i^*}^{\sigma_{i+1}^{n+\frac{1}{2}}} \frac{d\sigma}{\rho c_i^2} + \int_{\sigma_{i+1}^*}^{\sigma_{i+1}^*} \frac{d\sigma}{\rho c_i^2} + \int_{\sigma_{i+1}^*}^{\sigma_{i+1}^{n+\frac{1}{2}}} \frac{d\sigma}{\rho c_{i+1}^2} = \frac{\Delta t}{\Delta x} (v_{i+1}^n - v_i^n) \quad (41)$$

where $\sigma_i^*$ and $\sigma_{i+1}^*$ stand for the tensile yield stress in grid cells i and i + 1 respectively at time $t_n$ if the trial state did not satisfy the plastic criterion, or are equal to $\sigma_i^n$ (respectively $\sigma_{i+1}^n$) if the plastic criterion has been satisfied. Indeed, at one cell interface $i + 1/2$, the

<figure><img src="image_2.png" /><figcaption>**Figure 1:** Elastic-plastic integral curves passing by the states of the grid cells i and i + 1 at time $t_n$. This configuration corresponds to that both grid cells move from an elastic to an elastic-plastic state, so that plasticity progresses from both sides.</figcaption></figure>

plasticity can develop from its left side, from its right side, or from both. The computation of integrals involving the celerities $c_{i,i+1}$ in grid cells i and i + 1 depend on the type of hardening considered. Two cases can thus be considered.

If the hardening is linear (see equations (7) and (8)), the tangent modulus is constant and given by (9). Hence, equation (41) reads:

$$\sigma_{i+\frac{1}{2}}^{n+\frac{1}{2}} = \frac{1}{\frac{1}{\rho c_i^2} + \frac{1}{\rho c_{i+1}^2}} \left( \frac{\sigma_i^*}{\rho c_i^2} + \frac{\sigma_{i+1}^*}{\rho c_{i+1}^2} + \frac{1}{\rho c^2} (\sigma_i^n + \sigma_{i+1}^n - \sigma_i^* - \sigma_{i+1}^*) + \frac{\Delta t}{\Delta x} (v_{i+1}^n - v_i^n) \right) \quad (42)$$
where $c_i$ and $c_{i+1}$ denote the sound speed of grid cells $i$ and $i+1$. These can be equal to the elastic sound speed $c$ if the evolution is elastic, or $c_p$ the (constant) plastic sound speed ($=\sqrt{(d\sigma/d\epsilon)/\rho}$) if the evolution is elastic-plastic.

If the hardening is nonlinear as for (10), the tangent modulus depends on the cumulated plastic strain and thus on the stress. Hence, the computation of integrals involving celerities $c_{i,i+1}$ should be carried out according to a stress-dependent plastic sound speed. One can show that for a power law (10), the update formula (41) leads to a nonlinear equation written on the cell interface stress $\sigma_{i+1/2}^{n+1/2}$, which is a polynomial of order $(2m-1)/(m-1)$, $m$ being the exponent arising in equation (10).

However, it can prove more convenient to proceed otherwise than solving this nonlinear equation. Indeed, the second integral arising in equations (39) and (40) can be approximated by a trapezoidal rule:

$$ \int_{\sigma_i^*}^{\sigma_{i+1/2}^{n+1/2}} \frac{d\sigma}{\rho c^2} \approx \frac{\left(\sigma_{i+1/2}^{n+1/2} - \sigma_i^*\right)}{2} \left( \frac{1}{\rho c_p^2 (p_i^n)} + \frac{1}{\rho c_p^2 \left(p_{i+1/2}^{n+1/2}\right)} \right) \quad (43) $$

which is second order accurate and thus consistent with the Lax-Wendroff order of accuracy. In the above expression, $p_{i+1/2}^{n+1/2}$ is updated with a radial return algorithm [9], provided the elastic trial solution (30), the cumulated plastic strain $p_i^n$ in grid cell $i$ at time $t_n$, the initial tensile yield stress $\sigma_0$ and hardening parameters arising in (10). Since the hardening is nonlinear, the radial return algorithm still involve the solution of an algebraic nonlinear equation. However, the routines used for finite elements can be reused.

## 1.2. Lax-Wendroff second step

The second step of the Richtmyer algorithm is also based on a prediction-correction scheme. The trial state is obtained via the conservative formula (29):

$$ (\sigma_i^{n+1})^{\text{trial}} = \sigma_i^n + \frac{\rho c^2 \Delta t}{\Delta x} \left( v_{i+\frac{1}{2}}^{n+\frac{1}{2}} - v_{i-\frac{1}{2}}^{n+\frac{1}{2}} \right) \quad (44) $$

$$ v_i^{n+1} = v_i^n + \frac{\Delta t}{\Delta x} \left( \frac{\sigma_{i+\frac{1}{2}}^{n+\frac{1}{2}} - \sigma_{i-\frac{1}{2}}^{n+\frac{1}{2}}}{\rho} \right) \quad (45) $$

This trial state is also tested with respect to the plastic criterion within each adjacent grid cell, provided their respective plastic strain $(\epsilon^p)_{i,i+1}^n$ and cumulated plastic strain value $p_{i,i+1}^n$ as well as temperature $T_{i,i+1}^n$. If the plastic criterion is violated $f\left(\left(\sigma_i^{n+1}\right)^{\text{trial}}\right) > 0$, a plastic correction is required. This correction has to be a conservative update while taking care of the plastic threshold crossing. The formula (29) is thus adapted accounting for the fact that the tangent modulus appearing in (22)$_1$ depends on the stress. Writing the conservation law (21)$_1$ in its integral form on the control volume $[x_{i-1/2}, x_{i+1/2}] \times [t_n, t_{n+1}]$,
the tangent modulus may be moved to the denominator of the left-hand side, and the time integral changed into a path integral between stress states defined at times $t_n$ and $t_{n+1}$, that can be split accounting for the threshold:

$$ \frac{\Delta x}{\Delta t} \left( \int_{\sigma_i^n}^{\sigma_i^*} \frac{d\sigma}{\rho c^2} + \int_{\sigma_i^*}^{\sigma_i^{n+1}} \frac{d\sigma}{\rho c^2} \right) = v_{i+\frac{1}{2}}^{n+\frac{1}{2}} - v_{i-\frac{1}{2}}^{n+\frac{1}{2}} \quad (46) $$

The computation of the two integrals arising in the left-hand side is performed as discussed above.

Once the stress has been updated at time $t_{n+1}$, the plastic strain or the cumulated plastic strain is updated using the criterion (5) (or (18)). At last, the temperature is updated at the end of the time step through an implicit time discretization of the heat equation (23):

$$ T_i^{n+1} = T_i^n + \frac{1}{\rho C} \left( (\sigma_i^{n+1} - X_i^{n+1}) \Delta\epsilon_i^p - R_i^{n+1} \Delta p_i \right) \quad (47) $$

## II. High-order TVD methods

The high order Total Variation Diminishing (TVD) methods enable to meet both high order of accuracy in smooth regions and a high resolution of discontinuity without any spurious oscillations where shocks arise in the solution. 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. Following [1], the conservative formula (29) can be rewritten in term of flux-difference splitting, plus some additional flux to reach higher order of accuracy:

$$ \mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{\Delta x} \left( \mathcal{A}^+ \mathbf{U}_{i-\frac{1}{2}}^n + \mathcal{A}^- \mathbf{U}_{i+\frac{1}{2}}^n \right) - \frac{\Delta t}{\Delta x} \left( \tilde{\mathbf{F}}_{i+\frac{1}{2}}^n - \tilde{\mathbf{F}}_{i-\frac{1}{2}}^n \right) \quad (48) $$

where $\mathcal{A}^+ \mathbf{U}_{i-1/2}^n = \sum_{p=1}^m (\lambda_{i-1/2}^{(p)})^+ \mathcal{W}_{i-1/2}^{(p)}$ and $\mathcal{A}^- \mathbf{U}_{i+1/2}^n = \sum_{p=1}^m (\lambda_{i+1/2}^{(p)})^- \mathcal{W}_{i+1/2}^{(p)}$ are the right-going and left-going fluctuations respectively, provided $\mathbf{U}_{i-1/2}^n = \mathbf{U}_i^n - \mathbf{U}_{i-1}^n$, $\mathcal{W}_{i-1/2}^{(p)} = \alpha_{i-1/2}^{(p)} \mathbf{K}^{(p)}$ are the waves formed of the wave strength coefficient $\alpha_{i-1/2}^{(p)}$ weight- ing the right eigenvectors $\mathbf{K}^{(p)}$ of the jacobian matrix $\mathbf{A} = \partial\mathbf{f}/\partial\mathbf{u}$, and $(\lambda_{i-1/2}^{(p)})^{\pm}$ are the characteristic speeds (the eigenvalues of $\mathbf{A}$) travelling rightward or leftward respectively. In addition, $\tilde{\mathbf{F}}_{i+1/2}^n$ and $\tilde{\mathbf{F}}_{i-1/2}^n$ stand for limited additional correction fluxes designed so that the method achieves a high order of accuracy in smooth regions and a high resolution of discontinuity in rough ones. These limited fluxes are expressed as a function of limited waves $\tilde{\mathcal{W}}_{i-1/2}^{(p)} = \tilde{\alpha}_{i-1/2}^{(p)} \mathbf{K}^{(p)}$, $1 \le p \le m$:

$$ \tilde{\mathbf{F}}_{i-\frac{1}{2}}^n = \frac{1}{2} \sum_{p=1}^m |\lambda_{i-1/2}^{(p)}| \left( 1 - \frac{\Delta t}{\Delta x} |\lambda_{i-1/2}^{(p)}| \right) \tilde{\alpha}_{i-1/2}^{(p)} \mathbf{K}^{(p)} \quad (49) $$

These waves are limited based on an upwind ratio $\theta_{i-1/2}^{(p)}$ defined for the $p$th wave as:

$$ \theta_{i-1/2}^{(p)} = \frac{\| \mathcal{W}_{i-1/2}^{(p)} \|}{\| \mathcal{W}_{i-1/2}^{(p)} \|} = \frac{\alpha_{i-1/2}^{(p)}}{\alpha_{i-1/2}^{(p)}} \quad (50) $$
where the index *I* denotes the upwind interface of that located at $x_{i-1/2}$, that is *I* equals $i-1$ if $\lambda_{i-1/2}^{(p)} > 0$, $i+1$ if $\lambda_{i-1/2}^{(p)} < 0$. This ratio may be understood as a certain measure of the local regularity of the solution, the *p*-wave strength may thus be limited as a function of this measure:

$$ \tilde{\alpha}_{i-1/2}^{(p)} = \phi(\theta_{i-1/2}^{(p)}) \alpha_{i-1/2}^{(p)} \quad (51) $$

where $\phi(\theta_{i-1/2}^{(p)})$ is a limiting function. Among many others (see [1] and [10]), Lax-Wendroff is found for $\phi(\theta) = 1$, and for instance Superbee is obtained for $\phi(\theta) = \max(0, \min(1, 2\theta), \min(2, \theta))$.

As for the Lax-Wendroff method, implementing a thermo-elastic-plastic material also requires a special elastic-plastic Riemann solver. This is also based on a prediction-correction scheme. An elastic trial solution $(\sigma_{i+1/2})^{\text{trial}}$ is first computed at each interface, solution of an elastic Riemann problem consisting of two elastic (discontinuous) waves travelling at speeds $-c$ and $c$ (see [1]). Then, this trial stress state is tested against the yield criterion in both adjacent grid cells $i$ and $i+1$, being given their respective plastic strain $(\epsilon^p)_{i,i+1}^n$, cumulated plastic strain $p_{i,i+1}^n$, and temperature values $T_{i,i+1}^n$. In each of these grid cells, either the plastic yield criterion is satisfied, and then only one elastic wave occurs, or the criterion is violated, and a plastic correction should be carried out. This correction depends actually on the type of hardening accounted in the elastic-plastic constitutive model, which drives the type of waves arising in the Riemann problem.

## II.1. Linear hardening

For a linear hardening, a plastic discontinuous wave travelling at speed $\pm c_p$ toward grid cell $i+1$ or $i$ respectively, should be added to the elastic one if the elastic trial state does not satisfy the criterion in that grid cell. The plastic wave speed is computed with the constant and known tangent modulus (9) (or (19)). In all, the characteristic spectrum solution of the elastic-plastic Riemann problem may consist of (i) two elastic waves if the trial stress state satisfies the yield criterion in both adjacent grid cells, (ii) two elastic waves plus one left or right plastic discontinuous wave if the trial stress state has violated the yield criterion in the sole left of right grid cell adjacent to the interface, or (iii) two elastic waves plus two plastic discontinuous waves if the trial stress state violates the yield criterion in both grid cells $i$ and $i+1$, as shown in figure 2. This corresponds to the discrete application of the exact solution of the elastic-plastic Riemann problem [25]. In order to apply the finite volume scheme (48), it remains to compute the strength $\alpha_{i+1/2}^{(p)}$ of each of these waves. This is based on the decomposition of $U_{i+1/2}$ on the eigenbasis of the jacobian matrix $A = \partial f / \partial u$:

$$ U_{i+1/2}^n = U_{i+1}^n - U_i^n = \sum_{p=1}^{m} \alpha_{i+1/2}^{(p)} K^{(p)} \quad (52) $$

The strength of elastic waves in yielding grid cells are first computed since the stresses in areas $U_i^*$ and $U_{i+1}^*$ are known and equal to the respective tensile yield stresses. Then
<figure><img src="image_3.png" /><figcaption>**Figure 2:** Characteristic structure of the elastic-plastic Riemann problem with linear hardening in the case of plastic flow in both grid cells.</figcaption></figure>

following [3], the strengths of plastic waves are computed by forming a system of equations with the two remaining waves. For example, if plasticity occurs from both sides, the two plastic wave strengths are computed by solving the following linear system:

$$[\mathbf{K}^{(1P)}, \mathbf{K}^{(2P)}] \begin{bmatrix} \alpha_{i+1/2}^{(1P)} \\ \alpha_{i+1/2}^{(2P)} \end{bmatrix} = \mathbf{U}_{i+1}^n - \mathbf{U}_i^n - \alpha_{i+1/2}^{(1)} \mathbf{K}^{(1)} - \alpha_{i+1/2}^{(2)} \mathbf{K}^{(2)} \quad (53)$$

If only a right or left plastic discontinuous wave is added to the characteristic structure, a similar linear system is formed from both the plastic wave and the remaining elastic wave. The fluctuations $A^+U_{i+1/2}^n$ and $A^-U_{i+1/2}^n$ are then computed, and the unknowns $U_i^{n+1}$ are hence updated at time $t_{n+1}$ through the formula (48). The thermomechanical coupling is treated as discussed in section I.2.

## II.2. Nonlinear hardening

In the case of a nonlinear decreasing hardening material, centered rarefaction waves occur if plastic flow occurs [25] as shown in figure 3, within which the solution varies continuously between the head and the tail of the rarefaction fan. However, the elastic

<figure><img src="image_4.png" /><figcaption>**Figure 3:** Characteristic structure of the elastic-plastic Riemann problem for a decreasing hardening material in the case of plastic flow in both grid cells.</figcaption></figure>

discontinuity remains between areas of states $U_i^n, U_{i+1}^n$ and the head of each fan. The solution within the rarefaction fan $\tilde{u}(\xi)$ is a self-similar one, constant along any ray
$\xi = x/t$ equal to the $p^{th}$ characteristic wavespeed $\lambda_p$, $p = (1,2)$. Using the plastic compatibility equations (26), the solution within the $p^{th}$ rarefaction fan is solution of

$$v^{(p)} = v_{i,i+1}^{n} \pm \int_{\sigma_{i,i+1}^{n}}^{\sigma^{(p)}} \frac{d\sigma'}{\rho c(\sigma')} , \quad p = 2,1 \qquad (54)$$

However, the solution within each rarefaction fan is useless for a Godunov-type method that only requires the stationary ($x/t = 0$) solution $U^*$. Hence, an approximate elastic-plastic Riemann solver can advantageously be built in order to approximate $U^*$ whatever the remaining of the Riemann solution. The rarefaction waves are replaced by discontinuous elastic and plastic waves, as for the linear hardening case (see figure 2). The elastic waves are required to generate an elastic precursor, while the plastic ones should travel at the right speed in order to compute a correct approximation $\tilde{U}^*$ of $U^*$. The plastic wave speed is computed with a tangent modulus updated implicitly by means of a radial return algorithm [9] as for the Lax-Wendroff method (see section I.1). The states $U_i^*$ and $U_{i+1}^*$ are the yielding states corresponding to the states of heads of the rarefaction fans, solved explicitly as for the linear hardening case. Actually, this approximate Riemann solver amounts to match the elastic-plastic Riemann invariants on the time step (analogous to (42) written on half-time step), provided a rectangle method with an implicit assessment of the integrand to approximate (43). So the curved part of the elastic-plastic integral curves shown in figure 1 is approximated by a straight line. It is assumed that the numerical error generated is small enough if the elastic-plastic load step within the Riemann problem and therefore the time step of the global computation is sufficiently small. The discontinuous elastic and plastic waves are next limited based on the upwind ratio (50) to compute the additional correction fluxes (49).

## IV. Numerical examples

### I. Sudden unloading of strong discontinuous loading wave in linear hardening bar

We consider a semi-infinite bar made of a linear (isotropic) hardening material, in an initial natural state, suddenly loaded on its left side at time $t = 0$ with a constant tensile stress value $\sigma^*$ sufficiently large to deform the bar plastically. After time $t_u$, the applied load is suddenly released to zero. The analytical solution of this problem can be found in [25], is summarized in figure 4 and is given below

$$\begin{align*}
\sigma_1 &= \sigma_y; & \sigma_2 &= \sigma^* \\
\sigma_4 &= \sigma_2/\beta; & \beta &= \frac{c+c_p}{c-c_p} > 1 \\
\sigma_6 &= \frac{(\sigma^* - \sigma_y)}{2} \left(\frac{c}{c_p} - 1\right) - \frac{\sigma^*}{\beta} & & (55) \\
\varepsilon_2^p &= (\sigma^* - \sigma_y) \left(\frac{1}{H} - \frac{1}{E}\right) \\
\varepsilon_4^p &= (\sigma_4 - \sigma_y) \left(\frac{1}{H} - \frac{1}{E}\right)
\end{align*}$$
<figure><img src="image_5.png" /><figcaption>**Figure 4:** Analytical solution of a semi-infinite bar suddenly plastically loaded then unloaded on its end, plotted with two strain strong-discontinuities here.</figcaption></figure>

where *H* is the tangent modulus given by (9) (with *D* set at zero), *c* and *c*<sub>p</sub> are the elastic and plastic sound speeds. A stationary strong-discontinuous strain interface is formed at each unloading cycle, the rightward plastic wave being weakened until it disappears.

## I.1. Results for an elastic-plastic material

The numerical results obtained for a steel (see table 1 for numerical values of parameters) with the Lax-Wendroff and the second-order TVD Superbee finite volume methods are compared to the analytical solution of this test case and to the results obtained with the finite element method coupled with an explicit time integrator. P1-finite elements and a lumped mass matrix are used [7]. The constitutive equations are integrated with a

<table><caption>**Table 1:** Numerical values of parameters</caption><tr><td><i>E</i> = 2 · 10<sup>11</sup> Pa</td><td><i>σ</i><sub>0</sub> = 400 · 10<sup>6</sup> Pa</td><td><i>L</i> = 6 m</td></tr><tr><td><i>ρ</i> = 7800 kg.m<sup>-3</sup></td><td><i>A</i> = 10</td><td>number of grid cells/elements = 100</td></tr><tr><td><i>C</i> = 450 J.kg<sup>-1</sup>.K<sup>-1</sup></td><td><i>T</i><sub>0</sub> = 293 K</td><td><i>t</i><sub>u</sub> = 1.2 · 10<sup>-3</sup> s</td></tr><tr><td><i>Q</i> = 10 · 10<sup>9</sup> Pa</td><td><i>σ</i><sup>*</sup> = 900 · 10<sup>6</sup> Pa</td><td></td></tr></table>

radial return algorithm [9], and the solution is computed within a larger domain than that shown in figures 5 to mimic a semi-infinite bar. No additional viscosity is added to the explicit finite element solution for comparison purpose. For the finite volume solutions, transmittive boundary conditions have been set on the right side of the bar. Comparison is performed using the values at integration points for the finite element solution, consistently with centroid values of cells for finite volumes solutions. Both the stress and plastic strain fields are compared at instants *t*₁, *t*₂, *t*₃ and *t*₄ (see figure 4), in figures 5. When the bar is loaded, an elastic precursor travels at the elastic sound speed (see figure 5(a)), followed by a plastic wave. Since the Courant number has been set at one, the three numerical methods fit perfectly the analytical solution. However, the plastic wave travels slower, numerical oscillations appear for the finite element (the largest ones) and Lax-Wendroff solutions. For the latter one, this is due to the fact that
<figure><img src="image_6.png" /><figcaption>(a) Time $t_1 = 8.89 \cdot 10^{-4}$ seconds.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Time $t_2 = 1.48 \cdot 10^{-3}$ seconds.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) Time $t_3 = 1.78 \cdot 10^{-3}$ seconds.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d) Time $t_4 = 2.67 \cdot 10^{-3}$ seconds.</figcaption></figure>

**Figure 5:** Comparison of the analytical, Finite Element (FEM), Lax-Wendroff (LW) and Superbee (SB) stress and plastic strain fields at different times.
Lax-Wendroff is not a monotone method nor a TVD one. Moreover, the finite element solution overestimates the plastic strain value closed to the boundary, though it reduces to the right values when it departs. Superbee achieves the best resolution of the plastic wave on five cells here, and compute correctly the level of plastic strain.

Then, the prescribed load is released to zero, and an unloading wave pursues the forerunning plastic loading disturbance (see figure 5(b)). Strong oscillations appears at integration points within the unloaded area in the finite element solution. Many solutions can be used to reduce these oscillations, among others (i) unload in two time steps rather than in one, (ii) reaverage at integration points the stress moved at nodes using the finite element shape functions or (iii) adding some additional numerical viscosity. Superbee achieves a proper elastic unload, superposed with the analytical solution, while Lax-Wendroff is not so bad.

When the unloading wave catches up with the plastic wave, a discontinuity of the (plastic) strain appears at this point due to the strain history difference on both sides (see figure 5(c)). Thus, a stationary discontinuous interface is generated, as well as internal reflective waves, so that the plastic strain continues to propagate rightward, but with a smaller value. The resolution of these two plastic strain fronts is quite close for these three numerical methods.

A second unload disturbance is required here to stop the progession of plastic strain (see figure 5(d)). The stress front is poorly solved by the finite element method, Lax-Wendroff has also difficulties with the left front. Only Superbee gives acceptable results after several reflexions of plastic waves.

## I.2. *Results for a thermo-elastic-plastic material*

The thermomechanical coupling is here added, with numerical parameters listed in table 1. Since on the one hand a linear hardening is considered in this work, and on the other hand the thermal part influences the mechanical one through the sole yield stress, the plastic sound speed is not affected by the temperature here, but the stress and the plastic strain actually depend on the temperature. With the numerical values of parameters listed in table 1, the stress does not change much with respect to figure 5 at the beginning, though more after many wave reflexions. The change in plastic strain is much pronounced. The thermo-elastic-plastic (TEP) numerical solutions are plotted in figure 10, and the plastic strains superposed with the (isothermal) elastic-plastic (EP) analytical solution to observe the effect of the thermomechanical coupling. The decrease of the tensile yield stress through thermal softening leads to an increase of the plastic flow, and hence to an increase of temperature (figures 6(a) and 6(b)). Recall that adiabatic conditions have been assumed so that no heat conduction effects are here accounted for. This explains why the temperature and plastic strain profiles are similar.

## II. *Sudden unloading of centered plastic loading wave in a bar with a nonlinear isotropic hardening following a power law*

We consider now the problem presented in the previous section I, but with a nonlinear isotropic hardening following a power law (10). The parameters of that law are set at
<figure><img src="image_10.png" /><figcaption>**Figure 6:** Comparison of the analytical, Finite Element (FEM), Lax-Wendroff (LW) and Superbee (SB) plastic strain and temperature fields at different times.</figcaption></figure>

$B = 770 \cdot 10^6 \text{Pa}$ and $m = 0.557$. The analytical solution of that problem is recalled in [25], and can be particularized to a nonlinear isotropic hardening following a power law, the solution of which is summarized in the Lagrange diagram in figure 7. The

<figure><img src="image_11.png" /><figcaption>**Figure 7:** Sudden unloading of centered plastic loading waves.</figcaption></figure>

solution consists of a strong discontinuous unloading disturbance pursuing a set of weak-discontinuous plastic disturbances. The strong discontinuous unloading disturbance
<figure><img src="image_12.png" /><figcaption>(a) Time $t_1 = 1.13 \cdot 10^{-3}$ seconds.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Time $t_2 = 4.15 \cdot 10^{-3}$ seconds.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) Time $t_3 = 4.62 \cdot 10^{-3}$ seconds.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) Time $t_4 = 5.33 \cdot 10^{-3}$ seconds.</figcaption></figure>

**Figure 8:** Comparison of the Finite Element (FEM), Lax-Wendroff (LW), Superbee (SB) and analytical stress and plastic strain fields at different times for nonlinear isotropic hardening following a power law.
is then absorbed by the plastic loading zone at point $A_2$ since the prescribed stress has been set sufficiently high, then the elastic-plastic boundary is reduced to a weak-discontinuous unloading boundary, across which both the stress and the velocity are continuous. Unloading elastic disturbances thus travel back and forth, these waves weaken iteratively the plastic loading disturbances within segments $A_1A_2$, $B_1B_2$, $C_1C_2$, etc.

Figure 8 shows the numerical results obtained with the explicit Finite Element Method with P1-elements, the Lax-Wendroff and the TVD Superbee finite volume methods, compared to the analytical solution of this test case. These comparisons are performed at times $t_1$ to $t_4$ (see figure 7).

The same observations hold than these made for the linear hardening case I (i) during the loading stage (figure 8(a)), up to a longer plastic raise due to the nonlinear hardening, and (ii) during the first part of the unloading stage (figure 8(b)). During the unloading stage, no stationary plastic strain discontinuity occur now, it decreases continuously and reach a plateau on segments $A_2B_1$ and $B_2C_1$ (see figure 7). The three numerical solutions of the plastic strain are very close, though the plastic strain computed with the Superbee method seem a little early (fig. 8(d)), particularly in the second and third decreases associated to lines $B_1B_2$ and to the initial plastic loading stage respectively. This may be due to the approximate elastic-plastic Riemann solver used for nonlinear hardening, which accuracy is linked to the size of the local elastic-plastic increment at each time step. However, the stress computed with Superbee in the unloading area best fits the analytical solution, the Lax-Wendroff elastic solution being slightly shifted; actually a phase error arises due to its dispersive nature [1].

### III. Plane waves in a one-dimensional finite medium made of an elastic-plastic material with linear hardening and Riemann-type initial conditions

Let's consider a one-dimensional finite medium of length $L$ with free boundaries at its two ends, made of an elastic-plastic material with linear isotropic hardening. Riemann-type initial conditions are prescribed, the velocity is prescribed to $-\bar{v}$ in the first half of the medium $x \in [0, L/2[$, and to $\bar{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 plastic flow occurs:

$$ \bar{v} = 1.2 \frac{Y_H}{\rho c} \qquad (56) $$

where $Y_H = (\lambda + 2\mu)\sigma_0/2\mu$ denotes the Hugoniot elastic limit and $c$ is given by (25). The analytical solution of that problem can be developed using [25], and is summarized in figure 9. Expressions are detailed in Appendix A. The states involved within the time duration considered are plotted with the numerical values of table 1. 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 then elastically reflected at both free ends, then interact at the middle of the medium leading to a compressive reloading, first elastically, then plastically.
<figure><img src="image_16.png" /><figcaption>**Figure 9:** Analytical solution of the one-dimensional finite medium with Riemann-type initial conditions and a linear isotropic hardening, the states are computed with the numerical values of table 1.</figcaption></figure>

### III.1. Linear isotropic hardening

Figure 10 shows the numerical results obtained considering a linear isotropic hardening with the explicit Finite Element Method with P1-elements, the Lax-Wendroff and the second-order TVD Superbee finite volume methods, compared to the analytical solution of this test case. These comparisons are performed at times $t_1$ to $t_4$ (see figure 9).

The three numerical methods fit well the analytical stress states associated to elastic and plastic waves at time $t_1$ (fig. 10(a)), up to some numerical oscillations of the finite element solution. However, big discrepancies arise in the plastic strain field: the finite element solution really overestimates the plastic strain level and oscillations appear. Refining the mesh does not make the problem go away since these oscillations become more refined, so that no plateau could be distinguishable without the knowledge of the analytical solution. Lax-Wendroff also overestimates the value of this plateau, but more weakly and without oscillations. The best result is achieved by the Superbee method that describes very well the plastic strain plateau, and solves the plastic wavefront in five cells.

Once these waves are elastically reflected at both free ends (see fig. 10(b)), they interact at the middle location of the medium in such a way that plastic flow occur on both sides of the middle, and propagates leftward and rightward (states (8) and (8'), see fig. 9). Since an elastic-plastic compressive reloading occur, the plastic strain decreases locally. The three numerical methods fits fairly well the decrease of the plastic strain and the stress minimum level, up to some oscillations, and though the spatial discretization is here too coarse to capture properly the elastic precursors. The inner plastic waves cross so that two compressive elastic and plastic wave fronts propagate in each direction from the middle (fig. 10(d)). Though the stress states are well captured by the three methods, the inner plastic strain plateau is only solved correctly by the Superbee method, the two
<figure><img src="image_17.png" /><figcaption>**Figure 10:** Comparison of the Finite Element (FEM), Lax-Wendroff (LW), Superbee (SB) and analytical stress and plastic strain fields at different times of the one-dimensional finite medium submitted to Riemann-type initial conditions with a linear isotropic hardening.</figcaption></figure>
others underestimate the negative plastic strain value.

### III.2. Comparison between isotropic and kinematic hardenings

The previous example is well designed to show the influence of the type of hardening (isotropic or kinematic) since an elastic-plastic compressive reloading occurs after an elastic-plastic tensile loading. For comparison purpose, the kinematic hardening modulus $D$ is set at $2Q/3$, $Q$ being the isotropic hardening modulus, so that the same monotone behavior is achieved with both types of hardenings. Figures 10(a) and 10(b) are hence identical in both cases. Figures 11(a) and 11(b) show a comparison of numerical solutions obtained with the Superbee method at times $t_3$ and $t_4$ considering linear isotropic or kinematic hardenings, and the corresponding analytical solutions.

Though the stress states change very little (only the difference in elastic precursor states (7) and (7'') are visible), the plastic strain levels associated to the compressive elastic-plastic reloading differ significantly. The kinematic hardening leads to a greater plastic flow since the stress state has reached the opposite of the yield surface (for compression) before that corresponding to the isotropic hardening. Figure 11 also shows that the Superbee method allows to fit perfectly the stress and plastic strain levels in both cases, though a not so fine discretization (a hundred cells) has been considered here.

<figure><img src="image_18.png" /><figcaption>**Figure 11:** Comparison of the Superbee and analytical stress and plastic strain fields of the Riemann problem with linear isotropic (ISO) or kinematic (KIN) hardening at different times in a one-dimensional medium.</figcaption></figure>
### III.3. Extension to two-dimensional discretization

In order to demonstrate the extensibility to multidimensional problems of the previous numerical schemes, first results obtained with the Lax-Wendroff elastic-plastic scheme in the two-dimensional plane strain case are presented on this one-dimensional example involving plane waves.

Let's consider a rectangular domain $\Omega: (x, y) \in ]0, L[ \times ]0, H[$, made of an elastic-plastic material with linear isotropic hardening. Free traction conditions are set at the left and right sides, while zero shear stress and transverse velocity components are prescribed on the top and bottom sides, so that combined with the plane strain assumption (zero out-of-plane velocity and strain components), plane waves can be simulated. The aforementioned Riemann-type initial conditions are prescribed on the whole domain. The mesh is made of 200 $\times$ 3 quadrangular grid cells, and the same material properties than those of section III.1 are considered. Figure 12 shows for different times $t_1$ (fig. 12(a)) to $t_4$ (fig. 13(d)) the superposed plots of the longitudinal stress and plastic strain components computed with both the one-dimensional and two-dimensional Lax-Wendroff discretizations, for a Courant number set at 0.9. For comparison purpose, the analytical solution of the fields is also plotted on these graphs.

The same trends than these explained in section III.1 are observed, although the longitudinal plastic strain distribution shows a more pronounced non-physical peak at the middle of the one-dimensional medium (see fig. 12(a)) during the first stage of the simulation, that less appear in figure 10(a) since the mesh has been refined and the Courant number has been reduced to 0.9 to ensure numerical stability of the two-dimensional numerical solution. Observe especially the perfect matching of the one-dimensional and two-dimensional Lax-Wendroff numerical solutions on this test case, during both loading and unloading stages, showing the good behavior of the method in two dimensions as well.

## V. Conclusion

Two finite volume methods were presented for the numerical simulation of impacts on one-dimensional elastic-plastic solid media, both for linear and nonlinear hardenings, and both for bars and plane waves. First, the well-known Lax-Wendroff scheme was extended to elastic-plastic behaviour. Second, an approximate Riemann solver for elastic-plastic solids with nonlinear hardening was coupled with a high order TVD scheme based on flux difference splitting and the Superbee flux limiter. These methods were tested against the finite element method and analytical solutions on three test cases.

It has been shown on three test cases that the TVD scheme performs best, particularly when discontinuous solutions appear, capturing discontinuities in few cells without oscillations. The flux limiter also allows to capture accurately the stress and plastic strain plateau, even after many wave reflexions. Continuous solutions are also well described with the approximated Riemann solver.

Without additional viscosities, the finite element method is struggling with discontinuities, especially when elastic unloadings occur, but performs well in zones where the
<figure><img src="image_19.png" /><figcaption>**Figure 12:** Comparison between the Lax-Wendroff 1D and 2D results and the analytical solution on the longitudinal stress and plastic strain components.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 12:** Comparison between the Lax-Wendroff 1D and 2D results and the analytical solution on the longitudinal stress and plastic strain components.</figcaption></figure>
solution is continuous. Moreover, the plastic strains are overestimated close to impact area
and initial discontinuity. This can lead to significant discrepancies of residual stresses and
distortions in impact-type structural problems. Lax-Wendroff appears as an intermediate
with respect to the two previous methods. It slightly oscillates in plastic plateau, solve
discontinuities in a little more cells than Superbee, but still performs well in a global
manner. At last, Lax-Wendroff has been shown to be extensible to multidimensional
problems, that of the high order TVD scheme being the purpose of an ongoing work.
Appendix A. Analytical solution of the one-dimensional medium with free ends and Riemann-type initial conditions

The analytical expressions of the stress, velocity and plastic strain fields solution of the problem described in section III are given below for a linear isotropic hardening:

$$
\begin{align*}
v_0' &= \bar{v}; v_0 = -\bar{v}; \sigma_0 = \sigma_{0'} = 0 \\
v_1 &= \frac{Y_H}{\rho c} - \bar{v} = -v_1'; \sigma_1 = \sigma_{1'} = Y_H \\
v_2 &= 0; \sigma_2 = Y_H \left(1 - \frac{c_p}{c}\right) + \rho c \bar{v}; \epsilon_2^p = \frac{(\sigma_2 - Y_H)}{2\mu} \left(\frac{\lambda + 2\mu}{H} - 1\right) \\
v_3 &= \frac{2Y_H}{\rho c} - \bar{v} = -v_3'; \sigma_3 = -\sigma_{3'} = 0 \\
v_4 &= \frac{\sigma_2}{2\rho c} + \frac{v_3}{2} = -v_4'; \sigma_4 = \sigma_{4'} = \frac{\sigma_2 - \rho c v_3}{2} \\
v_5 &= 0; \sigma_5 = \sigma_4 - \rho c v_4 \\
v_6 &= \frac{\sigma_4}{\rho c} + v_4 = -v_6'; \sigma_6 = \sigma_{6'} = 0 \\
v_7 &= \frac{\sigma_7}{\rho c} + v_6 = -v_7''; v_{7'} = \frac{\sigma_5 - \sigma_{7'}}{\rho c} + v_5 = -v_7'' \tag{A.1} \\
\sigma_7 &= \sigma_{7'} = \sigma_{7''} = \sigma_{7'''} = \sigma_2 - \frac{\lambda + 2\mu}{\mu} Q \epsilon_2^p - 2Y_H \\
v_8 &= \frac{v_7 + v_{7'}}{2} = -v_{8'}; \sigma_8 = \sigma_{8'} = \sigma_7 + \frac{\rho c (v_{7'} - v_7)}{2}; \\
\epsilon_8^p &= \epsilon_{8'}^p = \epsilon_2^p + \frac{(\sigma_8 - \sigma_7)}{2\mu} \left(\frac{\lambda + 2\mu}{H} - 1\right) \\
v_9 &= 0; \sigma_9 = \sigma_7 - \rho c_p v_{7'}; \epsilon_9^p = \epsilon_2^p + \frac{(\sigma_9 - \sigma_7)}{2\mu} \left(\frac{\lambda + 2\mu}{H} - 1\right) \\
v_{10} &= -v_{10'} = \frac{\sigma_9 - \sigma_8}{2\rho c_p} + \frac{v_8}{2}; \sigma_{10} = \sigma_{10'} = \frac{\sigma_9 + \sigma_8 - \rho c_p v_8}{2}; \\
\epsilon_{10}^p &= \epsilon_{10'}^p = \epsilon_2^p + \frac{(\sigma_{10} - \sigma_7)}{2\mu} \left(\frac{\lambda + 2\mu}{H} - 1\right) \\
v_{11} &= 0; \sigma_{11} = \sigma_{10} - \rho c_p v_{10}; \epsilon_{11}^p = \epsilon_{10'}^p + \frac{(\sigma_{11} - \sigma_{10})}{2\mu} \left(\frac{\lambda + 2\mu}{H} - 1\right)
\end{align*}
$$

where H is the tangent modulus given by (19) with D = 0. In the case of a linear kinematic hardening, the following changes occur:

$$
\begin{align}
\sigma_7 &= \sigma_2 - 2Y_H \nonumber \\
\epsilon_i^p &= \frac{2\mu(\sigma_i + Y_H)}{\alpha(\lambda + 2\mu)}; \quad i = 8, 9, 10, 11 \tag{A.2}
\end{align}
$$

where $\alpha$ is defined by (32), and $H$ is defined by (19) with $Q = 0$.
[1] R.J. Leveque. *Finite volume methods for hyperbolic problems*. Cambridge University Press, 2002.

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

[3] Tiernan Fogarty. *Finite volume methods for acoustics and elasto-plasticity with damage in a heterogeneous medium*. PhD thesis, University of Washington, 2001.

[4] V. Psyk, D. Rich, B.L. Kinsey, A.-E. Tekkaya, and M. Kleiner. Electromagnetic metal forming - A review. *Journal of Materials Processing Technology*, 211:787–829, 2011.

[5] J.D. Thomas and N. Triantafyllidis. On electromagnetic forming processes in finitely strained solids: Theory and examples. *Journal of the Mechanics and Physics of Solids*, 57:1391–1416, 2009.

[6] T. Heuzé, A. Leygue, and G. Racineux. Parametric modeling of an electromagnetic compression device with the proper generalized decomposition. *International Journal of Material Forming*, pages 1–13, 2015.

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

[8] David J. Benson. Computational methods in lagrangian and eulerian hydrocodes. *Computer Methods in Applied Mechanics and Engineering*, 99:235–394, 1992.

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

[10] E.F. Toro. *Riemann solvers and numerical methods for fluid dynamics*. Springer, 2009.

[11] C.H. Lee, A.J. Gil, and J. Bonet. Development of a cell centred upwind finite volume algorithm for a new conservation law formulation in structural dynamics. *Computers and Structures*, 118:13–38, 2013.

[12] M.L. Wilkins. *Methods in computational physics*, volume 3, chapter Calculation of elastic-plastic flow, pages 211–263. Academic Press, 1964.

[13] J.A. Trangenstein and P. Collela. A higher-order godunov method for modeling finite deformation in elastic-plastic solids. *Communications in Pure Applied mathematics*, 47:41–100, 1991.

[14] G.H. Miller and P. Collela. A high-order eulerian godunov method for elastic-plastic flow in solids. *Journal of Computational Physics*, 167:131–176, 2001.

[15] D.J. Hill, D.I. Pullin, M. Ortiz, and D.I. Meiron. An eulerian hybrid weno centered-difference solver for elastic-plastic solids. *Journal of Computational Physics*, 229:9053–9072, 2010.
[16] P.T. Barton, D. Drikakis, and 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.

[17] A. Lopez Ortega, M. Lombardini, D.I. Pullin, and D.I. 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.

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

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

[20] M. Aguirre, A.J. Gil, J. Bonet, and C.H. Lee. An upwind vertex centred finite volume solver for lagrangian solid dynamics. *Journal of Computational Physics*, 300:387–422, 2015.

[21] S.K. Godunov and E.I. Romenskii. *Elements of continuum mechanics and conservation laws*. Kluwer Academic Plenum Publishers, New York, 2003.

[22] J. Lemaitre and J.L.Chaboche. *Mechanics of Solid Materials*. Cambridge University Press, 1994.

[23] F. Wang, J.G. Glimm, J.W. Grove, and B.J. Plohr. A conservative eulerian numerical scheme for elastoplasticity and application to plate impact problems. *Impact of computing science and engineering*, 5:285–308, 1993.

[24] B.J. Plohr and D.H. Sharp. A conservative formulation for plasticity. *Advances in applied mathematics*, 13:462–493, 1992.

[25] L. Wang. *Foundations of stress waves*. Elsevier, 2007.

[26] J. Hodowany, G. Ravichandran, A.J. Rosakis, and P. Rosakis. Partition of plastic work into heat and stored energy in metals. *Experimental Mechanics*, 40:113–123, 1999.

[27] B. Halphen and Q.S. Nguyen. Sur les matériaux standards généralisés. *Journal de mécanique*, 14(1):667–688, 1975. (in French).

[28] P. Germain, Q.S. Nguyen, and P. Suquet. Continuum Thermodynamics. *Journal of Applied Mechanics*, 50:1010–1020, 1983.