Interface Waves
in Pre-Stressed Incompressible Solids
Michel Destrade
Institut Jean Le Rond d’Alembert, CNRS/Université Pierre et Marie Curie, Paris, France
Abstract We study incremental wave propagation for what is seemingly the simplest boundary value problem, namely that constitued by the plane interface of a
semi-infinite solid. With a view to model loaded elastomers and soft tissues, we
focus on incompressible solids, subjected to large homogeneous static deformations.
The resulting strain-induced anisotropy complicates matters for the incremental
boundary value problem, but we transpose and take advantage of powerful techniques and results from the linear anisotropic elastodynamics theory. In particular
we cover several situations where fully explicit secular equations can be derived,
including Rayleigh and Stoneley waves in principal directions, and Rayleigh waves
polarized in a principal plane or propagating in any direction in a principal plane.
We also discuss the merits of polynomial secular equations with respect to more
robust, but less transparent, exact secular equations.
1 Introduction
The term “acousto-elastic effect” describes the interplay between the static deformation
of an elastic solid and the motion of an elastic wave. If both the deformation and the
motion are of infinitesimal amplitude, then all the governing equations are linearized, see
for instance the Chapter by Norris for examples and applications or the experimental
results of Pao et al. (1984). If both the deformation and the motion are of finite amplitude, then the resulting governing equations are highly nonlinear, and their resolution is
the subject of much research, see the Chapter by Fu for the weakly nonlinear theory and
the Chapter by Saccomandi for the fully nonlinear theory.
In between those two situations lies the theory of “small-on-large”, also known as the
theory of “incremental” motions, where the wave is an infinitesimal perturbation superimposed onto the large static homogeneous deformation of a generic hyperelastic solid.
There, the homogeneous character of the static deformation and the linear character of
the incremental equations of motion ensure that the calculations are valid for any strain
energy density (to be specified later for applications, if necessary). The next Section
of this Chapter briefly recalls the governing equations of incremental motions (see the
Chapter by Ogden for their derivation).
It turns out that many similarities can be drawn between the equations of incremental
motions and those of linear anisotropic elasticity, with the main difference that in the
latter case, the anisotropy is set once and for all for a given crystal whereas in the former
1
case, it is strain-induced and susceptible to great variations from one configuration to
another. Using the similarities, we may transpose the so-called Stroh formulation and
exploit its many results; on the other hand, when focussing on the differences, we may
highlight the influences of the pre-stress and of the choice of a strain-energy density
on the propagation of waves. In this Chapter, attention is restricted to waves at the
interface of pre-deformed, semi-infinite solids, in contact either with vacuum (Rayleigh
waves) or with another solid (Stoneley waves). With a view to model elastomers and
biological soft tissues, the solids are considered to be incompressible (mathematically,
this internal constraint lightens somewhat the expressions but does not prove essential
to the resolution).
Several situations are treated: principal wave propagation in Section 3, principal polarization in Section 4, and principal plane propagation in Section 5. The emphasis is
on deriving explicit secular equations in polynomial form, using some simple “fundamental equations” derived at the end of Section 2. Of course, as the setting gets more
and more involved, so does the search for a polynomial secular equation; eventually its
degree becomes too high for comfort and other techniques are required. The concluding
section (Section 6) discusses the pros and cons of such equations, as opposed to exact,
non-explicit, secular equations, free of spurious roots.
2 Basic equations
2.1
Finite deformation
Consider an isotropic, incompressible, hyperelastic solid at rest, characterized by a
mass density ρ and a strain energy function W . Then subject it to a large, static, homogeneous deformation (“the pre-strain”) carrying the particle at X in the undeformed
configuration to the position x in the deformed configuration.
Call F = ∂x/∂X the corresponding constant deformation gradient and B = F F t the
associated left Cauchy-Green strain tensor. This tensor being symmetric, the directions
of its eigenvectors are orthogonal; they are called the principal axes of pre-strain or in
short, the principal axes. Also, the eigenvalues of B are positive, λ21 , λ22 , λ23 , say, and
λ1 , λ2 , λ3 are called the principal stretches. Figure 1 shows how a unit cube with edges
aligned with the principal axes, is transformed by the pre-strain. Note that because the
solid is incompressible, its volume is preserved through any deformation so that here,
λ1 λ2 λ3 = 1.
(2.1)
The first two principal invariants of strain are defined as
I1 = tr B,
I2 = [(tr B)2 − tr (B 2 )]/2.
(2.2)
In the Cartesian coordinate system aligned with the principal axes, B is diagonal. Calling
e1 , e2 , e3 , the unit vectors in the x1 , x2 , x3 directions, respectively, we have
B = λ21 e1 ⊗ e1 + λ22 e2 ⊗ e2 + λ23 e3 ⊗ e3 ,
(2.3)
and the computation of I1 , I2 there gives
I1 = λ21 + λ22 + λ23 ,
I2 = λ21 λ22 + λ22 λ23 + λ23 λ21 .
2
(2.4)
Figure 1. Finite homogeneous deformation of a unit cube.
For an isotropic solid, W may be given as a function of the invariants: W = W (I1 , I2 )
or equivalently, as a symmetric function of the principal stretches: W = W (λ1 , λ2 , λ3 ),
according to what is most convenient for the analysis or according to how W has been
determined experimentally. With the first choice, the constant Cauchy stress (“the prestress”) necessary to maintain the solid in its state of finite homogeneous deformation
is
σ = −pI + 2(∂W/∂I1 + I1 ∂W/∂I2 )B − 2(∂W/∂I2 )B 2 ,
(2.5)
where p is a Lagrange multiplier due to the constraint of incompressibility (a yet arbitrary
constant scalar to be determined from initial and boundary conditions.) With the second
choice, the non-zero components of σ relative to the principal axes are written as
σi = −p + λi ∂W/∂λi ,
i = 1, 2, 3 (no sum).
(2.6)
The proof for the equivalence between (2.5) and (2.6) relies on the connections (2.4).
2.2
Incremental equations
Consider a half-space filled with an incompressible hyperelastic solid subject to a large
homogeneous deformation. We take the Cartesian coordinate system (x̂ 1 , x̂2 , x̂3 ) to be
oriented so that the boundary is at x̂2 = 0, and we study the propagation in the x̂1
direction of an infinitesimal interface wave in the solid.
This wave is inhomogeneous as it progresses in an harmonic manner in a direction
lying in the interface, while its amplitude decays with distance from the boundary.
We call u the mechanical displacement associated with the wave, and ṗ the increment
in the Lagrange multiplier p due to incompressibility. The incremental nominal stress
tensor s has components
sji = A0jilk uk,l + puj,i − ṗδij ,
3
(2.7)
where the comma denotes partial differentiation with respect to the coordinates x̂ 1 , x̂2 ,
x̂3 . Here, A0 is the fourth-order tensor of instantaneous elastic moduli, with components
A0jilk = Fjα Flβ
∂2W
= A0lkji .
Fiα Fkβ
(2.8)
Note that due to the symmetry above, A0 has in general 45 independent components.
In the principal axes coordinate system (x1 x2 x3 ) however, there are only 15 independent
non-zero components; they are (Ogden, 2001):
A0iijj = λi λj Wij ,
A0ijij = (λi Wi − λj Wj )λ2i /(λ2i − λ2j ),
i 6= j,
A0ijij = (A0iiii − A0iijj + λi Wi )/2,
i 6= j,
A0ijji = A0jiij = A0ijij − λi Wi ,
i 6= j,
λi 6= λj ,
λ i = λj ,
(2.9)
(no sums on repeated indexes here), where Wj = ∂W/∂λj and Wij = ∂ 2 W/(∂λi ∂λj ).
Finally, the governing equations are the incremental equations of motion and the
incremental constraint of incompressibility; they read
sji,j = ρ∂ 2 ui /∂t2 ,
uj,j = 0,
(2.10)
respectively.
Now everything is in place to solve an interface wave problem. We take u and ṗ in
the form
{u, ṗ} = {U (kx̂2 ), ikP (kx̂2 )}eik(n·x̂−vt) ,
(2.11)
where k is the wave number, U and P are functions of the variable kx̂2 only, n is the
unit vector in the direction of propagation, and v is the speed. Clearly by (2.7), s has a
similar form, say
s = ikS(kx̂2 )eik(n·x̂−vt) ,
(2.12)
where S is a function of kx̂2 only.
After substitution, it turns out that the incremental governing equations (2.10) can
be cast as the following first-order differential system (Chadwick, 1997),
ξ 0 = iN ξ,
where ξ = [U , t]t ,
(2.13)
the prime denotes differentiation with respect to the variable kx̂2 , and t are the tractions
acting on planes parallel to the boundary, with components tj = S2j . Here the matrix
N has the following block structure
N1
N2
,
(2.14)
N=
N 3 + ρv 2 I N t1
where N 1 , N 2 = N t2 , and N 3 = N t3 are square matrices. This is the so-called Stroh formulation. In effect, many of the results established thanks to the Stroh (1962) formalism
in linear anisotropic elasticity can formally be carried over to the context of incremental
dynamics in nonlinear elasticity, as shown by Chadwick and Jarvis (1979a), Chadwick
(1997), and Fu (2005a,b).
4
2.3
Resolution
The solution to the first-order differential system (2.13) is an exponential function in
kx̂2 ,
ξ(kx̂2 ) = eikqx̂2 ζ,
(2.15)
where ζ is a constant vector and q is a scalar. Then the following eigenvalue problem
emerges: N ζ = qζ. Its resolution is in two steps.
First, find the eigenvalues by solving the propagation condition,
det (N − qI) = 0,
(2.16)
for q, and keep those qj ’s which satisfy the decay condition. For instance, when the solid
fills up the x̂2 > 0 half-space, the decay condition is
=(q) > 0,
(2.17)
ensuring that the solution (2.15) is localized near the interface and vanishes away from
it. The penetration depth of the interface wave is clearly related to the magnitude of
=(q): the smaller this quantity is, the deeper the wave penetrates into the solid.
The propagation condition (2.16) is a polynomial in q with real coefficients and it
has only complex roots (Fu, 2005b), which come therefore in pairs of complex conjugate
quantities. Hence, half of all the roots to the propagation condition qualify as satisfying
the decay condition. Let ζ j be the eigenvector corresponding to the qualifying root qj .
Now proceed to the second step, which is to construct the general localized solution
to the equations of motion, as
P
ξ(kx̂2 ) = γj eikqj x̂2 ζ j ,
(2.18)
for some arbitrary constants γj . Then compute this vector at the interface x̂2 = 0 and
apply the boundary conditions. The vector ξ(0) is often decomposed as follows,
P
A
U (0)
γ,
(2.19)
= γj ζ j =
ξ(0) =
B
t(0)
where A and B are square matrices and γ is the vector with components γj . For instance,
the archetype of interface waves is the Rayleigh (1885) surface wave, which propagates
at the interface between a solid half-space and the vacuum, leaving the boundary free of
tractions. Mathematically, the corresponding boundary condition is that t(0) = Bγ = 0,
leading to
det B = 0.
(2.20)
This (complex) form of the secular equation is however not the optimal form, and it might
lead to unsatisfactory answers to the questions of existence and uniqueness of the wave
(see Barnett (2000) for an historical account of this point). From the Stroh formalism,
and its application to the present context, we learn that it is much more efficient to
work with the surface impedance matrix (Ingebrigsten and Tonning, 1969) than with the
matrices A and B; this matrix M is defined by
M = −iBA−1 .
5
(2.21)
It is Hermitian (Barnett and Lothe, 1985; Fu, 2005b) and so det M = −i(det B)/(det A)
is a real quantity and the secular equation for Rayleigh surface waves, written in the form
det M = 0,
(2.22)
is a real equation, in contrast to (2.20). Moreover, if there is a root to this equation in
the subsonic regime (where v is less than the speed of any bulk wave), then it is unique;
also, the existence of a root is equivalent to the existence of a surface wave.
Similar results also exist for other types of interface waves as seen in the course of this
Chapter. For instance, the boundary conditions for Stoneley (1924) interface waves are
that displacements and tractions are continuous across the boundary between two rigidly
bonded semi-infinite solids. Then ξ(0) = ξ ∗ (0) where the asterisk refers to quantities
for the solid in x̂2 6 0. Equivalently, Aγ = A∗ γ ∗ , Bγ = B ∗ γ ∗ , from which comes
[BA−1 A∗ − B ∗ ]γ ∗ = 0, leading to
det(M + M ∗ ) = 0,
(2.23)
the optimal form of the secular equation for Stoneley interface waves. Here M ∗ is the
surface impedance matrix for the solid in the x̂2 6 0 half-space, defined as
M ∗ = iB ∗ (A∗ )−1 .
2.4
(2.24)
Explicit secular equations
The derivation of a secular equation, preferably in the optimal form involving the
surface impedance matrix, is no sinecure in general. The problematic step lies in the
resolution of the propagation condition (2.16).
For principal wave propagation (x̂1 , x̂2 , x̂3 are aligned with the principal axes), the
propagation condition factorizes into the product of a term linear in q 2 and a term
quadratic in q 2 . Here we can compute the roots explicitly, keep the qualifying ones (see
(2.17)), and solve the boundary value problem in its entirety. Many problems falling in
this category have been solved over the years, and some are presented in Section 3.
For a non-principal wave with propagation direction and attenuation direction both
in a principal plane (the saggital plane (x̂1 x̂2 ) is a principal plane but x̂1 is not a principal
axis), the propagation condition factorizes into the product of a term linear in q 2 and a
term quartic in q. We treat this case in Section 4. Although it is possible to write down
formally the qualifying roots of the quartic (Fu, 2005a; Destrade and Fu, 2006; Fu and
Brookes, 2006), the formulas involved are cumbersome to interpret.
For a wave propagating in a principal plane but not in a principal direction (x̂2 is
aligned with a principal axis but neither x̂1 nor x̂3 are aligned with principal axes), the
propagation condition is a cubic in q 2 . We treat this case in Section 5. Now it is a
daunting task to find analytical expressions for the roots q satisfying the decay condition
(2.17).
Finally, for wave propagation in any other case, the propagation condition is a sextic
in q, unsolvable analytically according to Galois theory.
These observations suggest that, except in the case of principal waves, numerical
procedures are required in order to make progress. It is indeed the case that sophisticated
6
tools and efficient numerical recipes have been developed by Barnett and Lothe (1985),
Fu and Mielke (2002), and several others, with most satisfying results. However it is
also the case that some interface wave problems can be solved analytically, up to the
derivation of the secular equation in explicit polynomial form. The first steps in that
direction were taken by Currie (1979), and his advances were later refined by Taylor and
Currie (1981) and Taziev (1989), revisited by Mozhaev (1995) and by Ting (2004), and
extended by Destrade (2003).
The equations that turn out to be fundamental in the derivation of explicit polynomial
secular equations are
ˆ n ξ(0) = 0,
ξ(0) · IN
(2.25)
0
I
and n is an integer. Their derivation is most simple. First, it can be
where Iˆ =
I 0
shown by induction (Ting, 2004) that N n has a block structure similar to that of N ,
that is
#
"
(n)
(n)
N1
N2
n
(2.26)
N =
(n)t ,
K (n) N 1
(n)
with K (n) = K (n)t , N 2
(n)t
= N2
. It then follows that
#
"
(n)t
(n)
n
K
N
1
ˆ
IN
=
(n)
(n)
N1
N2
(2.27)
is symmetric for all n. Now take the scalar product of both sides of the governing
ˆ n ξ to get
equation (2.13) by IN
ˆ n ξ = iξ · IN
ˆ n+1 ξ;
ξ0 · IN
(2.28)
finally add its complex conjugate to this equality to end up with
ˆ n ξ + ξ · IN
ˆ n ξ 0 = 0,
ξ 0 · IN
(2.29)
and, by integration between the interface (at x̂2 = 0) and infinity (where U and t, and
thus ξ, vanish), arrive at (2.25).
For instance, the boundary condition for Rayleigh surface waves is that there are
no incremental tractions at the interface; thus ξ(0) = [U (0), 0]t , and the fundamental
equations (2.25) reduce to
U (0) · K (n) U (0) = 0.
(2.30)
3 Principal waves
Here we take (x̂1 , x̂2 , x̂3 ) to coincide with the principal axes (x1 , x2 , x3 ).
deformation is thus
x̂1 = λ1 X1 , x̂2 = λ2 X2 , x̂3 = λ3 X3 .
The pre(3.1)
Figure 2 summarizes the situation with respect to the waves’ characteristics near the
interface. Bear in mind that the wave analysis is linear and gives no indication about
the amplitude; moreover, a half-space has no characteristic length so that the secular
equation is non-dispersive and the wavelength remains undetermined.
7
x3
x1
x2
propagation
attenuation
Figure 2. Incremental wave propagation localized near the surface of a semi-infinite
deformed solid. The analysis does not give the amplitude nor the wavelength.
3.1
Governing equations
For principal waves, the fields (2.11) and (2.12) are independent of x̂3 = x3 , because
x̂2 = x2 and n · x̂ = x̂1 = x1 . Also, recall from (2.9) that the non-zero components of
A0 in the (x1 , x2 , x3 ) coordinate system of principal axes are
A01111 = λ21 W11 ,
λ−2
1 A01212
A01122 = λ1 λ2 W12 , A02222 = λ22 W22 ,
λ1 W 1 − λ 2 W 2
λ2 W 1 − λ 1 W 2
= λ−2
, A01221 =
λ1 λ2 ,
2 A02121 =
2
2
λ1 − λ 2
λ21 − λ22
(3.2)
and also A03333 , A01133 , A02233 , A01313 , A02323 , A03131 , A03232 , A01331 , and A02332 ,
whose expressions are not needed in this Section.
From these observations follows that the third equation of motion (2.10)3 reduces to
0
−S13 + iS23
= −ρv 2 U3 ,
(3.3)
where by (2.7),
iS23 = A02323 U30 .
iS13 = iA01313 U3 ,
(3.4)
Hence the movement along the x3 principal axis is governed by an equation which depends
only on U3 . For this equation, governing what is termed the anti-plane motion, we
take the trivial solution: U3 = 0, and we focus on the in-plane motion. According to
(2.10)1,2,4 , it is governed by
0
−S11 + iS21
= −ρv 2 U1 ,
0
−S12 + iS22
= −ρv 2 U2 ,
8
iU1 + U20 = 0,
(3.5)
where by (2.7),
iS11 = i(A01111 + p)U1 + A01122 U20 − P,
iS21 = A02121 U10 + i(A01221 + p)U2 ,
iS12 = (A01221 + p)U10 + iA01212 U2 ,
iS22 = iA01122 U1 + (A02222 + p)U20 − P.
(3.6)
We eliminate p in favour of the pre-stress: by (2.6) at j = 2, we have p = λ2 W2 − σ2 and
so by (3.2),
A01221 + p = A02121 − σ2 .
(3.7)
It then follows from the second equation above that
1
A02121 − σ2
U2 +
S21 ,
U10 = i −
A02121
A02121
(3.8)
and this constitutes the first line of the first-order system (2.13). The second line comes
from the incremental incompressibility constraint (3.5)3 as
U20 = i [−U1 ] .
(3.9)
0
0
Proceeding similarly for S21
, S22
, we find eventually that the governing equations are
indeed in the form (2.13), where ξ = [U1 , U2 , S21 , S22 ]t and −N 1 , N 2 , and −N 3 are
given by
γ21 − σ2
1
2(β12 + γ21 − σ2 )
0
0
0
(γ21 − σ2 )2 ,
(3.10)
,
γ21 , γ21
0
γ
−
12
1
0
0
0
γ21
respectively, where we used the following short-hand notations (no sums),
γij = A0ijij = λ2i λ−2
j γji ,
2βij = A0iiii + A0jjjj − 2A0iijj − 2A0ijji = 2βji .
(3.11)
or equivalently,
γij = (λi Wi − λj Wj )λ2i /(λ2i − λ2j ) = λ2i λ−2
j γji ,
2βij = λ2i Wii − 2λi λj Wij + λ2j Wjj + 2(λi Wj − λj Wi )λi λj /(λ2i − λ2j ) = 2βji .
3.2
(3.12)
Resolution
The propagation condition (2.16) reduces to a quadratic in q 2 ,
γ21 q 4 + (2β12 − ρv 2 )q 2 + γ12 − ρv 2 = 0.
Notice how σ2 , though present in N , does not appear explicitly in this equation.
9
(3.13)
Calling q12 , q22 , the roots of the quadratic, we have
q12 q22 =
γ12 − ρv 2
,
γ21
q12 + q22 = −
2β12 − ρv 2
.
γ21
(3.14)
The roots q1 , q2 of the biquadratic satisfying the decay condition (2.17) are in one of
the two following forms; either: q1 = iβ1 , q2 = iβ2 , where β1 > 0, β2 > 0, or: q1 = α + iβ,
q2 = −α + iβ, where β > 0. Whatever the case, q12 q22 > 0, q1 q2 < 0, and q1 + q2 is
a purely imaginary quantity. From the first inequality we deduce that (Dowaikh and
Ogden, 1990)
s
η=
γ12 − ρv 2
γ21
(3.15)
is a real quantity. From the second, and using the definitions of η and γij , we find
q1 q2 = −η,
(q1 + q2 )2 =
γ12 − 2β12
β12
− 2η − η 2 = λ21 λ−2
− 2η − η 2 .
2 −2
γ21
γ21
(3.16)
We compute the eigenvectors ζ 1 and ζ 2 of N corresponding to q1 and q2 as any
column of the matrix adjoint to N − q1 I and to N − q2 I, respectively. Choosing the
third column, we find
1
2
a
a
1
2
ζ = 1 , ζ = 2 ,
(3.17)
b
b
where
"
qj2 qj
a = −
,
γ21 γ21
j
#t
,
t
bj = −qj (qj2 − 1 + σ 2 ), qj2 (1 − σ 2 ) − η 2 ,
(3.18)
and σ 2 = σ2 /γ21 is a non-dimensional measure of the pre-stress.
We can now construct the A = [a1 |a2 ] and B = [b1 |b2 ] matrices, and the surface
impedance matrix M = −iBA−1 . It turns out to be
q1 + q 2
1 − σ2 − η
M = −iγ21
,
(3.19)
−(1 − σ 2 − η) (q1 + q2 )η
which is indeed Hermitian because q1 + q2 is a purely imaginary quantity and η is real.
For Rayleigh surface waves, the secular equation is (2.22), or here, using (3.16),
η 3 + η 2 + (2 − λ21 λ−2
2 +2
β12
− 2σ2 )η − (1 − σ 2 )2 = 0.
γ21
(3.20)
Dowaikh and Ogden (1990) established this form of the secular equation for principal
surface waves in pre-stressed incompressible solids, following other works by Hayes and
Rivlin (1961), Flavin (1963), Willson (1973a,b), Chadwick and Jarvis (1979a), Guz (2002,
review with an extensive bibliography), and many others. It is of course consistent with
Lord Rayleigh’s own analysis of surface waves in linear isotropic incompressible solids.
To check this, let the solid be un-stressed (σi = 0) and un-deformed (λi = 1); then η
10
p
reduces to 1 − ρv 2 /µ0 , where µ0 is the infinitesimal shear modulus; also, β12 = γ21
and η is the real root of η 3 + η 2 + 3η − 1 = 0, that is η ' 0.2956 giving ρv 2 /µ0 ' 0.9126,
as found by Rayleigh (1885).
For Stoneley interface waves, the secular equation is (2.23) where
∗
∗
γ (q + q2 ) − γ21
(q1∗ + q2∗ )
γ21 (1 − η) − γ21
(1 − η ∗ )
M + M ∗ = −i 21 1
. (3.21)
∗
∗
−γ21 (1 − η) + γ21
(1 − η ∗ ) γ21 (q1 + q2 )η − γ21
(q1∗ + q2∗ )η ∗
This equation was studied in great detail by Dowaikh and Ogden (1991) and by Chadwick
(1995). It is consistent with the analysis of Stoneley (1924) of interface waves in linear
isotropic incompressible solids. A remarkable feature of this secular equation for principal
Stoneley interface waves in deformed incompressible solids – first noted by Chadwick and
Jarvis (1979b) – is that the pre-stress σ2 does not appear explicitly in it, in contrast to
the equation for surface waves (3.20). This quantity, which is continuous across the
interface (σ2 = σ2∗ ), disappears in the addition of the two surface impedance matrices.
Of course it still plays an implicit role, in determining the pre-strain.
Dowaikh and Ogden (1990, 1991), Chadwick (1995), and Guz (2002, review) have
covered almost every aspect of principal interface wave propagation and more information
can be found in their respective articles. In the next Subsection we rapidly work out two
examples of surface waves.
3.3
Examples
First we present an example taken from the literature on elastomers, where the
Mooney-Rivlin strain energy function is often encountered. It is given by
W = D1 (λ21 + λ22 + λ23 − 3)/2 + D2 (λ21 λ22 + λ22 λ23 + λ23 λ21 − 3)/2,
(3.22)
where D1 and D2 are positive constants with the dimensions of a stiffness. The MooneyRivlin material enjoys special properties with respect to wave propagation (the neoHookean material, which corresponds to the special case D2 = 0, enjoys even more
special properties as is seen in Section 5.3). For instance, once subjected to a large
homogeneous pre-strain, it permits the propagation of bulk waves in every direction;
these waves can be infinitesimal, but also of arbitrary finite amplitude (Boulanger and
Hayes, 1992); they can be homogeneous plane waves but also inhomogeneous plane waves
(Destrade, 2000, 2002). The quantities (3.12) are also quite special; they are
γij = (D1 + D2 λ2k )λ2i ,
2βij = (D1 + D2 λ2k )(λ2i + λ2j ),
(3.23)
where k 6= i, j, and thus they satisfy
2βij = γij + γji .
(3.24)
These relationships mean that the biquadratic (3.13) factorizes to
(q 2 + 1)(q 2 + η 2 ) = 0,
(3.25)
and that the secular equation (3.20) reduces to
η 3 + η 2 + (3 − 2σ 2 )η − (1 − σ 2 )2 = 0.
11
(3.26)
Hence one qualifying root is q1 = i, whatever the values of the material constants D1
and D2 . The other root is q2 = iη. When there is no pre-stress normal to the boundary
(σ 2 = 0), then η is the real root of η 3 + η 2 + 3η − 1 = 0, that is η ' 0.2956 giving
ρv 2 = γ12 − γ21 η 2 = (D1 + D2 λ23 )(λ21 − 0.0874λ22),
(3.27)
a result first established by Flavin (1963). Here q1 = i and q2 ' 0.295i so that the
penetration depth of the surface wave is fixed and is completely independent of the prestrain and of the material parameters D1 and D2 . We say that the penetration depth is
universal relative to the class of Mooney-Rivlin materials.
The second example is taken from the biomechanics literature. From a series of
uniaxial tests on human aortic aneurysms, Raghavan and Vorp (2000) deduced that the
following strain energy density gave a satisfying fit with the data plots,
W = C1 (λ21 + λ22 + λ23 − 3) + C2 (λ21 + λ22 + λ23 − 3)2 ,
(3.28)
where, typically, C1 = 0.175 MPa, C2 = 1.9 MPa (Karduna et al. (1997) use the same
expression to model the response of passive myocardium.) A uniaxial pre-stress is σ 1 6= 0,
σ2 = σ3 = 0, leading through (2.6) to the following equi-biaxial pre-strain,
λ1 = λ,
λ2 = λ−1/2 ,
λ3 = λ−1/2 ,
(3.29)
where λ is calculated from
σ1 = 2(λ2 − λ−1 )[C1 + 2C2 (2 − 3λ−1 + λ−3 )].
(3.30)
Finally, using the following expressions for the relevant moduli,
γ21 = 2C1 λ−1 + 4C2 (2λ−1 − 3λ−2 + λ−4 ),
γ12 = λ3 γ21 ,
β21 = C1 (λ2 + λ−1 ) + 2C2 (4λ2 − 3λ − λ−1 − 3λ−2 + 3λ−4 ),
(3.31)
it is a simple matter to solve the secular equation (3.20) numerically and plot the variations of the squared wave speed, scaled with respect to the squared bulk wave speed
γ12 /ρ, with the pre-stretch λ. Figure 3b displays these variations; for comparison purposes, Figure 3a shows the variations of the scaled squared wave speed in the case of a
Mooney-Rivlin material in uniaxial stress; in that later case the graph is independent of
the material parameters D1 and D2 because by (3.27), ρv 2 /γ21 = 1 − 0.0874λ−3. The
dashed lines indicate the speed of Lord Rayleigh’s squared speed in the isotropic (no
pre-strain) case where λ = 1, ρv 2 /γ12 ' 0.9126. Notice how different the responses of
two solids are in that neighbourhood. Note also that for high compressive stretches, the
squared speeds eventually falls off to zero; this happens at λ ' 0.08741/3 ' 0.444 for all
Mooney-Rivlin materials, as shown by Biot (1963), and at λ ' 0.315 for the soft biological tissue model above. Beyond that critical compression stretch, v 2 < 0, leading to a
purely imaginary v, an amplitude which then grows exponentially with time according
to (2.11), and a breakdown of the linearized analysis. The search for critical compression
stretches is an extremely active area of research, clearly linked to the geometric stability
analysis of solids.
12
ρv2/γ1 2
ρv2/γ1 2
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
λ
0
0.2
0.4
0.6
0.8
1
1.2
λ
0
1.4
0.2
0.4
0.6
0.8
1
1.2
1.4
Figure 3. Variations of the scaled squared surface wave speed with the stretch in uniaxial
pre-stress (a) for any Mooney-Rivlin material and (b) for the solid with strain energy
density (3.28) where C1 = 0.175 MPa, C2 = 1.9 MPa.
4 Waves polarized in a principal plane
In this Section we study the case where x̂3 is aligned with the principal axis x3 but neither
x̂1 (propagation direction) nor x̂2 (attenuation direction) are aligned with principal axes.
The components Â0jilk in the (x̂1 x̂2 x̂3 ) coordinate system of the instantaneous elastic
moduli tensor A0 are related to the components A0jilk , given by (3.2), in the (x1 x2 x3 )
coordinate system of principal axes through the tensor transformations
cos Θ − sin Θ 0
Â0jilk = Ωjp Ωiq Ωlr Ωks A0pqrs , where Ω = sin Θ cos Θ 0 ,
(4.1)
0
0
1
and Θ is the angle between x1 and x̂1 .
In particular we find that the non-zero components in the forms Â013lk and Â023lk
are Â01313 , Â01323 , Â01331 , Â02323 , Â02313 , and Â02332 . The mechanical fields (2.11) and
(2.12) are independent of x̂3 = x3 because n · x̂ = x̂1 here. As a consequence, the third
equation of motion (2.10)3 reduces to
0
−S13 + iS23
= −ρv 2 U3 ,
(4.2)
where by (2.7),
iS13 = iÂ01313 U3 + Â01323 U30 ,
iS23 = iÂ02313 U3 + Â02323 U30 .
(4.3)
Hence the movement along the x3 principal axis is governed by an equation which depends
only on U3 , and for this anti-plane motion we take the trivial solution: U3 = 0.
The equations governing the in-plane motion have been derived in the case of a
general plane pre-strain by Fu (2005a) and solved for surface waves in the case of a prestrain consisting in a triaxial stretch followed by a simple shear by Destrade and Ogden
(2005). Instead of treating these cases again, we revisit the case relative to one of the
13
most important pre-strain fitting into the present context, that of finite simple shear,
presented originally by Connor and Ogden (1995) for surface waves.
Figure 4 sketches what happens to a unit cube when a solid is subject to the simple
shear of amount K,
x̂1 = X1 + KX2 ,
x̂2 = X2 ,
x̂3 = X3 .
(4.4)
Here the principal axes are x3 = X3 and x1 , x2 which make an angle ψ with X1 and
with X2 , respectively. That angle, and the corresponding principal stretches are (e.g.
Chadwick (1976)),
p
ψ = (1/2) tan−1 (2/K), λ1,2 = 1 + K 2 /4 ± K/2, λ3 = 1.
(4.5)
These relations highlight a major difference between this homogeneous pre-strain and
the triaxial pre-stretch (3.1): here the orientation of the principal axes with respect to
the plane interface changes as the magnitude of the pre-strain changes.
^
^
x3
x3
^
^
x1
x1
K< 0
K >0
^
x2
^
x2
Figure 4. Finite simple shear of amount K of a block near the interface. The x̂1 direction
is the direction of shear ; the (x̂1 x̂2 ) plane is the plane of shear ; the (x̂1 x̂3 ) plane is the
glide plane.
4.1
Governing equations for simple shear pre-strain
The deformation (4.4) is an example of plane strain. As we focus on two-partial
incremental waves in this Section, we may take advantage of formulas established by
Merodio and Ogden (2002) in a similar context. The components of the deformation
gradient tensor F and of the left Cauchy-Green strain tensor B for 2D pre-strain and
2D incremental motions, in the (x̂1 , x̂2 ) coordinate system (aligned with the (X1 , X2 )
system), are
1 + K2 K
1 K
.
(4.6)
,
B=
F =
K
1
0 1
14
For plane strain, λ3 = 1, so that in incompressible solids λ2 = λ−1
1 by (2.1). It follows
c (I1 ) by the
by (2.4) that I1 = I2 . Accordingly, we define the single-variable function W
identity
c (I1 ) = W (I1 , I1 ).
W
(4.7)
Then the 2D version of the constitutive equation (2.5) is
c1 B,
σ = −p̂I + 2W
(4.8)
c1 δik Bjl + 4W
c11 Bij Blk ,
Â0jilk = 2W
(4.9)
where p̂ is the Lagrange multiplier due to the incompressibility constraint. Also, the
components of A0 are (Merodio and Ogden, 2002),
c1 = W
c 0 (I1 ), W
c11 = W
c 00 (I1 ). With the help of (4.6)2 , we find that in the (x̂1 , x̂2 )
where W
coordinate system, the non-zero components Â0jilk relevant to in-plane motion are
c11 K(1 + K 2 ),
Â01112 = Â01211 = 4W
c1 K + 4W
c11 K(1 + K 2 ),
Â01121 = Â02111 = 2W
c11 (1 + K 2 ),
Â01122 = Â02211 = 4W
c1 (1 + K 2 ) + 4W
c11 K 2 ,
Â01212 = 2W
c1 K + 4W
c11 K,
Â01222 = Â02212 = 2W
c1 (1 + K 2 ) + 4W
c11 (1 + K 2 )2 ,
Â01111 = 2W
c11 K 2 ,
Â01221 = Â02112 = 4W
c1 + 4W
c11 K 2 ,
Â02121 = 2W
c11 K,
Â02122 = Â02221 = 4W
c1 + 4W
c11 .
Â02222 = 2W
(4.10)
Now the remaining incremental governing equations (2.10)1,2,4 reduce to
0
−S11 + iS21
= −ρv 2 U1 ,
0
−S12 + iS22
= −ρv 2 U2 ,
iU1 + U20 = 0,
(4.11)
where by (2.7),
iS11 = i(Â01111 + p̂)U1 + iÂ01112 U2 + Â01121 U10 + Â01122 U20 − iP̂ ,
iS12 = iÂ01112 U1 + iÂ01212 U2 + (Â01221 + p̂)U10 + Â01222 U20 .
iS21 = iÂ01121 U1 + i(Â01221 + p̂)U2 + Â02121 U10 + Â02122 U20 ,
iS22 = iÂ01122 U1 + iÂ01222 U2 + Â02122 U10 + (Â02222 + p̂)U20 − iP̂ ,
(4.12)
where P̂ is the increment of p̂. The pre-stress necessary to maintain the solid in the
static state of large simple shear is (4.8). In particular, the σ̂22 component along the x̂2
axis is found using (4.6) and (4.7) as
leading to the connection
c1 ,
σ̂22 = −p̂ + 2W
Â01221 + p̂ = Â02121 − σ̂22 .
15
(4.13)
(4.14)
Note that Connor and Ogden (1995) and Fu (2005a) keep σ2 rather than σ̂22 as a measure
of the pre-stress. It is the component of the pre-stress along the principal axis x2 , whose
orientation changes with the pre-strain, in contrast to σ̂22 , the component of the Cauchy
pre-stress tensor along the unchanged normal to the interface. As pointed out by Hussain
and Ogden (2000), it is σ̂22 which is continuous across the interface of two bonded sheared
solids.
Following the same procedure as in Section 3.1, we find the Stroh formulation of the
governing equations in the expected form (2.10), where ξ = [U1 , U2 , S21 , S22 ]t and −N 1
and N 2 , are given by
1
σ̂22
0
K
1
−
,
c1 + 4W
c11 K 2
c1 + 4W
c11 K 2 , 2W
(4.15)
2W
0
0
1
0
respectively, and −N 3 by
4.2
c1 − σ̂22 )
2(4W
c1 − σ̂22 )K
−(4W
Resolution
c1 − σ̂22 )K
−(4W
2
.
σ̂22
c1 K 2 + 2σ̂22 −
2W
2
c
c
2W1 + 4W11 K
(4.16)
The propagation condition (2.16) reduces to a quartic in q,
c1 + 2W
c11 K 2 )q 4 + 4(W
c1 + 2W
c11 K 2 )Kq 3
2(W
c1 (2 + K 2 ) − 4W
c11 K 2 (2 − K 2 ) − ρv 2 ]q 2
+ [2(W
c1 − 2W
c11 K 2 )Kq
+ 4(W
c1 (1 + K 2 ) + 4W
c11 K 2 − ρv 2 = 0. (4.17)
+ 2W
Notice how σ̂22 , though present in N , does not appear explicitly in this equation.
Quite surprisingly, there are two instances – both pointed out by Connor and Ogden
(1995, 1996) – where we can solve this quartic exactly and simply. The first instance
occurs for incremental deformations, when v = 0; it is applicable whatever the strain
energy density W might be. The reason for the simplicity of this resolution is made
apparent by the change of unknown from q to qe = q − K/2. Then the quartic becomes
c1 + 2W
c11 K 2 )e
c1 (4 − K 2 ) − 2W
c11 K 2 (4 + K 2 ) − ρv 2 ]e
2(W
q 4 + [W
q2
c1 + 2W
c11 K 2 )(4 + K 2 ) − 2ρv 2 ](4 + K 2 )/8 = 0, (4.18)
+ (ρv 2 )K qe + [(W
which is clearly a biquadratic at v = 0. The consequence is that any incremental static
problem can be solved in its entirety for sheared solids because the roots q are accessible
explicitly. This was to be expected though, because of an important theorem by Fu and
Mielke (2002) which states that “the buckling condition for a pre-stressed elastic halfspace is independent of the orientation of the surface as long as the surface normal remains
16
in the (x1 , x2 ) plane”; the buckling condition is what corresponds to the marginally
stable static solution obtained at v = 0; it is also called the wrinkling condition, or the
bifurcation criterion, or any other denomination associated with the onset of instability
in the linearised (incremental) theory.
The second instance where the quartic is easy to solve is when the solid is a MooneyRivlin material, see (3.22); this case is treated in Section 4.3.
In general however, the quartic is difficult (but not impossible, see Fu (2005a,b),
Destrade and Fu (2006), and the concluding Section) to solve analytically and other
methods, such as those relying on the fundamental equations (2.25), are required. For
the time being, we complete the picture with formal calculations.
Assuming the roots of the quartic have been computed, and calling them q1 , q2 , q 1 ,
q 2 , where q1 and q2 both satisfy the decay condition (2.17), we find that the eigenvectors
ζ 1 and ζ 2 associated with q1 and q2 , respectively, are
1
2
a
a
ζ1 = 1 , ζ2 = 2 ,
(4.19)
b
b
where
aj = [qj2 , −qj ]t ,
t
bj = qj γ̂21 (qj2 + Kqj − 1) + σ̂22 , −(γ̂21 − σ̂22 )qj2 + ν̂12 qj + γ̂12 − ρv 2 .
(4.20)
Here,
c1 + 4W
c11 K 2 ,
γ̂21 = 2W
c1 (2 + K 2 ) − 2W
c11 K 2 (2 − K 2 ),
β̂12 = W
c1 (1 + K 2 ) + 4W
c11 K 2 ,
γ̂12 = 2W
c1 K − 4W
c11 K 3 .
ν̂12 = 2W
(4.21)
Then we compute the A = [a1 |a2 ] and B = [b1 |b2 ] matrices, and eventually the surface
impedance matrix M = −iBA−1 , as
γ̂21 (K + q1 + q2 )
γ̂21 (1 + q1 q2 ) + σ̂22
(4.22)
M = −i ρv 2 − γ̂12
.
q1 + q 2
− γ̂21 − σ̂22 (ρv 2 − γ̂12 )
− ν̂12
q1 q2
q1 q2
Note that this matrix is indeed Hermitian; this is easy to show by using the quartic (4.17)
and the definitions (4.21) to uncover the identities:
q1 + q2 + q 1 + q 2 = −2K,
q1 q2 q 1 q 2 = (γ̂12 − ρv 2 )/γ̂21 ,
q1 q2 (q 1 + q 2 ) + q 1 q 2 (q1 + q2 ) = −2ν̂12 /γ̂21 .
These identities allow us to rewrite the surface impedance matrix as
"
#
γ̂21 (q1 + q2 − q 1 − q 2 )/2
γ̂21 (1 + q1 q2 ) − σ̂22
M = −i
.
−γ̂21 [q1 q2 (q 1 + q 2 ) − q 1 q 2 (q1 + q2 )]/2
−γ̂21 (1 + q 1 q 2 ) + σ̂22
17
(4.23)
(4.24)
Fu (2005a) derived the same form of the surface impedance matrix for the more general
case of any plane strain (any pre-strain where λ3 = 1), which includes the present case
of finite shear.
Notice also how the pre-stress σ̂22 is going to be explicitly present in the secular
equation for Rayleigh surface waves (2.22), but absent from the secular equation for
Stoneley interface waves (2.23). These secular equations remain implicit as long as the
roots q1 and q2 are not known. In the general case where the quartic (4.17) is not solvable
in a simple manner, we seek an explicit secular equation using the fundamental equations
(2.25).
Rayleigh surface waves. With the explicit expressions (4.15) and (4.16) for the blocks
of the matrix N , we can compute N −1 and N 2 . The lower left corner of these gives
in turn K (1) = N 3 + ρv 2 I and K (−1) , K (2) , respectively. Then the equations (2.30)
written at n = 1, −1, 2, yield the linear homogeneous system
(1)
(1)
(1)
K11
K12
K22
U1 (0)U 1 (0)
0
(−1)
(−1)
(−1)
(4.25)
K11
K12
K22 U1 (0)U 2 (0) + U 1 (0)U2 (0) = 0 .
(2)
(2)
(2)
0
U2 (0)U 2 (0)
K
K
K
11
12
22
The vanishing of the determinant of the 3 × 3 matrix on the left hand side is the explicit
polynomial secular equation for surface waves in a sheared incompressible semi-infinite
solid. It is a polynomial of degree 4 in ρv 2 . It is too long to reproduce in general but
easy to obtain (and solve numerically) with a computer algebra system. Here we present
its expression in the case where σ̂22 = 0. Then, K (1) and K (2) are given by
#
#
"
"
c1 − ρv 2 )K
c1 (4 − K 2 ) − 2ρv 2
c1
c1
2(4W
2W
ρv 2 − 8W
4W
, (4.26)
,
c1 (4 − K 2 ) − 2ρv 2
c1 K
c1
c1
2W
−8W
4W
ρv 2 + 2W
respectively, the components of K (−1) are
(−1)
K11
(−1)
K12
(−1)
K22
c1 − ρv 2 ),
= 2(2W
c1 (1 + K 2 ) − ρv 2 ]K,
= −2[2W
2
2
2
2
c
c
c1 (2 + K 2 )2 + ρv 2 ρv − 2W1 (5 + K ) − 4W11 K (1 + K ) ,
= 4W
c1 + 2W
c11 K 2
W
(4.27)
(up to an inessential common factor), and the secular equation is the quartic
!
c
5 + K2
4
3
2 W11
x − 5x + 8
x2
+K
c
4 + K2
W1
!
c
c1 + 2W
c11 K 2
1 + K2
W
2 W11
+
K
= 0, (4.28)
−8 4
x
+
8
c1
c1 (4 + K 2 )
4 + K2
W
W
where x is the following non-dimensional measure of the squared wave speed, x =
c1 (4 + K 2 )].
ρv 2 /[W
18
As stated above, the secular equation is also a quartic in the squared wave speed
when σ̂22 6= 0. For a given material, a given pre-stress, and a given shear, its numerical
resolution may yield more than one positive real root. If such is the case, then for
each corresponding speed, compute the roots to the quartic (4.17) and discard those
(supersonic) speeds which do not give two complex conjugate pairs of roots. Finally, find
which of the remaining speeds (if there is more than one) satisfies the secular equation
written in optimal form (2.22).
Stoneley shear-twin interface waves. For Stoneley interface waves, the fundamental equations (2.25) are not practical to derive secular equations in general, and we must
resort to other methods, such as those developed by Destrade and Fu (2006). The exception is the special case when each half-space is filled with Mooney-Rivlin materials,
because then the roots to the quartic (4.17) can be found easily, see Section 4.3.
We now focus on the possibility of propagating incremental waves at a shear-twin
interface. In this configuration, two solids, made of the same incompressible material, are
subject to equal and opposite shears, see Figure 5. The study of wave propagation at this
type of interface can have important repercussions in the non-destructive evaluation of a
twinned interface because this bimaterial can “simulate the finite (plastic) deformation
associated with a crystal twin” (Hussain and Ogden, 2000).
-K
x^3
x^1
K
x^2
Figure 5. Stoneley wave propagation at a shear-twin interface. Both half-spaces are
occupied with the same incompressible solid, one subject to a simple shear of amount
K, the other subject to a simple shear of amount −K. The analysis shows that Stoneley
waves cannot actually travel in the direction of shear.
For one half-space, the propagation condition is the quartic (4.17). For the other halfc1 and W
c11 do not change
space, the amount of shear is changed to its opposite, but W
2
signs because they are functions of K . It follows that if q1 and q2 are qualifying roots
of the propagation condition in one half-space, then −q1 and −q2 are qualifying roots in
the other half-space. Also, we see from (4.21) that γ̂21 , γ̂12 , and β̂21 in one half-space
19
are equal to their counterparts in the other half-space, but that the counterpart to ν̂ 12 is
−ν̂12 . Finally, σ̂22 is continuous across the interface. Then we conclude from (4.20) that
the counterparts to the matrices
B11 B12
A11 A12
,
(4.29)
,
B=
A=
B21 B22
A21 A22
in one half-space are the matrices (Destrade, 2003; Ting, 2005)
A11
A12
−B11 −B12
∗
∗
A =
,
B =
,
−A21 −A22
B21
B22
(4.30)
in the other half-space. In turn, this conclusion leads to the following sum of surface
impedance matrices
#
"
2M11
0
∗
,
(4.31)
M +M =
0
2M22
where M ∗ = iB ∗ (A∗ )−1 . Consequently, the exact secular equation for Stoneley sheartwin interface waves is, according to (2.23), (4.24), and (4.31) that either
=(q1 + q2 ) = 0,
or =[q1 q2 (q 1 + q 2 )] = 0.
(4.32)
However, as is easily proved, neither of these quantities can be zero when both q 1 and
q2 have positive imaginary parts. It follows that Stoneley waves cannot propagate in
the direction of sheared at a shear-twin interface, whatever the strain energy function is,
and whatever what the shear is. Note however that Hussain and Ogden (2000) show, in
their study of reflection and transmission of plane waves at shear-twin interface, that an
incident harmonic plane wave can give rise to an interfacial wave (but of a different type
than the Stoneley type).
4.3
Examples
Sheared Mooney-Rivlin solids. For the Mooney-Rivlin strain energy density (3.22)
we have
c1 = D/2,
c11 = 0,
W
W
(4.33)
where D = D1 + D2 is the shear modulus. The quartic propagation condition (4.17) then
factorizes to
p
(q 2 + 1)(q 2 + 2Kq + K 2 + η 2 ) = 0, where η = 1 − ρv 2 /D,
(4.34)
and its roots with positive imaginary parts are
q1 = i,
q2 = −K + iη.
(4.35)
We also have
γ̂21 = D,
γ̂12 = D(1 + K 2 ),
β̂12 = D(2 + K 2 ),
20
ν̂12 = DK,
(4.36)
and the surface impedance matrix of (4.24) reduces to
"
#
η+1
−K + i(η − 1) + iσ̂22 /D
M =D
.
−K − i(η − 1) − iσ̂22 /D
η2 + η + K 2
(4.37)
For Rayleigh surface waves in a sheared Mooney-Rivlin material, the secular equation
(2.22) is a cubic in η,
η 3 + η 2 + (3 + K 2 − 2σ̂22 /D)η − (1 − σ̂22 /D)2 = 0.
(4.38)
See Connor and Ogden (1995) for this equation with σ2 instead of σ̂22 , and Destrade and
Ogden (2005) for a generalization of this equation to a triaxial stretch followed by a shear.
See also Figure 6 for the variations of the squared scaled wave speed ρv 2 /D = 1 − η 2 with
the amount of shear K, for several values of the pre-stress σ̂22 . In particular, the plots
at σ̂22 = ±2D show that the Mooney-Rivlin material is unstable when subject only to a
hydrostatic pressure of that amount (that is, ρv 2 = 0 at K = 0) but regains stability as
soon as it is sheared (ρv 2 > 0 at K 6= 0) .
rv2/
0.5
1
1.5
0.8
2.0
0.6
-2.0
0.4
-1.5
0.2
K
– 3
– 2
– 1
0
1
2
3
Figure 6. Squared scaled surface wave speed ρv 2 /D as a function of the amount of shear
K for a Mooney-Rivlin material. For the solid plot, there is no pre-stress normal to the
boundary (σ̂22 = 0); for the dashed plots, the displayed number indicates the value of
the non-dimensional pre-stress σ̂22 /D. More plots are found in the article by Connor and
Ogden (1995).
As is clear from the roots (4.35), from the secular equation (4.38), and from Figure 6,
both the penetration depth and the speed are even functions of K. Also, once η is found
by solving the secular equation, we have ∂η/∂(K 2 ) = −η 2 /[2η 3 + η 2 + (1 − σ̂22 /D)2 ] < 0
and ∂=(q2 )/∂(K 2 ) = ∂η/∂(K 2 ) < 0, indicating that both the surface wave speed and
the penetration depth increase with the magnitude of the shear.
21
We note that having found an explicit expression (4.37) for the surface impedance
matrix allows us to solve the problem of Stoneley interface waves in its entirety, for halfspaces made of different Mooney-Rivlin solids, subject to different amounts of shear; see
Chadwick and Jarvis (1979b) for more discussion on this point.
Sheared Gent solids. As an example of a strain energy density for which the quartic
(4.17) is not easily solved, we now work with the following function,
I1 − 3
W = − 21 CJm ln 1 −
,
Jm
I1 < 3 + J m ,
(4.39)
proposed initially by Gent (1996) to describe strain-stiffening elastomers and since then
extensively used by Horgan and Saccomandi (2003, and references therein) to model
strain-stiffening soft biological tissues such as arteries. Here C > 0 is the infinitesimal
shear modulus and Jm > 0 is a constant, accounting for the limiting chain extensibility
of the solid. Hence the condition (4.39)2 limits the amount of shear by imposing
p
p
− Jm < K < Jm .
(4.40)
Turning our attention to the propagation of surface waves, we compute the following
quantities,
CJm
CJm
c1 =
c11 =
W
,
W
.
(4.41)
2
2(Jm − K )
2(Jm − K 2 )2
To fix the ideas, we take two Gent materials, one with Jm = 9.0, the other (stiffer) with
Jm = 1.0. Figure
7 shows the variations of the non-dimensional measure of the surface
p
wave speed ρv 2 /C with the amount of shear K. The bounds due to the inequalities
(4.40) are clearly visible, indicating that the solids become more and more rigid as their
limit of chain extensibility is approached. The speed is an even function of K and only
the K > 0 needs be displayed. For a boundary free of pre-stress (σ̂22 = 0), we solve
numerically the quartic secular equation (4.28) and find in general that there is two real
positive roots for ρv 2 ; one gives a supersonic speed and is discarded; the other gives a
speed for which the exact secular equation (2.22) is satisfied and it is kept. The Figure
also shows (dotted plots) the effect of the compressive pre-stress σ̂22 = −1.5C, which is to
slow the wave down in the small-to-moderate shear region; here again a quartic secular
equation is solved numerically and only one speed is kept.
5 Propagation in a principal plane
In this Section we consider an interface wave propagating in a principal plane, x2 = 0 say,
but not in a principal direction. We call θ the angle between the propagation direction
x̂1 and the principal axis x1 , see Figure 8. Although some results exist for non-principal
Stoneley waves (Destrade, 2005), the focus of this Section is on Rayleigh surface waves.
22
[ρv 2/C ]1/2
3
2.5
2
1.5
1
0.5
K
0
0.5
1
1.5
2
2.5
3
p
Figure 7. Scaled surface wave speed ρv 2 /C as a function of the amount of shear K for
two Gent materials. For the solid plots, there is no pre-stress normal to the boundary
(σ̂22 = 0); for the dashed plots, the pre-stress is σ̂22 = −1.5C. The vertical asymptotes
correspond to the maximal amount of shear, that is Kmax = ±3.0 for the Gent solid with
Jm = 9.0, and Kmax = ±1.0 for the other, stiffer, Gent solid with Jm = 1.0.
5.1
Governing equations
Hence we model this motion as
{u, ṗ, s} = {U (kx2 ), ikP (kx2 ), ikS(kx2 )}eik(cθ x1 +sθ x3 −vt) ,
(5.1)
where cθ = cos θ, sθ = sin θ. By projecting the governing equations in the coordinate axis
of the principal axes (where A0 has 15 independent non-zero components, see Section
2.2), we can write them in the Stroh form (2.13) as a homogeneous linear system of six
first-order differential equations, where ξ = [U1 , U2 , U3 , S21 , S22 , S23 ]t , see Destrade et al.
(2005) for details. Here the matrices −N 1 , N 2 , −N 3 are given by
0
cθ
0
cθ (γ21 − σ2 )/γ21
0
sθ (γ23 − σ2 )/γ23
0
sθ ,
0
1/γ21
0
0
23
0
0
0
0 ,
0 1/γ23
χ 0
0 ν
−κ 0
−κ
0 ,
µ
(5.2)
x3
^
x1
q
x1
l2
1
l3
1
1
x2
l1
Figure 8. Interface wave propagation when the boundary x2 = 0 is a principal plane of
pre-strain. The x̂1 direction is the direction of propagation, making an angle θ with the
principal direction of pre-strain x1 .
respectively, where
χ = 2c2θ (β12 + γ21 − σ2 ) + s2θ γ31 ,
ν = c2θ [γ12 − (γ21 − σ2 )2 /γ21 ] + s2θ [γ32 − (γ23 − σ2 )2 /γ23 ],
µ = c2θ γ13 + 2s2θ (β23 + γ23 − σ2 ),
κ = cθ sθ (β13 − β12 − β23 − γ21 − γ23 + 2σ2 ),
(5.3)
and the γij , βij are defined in (3.12).
Rogerson and Sandiford (1999) show that the propagation condition (2.16) is a cubic
in q 2 ,
γ21 γ23 q 6 − [(γ21 + γ23 )X − c1 ] q 4 + (X 2 − c2 X + c3 ) q 2 + (X − c4 )(X − c5 ) = 0, (5.4)
with X = ρv 2 and
c1 = (γ21 γ13 + 2β12 γ23 )c2θ + (γ23 γ31 + 2β23 γ21 )s2θ ,
c2 = (γ23 + γ13 + 2β12 )c2θ + (γ21 + γ31 + 2β23 )s2θ ,
c3 = (γ12 γ23 + 2β12 γ13 )c4θ + (γ21 γ32 + 2β23 γ31 )s4θ
+ [γ12 γ21 + γ13 γ31 + γ23 γ32 − (β13 − β12 − β23 )2 + 4β12 β23 ]c2θ s2θ ,
c4 = γ12 c2θ + γ32 s2θ ,
c5 = γ13 c4θ + 2β13 c2θ s2θ + γ31 s4θ .
(5.5)
24
5.2
Resolution for Rayleigh surface waves
Finding the eigenvectors ζ 1 , ζ 2 , ζ 3 of N corresponding to the qualifying roots q1 , q2 ,
q3 (say) is a lengthy task by hand, best left to a computer algebra system. In the end,
we find that the ζ i are of the forms:
1
2
3
a
a
a
1
2
3
ζ = 1 ,
ζ = 2 ,
ζ = 3 ,
(5.6)
b
b
b
where
a4 qi4 + a2 qi2 + a0
aj = −qj5 + b3 qj3 + b1 qj ,
d4 qj4 + d2 qj2 + d0
Here the quantities m and n are given by
m=
n=
h3 qj3 + h1 qj
bj = (ν − X)(qj4 + mqj2 + n) .
g3 qj3 + g1 qj
(5.7)
1
(γ21 − σ2 )2 2
(γ23 − σ2 )2 2
1
+ 2
cθ [η − X] +
+ 2
sθ [µ − X]
γ21
γ21 (ν − X)
γ23
γ23 (ν − X)
(γ21 − σ2 )(γ23 − σ2 )
cθ sθ , (5.8)
− 2κ
γ21 γ23 (ν − X)
(γ21 − σ2 )2 2 (γ23 − σ2 )2 2
−1
1+
cθ +
sθ (ν − X)
γ21
γ23
× [(µ − X)(η − X) − κ2 ]/(γ21 γ23 ). (5.9)
The expressions for the quantities ai , bi , di , hi , gi are too lengthy to reproduce but they
are easily obtained in a formal manner using a computer algebra system; as it turns out,
these constants are not needed in the secular equation for Rayleigh surface waves. Indeed,
the expression for the surface impedance matrix M is very lengthy, but its determinant
factorizes greatly. Hence we find that the exact secular equation for Rayleigh surface
waves (2.22) reduces to (Taziev, 1987; Destrade et al., 2005)
nωI − ωIII (m − ωII ) = 0,
(5.10)
where
ωI = −(q1 + q2 + q3 ),
ωII = q1 q2 + q2 q3 + q3 q1 ,
ωIII = −q1 q2 q3 .
(5.11)
This secular equation remains implicit as long as the roots q1 , q2 , q3 satisfying the decay
condition are not known. To compute them, we must find the wave speed, by using the
fundamental equations (2.30).
First, using the explicit expression (5.2) of the Stroh matrix N , we compute N −1
and N 3 , and in particular we find explicit expressions for the lower left blocks K (1) =
N 3 + ρv 2 I, K (−1) , K (3) . Next, we use the result (Destrade, 2005) that U (0) is in the
form
U (0) = U1 (0)[1, iα, β]T ,
(5.12)
25
where α, β are real numbers, to write the fundamental equations (2.30) at n = −1, 1, 3
as the following system of three equations,
(−1)
(−1)
(−1)
(−1)
−K11
K13
K33
K22
2β
(1)
(1)
(1)
(1)
(5.13)
K33
K22 β 2 = −K11 .
K13
(3)
(3)
(3)
(3)
α2
−K11
K13
K33
K22
Then, we solve this non-homogeneous system by Cramer’s rule to find
β 2 = ∆2 /∆,
2β = ∆1 /∆,
(5.14)
where ∆, ∆1 , and ∆2 are the following determinants,
(−1)
K13
(1)
∆ = K13
(3)
K13
(−1)
K33
(1)
K33
(3)
K33
(−1)
K22
(1)
K22 ,
(3)
K22
(−1)
−K11
(1)
∆1 = −K11
(3)
−K11
(−1)
K33
(1)
K33
(3)
K33
(−1)
K22
(1)
K22 ,
(3)
K22
(−1)
K13
(1)
∆2 = K13
(3)
K13
(−1)
−K11
(1)
−K11
(3)
−K11
(−1)
K22
(1)
K22 . (5.15)
(3)
K22
Finally we write down the compatibility of equations (5.14) as
∆21 − 4∆∆2 = 0,
(5.16)
which is the explicit secular equation for non-principal surface waves in deformed incompressible materials.
In general this equation is a polynomial of degree 12 in ρv 2 (Taziev, 1989; Destrade
et al., 2005), easy to solve numerically. Of the 12 possible roots, we keep those which are
real, positive, and give a subsonic speed (that is a speed for which the bicubic (5.4) has
three pairs of complex conjugate roots). We then test the remaining speeds against the
exact secular equation (5.10).
5.3
Examples
neo-Hookean solid. The neo-Hookean form of the strain energy density for an incompressible isotropic solid is a sub-case of the Mooney-Rivlin form, namely D2 = 0 in
(3.22) so that
W = D1 (λ21 + λ22 + λ23 )/2.
(5.17)
It leads to a stress-strain relationship (2.5) which is “linear” with respect to the left
Cauchy-green strain tensor,
σ = −pI + D1 B.
(5.18)
26
The neo-Hookean model is unable to capture neither qualitatively nor quantitatively
experimental data, over any range of deformations (Saccomandi, 2004). It is nonetheless a
highly popular model in the literature because it forms the basis of a statistical treatment
for the molecular description of rubber elasticity (Treloar, 1949).
With respect to wave propagation, it has more peculiar properties than the MooneyRivlin material, because here,
γij = D1 λ2i ,
2βij = D1 (λ2i + λ2j ),
(5.19)
which leads to great simplifications in the matrices N 1 , N 2 , N 3 of (5.2).
Now placing ourselves in the coordinates system (x̂1 , x2 , x̂3 ) attached to the wave
propagation, see Figure 8, we introduce the functions Ûi , Ŝ2i (i = 1, 2, 3), defined by
Ûi = Ωij Uj ,
Ŝ2i = Ωij S2j ,
cθ
where Ωij = 0
sθ
0 −sθ
1
0 .
0 cθ
(5.20)
With these functions, the governing equations (2.13) decouple the anti-saggital motion
[Û3 , Ŝ23 ] from its saggital counterpart (recall that the direction of propagation and the
normal to the interface define what is called the saggital plane.) For the latter motion
we find
0
0 t
[Û10 , Û20 , Ŝ21
, Ŝ22
] = iN [Û1 , Û2 , Ŝ21 , Ŝ22 ]t ,
(5.21)
with
where
0
−1
N =
ρv 2 − χ̂
0
−1 + σ 2
0
0
ρv 2 − ν̂
χ̂ = D1 (c2θ λ21 + s2θ λ23 + 3λ22 − 2σ 2 ),
1/(D1 λ22 )
0
0
−1 + σ 2
0
0
,
−1
0
ν̂ = D1 [c2θ λ21 + s2θ λ23 + λ22 (1 − σ 2 )2 ],
(5.22)
(5.23)
and σ 2 = σ2 /(D1 λ22 ) is a non-dimensional measure of the pre-stress. The associated
propagation condition is
(q 2 + 1)[λ22 q 2 + (c2θ λ21 + s2θ λ23 ) − ρv 2 /D1 ] = 0.
(5.24)
In other words, the situation is formally the same as that for principal waves in
Mooney-Rivlin material, see Section 3.3. The conclusion is that the speed is given by
ρv 2 = D1 (c2θ λ21 + s2θ λ23 − λ22 η 2 ),
(5.25)
where η is the real root of (3.26). Flavin (1963) established this result at σ 2 = 0. As he
also showed, the situation gets more complicated for Mooney-Rivlin solids, because the
wave is no longer plane polarized for a triaxial pre-stretch.
27
Mooney-Rivlin solid. The Mooney-Rivlin strain energy density (3.22) at D2 6= 0 does
not lead to a decoupling of the saggital motion from the anti-saggital motion. However,
as in Sections 3 and 4, we find that the propagation condition factorizes, here as
(q 2 + 1)(q 4 − Sq 2 + P ) = 0,
(5.26)
where
S=
1
1
+
γ21
γ23
X−
γ12
γ13
+
γ21
γ23
c2θ −
γ31
γ32
+
γ21
γ23
s2θ ,
P = (X − γ12 c2θ − γ32 s2θ )(X − γ13 c2θ − γ31 s2θ )/(γ21 γ23 ).
(5.27)
Flavin (1963) was the first to notice that q1 = i is an attenuation factor for nonprincipal waves in Mooney-Rivlin solids. Pichugin (2001) showed that a necessary and
sufficient condition for the factorization (5.26) to occur is that the relations (3.24) hold;
he also showed, completing earlier work by Willson (1973a), that another factorization of
the general bicubic (5.4) also occurs, this time for any strain energy function, when two of
the principal stretches of pre-strain are equal (equi-biaxial pre-strain). Finally note that
(q 2 + 1) always comes out as a factor in the propagation condition for inhomogeneous
waves in Mooney-Rivlin solids, whatever the direction of propagation is (not necessarily
in a principal plane as here), see Destrade (2002) for details.
Thanks to the factorization (5.26), we can actually compute explicitly the quantities
ωI , ωII , ωIII of (5.11). Indeed, we now have
q
√
√
q1 = i,
q 2 q3 = − P ,
q2 + q3 = i 2 P − S,
(5.28)
so that
ωI = −i(1 +
q
√
2 P − S),
−ωII =
√
P+
q
√
2 P − S,
√
ωIII = i P ,
leading to an explicit and exact form of the secular equation (5.10),
q
q
√
√
√
√
n 1 + 2 P − S + P m + P + 2 P − S = 0.
(5.29)
(5.30)
This result allows us to investigate the influence of pre-stress on surface wave propagation (see Destrade et al. (2005) for an example) and to address an important question in
the study of surface stability: how much can a Mooney-Rivlin half-space be compressed
before it buckles? In Section 3 we found the critical stretch in a principal direction (x1 ),
indicating the appearance of wrinkles parallel to x3 , see Figure 2 for a visualization.
However, could it be that wrinkles appeared earlier in the compression, in another direction? To answer this we take X(= ρv 2 ) = 0 (onset of instability) and solve (5.30) for
λ, for each value of θ.
For instance, take the case of the following plane pre-strain,
λ1 = λ,
λ2 = λ−1 ,
28
λ3 = 1,
(5.31)
imposed on a half-space made of the Mooney-Rivlin solid with material parameters
D1 = 2.0µ,
D2 = 0.8µ,
(5.32)
where µ has the dimension of a stiffness (the shear modulus of this Mooney-Rivlin solid is
(D1 + D2 )/2 = 1.4µ). Figure 9 displays the values of the critical stretch ratio, measured
in the principal direction x1 , for each angle θ and for several values of the pre-stress σ2 .
The solid line corresponds to σ2 = 0. At θ = 0 we find the critical compression ratio
λ0c (say) for wrinkles aligned along x3 ; hence the solid curve starts at λ0c = 0.544 as
expected from solving (3.27) with v = 0 and λ1 = λ0c = λ−1
2 . At θ 6= 0 we find that the
critical compression stretch is below λ0c ; it follows that λ0c is the absolute critical stretch
of compression for our example (5.31), (5.32). Of course, for Mooney-Rivlin solids other
than (5.32), or for pre-strains other than (5.31), or for solids other than Mooney-Rivlin
solids, we might end up with a different behaviour in compression. However the same
analysis can be brought in each case to its conclusion, with no additional difficulty.
lc
0.7
0.6
0.5
0.4
3.0
1.0
0.0
-1.0
-3.0
0.3
0.2
0
20
4 0
6 0
8 0
q
Figure 9. Critical stretch ratio of compression for a Mooney-Rivlin material as a function of the angle between the normal to the wrinkles and the principal axis of greatest
compression. The Mooney-Rivlin solid is subject to a finite plane strain compression
(λ3 = 1); its material parameters are D1 = 2.0µ, D2 = 0.8µ; the pre-stress normal to the
boundary is given by σ2 /µ = -3.0, -1.0, 0.0, 1.0, 3.0 (µ has the dimension of a stiffness).
Gent solid. We find that the elastic moduli (3.12) for Gent materials (4.39) are given
in general by
Jm
λ2 ,
Jm + 3 − λ21 − λ22 − λ23 i
"
#
2(λ2i − λ2j )2
Jm
2
2
2βij = C
λi + λ j +
.
Jm + 3 − λ21 − λ22 − λ23
Jm + 3 − λ21 − λ22 − λ23
γij = C
29
(5.33)
Now consider that a half-space made of a Gent material with Jm = 9.0 is subject to
a large shear such that the boundary is the plane of shear (in Section 4, the boundary
was the glide plane) and the boundary x2 = 0 is free of tractions. Then the principal
stretches are
p
p
λ1 = 1 + K 2 /4 + K/2,
λ2 = 1,
λ3 = 1 + K 2 /4 − K/2,
(5.34)
and of course, λ21 + λ22 + λ23 − 3 = K 2 .
Figure 10 shows the variations of the surface wave speed in the plane of shear, with
respect to the angle between the direction of propagation (x̂1 ) and the direction of shear
(X1 ). For a shear of amount K = 0.5, the principal axis of greatest stretch x1 makes
an angle 37.9◦ with the direction of shear, and the principal axis of smallest stretch x3
makes an angle 127.9◦ with the direction of shear. For K = 1, those two angles are 31.7◦
and 121.7◦ , respectively; for K = 2, those two angles are 22.5◦ and 112.5◦ , respectively.
Clearly, the surface wave reaches extremal values along those directions, indicating that
the principal directions can be determined acoustically.
[rv 2/C] 1/2
3
K = 2.0
2.5
2
K = 1.0
1.5
K = 0.5
1
0.5
q
0
20
40
60
80
100
120
140
160
180
Figure 10. Wave propagating
p in the plane of shear of a semi-infinite sheared Gent solid
with Jm = 9.0: scaled speed ρv 2 /C as a function of the angle between the propagation
direction and the direction of shear, for several values of the amount of shear.
6 Concluding remarks
This Chapter has presented several situations where it is possible to derive explicit secular
equations in exact form and in polynomial form, for incremental waves propagating
along the plane interface of one (or two) deformed semi-infinite hyperelastic solid (rigidly
30
bonded solids). The subject of interface waves in general is broad and includes many
other situations and geometries such as, to name but a few: the addition of a layer of finite
thickness (Love waves, Lamb waves, etc.), or of a fluid (Scholte waves for inviscid fluids,
Stoneley waves for viscous fluids, etc.), the consideration of an anisotropy due to families
of extensible fibres, of cylindrical or spherical coordinates, or of curved boundaries, and
so on. Reading the works by Ogden (2003, 2004) and by Guz (2002) and the references
therein gives a good overview of the vast panorama spanned by these other types of
interface wave problems.
This concluding Section expands on the reasons put forth to explain why it can be
advantageous at times to seek explicit secular equations rather than to turn to numerics
outright. Among other things, explicit secular equations in polynomial form
(i) are easy to solve numerically, with the greatest precision required;
(ii) are sometimes surprisingly simple and short, in spite of a complicated or even
unsolvable propagation condition;
(iii) lend themselves to simple asymptotic treatments, leading to approximate analytical
expressions for the wave speed;
(iv) account for all the solutions satisfying the propagation condition and the boundary
condition at the interface.
6.1
Numerics
Point (i) goes without saying because the numerical methods used to determine the
roots of a polynomial are safe and robust. Of course, at most only one root corresponds to
the actual solution and all the other roots must be discarded as being spurious. So, once
a likely root is found numerically (that is a root ρv 2 which is real and positive), it must
be checked that at that speed the exact secular equation is satisfied. This routine check
is simple enough to perform (solve the propagation condition, find the corresponding
impedance matrix, check that (2.22) (for Rayleigh waves) or that (2.23) (for Stoneley
waves) is satisfied).
Other numerical techniques to find the interface wave speed invoke the Barnett-LotheStroh integral formalism (see Ting (1996) for an account), computational algorithms for
eigenvalue problems Taylor (1981), or an algebraic Ricatti equation for the impedance
matrix (see the Chapter by Fu in this book and the references therein) to determine
M (v) for any numerical value of v; then v is varied in order to satisfy the exact secular
equation up to the desired precision. This however is not always an easy task, as seen in
the following example.
Consider the Mooney-Rivlin solid with material parameters (5.32), maintained in a
state of static pre-strain, with (Rogerson and Sandiford, 1999),
√
√
λ2 = 0.7,
λ3 = (λ1 λ2 )−1 ,
λ1 = 3.695,
σ2 = 0.8µ.
(6.1)
Take the interface between the solid and the vacuum to be the principal plane x2 = 0
and study the propagation of a Rayleigh surface wave in a direction close to x 3 (θ is close
to 90◦ in Figure 8.)
Along x3 we have a principal wave, travelling with a speed v(90)
p say, found from
Section 3.2. Here (Destrade et al., 2005) we find that v(90) ' 1.327 µ/ρ. On the other
31
hand, a shear (homogeneous) bulk wave linearly polarized along x1 travels with speed v1
say, given by ρv12 = γ31 , and a shear (homogeneous) bulk wave linearly polarized p
along
2
x2 travels with speed
v
say,
given
by
ρv
=
γ
.
Here
we
find
that
v
'
1.384
µ/ρ
2
32
1
2
p
and that v2 ' 0.995 µ/ρ, showing that the surface wave travels with a speed which is
intermediate between those of the shear bulk waves. That principal wave is two-partial
and polarized in the (x2 x3 ) plane, which is why it can afford to be faster than the shear
bulk wave polarized in the x1 direction. Also, it is isolated, in the sense that the transition
toward a surface wave propagating in a direction θ 6= 90◦ is abrupt, because this latter
wave is tri-partial and must therefore travel with a speed v(θ) say, which
p is less than the
speed
of
any
homogeneous
bulk
shear
wave,
in
particular
less
than
c4 /ρ and less than
p
c5 /ρ, where c4 and c5 are defined in (5.5). Figure 11 shows the variations of these 3
speeds in the neighbourhood of the x3 axis and makes those features clear.
1.5
_
r
._ v
m
1.4
1.3
principal surface wave
bulk shear waves
1.2
1.1
non-principal surface wave
1
78
80
82
84
86
88
90
q
Figure 11. Wave propagation in the x2 principal plane ofpsemi-infinite Mooney-Rivlin
solid (5.32) deformed by (6.1): variations of scaled speed ρv 2 /µ as a function of the
angle, near the x3 direction.
Now if our numerical method for the surface wave speed relies on using the known
principal wave speed as an “initial guess”, then it might run into difficulties here because
of the isolation of the principal wave. This special feature is characteristic of strongly
anisotropic crystals in linear elasticity. In incremental dynamics, it can occur at will,
simply by deforming the solid sufficiently to create a strong strain-induced anisotropy.
Also, as is clear from the Figure, the surface
p wave travels with a speed which is
extremely close to that of the bulk shear wave c4 /ρ. Hence at θ = 89.99◦, the speed
of the former is given by
p
ρv 2 /µ = 0.9948641879596778 . . .,
(6.2)
32
whilst the speed of the latter is given by
p
c4 /µ = 0.9948641879596860 . . ..
(6.3)
Consequently, if a numerical scheme for the surface wave speed relies on increasing v in
small steps until det M (v) = 0 is satisfied to any desired precision, then it will have to
be pushed to the 14th significant digit to insure that the speed does not correspond to a
bulk wave! Here the explicit polynomial secular gives only two positive roots, and it is a
routine check to verify that only (6.2) satisfies the exact secular equation.
6.2
Polynomials
For an interface wave propagating in any direction in a deformed solid, the propagation condition is a sextic in general. Such is for instance the case for a wave propagating
in any direction in the glide plane of the sheared block in Figure 4. Although a sextic is
unsolvable analytically, it is nonetheless possible to find a polynomial secular equation
for a Rayleigh surface wave, see Taziev (1989) or Ting (2004) for details. The polynomial turns out to be of degree 27 in ρv 2 . For Stoneley interface waves it does not seem
practical to look for a polynomial secular equation.
For an interface wave propagating in any direction in a principal plane of pre-deformation, the propagation condition is the bicubic (5.4). Although formulas exist for the
roots of a cubic, the analytical resolution of (5.4) proves to be complicated because it
is not known a priori whether the roots are purely imaginary or have a non-zero real
part and accordingly, which formula should be used. However, the polynomial explicit
secular equation for a Rayleigh surface wave is of degree 12 in ρv 2 as seen earlier. This
is also the case in the symmetry plane of an orthorhombic or monoclinic crystal. In
the symmetry plane of a cubic crystal, the degree of the polynomial comes down to 10,
as noted by Taylor (1981). For a two-partial surface wave, coupled to electrical fields
through piezoelectricity in 2mm crystals, the polynomial is also of degree 10 in ρv 2 for
metallized boundary conditions (there the propagation condition is also a bicubic, see
Collet and Destrade (2005)).
For a two-partial (non-principal) interface wave polarized in a principal plane of predeformation, the propagation condition is a quartic,
q 4 + d3 q 3 + d2 q 2 + d1 q + d0 = 0,
(6.4)
say (see (4.17) for the case of simple shear). In contrast to the case of the bicubic above,
we know here what the form of the roots should be, because the roots of a quartic are
either two pairs of complex conjugate numbers, or one double real root with one pair of
complex conjugate numbers, or four real roots. Here only the first scenario is acceptable
in order to satisfy the decay condition (2.17) and thus a single formula, always valid, is
required; it can be found in the textbooks. First introduce in turn the quantities r, s,
and h, defined by
r = d2 − 38 d23 ,
s = d1 − 21 d2 d3 + 81 d23 ,
h = d0 − 41 d1 d3 +
the quantities λ and φ ∈]0, π/3], defined by
3/2
1
1
12h + r2
,
φ = 31 arccos 2λ
λ = 27
33
2 3
27 r
1
2
16 d2 d3
−
+ s2 − 38 rh
3 4
256 d3 ,
,
(6.5)
(6.6)
and the quantities z1 , z2 , and z3 , defined by
z1 = 2λ1/3 cos(φ) − 32 r,
z2 = 2λ1/3 cos(φ + 2π/3) − 23 r,
z3 = 2λ1/3 cos(φ + 4π/3) − 32 r.
(6.7)
Then the qualifying roots are
√
√
√
p1 = 21 sign(s) z1 − 14 d3 + 21 i( −z2 + −z3 ),
√
√
√
p2 = − 21 sign(s) z1 − 41 d3 + 21 i( −z2 − −z3 ),
(6.8)
where sign(s) equals 1 if s is non-negative and −1 otherwise. These formulas are perfectly well handled by a computer algebra system, so that M , and ultimately the exact
secular equation, can be found. It then means that any interface wave problem can be
solved exactly. Nonetheless it might still be rewarding to look for the polynomial secular
equation, to check whether they turn out to be simple. For instance, the polynomial
secular equation for Rayleigh waves is a quartic in the squared speed when the solid is
sheared (Section 4.2), or tri-axially stretched and then sheared (Destrade and Ogden,
2005), or when the wave is polarized in the symmetry plane of a crystal (Currie, 1979).
Now consider a flexural wave travelling along the free edge of a thin (Love-Kirchhoff)
orthotropic plate where the edge makes an arbitrary angle with a principal axis of symmetry. Thompson et al. (2002) show that the corresponding propagation condition is
a quartic, but turn to a numerical resolution. However Fu (2003) shows that the polynomial secular equation is just a cubic in the squared speed. Finally consider the case
of a piezoacoustic Bleustein-Gulyaev surface wave in a rotated Y -cut about the Z axis
for 4 crystals. There also the propagation condition is a quartic, but the fundamental
equations reveal that the polynomial secular equation for metallized boundary conditions
is simply a quadratic in the squared speed (Collet and Destrade, 2004). Clearly it is a
worthy enterprise to unearth these polynomials rather than use numerics or the formulas
(6.5)-(6.8).
6.3
Approximate expressions
Once the polynomial secular equation is established, in the form P(v) = 0, say, it is a
straightforward matter to derive an analytical approximation for v. Calling v 0 an initial
approximation, we find in the first order that
v ' v0 − P(v0 )/P 0 (v0 ).
(6.9)
Of course we must choose v0 judiciously for an optimal expression.
Here we work out a simple example. We take a solid half-space subject to a hydrostatic
pressure only, so that λ1 = λ2 = λ3 ≡ λ say, and σ1 = σ2 = σ3 ≡ σ say. From
the incompressibility constraint (2.1), λ = 1 follows and we have a pre-stressed, but
unstrained, solid. It is thus isotropic and a surface wave propagates with the same speed
in every direction. To derive the secular equation, we specialize for instance (3.20) to the
34
isotropic case. We find that γ12 = γ21 = β12 = µ0 , the infinitesimal shear modulus, and
the exact secular equation is (Dowaikh and Ogden, 1990)
(6.10)
f (η) = η 3 + η 2 + (3 − 2σ)η − (1 − σ)2 = 0,
p
p
where σ ≡ σ/µ0 and η = 1 − ρv 2 /µ0 . Calling c ≡ ρv 2 /µ0 a dimensionless measure of
the wave speed, and multiplying (6.10) by f (−η), gives the polynomial secular equation
(Dowaikh and Ogden, 1990),
P(v) ≡ c6 − 4(2 − σ)c4 + 6(σ − 2)2 c2 + (σ 2 − 4)(σ − 2)2 = 0.
(6.11)
In the region where c is close to 1, the first order approximation (6.9) gives c ' 1 −
P(1)/P 0 (1), that is
21 − 28σ + 6σ 2 + 4σ 3 − σ 4
.
(6.12)
c'
22 − 32σ + 12σ 2
In particular we find that for an unstressed incompressible solid, c ' 21/22 = 0.9545 . . .,
which is a good approximation to the value found from the exact equation 0.9553. . . (Lord
Rayleigh, 1885). On Figure 12 we superpose the approximate and exact curves, and find
good agreement in the σ = −1 to σ = 1 region. Beyond that range, c departs from the
neighbourhood of 1, and a better choice must be made for the initial guess. This point
is not developed further here but it is clearly exposed in the paper by Mozhaev (1991).
1.0
c
0.9
0.8
0.7
0.6
σ
-1.5
-1.0
-0.5
0.5
1.0
1.5
Figure 12. Surface wave speed in an incompressible solid subject to hydrostatic pressure only.
Solid plot: exact calculation, dashed plot: approximate calculation.
6.4
Other waves
Taylor (1981) and Taziev (1989) both remarked that many of the roots to an explicit
secular equation in polynomial form have a physical origin, because they correspond
35
to a motion satisfying both the equations of motion and the boundary conditions at
the interface. In particular for the solid/vacuum interface, the polynomial provides the
parameters not only for the Rayleigh surface wave, but also for pseudo-surface waves,
for Brewster reflection of bulk waves, for the reflection of inhomogeneous waves, for
“organ-pipe” modes, etc.
One of the most important root is that corresponding to a leaky surface wave (if it
exists), whose energy is diffused slowly into the half-space, and for which the classical
methods evoked in Section 6.1 run into major difficulties.
Another application for this wealth of information is that many more wavefronts can
be found, in particular those observed in the neighbourhood of cusps, which correspond
to the interference of inhomogeneous plane waves (Huet, 2006). Figure 13 shows the
remarkable correspondence obtained between experiment and predictions, for wavefronts
propagating on the surface of a copper sample (cubic symmetry). Here the explicit secular
equation proves useful because it is easy to differentiate, as required for the derivation
of the curve.
1
0.5
–1
–0.5
0
0.5
1
–0.5
–1
Figure 13. Wavefronts on the surface of copper cut along a symmetry plane. Left: experimental results obtained by LASER impact (Huet, 2006). Right: predictions derived from the
secular equation in polynomial form, showing in black, the Rayleigh wave wavefront and in grey,
the wavefronts due to leaky waves and interference of inhomogeneous waves.
Bibliography
D.M. Barnett. Bulk, surface, and interfacial waves in anisotropic linear elastic solids.
International Journal of Solids and Structures, 37:45–54, 2000.
D.M. Barnett and J. Lothe. Free surface (Rayleigh) waves in anisotropic elastic halfspaces: the surface impedance method. Proceedings of the Royal Society of London,
Series A, 402:135–152, 1985.
36
M.A. Biot. Surface instability of rubber in compression. Applied Science Research, Series
A, 12:168–182, 1963.
Ph. Boulanger and M. Hayes. Finite-amplitude waves in deformed Mooney-Rivlin material. Quarterly Journal of Mechanics and Applied Mathematics, 45:575–593, 1992.
P. Chadwick. Continuum Mechanics. Allen & Unwin, 1976.
P. Chadwick. Interfacial and surface waves in pre-strained isotropic elastic media. ZAMP,
46:S51–S71, 1995.
P. Chadwick. The application of the Stroh formalism to prestressed elastic media. Mathematics and Mechanics of Solids, 2:379–403, 1997.
P. Chadwick and D.A. Jarvis. Surface waves in a pre-stressed elastic body. Proceedings
of the Royal Society of London, Series A, 366:517–536, 1979a.
P. Chadwick and D.A. Jarvis. Interfacial waves in a pre-strained neo-Hookean body.
Quarterly Journal of Mechanics and Applied Mathematics, 32:387–399, 1979b.
B. Collet and M. Destrade. Explicit secular equations for piezoacoustic surface waves:
Shear-horizontal modes. Journal of the Acoustical Society of America, 116:3432–3442,
2004.
B. Collet and M. Destrade. Explicit secular equations for piezoacoustic surface waves:
Rayleigh modes. Journal of Applied Physics, 98:054903, 2005.
P. Connor and R.W. Ogden. The effect of shear on the propagation of elastic surface
waves. International Journal of Engineering Science, 33:3432–3442, 1995.
P. Connor and R.W. Ogden. The influence of shear strain and hydrostatic stress on
stability and elastic waves in a layer. International Journal of Engineering Science,
34:375–397, 1996.
P.K. Currie. The secular equation for Rayleigh waves on elastic crystals. Quarterly
Journal of Mechanics and Applied Mathematics, 32:163–173, 1979.
M. Destrade. Finite-amplitude inhomogeneous plane waves in a deformed Mooney-Rivlin
material. Quarterly Journal of Mechanics and Applied Mathematics, 53:343–361, 2000.
M. Destrade. Small-amplitude inhomogeneous plane waves in a deformed Mooney-Rivlin
material. Quarterly Journal of Mechanics and Applied Mathematics, 55:109–126, 2002.
M. Destrade. Elastic interface acoustic waves in twinned crystals. International Journal
of Solids and Structures, 40:7375–7383, 2003.
M. Destrade. On interface waves in misoriented pre-stressed incompressible elastic solids.
IMA Journal of Applied Mathematics, 70:3–14, 2005.
M. Destrade and Y.B. Fu. The speed of interfacial waves polarized in a symmetry plane.
International Journal of Engineering Science, 44:26–36, 2006.
M. Destrade and R.W. Ogden. Surface waves in a stretched and sheared incompressible
elastic material. International Journal of Non Linear Mechanics, 40:241–253, 2005.
M. Destrade, M. Otténio, A.V. Pichugin, and G.A. Rogerson. Non-principal surface
waves in deformed incompressible materials. International Journal of Engineering
Science, 43:1092–1106, 2005.
M.A. Dowaikh and R.W. Ogden. On surface waves and deformations in a pre-stressed
incompressible elastic solid. IMA Journal of Applied Mathematics, 44:261–284, 1990.
37
M.A. Dowaikh and R.W. Ogden. Interfacial waves and deformations in pre-stressed
elastic media. Proceedings of the Royal Society of London, Series A, 433:313–328,
1991.
J.N. Flavin. Surface waves in pre-stressed Mooney material. Quarterly Journal of Mechanics and Applied Mathematics, 16:441–449, 1963.
Y.B. Fu. Existence and uniqueness of edge waves in a generally anisotropic elastic plate.
Quarterly Journal of Mechanics and Applied Mathematics, 56:605–616, 2003.
Y.B. Fu. An explicit expression for the surface-impedance matrix of a generally
anisotropic incompressible elastic material in a state of plane strain. International
Journal of Non Linear Mechanics, 40:229–239, 2005a.
Y.B. Fu. An integral representation of the surface-impedance tensor for incompressible
elastic materials. Journal of Elasticity, 81:75–90, 2005b.
Y.B. Fu and D.W. Brookes. An explicit expression for the surface-impedance tensor of a
compressible monoclinic material in a state of plane strain. IMA Journal of Applied
Mathematics, 71:434–445, 2006.
Y.B. Fu and A. Mielke. A new identity for the surface impedance matrix and its application to the determination of surface-wave speeds. Proceedings of the Royal Society
of London, Series A, 458:2523–2543, 2002.
A.N. Gent. A new constitutive relation for rubber. Rubber Chemistry and Technology,
69:59–61, 1996.
A.N. Guz. Elastic waves in bodies with initial (residual) stresses. International Applied
Mechanics, 38:23–59, 2002.
M.A. Hayes and R.S. Rivlin. Surface waves in deformed elastic materials. Archives for
Rational Mechanics and Analysis, 8:358–380, 1961.
C.O. Horgan and G. Saccomandi. A description of arterial wall mechanics using limiting
chain extensibility constitutive models. Biomechanics Modeling in Mechanobiology,
1:251–266, 2003.
G. Huet. Fronts d’Ondes Ultrasonores à la Surface d’un Milieu Semi-Infini Anisotrope:
Théorie des Rayons Réels et Complexes. PhD Thesis, Université de Bordeaux 1, 2006.
W. Hussain and R.W. Ogden. Reflection and transmission of plane waves at a shear-twin
interface. International Journal of Engineering Science, 38:1789–1810, 2000.
K.A. Ingebrigsten and A. Tonning. Elastic surface waves in crystals. Physical Review,
184:942–951, 1969.
A.R. Karduna, H.R. Halerpin, and F.C.P. Yin. Experimental and numerical analyses of
indentation in finite-sized isotropic and anisotropic rubber-like materials. Annals of
Biomedical Engineering, 25:1009–1016, 1997.
J. Merodio and R.W. Ogden. Material instabilities in fiber-reinforced nonlinearly elastic
solids under plane deformation. Archives of Mechanics, 54:525–552, 2002.
V.G. Mozhaev. Approximate analytical expressions for the velocity of Rayleigh waves
in isotropic media and on the basal plane in high-symmetry crystals. Soviet Physics
Acoustics, 37:1009–1016, 1991.
V.G. Mozhaev. Some new ideas in the theory of surface acoustic waves in anisotropic
media. In D.F. Parker and A.H. England, editors, IUTAM Symposium on Anisotropy,
Inhomogeneity and Nonlinearity in Solid Mechanics, pages 455–462. Kluwer, 1995.
38
R.W. Ogden. Elements of the theory of finite elasticity. In Y.B. Fu and R.W. Ogden, editors, Nonlinear Elasticity: Theory and Applications, pages 1–58. Cambridge
University Press, 2001.
R.W. Ogden. List of publications. Mathematics and Mechanics of Solids, 8, 9:449–450,
3–4, 442–443, 2003, 2004.
Y.-H. Pao, W. Sachse, and H. Fukuoka. Acoustoelasticity and ultrasonic measurements
of residual stresses. In W.P. Mason and R.N. Thurston, editors, Physical Acoustics,
Vol. 17, pages 61–143. Academic Press, 1984.
A.V. Pichugin. Asymptotic Models for Long Wave Motion in a Pre-Stressed Incompressible Elastic Plate. PhD Thesis, University of Salford, 2001.
M.L. Raghavan and D.A. Vorp. Toward a biomechanical tool to evaluate rupture potential
of abdominal aortic aneurysm: identification of a finite strain constitutive model and
evaluation of its applicability. Journal of Biomechanics, 33:475–482, 2000.
Lord Rayleigh. On waves propagated along the plane surface of an elastic solid. Proceedings of the London Mathematical Society, 17:4–11, 1885.
G.A. Rogerson and K.J. Sandiford. Harmonic wave propagation along a non-principal
direction in a pre-stressed elastic plate. International Journal of Engineering Science,
37:1663–1691, 1999.
G. Saccomandi. Phenomenology of rubber-like materials, CISM Lecture Notes 452. In
G. Saccomandi and R.W. Ogden, editors, Mechanics and Thermomechanics of Rubberlike Solids, pages 91–134. Springer, 2004.
R. Stoneley. Elastic waves at the surface of separation of two solids. Proceedings of the
Royal Society of London, 106:416–428, 1924.
A.N. Stroh. Some analytic solutions for rayleigh waves in cubic crystals. Journal of
Mathematics and Physics, 41:77–103, 1962.
D.B. Taylor. Surface waves in anisotropic media: the secular equation and its numerical
solution. Proceedings Royal Society of London, Series A, 376:265–300, 1981.
D.B. Taylor and P.K. Currie. The secular equation for Rayleigh waves on elastic crystals
ii: corrections and additions. Quarterly Journal of Mechanics and Applied Mathematics, 34:231–234, 1981.
R.M. Taziev. Bipartial surface acoustic waves. Soviet Physics Acoustics, 33:100–103,
1987.
R.M. Taziev. Dispersion relation for acoustic waves in an anisotropic elastic half-space.
Soviet Physics Acoustics, 35:535–538, 1989.
I. Thompson, I.D. Abrahams, and A.N. Norris. On the existence of flexural edge waves on
thin orthotropic plates. Journal of the Acoustical Society of America, 112:1756–1765,
2002.
T.C.T. Ting. Anisotropic Elasticity: Theory and Applications. University Press, 1996.
T.C.T. Ting. The polarization vector and secular equation for surface waves in an
anisotropic half-space. International Journal of Solids and Structures, 41:2065–2083,
2004.
T.C.T. Ting. The polarization vectors at the interface and the secular equation for
Stoneley waves in monoclinic bimaterials. Proceedings of the Royal Society of London,
Series A, 461:711–731, 2005.
39
L.R.G. Treloar. The Physics of Rubber Elasticity. Clarendon Press, 1949.
A.J. Willson. Surface and plate waves in biaxially-stressed elastic media. Pure and
Applied Geophysics, 102:182–192, 1973a.
A.J. Willson. Surface waves in uniaxially-stressed Mooney material. Pure and Applied
Geophysics, 112:352–364, 1973b.
40
View publication stats
Keep reading this paper — and 50 million others — with a free Academia account
Used by leading Academics
Alberto Pimpinelli
Université Clermont Auvergne
MARKUS MUELLER
Paul Scherrer Institute
Viorel Chihaia
Institute of Physical Chemistry "Ilie Murgulescu"
Khaled Dine
University of Dr moulay taher
