Scalar evolution equations for shear waves in incompressible solids: a simple derivation of the Z, ZK, KZK and KP equations

Michel Destrade, Alain Goriely, Giuseppe Saccomandi


We study the propagation of two-dimensional finite-amplitude shear waves in a nonlinear pre-strained incompressible solid, and derive several asymptotic amplitude equations in a simple, consistent and rigorous manner. The scalar Zabolotskaya (Z) equation is shown to be the asymptotic limit of the equations of motion for all elastic generalized neo-Hookean solids (with strain energy depending only on the first principal invariant of Cauchy–Green strain). However, we show that the Z equation cannot be a scalar equation for the propagation of two-dimensional shear waves in general elastic materials (with strain energy depending on the first and second principal invariants of strain). Then, we introduce dispersive and dissipative terms to deduce the scalar Kadomtsev–Petviashvili (KP), Zabolotskaya–Khokhlov (ZK) and Khokhlov–Zabolotskaya–Kuznetsov (KZK) equations of incompressible solid mechanics.

1. Introduction

The Zabolotskaya (Z) equation describes the propagation of finite-amplitude, two-dimensional, linearly polarized, shear waves in nonlinear solids. It can be presented in the following form (owing to Cramer & Andrews 2003): Embedded Image 1.1 where U=U(χ,η,τ) is a measure of the amplitude of the wave, χ is a spatial variable measured in a frame moving with the wave, η is a spatial variable transverse to the wave propagation, τ is a temporal variable and β characterizes the material properties of the solid. Zabolotskaya (1986) derived it via a perturbation scheme from the general equations of elastodynamics. It is a counterpart to the Zabolotskaya–Khokhlov (ZK) equation, itself describing the propagation of weakly nonlinear, weakly diffracting, two-dimensional sound beams (Zabolotskaya & Khokhlov 1969), Embedded Image 1.2 (Note that this form is also due to Cramer & Andrews (2003), and that the nonlinear parameter β is not the same in each equation.) The history of the ZK equation and its variants has recently been retraced by Rudenko (2010).

There are two major differences between the Z and ZK equations. One is obviously the cubic term in equation (1.1) versus the quadratic term in equation (1.2). The other is that the ZK equation governs the motion of the longitudinal displacement and is a scalar equation, while the Z equation governs the propagation of plane-polarized shear waves and is thus a vector equation, a priori. Nonetheless, recent papers (Enflo et al. 2006; Wochner et al. 2008; Pinton et al. 2010) have investigated the propagation of linearly polarized nonlinear shear waves in incompressible solids.

Here, we aim at finding what restrictions must be imposed, if at all, on the form of the strain energy W (say), in order to obtain a scalar Z equation. We also include the effects of pre-strain, similar to Cramer & Andrews (2003) in the compressible case. We point out that another problem of interest related to these model equations concerns the mathematical validation of the various approximations. This problem has been solved in the framework of the Navier–Stokes and Euler system by Rozanova-Pierrat (2009). Moreover, because we are mainly interested in understanding the physical underpinnings at the basement of various modelling assumptions, we adopt here a direct method for the derivation of the model equations, not a general reductive perturbation scheme such as the one proposed by Taniuti (1990), for example.

In §3, we establish that a scalar Z equation can indeed be generated for the whole class of the so-called generalized neo-Hookean solids, for which the strain-energy density does not depend on the second principal invariant of the Cauchy–Green strain tensor. We then try to generalize the study to the case of a generic fourth-order incompressible solid, for which W is expanded up to the fourth order in invariants of the Green strain tensor. We show in §3b that there, the Z equation is necessarily a vector equation that is, only plane-polarized—not linearly polarized—two-dimensional shear waves are possible.

When dissipative and dispersive effects are taken into account for the propagation of longitudinal waves, the ZK equation is generalized to Embedded Image 1.3 where β′ and β′′ are other material constants. This scalar equation covers the purely dissipative (at β′′=0) Khokhlov–Zabolotskaya–Kuznetsov (KZK) equation (Kuznetsov 1971) and the purely dispersive Kadomtsev–Petviashvili (KP) equation (Kadomtsev & Petviashvili 1970). In §4, we consider dissipative and dispersive effects in generalized neo-Hookean solids. We derive the counterpart to equation (1.3) for shear waves as Embedded Image 1.4 a scalar equation that covers the Z equation, the modified Khokhlov–Zabolotskaya–Kuznetsov (mKZK) equation and the modified Kadomtsev–Petviashvili (mKP) equation.

The aim of this paper is to provide a simple derivation of these scalar equations for incompressible materials. This quest is largely motivated by the pressing need to model acoustic-wave propagation in soft tissues and gels, where the combined effects of incompressibility, pre-strain, viscosity and dispersion cannot be ignored. This is a topic of fundamental relevance to ultrasound-based techniques of medical imaging such as elastography.

2. Basic equations

Consider a body of incompressible isotropic elastic material, initially at rest in a reference configuration Embedded Image, say. For the study of anti-plane shear motions, we choose a system of Cartesian coordinates (X1,X2,X3) to denote the position of a particle in the body in Embedded Image. We then call (x1,x2,x3) the Cartesian coordinates, aligned with (X1,X2,X3), corresponding to the current position x in the deformed configuration Embedded Image, say. We study the anti-plane motion Embedded Image 2.1 where λ is a positive constant—the axial pre-stretch—and u is a three-times continuously differentiable function—the shear motion.

We call F≡∂x/∂X the associated deformation gradient and B=FFT the left Cauchy–Green strain tensor. The first principal isotropic strain invariant, I1≡trB, is here Embedded Image 2.2 and the second principal isotropic strain invariant, Embedded Image, is equal to Embedded Image 2.3 where u1≡∂u/∂X1 and u2≡∂u/∂X2. In all generality, the strain-energy density W is a function of I1 and I2 only, and the Cauchy stress tensor σ is derived as (e.g. Ogden 1984) Embedded Image 2.4 where W1≡∂W/∂I1, W2≡∂W/∂I2 and p is a Lagrange multiplier introduced by the constraint of incompressibility (which imposes that Embedded Image at all times).

The equations of motion, in the absence of body forces, are Embedded Image 2.5 where ρ is the mass density, which remains constant owing to the incompressibility constraint. Here, the equations of motion reduce to the following system of partial differential equations (see Knowles (1976) for the static case): Embedded Image 2.6 and Embedded Image 2.7 where summation over the repeated index i=1,2 is assumed. The auxiliary quantity q is given by Embedded Image 2.8

3. The Z equation

(a) Results for generalized neo-Hookean materials

The system (2.6) consists of three partial differential equations for the two unknown functions p and u. Therefore, the system is overdetermined and it is not possible, in general, to find a general solution for arbitrary choices of the strain energy W. Progress is possible for some special forms of the strain-energy density, in particular for the Mooney–Rivlin form, which is linear in I1 and I2 (Adkins 1954), but only for static solutions (see Knowles (1976) for details). The dynamic situation for anti-plane motions is much more complex, as discussed by Hayes & Saccomandi (2004) in the Mooney–Rivlin case.

However, if we consider materials for which W=W(I1−3) only—the so-called generalized neo-Hookean materials—then system (2.6) reduces to a well-determined system, Embedded Image 3.1 where W′ is the derivative of W with respect to its argument (I1−3).

(i) Derivation of the Z equation

To derive the Z equation, we assume that W is at least of class Embedded Image. We assume that the primary propagation is in the X1-direction and introduce the following new displacement and new variables: Embedded Image 3.2 where ϵ is a small, non-dimensional parameter, giving a measure of the motion amplitude, the plus or minus sign corresponds to left or right running waves and α is a parameter to be adjusted to cancel all terms to order ϵ. (Note that the scalings used here mirror those of Zabolotskaya (1986), and give a consistent asymptotic expansion.) To obtain the Z equation, we expand equation (3.1), that is Embedded Image 3.3 to order ϵ3 (note that the subscript ‘i’ here denotes derivation with respect to Xi).

First, we compute the partial derivatives as follows: Embedded Image 3.4 which leads to Embedded Image 3.5

Second, we expand W′ in ϵ. Since Embedded Image 3.6 we have Embedded Image 3.7 where w(1)=W′(λ2+2λ−1−3) and w(2)=W′′(λ2+2λ−1−3). Third, we use the expansions for W′ and ui to compute each term in equation (3.3) up to order ϵ3. That is Embedded Image 3.8 Equation (3.3) now reads Embedded Image 3.9 Terms to order ϵ disappear by choosing α2=ρ/(2w(1)), which gives to order ϵ3, Embedded Image 3.10

Finally, taking a derivative of this equation with respect to τ and using the substitution U=vτ leads to the scalar Z equation Embedded Image 3.11 where β=3w(2)α5/ρ.

For example, consider the Yeoh strain energy described by Embedded Image 3.12 where μ>0 is the initial shear modulus and γ>0 is a nonlinear elastic constant. In this case, we find Embedded Image 3.13

Note that the apparent difference between equation (3.11) and the Z equation (1.1) is simply owing to the choice of independent variables. Here, we chose the same independent variables as Norris (1998), where τ is a retarded/advanced time, χ is the spatial variable in the propagation direction and η is the transverse spatial variable.

(ii) Special solutions

Following Cates & Crighton (1990) and Sionóid & Cates (1994), we perform the transformation Embedded Image 3.14 which reduces the Z equation (3.11) to an equation that, once integrated in ζ, yields Embedded Image 3.15 The solution of this one-dimensional equation does not provide a general solution to the two-dimensional equation (3.11); rather, it gives the solution to the particular Cauchy problems where initial data are given on the parabolic surfaces ζ=const.

A further change of function (from G to Embedded Image) and of variables (from χ to Embedded Image), Embedded Image 3.16 reduces equation (3.11) to the autonomous inviscid generalized Burgers equation Embedded Image 3.17 The general integral of this equation is Embedded Image, where F is an arbitrary function. This allows us to recover the following class of exact implicit solution for the Z equation (3.11): Embedded Image 3.18

It is quite remarkable that the transformation (3.14) should work for the Z equation (cubic nonlinearity) because it was initially intended for the KZK equation (quadratic nonlinearity). Through it, a huge class of exact solutions is generated for the Z equation.

General methods to obtain exact solutions for such equations include reduction methods such as the one based on Lie group analysis and its generalizations. Those methods have been systematically applied to the ZK equations with quadratic nonlinearities (e.g. Tajiri 1995; Bruzon et al. 2009), but it seems that there are very few results, regarding the case of a cubic nonlinearity, of interest to incompressible elastic materials.

(b) Results for general incompressible solids

We return to an important point evoked in §1: that the propagation of finite-amplitude transverse waves cannot be governed by a scalar Z equation (i.e. linearly polarized shear waves exist only in special solids, such as the generalized neo-Hookean materials).

To show this, we take solids with the following Rivlin strain energy: Embedded Image 3.19 (where C10, C01 and C20 are constants), which covers the most general model of fourth-order incompressible elasticity (Ogden 1974), Embedded Image 3.20 with the identifications (Destrade et al. 2010) Embedded Image 3.21 Here, E=(FTFI)/2 is the Green strain tensor and μ, A and D are the second-, third- and fourth-order elasticity constants, respectively.

With this strain energy, the first and second equations of (2.6) read Embedded Image 3.22 where Embedded Image and Embedded Image and Embedded Image are defined as in §2.

Guided by the results of the previous subsection, we know that we must perform the changes of function and variables defined by equation (3.2), if we aim at deriving the Z equation from the last equation in (2.6). Here, those changes give the first terms in the expansion of the two equations (3.22) as Embedded Image 3.23 Now the first of these tells us that Embedded Image must be of order ϵ2: Embedded Image, say. Then, equations (3.22) reduce further to Embedded Image 3.24 at order ϵ2. Clearly, these expressions for the derivatives of Q are not compatible.

We may thus conclude that linearly polarized waves are not possible in general incompressible elastic materials. Indeed, Horgan & Saccomandi (2003) explain that when an anti-plane shear deformation occurs in a general incompressible material, it is necessarily accompanied by secondary motions associated with in-plane deformations. This coupling increases dramatically the complexity of the calculations and is an obstruction to the propagation of linearly polarized two-dimensional shear waves. In general incompressible materials, it is necessary to take into account the in-plane motions coupled to the two-dimensional shear wave.

In fact, Wochner et al. (2008) acknowledge as much when they correct Enflo et al. (2006) and their study of linearly polarized nonlinear shear waves, and conclude that ‘Simply stated, quadratic nonlinearity prevents the existence of a pure linearly polarized shear wave beam’. However, they then state that ‘Cubic nonlinearity can be more pronounced than quadratic nonlinearity due to the quasi-planar nature of wavefronts in directional sound beams’, and go on to ‘consider only the effects of cubic nonlinearity’. But, ignoring quadratic terms in favour of cubic terms runs contrary to the very idea of asymptotic expansions. Upon examination of their equation of motion (eqn (11)), we see that the only way to make quadratic terms vanish is to take 4μ+A=0. According to equation (3.21) here, this is equivalent to saying that C01=0, making W in equation (3.20) belong to the class of generalized neo-Hookean materials.

4. Dissipation and dispersion

In our treatment of dissipative and dispersive effects, we use the following constitutive law for the Cauchy stress: Embedded Image 4.1 where Embedded Image 4.2 which separates the purely elastic part σE of the stress from a viscoelastic part σD.

More specifically, σE coincides with the Cauchy stress of an elastic Yeoh material with strain energy as in equation (3.12). We make this choice for convenience, but it is clear that once the asymptotic expansions are performed below, the Yeoh form recovers the entire class of generalized neo-Hookean solids, through the identifications Embedded Image 4.3

The viscoelastic part of the stress, σD, is the stress of a second-grade viscoelastic solid, where ν>0 and α1>0 are constants and A1 and A2 are the first two Rivlin–Ericksen tensors, defined by Embedded Image 4.4 with Embedded Image, the displacement gradient. Here, the superposed dot represents the material time derivative. The first part of σD is νA1=2νD, where D is the stretch tensor; that part describes dissipation by Newtonian viscosity effects, and coincides with the dissipative term in the Navier–Stokes theory. The second part of σD is Embedded Image; it describes dispersion by accounting for the intrinsic characteristic lengths of the solid (see Destrade & Saccomandi (2006) for details).

Note that the temporal evolution of the canonical energy E(t) for our solid (4.1) and (4.2) is given by Embedded Image, where Ωt is a given region occupied by the material at time t. It is independent of α1 (Fosdick & Yu 1996).

With respect to two-dimensional anti-plane motion, the derivation of the governing equation is not overly complicated by the additional terms. For ease of exposition, we consider that the solid is not pre-strained: λ=1. (It follows that once the asymptotic expansions are performed, the identifications (4.3) reduce to μ=2w(1), γ=2w(2)/w(1).)

We find that the equations qα=0(α=1,2) still apply here, while the last equation in (3.1) is changed to Embedded Image 4.5 In dimensionless form, it becomes Embedded Image 4.6 using the usual scalings.

When the solid is not dissipative (ν=0) and not dispersive (α1=0), the changes of function and variables (3.2) give a governing equation at order ϵ3, when α is chosen as Embedded Image. Here, the expansions of the dissipative and dispersive terms in equation (4.6) start as Embedded Image 4.7

Clearly, in order to be incorporated into the governing equation for the motion, the dissipative and the dispersive effects must be at least of the order Embedded Image and Embedded Image say, where Embedded Image and Embedded Image are constants. With that assumption, we find the following form of the equation of motion for velocity Embedded Image, Embedded Image 4.8 where Embedded Image 4.9 This equation covers not only the Z equation (at β′=β′′=0), but also the dissipative mZK equation (at β′′=0) and the dispersive mKP equation (at β′=0).

To the best of our knowledge, this is the first time that an mKP equation is derived in the continuum mechanics of homogeneous solids.

If we consider the diffusive version of equation (4.8), i.e. Embedded Image 4.10 and we seek solutions that propagate at an angle γ with respect to the τ-direction, it is possible to obtain an interesting reduction. To this end, let us consider solutions that depend only on χ and Embedded Image to reduce equation (4.10) to Embedded Image 4.11 This is the Nariboli–Lin-modified diffusive Burgers equation (Nariboli & Lin 1973). The plane-wave solution of this equation is computed as Embedded Image where V 0 is the speed of the plane wave. This solution is quite different from the plane-wave solution of the classical Burgers equation, which is expressed in terms of the hyperbolic tangent function.

If we consider the purely dispersive case and we rescale some of the variables and Embedded Image and Embedded Image, it is always possible to choose the scaling such that Embedded Image 4.12 the mKP equation. A connection between this equation and the modified Korteweg de Vries equation may be established with a similar transformation to the one used to reduce equation (4.10) to equation (4.11).

It is also possible to seek direct solutions of equation (4.12) in the form Embedded Image This ansatz produces the reduction Embedded Image 4.13 where c0 and c1 are integration constants.

If we consider solitary-wave solutions, the asymptotic boundary conditions require c0=c1=0 and the usual energy integral is given by Embedded Image where again c2 is an integration constant. It is clear that it must be c2=0 and if b2>0, an explicit solitary-wave solution is given by Embedded Image where ζ0 is an integration constant.

If we consider in equation (4.13) c0=0 and c1≠0, it is also possible to find cnoidal wave solutions in terms of elliptic functions.

5. Conclusions

We provide a clear and complete derivation of various model equations for shear waves in the framework of nonlinear incompressible elasticity. The rationale for investigating this type of problem is dictated by the use of imaging techniques of improved quality, such as latest generation elastography, to evaluate the nonlinear material properties of soft tissues.

We show that it is possible to derive a pure (scalar) Z equation only in the special case of generalized neo-Hookean materials. For general incompressible materials, the Z equation is coupled with a non-negligible transverse motion, an important issue that is often overlooked in the literature. Clearly, a rigorous evaluation of the validity of the perturbative expansion requires an extension of the existence, uniqueness and stability results that have been provided so far only for model equations with quadratic nonlinearity.

We point out that general structure of the various equations obtained here is as follows: Embedded Image where c is a suitable material parameter, is the second-order derivative (Laplace operator) in an orthogonal direction to the beam axis and Embedded Image is the equation governing plane-wave propagation (inviscid Burgers, Burgers or Korteweg de Vries equation). It would be an interesting problem to study the various transformation and group-analysis properties of such general equations to find exact and explicit solutions for various model equations following a general and complete scheme.


This work is supported by a Senior Marie Curie Fellowship awarded by the Seventh Framework Programme of the European Commission to the first author. This publication is based on work supported in part by award no. KUK-C1-013-04 , made by King Abdullah University of Science and Technology (KAUST) (A.G.), and also based in part upon work supported by the National Science Foundation under grant DMS-0907773 (A.G.).

  • Received October 2, 2010.
  • Accepted November 12, 2010.


View Abstract