## Abstract

This study describes a simple mathematical model for one-dimensional elastoplastic wave propagation in a metal in the regime where the applied stress greatly exceeds the yield stress. Attention is focused on the increasing ductility that occurs in the over-driven limit when the plastic wave speed approaches the elastic wave speed. Our model predicts that a plastic compression wave is unable to travel faster than the elastic wave speed, and instead splits into a compressive elastoplastic shock followed by a plastic expansion wave.

## 1. Introduction

Motivated by applications in the defence industry, we construct and analyse a mathematical model for the elastoplastic deformation of a metal subject to compressive stress considerably in excess of the yield stress. Specifically, we shall consider applied stress of the same order as the elastic moduli, giving rise to density variations of at least 10 per cent. For example, aluminium alloys, which have been well studied experimentally, have elastic moduli of the order 10–100 GPa, whereas the yield stress may be as low as 0.2 GPa. We will be working on a length scale of order 10^{−1} m, which is typical of many target specimens, and the elastic sound speeds are of order 10^{3} m s^{−1}. Typical particle velocities are also of order 10^{3} m s^{−1}, and it has been estimated that associated temperature rises are of order 10^{2} K.

Many phenomenological macroscopic models of the type we will consider have been proposed, which attempt to incorporate poorly understood microscopic phenomena empirically via adjustable parameters in the relevant constitutive assumptions. Although the resulting models are usually too complicated to be susceptible to any useful mathematical analysis, they underpin the hydrocodes that are used almost universally to predict macroscopic behaviour (Rakhmatulin 1945; Taylor 1946; Pack *et al.* 1948; Morland 1959; Von Kármán & Duwez 1959; Bland 1965; Johnson & Barker 1969; Prieto & Renero 1973; Swegle & Grady 1985; Steinberg & Lund 1989; Drumheller 1998; Molinari & Ravichandran 2004; Grady 2010).

Our objective is to extract from the literature a mathematical model that is both realistic and simple enough to be susceptible to systematic asymptotic and numerical analysis. It is fortunate for theoreticians that almost all the available experimental evidence concerns high velocity planar impact, at speeds of order 10^{3} m s^{−1}, on targets whose lateral extents greatly exceed their thicknesses. Such a target initially suffers only a uniaxial deformation, with the displacement varying only in the direction of impact. Within this restricted class of deformations, valuable experimental evidence comes from the results of Pack *et al.* (1948) and Johnson & Barker (1969). All these experiments reveal the presence of both classical elastic waves and plastic waves, as shown schematically in figure 1.

Figure 1 shows that, when the applied stress is below the yield stress, elastic waves, which could include shock waves, can propagate through the target with little decay. As the applied compressive stress increases above the yield stress then, provided the target is sufficiently thick, the leading elastic shock wave is observed to be the precursor for a slower plastic wave in which the flow variables are smooth functions of position and time, with the width of the plastic wave increasing with time (Johnson & Barker 1969).

However, when the applied stress becomes comparable to the elastic moduli, the plastic wave is observed to steepen markedly and possibly to catch up with the elastic wave. Also, it can propagate as a travelling wave over a relatively large distance of a few centimetres without significant change of shape. The associated dramatic increase in the metal ductility in this parameter regime is the main focus of the theoretical predictions we will be making in §3.

Concerning the basic physics modelling, the key phenomenon that we will consider is plastic flow at stresses well above the yield stress. It has been suggested that thermo-mechanical coupling is an important ingredient in this, especially in shock waves that are strong enough to engender stresses well in excess of the elastic moduli (Wallace 1981; Kalantar *et al.* 2005; Kadau *et al.* 2007; Hawreliak *et al.* 2008; Murphy *et al.* 2010). However, it is likely that dislocation nucleation and motion is the dominant microscale phenomenon (Davison & Graham 1979; Armstrong *et al.* 2009). Unfortunately, many different kinds of dislocation configurations are possible, such as classical slip systems involving dislocations, geometrically necessary or otherwise, partial dislocations and stacking faults, and twinned dislocations (Hirth & Lothe 1982; Armstrong *et al.* 2009). Hence, if we are to be able to create any coherent macroscopic mathematical model, it must be built on features that are common to averages over many such dislocation configurations.

Modelling such as that described in Johnson & Barker (1969), Clifton (1971) and Molinari & Ravichandran (2004) has paved the way for understanding violent elastoplastic flow on the macroscopic scales that are relevant to many practical applications. It is straightforward to write down conservation equations for mass and momentum in terms of the macroscopic material density, particle velocity and stress. Also, it is common to assume that the macroscopic stress is describable by a strain energy function of the macroscopic elastic strain. However, there is disagreement in the literature about how to distinguish the elastic and plastic components of the macroscopic displacement and stress (see Ilyushin 1963; Howell *et al.* 2009; Davison 2010, for discussion). These and related issues have also been approached from the rational mechanics viewpoint by many authors, including Willis (1969), Germain & Lee (1973) and Wang *et al.* (1995).

A key phenomenological hypothesis is that of Orowan (1954), which is the cornerstone of many influential theories of rate-dependent plasticity (often called viscoplasticity (Johnson & Barker 1969; Clifton 1971; Molinari & Ravichandran 2004)). The hypothesis is that, in a uniaxial deformation, the macroscopic plastic strain rate is proportional to the product of an average velocity of mobile dislocations and an average dislocation density. As shown in Davison (2010), further modelling hypotheses then need to be made concerning the dependence of the average dislocation velocity on the macroscopic shear stress and the dependence of the dislocation density on variables that may include the plastic strain history. One of the most important variables to emerge from such microscale modelling is the number density of mobile dislocations, which is a measure of the departure of the stress from the yield surface. Hence, we will adopt a simplified flow rule, in which the plastic strain rate is linearly proportional to the distance from the stress to the yield surface. This flow rule can be thought of as a relaxation law for the departure of the stress from the yield surface; as the time scale for this relaxation tends to zero, the model reduces to the simpler one of rate-independent plasticity, in which the yield stress is never exceeded. When our flow rule is combined with the laws of conservation of mass and momentum, we will be led to a system of three nonlinear partial differential equations and two nonlinear algebraic equations for the axial elastic and plastic strains, the material density and velocity, and the axial stress. The situation suggests analogies with the mathematical theories of combustion and relaxation in gas dynamics (Ockendon & Spence 1969).

In §2, we will derive our basic model for rate-dependent uniaxial compression. Then, in §3, we will study the predictions of our model for the waves that result from a violent compressive stress applied at the boundary of a one-dimensional target. This will lead us to a surprising prediction concerning the propagation of so-called over-driven plastic waves that travel at the elastic wave speed.

## 2. One-dimensional elastoplastic model

### (a) Governing equations

We focus our attention on one-dimensional uniaxial motion in which a material element initially at the Lagrangian position (*X*,*Y*,*Z*) is displaced to the Eulerian position (*x*(*X*,*t*),*Y*,*Z*) after a time *t*. Following, for example, Davison (2010) or Drumheller (1998), conservation of mass leads to the relation
2.1between the density *ρ*(*x*,*t*) and its initial value *ρ*_{0}, which we will assume to be constant. When we use the multiplicative decomposition proposed by Lee (1969) and others (see Germain & Lee 1973; who specifically address elastoplastic shock waves), the one-dimensional deformation gradient tensor reads
2.2where the elastic and plastic deformation gradients are
2.32.4

We suppose also that the stress field is uniaxial, so the Cauchy stress tensor takes the form ** τ**=diag(

*τ*

_{1},

*τ*

_{2},

*τ*

_{2}), where

*τ*

_{1}and

*τ*

_{2}are functions of

*X*and

*t*, and the momentum equation then reads 2.5The stress is assumed to be derivable from an energy density function , depending solely on the elastic strain, such that 2.6where the factor of 2 in the expression for

*τ*

_{2}arises from the assumed uniaxiality.

When considering the yield criterion to be applied, we focus on the simplest case of an isotropic, non-hardening material. In rate-independent plasticity, this leads to a yield condition of the form
2.7where *τ*_{Y} is the yield stress, assumed to be a prescribed constant, and *f* is the yield function. Noting that energy is dissipated by the material at a rate
2.8per unit volume, we make the common assumption that the stress field in the material when it is flowing plastically is such as to maximize the dissipation rate, while obeying the inequality (2.7). Hence, we obtain the associated flow rule
2.9where *Λ* is a Lagrange multiplier with the dimension of inverse time. Elastic behaviour corresponds to *Λ*=0 with strict inequality in (2.7). On the other hand, when the material has yielded, *Λ* must be non-negative to preclude negative dissipation. We can therefore extend (2.7) to the complementarity conditions
2.10

In practice, stresses significantly exceeding the yield stress are observed in experiments (Kalantar *et al.* 2005; Kadau *et al.* 2007; Hawreliak *et al.* 2008; Murphy *et al.* 2010). This motivates us to relax the yield condition (2.7) and to generalize (2.10) to a rate-dependent model, in which a constitutive relation is imposed between *Λ* and the yield function *f*(*τ*_{1},*τ*_{2}), with the property that *Λ*=0 whenever *f*<*τ*_{Y}. We will focus on the specific choice^{1}
2.11where *k* is a positive constant; the rate-independent complementarity conditions (2.10) are recovered in the limit .

To summarize, once we have specified constitutive relations for and *f*, the equations (2.1), (2.4)–(2.6), (2.9) and (2.11), in principle, give us a closed system for the nine unknowns *ρ*, *x*, , , , , *τ*_{1}, *τ*_{2} and *Λ*, all functions of *X* and *t*.

We will suppose that a characteristic time scale *T* is imposed by the boundary conditions and non-dimensionalize stress using the bulk elastic modulus *K*, which will shortly be seen to characterize the speed of small-amplitude plastic waves (Pack *et al.* 1948). We therefore write
2.12a
2.12bhenceforth dropping the primes on dimensionless variables to reduce clutter. There are only two dimensionless parameters that we choose to be
2.13and we will be particularly interested in the limit , modelling situations in which the imposed stress is far in excess of the yield stress.

To close our model, we need to select particular forms of the functions *f* and . For mathematical convenience, we use the von Mises yield criterion, which leads to the dimensionless yield function^{2}
2.14The flow rules (2.9) therefore read
2.15in dimensionless variables. A first integral of these equations yields
2.16so that the plastic deformation is incompressible, as expected when using the associated flow rule with a pressure-indifferent yield function.

For simplicity, we limit our attention to a mechanically linear constitutive relation modelled by the dimensional strain energy density
2.17where *λ* and *μ* are the usual Lamé constants (with *K*=*λ*+2*μ*/3), and the uniaxial elastic strain is given by
2.18In dimensionless form, (2.17) reads
2.19where
2.20is the dimensionless elastic wave speed, so that for materials with positive Poisson's ratio *ν*.

Of course, there are many other possible constitutive relations one might impose. In particular, in view of the large density changes that can occur, several authors have incorporated nonlinear dependence on the density, while assuming linearity with respect to the relatively small deviatoric elastic strain (Green & Naghdi 1965; Willis 1969; Wang *et al.* 1995). By contrast, in this study, we find that especially interesting behaviour occurs in a weakly nonlinear limit where the isotropic and deviatoric components are comparable.

We can now formulate a system of equations for the velocity *v*=∂*x*/∂*t*, *ρ*, , , *τ*_{1} and *τ*_{2} by writing (2.1) and (2.5) as
2.21and using (2.4), (2.6), (2.16) and (2.19) to give
2.22and
2.23Also (2.9) gives
2.24where
2.25

Equations (2.21)–(2.24) form the basis for the rest of the paper. In §3, we will analyse in detail their solution for the specific problem of unidirectional wave propagation in a semi-infinite slab of material. Before doing so, we first briefly demonstrate that they are consistent with other models in the literature in the rate-independent limit.

### (b) Rate-independent extreme plasticity

In the rate-independent limit , the flow rule (2.24) reduces to an algebraic relation between *τ*_{1} and *τ*_{2} and hence, through (2.23), to a relation between and . In principle, elimination between (2.22) and (2.23) then defines *τ*_{1} as a function of *ρ*. Therefore, as pointed out by Rakhmatulin (1945) and Ilyushin (1963), the equations of motion (2.21) reduce to the equations of one-dimensional barotropic flow. Some caution with this conclusion is warranted, however, because the effective equation of state between *τ*_{1} and *ρ* is in general multi-valued, and also because the limit is singular.

If in addition we take the extreme plasticity limit , in which the applied stress greatly exceeds the yield stress, then the yield condition (2.24) simply reduces to *τ*_{1}−*τ*_{2}=0, and (2.23) then leads to the relation
2.26Evidently there are multiple solution branches, as pointed out already, but the branch relevant to an initial undeformed state with is clearly . Thus, in this limit, the effective equation of state between the axial pressure and the density may be found explicitly and is given by
2.27

For density , this leads to a well-posed barotropic flow model with sound speed *c* given by
2.28It follows that *c* is unconditionally less than the dimensionless elastic wave speed *c*_{e}. This explains the appearance of an elastic precursor wave travelling ahead of a slower plastic wave, as sketched in figure 1 and further explored in §3. We observe also that the dimensionless wave-speed for small-amplitude perturbations from the initial unstressed state is given by *c*=1 when *ρ*=1. This corresponds to a dimensional plastic wave-speed of , as anticipated in (2.12).

## 3. Unidirectional wave propagation

### (a) Initial and boundary conditions

Motivated by the experimental evidence cited in §1, we focus on the problem of an initially unstressed semi-infinite slab of material occupying the half-space , whose face *X*=0 is displaced at a specified positive velocity *V* _{0}(*t*). This problem is of especial interest for violent impact testing, and our aim is to obtain qualitative agreement with figure 1 when we ignore reflections from the rear face of the target. The initial and boundary conditions to be imposed are therefore
3.1a
3.1bWe will treat *V* _{0}(*t*) as an *O*(1) function of time that corresponds to a dimensional boundary velocity comparable to the plastic wave-speed . We assume that *V* _{0}(*t*) is an increasing function which tends monotonically to a maximum value *V* _{m} as and that the time scale *T* has been chosen such that *V* _{0}(*t*)∼*t* as .

For 0<*ε*≪1, our model predicts that an elastic precursor wave travels in a narrow region of width *O*(*ε*) behind the leading characteristic *X*=*c*_{e}*t*, as shown in figure 2. As we cross this region, decreases and stays equal to unity. Hence, from (2.23), the maximum shear stress *τ*_{2}−*τ*_{1} increases to reach the yield stress and, by solving the linear elastic equations, it is easy to show that this happens on a characteristic *X*=*c*_{e}(*t*−*t*_{*}). Then, because the plastic wave speed is initially less than *c*_{e}, there is a region **C** of constant stress, velocity and density, in accordance with the ‘plateau’ in figure 1. In this region, both the purely elastic equations and the yield condition *τ*_{2}−*τ*_{1}=*ε* are satisfied, and it follows that the elastic strain components are given by and
3.2We deduce that the solution structure envisaged in figure 2 is possible only if the parameters *ε* and *c*_{e} satisfy the inequality
3.3which reduces to
3.4as .

The region **C** is bounded above by an elastoplastic characteristic of unit slope, beyond which the full elastoplastic model has to be solved in **P**. In this region, *τ*_{2}>*τ*_{1} and so from (2.15), we see that will decrease, which reflects the elastic response to the increase in the incompressible plastic deformation . In the rate-independent limit , we anticipate by analogy with compressive ‘piston’ problems in gas dynamics that the characteristics of the hyperbolic system (2.21) will converge to form a shock at some finite time *t*_{s}, as indicated in figure 2. In §3b, we derive a simplified weakly nonlinear model that captures all the key features anticipated in figure 2 while being more amenable to analytical and numerical solution.

### (b) The weakly nonlinear limit

In order to gain insights into the long-time evolution of these waves for either rate-dependent or rate-independent models, we carry out a weakly nonlinear analysis in which we assume that both *ε* and *c*_{e}−1 are small. With hindsight, we focus on the asymptotic limit where (*c*_{e}−1)=*O*(*ε*^{1/2});^{3} the subsequent analysis will reveal that, when (*c*_{e}−1)≫*ε*^{1/2}, elastic waves travel so fast that they cannot be overtaken by plastic waves and the scenario depicted in figure 2 persists indefinitely. Hence, we write *c*_{e}=1+*ε*^{1/2}*γ*, where *γ* is *O*(1): the inequality (3.4) requires that
3.5and for typical metals, we estimate that *γ* takes values approximately in the range 2.64 (aluminium alloy)–9 (copper). We anticipate that and *τ*_{2} will all be of *O*(*ε*^{1/2}) also.

We therefore consider the weakly nonlinear regime in which 3.6which implies that the stress components may then be expanded in the form 3.7and 3.8and hence 3.9

To lowest order in *ε*, the governing equations (2.21) and (2.24) therefore become
3.10
and
3.11where *α*=*ε*^{1/2}*a*. The initial and boundary conditions (3.1) give
3.12where *V* (*t*)=*ε*^{−1/2}*V* _{0}(*t*) is assumed to be monotonic increasing, concave and bounded, so that
3.13

The material remains elastic, with *ϕ*=0, until the yield condition *ϱ*(*ϱ*−3*γ*)=−1 is satisfied. This will occur at *t*=*t*_{*}, where,
3.14The lower bound (3.5) on *γ* ensures that *V* _{*} is real and lies in the range (0,1). We assume that the maximum imposed velocity *V* _{m} exceeds the value *V* _{*} so that the material will become plastic when , and *ϕ* will then be determined by the rate equation (3.11).

The theory of weakly nonlinear wave propagation (Whitham 1974) tells us that a distinguished limit emerges over a long time scale, where *t*=*O*(*ε*^{−1/2}). In terms of the travelling-wave variables
3.15we find that the governing equations (2.21) and (2.24) take the form
3.16
3.17
and
3.18

By matching with the unstrained material ahead of the travelling wave, we obtain the far-field behaviour
3.19We thus deduce from (3.16) and (3.17) that
3.20to lowest order in *ε*. Then, by adding (3.16) and (3.17) to eliminate the dominant term, we obtain the equation
3.21while (3.18) simplifies to
3.22Matching with the elastic solution as yields the initial condition
3.23Equations (3.21) and (3.22) form a classical hyperbolic system of partial differential equations for the scaled velocity and the plastic strain *ϕ*, with characteristics given by *τ*=constant and d*η*/d*τ*=2*ϕ*+*γ*. However, the rate-independent limit is singular.

So long as the material remains elastic, with *ϕ*=0, (3.21) reduces to
3.24and the solution is therefore
3.25

For *η*<*γτ*, the material is plastic and we have to solve the nonlinear system (3.21) and (3.22). Analytical progress is only possible in the rate-independent limit when (3.22) becomes
3.26which is a monotone increasing function for . We henceforth assume that *V* _{m}<3*γ*−1 so that (3.26) defines a one-to-one relation between and *ϕ*. It is now straightforward to solve
3.27using the method of characteristics, because *b* is also a monotonic increasing function of . The characteristics for a typical solution are shown in figure 3; this is the equivalent of figure 2 in the (*η*,*τ*)-plane. Here, we have used the parameter value *γ*=1 and applied a boundary velocity of the form
3.28with *V* _{m}=1.

We see in figure 3 that a shock forms at time
3.29Thereafter, as long as the shock remains at the **C**/**P** interface, it is tempting to simply assume the Rankine–Hugoniot condition for (3.27) in the form
3.30However, the shock conditions should properly be derived from the conservation laws (2.21), which imply that
3.31where the Eulerian shock position is given by
3.32Fortunately, when the expansions (3.6) are used in (3.31) and expanded to lowest order, they yield
3.33and
3.34the second of which is equivalent to (3.30) on using (3.26).

We can now continue the solution by using the condition (3.30) to determine the shock position while solving the hyperbolic equation (3.27) in the plastic region **P**. The result of such a calculation is shown in figure 3, with the shock position plotted as a heavy dashed curve.

We emphasize that this solution procedure works only for the rate-independent model where , and it is by no means certain that such solutions could be obtained from the full rate-dependent model ((3.21) and (3.22)) in the singular limit . We note, for example, that [*ϕ*]^{+}_{−} is in general non-zero across a shock solution to (3.27), which could not be the case when *α* is finite.

In addition, there is no obvious entropy inequality in our model and hence we may need to exclude solutions that violate causality. It may be shown that the Rankine–Hugoniot condition (3.30) guarantees that the characteristics of (3.27) are directed into the shock (as shown in figure 3), which suggests that such a shock should be stable. However, it could arise as the limit of a rate-dependent solution as only if causality were to be satisfied with respect to characteristics of the full system (3.21) and (3.22).

Finally, we find that (3.30) makes it possible for the shock to accelerate past the elastic wave-speed and overtake the **C**/**E** interface, in which case (3.27) does not apply on the right-hand side of what must now be considered an elastic/plastic free boundary. The conditions (3.31) still hold across this free boundary but now we can use only (3.26) to relate *ϕ* to in **P**, while *ϕ*=0 in **E**, which means that (3.30) must be replaced by
3.35and the solution must again be checked for causality.

It is easy to see that this situation will arise whenever the imposed velocity *V* (*t*) exceeds a critical value
3.36which lies between *V* _{*} and (3*γ*−1). We show an example in figure 4, again using a boundary velocity of the form (3.28) with now *γ*=4 and *V* _{m}=10>*V* _{c}≈7.92. With these parameter values, *t*_{*}≈0.0084 and the elastic region is therefore very narrow. We see a constant region between the elastic and plastic waves, as expected, and a shock that starts at *τ*≈0.1 before propagating rapidly into the constant region and eventually overtaking the elastic precursor. The corresponding solution is plotted in figure 5. Now we see a large-amplitude plastic wave that steepens into a shock and rapidly accelerates, eventually absorbing the relatively feeble elastic precursor. This behaviour is characteristic of experimentally observed over-driven waves (Kalantar *et al.* 2005; Kadau *et al.* 2007; Hawreliak *et al.* 2008; Murphy *et al.* 2010) in which dramatic steepening heralds the plastic wave overtaking the elastic precursor. We note that *V* _{c} increases as *γ* increases, making the occurrence of over-driven waves increasingly unlikely, and this motivates our focus on the parameter regime where *γ*=*O*(1).

### (c) Numerical solution

The solution of (3.22) for *ϕ*, subject to the matching condition as , is
3.37in terms of the *a priori* unknown , where refers to the quadratic function
3.38We can then write (3.21) as a partial integro-differential equation for . For computational convenience, we incorporate a viscous diffusive term to regularize any shocks that may arise and thus obtain the equation
3.39where 0<*δ*≪1. Equation (3.39) is to be solved subject to the initial condition (3.23) at *τ*=0.

In figure 6, we show numerical solutions of this model for *V* _{m}<*V* _{c} that are reminiscent of the time reversal of the curves in figure 1. The imposed velocity is again given by (3.28), and we use the parameter values *γ*=1 and *V* _{m}=1. We show snapshots of and *ϕ*(*η*,*τ*) versus *η* at *τ*=5 for two different values of *α*=1,10, as well as the rate-independent solution corresponding to , which is plotted as a dashed curve. We observe that, although is continuous for finite *α*, it has slope discontinuities at the interfaces between the undisturbed, elastic and constant regions.

For *η*>*γτ*−*t*_{*} (with *t*_{*}≈0.48), the material remains elastic, and the rate-dependent and rate-independent solutions are identical (modulo a small amount of artificial diffusion). Once the material has yielded (i.e. for *v*>*V* _{*}), we again observe a plateau in which and *ϕ*=0, followed by a plastic wave. In the rate-independent limit, this wave takes the form of a shock, across which both and *ϕ* are discontinuous. With these parameter values, for which *V* _{m}<*V* _{c} and the plastic shock is not over-driven, we see that the main effect of rate dependence is to smooth out this shock into a travelling wave whose width decreases with increasing *α*.

These observations lead us to seek wave solutions of (3.21) and (3.22) that travel with speed *β* relative to the (*η*,*τ*) frame, so that from (3.21),
3.40with
3.41a
3.41bwhere
3.42and *z*=*η*−*βτ*. By integrating and eliminating , we can soon show that
3.43in accordance with the conservation condition (3.35) when we put and . Hence, the smoothed shock travels more slowly than the effective elastic wave speed *γ* as long as *V* _{m}<*V* _{c}. However, a local analysis in the limit *V* _{m}↗*V* _{c} reveals the appearance of a slope singularity in *ϕ* and corresponding jump in as shown in figure 7. Hence, when *V* _{m}>*V* _{c}, there are no travelling waves satisfying (3.41) in which *ϕ* is continuous, suggesting that no over-driven elastic/plastic free boundaries can exist.

The situation is reminiscent of that in relaxing gas dynamics (Ockendon & Spence 1969), where the fact that the gas molecules cannot instantaneously return to their equilibrium configuration means that no compressive wave can propagate faster than the so-called frozen wave speed. This forces the formation of a partially dispersed shock which propagates at the frozen wave speed, resembling the dashed limiting solution shown in figure 7.

Armed with these insights, we now compute numerical solutions of the rate-dependent model ((3.37) and (3.39)) for the over-driven case in which *V* _{m}>*V* _{c}. In figure 8 we show the evolution of and *ϕ* with *α*=1 and *γ*=2, so that *V* _{*}≈0.1716 and *V* _{c}≈3.845. We use the expression (3.28) for the boundary velocity, here choosing *V* _{m}=4.5∈(*V* _{c},3*γ*−1).

As expected, the material is undisturbed ahead of the leading characteristic *η*=*γτ*, behind which a narrow elastic region is followed by an expanding constant region in which and *ϕ*=0. Behind this, we see a plastic wave that steepens to form a shock, across which is discontinuous, while *ϕ* suffers only a jump in its first derivative. Thereafter, the shock propagates at the elastic wave speed *γ*, and the value of immediately behind the shock increases, ultimately attaining a fixed value. Thus, we again have a partially dispersed shock, with transition back to the upstream values and *ϕ*=*ϕ*_{m} occurring via an unexpected smooth plastic expansion wave that lags behind the discontinuity. Consequently, a plateau opens up between the shock and the following plastic wave, in which takes an approximately constant value significantly in excess of *V* _{m}, while *ϕ* remains approximately equal to zero. (The small increase in *ϕ* across the shock is a diffusive artefact which decreases as the viscosity *δ* is reduced.)

This numerical experiment immediately suggests that the rate-independent results shown in figures 4 and 5 cannot be realized as large-*α* limits of solutions of the evolutionary rate-dependent model. The key new features revealed by figure 8 are, respectively, the existence of a smooth wave that lags behind the jump discontinuity in and the persistent continuity of *ϕ*. From (3.7) and (3.20), we see that the pressure is proportional to and hence the smooth wave observed in figure 8 is an expansion wave that tends to an expansion shock as . The continuity of *ϕ* corresponds to continuity of the transverse elastic strain , which is enforced by the rate equation (2.24). Ultimately, this is explained by the fact that the plastic strain is continuous across the leading shock, while the axial elastic strain suffers a jump discontinuity.

Because *ϕ* is continuous, the expression (3.34) for the shock speed simplifies to
3.44and a shock propagating into a region where and *ϕ*=0 must therefore propagate at the elastic wave-speed *γ*. This further demonstrates that rate-independent solutions containing an elastic/plastic free boundary are disallowed when rate-dependence is included, since they travel faster than *γ*.

The existence of the expansion wave is more difficult to explain analytically. However, figure 8 suggests that soon attains an almost constant value behind the **C/P** shock on which *dη*/*dτ*=*γ*. When we use this shock speed in (3.21) and (3.22), we find that the value of behind the shock satisfies
3.45so that, once exceeds 3*γ*/2, it rapidly tends to 1/*V* _{*} over a time scale of order *α*^{−1}. This is confirmed by inspection of figure 8, where the value of in the plateau region is found to be 1/*V* _{*}≈5.828.

This motivates us to think of the expansion wave as a travelling wave of width *O*(*α*^{−1}) as in (3.40), with (3.41a) replaced by
3.46The resulting profiles for *ϕ* and are plotted in figure 9 with *γ*=2 and different values of *V* _{m}>*V* _{c}≈3.845. The dashed curves showing the limiting case as *V* _{m}↘*V* _{c} are identical to those shown in figure 7, with a slope discontinuity in *ϕ* and a jump discontinuity in , and the amplitudes of the waves decrease as *V* _{m} increases from *V* _{c}.

As , the expansion wave can be considered as an expansion shock for the rate independent model (3.27), and it can easily be shown that the corresponding characteristics make this shock causal.

### (d) Summary

The results of this section put us in a position to propose a mathematically consistent scenario for elastoplastic waves under a compressive stress.

— For under-driven waves, the rate-dependent model ((3.21) and (3.22)) always has smooth solutions. As , these tend to weak solutions of the rate-independent model (3.27) that may have jump discontinuities satisfying the Rankine–Hugoniot condition (3.34) as in figure 6.

— When the applied compression is strong enough to induce an elastoplastic shock to travel at the elastic wave-speed (i.e. when

*V*_{m}>*V*_{c}), there is a rapid decompression behind that shock that propagates locally as a smooth travelling wave. As , this expansion wave steepens into an expansion shock that lags behind the leading elastoplastic shock, as in figure 8. Moreover, the expansion shock can be regarded as part of a weak solution of (3.27). Its stability is strongly suggested by the facts that it is a limit of a rate-dependent travelling wave, and that it is causal with respect to the characteristics of (3.27).

## 4. Discussion and conclusion

We have presented a systematic mathematical framework within which to describe one-dimensional compressive elastoplastic wave propagation when the ratio *ε* of yield stress to typical applied stress is small. Our principal new insight concerns the material response when the applied stress is so large that elastoplastic waves may catch up with the elastic precursor, leading to over-driven shocks. For ease of exposition, we have focused attention on a weakly nonlinear parameter regime in which the elastic and plastic acoustic wave speeds differ by a relative amount of *O*(*ε*^{1/2}) and wave steepening occurs over a relatively long time scale.

We find that over-driven shocks occur when the maximum applied velocity *V* _{m} exceeds a critical value *V* _{c}. For *V* _{m}<*V* _{c}, the effect of rate dependence is simply to smooth the plastic/plastic shock waves that result from nonlinear wave steepening in the rate-dependent limit. However, for *V* _{m}>*V* _{c}, the effect of rate dependence becomes much more dramatic. Now, even the rate-dependent response begins with a compressive elastoplastic shock that involves a jump discontinuity in the velocity, density and pressure (but not in the plastic strain). Moreover, this discontinuity moves at precisely the elastic wave speed and has a prescribed amplitude independent of the strength of the applied compressive stress. Behind the leading shock, the velocity takes a fixed value 1/*V* _{*} and returns to the imposed boundary velocity *V* _{m} via a plastic expansion wave, which steepens into an expansion shock in the rate-independent limit. This shock could not have been predicted using the rate-independent model alone. Unfortunately, all the readily available experimental data refer to impacts in which *V* _{m}<*V* _{c} and hence we are unaware of any experimental evidence of plastic expansion waves.

We emphasize that all these predictions depend strongly on our assumption of one-dimensionality. Also, all our detailed calculations depend on assumptions of a prescribed constant yield stress, mechanical linearity of the elastic response, and our simple linear rate equation for the deviation of the stress from the yield surface. However, our mathematical model appears to be quite robust to small changes in any of these assumptions. Hence, it is possible that our predictions may be relevant to, for example, the observation of confinement stresses in processes such as shaped charge penetration (Novokshanov & Ockendon 2006).

## Acknowledgements

We are grateful to AWE for initiating our interest in this class of problems and especially to Dr Chris Robinson and Dr John Curtis for many useful discussions. We are also indebted to Prof. Dennis Grady and Prof. Justin Wark for their helpful advice. We thank an anonymous referee whose careful reading and many useful comments have helped us to improve the accuracy and clarity of this manuscript. J.R.O. is grateful to the Leverhulme Trust for support that has facilitated the work contained in this paper.

## Footnotes

↵1 [

*z*]^{+}denotes the positive part of*z*, equal to zero when*z*<0 and equal to*z*when*z*≥0.↵2 Given our assumption of uniaxiality, any other isotropic yield function would reduce to (2.14) in the rate-independent limit, but not, in general, for rate-dependent problems.

↵3 In terms of the Lamé constants and the yield stress, this asymptotic regime corresponds to .

- Received May 4, 2012.
- Accepted August 6, 2012.

- This journal is © 2012 The Royal Society