Lax-Wendroff schemes for elastic-plastic solids

Two formulations of the Lax-Wendroff scheme are proposed in this paper for its extension to elastic-plastic solids. The unknown vector consists of the strain or stress components in addition to the velocity ones according to the chosen formulation. The Lax-Wendroff scheme is here implemented in its Richtmyer two-step version so that the projection of the elastic trial stresses onto the yield locus is performed twice per time step. It is shown that it allows to reduce the number of cells required for the approximation of discontinuous plastic waves with respect to a classical one-step strategy consisting in an a posteriori projection of the elastic trial stresses onto the yield locus. Elastic-plastic associated and non-associated constitutive models under small strains are considered in examples. In particular, dynamic ratchetting is simulated with the Armstrong-Frederick nonlinear kinematic hardening. Comparisons with respect to finite element numerical solutions show good agreements.


Introduction
The simulation of impacts on solids finds many applications in different domains such as crash or high speed forming processes. They are modeled by means of hyperbolic initial boundary value problems, whose solutions consist of possibly discontinuous waves that are important to be tracked in order to properly understand the physical phenomena occurring in impacted parts. Indeed for solid-type media, waves can propagate both continuous and discontinuous solutions, as well as large strains and irreversible phenomena. Further, wave paths provide information about the loading history undergone by each material particle in time, which in turn allows to compute the residual strains and stresses once the steady state is achieved.
The numerical simulation of impacts on dissipative solids has been and is again mainly performed with the classical finite element method coupled with the centered differences or Newmark finite difference schemes in time [1,2], which is implemented in many industrial codes. Indeed, the finite element method enables an easy management of nonlinear partial differential equations, especially for solid-type media for which history-dependent constitutive equations can be implemented by means of appropriate integration algorithms [3], and internal variables and 'stress' ones respectively. In the 'strain' formulation, the integration of the elastic-plastic constitutive equations is performed twice per time step, first at nodes, then cellwise. The 'stress' one amounts, for its first step, to solve a small system of equations at nodes to update the stresses, provided some effective elastic-plastic compliance moduli. Then, the stress update is performed cellwise during the Richtmyer second step and still accounts for the yield threshold. At last, the plastic strains and internal variables are updated iteratively by adjusting the elastic convex so that the updated stress state lies on its boundary, consistently with the plastic flow rule and hardening laws. The 'stress' formulation presented here is actually an extension to two space dimensional problems of that proposed in [17]. It should be noted that another extension of the Lax-Wendroff scheme to hypoelastic-plastic solids has already been proposed in [20]. However, only one projection onto the yield locus is performed per time step, and the formulations and strategies presented in this work are completely different.
The paper is organized as follows. First, the elastic-plastic initial boundary value problem is presented in Section 2. Next, the Lax-Wendroff scheme derived in its Richtmyer two-step version is presented in Section 3. The proposed 'strain' and 'stress' formulations are derived in this section. Section 4 is devoted to numerical applications. Several examples conducted in two-dimensional plane strain conditions with associated and non-associated elastic-plastic constitutive models are considered in the linearized geometrical framework, and comparisons with analytical and finite element numerical solutions are presented. In particular, an example of dynamic ratchetting is simulated with the Armstrong-Frederick [29,30] nonlinear kinematic hardening to show the possibilities of the proposed approach.

Conservation laws
We consider a three-dimensional continuum body Ω ∈ R 3 , of boundary ∂Ω, and current coordinates x = x i e i . Considering the linearized geometrical framework and isothermal setting, the conservation of linear momentum and the geometrical compatibility should be satisfied for the time interval of interest t ∈ τ =]0, T ]: which hold ∀(x, t) ∈ Ω×]0, T ], and where v, σ and ε refer to the velocity vector, the Cauchy stress and the linearized strain tensors respectively. Equation (2) can also be written in conservative form:ε where 1 is the second order identity tensor, and the operator is defined in the cartesian coordinate system as 1 v = δ ik v j e i ⊗ e j ⊗ e k , δ ik being the Kronecker delta symbol. The system formed by Equations (1) and (3) appears as a system of conservation laws, written as: where U and F denote the arrays of conserved quantities and associated fluxes respectively, defined as follows: More precisely, System (4) can be derived in full matrix form with cartesian coordinates as: with vector components such that where e i (1 ≤ i ≤ 3) refer to the cartesian basis vectors, and the redundant equations due to the symmetry of the strain tensor are not considered. The fluxes F depend on both the velocity and stress components, which also defines the following auxiliary vector: The closure of the above system of conservation laws (4) is performed by means of a set of constitutive equations, introduced thereafter. Besides, Systems (4) or (6) must also be supplemented with appropriate initial and boundary conditions.

Elastic-plastic constitutive model
The thermodynamical state of the material is described by a set of state variables [31,32] that consists of the linearized strain tensor ε and the temperature T , plus some internal state variables Z describing the evolution of the internal microstructure and stored energy due to plastic deformations and other irreversible processes. Under the small strains 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 some additional variables α I , 1 ≤ I ≤ N. From a Helmholtz free-energy density potential ψ(ε, T, Z) ≡ ψ(ε − ε p , T, α I, 1≤I≤N ), concave with respect to the temperature and convex with respect to the other variables, the state laws are derived as: where s and A I refer to the entropy density and the force conjugated to the variable α I respectively. Thus, the mechanical dissipation reads: where the vector Y denotes the thermodynamic forces conjugated to the internal variables Z, and the dot (˙ ) applied to the quantity stands for a time rate. A convex elastic domain is defined by means of the yield function f , and plastic flow occurs if the current state lies and remains on its boundary, defined by f = 0. Then, provided a dissipation pseudo-potential Φ(Y), convex with respect to the flux variables Y, positive (Φ(Y) ≥ 0), and containing the origin (Φ(0) = 0), the flow rules read:Ż whereλ denotes the plastic multiplier satisfying the complementarity conditionsλ ≥ 0, f (Y) ≤ 0, λ f = 0, the identity Φλ ≥ 0, and is determined by means of the consistency conditionḟ = 0. The convexity of the flow function Φ(Y) ensures that the mechanical dissipation is non-negative. The definition of both the flow and yield function Φ and f leads to the nonassociative framework for plasticity, which may be useful for instance for the modeling of nonlinear kinematic hardening. The associative case appears as a particular case of the former one, identifying both flow and yield functions.

Mesh notations
The computational domain is divided into N c non-overlapping polygonal cells. The nodes p = 1, . . . , N n are defined at cell vertices. Following the p − c notations introduced in [11,20], the set of nodes connected to a cell c is denoted by P(c), that of cells connected to a node p by C(p), that of edges linked to a cell c by E(c), and that of edges linked to a node p by E (p). The outward unit normal to a cell c at an edge e ∈ E(c) is denoted by n (c) e . The dual mesh consists of non-overlapping polygons surrounding a node p, so that the vertices of that polygon consist of the centroids of the cells C(p) connected to a node p, and the midpoints of the primal edges e ∈ E (p) connected to a node p; see Figure 1(a). Thus, the edges of dual cells consist of separator segments connecting the centroids of primal cells to the midpoints of the primal edges E (p). The overlap between the primal and dual cells c and p forms the subcell pc, which is a quadrangle. It is defined by the separators spc− and spc+, of lengths l spc− , l spc+ and outward unit normal n spc− , n spc+ respectively, and by the half of the edges e ∈ {p − p, pp + }. The length of these edges are denoted l e , and their outward unit normals are n pc− , n pc+ respectively. The volumes (actually areas in two-dimensions) of primal and dual cells are hence related to these of subcells by the following relations: Notice that a boundary node may be connected to two or three edges, so that its dual cell consists of one or two subcells respectively (see Figure 1(b)).

Approximation
Consider a polygon of volume V, one defines the following approximation within this domain: where U V denotes the unknown array, that is the integral average of the array of conserved variables within the polygon V, which can stand for the primal cell c or the dual one p thereafter. Let's consider the time increment [t n , t n+1 ], the time step being defined as ∆t = t n+1 − t n , for which the solution vector U n V is known at time t n . Then the set of conservation laws (4) written in integral form on the polygon V can be discretized as follows [20]: where l e is the length of an edge e ∈ E(c), and the time shift index α is either equal to 1/2 for the first step or to 1 for the second step. The Lax-Wendroff scheme defined in its Richtmyer two-step version consists in solving the discrete conservation laws (18) first on the dual mesh at half-time step (α = 1/2), then on the primal mesh at the end of the time step (α = 1).
3.3. Elastic-plastic 'strain' formulation 3.3.1. Lax-Wendroff first step Particularizing Equation (18) to a dual cell V p and setting α = 1/2, one gets the updated solution array at time t n + ∆t/2: where the fluxes F n c are evaluated such that F n c = F W n c , W n c refering to the integral average of the array W (8) on the primal cell c at time t n . Notice that for boundary nodes, boundary edges exist in E (p), so that their lengths and normals l e , n pc+ appear in place of separator contributions in Equation (19). The fluxes defined on these boundary edges are computed according to prescribed boundary conditions. Ghost cells are used and combined with an elastic Riemann solver to prescribe boundary conditions, and boundary fluxes are computed as the Godunov flux of the Riemann solution.
The unknown vector U p . First, the velocity of a node p is used to compute its displacement at time t n+1 : Second, the updated strain components are used for updating the stresses and internal variables at half-time step. To do so, the strain components, the thermodynamic forces as well as the set of internal variables are approximated in the dual cell V p at time t n by performing a weighted average of their respective values in the primal cells C(p) connected to a node p: Then, the set of elastic-plastic constitutive equations can be integrated according to any available integration algorithm [3], so that the update process can be formally written as: where T denotes the numerical integrator. In this work, the cutting-plane algorithm introduced in [33, 34, 3] is used. It consists, when a plastic loading occurs ( f (Y trial n+1/2 ) > 0), in linearizing the yield function f (k) n+1/2 at iteration k about the current values Y (k) n+1/2 . Dropping the index n + 1/2 for the moment, this reads: . Then, considering a freezed strain increment at a node p (given ∆ε p ), a plastic relaxation flow is computed consistently with the partition (9), the state laws (10) and the flow rules (14). The linearization and composition of the latters yield: Introducing (24)-(25) into (23) yields, accounting for the linearity of (24)-(25) with respect to (−δλ (k) ), the increment of the plastic multiplier between two iterations: Substitution of (26) into (24)- (25) gives the updated stress σ (k+1) and the thermodynamic forces , while its introduction in the flow rules (14) yields the updated plastic strain (ε p ) (k+1) and the internal variables α (k+1) I , 1 ≤ I ≤ N. Initial conditions for the above return mapping procedure are provided by the elastic predictor. The iterative process stops when the value of the yield function is lower than a given tolerance.

Lax-Wendroff second step
Let's now write Equation (18) on a primal cell c, and set α = 1. The updated solution at time t n+1 is given by the following conservative formula: where l e denotes the length of an edge e ∈ E(c). The fluxes F n+1/2 e defined at these edges are computed by averaging those defined at the end nodes p and p + of an edge e, such that: The fluxes F n+1/2 p and F n+1/2 p + are computed from nodal solution arrays defined at time t n + ∆t/2: where W n+ 1 2 p denotes the integral average of the array W (8) on the dual cell V p at half-time step. Next, the set of elastic-plastic constitutive equations is integrated between time t n and time t n+1 in a cell c according to the cutting-plane algorithm [33,34,3] (23)- (26). Therefore, the update process can be formally written as: Hence, Equations (19) and (27) define the Richtmyer two-step algorithm for the elastic-plastic 'strain' formulation, each being supplemented by a projection step of the elastic trial stresses onto the yield locus, (22) and (30) respectively, to compute the elastic-plastic behaviour.
3.4. Elastic-plastic 'stress' formulation 3.4.1. Quasi-linear form From the system of conservation laws (4), a quasi-linear form is conveniently written as follows [35]: where the auxiliary array W , the matrix A and the k th matrix B k are defined as: where 1 q denotes the identity tensor of order q, F k = F · e k and ∂ε/∂σ ≡ S stands for the elastic-plastic tangent compliance moduli if plastic flow occurs, or simply the elastic compliance tensor for an elastic evolution.

Lax-Wendroff first step
Since the first equation of (31) is identical to that of (4), integrating it over the discrete domain V p × t n , t n+ 1 2 yields the same update of the velocity field at a node p at time t n+1/2 as the one given by Equation (19): where l spc n spc = l spc+ n spc+ + l spc− n spc− . However, the stresses now occur in the second equation of (31), thus integrating it over the discrete domain V p × t n , t n+ 1 2 yields: Considering that the compliance moduli are homogeneous in each cell, and using the divergence theorem for the second term, it comes: where have been defined Provided the linearized geometrical framework, one gets: Equation (37) is the two space dimension analog of Equation (38) derived for one-dimensional problems in [17], which has been shown to be identical to the method of characteristics.
If an elastic constitutive response is assumed for the cells C(p), an elastic trial solution can be first computed for the dual cell V p at time t n+1/2 . Assuming a homogeneous medium, the elastic compliance moduli can be factorized in the computation of the left-hand-side of (37), so that the elastic trial stresses read: where C is the elastic stiffness tensor. Notice that Equation (38) can be deduced from the second equation of (19), associated with the 'strain' formulation, by means of multiplication by the elastic stiffness tensor. The yield criterion (13) is then computed with the elastic trial stresses (38) in each cell c ∈ C(p). If the yield criterion is satisfied in all cells connected to a node p, then the elastic trial stress state is the actual one, the evolution of the dual cell V p is elastic on the time increment t n , t n+ 1 2 . If the yield criterion is not satisfied in at least one cell of this patch, then an elasticplastic correction is required.
The cells C(p) may have three types of evolution during the time increment t n , t n+ 1 2 . A first subset C E (p) ⊂ C(p) consists of cells having an elastic evolution, for which the elastic trial stresses (38) satisfy the yield criterion. A second subset C EP (p) ⊂ C(p) consists of cells undergoing a transition from elasticity at time t n ((ε p ) n c = 0) to elastoplasticity at time t n+1/2 , for which the yield criterion computed with the elastic trial state is not satisfied at time t n+1/2 while these cells had an elastic evolution at time t n . A third subset C P (p) ⊂ C(p) consists of cells having an elastic-plastic evolution, for which the yield criterion computed with the elastic trial state is not satisfied at time t n+1/2 , while plastic flow occured at time t n . The computation of the left-hand-side of Equation (37) is thus performed accordingly, splitting the set of connected cells C(p) into the three subsets C E (p), C EP (p) and C P (p). Further, the integral pertaining to the cells C EP (p) should be split into two parts in order to account for the yield threshold, so that the stress state σ * c leading to the vanishing of the yield function f (σ * c , (A I,1≤I≤N ) n c ) = 0 appears in bounds of integrals. Hence, denoting by S(σ) the elastic-plastic tangent compliance moduli, and by S the elastic compliance moduli, the particularization of Equation (37) to the three above subsets yields: where the integrals pertaining to plastic flow have been approximated assuming constant elasticplastic tangent compliance moduli over the load increment: This approximation is valid if both the plastic strain increment is limited and the normal to the yield function varies little on the time increment t n , t n+ 1 2 . Since the time step ∆t is restricted by the Courant condition, this approximation is reasonable in pratice.
The left-hand-side of Equation (39) appears as some effective compliance moduli S p of the dual cell V p on the time increment, computed by a mixture rule. However, these effective compliance moduli cannot be explicitly inverted. Indeed, the elastic-plastic tangent compliance moduli S(σ) depend explicitly on the unit normal to the yield surface N(σ), whose eigenbasis depends on the plastic flow occurring in each cell connected to a node p, that may be different. Consequently, the computation of the stresses σ n+ 1 2 p has to be performed by solving a small linear system of equations at each node of the mesh, of the form: where B p stands for the right-hand-side of Equation (39).
3.6. Lax-Wendroff second step 3.6.1. Stress update The integration of the first equation of (31) over the discrete domain V c × [t n , t n + 1] gives the update of the velocity field in a cell c at time t n+1 , which can be shown to be identical to the first equation of (27). Next, the integration of the second equation of (31) over the same discrete domain V c × [t n , t n + 1] yields an equation similar to (34), but with different bounds of integration. Its discretization reads: where have been defined and v n+ 1 2 e is defined at an edge e according to (28), by averaging the velocities defined at the nodes p and p+ at time t n+1/2 : and C is the elastic stiffness tensor. Equation (45) can also be deduced from the second equation of formula (27) linked to the Lax-Wendroff second step of the 'strain' formulation, multiplying it by the elastic stiffness tensor. The trial elastic stresses (45) are then tested against the yield criterion (13) in this cell c. If the yield criterion is satisfied, then the elastic trial stresses are the actual ones, the evolution of that cell c is elastic during this step. Conversely, if the yield criterion is not satisfied, an elastic-plastic correction must be performed. The integral appearing in the left-hand-side of Equation (42) is thus computed according to the evolution of a cell c on the time increment [t n , t n+1 ]. If it undergoes a transition from elasticity at time t n to elastoplasticity at time t n+1 , the integral is split into two parts in order to account for the yield threshold, so that the stress state σ * c corresponding to the yield surface f (σ * c , (A I,1≤I≤N ) n c ) = 0 appears as an upper and a lower bound of the first and second integrals respectively. Provided that a similar approximation to (40) is done: Equation (42), premultiplied by the elastic-plastic tangent stiffness tensor C(σ), reads: The above expression can be further simplified in the case of isotropic elasticity and J 2 flow theory. The elastic-plastic tangent stiffness tensor C(σ) can be shown to be of the form [3, Equation (2.3.9)]: where β is a function of the elastic and hardening moduli, as well as hardening variables for certain constitutive models, and N is the unit normal to the yield surface. Introducing (48) into (47) gives the stress update for a cell undergoing an elastic-elastoplastic transition: where µ is the second Lamé parameter. Conversely, if the stress state in a cell c already lied on the yield surface at time t n , then this cell undergoes an elastic-plastic evolution on the time increment [t n , t n+1 ]. Provided an approximation similar to (46), Equation (42) gives:

Update of internal variables
The update of the plastic strains and additional internal variables appears for the 'stress' formulation as a stress-driven process, while the integration algorithms used for elastic-plastic constitutive equations are generally strain-driven [3], provided an increment of strain. The procedure thus amounts to adapt the size of the elastic convex and the location of its center (if both isotropic and kinematic hardenings are accounted for) so that the updated stress state ((49) or (50)) lies on its boundary, consistently with the plastic flow rule and hardening laws (14). Since elastic-plastic equations are nonlinear, the solution process is performed iteratively. One general way to proceed consists in modifying the cutting-plane algorithm [33] so that a stress-driven update of internal variables be performed. Thus the yield function f (k) is linearized at iteration k about current values of the sole additional thermodynamic forces A I 1 ≤ I ≤ N, the stresses being given: where δA (k) Then introducing Equations (25) into (51) and following the same procedure than that of Section 3.3.1 yields the increment of the plastic multiplier between two iterations: The iterative process stops when the value of the yield function is lower than a given tolerance. Notice that for the two-dimensional plane strain case, the out-of-plane stress component σ 33 has to be updated together with the internal variables. Therefore, the yield function (51) is also linearized with respect to that component. The increment of the plastic multiplier thus reads as (26) with the sole out-of-plane component of (24).

Plane waves with Riemann-type initial conditions
where Y H = (λ + 2µ)σ y /2µ denotes the Hugoniot elastic limit, c L is the speed of the elastic pressure waves c L = (λ + 2µ)/ρ, and λ, µ and ρ denote the Lamé parameters and the mass density respectively. Provided linear isotropic or kinematic hardenings, an analytical solution of that problem has been developed in [17]. This test case allows to generate a compressive plastic reloading after an initial tensile one, by means of wave interactions and reflexions at free ends. An elastic-plastic material with linear isotropic hardening is considered, associated with the von Mises yield criterion: where s is the deviatoric part of the Cauchy stress σ, C is the hardening modulus and p stands for the cumulated plastic strain. Table 1 summarizes the material properties and geometrical parameters used for the simulations. Comparison is here performed between (i) the Lax-Wendroff

Material parameters
Geometrical parameters solutions obtained with the 'stress' and 'strain' formulations, (ii) a 'projected' solution that consists of a Lax-Wendroff elastic computation (the first step is computed with (38), the second one with (45)) followed by a projection onto the yield locus by means of a radial return [14,3], (iii) a finite element solution obtained with the finite element code Cast3M [36] Figures 2 and 3 show for two different times the superposed plots of the longitudinal stress and plastic strain components computed with these different approaches. The solution first consists of two discontinuous elastic predictors travelling leftward and rightward (see Figure 2), captured in few cells by all the methods due to the chosen Courant number. The subsequent discontinuous plastic waves need much more cells to be captured. The numerical solutions obtained with the finite elements and with the Lax-Wendroff 'stress' and 'strain' formulations give results close to each other for the longitudinal stress. They allow to slightly improve the capture of the discontinuous plastic wave enabled by an a posteriori projection of a Lax-Wendroff elastic solution onto the yield locus, although the finite element solution shows oscillations on the stress 14 plateau. However, all the methods show a non-physical peak of plastic strain that appears at the middle of the medium. This non-physical peak can be explained by the combination of the initial discontinuous velocity profile and the fact that neither Lax-Wendroff nor FEM are monotone or Total Variation Diminishing [4] (TVD) schemes. Indeed, the use of a TVD scheme in [17] on that example allowed to remove that peak. The finite element and Lax-Wendroff 'strain' solutions slightly overestimate the plastic strain plateau, while the Lax-Wendroff 'stress' one underestimates it. The 'projected' numerical solution seems closer to that plateau on this example, although its plastic waves are less sharp than these obtained by the other methods. After reflexion of waves at the two ends, their interaction yields a compressive plastic reloading as shown in Figure 3, in which the same trend as depicted above is observed. The Lax-Wendroff 'stress' formulation gives slightly less accurate plastic strains on the compressive plastic plateau than the other methods. Results obtained from computations performed on coarser (100×3 grid cells) and finer meshes (500×3 grid cells) are respectively shown in Figures 4(a) and 4(b). Coarse and underresolved grids do not show any significant difference between a double projection onto the yield locus and a single one. However as the mesh is refined, a clear and net difference between these two procedures is observed in Figure 4(b): fewer cells are required for the Lax-Wendroff 'strain' and 'stress' solutions to capture the plastic discontinuity.

Partial impact on a plane
Let's now consider the bidimensional square domain shown in Figure 5, submitted to an impact on a part of its top face, by means of a step function of a pressure p. A symmetry condition is considered on the left side, and free boundaries are set at the bottom, right, and remaining part of top faces. This problem is treated in the bidimensional plane strain framework. For this test case, a comparison between the Lax-Wendroff numerical solutions and a Q1-finite element one is considered. The finite element numerical solution is obtained with the code Cast3M [36], computed with an explicit central difference time integrator, and uses a lumped mass matrix. Computations are performed with a mesh that consists of 100 × 100 regular quadrangular grid cells, and a Courant number set at 0.9. An elastic-plastic constitutive model with a linear isotropic hardening is also considered for this example with parameters listed in Table 1, to which are added these L L −pe 2 e 1 e 2 a Figure 5: Partial impact on a plane associated with the geometry and the loading of the partial impact test case, summarized in Table  2.
Geometry Loading a = 0.5m p = 2Y H L = 2m  Figures 6 and 7 show the comparison between the finite element and the Lax-Wendroff numerical solutions at two different times. These figures consist on the one hand of plots of numerical isovalues, in particular the normal stress σ 22 and the cumulated plastic strain ε eq ≡ p, and on the other hand of superposed plots of these fields along the symmetry line (left side). At time t = 1.88 × 10 −4 seconds (see Figure 6), the pressure and shear waves generated by the step pressure prescribed on a part of the top face propagate downward. An elastic precursor is followed by a discontinuous plastic wave, which propagate plastic strains. As can be seen in Figures 6(a) and 6(c), the Lax-Wendroff 'strain' and 'stress' formulations give numerical results close to these provided by finite elements, though the latter shows much spurious numerical oscillations. In particular, the three former solutions allow to capture well the discontinuous plastic wave, better than the Lax-Wendroff projected solution. This is also seen on the plastic strain front (see Figure 6(c)), though the overall isovalues shown in Figure 6(b) are comparable for all numerical solutions.
Once the initial waves have done several round trips (see Figure 7), the plastic strains do not evolve anymore in the structure. These have been much spread in the height of the domain, and the four numerical solutions look very close to each other, and should yield almost the same residual state.

Dynamic ratchetting
Let's consider an infinite bidimensional medium, containing a set of holes aligned vertically and equally spaced, as shown in Figure 8. This infinite medium is submitted to given oscillating pressure pulses propagating horizontally. This problem is reduced according to the symmetries of the system so that the computational domain shown in Figure 9 with simplified boundary conditions is considered for the analysis. An elastic-plastic material with a nonlinear kinematic hardening is considered, with the von Mises yield criterion: where X denotes the variable defining the center of the elastic convex, satisfying the Armstrong-Frederick law [29,30] where γ denotes an additional hardening parameter. The computational domain is submitted on its right side to non-symmetric alternating tension/compression tractions: where σ a , σ 0 and T refer respectively to the stress magnitude, the stress offset and a period defined as T = L/c L corresponding to the time required for an elastic pressure wave to travel the length of the domain. Provided these non-symmetric alternating stress pulses, the nonlinear kinematic hardening law is expected to produce plastic strains growing progressively with constant pitch, especially close to the hole due to the geometry of the domain considered. This phenomenon is known as ratchetting [32], and since it is triggered here by the propagation of waves, we call it dynamic ratchetting. The present test case does not involve any discontinuous solution, but rather aims at showing the ability of Lax-Wendroff schemes to simulate unconventional dynamic plastic flows with 'refined' elastic-plastic constitutive models.
Computations are performed in the bidimensional plane strain framework, and comparison is made between the Lax-Wendroff numerical solutions and an explicit Q1-finite element solution obtained with the code Cast3M [36] as for the previous example. Computations are performed with a mesh that consists of 4200 non-rectangular quadrangular grid cells, and a Courant number set at 0.55. The parameters used for the geometry and the loading are summarized in Table 3, while these pertaining to the material are the same as these listed in Table 1, to which should be added the value of the additional hardening parameter γ = 40 of the Armstrong-Frederick law.  Figure 10 shows some comparisons between the finite element and the Lax-Wendroff numerical solutions obtained with the 'strain' and 'stress' formulations at time t = 1.25 × 10 −3 seconds, that is after two cycles of loading/unloading close to the hole area. It can be observed that the isovalues of the longitudinal stress look very close for all numerical solutions, and predict it reaches a peak value at that time close to the top point of the hole. The computed cumulated plastic strain looks also quite close in the three numerical simulations, but shows however small differences close to the top point of the hole. This is better observed in Figure 11 where, at that point, plots of the time evolution of the plastic strain component ε p 11 , and of the stress component σ 11 as a function of this plastic strain component, are made. On the one hand it can be observed that the simulations conducted with the above modeling actually predict the occurence of a ratchet phenomenon as expected. On the other hand, slight differences in the computation of the ratchet appear between the finite element and the Lax-Wendroff numerical solutions. First, the ratchet pitch computed at that point by the finite elements appears slightly greater than these obtained by the two formulations of Lax-Wendroff, which also give results close to each other. Second, stresses computed with Lax-Wendroff are slightly larger than these obtained by finite elements.

Geometry Loading
However, if we look a little bit away from that corner point as shown in Figure 12, the agreement between numerical solutions looks better, the ratchet pitches computed by the finite 19 elements and the Lax-Wendroff 'stress' formulation perfectly match, while that computed by the 'strain' formulation is slightly lower. These differences may be attributed to the different ways of prescribing boundary conditions between the finite element and Lax-Wendroff schemes. However, the results obtained show the ability of the Lax-Wendroff scheme to compute complex plastic flows in dynamic elastic-plastic numerical simulations.

Conclusion
Two formulations were proposed in this work for the extension of the well-known Lax-Wendroff scheme to elastic-plastic solids in the linearized geometrical framework. These formulations are referred to as 'strain' and 'stress' ones since the strain or the stress components are respectively included in the unknown vector in addition to the velocity ones. The Lax-Wendroff scheme was here implemented in its Richtmyer two-step version so that the projection of the  elastic trial stresses onto the yield locus is performed twice per time step. It allows to reduce the number of cells required to properly capture discontinuous plastic waves with respect to a classical one-step strategy consisting in an a posteriori projection of the elastic trial stresses onto the yield locus. Besides, it has been shown on bidimensional plane strain examples the ability of these formulations to compute complex plastic flows like dynamic ratchetting in dynamic elastic-plastic numerical simulations.