t t

dm (giu— uig— ) dV (‚jσiju— ’ ui‚jσ— ).

= ’ +

dt dt

i i i ij

t0 V t0 V

(5.3.3)

The volume integral over the stress

dV(‚jσiju— ’ui‚jσ— ) = dV‚j(σiju— ’uiσ— ) ’ dV(σije— ’eijσ— ),

i ij i ij ij ij

V V V

(5.3.4)

and

σije— = cijklekle— = cklije— ekl = σ— ekl. (5.3.5)

ij ij ij kl

80 Linearised Elasticity and Viscoelasticity

Thus (5.3.3) can be written as

t

‚u—

‚ui —

ui ’ ui i

dVρ

‚t ‚t

V t0

t t

dm(giu— uig— ) dS(„iu— ’ ui„— ).

= ’ +

dt dt (5.3.6)

i i i i

t0 V t0 S

5.3.1 The reciprocal theorem

Take the time integration up to in¬nity and assume that u, u— are such that the

left-hand side of (5.3.6) vanishes, then

∞ ∞

dm giu— dS „iu— dm g— ui + dS „— u— ,

+ =

dt dt

i i i ii

’∞ V S ’∞ V S

(5.3.7)

and for a static ¬eld we can drop the time integral

dm giu— + dS „iu— = dm g— ui + dS „— u— . (5.3.8)

i i i ii

V S V S

This is Betti™s reciprocal theorem.

In applications of Betti™s reciprocal theorem (5.3.8) we take the required ¬eld to

be u and choose an appropriate starred ¬eld for which the solution is known and

which has as close a character as we can achieve to the desired ¬eld.

EXAMPLES: APPLICATIONS OF THE RECIPROCAL THEOREM

(i) In¬‚uence of self-gravitation

Self gravity is switched on in a homogeneous sphere “ ¬nd the decrease in radius:

At the surface r = a, „ = 0. Let g be the acceleration due to self gravity, then the

interior body force per unit mass is g = ’gx/a.

Choose u— = x implying uniform σ— = 3κ δij , then g— = 0 because the ¬eld is

ij

self-equilibrated. From (5.3.8)

gx

· xρ dV = 3κ 4πa2 u(a),

’

a

V

and hence

’gρ 4π 5

a = 3κ 4πa2 u(a),

a5

i.e., the decrease in radius is ’u(a) = ρga2 /(15κ). •

(ii) Loaded cylinder

A solid circular cylinder of height h and radius a is subjected to compressive loads L

applied via rigid plates in contact with its ends. Radial displacements over the plates are

opposed by a uniform frictional force „ (Figure 5.1). Find the reduction in height δh when

L is large enough to cause sliding everywhere over the ends.

5.3 Integral representations 81

L

„

A

x3

. h Figure 5.1. Con¬guration of loaded

a

cylinder and force distribution

L

Let x3 be measured axially from the centroid, and choose

u— = (…x1 , …x2 , ’x3 )p— /E,

corresponding to a uniaxial pressure p— , so g— = 0. Now also g = 0 in V, „— = „ = 0

over the curved surfaces. Thus, denoting an end by A,

„ · u— dA = „— · u dA,

the load L = ’ „3 dA, and so

a

p—

r„2πr dr = p— 1 πa2 δh.

1

2 HL ’… 2

E 0

Thus

h L 4…a„

δh = ’ .•

2

E πa 3h

5.3.2 The representation theorem

Consider a delta-function body force

gi = δikδ(x ’ ξ)δ(t ’ θ), (5.3.9)

acting in a volume V. The resulting displacement at x, t is speci¬ed by the Green™s

tensor with components

Gk(x, t; ξ, θ), (5.3.10)

i

82 Linearised Elasticity and Viscoelasticity

associated with the force in the k-direction at ξ, θ. Even when the properties of the

region vary Gk will depend on the time difference t ’ θ.

i

If we take gi as above and set

g— = δilδ(x ’ ζ)δ(σ ’ t), (5.3.11)

i

the displacement

u— = Gl(x, ’t; ζ, ’σ). (5.3.12)

i i

We assume that both Gk and Gl obey the same homogeneous boundary conditions

i i

on S, and then when we apply the reciprocal theorem (5.3.7) the surface integrals

vanish

∞

dV δikδ(x ’ ξ)δ(t ’ θ)Gl(x, ’t; ζ, ’σ)

dt i

’∞ V

∞

dV δikδilδ(x ’ ζ)δ(σ ’ t)Gk(x, t; ξ, θ),

= dt (5.3.13)

i

’∞ V

so that when ξ, ζ are in V

Gl (ξ, ’θ; ζ, ’σ) = Gk(ζ, σ; ξ, θ).

kξ lζ (5.3.14)

Equation (5.3.14) expresses Green™s tensor reciprocity.

Now set

g— = δikδ(x ’ ξ)δ(θ ’ t) (5.3.15)

i

in (5.3.6) and then

u— = Gk(x, ’t; ξ, ’θ), (5.3.16)

i i

and so

dVx{gi(x)Gk(x, ’t; ξ, ’θ) ’ ui(x)δ(x ’ ξ)δ(θ ’ t)}

0= dθ i

t0 V

dSx{„i(x)Gk(x, ’t; ξ, ’θ) ’ ui(x)Hk(x, ’t; ξ, ’θ)}

+ dθ i i

t0 S

dVx[‚tui(x)Gk ’ ui‚tGk]t=t0 ,

+ (5.3.17)

i i

V

where Hk is the traction associated with Gk.

i i

Now, using the reciprocity of the Green™s tensor we can obtain a representation

of the displacement as an integral over the boundaries, body forces and initial

conditions

dVx gi(x)Gi (ξ, t; x, s)

˜(ξ)uk(ξ, t) =

ξ ξ kξ

ds

t0 V

dSx{„i(x)Gi (ξ, t; x, s) ’ ui(x)Hi (ξ, t; x, s)}

+ kξ kξ

ds

t0 S

dVx[‚tui(x)Gi ’ ui‚tGi ]t=t0 .

+ (5.3.18)

k k

V

5.4 Elastic waves 83

where

˜(ξ) = 1,

ξ ξ in V,

= 0, otherwise .

5.4 Elastic waves

The equation of motion and the stress-strain relation for small disturbances in an

elastic medium, in the absence of body forces, can be combined to give

‚2ui

‚ ‚uk

cijkl =ρ 2 , (5.4.1)

‚xj ‚xl ‚t

as the governing differential equation for displacement u. Equation (5.4.1) controls

the spatial and temporal development of the displacement ¬eld and admits solutions

in the form of travelling waves. Consider then a plane wave travelling in the

anisotropic medium with frequency ω:

ui = Ui exp[iωpn · x ’ iωt]

(5.4.2)

= Ui exp[iωpnkxk ’ iωt].

n represents the direction of travel of the phase fronts and p the apparent slowness

(inverse of wave velocity) in that direction. On substituting this plane wave form

into (5.4.1) we obtain

ω2p2cijklnjnlUk = ’ω2ρUi, (5.4.3)

which constitutes an eigenvalue problem for the slowness p for waves travelling in

the direction n

[p2cijklnjnl ’ ρδik]Uk = 0. (5.4.4)

Equation (5.4.4) has three roots for p2 for each direction n, associated with

polarisations speci¬ed by the eigenvectors U(r)(n). The slownesses pr are

independent of frequency.

In a general anisotropic medium, the slownesses vary with direction and can give

quite complex slowness and wave surfaces as, e.g., those illustrated in Figure 5.2.

5.4.1 Isotropic media

For an isotropic medium, the eigenvalue equation for slowness takes the form

[p2(» + μ)nink + p2μδik ’ ρδik]Uk = 0, (5.4.5)

where we have used the isotropic form for the elastic modulus tensor cijkl.

(i) P waves

Consider a displacement ¬eld oriented along the propagation direction, i.e.,

uP = cn, (5.4.6)

84 Linearised Elasticity and Viscoelasticity

3 3

P SV

1

3

SH

1

Figure 5.2. Wave surfaces for a transversely isotropic (hexagonal) medium formed from

the envelope of plane wavefronts for a unit time of propagation.

then the slowness equation is

[(» + 2μ)p2 ’ ρ]ni = 0, (5.4.7)

so that we have a wave disturbance with slowness

a2 = ρ/(» + 2μ), (5.4.8)

and associated phase velocity

± = [(» + 2μ)/ρ]1/2. (5.4.9)

This longitudinal wave solution is called the P wave

uP = APn exp[iω(an · x ’ t)]. (5.4.10)

P

S

Figure 5.3. Elastic waves in a uniform isotropic medium: P waves with longitudinal mo-

tion and S waves with transverse (shear) motion.

5.4 Elastic waves 85

(ii) S waves

There is an alternative type of elastic wave motion in which the displacement is

transverse to the direction of motion. Consider a displacement ¬eld

uS = cs, s · n = 0,

with (5.4.11)

for which the slowness equation reduces to

[μp2 ’ ρ]si = 0, (5.4.12)

so that we have a wave disturbance with slowness

b2 = ρ/μ (5.4.13)

and associated phase velocity

β = [μ/ρ]1/2. (5.4.14)

This form of solution holds for any direction orthogonal to the direction of motion,

i.e., we have a degenerate eigenvalue problem from which we can choose two

orthogonal S wave vectors. It is conventional to choose one vector in the vertical

plane (denoted SV) and the other purely horizontal (denoted SH). This choice

simpli¬es the analysis of wave propagation in horizontally strati¬ed media, and

so we represent the S wave ¬eld as

uS = {BvsV + BHsH} exp[iω(bn · x ’ t)], (5.4.15)

where sV lies in the vertical plane through n and sH in the horizontal plane

[sV · n = sH · n = 0].

5.4.2 Green™s tensor for isotropic media

The wave disturbance at an observation point r produced by a delta-function point

force in the direction d at the origin in an unbounded isotropic medium can be

represented as

δ(t ’ r/±) δ(t ’ r/β)

4πu(r) = (d · ^)^ + ^ — (d — ^)

rr r r

r±2 rβ2

t

+(3(d · ^)^ ’ d) 3 [H(t ’ r/±) ’ H(t ’ r/β)],

rr (5.4.16)

r

where ^ is a unit vector in the direction r and H(t) is the Heaviside step function.

r

This result was ¬rst derived by Stokes in 1851 and can be deduced by many

different techniques, e.g., Fourier transform methods or the superposition of

potential solutions.

The disturbance consists of two delta functions in time, one longitudinal

travelling with the faster P-wave velocity ± and the other with transverse

displacement and the slower shear wave velocity β (this can be represented in

terms of SV and SH components if required). Between these two wavefronts there

86 Linearised Elasticity and Viscoelasticity

u

t

r /± r /β

Figure 5.4. The displacement from a point force in an unbounded medium with distinct P

and S arrivals linked by a ˜near-¬eld™ ramp disturbance.

is a ˜near-¬eld™ disturbance with intermediate character which decays more rapidly

with distance r (Figure 5.4).

The speci¬c form of the Green™s tensor Gk(x, t; 0, 0) for an observation point

i

with direction cosines γi is thus

t

Gk(x, t) = (3γiγj ’ δij) 3 [H(t ’ r/±) ’ H(t ’ r/β)]

i

r

1 1

+γiγj 2 δ(t ’ r/±) + (γiγj ’ δij) 2 δ(t ’ r/β). (5.4.17)

r± rβ

The radiation pattern for the far-¬eld terms which have a r’1 decay with distance

r away from the source takes a fairly simple form: the P waves depends on the

cosine of the inclination to the force direction and the S wave pattern on the sine of

the inclination (Figure 5.5).

More realistic sources can be obtained by using couples and dipoles whose

associated displacements can be derived by differentiation of the Green™s tensor.

The sense of differentiation is that of the displacement in the couple or dipole. The

double couple in Figure 5.5(c) is the equivalent force system for an in¬nitesimal

shear displacement on a surface aligned with one of the arms of the couple, and the

resulting radiation patterns are used in the interpretation of faulting mechanisms

from seismological observations.

5.4.3 Interfaces

The solution presented in the previous section is for elastic waves in an unbounded

medium, and more complicated behaviour results once material interfaces or a free

surface is introduced.

At a free surface, the requirement is that the traction should vanish. The surface

acts as a mirror for horizontally polarised S waves (SH); an incident plane SH wave

is re¬‚ected without change of amplitude and the SH component of the Green™s

5.4 Elastic waves 87

P S

(a)

1 1 1

2 2 2

3 3 3

(b)

1 1 1

2 2 2

3 3 3

(c)

1 1 1

2 2 2

3 3 3

Figure 5.5. Far-¬eld radiation patterns for P and S waves from (a) a single point force, (b)

a vertical dipole, (c) a double couple.

function within the half-space can be found by introducing a ¬ctitious mirrored

source lying above the free surface. The situation is more complex for vertically

polarised waves. An incident plane P or SV wave will be re¬‚ected from the free

surface with an equal inclination to the normal to the surface, in addition there will

be conversion between P and SV waves with the inclination of the converted wave

determined by Snell™s law

sin i sin j

= , (5.4.18)

±0 β0

where i is the inclination to the normal of the P wave with surface wavespeed ±0,

and j is the inclination of the SV wave with surface wavespeed β0. In addition

there is a special form of wave that satis¬es the vanishing traction condition called

a Rayleigh wave. This consists of coupled exponentially decaying P and S waves

travelling with a horizontal velocity about 0.9 of the shear wave velocity. A surface

source in a half-space excites the Rayleigh wave in addition to the P and SV body

88 Linearised Elasticity and Viscoelasticity

waves; an elegant treatment of the development of the Green™s functions for this

case is given by Hudson (1981).

At an interface between two different media the boundary conditions are that

(i) the displacement u is continuous,

(ii) the traction „(n) associated with the normal n to the interface is continuous.

These boundary conditions can only be satis¬ed by the coupling of P and S waves

at the interface. For a plane wave incident on a horizontal interface the P and SV

wave ¬elds are coupled, but SH is independent. Such incident plane waves are

both re¬‚ected and transmitted through the interface with inclinations to the vertical

dictated by Snell™s law, which may be viewed as a condition for the continuity of

phase along the interface.

A major application of linearised elastic wave theory is in seismology, and

a comprehensive development of seismic wave theory with applications to the

interpretation of observed seismograms can be found in Kennett (2001, 2002), see

also Chapter 11.

5.5 Linear viscoelasticity

As in (4.6.1) we assume that the stress tensor σ depends linearly on the history of

strain et so that

σij(x, t) = Cijpq{et }, (5.5.1)

pq

with ij, pq symmetry associated with the symmetries of the stress and strain tensor.

Creep and relaxation

Consider imposing a step-function stress cycle on a wire (i.e., a one-dimensional

problem)

σ(t) = σ0H(t), (5.5.2)

then the resulting strain for t > 0 can be written as

σ0

e(t) = [1 + ψ(t)], ψ(0) = 0, (5.5.3)

ED

where ED is the dynamic modulus and σ0/ED is an instantaneous elastic response.

The monotonically increasing function ψ(t) is known as the creep function. If ψ(t)

approaches a constant as t ’ ∞, then we have long-term solid behaviour, whereas

for materials like pitch whose ultimate behaviour is that of a (viscous) ¬‚uid ψ(t)

increases without limit.

We can extend this concept to a more general stress history, by considering the

superposition of many small steps and then taking the limit as the step interval tends

to zero, to obtain

t

E’1

e(t) = σ(t) + ds σ(s)ψ(t ’ s) ,

™ (5.5.4)

D

’∞

5.5 Linear viscoelasticity 89

linear viscoelastic liquid

ψ

E 0 e(t)

σ0 viscoelastic solid

1

elastic solid

Figure 5.6. Strain evolution as a function of time.

a convolution integral over stress rate. If the creep function is differentiable for

t ≥ 0, we may use the properties of the convolution to rewrite the strain as

t

™

E’1

e(t) = σ(t) + ds σ(s)ψ(t ’ s) . (5.5.5)

D

’∞

An alternative viewpoint is to consider a strain cycle for which

e(t) = e0H(t), (5.5.6)

with associated stress for t > 0

σ(t) = EDe0[1 ’ φ(t)], φ(0) = 0. (5.5.7)

φ(t) is termed the relaxation function and is a monotonically increasing function

with a monotonically decreasing slope. If the long-term behaviour φ(∞) = 0 we

have a liquid.

elastic

σ (t)

E0e0

partially relaxing (solid)

φ

completely relaxing (liquid)

t

Figure 5.7. Stress evolution as a function of time.

90 Linearised Elasticity and Viscoelasticity

For an arbitrary differentiable strain history we can represent the stress behaviour

as

t

σ(t) = ED e(t) ’ ds e(s)φ(t ’ s) .

™ (5.5.8)

’∞

Once again, if the relaxation function is differentiable for t ≥ 0, we can reorganise

the convolution to give

t

™

σ(t) = ED e(t) ’ ds e(s)φ(t ’ s) . (5.5.9)

’∞

Since (5.5.4), (5.5.8) have to represent the same mechanical behaviour we can

establish the relation between the creep and relaxation functions. Consider the

application of a constant unit stress, then

e(t) = E’1[1 + ψ(t)], (5.5.10)

D

and so from (5.5.8) for t > 0

t

™

1 = 1 + ψ(t) ’ dsψ(s)φ(t ’ s). (5.5.11)

0

Simplifying (5.5.11), we require

t

™

ψ(t) = ds ψ(s)φ(t ’ s), t > 0. (5.5.12)

0

Similarly, we ¬nd

t

™

φ(t) = ds φ(s)ψ(t ’ s), t > 0. (5.5.13)

0

These two coupled integral equations have their simplest solution in the Laplace

transform domain. We de¬ne, e.g.,

∞

¯ dt ψ(t)e’pt,

ψ(p) = (5.5.14)

0

and then (5.5.5) and (5.5.9) can be written as

¯

e(p) = E’1[1 + pψ(p)]σ(p),

¯ ¯ (5.5.15)

D

¯

σ(p) = ED[1 ’ pφ(p)]e(p).

¯ ¯ (5.5.16)

Since these relations are equivalent we require

¯ ¯

[1 + pψ(p)][1 ’ pφ(p)] = 1, (5.5.17)

which connects the Laplace transforms of ψ(t), φ(t). We can now de¬ne a

transform modulus

¯ ¯

E(p) = ED[1 ’ pφ(p)] = ED[1 + pψ(p)]’1, (5.5.18)

5.6 Viscoelastic behaviour 91

in terms of which the one-dimensional constitutive relation for linear viscoelasticity

can be written as

σ(p) = E(p)e(p).

¯ ¯ (5.5.19)

The formal analogy between the Laplace transform domain result (5.5.19) and the

equivalent elastic relation provides the formal basis for solving many viscoelastic

problems via the correspondence principle.

5.6 Viscoelastic behaviour

The relaxation function for the standard linear solid model introduced in Section

4.6 is

φ(t) = [1 ’ („R/„C)][1 ’ e’t/„R ], (5.6.1)

and the corresponding creep function is

ψ(t) = [(„C/„R) ’ 1][1 ’ e’t/„C ]. (5.6.2)

Although the standard linear solid model goes a long way towards providing a

description of a viscoelastic material it does not provide an adequate description

of any real material. More satisfactory representations are provided by joining

standard linear solid elements in series or parallel.

In the series case, the stress in each element is the same whilst the strain is the

sum of the strains in each element. The creep function is then a sum of terms of the

form (5.6.2)

f(„n)[1 ’ e’t/„n ],

ψ(t) = (5.6.3)

n

where f(„n) may be regarded as the contribution of a creep process with a

characteristic time „n.

For a set of parallel elements, the stress is cumulative and the relaxation function

can be built from a set of terms of the form (5.6.1)

g(„n)[1 ’ e’t/„n ].

φ(t) = (5.6.4)

n

In the most general case we can consider a spectrum of relaxation and creep

processes.

The creep function for the Burgers material takes the form

t

ψ(t) = ψU + ”ψ[1 ’ e’t/„] + , (5.6.5)

·M

where ψU is the reciprocal of the relevant instantaneous modulus, ”ψ describes the

extent of viscoelastic relaxation with a relaxation time „C and ·M is the viscosity

of the Newtonian viscous element. The equivalent Maxwell relaxation time for

the in¬‚uence of viscosity is „M = ·MψU. Once again more general behaviour

92 Linearised Elasticity and Viscoelasticity

can be produced by a sum or spectrum of creep and relaxation processes, and such

a combination of Burgers elements can be used as a versatile representation of

deformation in the mantle of the Earth.

Isotropic linear viscoelasticity

For an isotropic medium the tensor function Cijpq in (5.1.5) can be represented as

an isotropic tensor of the form

Cijpq = Lδijδpq + M(δipδjq + δiqδjp), (5.6.6)

where L and M are scalar functions, and so

σij(x, t) = δijL{et } + 2M{et }, (5.6.7)

kk ij

and we can separate the dependence of the dilatation component from the deviatoric

stress and strain.

By analogy with the equations of isotropic elasticity we can write the

stress“strain relation as

σij(x, t) = »0δijekk(x, t) + 2μ0eij(x, t)

t t

™ ™

+ ds δijekk(x, s)R»(t ’ s) + ds 2eij(x, s)Rμ(t ’ s), (5.6.8)

0 0

in terms of relaxation functions R», Rμ and a representation analogous to (5.5.9)

for the one-dimensional case.

An equivalent representation to (5.6.8) expressing the strain as a function of the

stress history can be made in terms of creep functions C», Cμ. These are related

to the relaxation functions via coupled integral equations similar to (5.5.13). The

correspondence principle in the Laplace transform domain can be extended to these

three-dimensional problems and provides a means of solving many viscoelastic

problems.

5.7 Damping of harmonic oscillations

We consider a harmonic strain cycle of the form

eij = Eije’iωt, (5.7.1)

and then the corresponding stress is

σijeiωt = »0δijEkk + 2μ0Eij

t

™ ™

ds [δijEkkR»(t ’ s) + 2EijRμ(t ’ s)]eiω(t’s),

+ (5.7.2)

’T

where t = ’T is the time at which the disturbance starts. Now let T ’ ∞ and

introduce the Fourier transforms of the time derivatives of the relaxation functions

∞ ∞

™ ™

iωt

dt Rμ(t)eiωt.

»1(ω) = dt R»(t)e , μ1(ω) = (5.7.3)

0 0

5.7 Damping of harmonic oscillations 93

We can express the stress as

σij = {[»0 + »1(ω)]δijEkk + 2[μ0 + μ1(ω)]Eij}e’iωt, (5.7.4)

with a stress“strain relation in the same form as for an elastic medium but with

frequency-dependent complex moduli

¯

» = »0 + »1(ω), μ = μ0 + μ1(ω).

¯ (5.7.5)

We may therefore employ the descriptions of P and S waves developed in section

5.4, but have to allow for the frequency dependent velocities and consequent energy

loss.

The rate of energy dissipation per unit volume due to deformation of the medium

is from (4.2.11)

‚vj

D

ρ U = σij ≈ σijeij,

™ (5.7.6)

Dt ‚xi

in the linear approximation. The average work dissipated per unit time and volume

is thus

¯^ ^ ij

¯ — iωt

= ’ 1 ω Im[κ””— + 2μEij E— ],

1

W= 2 Re[iωσijEije ] ¯ (5.7.7)

2

^

where the dilatation ” = Ekk, and the deviatoric strain Eij = Eij ’ 1 ”δij. The

3

complex bulk modulus

¯ ¯ 3¯

κ = » + 2 μ = κ0 + κ1(ω). (5.7.8)

A convenient measure of the rate of energy dissipation is the loss factor Q’1(ω)

which may be de¬ned as

¯

”W(ω)

Q’1 = , (5.7.9)

2πW0(ω)

¯

where ”W(ω) is the energy loss in a cycle at frequency ω and W0(ω) is the

˜elastic™ energy stored in the oscillation, i.e., the sum of the strain and kinetic

energy associated with just the instantaneous elastic moduli. For purely deviatoric

disturbances

^ ^ ij

’ 1 Im μ(ω)EijE—

¯ ’Im μ1(ω)

2

’1

Qμ (ω) = = . (5.7.10)

1 ^ ^— μ0

EijEij

2

Similarly for purely dilatational disturbances

’Im κ1(ω)

Q’1(ω) = ; (5.7.11)

κ

κ0

normally loss in pure dilatation is much less signi¬cant than loss in shear and

Q’1 Q’1.

κ μ

In general the loss factor for dilatation or shear can be represented as

’Im m1(ω)

Q’1(ω) = , (5.7.12)

m

m0

94 Linearised Elasticity and Viscoelasticity

where m is the appropriate complex modulus. Thus, for the standard linear solid

model,

„’1 ’ iω

C

m(ω) = m2 (5.7.13)

„’1 ’ iω

R

ω2 + „’1„’1 ’ iω(„’1 ’ „’1)

CR R C

= m2 , (5.7.14)

ω2 + „’2

R

and the instantaneous modulus is m2, and thus the loss factor

„’1 ’ „’1

Q’1(ω) R C

=ω . (5.7.15)

m

ω2 + „’2

R

By superposing a spectrum of relaxation times, a wide range of frequency

behaviour can be simulated.

The dissipation via Q’1(ω) is accompanied by a frequency dependent

correction to the real part of the complex modulus and hence to the apparent

wavespeed

ω2 + „’1„’1

CR

Re m(ω) = m2 . (5.7.16)

ω2 + „’2

R

This property is a consequence of a causal dissipation mechanism. Since the

™

relaxation contributions depend only on the past history of the strain Rm(t)

vanishes for t < 0, so that the transform m1(ω) must be analytic in the upper

half-plane. In consequence the real and imaginary parts of m1(ω) are Hilbert

transforms of each other,

∞ ∞

ω Q’1(ω )

1 Im m1(ω ) 2m0

Re m1(ω) = P =’ ,

dω P dω

ω 2 ’ ω2

π ω ’ω π

’∞ ’∞

(5.7.17)

where P denotes the Cauchy principal value of the integral.

6

Continua under Pressure

The linearised development of elasticity in the previous chapter can be applied

about a stress-free state and has therefore direct application to the shallower parts

of the Earth. As we shall see, we can make an incremental treatment about

the stress-state at depth, which will be dominantly hydrostatic, and derive the

incremental elastic moduli for small disturbances such as those created by the

passage of seismic waves.

It is dif¬cult in experimental studies to achieve the pressures and temperatures

appropriate to a direct study of the Earth. For the lower mantle, in particular,

recourse needs to be made to extrapolation from available conditions. A convenient

summary of the property of Earth materials at depth is provided by the equation of

state, which is commonly expressed as a relation between the pressure, the speci¬c

volume and the temperature. This simpli¬ed approach ignores non-hydrostatic

stresses and any non-isotropic response to the assumed hydrostatic pressure. Finite

strain representations are used in the Birch“Murnaghan equations to relate the

conditions at depths to the properties at zero pressure.

Long-term response to large non-hydrostatic stresses leads to irreversible

(plastic) deformation, and so over geological time materials are expected to have

yielded to such stresses, leaving just the hydrostatic load at depth. Non-hydrostatic

stresses are, however, important at shallow depth where they approach the

magnitude of the pressure and are responsible for the occurrence of earthquakes.

Although materials within the Earth are undoubtedly anisotropic on a

microscopic scale, the properties of the polycrystalline aggregates averaged over a

typical seismic wavelength are close to isotropic except towards the top and bottom

of the mantle.

6.1 Effect of radial strati¬cation

A convenient simpli¬cation when considering the properties of the Earth at depth

is to concentrate on the variations as a function of radius and to regard any

three-dimensional variation as a perturbation around this state. This enables the

use of a simple equation of state in terms of the hydrostatic pressure as a function

95

96 Continua under Pressure

of radius, which depends directly on knowledge of the density distribution within

the Earth.

6.1.1 Hydrostatic pressure

The predominant stress-¬eld within the Earth is imposed by the hydrostatic

pressure induced by gravitational effects. The gravitational force is directed

radially inwards and so g(r) = ’g^r, with the acceleration due to gravity g given

e

by

r

4πG

ds s2ρ(s),

g= 2 (6.1.1)

r 0

in terms of the radial distribution of density ρ(r). For hydrostatic balance, the

pressure is determined by

re

‚p

+ ρg = 0, i.e., p(r) = ds ρ(s)g(s), (6.1.2)

‚r r

with the integration taken from the free surface at re , where p(re ) = 0, down to

radius r.

16 400

14

12 300

g

10

[Mg/m3], [m/s2]

[GPa]

8 200

6

ρ

4 100

p

2

0 0

0 1000 2000 3000 4000 5000 6000

Depth [km]

Figure 6.1. Pressure p, gravitational acceleration g for the density distribution associated

with the AK135M model.

6.1 Effect of radial strati¬cation 97

The behaviour of the pressure and the gravitational acceleration within the Earth

is illustrated in Figure 6.1, using the density distribution ρ for the AK135 model

(Figure 1.6) modi¬ed to remove the low density zone in the upper mantle. Both the

pressure and the gravitational acceleration are continuous across the discontinuities

in density associated with the phase transitions in the upper mantle, the core“mantle

boundary and the boundary between the inner and outer core.

The gravitational acceleration grows through the core as we move from the centre

of the Earth outwards and reaches its maximum value at the core“mantle boundary,

because the mantle material is lighter and the increased volume associated with

spherical shells closer to the surface of the Earth does not compensate for the

lowered density. The secondary maximum at the base of the transition zone again

re¬‚ects the step down to lower densities in the outer layers.

The pressure increases steadily with depth with slight changes in slope at each of

the discontinuities in density. The pressure tends to a constant value at the centre of

the Earth of about 370 GPa; the reduction in slope for the inner core is associated

with the steady decrease of g towards zero as the centre of the Earth is approached.

6.1.2 Thermodynamic relations

Most equations of state can be regarded as statements of the conservation of energy.

In terms of the total energy U of a system, the ¬rst law of thermodynamics can be

written as

dU(S, V) = T dS ’ p dV, (6.1.3)

as a function of temperature T , entropy S, pressure p and volume V. Because

we want to consider the thermal state of a material it is most convenient from a

theoretical viewpoint to work in terms of temperature T and volume V, since these

can be related directly to the reference state, whereas stresses are most naturally

introduced in the deformed state (Chapter 3). The ¬rst law can then be written in

terms of the Helmholtz free energy F = U ’ T S:

dF(T, V) = ’S dT ’ p dV. (6.1.4)

In experiments on earth materials, it is usually the pressure rather than the volume

which is adjusted, and we have noted in section 4.2 that propagation is adiabatic

rather than isothermal. From (6.1.4)

‚F

p= . (6.1.5)

‚V T

To provide a key to the mineral physics literature we will use the notation K rather

than κ for the bulk modulus and G for the shear modulus instead of μ. The bulk

modulus (incompressibility) K is de¬ned as

‚p ‚p

KT = ’V , KS = ’V . (6.1.6)

‚V ‚V

T S

98 Continua under Pressure

KT applies to isothermal conditions in slow laboratory experiments and KS to the

adiabatic conditions in transient wave propagation.

A related quantity is the ˜seismic parameter™

KS ‚p

¦= = . (6.1.7)

ρ ‚ρ S

For an isotropic material, the bulk-sound speed

√

(±2 ’ 4 β2),

φ= ¦= (6.1.8)

3

where ± is the P wavespeed (5.4.9) (also designated vp) and β is the S wavespeed

(5.4.14) (also designated vs).

Regarding pressure p as a function of V, T , S

‚p ‚p ‚p ‚T

= ’ , (6.1.9)

‚V ‚V ‚T ‚V

T S V S

and so the difference between the isothermal and adiabatic moduli is

‚p ‚T

K T ’ KS = V . (6.1.10)

‚T ‚V

V S

The volume expansion coef¬cient ±th

1 ‚V ‚ ln V

±th = = ; (6.1.11)

V ‚T ‚T

p p

and so we can use the chain rule for partial derivatives to write

‚p 1 ‚V ‚p

=’ V = ±th KT . (6.1.12)

‚T V ‚T ‚V

V p T

The adiabatic derivative of temperature with respect to speci¬c volume can be

written as

‚T T

= ’γth , (6.1.13)

‚V V

S

where γth is the thermodynamic Gr¨ neisen parameter. The speci¬c heats at

u

constant volume CV and constant pressure Cp are given by

‚S ‚S

CV = T , Cp = T , (6.1.14)

‚T ‚T

V p

where S is the speci¬c entropy.

We are therefore able to express the Gr¨ neisen parameter in the form

u

±th VKS ±th VKT

γth = = . (6.1.15)

Cp CV

6.1 Effect of radial strati¬cation 99

Note also that

‚T T

= γth . (6.1.16)

‚p KS

S

In terms of these new variables we can express the difference between the

isothermal and adiabatic moduli as

KT ’ KS = ’±th KT Tγth , (6.1.17)

and so we ¬nd that

KS

= 1 + γth ±th T. (6.1.18)

KT

The behaviour of the bulk and shear moduli in the Earth with depth, inferred

from the properties of seismic waves, is illustrated in Figure 6.2 for the AK135

model.

16 16

14 14

12 12

K

10 10

[1011 Pa]

8 8

[1011 Pa]

6 6

K

4 4

G p

G

2 2

0 0

0 1000 2000 3000 4000 5000 6000

Depth [km]

Figure 6.2. The behaviour of the bulk modulus K and shear modulus G as a function of

depth for the AK135M model, compared to the variation in pressure p.

The effect of compression with depth plays a major role in determining the

physical properties, but the increase in temperature with depth also plays a

signi¬cant role. Seismic wave propagation will be close to isentropic conditions

and hence measurements are of the adiabatic modulus KS. For both the bulk of the

100 Continua under Pressure

lower mantle and the core, where the composition is expected to be homogeneous,

the bulk modulus K is approximately proportional to pressure p.

6.2 Finite strain deformation

As we have seen in Section 4.3, the elastic properties of a material can be directly

related to a local Lagrangian treatment about a suitable reference state. However,

until recently, the appropriate conditions for much of the Earth™s mantle were not

accessible to experimental studies. In particular, even where measurements were

carried out at high pressure, the temperature was close to ambient.

Estimates of the material behaviour in the Earth have therefore been made by

using a systematic extrapolation from the low-pressure regime of most experiments

to likely Earth conditions. A suitable equation of state for mineral systems can be

derived by examining the in¬‚uence of large ¬nite strain, suf¬cient to change the

volume of a unit cell from its zero-pressure value to that at depth. The aim of

such studies is to try to understand the distribution of physical properties revealed

from seismological work, as in Figure 6.2, in terms of the appropriate mineral

assemblages and temperature regimes.

We recall from Section 2.2.5 the de¬nition of the Green strain in the Lagrangian

frame

‚ui ‚uj ‚uk ‚uk

E = 1 (FT F ’ I), Eij = 1

2 ‚ξ + ‚ξ + ‚ξ ‚ξ , (6.2.1)

2

j i i j

and the Cauchy strain in the Eulerian frame

‚ui ‚uj ‚uk ‚uk

e = 1 (I ’ (FFT )’1), 1

eij = + ’ . (6.2.2)

2 2 ‚xj ‚xi ‚xi ‚xj

The volume transformation under ¬nite strain can be expressed either in terms

of the Green strain in a Lagrangian viewpoint,

2

V

= J2 = det[FT F] = det[I + 2E], (6.2.3)

V0

or alternatively in terms of the Cauchy strain in an Eulerian perspective,

2

V0

= J’2 = det[FT F]’1 = det[2e ’ I]. (6.2.4)

V

Consider hydrostatic pressure on a cubic or isotropic material. Then in the

Lagrangian formulation the Green strain tensor takes the form

det[I + 2E] = (1 + 2·)3,

Eij = ·δij, (6.2.5)

and so, from (6.2.3),

2/3 2/3

V ρ0

1 1

·= ’1 = ’1 ; (6.2.6)

2 2

V0 ρ

6.2 Finite strain deformation 101

this scalar · varies from 0 to ’ 1 as the pressure varies from 0 to ∞.

2

Alternatively, in the Eulerian formulation, the Cauchy strain is given by

det[2e ’ I] = (2 ’ 1)3,

eij = δij, (6.2.7)

and so

2/3 2/3

V0 ρ

’1 ’1

= ’1 = ’1 . (6.2.8)

2 2

V ρ0

In terms of the strain variables · and , the volume change

V0 ρ

= (1 + 2·)’3/2 = (1 ’ 2 )3/2.

= (6.2.9)

V ρ0

In an elastic material the strain energy W can be equated with the

speci¬c Helmholtz free energy F/ρ (we neglect thermal effects in either slow

near-isothermal deformations or the nearly isentropic incremental deformation

associated with the passage of seismic waves). The Helmholtz free energy can

then be expressed as a function of strain, e.g., through its dependence on the Green

strain E. With a generalisation of (6.1.4) for speci¬c quantities

ρ dW = ’ρS dT + ρ0σij dAij, (6.2.10)

in terms of the displacement gradient A = F ’ I (5.1.1). The Cauchy stress tensor

σij is thus

ρ ‚F

σij = ; (6.2.11)

ρ0 ‚Aij

which is equivalent to our previous expression (4.3.10)

‚W

det F σij = Fim Fjn, (6.2.12)

‚Emn

The effective elastic moduli Cijkl are de¬ned by

ρ ‚2W

‚σij

Cijkl = = ’ σijδkl, (6.2.13)

‚Akl ρ0 ‚Aij‚Akl

since ‚ρ/‚Akl = ’ρδkl. The differentiation may be carried out either isothermally

or adiabatically to give the appropriate moduli. In terms of the dependence on

strain, the effective moduli in the case of hydrostatic pre-stress take the form

‚2W

FjnFlq ’ pδkl,

det F Cijkl = FimFkp (6.2.14)

ij

‚Emn‚Epq

where

δkl = δikδjl + δilδjk ’ δijδkl. (6.2.15)

ij

102 Continua under Pressure

6.3 Expansion of Helmholtz free energy and equations of state

A commonly used approach to handling ¬nite strain effects is to expand the

Helmholtz free energy F as a Taylor series in the elastic strain from the

zero-pressure state. This yields a power series expansion in terms of a suitable

measure of ¬nite strain.

The level of strain required to attain the required density at depth in the Earth

is then employed to determine the pressure and the other desired properties of the

medium in the strained, i.e. pressured, state. For the Lagrangian viewpoint we

make a power-series expansion in the strain · from (6.2.6),

∞

Ai·i.

F = V0 (6.3.1)

i=0

Whereas for the Eulerian viewpoint (Birch“Murnaghan) the expansion is in terms

of the strain from (6.2.8),

∞

Bi i.

F = V0 (6.3.2)

i=0

The isothermal thermodynamic relations

‚F ‚p ‚K

p=’ , KT = ’V , KT = (6.3.3)

‚V ‚V ‚T

T T S

enable us to identify the coef¬cients in these series expansions.

For the Lagrangian expansion the coef¬cients in the expansion (6.3.1) for the

Helmholtz free energy F = V0 ∞ Ai·i are

i=0

A0 = A1 = 0, (6.3.4)

A2 = 9 K0, (6.3.5)

2

A3 = ’ 9 K0K0, (6.3.6)

2

A4 = 9 K0[K0K0 + K0(K0 + 1) ’ 1 ], (6.3.7)

2 9

where the subscript 0 indicates evaluation at zero pressure and ambient temperature

T0. The Lagrangian equation of state is then

’1/3 ∞

‚F ‚· V

(i + 1)Ai+1·i,

= ’1

pL = ’ (6.3.8)

3

‚· ‚V V0

i=1

since

’1/3

‚· V

= 1 (1 + 2·)’1/2 = 1

. (6.3.9)

3 3

‚V V0

Thus the pressure pL is represented as

’1/3

V

· ’ 3 K0·2 + · · · .

pL = ’3K0 (6.3.10)

2

V0

6.3 Expansion of Helmholtz free energy and equations of state 103

and the isothermal modulus can be found in terms of · from (6.1.6).

For the Eulerian viewpoint the coef¬cients in the expansion (6.3.2) for the

Helmholtz free energy F = V0 ∞ Bi i can be expressed as

i=0

B0 = B1 = 0, (6.3.11)

B2 = A2 = 9 K0, (6.3.12)

2

B3 = ’ 9 K0(K0 ’ 4), (6.3.13)

2

B4 = 9 K0[K0K0 + K0(K0 ’ 7) ’ 143

9 ], (6.3.14)

2

and the equation of state for the pressure takes the form

5/3

V 2

’ 3 (K0 ’ 4)

pE = ’3K0 + ··· . (6.3.15)

2

V0

When we reinstate the full dependence on by expressing the volume change in

terms of the strain, we ¬nd

pE = ’(1 ’ 2 )5/2 C1 ’ C2 2/2 + C3 3/6 , (6.3.16)

where

C1 = 2 B2 = 3K0, C2 = ’2B3, C3 = 6B4. (6.3.17)

3

The corresponding expression for the isothermal bulk modulus is

KT = 1 (1 ’ 2 )5/2 C1 + (C2 ’ 7C1) + (C3 ’ 9C2) 2 3

’ 11C3 (6.3.18)

3

= (1 ’ 2 )5/2 K0 + K1 + K2 2

+ ··· . (6.3.19)

For moderate changes in volume 1.0 > V/V0 > 0.9 the predictions of the

Lagrangian and Eulerian equations of state are essentially equivalent, but for larger

strains the pressure estimate pE increases more rapidly than pL.

Following Birch (1952) the preferred formulation for studies of the Earth is that

based on the Eulerian viewpoint because the in¬‚uence of the third-order term in

strain is much reduced. We note that the third-order coef¬cient B3 depends on the

pressure derivative of the bulk modulus through K0 ’ 4. For most minerals the

pressure derivative K0 lies in the range 3 to 6 and is frequently close to 4 and so

B3, C2 are small, with only a limited in¬‚uence from uncertainties in measurement.

The third order coef¬cient in the equivalent Lagrangian expression depends solely

on K0, and has a stronger in¬‚uence.

The third-order truncation of the Birch“Murnaghan equation of state has the

advantage that no knowledge is required of second pressure derivatives of physical

properties that are more dif¬cult to measure than the ¬rst derivatives, and hence are

rather poorly known. The fourth-order representation is nevertheless more ¬‚exible

because of the larger number of parameters.

104 Continua under Pressure

The isothermal seismic parameter

KT

¦T = (6.3.20)

ρ

1

(1 ’ 2 ) C1(1 ’ 7 ) + C2( ’ 9 2/2) + C3( 2

’ 11 3/6) ;

=

3ρ0

the grouping of the terms emphasises that the coef¬cients of the powers of are

not all independent.

With the aid of an Eulerian representation comparable to (6.2.14) we can develop

representations for the wavespeeds of both P and S waves in the form

1

±2 = v2 = 2

(1 ’ 2 ) L0 + L1 + L2 + ··· , (6.3.21)

P

ρ0

1

β2 = v2 = 2

(1 ’ 2 ) G0 + G1 + G2 + ··· , (6.3.22)

S

ρ0

where

L0 = K0 + 4 G0, L1 = K1 + 4 G1, etc. (6.3.23)

3 3

The ¬rst and second pressure derivatives of the shear and bulk moduli at zero

pressure are then given by

G0 = 1 (5G1 ’ G2)/K0,

3

K0 = 1 (5L1 ’ L2)/K0 ’ 4 G0,

3 3

(6.3.24)

5)G1]/K2,

1

G0 = 9 [5G3 + 3(K0 ’ 4)G2 ’ 5(3K0 ’ 0

K0 = 1 [5L3 + 3(K0 ’ 4)L2 ’ 5(3K0 ’ 5)L1]/K2 ’ 4 G0.

0

9 3

The coef¬cients for the contributions from the bulk modulus can be related to

quantities which can be measured, or estimated, from mineral physics experiments.

However, the corresponding results for the shear modulus are rather poorly

known, it is only recently that measurements of shear properties under conditions

appropriate to the Earth™s transition zone have been accomplished and the

derivatives are weakly constrained. Ab initio quantum mechanical calculations of

mechanical properties are also easier to perform for the bulk modulus, because

allowance does not need to made for anelastic thermal effects (see Chapter 9).

The ¬nite strain theory does not, by itself, provide a relation between the bulk

and shear moduli. Current seismological models are consistent with a dependence

of the shear modulus G on pressure p and bulk modulus K of the form

G = aK ’ bp. (6.3.25)

This behaviour is ilustrated in Figure 6.3 for the AK135 model, by the linear

segments associated with different parts of the mantle. The slopes are quite

similar, except for the zone just below the 660 km discontinuity and the highly

heterogeneous D region at the base of the mantle.

6.4 Incremental stress and strain 105

12

10

210

8

410

6