Keywords: cell-centered Lagrangian hydrodynamics, node-centered solver, Godunov-type method, elastic-plastic flow, high-order finite volume methods, multi-dimensional unstructured mesh, Generalized Riemann Problem, Geometric Conservation Hydrodynamique lagrangienne, schéma centré cellule, Méthode de Godunov, Ecoulements elasti-plastiques, schémas volumes fini d'ordre élevé, maillages non structurés, problde Rieman

In this paper, we describe a cell-centered Lagrangian scheme devoted to the numerical simulation of solid dynamics on two-dimensional unstructured grids in planar geometry. This numerical method, utilizes the classical elastic-perfectly plastic material model initially proposed by Wilkins [M.L. Wilkins, Calculation of elastic-plastic flow, Meth. Comput. Phys. (1964)]. In this model, the Cauchy stress tensor is decomposed into the sum of its deviatoric part and the thermodynamic pressure which is defined by means of an equation of state. Regarding the deviatoric stress, its time evolution is governed by a classical constitutive law for isotropic material. The plasticity model employs the von Mises yield criterion and is implemented by means of the radial return algorithm. The numerical scheme relies on a finite volume cell-centered method wherein numerical fluxes are expressed in terms of sub-cell force. The generic form of the sub-cell force is obtained by requiring the scheme to satisfy a semi-discrete dissipation inequality. Sub-cell force and nodal velocity to move the grid are computed consistently with cell volume variation by means of a node-centered solver, which results from total energy conservation. The nominally second-order extension is achieved by developing a two-dimensional extension in the Lagrangian framework of the Generalized Riemann Problem methodology, introduced by Ben-Artzi and Falcovitz [M. Ben-Artzi and J. Falcovitz, Generalized Riemann Problems in Computational Fluid Dynamics, Cambridge Monographs on Applied and Computational Mathematics, 2003]. Finally, the robustness and the accuracy of the numerical scheme are assessed through the computation of several test cases. .

Un schéma du second ordre de type lagrangien pour la simulation des écoulements élasto-plastiques sur maillaes non structurés bidimensionnels Résumé : In this paper, we describe a cell-centered Lagrangian scheme devoted to the numerical simulation of solid dynamics on two-dimensional unstructured grids in planar geometry. This numerical method, utilizes the classical elastic-perfectly plastic material model initially proposed by Wilkins [M.L. Wilkins, Calculation of elastic-plastic flow, Meth. Comput. Phys. (1964)]. In this model, the Cauchy stress tensor is decomposed into the sum of its deviatoric part and the thermodynamic pressure which is defined by means of an equation of state. Regarding the deviatoric stress, its time evolution is governed by a classical constitutive law for isotropic material.

The plasticity model employs the von Mises yield criterion and is implemented by means of the radial return algorithm. The numerical scheme relies on a finite volume cell-centered method wherein numerical fluxes are expressed in terms of sub-cell force. The generic form of the sub-cell force is obtained by requiring the scheme to satisfy a semi-discrete dissipation inequality. Sub-cell force and nodal velocity to move the grid are computed consistently with cell volume variation by means of a node-centered solver, which results from total energy conservation. The nominally second-order extension is achieved by developing a two-dimensional extension in the Lagrangian framework of the Generalized Riemann Problem methodology, introduced by Ben-Artzi and Falcovitz [M. Ben-Artzi and J. Falcovitz, Generalized Riemann Problems in Computational Fluid Dynamics, Cambridge Monographs on Applied and Computational Mathematics, 2003]. Finally, the robustness and the accuracy of the numerical scheme are assessed through the computation of several test cases.

Introduction

This paper aims at describing a high-order cell-centered Lagrangian scheme for elastic-plastic solids in two-dimensional planar geometry. We point out that this numerical method is the non trivial extension to solid dynamics of the EUCCLHYD (Explicit Unstructured Cell-Centered Lagrangian HYDrodynamics) scheme which has been initially derived for solving gas dynamics equations written under the Lagrangian formalism, refer to [START_REF] Maire | A high-order one-step sub-cell force-based discretization for cell-centered lagrangian hydrodynamics on polygonal grids[END_REF][START_REF] Maire | Contribution to the numerical modeling of inertial confinement fusion[END_REF][START_REF] Maire | A unified sub-cell force-based discretization for cell-centered lagrangian hydrodynamics on polygonal grids[END_REF][START_REF] Maire | A cell-centered Lagrangian scheme for twodimensional compressible flow problems[END_REF]. Let us recall that this formalism is characterized by the fact that the computational grid is attached to the material motion. Thus, the present numerical scheme belongs to the class of finite volume methods on moving mesh. Moreover, the scheme is referred as to cell-centered since the primarily unknowns are located at the cell center.

Here, we study the discretization of the isotropic elastic-plastic model that has been initially introduced by Wilkins [START_REF] Wilkins | Methods in Computational Physics, volume 3, chapter Calculation of Elastic-Plastic Flow[END_REF]. In this seminal paper, Wilkins has derived a staggered scheme, wherein the kinematical variables are located at the node whereas the thermodynamical variables are defined at the cell center. This numerical method solves the specific internal energy equation and requires the use of an artificial viscosity to ensure a proper dissipation of kinetic energy into internal energy through irreversible transformation such as shock waves. The elastic-plastic model consists of the usual conservation laws of the volume, the momentum and the total energy. The Cauchy stress tensor is decomposed into a spherical part, the thermodynamic pressure, and a deviatoric part, the deviatoric stress. The pressure is expressed by means of an equation of state whereas the deviatoric stress is described by means of a constitutive law written under a rate form. This constitutive law expresses the Jaumann rate of the deviatoric stress in terms of the strain rate tensor. Let us note that the Jaumann rate is a specific rate which ensures that the constitutive law satisfies the principle of material objectivity otherwise named principle of material frame-indifference, refer to [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF]. Namely, the constitutive law should not depend on whatever external frame of reference is used to describe it. The plasticity modeling relies on the classical additive decomposition of the strain rate tensor into an elastic part and a plastic deviatoric part. The plastic strain rate tensor is defined by means of a plastic flow rule which is of Prandtl-Reuss type [START_REF] De Souza Neto | Computational Methods for Plasticity: Theory and Applications[END_REF]. The plastic yielding is characterized by the von Mises yield criterion which is particularly appropriate for metals [START_REF] Lemaitre | Mécanique des matériaux solides[END_REF].

The previous elastic-plastic model relies on the use of a hypoelastic-based constitutive formulation. In this framework, standard infinitesimal elastoplasticity models are extended to the finite strain range by formulating the constitutive law in terms of frame invariant stress rates. In spite of its simplicity, this approach arises two main difficulties. First, the definition of objective stress rates is somewhat arbitrary and many different stress rates have been proposed, refer to [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF]. Moreover, the discretization of these objective stress rates must be performed with great care to ensure that the constitutive law at the discrete level still satisfies the principle of material objectivity. This important property refers to as the incremental objectivity. Let us also note that an objective stress rate such as the Jaumann rate leads to an evolution equation for the stress which cannot be written under conservative form for multi-dimensionnal flows. This flaw renders the mathematical analysis of discontinuous solutions questionable as noticed in [START_REF] Gavrilyuk | Modelling wave dynamics of compressible elastic materials[END_REF]. Second, hypoelastic-based models are characterized by a dissipative behavior within the elastic regime. In other words, reversible transformations of an elastic material are characterized by a non-zero rate of entropy. Let us point out that these difficulties do not prevent from developing numerical methods for simulating hypoelastic-based models. For instance, Trangenstein and Collela [START_REF] Trangenstein | A High-Order Godunov Method for Modeling Finite Deformation in Elastic-Plastic Solids[END_REF], have developed a first-order system of conservation laws for finite deformations in solids. They describe the characteristic structure of this system for both Lagrangian and Eulerian formalisms and provide a high-order Godunov method for problems involving finite deformation and plasticity. They also investigate the hyperbolic structure of models written in Inria rate-form which are quite close to the Wilkins model. They show that for Jaumann rate-based models the characteristic speeds are always real within the elastic regime. In [START_REF] Howell | A free-Lagrange augmented Godunov method for the simulation of elastic-plastic solids[END_REF], Howell and Ball have proposed a Free-Lagrange type numerical method for the simulation of elastic-plastic solids. Their scheme, which employs the Wilkins model, relies on a time-operator splitting procedure wherein the hydrodynamic fluxes are considered separately from the strength fluxes. The numerical method for the hydrodynamic part is a classical finite volume scheme for gas dynamics equations written under the Lagrangian formalism. It consists of a Godunov-type method in which the numerical fluxes are evaluated by solving one-dimensional approximated Riemann problems in the direction of the unit outward normal to the cell faces. The free-Lagrange mesh, which consists of a Voronoi tessellation, and cell-centered formulation also allows for general topology and arbitrary and automatically updated connectivity.

In contrast to the above hypoelastic-based formulations, hyperelastic-based formulations rely on a constitutive law which is thermodynamically consistent and satisfies automatically the principle of material objectivity. In this approach, constitutive laws are derived following the Colemann-Noll procedure [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF] which takes as a basic hypothesis the requirements that all constitutive processes be consistent with the second law of thermodynamics and the principle of frame-indifference. Being given a free-energy potential which depends on the deformation gradient, this requirement yields a constitutive law for which the Cauchy stress tensor derives directly from the free-energy function, thus leading to a system of equations written under conservative form. These good theoretical properties allow the construction of numerical methods based on approximate Riemann solvers and explain the high level of attention paid in the recent literature to the numerical schemes devoted to the discretization of hyperelastic models written either under Eulerian or Lagrangian formalism. The interesting reader might refer to [START_REF] Miller | A high-order Eulerian Godunov method for elastic-plastic flow in solids[END_REF][START_REF] Gavrilyuk | Modelling wave dynamics of compressible elastic materials[END_REF][START_REF] Kluth | Discretization of hyperelasticity on unstructured mesh with a cell-centered lagrangian scheme[END_REF][START_REF] Hill | An Eulerian hybrid WENO centered-difference solver for elastic-plastic solids[END_REF][START_REF] Favrie | Diffuse interface model for compressible fluid-compressible elastic-plastic solid interaction[END_REF] and the references therein.

However, we want to emphasize that, in spite of the previously mentioned drawbacks, hypoelastic models are able to reproduce easily laboratory data where the data experiments take the form of measurements of changes in stress with respect to changes in strain. That is why they are extensively used by many researchers and engineers and are currently used in many commercial computational codes such as LS-DYNA, refer to http://www.lstc.com/products/ls-dyna.

The numerical scheme described in this paper is based on a compatible cell-centered discretization for which the Geometric Conservation Law is rigorously satisfied regardless the shape of the cell. That is, this amounts to derive a discrete divergence operator with respect to the nodal velocities. It ensures the compatibility between the time rate of change of the cell volume and the displacement of its vertices according to the trajectory equation. Bearing this in mind, we define a mimetic discretization of the velocity gradient tensor that allows to deduce compatible discretization of the strain rate and spin tensors with respect to the nodal velocities.

Momentum and total energy equations are discretized by means of sub-cell forces. More precisely, in each cell, momentum flux is written as the summation of the sub-cell forces that are acting at each vertex of the cell. We note that the sub-cell force concept has been firstly introduced in staggered discretization [START_REF] Caramana | The construction of compatible hydrodynamics algorithms utilizing conservation of total energy[END_REF] to derive compatible conservative staggered schemes. In the present work, the sub-cell force corresponds to the integral of the stress tensor divergence over a sub-cell, knowing that a sub-cell is the quadrilateral obtained by joining the cell center, a particular vertex and the midpoints of the two edges impinging at this vertex. The total energy flux is simply deduced from the momentum flux by dot-multiplying sub-cell force by its corresponding vertex velocity. To achieve the discretization, it remains to construct an approximation of the sub-cell force and to compute the vertex velocity. The former task is completed by deriving an expression of the sub-cell force that is consistent with thermodynamics. Namely, after computing the time rate of change of specific internal energy within a cell using the semi-discrete conservation equations, we deduce a general way of writing the sub-cell force so that a cell dissipation inequality is satisfied. In this manner, the sub-cell force writes as a stress tensor contribution plus a tensorial viscous contribution, which is proportional to the difference between the cell-centered and the vertex velocities. To satisfy the dissipation inequality, the local sub-cell tensor involved in the viscous part of the sub-cell force must be positive definite. This tensor is the cornerstone of the scheme. A particular expression of this tensor is proposed by employing an approximation of the acoustic tensor deduced from a plane-wave analysis of the Wilkins model. We note that this tensor, in the hydrodynamic limit, allows to recover the scheme initially introduced in [START_REF] Maire | A cell-centered Lagrangian scheme for twodimensional compressible flow problems[END_REF], and to make the link with a node-centered approximate Riemann solver. The vertex velocity is computed by invoking the global conservation of the total energy. This last statement amounts to write that the summation of the sub-cell forces over the sub-cells surrounding a vertex is equal to zero. This balance equation also ensures the momentum conservation and leads to a node-centered solver that uniquely determines the nodal velocity.

The high-order extension of this cell-centered scheme is achieved using a one-step time discretization, wherein the fluxes are computed by means of a Taylor expansion. The time derivatives of the fluxes are obtained through the use of a node-centered solver [START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF][START_REF] Maire | A high-order one-step sub-cell force-based discretization for cell-centered lagrangian hydrodynamics on polygonal grids[END_REF] which can be viewed as a two-dimensional extension of the Generalized Riemann Problem (GRP) methodology introduced by Ben-Artzi and Falcovitz [START_REF] Ben-Artzi | Generalized Riemann problems in Computational Fluids Dynamics[END_REF].

The robustness and the accuracy of this finite volume cell-centered scheme against several relevant test cases are demonstrated.

The remainder of the paper is structured as follows: the main features of the elastic-plastic model are presented in Section 2. The first-order spatial discretization is derived in Section 3 employing the concept of sub-cell force. The nodal velocity required to move the mesh is computed by means of a node-centered solver which is introduced in Section 4. In Section 5, the second-order space and time discretization, which relies on the Generalized Riemann Problem methodology, is thoroughly exposed. To assess the robustness and the accuracy of the numerical method some numerical experiments are reported in Section 6. Concluding remarks and perspectives are given in Section 7. We recall the definition and the main properties of the exponential of a tensor in Appendix A. Finally, in Appendix B, a plane-wave analysis of the Wilkins model shows that this model, in the elastic regime, has characteristic speeds which are always real, thus ensuring the hyperbolic feature of this model.

2 Elastic-plastic model

Governing equations

Following Wilkins [START_REF] Wilkins | Methods in Computational Physics, volume 3, chapter Calculation of Elastic-Plastic Flow[END_REF], the multi-dimensional elastic-plastic flow equations in Lagrangian form write

ρ d dt ( 1 ρ ) -∇ • V = 0, (1a)
ρ d dt V -∇ • T = 0, (1b)
ρ d dt E -∇ • (TV ) = 0. ( 1c
)
Here, d dt denotes the material derivative, and ρ, V , E and T are the mass density, velocity, specific total energy and the Cauchy stress tensor which is symmetric, i.e., T t = T, where the superscript t denotes the transpose operator. The stress tensor is decomposed into a hydrostatic component P and a deviator component S by setting

T = -P I d + S, (2)
Inria where, P denotes the hydrodynamic pressure, I d is the identity tensor and S is the deviatoric part of the stress tensor, characterized by tr S = 0, where tr denotes the trace operator. Let us notice that the deviatoric stress tensor is also symmetric, that is S t = S. System (1) requires a closure to define expressions of the pressure and the deviatoric stress. This will be achieved in the next paragraph by introducing constitutive equations.

Constitutive equations

The thermodynamic closure of system ( 1) is achieved by means of an equation of state and an incremental constitutive law which describe elastic, elastic-plastic and hydrodynamic regimes.

An appropriate yield criteria is also included in the latter two regimes.

Equation of state

Following [START_REF] Wilkins | Methods in Computational Physics, volume 3, chapter Calculation of Elastic-Plastic Flow[END_REF], the equation of state for the hydrodynamic component is introduced by expressing the pressure in terms of the density and the specific internal energy under the form

P = P(ρ, ε), (3)
where ε = E -1 2 V 2 denotes the specific internal energy.

Incremental constitutive law

The evolution equation for the deviatoric stress, S, is given by an incremental constitutive law which applies to elastic perfectly plastic material. First, let us introduce some notation that will be useful for the expression of the constitutive equation. We denote by D the strain rate tensor, i.e., the symmetric part of the velocity gradient, defined by

D = 1 2 ∇V + (∇V ) t . (4)
We remark that the trace of the strain rate tensor is equal to the divergence of the velocity field, i.e., tr

(D) = ∇ • V .
In what follows, we make use of the conventional theory of plasticity which is applicable to a wide class of materials under the assumption of small deformations, refer to [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF]. In this framework, the strain rate tensor, D, admits the additive decomposition

D = D e + D p , (5)
where D e is the elastic strain rate and D p is the plastic strain rate. D e represents the rate of deformation of the underlying microscopic structure of material, whereas D p represents the rate of the local deformation due to the formation and motion of dislocations through that structure.

Finally, consistent with the general observation that the flow of dislocations does not induce changes in volume, we assume that the plastic strain rate is deviatoric

tr(D p ) = 0. ( 6
)
Moreover, we assume the dissipation inequality [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF] S : D p ≥ 0, [START_REF] Bower | Applied Mechanics of Solids[END_REF] where : denotes the inner products of tensors defined by R : Q = tr(R t Q) for any arbitrary pair of tensors R and Q. The dissipation inequality [START_REF] Bower | Applied Mechanics of Solids[END_REF], which characterizes the energy dissipated during plastic flow, will be useful for deriving our numerical scheme. Bearing this notation in mind, the incremental constitutive law for the deviatoric stress writes as

d dt S = 2µ(D 0 -D p ) -(SW -WS). (8)
Here, µ is the Lamé coefficient which is material dependent, D 0 denotes the deviatoric part of the strain rate tensor, D, which reads

D 0 = D - 1 3 tr(D)I d .
Finally, W corresponds to the antisymmetric part of the velocity gradient

W = 1 2 ∇V -(∇V ) t . (9)
The deviatoric stress, S, is subjected to the plastic boundedness inequality

f =| S | - 2 3 Y 0 ≤ 0. ( 10
)
Here, f is the yield function, | S |= (S : S) is the magnitude of the deviatoric stress and Y 0 is the flow resistance (yield strength) which is a constant in the case of elastic perfectly plastic materials. In this case, the yield stress level does not depend in any way on the degree of plastification. Thus, the yield surface remains fixed regardless of any deformation process the material may experience. This inequality corresponds to the von Mises yield condition. Finally, the plastic strain rate is determined through the equation

D p = χ(N p : D)N p . (11)
In the space of symmetric, deviatoric tensors,

N p = S | S |
represents the plastic flow direction according to the terminology introduced in [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF]. Moreover, whenever the plastic flow direction N p is mentioned, it is assumed that S = 0. In formula [START_REF] Després | A geometrical approach to non-conservative shocks and elastoplastic shocks[END_REF], χ is the switching parameter defined by

χ = 0 if f < 0 or if f = 0 and (N p : D) ≤ 0, 1 if f = 0 and (N p : D) > 0. ( 12
)
Let us notice that the case f < 0 corresponds to a behavior in the elastic range whereas the case f = 0 corresponds to the following behaviors depending on the sign of (N p : D)

• (N p : D) < 0: elastic unloading,

• (N p : D) = 0: neutral loading,

• (N p : D) > 0: plastic loading.

Inria

A simple computation shows that N p represents also the unit outward normal to the yield surface f = 0, i.e., N p = ∂f ∂S , where ∂f ∂S denotes the gradient of the scalar valued function, f , with respect to the tensor S.

Let us check the consistency of the flow rule [START_REF] Després | A geometrical approach to non-conservative shocks and elastoplastic shocks[END_REF] with respect to the plastic boundedness inequality [START_REF] Desideri | Compressible flow using unstructured grids[END_REF] Gathering the previous results leads to

S : dS dt = 2µ(1 -χ)(S : D), (13)
where the switching parameter, χ, is defined by [START_REF] Donea | Encyclopedia of Computational Mechanics[END_REF]. Finally, the time rate of change of the yield function writes

df dt = 2µ(1 -χ)(N p : D), (14)
where the plastic flow direction is given by N p = S |S| . If the yield condition is satisfied, i.e., f = 0, then by [START_REF] Donea | Encyclopedia of Computational Mechanics[END_REF] we obtain the following alternative: either (N p : D) ≤ 0, which corresponds to elastic unloading or neutral loading, hence χ = 0 and df dt ≤ 0 or (N p : D) > 0, which corresponds to plastic loading, hence χ = 1 and df dt = 0. Let us emphasize that a more detailed treatment of this question is beyond the scope of the present paper. For a complete description of the plasticity modeling using a rigorous mathematical framework, the interested reader should refer to [START_REF] Simó | Computational Inelasticity[END_REF].

Comment 1. The previous elastic perfectly plastic model can be easily enriched to take into account strain-hardening, a phenomenon characterized by an increase in flow resistance with plastic deformation, refer for instance to [START_REF] Bower | Applied Mechanics of Solids[END_REF]. In such an approach, the flow resistance Y depends on a hardening variable which is usually defined to be the accumulated plastic strain defined by

e p (t) = t 0 | D p (τ ) | dτ.
In this case, the flow rule [START_REF] Després | A geometrical approach to non-conservative shocks and elastoplastic shocks[END_REF] is modified to become

D p = χβ(e p )(N p : D)N p , ( 15
)
where β is the stiffness ratio given by

β(e p ) = 2µ 2µ + H(e p ) ,
and H is the hardening modulus defined by H(e p ) = dY de p . The switching parameter, χ, is given as previously by [START_REF] Donea | Encyclopedia of Computational Mechanics[END_REF].

The Jaumann derivative

The last term in the incremental stress-strain relationship [START_REF] Burton | Multidimensional Discretization of Conservation Laws for Unstructured Polyhedral Grids[END_REF] ensures that the constitutive law obeys the principle of material frame independence [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF]. This principle, also known as the principle of material objectivity, states that the material response is independant of the observer. This means that the tensor-valued functions that relate stress to deformation measure must transform correctly under a change of basis and change of origin for the coordinate system. More precisely, under the change of frame characterized by the rotation, Q(t), the principle of material objectivity requires that the deviatoric stress tensor, S, transforms according to the rule

S(t) → Q(t)S(t)Q t (t).
Bearing this in mind, if DS Dt denotes a generic stress rate, it is said to be objective if it transforms according to the same rule, i.e.,

D Dt S → Q(t) DS Dt Q t (t).
It is well known that the material derivative of the deviatoric stress, d dt S, is not an objective stress rate. To ensure the objectivity of the formulation of the constitutive law expressed in terms of the stress rate, it is crucial to employ an objective stress rate knowing that it exists many different objective stress rates, refer to [START_REF] De Souza Neto | Computational Methods for Plasticity: Theory and Applications[END_REF]. Here, we shall employ the Jaumann rate or Jaumann derivative which reads

S = d dt S + SW -WS. (16)
The use of this objective stress rate leads to the formulation of a stress rate strain rate constitutive law, [START_REF] Burton | Multidimensional Discretization of Conservation Laws for Unstructured Polyhedral Grids[END_REF], which satisfies the principle of material objectivity in the sense that has been previously given, refer to [START_REF] De Souza Neto | Computational Methods for Plasticity: Theory and Applications[END_REF][START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF] for a detailed demonstration. At this point it is interesting to give an interpretation of the Jaumann stress rate that will be usefull for the subsequent time discretization. To this end, we follow the presentation introduced in [START_REF] Trangenstein | A High-Order Godunov Method for Modeling Finite Deformation in Elastic-Plastic Solids[END_REF]. Let Ω be the unique tensor solution to the initial valued problem

d dt Ω = WΩ, Ω(0) = Id. (17)
Applying the transpose operator to the above equation and recalling that W is skew symmetric, i.e., W t = -W, leads to

d dt Ω t = -Ω t W, Ω t (0) = Id. (18)
Now, we compute the time derivative of Ω t Ω employing ( 17) and ( 18) to obtain

d dt (Ω t Ω) = d dt (Ω t )Ω + Ω t d dt (Ω) = -Ω t WΩ + Ω t WΩ =0.
Hence, Ω t Ω(t) = Ω t Ω(0) = Id for all t > 0, thus Ω(t) is an orthogonal tensor. More precisely, Ω(t) is a rotation since det[Ω(t)] = 1. Bearing this in mind, let us recall that ΩSΩ t is the image of the deviatoric stress, S, by the rotation Ω. Accordingly, Ω t SΩ is the pre-image of S by the rotation Ω. Finally, utilizing the previous results, it is straightforward to compute the time derivative of Inria the pre-image of the deviatoric stress to get

d dt (Ω t SΩ) = d dt (Ω t )SΩ + Ω t d dt (S)Ω + Ω t S d dt (Ω) = -Ω t WSΩ + Ω t d dt (S)Ω + Ω t SWΩ.
By rotating the above equation and comparing it to the definition of the Jaumann stress rate (16) leads to

S = Ω d dt (Ω t SΩ) Ω t . ( 19
)
This shows that the Jaumann rate of the deviatoric stress, S, is the image by the rotation, Ω, of the time derivative of the pre-image of S by the same rotation. This formula will be used to derive the time discretization of the constitutive law in the elastic regime, refer to Section 5.2.1.

Some comments about thermodynamics

The previous hypo-elastic model has been successfully employed over the past decades since the seminal work of Wilkins [START_REF] Wilkins | Methods in Computational Physics, volume 3, chapter Calculation of Elastic-Plastic Flow[END_REF]. However, it suffers from some deficiencies since it does not rely on rigorous thermodynamics basis, even in the elastic regime. Namely, the Cauchy stress tensor does not derive from a thermodynamical potential but its deviatoric part is the solution of an evolution equation. Indeed, for reversible elastic transformations, the model does not follow an isentropic evolution. Therefore, this model is thermodynamically inconsistent as it has been noticed in [START_REF] Gavrilyuk | Modelling wave dynamics of compressible elastic materials[END_REF]. From a theoretical point of view, this flaw renders quite difficult the mathematical study of this model. To investigate further its thermodynamical behavior, let us compute the time rate of change of the specific entropy, η. To this end, let us recall the fundamental Gibbs relation that defines specific entropy

θdη = dε + P d( 1 ρ ), (20)
where θ > 0 denotes the temperature. Let us assume that the specific internal energy, ε, is given by ε = E -1 2 V 2 . The time rate of change of specific entropy is obtained by means of

ρθ dη dt =ρ dE dt -ρV • d dt V + P ρ d dt ( 1 ρ ) =∇ • (TV ) -(∇ • T) • V + P ∇ • V thanks to (1)
=T : ∇V + P tr(D)

=T : D + P tr(D).

Here, we have made the assumption that the flow variables were sufficiently smooth. We have also employed the identity [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF] ∇

• (QA) = A • ∇ • (Q t ) + Q t : ∇A,
where Q is a tensor and A a vector. Finally, recalling that the stress tensor is decomposed as T = -P Id + S leads to the following equation for the time rate of change of entropy

ρθ dη dt = S : D, ( 21
)
where D is the strain rate tensor and S is the deviatoric stress solution of [START_REF] Burton | Multidimensional Discretization of Conservation Laws for Unstructured Polyhedral Grids[END_REF]. The right-hand side of ( 21) is the work rate of the deviatoric stress which only cancels for uniform elastic flows or pure hydrodynamic flows, i.e. when S = 0. We conclude that even in the elastic regime for smooth flows the present model is characterized by a non-zero entropy production which is not correct from a thermodynamical point of view. To correct this flaw, one may modify the definition of specific internal energy to take into account the elastic energy related to shear stress [START_REF] Kulikowskii | Mathematical aspects of numerical solutions of hyperbolic systems[END_REF] ε e = 1 4µρ (S : S).

Using [START_REF] Dukowicz | Accurate conservative remapping (rezoning) for arbitrary lagrangian-eulerian computations[END_REF], it is clear that in the elastic range

d dt (S : S) = 4µ(S : D), (23)
hence, using the definition of the elastic energy related to shear stress ( 22)

1 4µ d dt (4µρε e ) = (S : D).
Thus, if we set In the elastic range, by virtue of [START_REF] Desideri | Compressible flow using unstructured grids[END_REF], the first term in the right-hand side of the above equation satisfies the following inequality 0 ≤ (S : S)

E = ε + ε e + 1 2 V 2 ,
4µ 2 ρ ≤ 1 6ρ Y 0 µ 2 .
For most solids the shear modulus, µ, is two or three orders of magnitude larger than the yield strength, Y 0 . For instance, for aluminum we have µ = 27.6 10 9 Pa and Y 0 = 300 10 6 Pa [START_REF] Howell | A free-Lagrange augmented Godunov method for the simulation of elastic-plastic solids[END_REF]. Moreover, the density of a solid material is a rather high number, for instance for aluminum the initial density is ρ 0 = 2785 kg m -3 . Thus, provided that d dt (µρ) remains bounded, the right-hand side of the entropy equation ( 23) should be very small.

To provide a rigorous mathematical analysis of this model, Gavrilyuk et al. [START_REF] Gavrilyuk | Modelling wave dynamics of compressible elastic materials[END_REF] have proposed to modify it, invoking the fact that in absence of shock waves elastic transformations are reversible, the model has to be isentropic. To this end, they added the following extra evolution equation

d dt (µρ) = 0. ( 24
)
This means that the product of the density of the material by its shear modulus is conserved along the particle paths. For the same reasons and the same goals, Després [START_REF] Després | A geometrical approach to non-conservative shocks and elastoplastic shocks[END_REF] has proposed the stronger assumption that µρ remains constant, which is compatible with [START_REF] Li | Remark on the generalized Riemann problem method for compressible fluid flows[END_REF]. Let us note that, from a physical point of view, these assumptions are neither completely justified nor satisfying. Finally, we point out that Le Tallec [START_REF] Tallec | Viscous and elastoplastic constitutive laws for hydrocodes: an energetic perspective[END_REF] has proposed an alternative approach for which the Inria derivation of large strain elastic-plastic constitutive laws is performed by means of a thermodynamically consistent framework. The small strain limit of this model is characterized by an incremental constitutive law which describes the time evolution of the deviatoric elastic strain. The deviatoric stress is expressed by means of the Kirchoff stress tensor, i.e., ρ0 ρ T, where ρ 0 denotes the initial density.

In the present work, we will keep the original form of the model defined by the system of equations [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF]. Namely, we define the internal energy as ε = E -1 2 V 2 . This means that we are neglecting the elastic energy related to shear stress, ε e . However, we want to emphasize that this model, in spite of its thermodynamical inconsistency, is largely used in engineering sciences not only because it allows to incorporate easily plasticity effects but also because of its accuracy in the domain of solid dynamics for metals. Moreover, the aim of this paper is not to discuss the weakness of this model but rather to describe a new numerical method devoted to its discretization.

We conclude this paragraph by exhibiting an inequality which will be useful for the construction of the numerical scheme. First, let us recall that the plasticity model previously described relies on the dissipation inequality [START_REF] Bower | Applied Mechanics of Solids[END_REF] which reads S : D p ≥ 0, where D p is the plastic strain rate. For sufficiently smooth flow variables, we claim that the previous elastic plastic model satisfies

ρ d dt ε -T : D e ≥ 0, (25)
where D e is the elastic strain rate. To prove this inequality, we compute the time rate of change of specific internal energy. To this end, we write the time rate of change of kinetic energy by dot multiplying the momentum equation (1b) by V and then subtracting it to the total energy equation (1c) leads to

ρ d dt ε -T : D = 0. ( 26
)
Now, substituting the additive decomposition (5) of the strain rate tensor into the previous equation yields

ρ d dt ε -T : D e = S : D p .
Here, we have also used the decomposition of the Cauchy stress tensor (2) and the assumption that the plastic strain rate, D p , is deviatoric, i.e., tr(D p ) = 0. Finally, using the dissipation inequality [START_REF] Bower | Applied Mechanics of Solids[END_REF] leads to the inequality [START_REF] Loubère | Volume consistency in a staggered Lagrangian hydrodynamics scheme[END_REF], which ends the proof.

Spatial discretization

This section aims at providing a spatial discretization of system (1) using the concepts of sub-cell force and compatible discretization which have been initially introduced in [START_REF] Caramana | The construction of compatible hydrodynamics algorithms utilizing conservation of total energy[END_REF] for Lagrangian hydrodynamics in the framework of staggered discretization.

Control volume form of the elastic-plastic equations

Let D be an open set of R 2 , filled with an elastic plastic material and equipped with the orthonormal frame (0, x, y) and the orthonormal basis (e x , e y ) which is naturally completed by the unit vector e z = e x × e y . We want to develop a finite volume discretization of the elastic-plastic equations [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF]. To this end, we first recall the control volume formulation of these equations. This formulation, which consists in following a control volume moving with the material velocity, writes

d dt ω(t) dv - ∂ω(t) V • n ds = 0, (27a)
d dt ω(t) ρV dv - ∂ω(t)
Tn ds = 0, (27b)

d dt ω(t) ρE dv - ∂ω(t) (TV ) • n ds = 0. ( 27c
)
Here, ω(t) is the moving control volume and ∂ω(t) its boundary surface, whose unit outward normal is n. Equations (27a)-(27c) express the conservation of volume, momentum and total energy. Obviously, the mass of the control volume ω(t) is conserved

d dt ω(t) ρ dv = 0. ( 28
)
We note that first equation of system [START_REF] Luttwak | Slope limiting for vectors: A novel vector limiting algorithm[END_REF] expresses the time rate of change of the volume of the fluid and is often named Geometric Conservation Law (GCL). It is strongly linked to the trajectory equation

d dt x = V (x(t), t), x(0) = X, (29)
where x = x(t) is the position vector of a material point at time t > 0, which was initially located at X. Let us remark that the integral form of the conservation laws [START_REF] Luttwak | Slope limiting for vectors: A novel vector limiting algorithm[END_REF] has been obtained utilizing the following transport relation

d dt ω(t) ρφ dv = ω(t) ρ dφ dt dv, (30)
where φ = φ(x, t) denotes a physical variable. Transport relation [START_REF] Maire | A high-order one-step sub-cell force-based discretization for cell-centered lagrangian hydrodynamics on polygonal grids[END_REF] is a consequence of the Reynolds transport relation [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF] and mass conservation. Bearing this in mind, the integral form of the incremental stress-strain relationship results from multiplying (8) by the mass density and integrating over the control volume ω(t)

d dt ω(t) ρS dv - ω(t) 2ρµ(D 0 -D p ) dv + ω(t) ρ(SW -WS) dv = 0. ( 31
)

Compatible cell-centered discretization

We discretize the previous set of equations over a domain which is paved using a set of polygonal cells without gaps or overlaps. Let us introduce the notations and the assumptions that are necessary to develop our cell-centered discretization.

Notation and assumptions

Each polygonal cell is assigned a unique index c and is denoted by ω c (t). A generic point (vertex) is labeled by the index p, its corresponding position vector is x p . For a cell c, we introduce the set P(c) which is the counterclockwise ordered list of points of cell c. Conversely, for a given point p, we introduce the set C(p) containing the cells that surround point p. Being given p ∈ P(c), p -

Inria p + n p -p p - n pp + p l pp + l p -p ω c
m c d dt ( 1 ρ c ) - ∂ωc(t) V • n ds = 0, (32a)
m c d dt V c - ∂ωc(t)
Tn ds = 0, (32b)

m c d dt E c - ∂ωc(t) (TV ) • n ds = 0. ( 32c
)
Here, ∂ω c (t) is the boundary of the cell ω c (t) and m c denotes the mass of the polygonal cell, which is constant according to [START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF]. For a flow variable φ, we have defined its mass averaged value over the cell ω c (t) as

φ c = 1 m c ωc(t) ρφ dv.
We notice that the first equation corresponds to the GCL since m c = ρ c v c where v c is the volume of cell c. We have obtained a set of semi-discrete evolution equations for the primary variables ( 1 ρc , V c , E c ) whose thermodynamic closure is given by the equation of state

P c = P (ρ c , ε c ), where ε c = E c - 1 2 V 2 c . (33)
The discrete version of the incremental constitutive law and the associated flow rule read

d dt S c = 2µ c (D 0,c -D p c ) -(S c W c -W c S c ), (34a)
f c =| S c | - 2 3 Y 0 ≤ 0, (34b)
D p c = χ c (N p c : D c )N p c . (34c)
Here, D p c denotes the plastic strain, and N p c = Sc |Sc| is the plastic flow direction and χ c the switching parameter defined within cell c by

χ c = 0 if f c < 0 or if f c = 0 and (N p c : D c ) ≤ 0, 1 if f c = 0 and (N p c : D c ) > 0. ( 35
)
Let us recall that D 0,c denotes the deviatoric part of the strain rate tensor, D c , that is

D 0,c = D c - 1 3 tr(D c )I d .
The motion of the grid is ruled by the semi-discrete trajectory equation written at each point

d dt x p = V p (x p (t), t), x p (0) = X p . (36)
To achieve the discretization, it remains not only to compute the numerical fluxes related to volume, momentum and total energy but also to determine the nodal velocity to compute the grid motion. We have also to provide a discretization for the piecewise constant strain rate and spin tensors D c and W c .

Geometric Conservation Law compatibility

Here, the expression GCL compatibility means that we are deriving a discrete divergence operator for the volume equation (1a) by requiring consistency of the divergence of the velocity field with the time rate of change of volume of a cell. Noticing that m c = ρ c v c , Eq. (32a) rewrites

d dt v c = ∂ωc(t) V • n ds. ( 37
)
This equations states that the time rate of change of the cell volume v c must be equal to the volume swept by the element boundary during its displacement with the flow velocity. Assuming that the volume v c (t) in the left-hand side can be computed exactly, this amounts to requiring the exact computation of the volume flux in the right-hand side also. In this manner, we obtain a compatible discretization of the volume flux. Using the triangulation displayed in Fig. 2, the polygonal cell volume writes

v c (t) = 1 2 p∈P(c) x p (t) × x p + (t) • e z . (38)
We remark that the volume is expressed as a function of the position vectors of its vertices.

Using the chain rule derivative, time differentiation of Eq. ( 38) yields

d dt v c = p∈P(c) ∇ xp v c • d dt x p ,
where ∇ xp v c is the gradient of the cell volume according to the position vector x p . This gradient is computed directly through the use of ( 38) We define the normal vector at corner pc by setting

∇ xp v c = 1 2 ∇ xp (x p + -x p -) × e z • x p = 1 2 (x p + -x p -) × e z . Inria p + n p -p p - O n pp + U p ω c p l pp + l p -p l pc n pc
l pc n pc = ∇ xp v c = 1 2 (x p + -x p -) × e z , (39)
where the length, l pc , is determined knowing that n 2 pc = 1. We notice that this corner vector can be expressed in terms of the two outward normals impinging at node p as

l pc n pc = 1 2 (l p -p n p -p + l pp + n pp + ), (40)
where l p -p and n p -p are the length and the unit outward normal related to the edge [p -, p] whereas l pp + and n pp + are the length and the unit outward normal related to the edge [p, p + ], refer to Fig. 1. It turns out that the corner vector l pc n pc is the fundamental geometric object that allows to define uniquely the time rate of change of cell volume as

d dt v c = p∈P(c) l pc n pc • V p . (41)
Here, we have used the trajectory equation [START_REF] Menikoff | Equations of state and fluid dynamics[END_REF], i.e. d dt x p = V p . The corner vector l pc n pc satisfies the fundamental geometrical identity p∈P(c)

l pc n pc = 0, (42)
which is equivalent, using [START_REF] Sambasivan | Exploration of new limiter schemes for stress tensors in Lagrangian and ALE hydrocodes[END_REF], to the well known result that the summation of the outward normals to a closed polygonal contour is equal to zero. This also involves that the volume of a Lagrangian cell in Cartesian geometry moving in a uniform flow does not change. We claim that with this purely geometric derivation we have completely defined the volume flux. Moreover, this definition is by construction compatible with mesh motion. We also remark that this result can be used to derive the discrete divergence operator, DIV c over cell c as

DIV c (V ) = 1 v c ∂ωc(t) V • n ds.
Now, combining [START_REF] Miller | A high-order Eulerian Godunov method for elastic-plastic flow in solids[END_REF] and ( 41) we get

DIV c (V ) = 1 v c d dt v c = 1 v c p∈P(c) l pc n pc • V p . (43)
Let us remark that property [START_REF] Tallec | Viscous and elastoplastic constitutive laws for hydrocodes: an energetic perspective[END_REF] shows that the null space of the discrete divergence operator is spanned by the set of constant vectors. We recover the compatible discretization of the divergence operator currently used in the derivation of staggered grid Lagrangian schemes for compressible hydrodynamics [START_REF] Burton | Multidimensional Discretization of Conservation Laws for Unstructured Polyhedral Grids[END_REF][START_REF] Whalen | Algebraic limitations on two dimensionnal hydrodynamics simulations[END_REF][START_REF] Caramana | The construction of compatible hydrodynamics algorithms utilizing conservation of total energy[END_REF].

Comment 2. We want to mention that the above discrete divergence operator could have been derived utilizing an alternative approach which relies on the approximation of the surface integral, ∂ωc(t) V • n ds, wherein a piecewise linear distribution of the velocity over each edge is assumed; this piecewise linear distribution being determined from the knowledge of the nodal velocity. Let us notice that this approach has been introduced a long time ago in [START_REF] Noh | Methods in Computational Physics, volume 3, chapter CEL: A Time-dependent Two-Space-Dimensional Coupled Eulerian-Lagrangian Code[END_REF] to compute the difference approximation for partial derivatives over an arbitrary region.

Using the same approach, it is straightforward to define the discrete gradient operator, GRAD c , over cell c. Let Q be a scalar valued function over cell c, we define it as

GRAD c (Q) = 1 v c ∂ωc(t)
Qn ds.

The corresponding practical formula writes

GRAD c (Q) = 1 v c p∈P(c) Q p l pc n pc . (44)
where Q p is the value taken by Q at the node p. We claim that the above formula is exact for linear functions. To demonstrate this result, let us introduce the useful geometrical identity

p∈P(c) l pc (n pc ⊗ x p ) = v c I d . (45)
This result can be demonstrated directly using simple algebra [START_REF] Maire | Contribution to the numerical modeling of inertial confinement fusion[END_REF] or applying the divergence theorem as follows [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF]

v c I d = ωc ∇x dv = ∂ωc x ⊗ n ds.
Bearing this in mind, we can demonstrate that ( 44) is exact for linear functions. To this end, let us consider the linear function, Q = A + B • x, where A and B are a constant scalar and a constant vector. Obviously the nodal value of

Q reads Q p = A + B • x p . Substituting this expression into (44) leads to v c GRAD c (Q) = p∈P(c) (A + B • x p )l pc n pc = A p∈P(c) l pc n pc + p∈P(c) l pc (n pc ⊗ x p )B, recalling that (u ⊗ v)w = u(v • w) = v c I d B, thanks to (42) and (45) = v c B.
Inria

Strain rate and spin tensors compatible discretization

We construct a compatible cell-centered discretization for the strain rate and spin tensors. To this end we start by constructing a compatible discretization for the velocity gradient. This amounts to define a discrete gradient operator over the polygonal cell c using the divergence theorem as follows

GRAD c (V ) = 1 v c ∂ωc V ⊗ n ds.
This definition leads to the following practical formula for the discrete gradient operator over a polygonal cell

GRAD c (V ) = 1 v c p∈P(c) l pc (V p ⊗ n pc ). ( 46
)
We claim that this formula is exact for linear velocity field. To demonstrate this, let us introduce the following linear velocity field, V (x) = A + Bx, where A and B are respectively a constant vector and a second order constant tensor. Substituting this linear velocity field into (46) yields

v c GRAD c (V ) = p∈P(c) l pc (A ⊗ n pc ) + l pc (Bx p ⊗ n pc ) = p∈P(c) Bl pc (x p ⊗ n pc ), thanks to (42) =Bv c I d .
Finally, the geometric identity (45) leads to GRAD c (V ) = B, which ends the proof. This shows that the present discretization is able to capture exactly rigid body motion.

Recalling that the strain rate tensor and the spin tensor are respectively the symmetric and anti-symmetric part of the velocity gradient, it is straightforward to define their discrete counterpart as

D c = 1 2 GRAD c (V ) + GRAD t c (V ) , W c = 1 2 GRAD c (V ) -GRAD t c (V ) .
Therefore, the cell-averaged values of the strain rate tensor and the spin tensor are given by

D c = 1 v c p∈P(c) 1 2 l pc (V p ⊗ n pc + n pc ⊗ V p ) , (47)
W c = 1 v c p∈P(c) 1 2 l pc (V p ⊗ n pc -n pc ⊗ V p ) . (48)
Knowing that tr(A ⊗ B) = A • B, where A and B are two vectors and taking the trace of Eq. ( 47) we obtain

tr(D c ) = 1 v c p∈P(c) l pc n pc • V p .
The right-hand side of the above equation is nothing but the discrete divergence operator of the velocity field, that is tr(D c ) = DIV c (V ), which is the discrete counterpart of the identity tr D = ∇ • V . We conclude by claiming that we have construct a compatible discretization of the strain rate tensor. We also point out that the averaged spin tensor W c is traceless, that is tr(W c ) = 0.

Comment 3. It is important to notice that for a uniform flow, i.e. V p = V 0 , ∀ p, the strain rate tensor and the spin tensor are equal to zero thanks to the geometrical identity (42).

Compatible discretization of a tensor identity

In this paragraph, in view of the subsequent spatial discretization, we shall investigate the compatible discretization of the following tensor identity

∇ • (QA) = A • ∇ • (Q t ) + Q t : ∇A, ( 49
)
where Q is an arbitrary tensor and A an arbitrary vector. This identity gives the counterpart for fields, of the standard product rule for scalar functions of scalar variable. Let us note that this identity, when applied to the stress tensor, Q = T and the velocity field A = V , enables to derive the specific internal energy equation by subtracting the kinetic energy equation to the total energy equation. Knowing that the kinetic energy equation results from dot multiplying the momentum equation by the velocity vector. Here, we complete the work initiated in the last two paragraphs, by extending the compatible discretization of the divergence operator to the case of tensor fields. This amounts to define a discrete divergence operator over the polygonal cell c using the divergence theorem as follows

DIV c (Q) = 1 v c ∂ωc Qn ds.
This definition leads to the following practical formula for the discrete divergence operator over a polygonal cell

DIV c (Q) = 1 v c p∈P(c) l pc Q pc n pc . (50)
We refer to Q pc as the pointwise value of Q at the node p of cell c. Let us consider the product QA, which is a vector, and let us define its pointwise value at node p by (QA) p = Q pc A p where A p is the pointwise value of the vector A at node p. Employing the definition of the discrete divergence operator for a vector field (43) leads to

DIV c (QA) = 1 v c p∈P(c) l pc (Q pc A p ) • n pc .
Now, combining (50), the above formula and the discrete gradient operator for vector fields [START_REF] Whalen | Algebraic limitations on two dimensionnal hydrodynamics simulations[END_REF] yields

v c DIV c (QA) -A c • DIV c (Q t ) -Q t c : GRAD c (A) = p∈P(c) l pc [(Q pc -Q c )(A p -A c )] • n pc .
(51) Here, Q c and A c denote the mean values of Q and A over cell c. This result has been obtained using the identity Q : (A ⊗ B) = A • (QB), where B is an arbitrary vector and also the geometric identity [START_REF] Tallec | Viscous and elastoplastic constitutive laws for hydrocodes: an energetic perspective[END_REF]. Formula (51) consists of the discrete counterpart of the identity (49). The righthand side can be viewed as a residual term, which cancels only for uniform tensor and vector fields.

Sub-cell force-based discretization

It remains to discretize momentum and total energy equations. To this end, we utilize the fundamental concept of sub-cell force initially introduced in [START_REF] Caramana | The construction of compatible hydrodynamics algorithms utilizing conservation of total energy[END_REF] in the context of staggered discretization. Let us introduce some notations. Being given a polygonal cell, ω c , for each vertex p ∈ P(c), we define the sub-cell ω pc by connecting the centroid of ω c to the midpoints of edges [p -, p] and [p, p + ] impinging at node p, refer to Fig. 3. In two dimensions the sub-cell, as just defined, is always a quadrilateral regardless of the type of cells that compose the underlying grid. Using the sub-cell definition, cell ω c and its boundary ∂ω c can be decomposed as Here, we have made the assumption that the velocity field over the outer sub-cell boundary, ∂ω pc ∩ ∂ω c is constant and equal to the nodal velocity, V p . Finally, total energy equation (32c

ω c = p∈P(c) ω pc , ∂ω c = p∈P(c) ∂ω pc ∩ ∂ω c . (52)
) rewrites m c d dt E c = p∈P(c) F pc • V p .
This leads to the mechanical interpretation that the time variation of total energy results from the summation over the sub-cell of the rate of work done by sub-cell force F pc applied on each vertex p. Combining the previous results, system [START_REF] Maire | A cell-centered Lagrangian scheme for twodimensional compressible flow problems[END_REF] becomes

m c d dt ( 1 ρ c ) - p∈P(c) l pc n pc • V p = 0, (54a)
m c d dt V c - p∈P(c) F pc = 0, (54b)
m c d dt E c - p∈P(c) F pc • V p = 0. (54c)
The cell-centered discrete unknowns ( 1 ρc , V c , E c ) satisfy a system of semi-discrete evolution equations wherein the numerical fluxes are expressed as functions of the nodal velocity, V p , and the sub-cell force F pc . Let us recall that thermodynamic closure is given by the equation of state [START_REF] Maire | A second-order cell-centered Lagrangian scheme for twodimensional compressible flow problems[END_REF], the incremental constitutive law by [START_REF] Maire | Multi-scale Godunov-type method for cell-centered discrete Lagrangian hydrodynamics[END_REF] and the grid motion is governed by the semi-discrete trajectory equation [START_REF] Menikoff | Equations of state and fluid dynamics[END_REF].

We conclude this section by giving an interpretation of the above sub-cell force-based discretization in terms of the discrete divergence operator previously introduced. Applying the definitions ( 43), (50) of the discrete divergence operators for vectors and tensors to the system (32) directly yields

m c d dt ( 1 ρ c ) -v c DIV c (V ) = 0, (55a)
m c d dt V c -v c DIV c (T) = 0, (55b)
m c d dt E c -v c DIV c (TV ) = 0. ( 55c
)
Thanks to (50), the discrete divergence of the stress tensor reads

DIV c (T) = 1 v c p∈P(c) l pc T pc n pc ,
where T pc is the pointwise of T at point p of cell c. Bearing this in mind, the comparison between (55b) and (54b) shows that the sub-cell force, F pc , and the pointwise stress tensor, T pc , are linked by

F pc = l pc T pc n pc . ( 56
)

Inria

To complete the discretization, it remains to compute the nodal velocity and construct an approximation of the sub-cell force. This goal will be achieved by investigating the properties of the scheme regarding its approximate thermodynamic consistency and its conservation for momentum and total energy.

Thermodynamic consistency of the semi-discrete scheme

We derive a general form of the sub-cell force requiring that the semi-discrete scheme (54) satisfies a semi-discrete version of the dissipation inequality [START_REF] Loubère | Volume consistency in a staggered Lagrangian hydrodynamics scheme[END_REF]. To this end, we compute the semi-discrete version of the time rate of change of specific internal energy within cell c. First, recalling that the specific internal energy within cell c is given by

ε c = E c -1 2 V 2 c , we get m c d dt ε c = m c d dt E c -m c V c • d dt V c .
Bearing this in mind, we dot multiply the momentum equation (55b) by V c and subtract it to the total energy equation (55c) to obtain

m c d dt ε c -v c [DIV c (TV ) -V c • DIV c (T)] = 0.
Now, by means of the discrete tensor identity (51), we rewrite the above equation as

m c d dt ε c -v c T c : GRAD c (V ) = p∈P(c) l pc [(T pc -T c )n pc ] • (V p -V c ).
Here, we have used the symmetry of the Cauchy stress tensor, i.e., T t = T. Finally, recalling that the discrete strain rate tensor, D c , is defined in terms of the discrete velocity gradient tensor by

D c = 1 2 [GRAD c (V ) + GRAD t c (V )], leads to m c d dt ε c -v c (T c : D c ) = p∈P(c) l pc [(T pc -T c )n pc ] • (V p -V c ). ( 57
)
Thanks to the additive decomposition of the strain rate into elastic and plastic parts, i.e., D c = D e c + D p c and recalling that the plastic strain rate tensor is deviatoric, the previous equation rewrites as

m c d dt ε c -v c (T c : D e c ) = v c (S c : D p c ) + p∈P(c) (F pc -l pc T c n pc ) • (V p -V c ), (58)
where S c is the deviatoric stress tensor. Here, we have introduced the sub-cell force in the residual term at the right-hand side by means of (56). Let us remark that the plastic strain rate tensor is defined by the flow rule (34c). We want our discretization to be thermodynamically consistent in the sense that the semi-discrete version of the scheme must satisfy the dissipation inequality

m c d dt ε c -v c (T c : D e c ) ≥ 0. ( 59
)
Let us notice that this dissipation inequality is nothing but the discrete counterpart of [START_REF] Loubère | Volume consistency in a staggered Lagrangian hydrodynamics scheme[END_REF]. Knowing that by definition (34c), the plastic work rate, (S c : D p c ) is non negative, a sufficient condition to obtain this consists in defining the sub-cell force so that the residual term in the right-hand side of (58) remains always non negative. To this end, we set

F pc = l pc T c n pc + M pc (V p -V c ), ( 60
)
where M pc is a sub-cell-based 2 × 2 matrix which must satisfy the following properties

RR n°7975 • M pc is positive semidefinite, i.e. M pc V • V ≥ 0, ∀ V ∈ R 2 .
• M pc has a physical dimension equal to a length times a mass density times a velocity, i.e.

[M pc ] = L ρ U .

Moreover, its definition must be compatible with the principle of material frame-indifference [START_REF] Botsis | Mécanique des milieux continus[END_REF], namely it should not depend on whatever external frame of reference is used to describe it. In a nutshell, M pc is invariant by translation and transforms as M pc = RM pc R t for a rigid rotation R. Substituting (60) into (58) leads to

m c d dt ε c -v c (T c : D e c ) = v c (S c : D p c ) + p∈P(c) M pc (V p -V c ) • (V p -V c ). ( 61
)
Knowing that the (S c : D p c ) ≥ 0 and since the right-hand side is a positive semidefinite quadratic form, the dissipation inequality (59) is satisfied. We remark that the numerical dissipation rate within cell c is directly governed by the positive semidefinite sub-cell matrix M pc and the velocity jump between the nodal and the cell-centered velocity, ∆V pc = V p -V c .

We conclude this paragraph by evaluating the time rate of change of entropy associated to our semi-discrete scheme. Thanks to the Gibbs formula, i.e., θdη = dε + P d( 1ρ ), the time rate of change of entropy within cell c writes

m c θ c d dt η c = m c d dt ε c + P c d dt 1 ρ c ,
where θ c > 0 denotes the mean temperature within cell c. Replacing the time rate of change of specific internal energy by means of (61) in the above equation yields

m c θ c d dt η c = v c (T c : D c ) + P c v c tr(D c ) + p∈P(c) M pc (V p -V c ) • (V p -V c ).
Here

m c θ c d dt η c = v c (S c : D c ) + p∈P(c) M pc (V p -V c ) • (V p -V c ). ( 62
)
Knowing that the second term in the right-hand side of the above equation is always non-negative, the specific entropy satisfies the following inequality within cell c

m c θ c d dt η c ≥ v c (S c : D c ). ( 63
)
This inequality which characterizes the entropy production of our semi-discrete scheme is the discrete counterpart of the entropy production equation ( 21) derived in Section 2.3.

Comment 4. At this point it is interesting to investigate the limit of the previous discretization for a material characterized by a pure hydrodynamic behavior. In this particular case, the Cauchy stress tensor collapses to its spherical part, that is, S = 0, T = -P Id and the system (1) collapses to the system of gas dynamics equations. The plastic strain rate is equal to zero and the sub-cell force expression becomes

F pc = -l pc P c n pc + M pc (V p -V c ),

Inria

where P c denotes the piecewise constant pressure within cell c. We point out that with this expression of the sub-cell force, we recover the sub-cell force-based discretization that has been derived in [START_REF] Maire | Contribution to the numerical modeling of inertial confinement fusion[END_REF]. Furthermore, the dissipation inequality transforms into

m c d dt ε c + P c v c tr(D c ) ≥ 0.
Noticing that tr(D c ) = DIV c (V ) and substituting the volume conservation equation (55a) allows to rewrite the previous inequality as

m c d dt ε c + m c P c d dt ( 1 ρ c ) ≥ 0.
This last equation is nothing but the entropy inequality associated to the semi-discrete system of the gas dynamics equations. We note also that in the pure hydrodynamic limit, (61) collapses to the time rate of change of specific entropic

m c d dt ε c + m c P c d dt ( 1 ρ c ) = p∈P(c) M pc (V p -V c ) • (V p -V c ).
We conclude by claiming that the present spatial discretization is a natural extension of the numerical scheme that has been developed initially for Lagrangian hydrodynamics [START_REF] Maire | A high-order one-step sub-cell force-based discretization for cell-centered lagrangian hydrodynamics on polygonal grids[END_REF].

Summary of the semi-discrete scheme

Having determined the generic sub-cell force form consistent with the dissipation inequality (59), we present a summary of the semi-discrete equations that govern the evolution of the primary variables

( 1 ρc , V c , E c ) m c d dt ( 1 ρ c ) - p∈P(c) l pc n pc • V p = 0, m c d dt V c - p∈P(c) F pc = 0, m c d dt E c - p∈P(c) F pc • V p = 0.
Here, the sub-cell force reads

F pc = l pc T c n pc + M pc (V p -V c ),
where the piecewise stress tensor is decomposed into a hydrostatic component and a deviator component as T c = -P c I d + S c . The closure of the above system is ensured by means of the following equation of state

P c = P (ρ c , ε c ), where c = E c - 1 2 V 2 c .
Moreover, the behavior of the elastic plastic material is described by the following ingredients:

d dt S c = 2µ c (D 0,c -D p c ) -(S c W c -W c S c ), f c =| S c | - 2 3 Y 0 ≤ 0, D p c = χ c (N p c : D c )N p c , where N p c = S c | S c | .
Here, D c , D 0,c and W c denote respectively the piecewise constant strain rate tensor, its deviatoric part and the spin tensor which are discretized in terms of the nodal velocity as follows

D c = 1 v c p∈P(c) 1 2 l pc (V p ⊗ n pc + n pc ⊗ V p ) , D 0,c = 1 v c p∈P(c) 1 2 l pc (V p ⊗ n pc + n pc ⊗ V p ) - 1 3v c   p∈P(c) l pc n pc • V p   Id, W c = 1 v c p∈P(c) 1 2 l pc (V p ⊗ n pc -n pc ⊗ V p ) .
We have derived a compatible semi-discrete cell-centered Lagrangian scheme by means of a subcell force formalism. This scheme satisfies a semi-discrete local dissipation inequality in the sense we have defined in Section 2.3. To achieve the discretization it remains to compute the nodal velocity V p and also to provide an expression for the sub-cell matrix M pc which is the cornerstone of the scheme.

4 Node-centered solver for the grid velocity

Total energy and momentum conservation

We investigate the conservation of our cell-centered scheme regarding total energy and momentum. Total energy over the whole grid is defined as Before proceeding further, we discretize the right-hand side. To this end, let us introduce some specific notations, assuming that the boundary is a closed contour. Let p be a node located on the boundary ∂D, we denote by p -and p + the previous and next points on the boundary with respect to p in the counterclockwise ordered list of points located on ∂D. The curvilinear boundary ∂D is discretized using the decomposition ∂D = ∪ p ∂D p . Here,

E(t) = c m c E c (t),
∂D p = [i -, p] ∪ [p, i + ],
where i ± is the midpoint of the segment [p ± , p]. Using this decomposition, boundary term contribution is discretized as

∂D TV • n ds = p∈∂D ∂Dp TV • n ds = p∈∂D ∂Dp Tn ds • V p .
The term between parentheses in the right-hand side of the second line represents a corner force that acts from the exterior boundary onto boundary points. Then, it is natural to set

F p = ∂Dp Tn ds, (65)
Inria where F p is the boundary corner force acting onto point p. Combining the previous notations and substituting the specific total energy equation (54c) into (64) yields the balance of total energy over the entire domain c p∈P(c)

F pc • V p = p∈∂D F p • V p .
Interchanging the order in the double sum in the left-hand side yields

p   c∈C(p) F pc   • V p = p∈∂D F p • V p ,
where C(p) is the set of cells surrounding point p. Finally, left-hand side of the above equation is divided into two parts depending on the points location

p∈D o   c∈C(p) F pc   • V p + p∈∂D   c∈C(p) F pc   • V p = p∈∂D F p • V p , (66)
where D o is the interior of the domain D. Knowing that the total energy balance (66) must hold whatever point velocity is, total energy conservation is ensured if and only if

∀ p ∈ D o , c∈C (p)
F pc = 0, (67a)
∀ p ∈ ∂D, c∈C (p)
F pc = F p . (67b)
It remains to check that the previous conditions also lead to momentum conservation. Let Q denotes the total momentum over the entire domain, i.e. Q = c m c V c . We compute its time rate of change

d dt Q = c m c d dt V c = c p∈P(c) F pc , thanks to (54b) = p∈D o c∈C(p) F pc + p∈∂D c∈C(p)
F pc , by interchanging the double sum

= p∈∂D F p , thanks to (67) = p∈∂D ∂Dp
Tn ds, thanks to (65).

We conclude that up to the boundary terms contribution, momentum is conserved over the entire domain. Hence, conditions (67a) and (67b) turn out to ensure not only total energy but also momentum conservation. Moreover, as we shall show it in next section, they also provide a vectorial equation that enables us to determine the nodal velocity.

Comment 5. The fundamental equation (66) can be interpreted as discrete variational formulation wherein the point velocity V p is a test function. It can also be viewed as a principle of virtual work [16].

Generic formulation of the nodal solver

Remembering that the general sub-cell force form reads

F pc = l pc T c n pc + M pc (V p -V c )
, where M pc is a 2 × 2 symmetric positive semidefinite matrix, and using the conservation condition (67), we are now in position to write the 2 × 2 system that solves the nodal velocity V p . In the general case this system writes

∀ p ∈ D o , M p V p = c∈C(p) (-l pc T c n pc + M pc V c ) , (68a)
∀ p ∈ ∂D, M p V p = c∈C(p) (-l pc T c n pc + M pc V c ) + F p , (68b)
where M p denotes the 2 × 2 node-centered matrix defined as

M p = c∈C(p) M pc . (69)
We emphasize that we have divided the nodal velocity determination into two cases depending on the node location. As M pc is positive semidefinite, M p also shares the same property. To enforce the solvability of Eq. (68a) and Eq. (68b), we assume that the matrix M pc is positive definite. This yields to a M p matrix which is always invertible. Therefore, provided that the M pc matrix is defined, the nodal velocity, V p is always uniquely determined by inverting equations (68a) and (68b). We recall that M pc has the physical dimension of a length times a density times a velocity. If M pc does not depend on the nodal velocity, hence (68a) and (68b) are linear equations whose solutions are easily obtained. Conversely, if M pc depends explicitly on the nodal velocity, one has to solve non-linear equations by using an iterative method such as fixed point algorithm. In this latter case, the invariance by translation requires that M pc is expressed as a function of the difference between the cell velocity and the nodal velocity.

Practical issue related to boundary conditions

In this section, we describe briefly boundary conditions implementation for the elastic-plastic model [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF]. In Lagrangian formalism, this task is quite simple as we have to consider only two types of natural boundary conditions. Namely, on the boundary of the domain, ∂D, either the traction or the velocity is prescribed. We also describe the wall boundary condition, which is a usefull boundary condition devoted to the numerical simulation of projectile impact on a rigid wall. We present a boundary conditions implementation which is fully compatible with the node-centered solver previously developed. To this end, let us consider a generic point p located on the boundary, we denote p -and p + its previous and next neighbor with respect to p in the counterclockwise ordered list of boundary points. Without loss of generality we make the assumption that ∂D is a closed contour. The two outward normals to the edges located on the boundary and connected to point p are l - p n - p and l + p n + p , where l ± p denotes one half of the length of the segment [p, p ± ], refer to Fig. 4. These outward normals are linked to the corner vectors l pc n pc as follows c∈C(p)

l pc n pc = l - p n - p + l + p n + p . (70)
Once more, this is due to the fact that the dual cell contour is closed (red dashed line in Fig. 4). Now, let us investigate the two following cases.

Inria 000000000000000000000000 000000000000000000000000 000000000000000000000000 000000000000000000000000 000000000000000000000000 000000000000000000000000 000000000000000000000000 000000000000000000000000 000000000000000000000000 000000000000000000000000 000000000000000000000000 000000000000000000000000 • Prescribed traction: the boundary corner force acting onto point p defining by (65) writes

F p = l - p t -, p + l + p t +,
p , where t ±, p are the prescribed tractions on the two half-edges impinging at point p, that is Tn ± p = t ±, p . Substituting this expression of the boundary corner force into (68b) leads to

M p V p = c∈C(p) (-l pc T c n pc + M pc V c ) + l - p t -, p + l + p t +, p . (71)
Due to (70), we note that this formula preserves uniform fluid flows. Namely, if (T c , V c ) = (T 0 , V 0 ), ∀ c, and t ±, p = T 0 n ± p , ∀ p, then V p = V 0 .

• Prescribed velocity: this boundary condition is particularly easy to implement. Indeed, being given the velocity field distribution, V (x), on the boundary, we prescribe the nodal velocity at point p as V p = V (x p ).

• Wall boundary condition: let V ±, * p denotes the prescribed normal velocity on sides of point p, that is (V • n ± p ) = V ±, * p . We point out that most of the time, the rigid wall is at rest and characterized by a planar interface, thus V ±, * p = 0 and the unit outward normals, n ± p , are colinear. Let us note that the nodal velocity, V p , satisfies (68b)

M p V p = c∈C(p) (-l pc T c n pc + M pc V c ) + F p ,
where F p is an unknown force acting at node p. We point out that the nodal velocity is constrained by the normal velocity boundary conditions as follows

l - p n - p + l + p n + p • V p = l - p V -, * p + l + p V +, * p . ( 72
)
It remains to determine the expression of the boundary force F p . To this end, we first recall that the expression of the force element that acts on the boundary surface element nds reads dF = Tnds, where n is the unit outward normal to the boundary. Using the decomposition of the Cauchy stress tensor, T = -P Id + S, and splitting the traction between its normal and tangent components yields

dF = {[-P + (Sn • n)] n + [Id -(n ⊗ n)] (Sn)} ds.
Let us point out that [Id -(n ⊗ n)] is the orthogonal projector onto the plane tangent to the boundary surface, hence for any vector λn, where λ is an arbitrary real, we get [Id -(n ⊗ n)](λn) = 0. Now, we compute the power element, dP, of the above force element resulting from a prescribed velocity of the form

V = V n dP =dF • V = {[-P + (Sn • n)] n + [Id -(n ⊗ n)] (Sn)} • V nds = [-P + (Sn • n)] V .
This shows that only the normal component of the traction contributes to the power generated by a prescribed normal velocity field. Bearing this in mind, and recalling that we want to ensure total energy conservation (66), we propose the following expression for the boundary corner force, F p ,

F p = (l - p n - p + l + p n + p )Π p ,
where Π p stands for an averaged normal traction acting onto point p. Substituting this expression into (68b) and making use of (72) leads to

M p V p -(l - p n - p + l + p n + p )Π p = c∈C(p) -l pc T c n pc + M pc V c . (73)
Thanks to the above equation, we express the nodal velocity in terms of the normal traction Π p , then dot multiplying this resulting equation by the corner normal N p = l - p n - p + l + p n + p , and using the boundary condition (72), we finally get the equation satisfied by the normal traction, Π p

(M -1 p N p • N p )Π p = l - p V -, * p + l + p V +, * p - c∈C(p) M -1 p (-l pc T c n pc + M pc V c ) • N p . (74)
It is clear that this equation has always a unique solution provided that M -1 p is positive definite. Computing Π p by means of (74), it remains to obtain V p by solving

M p V p = c∈C(p) (-l pc T c n pc + M pc V c ) + Π p N p ,
At this point, we can conclude that the sub-cell force formalism provides a general framework that allows to construct a family of compatible cell-centered schemes that share good physical properties such as thermodynamic consistency and conservation. The numerical fluxes and the nodal velocity are computed in a compatible manner by means of a node-centered solver which enables us to consistently derive boundary conditions. The key point in designing these schemes lies in constructing the corner matrix M pc . We shall investigate this latter point in the next paragraph, by presenting one example of such a construction and making the link with approximate Riemann solvers.

Inria n - pc p - n + pc p + l + pc l - pc (ρ c , V c , T c ) ω c T - pc p V p l pc n pc = l - pc n - pc + l + pc n + pc T + pc

Definition of the corner matrix

The definition of the corner matrix is based on a methodology that has been developed initially in [START_REF] Maire | A cell-centered Lagrangian scheme for twodimensional compressible flow problems[END_REF] and revisited in [START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF][START_REF] Maire | Multi-scale Godunov-type method for cell-centered discrete Lagrangian hydrodynamics[END_REF] for pure Lagrangian hydrodynamics. Here, we present an extension of this methodology to the case of elastic flows which relies on the plane-wave analysis performed in Appendix B. Let us recall that the abstract formulation of the sub-cell force reads

F pc = l pc T c n pc + M pc (V p -V c ),
where M pc must be a 2 × 2 positive semidefinite matrix that remains to express. To this end, we set

M - pc = ρ c l - pc a L,c (n - pc ⊗ n - pc ) + a T,c Id -(n - pc ⊗ n - pc ) , (75a)
M + pc = ρ c l + pc a L,c (n + pc ⊗ n + pc ) + a T,c Id -(n + pc ⊗ n + pc ) , (75b)
M pc = M - pc + M + pc . (75c)
where l - pc n - pc and l + pc n + pc are the outward normals to edges [p -, p], [p, p + ]. Let us note that l ± pc are one half of the length of the corresponding edges, refer to Fig. 5. In the above equations, the coefficient a L,c and a T,c are non negative scalars which are defined by

a 2 L,c = a 2 c + 4 3 µ c ρ c , (76a)
a 2 T,c = µ c ρ c , (76b)
where µ c is the shear modulus of the material and a c is the thermodynamic sound speed computed from the equation of state of the material as

a 2 c = ∂P ∂ρ ε + P ρ 2 ∂P ∂ε ρ .
Let us point that a L,c and a T,c are respectively the longitudinal and transverse elastic wave speeds. We observe that the structure of the matrix M pc relies on the acoustic tensor, H, which has been obtained in Appendix B after performing a plane-wave analysis of the Wilkins model. This acoustic tensor reads

H(n) = a 2 L (n ⊗ n) + a 2 T [Id -(n ⊗ n)] ,
where n is the unit vector which defines the wave direction. It is clear that a 2 L is the eigenvalue of H associated to the eigenvector n, whereas a 2 T is the eigenvalue of H related to the eigenvector τ perpendicular to n. The matrix 1 ρcl ± pc M ± pc has been defined as being the square root of the acoustic tensor H(n ± pc ). We point out that this choice is an approximation of the full acoustic tensor of the Wilkins model which has been derived in Appendix B. The underlying approximation relies on the assumption Y 0 µ 1 which allows to neglect the contribution of the Jaumann derivative in the expression of the acoustic tensor. In a future work, we plan to investigate the sensitivity of the numerical scheme to this approximation by defining a corner matrix which relies directly on the full acoustic tensor of the Wilkins model. We also note that the corner matrix we have proposed is quite similar to one of the corner matrices that has been derived in [START_REF] Kluth | Discretization of hyperelasticity on unstructured mesh with a cell-centered lagrangian scheme[END_REF] in the framework of hyperelasticity. The corner matrix defined by ( 75) is symmetric definite positive and corresponds to the extension of the corner matrix initially introduced in the case of Lagrangian hydrodynamics [START_REF] Maire | A unified sub-cell force-based discretization for cell-centered lagrangian hydrodynamics on polygonal grids[END_REF][START_REF] Maire | A high-order one-step sub-cell force-based discretization for cell-centered lagrangian hydrodynamics on polygonal grids[END_REF][START_REF] Maire | Contribution to the numerical modeling of inertial confinement fusion[END_REF].

We conclude this paragraph by expressing the sub-cell force in terms of nodal stresses. Substituting the definition (75) of the corner matrix into the sub-cell force and using the geometric identity l pc n pc = l - pc n - pc + l + pc n + pc leads to the following expression of the sub-cell force

F pc = l - pc T - pc n - pc + l + pc T + pc n + pc , (77)
where T - pc n - pc and T + pc n + pc are the two nodal traction forces defined at the corner pc by

(T - pc -T c )n - pc = ρ c a L,c (n - pc ⊗ n - pc ) + a T,c Id -(n - pc ⊗ n - pc ) (V p -V c ), (T + pc -T c )n + pc = ρ c a L,c (n + pc ⊗ n + pc ) + a T,c Id -(n + pc ⊗ n + pc ) (V p -V c
). In the case of a hydrodynamic flow, the stress tensor collapses to T = -P I d , the shear modulus is equal to zero and thus the wave speeds become a L,c = a c , a T,c = 0. Thus, the nodal traction forces transform into pressure forces P - pc n - pc and P + pc n + pc respectively defined by

(P - pc -P c )n - pc = -ρ c a c (n - pc ⊗ n - pc )(V p -V c ), (P + pc -P c )n + pc = -ρ c a c (n + pc ⊗ n + pc )(V p -V c
). These last equations correspond to the half Riemann invariants written in the direction of the outward unit normals n - pc and n + pc . We conclude this section by pointing out that the choice of the corner matrix described above is certainly not the optimal one from the point of view of numerical dissipation. Another possible choice simply consists in re-using the corner matrix that has been initially developed in the context of pure hydrodynamic flow. Let us also note that this choice which is very similar to what is done in the framework of staggered discretization. Namely, for most of the staggered discretizations utilized in classical hydrocodes, the artificial viscosity, employed to ensure the dissipation of kinetic energy into internal energy through shock waves, remains the same regardless the material modeling, either hydrodynamic or elastic-plastic.

High-order extension based on the acoustic GRP method

We aim at deriving the high-order extension of the compatible cell-centered discretization presented in the previous section. To construct a high-order extension, many methods are available.

Inria

The most obvious one, consists in defining a monotonic piecewise linear reconstruction of the pressure and the velocity by means of a slope limiter. The nodal velocity is computed through the use of the node-centered solver wherein the input data are defined as the nodal extrapolated values of the pressure and the velocity. The time discretization is simply based on a twosteps Runge-Kutta procedure. Such a methodology has been successfully developed in [START_REF] Maire | A second-order cell-centered Lagrangian scheme for twodimensional compressible flow problems[END_REF][START_REF] Maire | Multi-scale Godunov-type method for cell-centered discrete Lagrangian hydrodynamics[END_REF]. However, this approach is rather expensive since it needs a two-step integration in time. This point becomes particularly crucial when coupling the hydrodynamic scheme with more complex physics. In particular, we want to avoid as much as possible the evaluation of the thermodynamic coefficients, that is we want to minimize the equation of state call. For this reason, we prefer to use a one-step time integrator based on the GRP (Generalized Riemann problem) method of Ben-Artzi and Falcovitz [START_REF] Ben-Artzi | A second-order Godunov-type scheme for compressible fluid dynamics[END_REF][START_REF] Ben-Artzi | An upwind second-order scheme for compressible duct flows[END_REF][START_REF] Ben-Artzi | Application of the "Generalized Riemann Problem" method to 1-D compressible flows[END_REF][START_REF] Ben-Artzi | Generalized Riemann problems in Computational Fluids Dynamics[END_REF][START_REF] Li | Remark on the generalized Riemann problem method for compressible fluid flows[END_REF]. This method consists in solving the higher-order Riemann problem with piecewise linear polynomials, whereby the approximate solution is given as a time power series expansion right at the interface, thus providing a numerical flux for a high-order Godunov-type method. Here, we are using the acoustic version of the GRP method. This approximation provides a framework in which the solution of the GRP is simple to compute and easy to handle. In the case of one-dimensional Lagrangian hydrodynamics, this method has been completely derived in the monograph [START_REF] Ben-Artzi | Generalized Riemann problems in Computational Fluids Dynamics[END_REF]. In what follows, we present the non-trivial extension of the acoustic GRP methodology to our two-dimensional Lagrangian scheme for elastic plastic flows.

Second-order time discretization

Let us introduce the time discretization of the semi-discrete system (54). All the physical and geometric variables are assumed to be known at the beginning of time step t n and we denote them using the superscript n. Their updated values at time t n+1 = t n + ∆t, where ∆t is the current time step, are obtained by means of the following system

m c 1 ρ n+1 c - 1 ρ n c - p∈P(c) t n+1 t n (l pc n pc )(t) • V p (t) dt = 0, (78a)
m c V n+1 c -V n c - p∈P(c) t n+1 t n F pc (t) dt = 0, (78b)
m c E n+1 c -E n c - p∈P(c) t n+1 t n F pc (t) • V p (t) dt = 0. (78c)
Here, the sub-cell force F pc is a time-dependent function which writes

F pc (t) = T c (t)(l pc n pc )(t) + M pc [V p (t) -V c (t)],
where the corner matrix, M pc is assumed to be positive definite. The motion of the grid is governed by the discrete trajectory equation

x n+1 p = x n p + t n+1 t n V p (t) dt, x 0 p = X p . (79)
We want to perform a one-step time integration of the numerical fluxes, which is at least secondorder accurate. To this end, we make use of the following Taylor expansions of the nodal velocity and the sub-cell force

∀t ∈ [t n , t n+1 ],        V p (t) = V p (t n ) + (t -t n ) d dt V p (t n ) + O((t -t n ) 2 ), F pc (t) = F pc (t n ) + (t -t n ) d dt F pc (t n ) + O((t -t n ) 2 ).
Employing these Taylor expansions, integrals over the time interval [t n , t n+1 ] of the nodal velocity and the sub-cell force are approximated by

t n+1 t n V p (t) dt = ∆t V n p + ∆t 2 d dt V p n , t n+1 t n F pc (t) dt = ∆t F n pc + ∆t 2 d dt F pc n ,
where d dt V p n and d dt F pc n denote the time derivative of the nodal velocity and the sub-cell force evaluated at time t n . Let us introduce the time centered values

V n+ 1 2 p = V n p + ∆t 2 d dt V p n , F n+ 1 2 pc = F n pc + ∆t 2 d dt F pc n .
The first equation of the above system shows that the second-order time discretization of the trajectory equation writes

x n+1 p = x n p + ∆tV n+ 1 2 p , x 0 p = X p .
Now, being given the time-centered value of the nodal velocity, observe that the position vector of point p can be parameterized over the time interval [t n , t n+1 ] as

x p (t) = x n p + (t -t n )V n+ 1 2 p , ∀t ∈ [t n , t n+1 ].
Let us also recall that the corner vector can be expressed in terms of the position vector of the vertex of cell c as

(l pc n pc )(t) = 1 2 x p + (t) -x p -(t) × e z ,
where x p + and x p -are the position vectors of the next and previous points with respect to point p in the counterclockwise ordered list of points of cell c. This shows that the corner vector, (l pc n pc )(t), has a linear dependency with respect to time. Thus, the time-centered approximation of the corner vector is computed exactly as being

(l pc n pc ) n+ 1 2 = 1 ∆t t n+1 t n (l pc n pc )(t) dt = (l pc n pc )(t n + ∆t 2 ). (80)
Therefore, the sub-cell contribution to the volume flux is computed as

t n+1 t n (l pc n pc )(t) • V p (t) dt = ∆t(l pc n pc ) n+ 1 2 • V n+ 1 2 p .

Inria

Gathering the previous results yields the following formal second-order time discretization of system (78)

m c 1 ρ n+1 c - 1 ρ n c -∆t p∈P(c) (l pc n pc ) n+ 1 2 • V n+ 1 2 p = 0, (81a)
m c V n+1 c -V n c -∆t p∈P(c) F n+ 1 2 pc = 0, (81b)
m c E n+1 c -E n c -∆t p∈P(c) F n+ 1 2 pc • V n+ 1 2 p = 0. ( 81c
)
The grid displacement is computed using the discrete trajectory equation

x n+1 p = x n p + ∆tV n+ 1 2 p , x 0 p = X p . (82)
Here, the time-centered variables nodal fluxes are given by

(l pc n pc ) n+ 1 2 = (l pc n pc )(t n + ∆t 2 ), (83a)
V

n+ 1 2 p = V n p + ∆t 2 d dt V p n , (83b)
F n+ 1 2 pc = F n pc + ∆t 2 d dt F pc n . (83c)
Comment 6. We show in [START_REF] Maire | Contribution to the numerical modeling of inertial confinement fusion[END_REF] that the previous time discretization fulfills the discrete geometric conservation law [START_REF] Donea | Encyclopedia of Computational Mechanics[END_REF]. This means that the zone volume that is computed directly from the coordinates of its vertices is rigorously equal to the zone volume that is deduced from solving the discrete volume conservation equation (81a). This result is of first importance and shows the consistency of our cell-centered discretization. We point out that this property is not straightforward to achieve for usual two-dimensional staggered discretizations unless a special procedure is applied [START_REF] Loubère | Volume consistency in a staggered Lagrangian hydrodynamics scheme[END_REF].

To achieve the above second-order time discretization it remains to compute the time derivatives of both nodal velocity and sub-cell force. This will be the main task of a next paragraph invoking the conservation principle of total energy. We have also to present the time discretization of the incremental constitutive law. This task is investigated in the next paragraph.

Time discretization of the constitutive law

Here, we present the time discretization of the evolution equation that governs the deviatoric stress

d dt S c = 2µ c (D 0,c -D p c ) -(S c W c -W c S c ),
where D 0,c denotes the deviatoric part of the strain rate tensor and D p c is the plastic strain rate defined by the associated flow rule (34b), (34c) and [START_REF] Margolin | Numerical simulation of plasticity at high strain rate[END_REF]. The cell-averaged strain rate tensor, D c and spin tensor, W c are given respectively by the mimetic discretizations [START_REF] Wilkins | Methods in Computational Physics, volume 3, chapter Calculation of Elastic-Plastic Flow[END_REF] and [START_REF] Ya | Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena[END_REF]. We point out that being given the nodal velocity field, V p , the above equation is nothing but a differential equation for the cell-averaged deviatoric stress, S c . We start by studying the time discretization of the elastic part of the incremental constitutive law.

Time discretization of the elastic part

We study the time discretization of the following initial value problem over the time interval

[t n , t n+1 ] Sc = 2µ c D 0,c , S c (t n ) = S n c , ( 84
)
where S is the Jaumann derivative of S defined by

S = d dt S + SW -WS.
Let us recall that this corotational rate satisfies the property of frame-indifference, moreover for rigid body rotation, the Jaumann derivative of S is equal to zero, refer to Section 2.2.3. We also point out that the Jaumann derivative preserves the material derivative of the squared norm of the deviatoric stress

(S c : S c ) = d dt (S c : S c ). ( 85
)
This important property is due to the fact that S c :

(S c W c -W c S c ) = 0.
In what follows, we shall derive a time discretization that preserves this property, knowing that the magnitude of the deviatoric stress is a fundamental ingredient in the plastic flow rule, which should not be modified by spurious discretization errors.

Let us note that being given the nodal velocity field and the geometry at t n+ 1 2 = 1 2 (t n +t n+1 ), we easily define the time-centered discretization of the strain rate tensor, D

n+ 1 2 c
, and the spin tensor, W n+ 1 2 by means of ( 47) and [START_REF] Ya | Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena[END_REF]. The deviatoric strain rate tensor, D

n+ 1 2 0,c
, is also easily obtained from the strain rate tensor by subtracting to it its trace. Let us integrate the differential equation ( 84) assuming that µ c , D 0,c and W c are constant over the time step [t n , t n+1 ]. More precisely, we assume that the deviatoric strain rate and the spin tensor are defined by their time-centered values: D . Now, we are going to transform the differential equation (84) utilizing the interpretation of the Jaumann rate in terms of rotation that has been presented in Section 2.2.3. More precisely, we have shown that the Jaumann rate reads

S = Ω d dt (Ω t SΩ) Ω t ,
where the rotation, Ω, is the unique solution to the initial value problem

d dt Ω = WΩ, Ω(0) = Id.
Therefore, the semi-discrete constituve law (84) rewrites

d dt (Ω t c S c Ω c ) = 2µ c Ω t c D 0,c Ω c . (86)
Here, Ω c denotes the solution to the initial value problem

d dt Ω c = W n+ 1 2 c Ω c , ∀t ∈ [t n , t n+1 ], Ω(t n ) = Id.
As it is shown in Appendix A, the unique solution to this initial value problem is expressed in terms of the tensor exponential of the discrete spin tensor, W

n+ 1 2 c
,

Ω c (t) = exp (t -t n )W n+ 1 2 c , ∀ ∈ [t n , t n+1 ]. ( 87
)

Inria

The time integration of (86) over the time interval [t n , t n+1 ] yields

Ω t c (t n+1 )S n+1 Ω c (t n+1 ) -S n = 2µ c t n+1 t n Ω t c (t)D n+ 1 2 0,c Ω c (t) dt.
Finally, substituting the definition (87) of the rotation into the above equation leads to S n+1exp(∆tW

n+ 1 2 c
)S n exp(-∆tW

n+ 1 2 c ) = 2µ c ∆t 0 exp(tW n+ 1 2 c
)D

n+ 1 2 0,c exp(-tW n+ 1 2 c ) dt, ( 88
)
where ∆t = t n+1t n denotes the time step. Let us emphasize that (88) is the exact solution of the initial value problem (84) granted that the spin tensor and the deviatoric strain rate tensor are constant over the time step [t n , t n+1 ]. Finally, approximating the integral at the right-hand side by the classical midpoint quadrature yields the second-order accurate numerical scheme

S n+1 c -exp(∆tW n+ 1 2 c
)S n c exp(-∆tW

n+ 1 2 c ) = 2µ c ∆t exp( ∆t 2 W n+ 1 2 c
)D

n+ 1 2 0,c exp(- ∆t 2 W n+ 1 2 c
). ( 89)

From a practical point of view, the numerical evaluation of the tensor exponential function is performed by means of the second-order accurate Padé approximant. It has been introduced in Appendix A and we recall its expression hereafter

Q[A](t) = Id - t 2 A -1 Id + t 2 A , ( 90
)
where A is an arbitrary skew symmetric tensor. Let us recall that Q[A](t) is a rotation which approximates the tensor exponential of A up to second-order in the vicinity of t = 0, i.e., exp(tA) = Q[A](t) + O(t 3 A 3 ). With this notation, the numerical scheme (89) turns into

S n+1 c -Q[W n+ 1 2 c ](∆t)S n c Q t [W n+ 1 2 c ](∆t) = 2µ c ∆tQ[W n+ 1 2 c ]( ∆t 2 )D n+ 1 2 0,c Q t [W n+ 1 2 c ]( ∆t 2 ), ( 91
)
where Q is defined by (90). Therefore, we have constructed a second-order accurate time discretization of the constitutive law in the elastic regime. Moreover, we claim that this scheme preserves the magnitude of the deviatoric stress in the case of a rigid body rotation. In this particular case, it is clear that the deviatoric strain is equal to zero, hence the updated elastic deviatoric stress reads as

S n+1 c = Q[W n+ 1 2 c ](∆t)S n c Q t [W n+ 1 2 c
](∆t). This means that S n+1 c is nothing but the image of S n c by the rotation Q[W

n+ 1 2 c ](∆t), thus | S n+1 c |=| S n c |.
Finally, we point out that this numerical scheme belongs to the class of incrementally objective algorithms, that is, the discrete scheme preserves the crucial property of material objectivity, refer to [START_REF] De Souza Neto | Computational Methods for Plasticity: Theory and Applications[END_REF][START_REF] Simó | Computational Inelasticity[END_REF].

Time discretization of the plastic part

Let us remark that equation (91) corresponds to the time discretization of the incremental constitutive law after having frozen the plastic flow contribution, i.e., D p c . Bearing this in mind, we denote by S n+1,tr c , the value of the deviatoric stress at time t n+1 given by (91), in the sense that it is a predicted trial elastic state. This intermediate value requires to be corrected by means of a corrector step to take into account plasticity. To this end, we assess the value of the yield (91) is plastic and the trial state should be corrected in order to satisfy the plastic consistency condition, f (S n+1 ) = 0. To this end, we will use the radial return mapping algorithm originally proposed by Wilkins [START_REF] Wilkins | Methods in Computational Physics, volume 3, chapter Calculation of Elastic-Plastic Flow[END_REF]. First, let us recall that the evolution of the deviatoric stress tensor during the plastic phase is ruled by

function, f (S) =| S | -2 3 Y 0 . If f (S n+1
d dt S c = -2µ c D p c ,
Knowing from the flow rule (34c) that the plastic strain rate, D p c , is colinear to the unit normal, N p c = Sc |Sc| , that characterizes the plastic flow direction, we write it as D p c = Λ c N p , where Λ c is a real valued coefficient named the plastic multiplier. Time integration of the previous equation, assuming the plastic strain rate remains constant over the time step, leads to

S n+1 c -S n+1,tr c = -2µ c ∆t Λc N p,n+1 , (92)
where Λc =

| S n+1 c | N p,n+1 -| S n+1,tr c | N p,n+1,tr = -2µ c ∆t Λc N p,n+1 . ( 93
)
This clearly shows that the trial unit normal vector, N p,n+1,tr , which is known from the trial state, is in the direction of the actual unit normal vector, N p,n+1 , which is not known and should be found. This property shows that the plastic corrector step is performed in the radial direction to the yield surface and the trial unit vector remains unchanged during the plastic update. Thus, the trial unit normal vector defines the direction of the return to the actual yield surface defined by f (S n+1

S n+1,tr c -Q[W n+ 1 2 c ](∆t)S n c Q t [W n+ 1 2 c ](∆t) = 2µ c ∆tQ[W n+ 1 2 c ]( ∆t 2 )D n+ 1 2 0,c Q t [W n+ 1 2 c ]( ∆t 2 ),
where the rotation, Q, is given by

Q[A](t) = Id - t 2 A -1 Id + t 2 A ,
for any arbitrary skew symmetric tensor A. At this point it is interesting to mention that the method of radial return described above is an explicit first-order approximation to the actual plastic stress strain equation. This approximation is not always accurate particularly when the deviator strain rate is comparable in magnitude to the total dilatational strain rate as it has been noticed in [START_REF] Margolin | Numerical simulation of plasticity at high strain rate[END_REF]. In this reference, the authors propose a methodology to cure this flaw by constructing an exact solution to the elastic-plastic equation assuming that the strain rate tensor is constant over a time step. In the future we plan to investigate this approach and to adapt it to our numerical scheme.

Total energy and momentum conservation

Let us investigate total energy conservation over the entire domain at the discrete level. To this end, we transpose at the discrete time level, the reasoning developed in Section 4.1. Recalling that total energy over the whole grid at time t n is defined by E n = c m c E n c , its conservation over the time interval [t n , t n+1 ] writes as

E n+1 -E n = t n+1 t n ∂D (TV ) • n ds dt,
where the right-hand side expresses the time rate of pressure work on the boundary, ∂D, of the domain, D, occupied by the material. By definition of total energy, this last equation rewrites

c m c (E n+1 c -E n ) = t n+1 t n ∂D (TV ) • n ds dt. ( 95
)
The discretization of the right-hand side of the above equation follows the methodology initially introduced in Section 4.1. Therefore, substituting the specific total energy equation (81c) into (95) leads to

c p∈P(c) F n+ 1 2 pc • V n+ 1 2 p = p∈∂D F ,n+ 1 2 p • V n+ 1 2 p . ( 96
)
where

F ,n+ 1 2 p
is a time-centered evaluation of the prescribed boundary force acting onto point p defined by (65). This boundary force is expressed as

F ,n+ 1 2 p = F n, p + ∆t 2 d dt F p n . (97)
Here, F n, p denotes the value of the prescribed boundary force at time t n whereas d dt F p n represents its time derivative at time t n . Let us emphasize that these quantities are known since they correspond to the prescribed boundary conditions. Interchanging the order of summation in the left-hand side of (96) yields

p ( c∈C(p) F n+ 1 2 pc ) • V n+ 1 2 p = p∈∂D F ,n+ 1 2 p • V n+ 1 2 p , RR n°7975
where C(p) is the set of cells surrounding point p. Splitting the left-hand side of the above equation into two parts depending on the points location leads to

p∈D o ( c∈C(p) F n+ 1 2 pc ) • V n+ 1 2 p + p∈∂D ( c∈C(p) F n+ 1 2 pc ) • V n+ 1 2 p = p∈∂D F ,n+ 1 2 p • V n+ 1 2 p , ( 98
)
where D o is the interior of the domain D. The total energy conservation principle requires that the above equation must hold regardless the value of the nodal velocity. This amounts to state the following necessary and sufficient conditions to have total energy conservation

∀ p ∈ D o , c∈C(p) F n+ 1 2 pc = 0, (99a)
∀ p ∈ ∂D, c∈C(p) F n+ 1 2 pc = F ,n+ 1 2 p . ( 99b
)
Now, substituting the expression of the time-centered sub-cell force, F

n+ 1 2 pc = F n pc + ∆t 2 d dt F pc n
and invoking the fact that (99a) and (99b) must be satisfied regardless of the time step value leads to the following necessary and sufficient conditions:

• For the sub-cell force

∀ p ∈ D o , c∈C (p)
F n pc = 0, (100a)
∀ p ∈ ∂D, c∈C (p)
F n pc = F ,n p . (100b)
• For the time derivative of the sub-cell force

∀ p ∈ D o , c∈C (p)
d dt F pc n = 0, (101a)
∀ p ∈ ∂D, c∈C (p)
d dt F pc n = d dt F p ,n . (101b)
Condition (100) corresponds to the sub-cell forces balance at point p, whereas condition (101) represents the balance of the time derivative of the sub-cell forces acting onto point p. We claim that these conditions also provide momentum conservation over the entire domain, the proof is left to the reader. Recalling that the sub-cell force reads as

F pc (t) = T c (t)(l pc n pc )(t) + M pc [V p (t) -V c (t)]
, the two conditions (100) and (101) allow us to construct respectively a node-centered solver for computing the nodal velocity at time t n and a node-centered solver for computing the nodal acceleration at the same time.

Node-centered solver for the grid velocity at time t n

We present the node-centered solver that allows to compute the nodal velocity at time t n . For the high-order extension, this solver relies on the condition (100) wherein F n pc is the nodal extrapolation of the sub-cell force given by

F n pc = T c (x n p ) [(l pc n pc ) n ] + M n pc V n p -V c (x n p ) . (102)

Inria

Here, T c = T c (x) and V c = V c (x) denote the piecewise monotonic linear reconstruction of the Cauchy stress tensor and the velocity field within the cell ω c . More precisely, T c = -P c Id + S c , where P c and S c are piecewise monotonic linear reconstructions of the pressure and the deviatoric stress within cell c. These piecewise linear reconstruction will be described in Section 5.7. Substituting the expression of the nodal extrapolation of the sub-cell force (102) into the conservation condition (100) leads to the 2 × 2 system that solves the nodal velocity

V n p ∀ p ∈ D o , M n p V n p = c∈C(p) -T c (x n p ) [(l pc n pc ) n ] + M n pc V c (x n p ) , (103a)
∀ p ∈ ∂D, M n p V n p = c∈C(p) -T c (x n p ) [(l pc n pc ) n ] + M n pc V c (x n p ) + F ,n p , (103b)
where M n p denotes the 2 × 2 node-centered matrix defined as

M n p = c∈C(p) M n pc . (104)
The nodal velocity determination has been divided into two cases depending on the node location. Noticing that M n p is positive definite, V n p is uniquely defined by solving systems (103a) and (103b).

Node-centered solver for computing the acceleration of the grid velocity at time t n

In this paragraph we aim at describing the node-centered solver which allows to compute the acceleration of the grid velocity at time t n . To this end, we shall use the conservation condition (101) that has been derived in the last paragraph. Knowing that the sub-cell force reads

F pc (t) = T c (t)(l pc n pc )(t) + M pc [V p (t) -V c (t)],
we time differentiate it using the chain rule

d dt F pc = d dt T c (l pc n pc ) + M pc d dt V p - d dt V c , (105)
where d dt T c and d dt V c denotes the cell-centered time derivatives of the Cauchy stress tensor and the velocity field. We shall describe later how to compute these quantities. Here, we have only taken into account the time variation of the Cauchy stress tensor, the nodal velocity and the cell-centered velocity, whereas we have neglected the time variation of the corner normal and the corner matrix. We have checked that this assumption did not decrease the overall accuracy of the numerical scheme. On the other hand, it considerably simplifies the construction of the nodecentered solver for the acceleration by rendering its structure similar to that of the node-centered solver for the velocity.

Substituting the expression of the time derivative of the sub-cell force (105) evaluated at time t n into the conservation condition (101) leads to the 2 × 2 system satisfied by the nodal

acceleration d dt V p n ∀ p ∈ D o , M n p d dt V p n = c∈C(p) - d dt T n c [(l pc n pc ) n ] + M n pc d dt V n c , (106a)
∀ p ∈ ∂D, M n p d dt V p n = c∈C(p) - d dt T n c [(l pc n pc ) n ] + M n pc d dt V n c + d dt F p ,n . (106b)
where (∇ • V ) c denotes the cell-centered velocity divergence obtained by means of the piecewise monotonic linear reconstruction of the velocity field within cell c. Let us note that the cellaveraged value of the strain rate tensor is given by the mimetic discretization [START_REF] Wilkins | Methods in Computational Physics, volume 3, chapter Calculation of Elastic-Plastic Flow[END_REF]. Last, we compute the cell-centered time derivative of the deviatoric stress by means of the incremental constitutive law which reads

d dt S c = 2µ c (D 0,c -D p c ) -(S c W c -W c S c ),
where D 0,c denotes the deviatoric part of the strain rate tensor and D p c is the plastic strain rate defined by the associated flow rule (34b), (34c) and [START_REF] Margolin | Numerical simulation of plasticity at high strain rate[END_REF]. The cell-averaged spin tensor, W c is given by the mimetic discretization (48).

Piecewise monotonic linear reconstruction

Here, we need to perform piecewise linear monotonic reconstruction of the pressure, the velocity and the deviatoric stress tensor on general two-dimensional unstructured grids. To this end, we employ a classical least squares procedure [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF], followed by a slope limitation procedure to enforce monotonicity [START_REF] Desideri | Compressible flow using unstructured grids[END_REF][START_REF] Dukowicz | Accurate conservative remapping (rezoning) for arbitrary lagrangian-eulerian computations[END_REF].

Piecewise linear reconstruction by means of a least squares approach

Let W ≡ W (x) denotes a fluid variable (pressure, velocity or deviatoric stress components), which has a piecewise linear representation within cell c defined by

W c (x) = W c + (∇W ) c • (x -x c ) , (111)
where W c is the mean value of W in cell c and (∇W ) c is the gradient of W that we are looking for. To ensure a conservative reconstruction, we define x c as being the cell centroid, that is x c = 1 vc ωc x dv. Then, the gradient in (111) is computed by imposing that

W c (x d ) = W d for d ∈ C(c),
where C(c) is the set of the neighboring cells of cell c. This problem is generally over-determined and thus this gradient is obtain by using a least squares procedure. Therefore, it is the unique solution of the following minimization problem

(∇W ) c = argmin G∈R 2 d∈C(c) [W d -W c -G • (x d -x c )] 2 .
A straightforward computation shows that this solution writes

(∇W ) c = M -1 c d∈C(c) (W d -W c ) (x d -x c ) , (112)
where M c is the 2 × 2 matrix which reads

M c = d∈C(c) (x d -x c ) ⊗ (x d -x c ) .
We notice that M c is symmetric positive definite and thus always invertible. The main feature of this least squares procedure is that it is valid for any type of unstructured mesh and moreover it preserves the linear fields.

Using this approach, we can construct the following piecewise linear reconstruction of the pressure, the deviatoric stress and the velocity field within cell c

P c (x) = P c + (∇P ) c • (x -x c ), S c (x) = S c + (∇S) c (x -x c ), V c (x) = V c + (∇V ) c (x -x c ).
Here, (∇P ) c is the cell-centered pressure gradient, which is a vector, (∇S) c is the cell-centered deviatoric stress gradient, which is a third order tensor, and (∇V ) c is the cell-centered velocity gradient, which is a second order tensor. Let us point out that for vectors and tensors the least squares procedure is employed componentwise.

Monotonicity

To preserve monotonicity, we limit the value that the gradient is allowed to take. For each cell, we introduce the slope limiter φ c ∈ [0, 1] and the limited reconstructed field as

W lim c (x) = W c + φ c (∇W ) c • (x -x c ) , (113)
where (∇W ) c denotes the approximate gradient given by (112). The coefficient φ c is determined by enforcing the following local monotonicity criterion

W min c ≤ W lim c (x) ≤ W max c , ∀x ∈ c, (114)
where we have set

W min c = min min d∈C(c) {W d }, W c and W max c = max max d∈C(c) {W d }, W c .
Since the reconstructed field is linear we note that it is sufficient to enforce the following conditions at any point p ∈ P(c)

W min c ≤ W lim c (x p ) ≤ W max c , (115)
so that the quantity W in the cell c does not lie outside the range of the average quantities in the neighboring cells. Thanks to this formula, we can define the slope limiter as

φ c = min p∈P(c) φ c,p knowing that φ c,p =                          ν W max c -W c W c (x p ) -W c if W c (x p ) -W c > 0, ν W min c -W c W c (x p ) -W c if W c (x p ) -W c < 0, 1 if W c (x p ) -W c = 0.
Here, ν denotes a real valued function characterizing the limiter. By setting ν(x) = min(1, x) we recover a limiter which consists of the multi-dimensional extension [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF] of the van Leer's classical method. We can also define a smoother -in the sense that it is more differentiable-limiter by

Inria setting ν(x) = x 2 +2x
x 2 +x+2 . This limiter has been introduced by Venkatakrishnan [START_REF] Venkatakrishnan | Convergence to steady state solutions of the Euler equations on unstructured grids with limiters[END_REF] in order to improve the convergence towards steady solutions for the Euler equations.

We pursue by giving indications about the limitation procedure applied to the velocity field and the deviatoric stress. We point out that the previous limiting algorithm has been developed for scalar variables. For vectors and tensors, limiting is usually applied separately to each component. However, such a procedure is frame dependent and thus leads to rotational symmetry distortion. Namely, component limiters do not preserve symmetry since a rotation of the coordinate axis produces different results. This drawback is crucial in the framework of Lagrangian hydrodynamics since we are dealing with moving mesh discretizations which are particularly sensitive about symmetry loss. To correct this flaw, we have to construct a limiting procedure which is frame invariant for vectors and tensors. We note that for vectors it is possible to use the Vector Image Polygon (VIP) methodology derived by Luttwak and Falcovitz in [START_REF] Luttwak | Slope limiting for vectors: a novel vector limiting algorithm[END_REF][START_REF] Luttwak | Slope limiting for vectors: A novel vector limiting algorithm[END_REF]. Recently, a frame invariant limiting procedure for symmetric tensors, which relies on the principal invariants, has been proposed in [START_REF] Sambasivan | Exploration of new limiter schemes for stress tensors in Lagrangian and ALE hydrocodes[END_REF].

Here, we shall employ for vectors the limiting procedure that has been described in [START_REF] Maire | A high-order one-step sub-cell force-based discretization for cell-centered lagrangian hydrodynamics on polygonal grids[END_REF]. This method consists in defining in each cell an intrinsic orthonormal basis which relies for instance on the natural direction of the flow computed from the eigenvectors of the strain rate tensor. Then, the velocity field is projected onto this local basis and the above scalar limiting procedure is applied for each projected components of the velocity field. We obtain one slope limiter for each local component of the velocity field. Finally, we transform this pair of limiters back into the Cartesian coordinates and compute the limited velocity gradient. We claim that this limitation procedure is frame invariant and thus preserves rotational symmetry. Regarding the tensors, we apply a straightforward extension of the above method by computing the slope limiters corresponding to the local components of the tensor projected onto the aforementioned local orthonormal basis.

Finally, let us note that the above limitation procedure does not ensure that the extrapolated deviatoric stress, S c (x p ), fulfills the yield criterion. To correct this flaw, we apply the radial return algorithm in the following manner. Let us denote S tr c (x p ) the nodal extrapolated value of the deviatoric stress, resulting from the piecewise monotonic linear reconstruction. The superscript tr refers to as a trial state that requires to be corrected to take into account plasticity. To this end, we assess the value of its corresponding yield function, i.e., f (S) =| S | -2 3 Y 0 . If f [ S tr c (x p )] ≤ 0, the trial nodal extrapolated deviatoric stress remains below the yield strength and thus is not corrected. On the other hand, if f [ S tr c (x p )] > 0, the trial nodal extrapolated deviatoric stress is plastic and it must be corrected to satisfy the plastic consistency condition f [ S c (x p )] = 0. To this end, we define the corrected nodal extrapolated value as

S c (x p ) = 2 3 Y 0 S tr c (x p ) | S tr c (x p ) | . ( 116
)
6 Numerical results

We present some numerical results to assess the accuracy and the robustness of our numerical scheme. In all these test cases the materials under consideration are characterized by a Mie-Grüneisen equation of state. For sake of completeness, we recall hereafter the main feature related to this type of equation of state.

Mie-Grüneisen equation of state

Being given the specific internal energy, ε, and the specific volume, V = 1 ρ , the Mie-Grüneisen equation of state expresses the pressure, P , as follows

P (V, ε) = a 2 0 (V 0 -V ) [V 0 -s(V 0 -V )] 2 + Γ(V ) V ε - 1 2 a 0 (V 0 -V ) V 0 -s(V 0 -V ) 2 .
(117)

Here, V 0 1 ρ0 , denotes the specific volume of the unstressed material, a 0 and s are coefficient that relate the shock speed U s and the particle velocity u p which is related to the shock velocity by the empirically obtained linear relation ship U s = a 0 + su p , refer to [START_REF] Menikoff | Equations of state and fluid dynamics[END_REF]. The Grüneisen parameter, Γ, is given by

Γ = V ∂P ∂ε V = Γ 0 V V 0 .
For practical numerical computations, we recast this equation of state under the form

P (ρ, ε) = ρ 0 a 2 0 f (η) + ρ 0 Γ 0 ε, f (η) = (η -1)[η -1 2 Γ 0 (η -1)] [η -s(η -1)] 2
, where η = ρ ρ 0 .

For a generic equation of state written under the form, P = P (ρ, ε) the isentropic sound speed, a, is defined by

a 2 = ∂P ∂ρ ε + P ρ 2 ∂P ∂ε ρ .
Employing the definition of the Mie-Grüneisen equation of state leads to the practical formula

a 2 = a 2 0 f (η) + Γ 0 ρ 0 P η 2 , f (η) = η + (s -Γ 0 )(η -1) [η -s(η -1)] 3 .

Elastic-plastic piston-like problem

In this paragraph, we assess the accuracy of our numerical scheme against an elastic-plastic piston-like problem. This is a one-dimensional flow characterized by an analytical solution.

Before showing the numerical results, let us recall briefly the main ingredients used to derive the analytical solution. For a one-dimensional flow aligned with the x direction, the elastic-plastic system (1) becomes

ρ d dt ( 1 ρ ) - ∂u ∂x = 0, (118a)
ρ du dt - ∂T ∂x = 0, (118b)
ρ dE dt - ∂ ∂x (T u) = 0. (118c)
Here, the first component of the Cauchy stress tensor reads T = -P + S xx , where P is the pressure given by the Mie-Grüneisen equation of state and S xx is the first component of the Inria deviatoric stress. The incremental constitutive law in the elastic regime writes

d dt S xx = 4 3 µ ∂u ∂x , (119a)
d dt S yy = - 2 3 µ ∂u ∂x , (119b)
d dt S zz = - 2 3 µ ∂u ∂x . (119c)
Let us note that for a one dimensional flow, the material derivative coincides with the Jaumann derivative. The plasticity is taken into account by means of the yield criterion f =| S | -2 3 Y 0 ≤ 0. Assuming a deviatoric stress free initial state, the solution of (119) satisfies S yy = S zz = - Knowing that S xx (x, 0) = 0, the above equation is straightforwardly solved and we get

S xx (x, t) = - 4 3 µ log ρ ρ 0 , (120)
where ρ 0 is the initial density. Knowing that the elastic regime is bounded by the yield criterion, | S xx | -2 3 Y 0 ≤ 0, leads to the following inequality for the density exp -

Y 0 2µ ≤ ρ ρ 0 ≤ exp Y 0 2µ . ( 121
)
At this point, it is natural to introduce

ρ - e = ρ 0 exp - Y 0 2µ , ρ + e = ρ 0 exp Y 0 2µ .
Here, ρ - e and ρ + e are respectively the minimum and the maximum density related to an elastic rarefaction and an elastic compression. With this notation, we are able to express the first deviatoric stress component in terms of the density

S xx (ρ) =                -2 3 Y 0 if ρ ≥ ρ + e , -4 3 µ log ρ ρ0 if ρ - e ≤ ρ ≤ ρ + e , 2 3 Y 0 if ρ ≤ ρ - e (122)
Using this analytical expression of the first component of the deviatoric stress, the first component of the Cauchy stress tensor reads T = -P (ρ, ε)+S xx (ρ). Thus, the one-dimensional system (118) has a complete thermodynamic closure. Therefore, it is possible to solve the Riemann problem corresponding to this system. For instance, to investigate the admissible shock waves related to this system, it suffices to write the classical Rankine-Hugoniot equations. For a piston-like problem, being given the piston velocity, it is straightforward to solve the Rankine-Hugoniot equations in the elastic regime, i.e., ρ ∈ [ρ 0 , ρ + e ], to get all the physical variables behind the elastic shock wave. In the elastic-plastic regime, the analysis is more subtle due to the fact that the Hugoniot curve is characterized by a slope discontinuity at the density ρ = ρ + e . This slope discontinuity leads to a classical two-wave structure characterized by a leading elastic shock wave, named elastic precursor, followed by a plastic shock wave, refer to [START_REF] Ya | Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena[END_REF]. For a detailed analysis of this type of analytical solution, the interested reader might refer to [START_REF] Udaykumar | An Eulerian method for computation of multimaterial impact with ENO shock-capturing and sharp interfaces[END_REF]. Now, let us describe the data which characterize our elastic-plastic piston-like problem. The material under consideration is made of copper and its equation of state is given by the Mie-Grüneisen model (117) characterized by the following parameters: ρ 0 = 8930 kg/m 3 , a 0 = 3940 m/s, Γ 0 = 2 and s = 1.49. The incremental elastic-plastic constitutive law is characterized by the shear modulus, µ = 45 10 9 Pa, and the yield strength, Y 0 = 90 10 6 Pa. The computational domain is defined by (x, y) ∈ [0, 1] × [0, 0.1] and paved with (n x × n y ) cells, where n y = 2 and n x = 100, 200, 400. The initial data are prescribed by (P 0 , V 0 ) = (10 5 Pa, 0e x ), where P 0 and V 0 denote respectively the initial pressure and the initial velocity. At the left boundary, we impose a velocity boundary condition given by V piston = 20 m/s, whereas at the right boundary we prescribe a wall boundary condition. This test case is run until the stopping time t = 150 10 -6 s. We have respectively displayed in Figure 6 and in Figure 7 the density, the pressure, the velocity and the deviatoric stress. We have performed a convergence analysis using a sequence of three grids characterized respectively by ∆x = 10 -2 , 5 10 -3 , 2.5 10 -3 . The results are evidently converging to the analytical solution, in terms of value and the shock locations as well for the leading elastic shock as for the plastic shock wave.

Wilkins problem: flying aluminium plate striking a target plate

This one-dimensional test case, which has been initially proposed by Wilkins in [START_REF] Wilkins | Methods in Computational Physics, volume 3, chapter Calculation of Elastic-Plastic Flow[END_REF], consists in computing the flow resulting from a flying aluminium plate that strikes a target aluminium plate. Let us remark that this problem has been also studied in the following references [START_REF] Favrie | Diffuse interface model for compressible fluid-compressible elastic-plastic solid interaction[END_REF][START_REF] Miller | A high-order Eulerian Godunov method for elastic-plastic flow in solids[END_REF][START_REF] Hill | An Eulerian hybrid WENO centered-difference solver for elastic-plastic solids[END_REF]] by means of Eulerian numerical schemes. The equation of state is given by the Mie-Grüneisen model (117) characterized by the following parameters: ρ 0 = 2785 kg/m 3 , a 0 = 5328 m/s, Γ 0 = 2 and s = 1.338. The incremental elastic-plastic constitutive law is characterized by the shear modulus, µ = 27.6 10 9 Pa, and the yield strength, Y 0 = 300 10 6 Pa. The computational domain is defined by (x, y) ∈ [0, 50 10 -3 ] × [0, 5 10 -3 ] and paved with (500 × 2) cells, which corresponds to x-spacing used in [START_REF] Wilkins | Methods in Computational Physics, volume 3, chapter Calculation of Elastic-Plastic Flow[END_REF]. At time t = 0, the flying plate is located in the interval x ∈ [0, 5 10 -3 ], that is, it has a width of 5 10 -3 m. The initial data are prescribed by (P 0 , V 0 ) = (10 -6 Pa, 800e x ), for x ∈ [0, 5 10 -3 ], (P 0 , V 0 ) = (10 -6 Pa, 0), for x ∈ [5 10 -3 , 50 10 -3 ],

where P 0 and V 0 denote respectively the initial pressure and the initial velocity. We observe that the flying plate is characterized by a initial velocity of 800 m/s. We impose a free boundary condition at the left side and a wall boundary condition at the right side. This test case is run until the stopping time t = 5 10 -6 s. We have respectively displayed in Figure 8 and in Figure 9 the density, the pressure, the velocity and the deviatoric stress. At initial time, left-and right-facing shocks propagate from the interface located at x = 5 10 -3 . Due to the plasticity, the right-facing shock leads to two right-facing shocks: an elastic precursor followed by a plastic shock. These right-facing waves are clearly visible in Figures 8(a), 8(b) and 9(a). Similar left-facing shock waves are also created. However, when, these left-facing shocks reach the free surface, they are reflected into right-facing rarefaction waves which begin to overtake the right-facing shocks. We observe in Figure 8 that the elastic rarefaction wave is followed by a plastic rarefaction wave. Let us note that the numerical results obtained with the present numerical are in good agreement with the results obtained in [START_REF] Favrie | Diffuse interface model for compressible fluid-compressible elastic-plastic solid interaction[END_REF][START_REF] Miller | A high-order Eulerian Godunov method for elastic-plastic flow in solids[END_REF].

Collapse of a thick-walled cylindrical beryllium shell

This test problem has been initially proposed by Howell and Ball in [START_REF] Howell | A free-Lagrange augmented Godunov method for the simulation of elastic-plastic solids[END_REF]. It consists in computing the collapse of a cylindrical beryllium shell undergoing a initial radial velocity field directed towards its center. The initial kinetic energy of the shell is entirely converted into internal energy through plastic dissipation. The final state of the shell at the end of the collapse is characterized by its inner and outer stopping radius which can be expressed analytically in terms of the initial conditions. The equation of state is given by the Mie-Grüneisen model (117) characterized by the following parameters: ρ 0 = 1845 kg/m 3 , a 0 = 12870 m/s, Γ 0 = 2 and s = 1.124. The incremental elastic-plastic constitutive law is characterized by the shear modulus, µ = 151.9 10 9 Pa, and the yield strength, Y 0 = 330 10 6 Pa. The computational domain is defined in polar coordinates by (r, φ)

∈ [R i , R o ] × [0, π 2 ]
, where R i = 80 10 -3 m is the initial inner radius and and R o = 100 10 -3 m is the initial outer radius. To perform a convergence Inria analysis we have paved the initial computational domain we a sequence of three grids which are characterized respectively by (10 × 8), (20 × 16) and (40 × 32) uniformly spaced cells. The initial data are prescribed by (P 0 , V 0 ) = (10 -6 Pa, -V 0 Ri r e r ) where V 0 denotes the magnitude of the initial radial velocity field, and r = x 2 + y 2 is the radius of a generic point whose Cartesian coordinates are (x, y). We prescribe free boundary conditions, i.e., P = P 0 at the inner and outer boundary of the shell whereas we impose symmetry boundary conditions at the edges characterized by φ = 0 and φ = Π 2 . Howell and Ball have derived an analytical solution of this problem by writing the onedimensional Lagrangian governing equations for cylindrical symmetry. Then, assuming that the shell material is incompressible, they compute the stresses to produce yielding of the cylinder. Finally, writing the conservation of energy, that is, equating the plastic work done during the collapse to the total initial kinetic energy, leads to an analytical relationship between the magnitude of the initial velocity and the inner and outer stopping radius. Let us point out that this relationship depends solely on the geometrical configuration, the yield strength, Y 0 , and the density. For this reason, although we are not using the same equation of state than the one described in [START_REF] Howell | A free-Lagrange augmented Godunov method for the simulation of elastic-plastic solids[END_REF], we can utilize directly their analytical solution. For V 0 = 417.1 m/s, Howell and Ball have computed the following values for the inner and outer stopping radius: R i = 50 10 -3 m and R o = 78.1 10 -3 m.

To assess the numerical convergence of our scheme, the time history of the inner and outer radius of the shell are plotted versus their analytical values in Figure 10(a) for the three aforementioned computational grids. We obtain a rather good agreement between the theoretical values and the numerical ones. Moreover, we clearly observe the numerical convergence towards the analytical values. We have displayed in Figure 10(b) the time history of the kinetic energy and the internal energy of the collapsing shell for the three computational grids used. We observe that the kinetic energy is converted into internal energy through plastic dissipation while the summation of kinetic and internal energies remains evidently equal to a constant total energy. Finally, we have displayed the finest computational grid at initial and final time in Figure 11(a) and Figure 11(b). The final time being the time at which 0.99 of the initial kinetic energy has been converted into internal energy. We observe that our numerical scheme ensures a rather good symmetry preservation. We remark that the shell has thickened during its collapse due to the radial convergence. We have also displayed in Figure 12(a), the color map of the normalized plasticity threshold given by 3 2 √ (S:S) Y 0 . The density color map is plotted in Figure 12(b). These color maps are plotted at time t = 130 10 -6 s, which is slightly greater than the final time. Accordingly, we observe that the shell is coming back to the elastic regime.

Two-dimensional projectile impact

This problem which consists of the impact of a two-dimensional aluminium projectile against a rigid wall, has been initially proposed and simulated in [START_REF] Kluth | Discretization of hyperelasticity on unstructured mesh with a cell-centered lagrangian scheme[END_REF] by means of an hyperelastic model. It is the counterpart in planar geometry of the famous Taylor impact problem which employs cylindrical rod. Here, we use the same data than in [START_REF] Kluth | Discretization of hyperelasticity on unstructured mesh with a cell-centered lagrangian scheme[END_REF], recalling that the models used are quite different. The initial projectile is a rectangular plate of length L = 5 m and thickness H = 1 m. The material under consideration is made of aluminium and is described by means of the Mie Grüneisen equation of state. The values of the parameters are the same than those described in Section 6.3. The initial computational domain is [x, y] ∈ [0, 5] × [-0.5, 0.5]. The initial velocity is given by V 0 = -150e x m/s. The boundary conditions are free traction boundary conditions for all the boundary of the domain excepted the left boundary which is a wall boundary condition. We point out that this problem has no analytical solution neither Inria experimental results. However, it is a good test to assess the robustness and the accuracy of our numerical method. To perform a convergence analysis, we run this problem using a sequence of three computational grids characterized by 50 × 10, 100 × 20 and 200 × 40 uniformly spaced cells. The stopping time is t = 5 10 -3 s. During the impact, the kinetic energy is entirely converted into internal energy through plastic dissipation as it can be observed in Figure 14(a). We have plotted the time history of the length of the projectile for the three computational grids. It turns out that this length tends towards a limit which is approximately equal to L final = 4.62 m. We have also displayed in Figure 13(a) the computational grid at the final time and in Figure 13(b) the color map of the plasticity threshold which is defined by 3 2 √ (S:S) Y 0 . To assess the ability of our numerical scheme to cope with unstructured grids, we have run this problem employing a general polygonal grid. This computational grid, which is displayed in Figure 15(a), consists of 2000 polygonal cells. We observe in Figure 15(b) that the polygonal grid at the stopping time is still valid even in the vicinity of the wall. Let us emphasize that the lack of symmetry of this grid with respect to the axis y = 0 is uniquely due to the lack of symmetry of the initial grid. We have also plotted in Figure 14, the length, the kinetic energy and the internal energy of the rod versus time. Let us note that the numerical results obtained using the polygonal grid are superimposed with the finer Cartesian grid. Finally, plasticity threshold map is displayed in Figure 16 at the stopping time. We note that the plasticity threshold map obtained using the unstructured grid is quite similar to the one displayed in Figure 13(b) which has been obtained using the finest structured grid. These numerical results demonstrate that our numerical scheme is able to handle general unstructured grids as well as classical Cartesian grids. Moreover, the quality of the results remains the same.

Inria 7 Conclusion

We have presented a cell-centered Lagrangian scheme for solving elastic-plastic flows on twodimensional unstructured grids in planar geometry. The underlying material strength model relies on a hypo-elastic constitutive law wherein the Jaumann rate of the deviatoric stress tensor is expressed by means of the deviatoric strain rate tensor. The plasticity is taken into account by means of the von Mises yield criterion for perfectly plastic material. It has been implemented employing the classical Wilkins radial return algorithm [START_REF] Wilkins | Methods in Computational Physics, volume 3, chapter Calculation of Elastic-Plastic Flow[END_REF]. The present numerical scheme possesses several attractive features. First, it is characterized by a compatible cell-centered discretization which ensures the compatibility between the geometric conservation law and the mesh motion. Second, numerical fluxes are derived so that a dissipation inequality is satisfied at the semi-discrete level. Their computation is performed by means of a node-centered solver which results from the total energy conservation. Finally, the second-order extension is achieved by extending to the two-dimensional Lagrangian formalism, the Generalized Riemann Problem methodology initially introduced by Ben-Artzi and Falcovitz [START_REF] Ben-Artzi | Generalized Riemann problems in Computational Fluids Dynamics[END_REF].

In the future, we intend to improve the formulation of the node-centered solver by deriving and studying more sophisticated corner viscous matrices. We also plan to pursue the numerical validation by performing code to code comparison using more demanding two-dimensional test cases. Moreover, in a forthcoming paper, we shall derive the cylindrical version of the present scheme, addressing the difficult issue related to the problem of symmetry preservation.

• For any tensor A, the tensor exp(tA) is invertible and det[exp(tA)] = exp[t tr(A)] > 0.

(126)

• For any skew symmetric tensor, i.e., A t = -A, the tensor exponential, exp(tA) is a rotation.

It is interesting for the subsequent developments to recall the demonstration of the second property. Let us introduce Z(t) = exp(tA). First, using the first property and noticing that tr(A) = 0 yields immediately det[Z(t)] = 1. By definition, Z(t) satisfies the following equation

d dt Z = AZ(t), Z(0) = Id.
Similarly, Z t satisfies also d dt Z t = Z t (t)A t = -Z t (t)A.

Using these results, we compute

d dt (Z t Z)(t) = dZ t dt Z + Z t dZ dt = -Z t AZ + Z t AZ =0.
Hence, for all time t > 0, we have (ZZ t )(t) = (ZZ t )(0) = Id, which ends the proof.

A.2 A second-order accurate approximation of a tensor exponential

In this paragraph, we propose to construct a simple approximation of the tensor exponential, exp(tA) in the case where A is skew symmetric. To this end, we start by deriving an approximate solution of the ordinary differential equation whose the tensor exponential is the unique solution, that is dZ dt = AZ.

Let t n denotes an arbitrary time, we construct an approximate solution by integrating the above equation over the time interval [t n , t n + ∆t], where ∆t > 0 denotes the time step. Assuming that A is constant over [t n , t n + ∆t] leads to

Z(t n + ∆t) -Z(t n ) = A t n +∆t t n Z(t) dt.
Now, introducing the classical midpoint rule leads to the numerical scheme

Z n+1 -Z n = ∆t 2 A(Z n+1 + Z n ),
where Z n+1 and Z n denotes the approximations of Z(t n + ∆t) and Z(t n ). Finally, we get

Z n+1 = Id - ∆t 2 A -1 Id + ∆t 2 A Z n . ( 127
)

Inria

Here, (Id -∆t 2 A) is an invertible matrix since det(Id + Ω) = 1 + 1 2 | Ω | 2 for any arbitrary skew symmetric tensor Ω, refer to [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF]. We claim that (127) represents a second-order accurate time discretization of the differential equation (123). Thus, the tensor Q(t) given by

Q(t) = Id - t 2 A -1 Id + t 2 A , (128)

B Plane-wave analysis of the Wilkins model in the elastic regime

In this section, we aim at investigating the hyperbolicity of the system of equations [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF]. To this end we shall seek plane-wave solutions.

B.1 Hyperbolicity and plane-wave analysis

For sake of completeness, we briefly recall the definition of hyperbolicity for a system of conservation law in several space variables [START_REF] Godlewski | Hyperbolic Systems of Conservation Laws[END_REF]. Let Ω be an open subset of R p , and let F j , j = 1 . . . d, be d smooth functions from Ω into R p . The general form of a system of conservation laws writes

∂ ∂t U + d j=1 ∂ ∂x j F j (U ) = 0. (129)
Here, U = U (x, t) is the vector of unknowns, x ∈ R d and t > 0 denote the position vector and the time. The function F j is called the flux function and we denote by B j = B j (U ) its Jacobian matrix, defined as follows

B j = ∇ U F j (U ),
where ∇ U denotes the gradient operator with respect to the components of U . With this notation, the system (129) reads

∂ ∂t U + d j=1 B j (U ) ∂ ∂x j U = 0. ( 130
)
This equation corresponds to the non conservative form of the system (129). We are now in position to recall the definition of hyperbolicity. The system (129) is called hyperbolic if, for any U ∈ Ω and any ω ∈ R d , ω = 0, the matrix

B(U , ω) = d j=1 ω j B j (U ) (131)
has p real eigenvalues and p linearly independent corresponding eigenvectors. Now, let us show the links between the notion of hyperbolicty and the plane-wave analysis. The plane-wave analysis consists in looking for solution of system (129) which writes under the form U = U (x • nct). Here, n, is a unit vector of R d which characterizes the wave direction and c ∈ R is the wave speed. We introduce ξ = x • nct and denote by U the derivative of U with respect to ξ. Applying the chain rules leads to

∂ ∂t U = -cU , ∂ ∂x j U = n j U ,
where n j denotes the jth component of n. Substituting the previous formulas into the system written under the non conservative form (130) leads to

-cU + d j=1 n j B j (U )U = 0.
Utilizing the definition (131), we observe that this system rewrites

[B(U , n) -cId] U = 0.
This equation is nothing but the eigenvalue problem which corresponds to the study of the hyperbolicity of system (129). More precisely, the plane wave characterized by (n, c) exists if and only if c is a real eigenvalue of B. This shows the links between the plane-wave analysis and the hyperbolic feature of a system of conservation laws.

B.2 Study of the Wilkins model

Before proceeding any further, assuming a sufficient smoothness of the flow variables, the nonconservative form of the Wilkins model (1) writes

d dt P + ρa 2 ∇ • V - β ρ (S : D) = 0, (132a)
d dt V + 1 ρ ∇P - 1 ρ ∇ • S = 0, (132b)
d dt S -2µD 0 + SW -WS = 0. ( 132c
)
For sake of completeness, let us recall that D and W are respectively the symmetric and the skew symmetric part of the velocity gradient tensor:

D = 1 2 ∇V + (∇V ) t , W = 1 2 ∇V -(∇V ) t .
D is the strain rate tensor and W is the spin tensor. The deviatoric part of the strain rate tensor, D, is defined by

D 0 = D - 1 3 tr(D)Id.

Inria

Finally, the equation of state being written under the form P = P (ρ, ε), where ε is the specific internal energy, we write the thermodynamical coefficient α and β as

α = ∂P ∂ρ ε , β = ∂P ∂ε ρ .
This allows to express the thermodynamical sound speed as

a 2 = α + P ρ 2 β.
Let us point out that in the above system the first equation which describes the pressure evolution has been derived in Section 5.6. Let ϕ = ϕ(x, t) be a generic flow variable, we are looking for plane-wave solutions of the above system and written under the form ϕ = ϕ(x • nct). Here, n is a unit vector of R 3 which characterizes the direction of the wave and c is the wave velocity. Setting ξ = x • nct, we denote by ϕ (ξ) the derivative of ϕ with respect to ξ. Now, using the previous notation and the chain rule derivative leads to the following useful formulas

∇P = P n, (133a) ∇ • V = V • n, ∇ × V = n × V , (133b) ∇V = V ⊗ n, ∇ • S = S n. ( 133c
)
Recalling that the material derivative reads d dt ϕ = ∂ϕ ∂t + V • ∇ϕ yields the practical formulas

d dt P = -(c -V • n)P , (134a)
d dt V = -(c -V • n)V , (134b)
d dt S = -(c -V • n)S . (134c)
It is convenient for the subsequent computations to introduce w = c -V • n which is the relative speed of the wave with respect to the flow velocity projected onto the wave direction. We are now in position to express the strain rate tensor, its deviatoric part and the spin tensor in terms of V and n by employing (133)

D = 1 2 V ⊗ n + n ⊗ V , (135a)
D 0 = 1 2 V ⊗ n + n ⊗ V - 1 3 (V • n)Id, (135b)
W = 1 2 V ⊗ n -n ⊗ V . (135c)
We note in passing that the wave is longitudinal if and only if ∇ × V = 0 and transverse if and only if ∇ • V = 0. Substituting (133), ( 134) and (135) into system (132) yields

-wP + ρa 2 (V • n) - β ρ (V • Sn) = 0, (136)
-wV + 1 ρ P n - 1 ρ S n = 0, (137)
-wS -µ V ⊗ n + n ⊗ V - 2 3 (V • n)Id + 1 2 S(V ⊗ n) -S(n ⊗ V ) -(V ⊗ n)S + (n ⊗ V )S = 0. (138)

RR n°7975

Here, to compute the third term in the left-hand side of (136), we make use of the identity L : (a ⊗ b) = a • Lb. Multiplying (136) and (138) by n ρ and using the identity (a ⊗ b)c = a(b • c) leads to

w ρ P n =a 2 (n ⊗ n)V - β ρ 2 (n ⊗ Sn)V , (139a)
- w ρ S n = µ ρ [Id -(n ⊗ n)] V + 4 3 µ ρ (n ⊗ n)V - 1 2ρ [S -(Sn ⊗ n) -(Sn • n)Id + (n ⊗ Sn)] V . (139b)
Substituting (139a) and (139b) into (137) and assuming w = 0 leads to

w 2 V = H + C - β ρ 2 (n ⊗ Sn) V , (140)
where the tensors H and C are defined by

H = (a 2 + 4 3 µ ρ )(n ⊗ n) + µ ρ [Id -(n ⊗ n)] , (141a)
C = 1 2ρ [(Sn ⊗ n) -S + (Sn • n)Id -(n ⊗ Sn)] . (141b)
Let us point out that the case w = 0 corresponds to a wave speed given by c = V • n, that is, a plane wave which propagates with the flow velocity. Finally, let us define the acoustic tensor related to the Wilkins model as

A = H + C - β ρ 2 (n ⊗ Sn). ( 142
)
This tensor is decomposed into the three following terms:

• The tensor H which is symmetric positive definite and corresponds to the classical elastodynamic part of the model, refer to [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF].

• The tensor C which is non symmetric and comes from the Jaumann rate.

• The tensor -β ρ 2 (n ⊗ Sn) that corresponds to the dissipative feature of the model.

The plane-wave analysis consists in investigating the propagation condition

AV = w 2 V . ( 143
)
We claim that the model is hyperbolic if and only if the acoustic tensor, A, has real non-negative eigenvalues.

First, let us investigate the simple case for which we ignore the two last terms in the definition of A. Namely, A = H and it collapses to the acoustic tensor which characterizes the classical equations of elastodynamics. In this case, a L = (a 2 + 4 is the transverse sound speed. We point out that a 2 L and a 2 T are respectively the eigenvalues of the acoustic tensor, H, associated to the eigenvectors n and τ , where τ is a unit vector in the plane perpendicular to n. We note that in R 3 , a 2

T has a multiplicity of 2.

Inria Q t AQ =   An • n Aτ 1 • n Aτ 2 • n An • τ 1 Aτ 1 • τ 1 Aτ 2 • τ 1 An • τ 2 Aτ 1 • τ 2 Aτ 2 • τ 2   .
Employing (144a) and (144b) yields

Q t AQ =       a 2 L -β ρ 2 (Sn • n) -( 1 ρ + β ρ 2 )(Sn • τ 1 ) -( 1 ρ + β ρ 2 )(Sn • τ 2 ) 0 a 2 T + 1 2ρ [(Sn • n) -(Sτ 1 • τ 1 )] -1 2ρ (Sτ 2 • τ 1 ) 0 -1 2ρ (Sτ 1 • τ 2 ) a 2 T + 1 2ρ [(Sn • n) -(Sτ 2 • τ 2 )]      
We observe that the right lower block in the above matrix is a 2×2 symmetric matrix which admits two real eigenvalues. It is clear that these eigenvalues are the two remaining real eigenvalues of the acoustic tensor. Let us study the sign of these eigenvalues. To this end, we introduce the matrix

M = [a 2 T + 1 2ρ (Sn • n)]Id - 1 2ρ (Sτ 1 • τ 1 ) (Sτ 2 • τ 1 ) (Sτ 1 • τ 2 ) (Sτ 2 • τ 2 ) .
We point out that the matrix M is symmetric since S = S t . We note in passing that the 2 × 2 matrix in the right-hand side of the above equation is similar to S since it reads Q t 2 SQ 2 , where Q 2 is the 3 × 2 orthogonal matrix defined by Q 2 = [τ 1 , τ 2 ]. Bearing this in mind, we are going to localize the eigenvalues of M by studying (Mν • ν) where ν is a unit vector of R 2 . Using the definition of M leads to We conclude by claiming that the two remaining eigenvalues are real and non-negative granted that Y 0 µ 1. This shows that the acoustic tensor has three real non-negative eigenvalues and thus the Wilkins model is hyperbolic. We point out that it is also possible to study the effect of the plasticity on the hyperbolicity of the model. To this end, one might use the methodology developed by Trangenstein and Collela in [START_REF] Trangenstein | A High-Order Godunov Method for Modeling Finite Deformation in Elastic-Plastic Solids[END_REF].

(Mν • ν) = a 2 T + 1 2ρ (Sn • n) - 1 2ρ (SQ 2 ν • SQ 2 ν), (148)

  by computing the time rate of change of the yield function f . To this end, let us remark that since f =| S | -2 3 Y 0 , we have df dt = 1 |S| (S : dS dt ). Now, taking the inner product of (8) by S yields S : dS dt = 2µ [(D 0 : S) -(D p : S)] . Here, we have used the fact that (SW -WS) : S = 0 since S t = S and W t = -W. Recalling that tr(S) = 0 leads to (D 0 : S) = (D : S) and substituting the flow rule in the second term of the right-hand side of the above equation yields (D p : S) =χ S | S | : D S | S | : S , thanks to (11) =χ(S : D), since (S : S) =| S | 2 .

Figure 1 :

 1 Figure 1: Notation for a polygonal cell.

Figure 2 :

 2 Figure 2: Triangulation of the polygonal cell ω c .

Figure 3 :

 3 Figure 3: Sub-cell ω pc related to polygonal cell ω c and point p.

  This decomposition allows to rewrite the momentum flux as a summation of contributions coming from each sub-cell boundary ∂ωc Tn ds = p∈P(c) ∂ωpc∩∂ωc Tn ds. Hence, we define the sub-cell force related to cell c and point p F pc = ∂ωpc∩∂ωc Tn ds. (53) This definition enables us to rewrite momentum equation (32b) as m c d dt V c = p∈P(c) F pc , which is the Newton law applied to a particle of mass m c moving with velocity V c . We also use the sub-cell-based decomposition to write the total energy flux ∂ωc TV • n ds = p∈P(c) ∂ωpc∩∂ωc TV • n ds. The contribution of the sub-cell, ω pc to the total energy flux is expressed in terms of sub-cell force, F pc , by means of the approximation ∂ωpc∩∂ωc TV • n ds = ∂ωpc∩∂ωc Tn ds • V p since T t = T = F pc • V p , using sub-cell force definition (53).

  , we have employed (S c : D p c ) = (T c : D p c ) since tr(S c ) = 0 and D c = D e c + D p c . We have also expressed the time rate of change of the specific volume by means of (55a) and used the fact that DIV c (V ) = tr(D c ). Finally, introducing the decomposition of the Cauchy stress tensor, T c = -P c Id + S c , leads to

  its conservation involves d dt E = ∂D TV • n ds, where the right-hand side expresses the time rate of pressure work on the boundary, ∂D, of the domain, D, occupied by the fluid. By definition of total energy, this last equation rewrites c m c d dt E c = ∂D TV • n ds. (64)

Figure 4 :

 4 Figure 4: Fragment of a polygonal grid in the vicinity of point p located on the boundary ∂D.

Figure 5 :

 5 Figure 5: Notation related to the corner of a polygonal cell.

,tr c )

 c ≤ 0, the predicted trial elastic state given by (91) remains below the yield strength and thus the updated deviatoric stress, S n+1 c , is the same as the predicted value, i.e., S n+1 c = S n+1,tr c . On the other hand, if f (S n+1,tr c ) > 0, the incremental step RR n°7975

c) = 0 . 3 2 0 2 c

 0322 Bearing this in mind, equation (93) leads to | S n+1 c | -| S n+1,tr c |= -2µ c ∆t Λc . Now, recalling that the actual state is located on the yield surface leads to this expression of the plastic multiplier into (92) yields the value of the updated deviatoric stress Summary of the time discretization of the constitutive law Being given the initial deviatoric stress, S n c , the time centered values of the deviatoric strain, D n+ 1 ,c , and the spin tensor, W n+ 1 , then we compute the final value of the deviatoric stress by means of the following algorithm over the time interval [t n , t n+1 ] 1. Elastic prediction

  function, f , reads f (S) =| S | -2 3 Y 0 .

1 2 S

 2 xx for all time t > 0. Bearing this in mind, the magnitude of the deviatoric stress reads | S |= 3 2 | S xx |. Hence, the yield criterion becomes | S xx | -2 3 Y 0 ≤ 0. Finally, substituting (118a) into (119a) yields d

Figure 6 :

 6 Figure 6: Density and pressure resulting from the elastic-plastic piston-like problem at time t = 150 µs.

Figure 7 :

 7 Figure 7: Velocity and deviatoric stress resulting from the elastic-plastic piston-like problem at time t = 150 µs.

Figure 8 :

 8 Figure 8: Density and pressure from a flying plate striking a target plate at time t = 5 µs.

  Velocity in m/s.

  Deviatoric stress in Pa.

Figure 9 :

 9 Figure 9: Velocity and deviatoric stress from a flying plate striking a target plate at time t = 5 µs.

  (b) Balance of kinetic, internal and total energy.

Figure 10 :

 10 Figure 10: Collapse of a thick-walled cylindrical beryllium shell. Convergence analysis. Inria

Figure 11 :

 11 Figure 11: Collapse of a thick-walled cylindrical beryllium shell. Computational grid characterized by (40 × 32) uniformly spaced cells.

Figure 12 :

 12 Figure 12: Collapse of a thick-walled cylindrical beryllium shell. Color maps of the density and the plasticity threshold at time t = 130 10 -6 s. Computational grid characterized by (40 × 32) uniformly spaced cells.

Figure 13 :

 13 Figure 13: Two-dimensional projectile impact. Computational grid and color map of the plasticity threshold at time t = 5 10 -3 s. The computational grid is characterized by (200 × 40) uniformly spaced cells.

Figure 14 :

 14 Figure 14: Two-dimensional projectile impact. Convergence analysis.

Figure 15 :

 15 Figure 15: Two-dimensional projectile impact. Computational unstructured grids made of 2000 polygonal cells.

Figure 16 :

 16 Figure 16: Two-dimensional projectile impact. Plasticity threshold map at the stopping time.

where | Q 2 ν

 2 |= 1 since Q 2 is orthogonal. Let s min and s max be the smallest and the biggest eigenvalue of S, it is clear that (Sn • n) and (SQ 2 ν • SQ 2 ν) lie in the interval [s min , s max ]. Thus, (Mν • ν) satisfies the following inequality a 2 T -s 2ρ ≤ (Mν • ν) ≤ a 2 T + s 2ρ , where s = s maxs min . Accordingly, the eigenvalues of M lie in the interval [ 2µ-s 2ρ , 2µ+s 2ρ ]. Due to the von Mises yield criterion s ∈ [-2 2 3 Y 0 , 2 2 3 Y 0 ] and 2µs ≥ 2(µ -2 3 Y 0 ) > 0 since for the materials under consideration Y 0 µ 1.

  we recompute the time rate of change of specific entropy as

	ρθ	dη dt	=(S : D) -ρ	dε e dt	,
			=ρ	dε e dt	+	ε e µ	d dt	(µρ) -ρ	dε e dt	,
			=	(S : S) 4µ 2 ρ	d dt	(µρ).

  represents a second-order accurate approximation of the exponential tensor, exp(tA), i.e., exp(tA) = Q(t) + O(A 3 t 3 ) in the vicinity of t = 0. Let us remark that this approximation corresponds to the second Padé approximant of the exponential tensor. We conclude this paragraph by showing that Q(t) is a rotation for all t > 0. First, let us notice that

	det Q(t) = =	det(Id + t 2 A) 2 A) det(Id -t det(Id + t 2 A) det[(Id + t 2 A) t ]	since A t = -A
	=1.

Second, since A is skew symmetric then

Q t = (Id -t 2 A)(Id + t 2 A) -1 . Moreover, knowing that (Id + t 2 A) and (Id -t 2 A) are commuting leads to Q(t)Q t (t) = Id. Thus, Q(t) is a rotation.

RR n°7975

351, Cours de la Libération Bâtiment A 29

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgements

R. Abgrall has been funded in part by the EU ERC Advanced grant "ADDECCO" #226616.

Let us remark that the structure of the above systems is quite similar to the structure of the systems derived for computing the nodal velocity. To achieve the computation of the nodal acceleration, it remains to determine an expression of the cell-centered time derivatives of the Cauchy stress tensor and the velocity field.

Expression of the cell-centered time derivatives

Being given a piecewise monotonic linear reconstruction of the pressure, the deviatoric stress tensor and the velocity field within the cell c, we compute the cell-centered time derivatives of this physical variables by means of a Lax-Wendroff procedure. This procedure consists in expressing the time derivatives of the physical variables in terms of the spatial derivatives by means of the conservation laws satisfied by these variables. Let us recall that we need to compute

dt V is obtained thanks to the momentum equation which reads

Hence, the cell-centered time derivative of the velocity field writes

Here, (∇P ) c and (∇ • S) c denote respectively the cell-centered pressure gradient and the cellcentered deviatoric stress divergence computed by means of the piecewise monotonic linear reconstruction of the pressure and the deviatoric stress within cell c. Second, we compute d dt P recalling that the pressure is obtained by means of the equation of state under the form P = P (ρ, ε). Bearing this in mind and introducing the thermodynamic coefficients α = ∂P ∂ρ ε and β = ∂P ∂ε ρ yields

Assuming a sufficient smoothness for the flow variables, the time derivatives of the density and specific internal energy read

where D is the strain rate tensor. Now, replacing the above equations into (108) leads to

Here, we have used the decomposition of the Cauchy stress tensor T = -P Id + S. Noticing that a 2 = α + P ρ 2 β is nothing but the square of the thermodynamic sound speed, we arrive at the final expression

Finally, the cell-centered time derivative of the pressure writes

Inria

A The exponential of a tensor

A.1 Definition and main properties

The aim of this section is to recall briefly the definition and the main features of the exponential of a tensor following the presentation given in [START_REF] De Souza Neto | Computational Methods for Plasticity: Theory and Applications[END_REF][START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF]. To this end, let A be a given constant tensor and let us consider the initial value problem

for a tensor function Z(t). The above problem is nothing but a tensor ordinary differential equation for which the classical Cauchy Lipschitz theorem applies. Therefore, this problem has a unique solution which writes Z(t) = exp(tA) . This defines the tensor exponential function.

Explicitly it can be expressed by means of its Taylor series representation, that is for a generic tensor

Let us remark that the transpose of the tensor exponential, i.e., [exp(A)] t is defined as being the unique solution of the initial value problem

The tensor exponential satisfies the two following important properties:

RR n°7975

Bearing this in mind, let us perform the analysis of the complete case for which A is given by (142). We shall demonstrate that A admits 3 real non-negative eigenvalues. To this end, let us compute An and Aτ , where τ is a unit vector which belongs to the plane perpendicular to n

We observe that n is an eigenvector of A associated to the eigenvalue a 2 L -β ρ 2 (Sn•n). Let us prove that this eigenvalue is non-negative. Recalling that β = ( ∂P ∂ε ) ρ , for the Mie-Grüneisen equation of state, refer to Section 6.1, we obtain β = ρ 0 Γ 0 where ρ 0 is the density of the unstressed material and Γ 0 is the Grüneisen parameter. The deviatoric stress tensor, S, being symmetric, it admits the following spectral decomposition

where s i and ν i are the ith eigenvalue and eigenvector of S. We also point out that (ν 1 , ν 2 , ν 3 ) is an orthonormal basis of R 3 . Decomposing n over this basis as n =

Knowing that

which characterizes the elastic regime. Thus, the eigenvalue a 2 L -β ρ 2 (Sn•n) satisfies the following inequality

As we have already mentionned in Section 2.3, for most solids the shear modulus, µ, is two or three orders of magnitude larger than the yield strength, Y 0 , hence Y 0 µ 1. Bearing this in mind and granted that the density ratio remains bounded we claim that the eigenvalue a 2 L -β ρ 2 (Sn • n) is non-negative. Now, we demonstrate that A admits two more real non-negative real eigenvalues. To this end, let us introduce (τ 1 , τ 2 ) an orthonormal basis of the plane perpendicular to n. It is clear that (n, τ 1 , τ 2 ) is an orthonormal basis of R 3 . Moreover, the matrix Q whose columns are respectively n, τ 1 and τ 2 , i.e., Q = [n, τ 1 , τ 2 ], is an orthogonal matrix since Q t Q = Id. Thus, the acoustic tensor, A, is similar to Q t AQ where