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

# Lax-Wendroff schemes for elastic-plastic solids

Thomas Heuzé

► To cite this version:

Thomas Heuzé. Lax-Wendroff schemes for elastic-plastic solids. Journal of Computational Physics, 2019, 396, pp.89-105. 10.1016/j.jcp.2019.06.050. hal-02279637

HAL Id: hal-02279637

https://hal.science/hal-02279637

Submitted on 30 Sep 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Lax-Wendroff schemes for 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

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.

**Keywords:** Elastic-plastic solids, Lax-Wendroff, Associated and non-associated plasticity

## 1. 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

Email address: thomas.heuze@ec-nantes.fr (Thomas Heuzé)
are stored at integration points in each element. However, it is first well-known that explicit finite element schemes can exhibit numerical noise near sharp solutions that can be hard to remove with artificial viscosity without loss of accuracy. Second, and especially for elastic-plastic simulations, the stresses and plastic strains converge at a rate that is one order lower than that of the displacement components.

On another side, the finite volume methods initially developed for the simulation of gas dynamics [4, 5], have gained recently more and more interest for problems involving impacts on solid media (see e.g. [6, 7, 8, 9, 10, 11, 12, 13]). This family of methods is constructed based on the characteristic structure of hyperbolic equations, either introduced explicitly by means of the solution of a Riemann problem at cell interfaces for upwind cell-centered schemes, or implicitly for centered schemes through an average process to define the updated quantities. One interesting feature of these schemes is that the same order of convergence is achieved for all components of the unknown vector, and in particular for both the velocity and stress fields in solid mechanics [9]. Next, some schemes permit to control locally the amount of numerical viscosity introduced as a function of the local regularity of the solution, so as to deliver non-oscillatory numerical solutions while preserving a high order of accuracy in more regular zones [4, 5].

Since the early work of Wilkins [14] and Trangenstein et al. [15], several authors have proposed many ways to simulate impacts on dissipative solid media with this class of methods, both for rate-independent and rate-dependent solid behaviours. Eulerian [16, 7, 6, 10, 13] as well as Lagrangian [11, 12, 17] approaches have been proposed, each pulling its own profits according to the level of strain reached, mesh entanglement and refined history-dependent constitutive equations. One key point that should be emphasized is that for rate-independent plasticity, no elastic-plastic Riemann solvers exist for two or three-dimensional stress states, even approximate. Only that of Lin & Ballman [18] was introduced for combined tension-torsion loadings, and is based on the analytical solution provided by Clifton [19]. Anyway these solvers are generally considered too complex or costly so that upwind schemes are usually coupled to linear elastic/acoustic Riemann solvers to define an elastic trial solution, and then either it is projected onto the yield locus for rate-independent constitutive behaviour [11, 20] or a relaxation step associated with a splitting scheme is performed for rate-dependent solids. The latter is often called the Maxwell-type relaxation approach [21, 7, 6, 10, 13, 8] which actually refers to an adapted version of Perzyna's elastic-viscoplastic solids [22], and converge to the rate-independent solution when the relaxation parameter tends to zero, provided the solution scheme employed to solve the source term is asymptotic preserving and accurate [23, 24]. However, both approaches generally yield a poor solution of discontinuous plastic waves (e.g. plane wave with linear hardening), which are smeared over a lot of cells even though the Courant number is set at one since the elastic precursor travels faster. Conversely, they behave better in more general cases for which plastic waves are simple waves (for nonlinear hardenings or non plane wave propagation).

One way to go around the former approaches while still avoiding elastic-plastic Riemann solvers can be done by looking at the well-known Lax-Wendroff scheme [25], that has two interesting features for the computation of elastic-plastic solids. First, it belongs to the class of centered schemes [26] so that no Riemann solver is required. Second, its Richtmyer two-step implementation [27, 28] allows to perform the projection of the elastic trial stresses onto the yield locus twice per time step, which can improve the solution of plastic waves in case these are discontinuous. The purpose of this work is to develop this idea. Therefore, two formulations of the Lax-Wendroff scheme are proposed in this paper for its extension to the computation of elastic-plastic solids. The unknown vector consists either of the strain or stress components in addition to the velocity ones according to the chosen formulation, hereafter denoted as 'STRAIN'
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.

## 2. Elastic-plastic initial boundary value problem

### 2.1. Conservation laws

We consider a three-dimensional continuum body $\Omega \in \mathbb{R}^3$, of boundary $\partial\Omega$, and current coordinates $\mathbf{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 \in \tau = ]0, T]$:

$$ \rho \dot{\mathbf{v}} - \nabla \cdot \mathbf{\sigma} = \mathbf{0} \qquad (1) $$

$$ \dot{\mathbf{\epsilon}} - \frac{1}{2}(\nabla \mathbf{v} + (\nabla \mathbf{v})^T) = \mathbf{0} \qquad (2) $$

which hold $\forall (\mathbf{x}, t) \in \Omega \times ]0, T]$, and where $\mathbf{v}$, $\mathbf{\sigma}$ and $\mathbf{\epsilon}$ refer to the velocity vector, the Cauchy stress and the linearized strain tensors respectively. Equation (2) can also be written in conservative form:

$$ \dot{\mathbf{\epsilon}} - \nabla \cdot \left( \frac{\mathbf{v} \otimes \mathbf{1} + \mathbf{1} \otimes \mathbf{v}}{2} \right) = \mathbf{0} \qquad (3) $$

where $\mathbf{1}$ is the second order identity tensor, and the operator $\otimes$ is defined in the cartesian coordinate system as $\mathbf{1} \otimes \mathbf{v} = \delta_{ik} v_j e_i \otimes e_j \otimes e_k$, $\delta_{ik}$ being the Kronecker delta symbol. The system formed by Equations (1) and (3) appears as a system of conservation laws, written as:

$$ \frac{\partial \mathcal{U}}{\partial t} + \nabla \cdot \mathcal{F} = 0 \qquad (4) $$
where $\mathcal{W}$ and $\mathcal{F}$ denote the arrays of conserved quantities and associated fluxes respectively, defined as follows:

$$ \mathcal{W} = \left\{ \begin{matrix} \mathbf{v} \\ \boldsymbol{\epsilon} \end{matrix} \right\}; \quad \mathcal{F} = \left\{ \begin{matrix} -\frac{\boldsymbol{\sigma}}{\rho} \\ -\left( \frac{\mathbf{v} \otimes \mathbf{1} + \mathbf{1} \otimes \mathbf{v}}{2} \right) \end{matrix} \right\} \qquad (5) $$

More precisely, System (4) can be derived in full matrix form with cartesian coordinates as:

$$ \frac{\partial \mathcal{W}}{\partial t} + \sum_{k=1}^{3} \frac{\partial \mathcal{F}_k}{\partial x_k} = 0 \qquad (6) $$

with vector components such that

$$ \mathcal{U} = \left\{ \begin{matrix} \mathbf{v} \\ \{\boldsymbol{\epsilon} \cdot \mathbf{e}_i\}_{1 \le i \le 3} \end{matrix} \right\}; \quad \mathcal{F}_k = \left\{ \begin{matrix} -\frac{\boldsymbol{\sigma} \cdot \mathbf{e}_k}{2} \\ -\left( \frac{v_p \delta_{ik} + v_i \delta_{pk}}{2} \right)^{\rho} \mathbf{e}_p \end{matrix} \right\}_{1 \le i \le 3} = \mathcal{F} \cdot \mathbf{e}_k \qquad (7) $$

where $\mathbf{e}_i$ ($1 \le i \le 3$) refer to the cartesian basis vectors, and the redundant equations due to the symmetry of the strain tensor are not considered. The fluxes $\mathcal{F}$ depend on both the velocity and stress components, which also defines the following auxiliary vector:

$$ \mathcal{W} = \left\{ \begin{matrix} \mathbf{v} \\ \boldsymbol{\sigma} \end{matrix} \right\} \qquad (8) $$

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.

## 2.2. 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 $\boldsymbol{\epsilon}$ and the temperature $T$, plus some internal state variables $\mathbf{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 $\boldsymbol{\epsilon}$ is additively decomposed into an elastic strain ($\boldsymbol{\epsilon}^e$) and a plastic strain ($\boldsymbol{\epsilon}^p$):

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

The set of internal state variables $\mathbf{Z}$ consists here of the plastic strain $\boldsymbol{\epsilon}^p$ plus some additional variables $\alpha_I$, $1 \le I \le N$. From a Helmholtz free-energy density potential $\psi(\boldsymbol{\epsilon}, T, \mathbf{Z}) \equiv \psi(\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^p, T, \alpha_I, 1 \le I \le N)$, concave with respect to the temperature and convex with respect to the other variables, the state laws are derived as:

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

where $s$ and $A_I$ refer to the entropy density and the force conjugated to the variable $\alpha_I$ respectively. Thus, the mechanical dissipation reads:

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

with

$$ \mathbf{Y} = \{\boldsymbol{\sigma}, A_I\}; \quad \mathbf{Z} = \{\boldsymbol{\epsilon}^p, -\alpha_I\} \qquad (12) $$
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

$$ \mathcal{C} = \{ \mathbf{Y} | f(\mathbf{Y}) \le 0 \} \qquad (13) $$

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:

$$ \dot{\mathbf{Z}} = \lambda \frac{\partial \Phi}{\partial \mathbf{Y}} ; \quad \lambda \ge 0 \qquad (14) $$

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 *f*̇ = 0. The convexity of the flow function Φ(**Y**) ensures that the mechanical dissipation

$$ \mathcal{D} = \lambda \mathbf{Y} \cdot \frac{\partial \Phi}{\partial \mathbf{Y}} \ge 0 \qquad (15) $$

is non-negative. The definition of both the flow and yield function Φ and *f* leads to the non-associative framework for plasticity, which may be useful for instance for the modeling of non-linear kinematic hardening. The associative case appears as a particular case of the former one, identifying both flow and yield functions.

### 3. Lax-Wendroff numerical schemes

#### 3.1. Mesh notations

The computational domain is divided into $N_c$ non-overlapping polygonal cells. The nodes $p = 1, \dots, 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 $\mathcal{P}(c)$, that of cells connected to a node $p$ by $\mathcal{C}(p)$, that of edges linked to a cell $c$ by $\mathcal{E}(c)$, and that of edges linked to a node $p$ by $\mathcal{E}(p)$. The outward unit normal to a cell $c$ at an edge $e \in \mathcal{E}(c)$ is denoted by $\mathbf{n}_e^{(c)}$.

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 $\mathcal{C}(p)$ connected to a node $p$, and the midpoints of the primal edges $e \in \mathcal{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 $\mathcal{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 $\mathbf{n}_{spc-}$, $\mathbf{n}_{spc+}$ respectively, and by the half of the edges $e \in \{p^-p, pp^+\}$. The length of these edges are denoted $l_e$, and their outward unit normals are $\mathbf{n}_{pc-}$, $\mathbf{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:

$$ V_c = \sum_{p \in \mathcal{P}(c)} V_{pc} ; \quad V_p = \sum_{c \in \mathcal{C}(p)} V_{pc} \qquad (16) $$

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)).
<figure><img src="image_3.png" /><figcaption>Figure 1: The dual control volumes for interior and boundary nodes are shown in red background, and the edges $e \in \mathcal{E}(p)$ connected to a node $p$ in red. The normals to the dual cell boundaries are in blue, while the normals to primal edges are in black. The reader is refered to the web version of this article for interpretation of colours.</figcaption></figure>

## 3.2. Approximation

Consider a polygon of volume $V$, one defines the following approximation within this domain:

$$ U_V = \frac{1}{V} \int_V \mathcal{W} dV \qquad (17) $$

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 $\Delta t = t_{n+1} - t_n$, for which the solution vector $U_V^n$ 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]:

$$ \frac{V(\mathbf{U}_V^{n+\alpha} - \mathbf{U}_V^n)}{\alpha \Delta t} = - \sum_{e \in \mathcal{E}(c)} l_e \mathbf{F}_e^{n+\alpha - \frac{1}{2}} \cdot \mathbf{n}_e^{n+\alpha - \frac{1}{2}} \qquad (18) $$

where $l_e$ is the length of an edge $e \in \mathcal{E}(c)$, and the time shift index $\alpha$ 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 ($\alpha = 1/2$), then on the primal mesh at the end of the time step ($\alpha = 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 $\alpha = 1/2$, one gets the updated solution array at time $t_n + \Delta t/2$:

$$V_p \mathbf{U}_p^{n+\frac{1}{2}} = \sum_{c \in C(p)} \left( V_{pc} \mathbf{U}_c^n - \frac{\Delta t}{2} \mathbf{F}_c^n \cdot (l_{spc-c} \mathbf{n}_{spc-c} + l_{spc+c} \mathbf{n}_{spc+c}) \right) \quad (19)$$

where the fluxes $\mathbf{F}_c^n$ are evaluated such that $\mathbf{F}_c^n = \mathcal{F}(\mathbf{W}_c^n)$, $\mathbf{W}_c^n$ referring to the integral average of the array $\mathcal{W}$ (8) on the primal cell $c$ at time $t_n$. Notice that for boundary nodes, boundary edges exist in $\mathcal{E}(p)$, so that their lengths and normals $l_e$, $\mathbf{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 $\mathbf{U}_p^{n+\frac{1}{2}}$ consists of both the velocity $\mathbf{v}_p^{n+\frac{1}{2}}$ and the strain components $\boldsymbol{\epsilon}_p^{n+\frac{1}{2}}$. First, the velocity of a node $p$ is used to compute its displacement at time $t_{n+1}$:

$$\mathbf{u}_p^{n+1} = \mathbf{u}_p^n + \Delta t \mathbf{v}_p^{n+\frac{1}{2}} \quad (20)$$

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$:

$$(\boldsymbol{\epsilon})_p^n = \frac{\sum_{c \in C(p)} V_{pc} \boldsymbol{\epsilon}_c^n}{V_p}; \quad \mathbf{Y}_p^n = \frac{\sum_{c \in C(p)} V_{pc} \mathbf{Y}_c^n}{V_p}; \quad \mathbf{Z}_p^n = \frac{\sum_{c \in C(p)} V_{pc} \mathbf{Z}_c^n}{V_p} \quad (21)$$

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:

$$\{ \mathbf{Y}_p^{n+\frac{1}{2}}, \mathbf{Z}_p^{n+\frac{1}{2}} \} = \mathcal{T}((\boldsymbol{\epsilon})_p^{n+\frac{1}{2}}, \mathbf{Y}_p^n, \mathbf{Z}_p^n) \quad (22)$$

where $\mathcal{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(\mathbf{Y}_{n+1/2}^{\text{trial}}) > 0$), in linearizing the yield function $f_{n+1/2}^{(k)}$ at iteration $k$ about the current values $\mathbf{Y}_{n+1/2}^{(k)}$. Dropping the index $n + 1/2$ for the moment, this reads:

$$f \approx f(\mathbf{Y}^{(k)}) + \frac{\partial f}{\partial \mathbf{Y}}|^{(k)} \cdot \delta \mathbf{Y}^{(k)} = 0 \quad (23)$$

where $\delta \mathbf{Y}^{(k)} = \mathbf{Y}^{(k+1)} - \mathbf{Y}^{(k)}$. Then, considering a frozen strain increment at a node $p$ (given $\Delta \boldsymbol{\epsilon}_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 latter yield:

$$\delta \boldsymbol{\sigma}^{(k)} = \rho \frac{\partial \psi}{\partial \delta \boldsymbol{\epsilon}^e} \circ \left( -\delta \lambda^{(k)} \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} \right)^{(k)} \quad (24)$$

$$\delta A_I^{(k)} = \rho \frac{\partial \psi}{\partial \delta \alpha_I} \circ \left( -\delta \lambda^{(k)} \frac{\partial \Phi}{\partial A_I} \right)^{(k)} \quad 1 \le I \le N \quad (25)$$
Introducing (24)-(25) into (23) yields, accounting for the linearity of (24)-(25) with respect to $(-\delta\lambda^{(k)})$, the increment of the plastic multiplier between two iterations:

$$ \delta\lambda^{(k)} = \frac{f^{(k)}}{\rho \left|\frac{\partial\psi}{\partial\delta\epsilon^s} \circ \frac{\partial\Phi}{\partial\sigma}\right|^{(k)} + \sum_{l=1}^{N} \rho \left|\frac{\partial\psi}{\partial\delta\alpha_l} \circ \frac{\partial\Phi}{\partial A_l}\right|^{(k)}} ; \quad \lambda^{(k+1)} = \lambda^{(k)} + \delta\lambda^{(k)} \quad (26) $$

Substitution of (26) into (24)-(25) gives the updated stress $\sigma^{(k+1)}$ and the thermodynamic forces $A_l^{(k+1)}$, while its introduction in the flow rules (14) yields the updated plastic strain $(\epsilon^p)^{(k+1)}$ and the internal variables $\alpha_l^{(k+1)}$, $1 \le I \le 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.

### 3.3.2. Lax-Wendroff second step

Let's now write Equation (18) on a primal cell $c$, and set $\alpha = 1$. The updated solution at time $t_{n+1}$ is given by the following conservative formula:

$$ \mathbf{U}_c^{n+1} = \mathbf{U}_c^n - \frac{\Delta t}{V_c} \sum_{e \in \mathcal{E}(c)} (l_e \mathbf{F}_e^{n+1/2} \cdot \mathbf{n}_e^{(c)}) \quad (27) $$

where $l_e$ denotes the length of an edge $e \in \mathcal{E}(c)$. The fluxes $\mathbf{F}_e^{n+1/2}$ defined at these edges are computed by averaging those defined at the end nodes $p$ and $p^+$ of an edge $e$, such that:

$$ \mathbf{F}_e^{n+1/2} = \frac{\mathbf{F}_p^{n+1/2} + \mathbf{F}_{p^+}^{n+1/2}}{2} \quad (28) $$

The fluxes $\mathbf{F}_p^{n+1/2}$ and $\mathbf{F}_{p^+}^{n+1/2}$ are computed from nodal solution arrays defined at time $t_n + \Delta t/2$:

$$ \mathbf{F}_p^{n+1/2} = \mathcal{F}(\mathbf{W}_p^{n+1/2}) ; \quad \mathbf{F}_{p^+}^{n+1/2} = \mathcal{F}(\mathbf{W}_{p^+}^{n+1/2}) \quad (29) $$

where $\mathbf{W}_p^{n+1/2}$ denotes the integral average of the array $\mathcal{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:

$$ \{\mathbf{Y}_c^{n+1}, \mathbf{Z}_c^{n+1}\} = \mathcal{F}(\boldsymbol{\epsilon}_c^{n+1}, \mathbf{Y}_c^n, \mathbf{Z}_c^n) \quad (30) $$

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]:

$$ \mathcal{A} \frac{\partial \mathcal{W}}{\partial t} + \sum_{k=1}^{3} \mathcal{B}_k \frac{\partial \mathcal{W}}{\partial x_k} = 0 \quad (31) $$
where the auxiliary array $\mathcal{W}$, the matrix $\mathcal{A}$ and the $k^{th}$ matrix $\mathcal{B}_k$ are defined as:

$$ \mathcal{W} = \begin{Bmatrix} \mathbf{v} \\ \boldsymbol{\sigma} \end{Bmatrix}; \quad \mathcal{A} = \frac{\partial \mathcal{W}}{\partial \mathcal{W}} = \begin{bmatrix} \mathbf{1}_2 & \mathbf{0}_3 \\ \mathbf{0}_3 & \frac{\partial \boldsymbol{\varepsilon}}{\partial \boldsymbol{\sigma}} \end{bmatrix}; \quad \mathcal{B}_k = \frac{\partial \mathcal{F}_k}{\partial \mathcal{W}} = - \begin{bmatrix} \mathbf{0}_2 & \frac{\mathbf{1}_4 \cdot \mathbf{e}_k}{\rho} \\ \mathbf{1}_4 \cdot \mathbf{e}_k & \mathbf{0}_4 \end{bmatrix} \quad (32) $$

where $\mathbf{1}_q$ denotes the identity tensor of order $q$, $\mathcal{F}_k = \mathcal{F} \cdot \mathbf{e}_k$ and $\partial\boldsymbol{\varepsilon}/\partial\boldsymbol{\sigma} \equiv \mathcal{S}$ stands for the elastic-plastic tangent compliance moduli if plastic flow occurs, or simply the elastic compliance tensor for an elastic evolution.

### 3.5. Lax-Wendroff first step

Since the first equation of (31) is identical to that of (4), integrating it over the discrete domain $V_p \times [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):

$$ V_p \mathbf{v}_p^{n+1/2} = \sum_{c \in C(p)} \left( V_{pc} \mathbf{v}_c^n + \frac{\Delta t}{2} \frac{l_{spc}}{\rho} \boldsymbol{\sigma}_c^n \cdot \mathbf{n}_{spc} \right) \quad (33) $$

where $l_{spc}\mathbf{n}_{spc} = l_{spc+}\mathbf{n}_{spc+} + l_{spc-}\mathbf{n}_{spc-}$. However, the stresses now occur in the second equation of (31), thus integrating it over the discrete domain $V_p \times [t_n, t_{n+1/2}]$ yields:

$$ \int_{t_n}^{t_{n+1/2}} \int_{V_p} \left[ \mathcal{S} : \dot{\boldsymbol{\sigma}} - \nabla \cdot \left( \frac{\mathbf{v} \otimes \mathbf{1} + \mathbf{1} \boxtimes \mathbf{v}}{2} \right) \right] dV dt = \mathbf{0} \quad (34) $$

Considering that the compliance moduli are homogeneous in each cell, and using the divergence theorem for the second term, it comes:

$$ \int_{t_n}^{t_{n+1/2}} \sum_{c \in C(p)} V_{pc} \mathcal{S}_c : \dot{\boldsymbol{\sigma}}_c dt = \frac{\Delta t}{2} \sum_{c \in C(p)} \int_{[spc-\cup spc+]} \left( \frac{\mathbf{v}_c^n \otimes \mathbf{n} + \mathbf{n} \otimes \mathbf{v}_c^n}{2} \right) dl \quad (35) $$

where have been defined

$$ \boldsymbol{\sigma}_c = \frac{1}{V_{pc}} \int_{V_{pc}} \boldsymbol{\sigma} dV; \quad \mathbf{v}^n = \frac{1}{\Delta t} \int_{t_n}^{t_{n+1/2}} \mathbf{v} dt \quad (36) $$

Provided the linearized geometrical framework, one gets:

$$ \sum_{c \in C(p)} V_{pc} \int_{\boldsymbol{\sigma}_c^n}^{\boldsymbol{\sigma}_p^{n+1/2}} \mathcal{S}_c : d\boldsymbol{\sigma} = \frac{\Delta t}{2} \sum_{c \in C(p)} l_{spc} \left( \frac{\mathbf{v}_c^n \otimes \mathbf{n}_{spc} + \mathbf{n}_{spc} \otimes \mathbf{v}_c^n}{2} \right) \quad (37) $$

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:

$$ V_p \boldsymbol{\sigma}_p^{\text{trial}} = \sum_{c \in C(p)} \left( V_{pc} \boldsymbol{\sigma}_c^n + \frac{\Delta t}{2} l_{spc} C_c : \mathbf{v}_c^n \otimes \mathbf{n}_{spc} \right) \quad (38) $$
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 \in 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+\frac{1}{2}}]$. If the yield criterion is not satisfied in at least one cell of this patch, then an elastic-plastic correction is required.

The cells $C(p)$ may have three types of evolution during the time increment $[t_n, t_{n+\frac{1}{2}}]$. A first subset $C^E(p) \subset 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) \subset C(p)$ consists of cells undergoing a transition from elasticity at time $t_n$ ($\dot{\varepsilon}^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) \subset 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 $\sigma_c^*$ leading to the vanishing of the yield function $f(\sigma_c^*, (A_{I,1} \le I \le N})^n_c) = 0$ appears in bounds of integrals. Hence, denoting by $S(\sigma)$ the elastic-plastic tangent compliance moduli, and by $S$ the elastic compliance moduli, the particularization of Equation (37) to the three above subsets yields:

$$ 
\begin{aligned}
& \left( \sum_{c \in C^E(p)} V_{pc} S_c + \sum_{c \in C^{EP}(p)} V_{pc} S_c(\sigma_c^*) + \sum_{c \in C^P(p)} V_{pc} S_c(\sigma_c^n) \right) : \sigma_p^{n+\frac{1}{2}} \\
& = \sum_{c \in C^E(p)} V_{pc} S_c : \sigma_c^n + \sum_{c \in C^{EP}(p)} V_{pc} S_c(\sigma_c^*) : \sigma_c^* - \sum_{c \in C^{EP}(p)} V_{pc} S_c : (\sigma_c^* - \sigma_c^n) \\
& + \sum_{c \in C^P(p)} V_{pc} S_c(\sigma_c^n) : \sigma_c^n + \frac{\Delta t}{2} \sum_{c \in C(p)} l_{spc} \left( \frac{\mathbf{v}_c^n \otimes \mathbf{n}_{spc} + \mathbf{n}_{spc} \otimes \mathbf{v}_c^n}{2} \right)
\end{aligned} 
\quad (39) $$

where the integrals pertaining to plastic flow have been approximated assuming constant elastic-plastic tangent compliance moduli over the load increment:

$$ \int_{\sigma_c^*}^{\sigma_p^{n+\frac{1}{2}}} \mathbf{S}_c(\sigma) : d\sigma \approx \mathbf{S}_c(\sigma_c^*) : (\sigma_p^{n+\frac{1}{2}} - \sigma_c^*) \quad (40) $$

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+\frac{1}{2}}]$. Since the time step $\Delta t$ is restricted by the Courant condition, this approximation is reasonable in practice.

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(\sigma)$ depend explicitly on the unit normal to the yield surface $N(\sigma)$, 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 $\sigma_p^{n+\frac{1}{2}}$ has to be performed by solving a small linear system of equations at each node of the mesh, of the form:

$$S_p : \sigma_p^{n+\frac{1}{2}} = \mathbf{B}_p \qquad (41)$$

where $\mathbf{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 \times [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 \times [t_n, t_n + 1]$ yields an equation similar to (34), but with different bounds of integration. Its discretization reads:

$$\int_{\sigma_c^n}^{\sigma_c^{n+1}} S_c : d\sigma = \frac{\Delta t}{V_c} \sum_{e \in \mathcal{E}(c)} l_e \left( \frac{\mathbf{v}_e^{n+\frac{1}{2}} \otimes \mathbf{n}_e^{(c)} + \mathbf{n}_e^{(c)} \otimes \mathbf{v}_e^{n+\frac{1}{2}}}{2} \right) \qquad (42)$$

where have been defined

$$\sigma_c = \frac{1}{V_c} \int_{V_c} \sigma dV \quad ; \quad \mathbf{v}^{n+\frac{1}{2}} = \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} \mathbf{v} dt \qquad (43)$$

and $\mathbf{v}_e^{n+\frac{1}{2}}$ 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}$:

$$\mathbf{v}_e^{n+\frac{1}{2}} = \frac{\mathbf{v}_p^{n+\frac{1}{2}} + \mathbf{v}_{p+}^{n+\frac{1}{2}}}{2} \qquad (44)$$

An elastic trial solution can be first computed for a cell c. In that case, Equation (42) reads:

$$\sigma_c^{\text{trial}} = \sigma_c^n + \frac{\Delta t}{V_c} C_c : \sum_{e \in \mathcal{E}(c)} l_e \mathbf{v}_e^{n+\frac{1}{2}} \otimes \mathbf{n}_e^{(c)} \qquad (45)$$

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 $\sigma_c^*$ corresponding to the yield surface $f(\sigma_c^*, (A_{I,1} \le I \le N})_c^n = 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:

$$\int_{\sigma_c^*}^{\sigma_c^{n+1}} \mathbf{S}(\sigma) : d\sigma \approx \mathbf{S}(\sigma_c^*) : (\sigma_c^{n+1} - \sigma_c^*), \qquad (46)$$
Equation (42), premultiplied by the elastic-plastic tangent stiffness tensor C(σ), reads:

$$ \sigma_c^{n+1} = C_c(\sigma_c^*) : S_c : \sigma_c^n + (1_4 - C_c(\sigma_c^*) : S_c) : \sigma_c^* \\ + \frac{\Delta t}{V_c} C_c(\sigma_c^*) : \sum_{e \in \mathcal{E}(c)} l_e \mathbf{v}_e^{n+\frac{1}{2}} \otimes \mathbf{n}_e^{(c)} \quad (47) $$

The above expression can be further simplified in the case of isotropic elasticity and J₂ flow theory. The elastic-plastic tangent stiffness tensor C(σ) can be shown to be of the form [3, Equation (2.3.9)]:

$$ C(\sigma) = C - \beta \mathbf{N} \otimes \mathbf{N} \quad (48) $$

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:

$$ \sigma_c^{n+1} = \sigma_c^n + \frac{\beta}{2\mu} \mathbf{N}^* \otimes \mathbf{N}^* : (\sigma_c^* - \sigma_c^n) + \frac{\Delta t}{V_c} C_c(\sigma_c^*) : \sum_{e \in \mathcal{E}(c)} l_e \mathbf{v}_e^{n+\frac{1}{2}} \otimes \mathbf{n}_e^{(c)} \quad (49) $$

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:

$$ \sigma_c^{n+1} = \sigma_c^n + \frac{\Delta t}{V_c} C_c(\sigma_c^n) : \sum_{e \in \mathcal{E}(c)} l_e \mathbf{v}_e^{n+\frac{1}{2}} \otimes \mathbf{n}_e^{(c)} \quad (50) $$

### 3.6.2. 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:

$$ f \approx f(\sigma^{n+1}, A_I^{(k)}) + \sum_{I=1}^{N} \left. \frac{\partial f}{\partial A_I} \right|^{(k)} \cdot \delta A_I^{(k)} = 0 \quad (51) $$

where δA_I^{(k)} = A_I^{(k+1)} - A_I^{(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:

$$ \delta\lambda^{(k)} = \frac{f^{(k)}}{\sum_{I=1}^{N} \rho \frac{\partial\psi}{\partial\delta\alpha_I} \circ \frac{\partial\Phi}{\partial A_I}^{(k)}} ; \quad \lambda^{(k+1)} = \lambda^{(k)} + \delta\lambda^{(k)} \quad (52) $$
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 $\sigma_{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).

# 4. Applications

## 4.1. Plane waves with Riemann-type initial conditions

Let's consider a rectangular domain $\Omega : (x, y) \in ]0, L[ \times ]0, H[$, made of an elastic-plastic material. Traction free conditions are set at the left and right sides, while zero shear stress and normal 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. Riemann-type initial conditions are prescribed: the longitudinal component of 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 stresses are considered to be zero everywhere initially. The prescribed velocity is set so that plastic flow occurs:

$$ \bar{v} = 2 \frac{Y_H}{\rho c_L} \qquad (53) $$

where $Y_H = (\lambda + 2\mu)\sigma_y/2\mu$ denotes the Hugoniot elastic limit, $c_L$ is the speed of the elastic pressure waves $c_L = \sqrt{(\lambda + 2\mu)/\rho}$, and $\lambda, \mu$ and $\rho$ 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:

$$ f = \sigma_{\text{eq}} - C p - \sigma_y ; \quad \sigma_{\text{eq}} = \sqrt{\frac{3}{2} \mathbf{s} : \mathbf{s}} ; \quad \mathcal{C} = \{ \sigma | f(\sigma) \le 0 \} \qquad (54) $$

where $\mathbf{s}$ is the deviatoric part of the Cauchy stress $\sigma$, $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

<table><caption>Table 1: Parameters</caption><thead><tr><th colspan="2">Material parameters</th><th>Geometrical parameters</th></tr></thead><tbody><tr><td>$E = 200\text{GPa}$</td><td>$C = 10\text{GPa}$</td><td>$L = 6\text{m}$</td></tr><tr><td>$\nu = 0.3$</td><td>$\sigma_y = 400\text{MPa}$</td><td>$H = 0.08\text{m}$</td></tr><tr><td>$\rho = 7800\text{kg.m}^{-3}$</td><td></td><td></td></tr></tbody></table>

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] and computed
with an explicit central difference time integrator, Q1 elements and a lumped mass matrix [1],
and (iv) the analytical solution. Computations are performed with a mesh made of 200 × 3
quadrangular grid cells and a Courant number set at 0.9.

<figure><img src="image_4.png" /><figcaption>Figure 2: Comparison at time 4.37 × 10⁻⁴ seconds of the Lax-Wendroff solutions ('STRESS' and 'STRAIN' formulations), an elastic Lax-Wendroff solution projected onto the yield locus, a finite element solution and the analytical stress and plastic strain fields.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Comparison at time 1.38 × 10⁻³ seconds of the Lax-Wendroff solutions ('STRESS' and 'STRAIN' formulations), an elastic Lax-Wendroff solution projected onto the yield locus, a finite element solution and the analytical stress and plastic strain fields.</figcaption></figure>

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

<figure><img src="image_6.png" /><figcaption>Figure 4: Comparison at time $4.37 \times 10^{-4}$ seconds of the Lax-Wendroff solutions ('STRESS' and 'STRAIN' formulations), an elastic Lax-Wendroff solution projected onto the yield locus and the analytical stress for coarser and finer meshes.</figcaption></figure>

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.

## 4.2. 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 \times 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
<figure><img src="image_7.png" /><figcaption>Figure 5: Partial impact on a plane</figcaption></figure>

associated with the geometry and the loading of the partial impact test case, summarized in Table 2.

<table><caption>Table 2: Parameters of the partial impact test case</caption><thead><tr><td>Geometry</td><td>Loading</td></tr></thead><tbody><tr><td>a = 0.5m</td><td rowspan="2">p = 2Y<sub>H</sub></td></tr><tr><td>L = 2m</td></tr></tbody></table>

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 $\sigma_{22}$ and the cumulated plastic strain $\varepsilon_{eq} \equiv p$, and on the other hand of superposed plots of these fields along the symmetry line (left side). At time $t = 1.88 \times 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.
<figure><img src="image_8.png" /><figcaption>Figure 6: Comparison at time $t = 1.88 \times 10^{-4}$ seconds of the Lax-Wendroff solutions ('STRESS' and 'STRAIN' formulations), an elastic Lax-Wendroff solution projected onto the yield locus and a finite element solution.</figcaption></figure>

### 4.3. 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:

$$f = \sigma_{\text{eq}}(\xi) - \sigma_y ; \quad \sigma_{\text{eq}} = \sqrt{\frac{3}{2}\xi : \xi} ; \quad \xi = s - X ; \quad \mathcal{C} = \{(\sigma, X) | f \le 0\} \quad (55)$$
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: Comparison at time $t = 5 \times 10^{-3}$ seconds of the Lax-Wendroff solutions ('STRESS' and 'STRAIN' formulations), an elastic Lax-Wendroff solution projected onto the yield locus and a finite element solution.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8: Infinite bidimensional medium, containing a set of holes, submitted to given oscillating pressure pulses propagating horizontally.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9: Computational domain and boundary conditions</figcaption></figure>

where X denotes the variable defining the center of the elastic convex, satisfying the Armstrong-Frederick law [29, 30]:

$$ \dot{\mathbf{X}} = \frac{2}{3} C \dot{\boldsymbol{\epsilon}}^p - \gamma \mathbf{X} \dot{\boldsymbol{p}} \quad (56) $$
where $\gamma$ denotes an additional hardening parameter. The computational domain is submitted on its right side to non-symmetric alternating tension/compression tractions:

$$ \sigma_{\infty}(t) = \sigma_a \sin\left(\frac{\pi t}{2T}\right) + \sigma_0 \qquad (57) $$

where $\sigma_a$, $\sigma_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 $\gamma = 40$ of the Armstrong-Frederick law.

<table><caption>Table 3: Parameters of the dynamic ratchetting test case</caption><thead><tr><td>Geometry</td><td>Loading</td></tr></thead><tbody><tr><td>r = 0.5m</td><td>σ<sub>a</sub> = 3Y<sub>H</sub></td></tr><tr><td>L = 1m</td><td>σ<sub>0</sub> = Y<sub>H</sub>/2</td></tr></tbody></table>

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 \times 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 $\epsilon_{11}^p$, and of the stress component $\sigma_{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 occurrence 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.

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
<figure><img src="image_13.png" /><figcaption>Figure 10: Comparison at time $t = 1.25 \times 10^{-3}$ seconds of the Lax-Wendroff solutions ('STRESS' and 'STRAIN' formulations) and a finite element solution.</figcaption></figure>

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.

## 5. 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
<figure><img src="image_14.png" /><figcaption>Figure 11: Plots at the top point of the hole ($x_1 = 0, x_2 = r$) of the time evolution of the plastic strain component $\epsilon_{11}^p$ (at the left), and that of the stress component $\sigma_{11}$ as a function of this plastic strain component (at the right).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 12: Plots at the location ($x_1 = r/5, x_2 = r$) of the time evolution of the plastic strain component $\epsilon_{11}^p$ (at the left), and that of the stress component $\sigma_{11}$ as a function of this plastic strain component (at the right).</figcaption></figure>

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.

[1] T. Belytschko, W. Liu, B. Moran, Nonlinear finite elements for continua and structures, Wiley, 2000 (2000).

[2] D. J. Benson, Computational methods in Lagrangian and Eulerian hydrocodes, Computer Methods in Applied Mechanics and Engineering 99 (1992) 235–394 (1992).

[3] J. Simo, T. Hughes, Computational inelasticity, Springer, 1997 (1997).

[4] R. Leveque, Finite volume methods for hyperbolic problems, Cambridge University Press, 2002 (2002).

[5] E. Toro, Riemann solvers and numerical methods for fluid dynamics, Springer, 2009 (2009).

[6] P. Barton, D. Drikakis, E. Romenskii, A high-order Eulerian Godunov method for elastic-plastic flow in solids, International Journal for Numerical Methods in Engineering 81 (2010) 453–484 (2010).

[7] D. Hill, D. Pullin, M. Ortiz, D. Meiron, An Eulerian hybrid WENO centered-difference solver for elastic-plastic solids, Journal of Computational Physics 229 (2010) 9053–9072 (2010).

[8] N. Favrie, S. Gavrilyuk, Mathematical and numerical model for nonlinear viscoplasticity, Philosophical Transactions of the Royal Society A 369 (2011) 2864–2880 (2011).

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

[10] A. L. Ortega, M. Lombardini, D. Pullin, D. Meiron, Numerical simulation of elastic-plastic solid mechanics using
an Eulerian stretch tensor approach and HLLD Riemann solver, Journal of Computational Physics 257 (2014) 414–441 (2014).

[11] P. Maire, R. Abgrall, J. Breil, R. Loubère, 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 (2013) 626–665 (2013).

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

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

[14] M. Wilkins, Methods in computational physics, Vol. 3, Academic Press, 1964, Ch. Calculation of elastic-plastic flow, pp. 211–263 (1964).

[15] J. Trangenstein, P. Collela, A higher-order Godunov method for modeling finite deformation in elastic-plastic solids, Communications in Pure Applied mathematics 47 (1991) 41–100 (1991).

[16] G. Miller, P. Collela, A high-order eulerian godunov method for elastic-plastic flow in solids, Journal of Computational Physics 167 (2001) 131–176 (2001).

[17] T. Heuzé, Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elasticplastic solids, Journal of Computational Physics 346 (2017) 369–388 (2017).

[18] X. Lin, J. Ballman, A Riemann solver and a second-order Godunov method for elastic-plastic wave propagation in solids, International Journal of Impact Engineering 13 (1993) 463–478 (1993).

[19] R. Clifton, An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube, Tech. rep., Brown university Providence R.I, Division of Engineering (1966).

[20] D. Fridrich, R. Liska, B. Wendroff, Cell-centered lagangian Lax-Wendroff HLL hybrid method for elastic-plastic flows, Computers and Fluids 157 (2017) 164–174 (2017).

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

[22] P. Perzyna, Fundamental problems in viscoplasticity, Advances in Applied Mechanics 9 (1966) 243–277 (1966).

[23] T. Heuzé, Simulation of impacts on elasticviscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes, Advanced Modeling and Simulation in Engineering Sciences 5 (1) (2018) 9 (2018).

[24] M. Dumbser, C. Enaux, E. Toro, Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws, Journal of computational physics 227 (2008) 3971–4001 (2008).

[25] P. Lax, B. Wendroff, Systems of conservation laws, Communications on Pure and Applied mathematics 13 (2) (1960) 217–237 (1960).

[26] E. Toro, S. Billett, Centred TVD schemes for hyperbolic conservation laws, IMA Journal of Numerical Analysis 20 (2000) 47–79 (2000).

[27] R. Richtmyer, A survey of difference methods for non-steady fluid dynamics, Tech. rep., National Center for Atmospheric Research, http://dx.doi.org/10.5065/D67P8WCQ (1962).

[28] G. Zwas, On Two Step Lax-Wendroff Methods in Several Dimensions, Numerische Mathematik 20 (1973) 350–355 (1973).

[29] C. Frederick, P. Armstrong, A mathematical representation of the multiaxial Bauchinger effect, Tech. rep., Central Electricity Generating Board, Berkeley, UK (1966).

[30] C. Frederick, P. Armstrong, A mathematical representation of the multiaxial Bauschinger effect, Materials at High Temperatures 24 (2007) 1–26 (2007).

[31] P. Germain, Q. Nguyen, P. Suquet, Continuum Thermodynamics, Journal of Applied Mechanics 50 (1983) 1010–1020 (1983).

[32] J. Lemaitre, J.L.Chaboche, Mechanics of Solid Materials, Cambridge University Press, 1994 (1994).

[33] J. Simo, M. Ortiz, A unified approach to finite deformation elastoplastic analysis based on the use of hyperelastic constitutive equations, Computer methods in applied mechanics and engineering 49 (2) (1985) 221–245 (1985).

[34] M. Ortiz, J. Simo, An analysis of a new class of integration algorithms for elastoplastic constitutive relations, International Journal for Numerical Methods in Engineering 23 (3) (1986) 353–366 (1986).

[35] J. Trangenstein, Numerical solution of hyperbolic partial differential equations, Cambridge, 2009 (2009).

[36] Cast3M, user's manual. www.cast3M-cea.fr (2018).