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

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.


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 elasticviscoplastic 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 fluxdifference 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 conservations 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.

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 ε, the temperature 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 processes. 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 ): 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: 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),ṗ denotes the rate of cumulated plastic strain defined for a bar asṗ = |ε p |. The tensile yield stress consists of two parts, the first one associated to the initial tensile yield stress σ 0 , submitted to thermal softening assumed here to decrease it linearly with temperature by means of a decrease coefficient A, T 0 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: where D and Q are hardening moduli leading to the constant tangent modulus 4 or a nonlinear isotropic hardening following a power law: 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 σ/dε 2 ≤ 0). For plane waves, the strain state remains unidimensional while the stress state is threedimensional, denoting by x, r, θ the axial, radial and hoop coordinates respectively, these are expressed as: According to the plastic strain incompressibility, the plastic strain and the back stress tensors read: Hence, the equations (2), (3), (4) and (5) are replaced by the following set of equations: where σ x ≡ σ is the axial stress component and σ r (= σ θ ) the radial one, λ and µ are the Lamé's parameters. Besides, one can show that the tangent modulus is changed into 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 ∂ε ∂t combined with the incremental stress-strain relationship, can be written as the following system of first-order partial differential equations: ∂u ∂t 5 within a domain of length L and a duration t end , the unknown vector u and the flux f(u) are defined as where ρ is the mass density and σ 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 (λ + 2µ) (λ and µ 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: where C is the specific thermal heat capacity, and the mechanical dissipation is given by

III. Characteristic analysis
The jacobian matrix associated to the system (21), given by: admits real eigenvalues function of the modulus H. If elastic loading ( f < 0) or unloading ( f = 0 andḟ < 0) is considered, classical elastic characteristic line equations and compatibility equations on characteristic lines are solution: where H is equal to E or (λ + 2µ) depending on whether bar or plane wave is considered, c is the elastic celerity. If elastic-plastic loading is considered ( f =ḟ = 0), plastic characteristic line equations and plastic compatibility equations are solution [25]: where H is now the tangent modulus dσ/dε. 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 ∆x for one-dimensional problems, and defining the approximation U i of 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 u at the midpoint in time t n+1/2 = t n + ∆t/2 (n referring to the time step number) and at each interface by: where U n i , U n i+1 and F n i , F n i+1 denote respectively the approximation of u and f(u) in grid cells i and i + 1 at time t n . The second step requires to evaluate the flux at this point (28) in order to update the unknowns at time t n+1 with the following conservative scheme: For implementing a thermo-elastic-plastic material, a dedicated elastic-plastic Riemann solver should be considered. At each step of the Richtmyer algorithm, a predictioncorrection 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: where c stands for the elastic sound speed defined by equation (25), generally set at ∆t/∆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) n i,i+1 , respectively linked to the plastic strain (ε p ) n i,i+1 and cumulated plastic strain values p n i,i+1 through the evolution laws (7) and (8) or (10), and temperature values T n i,i+1 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: Two cases thus arise for each grid cell: either the plastic criterion is satisfied, and hence the evolution is actually elastic 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[ Introducing the incremental stress-strain relationship yields: 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: where Defining it comes the update formula of the first step of Lax-Wendroff for an elastic-plastic solid: 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 The matching of (39) and (40) yields (38). Accounting for the plastic threshold, it gives: where σ * i and σ * 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 σ n i (respectively σ n i+1 ) if the plastic criterion has been satisfied. Indeed, at one 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.
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: 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 (= (dσ/dε)/ρ) 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 σ n+1/2 i+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: which is second order accurate and thus consistent with the Lax-Wendroff order of accuracy. In the above expression, p n+1/2 i+1/2 is updated with a radial return algorithm [9], provided the elastic trial solution (30), the cumulated plastic strain p n i in grid cell i at time t n , the initial tensile yield stress σ 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.

I.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): This trial state is also tested with respect to the plastic criterion within each adjacent grid cell, provided their respective plastic strain (ε p ) n i,i+1 and cumulated plastic strain value p n i,i+1 as well as temperature T n i,i+1 . If the plastic criterion is violated f σ n+1 i trial > 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 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: 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):

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: i+1/2 are the right-going and left-going fluctuations respectively, provided are the waves formed of the wave strength coefficient α (p) i−1/2 weighting the right eigenvectors K (p) of the jacobian matrix A = ∂f/∂u, and (λ (p) i−1/2 ) ± are the characteristic speeds (the eigenvalues of A) travelling rightward or leftward respectively. In addition,F n i+1/2 andF n i−1/2 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 wavesW These waves are limited based on an upwind ratio θ (p) i−1/2 defined for the p th wave as: where the index I denotes the upwind interface of that located at x i−1/2 , that is I equals 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:α i−1/2 ) is a limiting function. Among many others (see [1] and [10]), Lax-Wendroff is found for φ(θ) = 1, and for instance Superbee is obtained for φ(θ) = max(0, min(1, 2θ), min(2, θ)).
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 predictioncorrection scheme. An elastic trial solution (σ i+1/2 ) 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 (ε p ) n i,i+1 , cumulated plastic strain p n i,i+1 , and temperature values T n i,i+1 . In each of these grid cells, either the plastic yield criterion is satified, 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 elasticplastic 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 ±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 α (p) i+1/2 of each of these waves. This is based on the decomposition of U i+1/2 on the eigenbasis of the jacobian matrix A = ∂f/∂u: 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 2: Characteristic structure of the elastic-plastic Riemann problem with linear hardening in the case of plastic flow in both grid cells.
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: 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 n i+1/2 and A − U n i+1/2 are then computed, and the unknowns U n+1 i 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 discontinuity remains between areas of states U n i , U n i+1 and the head of each fan. The solution within the rarefaction fanũ(ξ) is a self-similar one, constant along any ray ξ = x/t equal to the p th characteristic wavespeed λ p , p = (1, 2). Using the plastic compatibility equations (26), the solution within the p th rarefaction fan is solution of 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Ũ * 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 coresponding 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).

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 σ * 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 where H is the tangent modulus given by (9) (with D set at zero), c and c p are the elastic and plastic sound speeds. A stationary strong-discontinuous strain interface is formed at each unloading cycle, the righward 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 E = 2 · 10 11 Pa σ 0 = 400 · 10 6 Pa L = 6 m ρ = 7800 kg.m −3 A = 10 number of grid cells/elements = 100 C = 450 J.kg −1 .K −1 T 0 = 293 K t u = 1.2 · 10 −3 s Q = 10 · 10 9 Pa σ * = 900 · 10 6 Pa Table 1: Numerical values of parameters 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 1 , t 2 , t 3 and t 4 (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 appears for the finite element (the largest ones) and Lax-Wendroff solutions. For the latter one, this is due to the fact that      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   B = 770 · 10 6 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 solution consists of a strong discontinuous unloading disturbance pursuing a set of weakdiscontinuous plastic disturbances. The strong discontinuous unloading disturbance     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 weakdiscontinuous 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 1 A 2 , B 1 B 2 , C 1 C 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 2 B 1 and B 2 C 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 1 B 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. Riemanntype initial conditions are prescribed, the velocity is prescribed to −v in the first half of the medium x ∈ [0, L/2[, and tov in the second half x ∈]L/2, L], while the stress is considered to be zero everywhere initially. The prescribed velocity is set so that plastic flow occurs: where Y H = (λ + 2µ)σ 0 /2µ 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 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).

III.1. Linear isotropic hardening
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

22
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.

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 Ω : (x, y) ∈]0, L[×]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), planes waves can be simulated. The aforementioned Riemann-type initial conditions are prescribed on the whole domain. The mesh is made of 200 × 3 quadrangular grid cells, and the same material properties than these 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, althought 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 tests 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 were the 24 (a) Time t 1 = 3.58 · 10 −4 seconds.    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.