## Abstract

A macroscopic model describing elastic–plastic solids is derived in a special case of the internal specific energy taken in separable form: it is the sum of a hydrodynamic part depending only on the density and entropy, and a shear part depending on other invariants of the Finger tensor. In particular, the relaxation terms are constructed compatible with the von Mises yield criteria. In addition, Maxwell-type material behaviour is shown up: the deviatoric part of the stress tensor decays during plastic deformations. Numerical examples show the ability of this model to deal with real physical phenomena.

## 1. Introduction

In finite deformation dynamics, two types of models of elastic–plastic solids exist. The first class is the class of hypoelastic models. This type of model is commonly used at engineering level and is implemented in many numerical codes such as OURANOS of the Direction Générale de l’Armement and the Commissariat à l’Energie Atomic (France) or a shock wave physics computer code (called CTH) developed at the Sandia National Laboratory (USA). For such models, the Maxwell relaxation equation for the deviatoric part of the stress tensor is used [1,2]. Compatibility of these models with the objectivity principle and with general thermodynamical principles is not straightforward.

The second class of models is the class of hyperelastic models ([3–10] and others). The model is hyperelastic if the stress tensor is defined in terms of a stored energy function. These models are, by construction, in conservative form. They are also objective and thermodynamically consistent. In papers by Plohr & Sharp [5], Walter *et al.* [7] and Miller & Collela [8], plastic deformations were defined, and new conservation laws for these new variables were proposed to close the system. Also, an equation on a hardening parameter was proposed to deal with evolution of the von Mises yield limit. Kluth & Desprès [11] proposed a simplified plasticity model by introducing a new potential inspired by the one used for the description of phase transitions. The total potential is ‘glued’ from two convex potentials: one describes elastic transformations and the other one describes plastic transformations. The boundary where the elastic–plastic transition occurs is compatible with the von Mises criterion. This model is isentropic for smooth solutions and has a restricted domain of application owing to its reversibility. An important class of hyperbolic type models describing the plastic behaviour of materials under large stresses (or the equivalent for small values of yield strength) was proposed by Godunov [3], Godunov & Romenskii [9] and Godunov & Peshkov [12]. This model describes, in particular, the evolution of the effective (elastic) deformation tensor in terms of which the stress tensor is defined. The relaxation equation for this tensor is compatible with the mass conservation law. At the end of the relaxation process, the stress tensor becomes spherical (the solid evolves as a fluid for large pressures).

We propose an extension of the last approach to the situation not necessarily involving small values of the yield strength. The internal specific energy is taken in separable form: it is the sum of the hydrodynamic part (depending only on the density and entropy) and the shear part (depending on other invariants of the Finger tensor). The relaxation terms of the model are built in agreement with the following requirements:

— the mass is conserved and does not evolve during the relaxation process,

— entropy is created when plasticity occurs (irreversibility of plastic deformations),

— the intensity of shear stresses decreases during the relaxation (Maxwell-type model), and

— the von Mises yield limit is reached at the end of the relaxation process.

This paper is organized as follows. In §2, the hyperelasticity model is summarized. In §3, the relaxation terms are added to model plastic deformations. These terms are compatible with the requirements mentioned above. Numerical results are presented in §§4 and 5.

## 2. Governing equations of isotropic elastic solids

The governing equations of *isotropic elastic* solids can be written in the following form [3,8,9,13]:
2.1
Here, are the columns of , is the deformation gradient (we call them *cobasis vectors* because they are the gradients of the Lagrangian coordinates), is the solid density, *ρ*_{0} is the reference density depending only on the Lagrangian coordinates, is the Finger tensor, **v** is the velocity field, *σ* is the stress tensor, is the specific total energy, is the specific internal energy depending only on invariants of and *η* is the specific entropy. The entropy is conserved along trajectories of smooth solutions of the system. The symmetric stress tensor is given by
where *p* is the hydrodynamic pressure and **S** is the deviatoric part of the stress tensor. Such a model of isotropic elastic bodies (2.1) is quite conventional. However, it should be closed by giving the equation of state . We take the energy in separable form,
2.2
The elastic part of the internal energy *ε*^{e}(**g**) depends only on **g**. The tensor **g** has a unit determinant, so it is unaffected by the volume change. The idea to take the arguments of the internal energy in the form *ρ*, *η* and **g** was proposed by Gouin & Debieve [14] (see also [10,15]), but for the dependence of the energy on the right Cauchy–Green tensor, . The choice of the Finger tensor is more natural for the Eulerian description of solids. The stress tensor is then
The shear part of the energy *ε*^{e} has no influence on the pressure *p*. The pressure is determined only by the hydrodynamic part *ε*^{h}. For applications, the hydrodynamic part of the energy *ε*^{h}(*ρ*,*s*) is taken in the form of stiffened gas EOS equation of state,
2.3
and the elastic energy *ε*^{e}(**g**) is taken as
2.4
where *μ* is the shear modulus. In the case of small deformations, we obtain the classical Hooke law from equation (2.2). For large solicitations, we recover the fluid behaviour (2.3). The choice (2.3) and (2.4) guarantees, in particular, the hyperbolicity of equations (2.1) in the one-dimensional case [16], and hence the well-posedness of the corresponding Cauchy problem. The stress tensor will then be
2.5
The first equation of (2.1) can also be written in non-conservative form
2.6
or replaced by an equivalent system of equations coupling the evolution of the vectors ,
2.7

If *ω*^{α}=0 initially, it will always be zero for any time. This form allows us to consider *ω*^{α} as a source term in the first equation. Obviously, this presentation is trivial in hyperelasticity.

## 3. Governing equations of elastic–plastic solids

The governing equations of *elastic–plastic* solids were derived in Godunov [3] and Godunov & Romenskii [9], Godunov & Peshkov [12] and Barton *et al.* [17]. We propose to take them in the following generic form:
3.1
The notations **e**^{β} are for the *elastic local cobasis*. In addition, we will define the *elastic Finger tensor* (it is also called the *effective elastic deformation tensor* by Godunov [3] and Godunov & Romenskii [9]), and the *elastic deformation gradient* **F** related to **e**^{β} by the relation **F**^{−T}=(**e**^{1},**e**^{2},**e**^{3}). The *real Finger tensor* is obtained in the following way. We calculate the trajectories of the solid body,
where **X** are the Lagrangian coordinates. Then, we calculate the real deformation gradient,
and finally the real Finger tensor,
The real Finger tensor verifies the equation
3.2
The unknown symmetric tensor **R**=**R**^{T} should be determined in accordance with some basic principles formulated below. The evolution equation for the elastic Finger tensor can be obtained from the evolution equations for **e**^{β},
3.3
To close this generic system (3.1), the relaxation term **R** and the relaxation time *τ*_{rel}>0 should be determined in such a way that

— they should be compatible with the mass conservation law,

— they should be compatible with the entropy inequality,

— the intensity of shear stresses decays during the relaxation (Maxwell-type model), and

— they should be compatible with the von Mises yield criterion.

### Remark 3.1

The form of equations used, for example, in Barton *et al.* [17] is different from the present paper. In particular, they use the evolution equation for the effective deformation gradient **F** supplemented by the evolution equation for the vector div(*ρ***F**) (here the divergence of a matrix is a vector, each component of which is the divergence of the matrix column). The fact that this vector is zero or not allows us to say whether the solid state is elastic or plastic. This is a sort of ‘indicator’ of a plastic state. A more general approach consists of using the Burgers tensor rot(**F**^{−T}) as an indicator of the plastic state. Let us finally remark that adding *ω*^{β}=(**e**^{β}) as unknowns, we can rewrite the evolution equations for **e**^{β} in conservative form analogous to equations (2.7),
and
An analogous conservative form can also be written for div(*ρ***F**) [17]. However, a non-conservative form is preferred for numerical purposes.

### (a) Compatibility with the mass conservation law

Let *ρ*_{0} be a function conserving along trajectories
It follows from the definition of density and equation (3.3) that
Hence, for mass conservation, it is *necessary* that

### (b) Compatibility with the entropy inequality

The energy equation
is equivalent to the following one:
Using the relaxation equation for **G** and the definition of the stress tensor, we get
because
Hence,
if
*A sufficient condition* for that is
where *a* is a positive scalar function, *a*>0. Let us also remark that with such a choice, the tensor **GR**=−*a***GS** is symmetric because, for isotropic materials, **S** is a second-order polynomial of **G**.

### (c) Compatibility with the von Mises criterion of yielding and the Lyapunov function

The von Mises criterion says that the material starts yielding when the corresponding *yield function*
becomes positive. The surface *f*(**S**)=0 is the *yield surface* and *Y* is the *yield strength*. If the yield function is negative, the material is elastic, and we use in the domain of elasticity.

Consider now the domain where the yield function is positive. The idea is to relax the deformations in such a way that at the end of the relaxation process, the yield surface is recovered. Consider tensor **G** in principal axes,
The density is then
The stress tensor will be
and
Consider the *singular value decomposition* of **F**^{−1}:**F**^{−1}=**UKV**^{T}, where **U** and **V** are orthogonal matrices and **K** is the diagonal matrix (e.g. [9]). The vectors **e**^{β} are the columns of the matrix (**F**^{T})^{−1}=**VKU**^{T}. The eigenvalues *κ*_{α} of the matrix **G**=**F**^{−T} **F**^{−1} are related to the eigenvalues *k*_{α} of the diagonal matrix **K** (or, equivalently, the singular values of (**F**^{T})^{−1}) by the following obvious formula:
Let us suppose that the matrices **U** and **V** do not vary during the relaxation process. This hypothesis was proposed by Godunov [3], Godunov & Romenskii [9] and Godunov & Peshkov [12]. It is related to the fact that the description of the Maxwell-like models in terms of deformation is not direct. So, we have to decompose the relaxation process into two parts: one of them is related to the geometry through the orthogonal matrices **U** and **V** and the other one to the physics through the eigenvalues of **K** (or **G**). The geometry varies only if the elastic deformations occur, while it is ‘frozen’ during the relaxation process. Since we know at the end of the relaxation process the eigenvalues *κ*_{α} (or *k*_{α}), we can easily reconstruct the local cobasis **e**^{β} since the matrices **U** and **V** are unchanged. The reconstructed local cobasis is not necessarily rotation free. If space variations do not occur during the relaxation process, the derivative *D*/*Dt* should be replaced by the partial derivative with respect to time,
and we will simply write
For each numerical cell, we want to solve the following relaxation equation:
3.4
since in the isotropic case, the tensor **GS** is symmetric. Or, in terms of eigenvalues of **G**,
3.5
Obviously, equation (3.5) (or, equivalently, equation (3.4)) admit the first integral, which is just the mass conservation:
The fact that the density is conserved during plastic deformations is quite conventional [3,18]. We will denote
We introduce new variables
In these variables, the relaxation equations are
where
3.6

#### (i) Lyapunov function

The aim of this section is to show that
is a Lyapunov function for our system. Replacing equation (3.6), we have
Differentiating it with respect to time, we obtain
Obviously, the time derivative of *L* is negative in the vicinity of the equilibrium where all *β*_{α}=1. A numerical study shows that this time derivative is indeed negative for any *β*_{α}>0, such that *β*_{1}*β*_{2}*β*_{3}=1.

Since only the ratio *a*/*τ*_{rel} is important, we choose *a* as
and the relaxation time in the form
3.7
where *τ*_{0} is a characteristic constant time and *n*>0. In particular, with such a choice of the relaxation time, the trajectories are attracted to the yield surface in finite time if 1>*n*>0, and in infinite time if *n*≥1.

In severe conditions, the deviatoric part of the stress tensor can be neglected, and we can take *Y* =0. In this particular case, formulae for the relaxation time can be found in Godunov *et al.* [19].

#### (ii) Simplified relaxation equation

Finally, the governing equations in terms of *β*_{α} are
3.8
Even if the ordinary differential equations for *β*_{α} can easily be solved numerically, it would be useful to have simple algebraic relations (exact solutions) of the time dependence of *β*_{α}, to make numerical codes more efficient. Here, we present such a simplification in the case of small deformations. Let *β*_{α}=1+*γ*_{α}, where *γ*_{α} is small. We can now estimate the principal eigenvalues of the deviatoric part of the stress tensor,

Equation (3.8) can be replaced by the approximate one
or, equivalently,
3.9
However, we will not suppose that the relaxation time is constant, we will choose it in the form
3.10
Then, equations (3.9) and (3.10) admit a time-dependent first integral,
The equations for *S*_{α} (3.9) become linear,
where
The solution of these equations is
3.11
This explicit formula reflects quite well the qualitative behaviour of the deviatoric part of the stress tensor. In particular, if the characteristic scale Δ*t* is large with respect to the relaxation scale *τ*_{0}, the final values of *S*_{α} can be estimated as
This asymptotic result can be used for the construction of efficient Riemann solvers in the case where the temporal time step is much larger than the characteristic relaxation time.

### Remark 3.2

To model the solid behaviour at large pressures where the solid behaves like a fluid, i.e. the stress tensor becomes spherical at the end of the relaxation process, Godunov [3] has proposed to take the relaxation equations for the singular values of the effective deformation gradient **F**. The singular value decomposition of **F** has been used, **F**=**VDU**^{T}, where **D**=**K**^{−1} (**K** is a diagonal matrix, and **U** and **V** are orthogonal matrices). For eigenvalues *d*_{i}, *i*=1,2,3, of the diagonal matrix
the following equation has been used:
where *τ*>0 is a relaxation time. This is compatible with the conservation of mass. The same relaxation equation was also used in Bonnetier *et al.* [4] and Merzhievsky & Resnyansky [6], with a modification taking into account the von Mises criteria. Our relaxation equation is different and specific to the choice of the equation of state in separable form (2.2). Under some natural constraints, this choice guarantees the hyperbolicity of the governing equations in the whole domain of parameters (e.g. [16]). We use the relaxation system for the eigenvalues of **G**, and show that the resulting equations are compatible, not only with the mass conservation law, but also the entropy inequality and the von Mises yield condition. It allows it to have a larger field of applications.

### Remark 3.3

Let us show that in the limit of small variations, the relaxation equations for the elastic deformation tensor are similar to conventional relaxation equations of viscoplasticity. Indeed, for small deformations, the difference between gives us the plastic deformation tensor, (for definitions, see Lemaitre & Chaboche [20]), Hence, the equation for the plastic deformation tensor following equations (3.2) and (3.3) is This is a classical equation in viscoplasticity [20].

## 4. Numerical results

### (a) Impact test problem

A higher order Godunov method has been used for the computations [21]. Details of the method are presented, in particular, in Gavrilyuk *et al.* [13] and Favrie *et al.* [16]. Some important modifications have been carried out for an efficient treatment of relaxation terms and the details will appear in Favrie & Gavrilyuk [22].

The problem of symmetric impact is considered under three different relative impact velocities: 100, 300 and 600 m s^{−1}. The materials have, in all these tests, the same initial pressure (*p*=0.1 MPa), the same initial density (*ρ*=7850 kg m^{−3}) and zero initial shear stress. Impacts of the two materials occur at *x*=0.5 m.

In figure 1, the output is shown at time *t*=0.1 ms. The 100 m s^{−1} impact (elastic) is represented by a dashed line, the 300 m s^{−1} impact by a thick line and the 600 m s^{−1} impact by a thin line. For a lower velocity impact, the plasticity limit is not reached. For a higher velocity impact, we notice the propagation of an elastic precursor followed by a plastic wave.

### (b) Quasi-static experiment

We study a quasi-static tensile experiment where a load–unload cycle is imposed (figure 2). Let an immobile rigid wall be on the right, and on the left, we consider a piston moving at a constant velocity of −5 m s^{−1} for 30 ms (), then we unload with a speed of 5 m s^{−1} for 50 ms (), and then we reload at a constant velocity of −5 m s^{−1} () for 20 ms. The sample length is 1 m. We take the relaxation time to be very small (of the order of 10^{−20} s) so that, under such conditions, the results do not depend on the loading speed. This is why we call this test ‘quasi-static’. The experimental configuration is shown in figure 2. In figure 3, we represent the total stress *σ*_{11}=−*p*+*S*_{11} and the deviatoric part of the stress tensor *S*_{11} as a function of strain (dimensionless piston position) during this load–unload cycle. In states 1, 3 and 5, the yield limit is attained, and the deviatoric part of the stress tensor does not vary anymore. When this limit is reached, the pressure increases (or decreases) as in a compressible fluid.

## 5. Diffuse interface model

In Favrie *et al.* [16], a diffuse interface model has been constructed for description of a fluid–elastic solid interaction. The dynamics of interfaces separating fluid and solid as well as the corresponding boundary conditions have been naturally included in such a formulation. An analogous ‘multi-phase’ formulation can be derived for the interaction of elastic–plastic solids with an ideal fluid. A model derivation will be published in a forthcoming paper. However, we would like to illustrate the ability of such a multi-phase formulation for the treatment of complex physical problems. Let us consider an impact of a solid copper projectile on a solid copper plate (figure 4). Such a problem has already been studied in an elastic case in Favrie *et al.* [16]. The copper projectile has an initial velocity *u*=500 m s^{−1}. This projectile is a square of 0.1 m length. The plate is 0.5 m long and 0.1 m wide. The contact between the plate and the projectile has already been made at the initial time. The other part of the computational domain contains air at atmospheric conditions. The computational domain is 0.7 m long and 0.7 m high. The mesh contains 1000 cells along the *x*- and *y*-directions. The parameters of the copper are as follows: *γ*_{s}=4.22, Pa, *μ*=9.2×10^{10} Pa, *Y* =2.49 GPa , *ρ*=8900 kg m^{−3}. The air is considered as an ideal gas with *γ*_{g}=1.4. A small amount of air (*α*_{g}=10^{−4}) is present in the solid. The presence of air allows to the shock-induced void nucleation and the appearance of new interfaces (cracks) in solids to be observed. The dynamics is shown at instant 249 μs in figure 5 (the elastic case is shown in figure 5*a*, and the elastic plastic case is shown in figure 5*b*).

A quantitative comparison of numerical solutions and conventional benchmark tests (the problem of a cylindrical rod that impacts a rigid surface, or perforation of thick plates by projectiles) are in progress.

## 6. Conclusion

An extension of a hyperelastic model is proposed to deal with elasto-plasticity. The model is compatible with the von Mises yield criteria and the entropy inequality. Numerical solutions of the model are in good qualitative agreement with the physics of plastic deformations. In particular, the model is able to predict the splitting of a signal into an elastic precursor and a plastic wave. In addition, the quasi-static experiments clearly show the irreversibility of the model. Perspectives to the description of solid–fluid interfaces have been discussed and some preliminary results showing the possibility of the model to generate new interfaces (cracks) have been presented.

## Acknowledgements

The authors thank Richard Saurel and Pierre Suquet for fruitful discussions about modelling and numerics. N.F. has largely benefited from a discussion with Prof. S. K. Godunov and Dr I. Peshkov during his post-doctoral study at Sobolev Institute of Mathematics, Novosibirsk. The authors have been partially supported by grant N 07.34.048 from the DGA.

## Footnotes

One contribution of 13 to a Theme Issue ‘The mathematical challenges and modelling of hydroelasticity’.

- This journal is © 2011 The Royal Society