Reciprocity theorems for diﬀusion,ﬂow and waves
∗
Kees Wapenaar and Jacob Fokkema
Department of Geotechnology,Delft University of Technology,
P.O.Box 5048,2600 GA Delft,The Netherlands
(Dated:May 20,2010)
Diﬀusion,ﬂow and wave phenomena can each be captured by a uniﬁed diﬀerential equation in
matrixvector form.This equation forms the basis for the derivation of uniﬁed reciprocity theorems
for diﬀusion,ﬂow and wave phenomena.
PACS numbers:
I.INTRODUCTION
Diﬀusion,ﬂow and wave phenomena can each be cap
tured by the following diﬀerential equation in matrix
vector form
A
Du
Dt
+Bu +D
x
u = s,(1)
where u = u(x,t) is a vector containing space and time
dependent ﬁeld quantities,s = s(x,t) is a source vec
tor,A = A(x) and B = B(x) are matrices containing
spacedependent material parameters and D
x
is a ma
trix containing the spatial diﬀerential operators ∂/∂x
1
,
∂/∂x
2
and ∂/∂x
3
.Finally,D/Dt denotes the material
time derivative,deﬁned as
D
Dt
=
∂
∂t
+v¢ ∇=
∂
∂t
+v
k
∂
∂x
k
,
where ∂/∂t denotes the time derivative in the reference
frame and v = v(x) is the spacedependent ﬂow veloc
ity of the material;v
k
denotes the kth component of v.
Throughout this paper the summation convention ap
plies to repeated subscripts;lowercase Latin subscripts
run from1 to 3.The vectors and matrices in equation (1)
are further deﬁned in the appendices for diﬀusion (Ap
pendix A),acoustic wave propagation in moving ﬂuids
(Appendix B),momentumtransport (Appendix C),elas
todynamic wave propagation in solids (Appendix D) elec
tromagnetic diﬀusion and wave propagation (Appendix
E),and coupled elastodynamic and electromagnetic wave
propagation in porous solids (Appendix F).In this pa
per we use equation (1) as the basis for deriving uniﬁed
reciprocity theorems for these phenomena.In general,a
reciprocity theorem interrelates the quantities that char
acterize two admissible physical states that could occur
in one and the same domain [1].One can distinguish be
tween convolution type and correlation type reciprocity
theorems [2].Generally speaking,these two types of reci
procity theorems ﬁnd their applications in forward and
inverse problems,respectively.Both types of reciprocity
theorems will be derived for the ﬁeld vector u.
∗
This paper was originally published in 2004 in A.S.M.E.Journal
of Applied Mechanics,Vol.71,145150.This new version con
tains improved appendices,but the body of the paper remained
unchanged (except that the references to the appendices have been
updated and some printing errors have been corrected).
II.THE DIFFERENTIAL EQUATION IN THE
FREQUENCY DOMAIN
Reciprocity theorems can be derived in the time do
main,the Laplace domain and the frequency domain [3].
Here we only consider the frequency domain.We deﬁne
the Fourier transform of a timedependent function f(t)
as
ˆ
f(ω) =
R
∞
−∞
f(t) exp(−jωt)dt,where j is the imagi
nary unit and ω denotes the angular frequency.We apply
the Fourier transform to all terms in equation (1),under
the assumption that this equation is linear in u.Hence,
we only consider those cases in which the ﬁeld quantities
in u do not appear in any of the matrices or operators in
equation (1).In particular,this is why the term Du/Dt
in the momentum transport equation (C15) is replaced
by ∂u/∂t in (C20).Transforming equation (1) to the
frequency domain yields
A
³
jω +v ¢ ∇
´
ˆu +Bˆu +D
x
ˆu =ˆs,(2)
where ˆu = ˆu(x,ω) is the space and frequency dependent
ﬁeld vector and ˆs = ˆs(x,ω) is the space and frequency
dependent source vector.The term v ¢ ∇ should be
dropped for linearized momentum transport (Appendix
C) as well as for wave phenomena in nonmoving me
dia (Appendices D − F).Finally we remark that in a
number of cases matrix B contains temporal convolution
kernels in the time domain (Appendix B) or,equivalently,
complex frequencydependent material parameters in the
frequency domain (Appendices B and F).
III.MODIFICATION OF GAUSS’
DIVERGENCE THEOREM
The reciprocity theorem will be derived for a volume
V enclosed by surface ∂V with outward pointing normal
vector n.Note that ∂V does not necessarily coincide with
a physical boundary.Gauss’ divergence theorem plays a
central role in the derivation.For a scalar ﬁeld a(x),this
theorem reads
Z
V
∂a(x)
∂x
i
d
3
x =
I
∂V
a(x)n
i
d
2
x,(3)
where n
i
denotes the ith component of n.In this section
we will modify this theorem for the diﬀerential operator
2
matrix D
x
appearing in equations (1) and (2).Note that
D
x
= D
Tx
for all forms of D
x
appearing in the appendices
(here superscript
T
denotes matrix transposition only;it
does not denote operator transposition).Let D
IJ
denote
the operator in row I and column J of matrix D
x
.The
symmetry of D
x
implies D
IJ
= D
JI
.We deﬁne a ma
trix N
x
which contains the components of the normal
vector n,organized in a similar way as matrix D
x
,see
the appendices for details.Hence,if N
IJ
denotes the
element in row I and column J of matrix N
x
,we have
N
IJ
= N
JI
.For example,for matrices D
x
and N
x
in
equations (A4) and (A6) we have D
12
= D
21
= ∂/∂x
1
and N
12
= N
21
= n
1
.If we now replace the scalar ﬁeld
a(x) by a
I
(x)b
J
(x),we may generalize equation (3) to
Z
V
D
IJ
[a
I
(x)b
J
(x)] d
3
x =
I
∂V
a
I
(x)b
J
(x)N
IJ
d
2
x,(4)
where the summation convention applies for repeated
capital Latin subscripts (which may run from 1 to 4,6,9
or 19,depending on the choice of operator D
x
).Applying
the product rule for diﬀerentiation and using the symme
try property D
IJ
= D
JI
,we obtain for the integrand in
the lefthand side of equation (4)
D
IJ
(a
I
b
J
) = a
I
D
IJ
b
J
+(D
JI
a
I
)b
J
= a
T
D
x
b +(D
x
a)
T
b,(5)
where a and b are vector functions,containing the scalar
functions a
I
(x) and b
J
(x),respectively.Rewriting the
integrand in the righthand side of equation (4) in a sim
ilar way,we thus obtain
Z
V
[a
T
D
x
b +(D
x
a)
T
b] d
3
x =
I
∂V
a
T
N
x
bd
2
x.(6)
Finally we consider a variant of this equation.We replace
a by Ka,where Kis a diagonal matrix with the following
property
D
x
K= −KD
x
,(7)
see the appendices for details.With this replacement,
equation (6) becomes
Z
V
[a
T
KD
x
b −(D
x
a)
T
Kb] d
3
x =
I
∂V
a
T
KN
x
bd
2
x.
(8)
IV.RECIPROCITY THEOREM OF THE
CONVOLUTION TYPE
We consider two physical states in volume V.The
ﬁeld quantities,the material parameters,the ﬂow veloc
ity as well as the source functions may be diﬀerent in both
states and they will be distinguished with subscripts A
and B (of course the summation convention does not ap
ply for these subscripts).We substitute ˆu
A
and ˆu
B
for
a and b in equation (8),apply equation (2) for states A
and B and use the symmetry properties
A
T
A
K= KA
A
and B
T
A
K= KB
A
(9)
(see the appendices).This yields
I
∂V
ˆu
T
A
KN
x
ˆu
B
d
2
x =
Z
V
h
ˆu
T
A
Kˆs
B
−ˆs
T
A
Kˆu
B
i
d
3
x
+
Z
V
ˆu
T
A
K
£
jω(A
A
−A
B
) +(B
A
−B
B
)
¤
ˆu
B
d
3
x +
Z
V
h
¡
(v
A
¢ ∇)ˆu
A
¢
T
KA
A
− ˆu
T
A
KA
B
(v
B
¢ ∇)
i
ˆu
B
d
3
x.
(10)
The ﬁrst term in the last integral can be written as
¡
(v
A
¢ ∇)ˆu
A
¢
T
KA
A
ˆu
B
= ∇¢
¡
v
A
ˆu
T
A
KA
A
ˆu
B
¢
−
ˆ
u
T
A
K
∂(v
i,A
A
A
)
∂x
i
ˆ
u
B
−
ˆ
u
T
A
KA
A
(v
A
¢ ∇)
ˆ
u
B
.(11)
For diﬀusion (Appendix A) the term∂(v
i,A
A
A
)/∂x
i
van
ishes on account of the equation of continuity.For
acoustic wave propagation this term (with v
i,A
replaced
by v
0
i,A
,Appendix B) is negligible in comparison with
the spatial derivatives of the wave ﬁelds ˆu
A
and ˆu
B
.
For the other situations considered in the appendices,
v
i,A
is taken equal to zero.Hence,the term contain
ing ∂(v
i,A
A
A
)/∂x
i
will be dropped.Substituting the
remainder of the righthand side of equation (11) into
equation (10) and applying the theorem of Gauss for the
term containing the divergence operator,yields
I
∂V
ˆ
u
T
A
KN
x
ˆ
u
B
d
2
x =
Z
V
h
ˆ
u
T
A
K
ˆ
s
B
−
ˆ
s
T
A
K
ˆ
u
B
i
d
3
x
+
Z
V
ˆu
T
A
K
£
jω(A
A
−A
B
) +(B
A
−B
B
)
¤
ˆu
B
d
3
x
−
Z
V
ˆu
T
A
K
£
A
A
(v
A
¢ ∇) +A
B
(v
B
¢ ∇)
¤
ˆu
B
d
3
x
+
I
∂V
¡
ˆu
T
A
KA
A
ˆu
B
¢
v
A
¢ nd
2
x.(12)
This is the uniﬁed reciprocity theorem of the convolution
type (we speak of convolution type,because the multipli
cations in the frequency domain correspond to convolu
tions in the time domain).It interrelates the ﬁeld quan
tities (contained in ˆu
A
and ˆu
B
),the material parameters
(contained in A
A
,B
A
,A
B
and B
B
),the ﬂow velocities
3
(v
A
and v
B
) as well as the source functions (contained
in
ˆ
s
A
and
ˆ
s
B
) of states A and B.The lefthand side is a
boundary integral which contains a speciﬁc combination
of the ﬁeld quantities of states A and B at the bound
ary of the volume V.The ﬁrst integral on the righthand
side interrelates the ﬁeld quantities and the source func
tions in V.The second integral contains the diﬀerences
of the medium parameters in both states;obviously this
integral vanishes when the medium parameters in both
states are identical.The third integral on the righthand
side contains the ﬂow velocities in V;this integral van
ishes when the medium parameters in both states are
identical and the ﬂow velocities in both states are oppo
site to each other.The last integral on the righthand
side is a boundary integral containing the normal com
ponent of the ﬂow velocity in state A;it vanishes when
this ﬂow velocity is tangential to the boundary ∂V.De
pending on the type of application,states A and B can
be both physical states,or both mathematical states (e.g.
Green’s states),or one can be a physical state and the
other a mathematical state (the latter situation leads to
representation integrals).For further discussions on con
volution type reciprocity theorems in diﬀerent ﬁelds of
application we refer to Lyamshev [4],De Hoop and Stam
[1],Fokkema and Van den Berg [3],Allard et al.[5],
Pride and Haartsen [6] and Belinskiy [7].
V.RECIPROCITY THEOREM OF THE
CORRELATION TYPE
We substitute
ˆ
u
∗
A
and
ˆ
u
B
for a and b in equation
(6),where
∗
denotes complex conjugation.Following the
same procedure as in the previous section,using the sym
metry property
A
H
A
= A
A
,(13)
where
H
denotes complex conjugation and transposition,
we obtain
I
∂V
ˆu
H
A
N
x
ˆu
B
d
2
x =
Z
V
h
ˆu
H
A
ˆs
B
+ˆs
H
A
ˆu
B
i
d
3
x
+
Z
V
ˆu
H
A
£
jω(A
A
−A
B
) −(B
H
A
+B
B
)
¤
ˆu
B
d
3
x
+
Z
V
ˆu
H
A
£
A
A
(v
A
¢ ∇) −A
B
(v
B
¢ ∇)
¤
ˆu
B
d
3
x
−
I
∂V
¡
ˆ
u
H
A
A
A
ˆ
u
B
¢
v
A
¢ nd
2
x.(14)
This is the uniﬁed reciprocity theorem of the correlation
type (we speak of correlation type,because the multi
plications in the frequency domain correspond to cor
relations in the time domain).The term ˆu
H
A
contains
‘backpropagating’ ﬁeld quantities in state A,[2].When
we compare this reciprocity theorem with equation (12),
we observe that,apart from the complex conjugation,
the diagonal matrix K is absent in all integrals and
that some plus and minus signs have been changed.In
particular,the term (B
A
− B
B
) has been replaced by
(B
H
A
+ B
B
),which means that the second integral on
the righthand side no longer vanishes when the medium
parameters in both states are identical.Moreover,the
term
£
A
A
(v
A
¢ ∇) +A
B
(v
B
¢ ∇)
¤
has been replaced by
£
A
A
(v
A
¢ ∇) −A
B
(v
B
¢ ∇)
¤
,which means that the third
integral on the righthand side vanishes when the medium
parameters contained in matrix A as well as the ﬂow ve
locities are identical in both states.For a discussion on
the application of correlation type reciprocity theorems
to inverse problems we refer to Fisher and Langenberg
[8] and De Hoop and Stam [1].
VI.CONCLUSIONS
We have formulated a general diﬀerential equation in
matrixvector form [equation (1)],which applies to diﬀu
sion (Appendix A),acoustic wave propagation in moving
ﬂuids (Appendix B),momentum transport (Appendix
C),elastodynamic wave propagation in solids (Appendix
D) electromagnetic diﬀusion and wave propagation (Ap
pendix E) and coupled elastodynamic and electromag
netic wave propagation in ﬂuidsaturated porous solids
(Appendix F).For linear phenomena (which excludes
nonlinear momentum transport) we have transformed
the general equation from the time domain to the fre
quency domain [equation (2)].Based on this general
equation as well as the symmetry properties (7),(9) and
(13) we have derived uniﬁed reciprocity theorems of the
convolution type [equation (12)] and of the correlation
type [equation (14)],respectively.
Acknowledgements
We thank our colleague David Smeulders for his advice
on the theoretical aspects of wave propagation in porous
media.This work is part of the research programme
of the Netherlands research centre for Integrated Solid
Earth Science (ISES).
APPENDIX A:MASS DIFFUSION
The equation of continuity for species k in a mixture
of ﬂuids reads
ρ
DY
(k)
Dt
+
∂J
(k)
j
∂x
j
= ˙ω
(k)
,(A1)
where Y
(k)
is the mass fraction of species k,J
(k)
j
its mass
ﬂux relative to the mixture,ρ is the mass density of the
mixture and ˙ω
(k)
the mass production rate density of
4
species k (due to chemical reactions).Fick’s ﬁrst law of
diﬀusion reads
J
(k)
j
+ρD
(k)
∂Y
(k)
∂x
j
= 0,(A2)
where D
(k)
is the diﬀusion coeﬃcient for species k.Elimi
nating J
(k)
j
fromequations (A1) and (A2) yields the mass
diﬀusion equation,according to
ρ
DY
(k)
Dt
−
∂
∂x
j
³
ρD
(k)
∂Y
(k)
∂x
j
´
= ˙ω
(k)
.(A3)
On the other hand,equations (A1) and (A2) can be com
bined to yield equation (1),with
u =
Y
(k)
J
(k)
1
J
(k)
2
J
(k)
3
,s =
˙ω
(k)
0
0
0
,D
x
=
0
∂
∂x
1
∂
∂x
2
∂
∂x
3
∂
∂x
1
0 0 0
∂
∂x
2
0 0 0
∂
∂x
3
0 0 0
,
(A4)
A=
ρ 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
and B =
0 0 0 0
0
1
ρD
(k)
0 0
0 0
1
ρD
(k)
0
0 0 0
1
ρD
(k)
.
(A5)
Matrices N
x
and K,appearing in the modiﬁed diver
gence theorems (6) and (8),read
N
x
=
0 n
1
n
2
n
3
n
1
0 0 0
n
2
0 0 0
n
3
0 0 0
and K=
1 0 0 0
0 −1 0 0
0 0 −1 0
0 0 0 −1
.
(A6)
Note that matrices D
x
,A and B obey equations (7),(9)
and (13).Other diﬀusion phenomena can be formulated
in a similar way.
APPENDIX B:ACOUSTIC WAVE
PROPAGATION IN A MOVING FLUID
The linearized equation of motion in a moving ﬂuid
reads
ρ
Dv
i
Dt
+b
v
∗ v
i
+
∂p
∂x
i
= f
i
,(B1)
with
D
Dt
=
∂
∂t
+v
0
k
∂
∂x
k
,(B2)
where p is the acoustic pressure,v
i
the particle veloc
ity associated to the acoustic wave motion (which is to
be distinguished fromthe ﬂow velocity v
0
k
in the operator
D/Dt),ρ the mass density of the mediumin equilibrium,
f
i
the volume density of external force and b
v
a causal
loss function (∗ denotes a temporal convolution).The
spatial variations of the ﬂow velocity are assumed small
in comparison with those of the particle velocity of the
acoustic wave ﬁeld,i.e.∂v
0
i
/∂x
j
<< ∂v
i
/∂x
j
(this as
sumption can be relaxed,but then the equations become
more involved [9]).The linearized stressstrain relation
reads
κ
Dp
Dt
+b
p
∗ p +
∂v
i
∂x
i
= q,(B3)
where κ is the compressibility,q the volume injection rate
density and b
p
a causal loss function.
Eliminating v
i
from equations (B1) and (B3) for the
lossless situation (b
v
= b
p
= 0),yields the acoustic wave
equation,according to
D
Dt
³
κ
Dp
Dt
´
−
∂
∂x
i
³
1
ρ
∂p
∂x
i
´
= −
∂
∂x
i
³
f
i
ρ
´
+
Dq
Dt
.(B4)
On the other hand,equations (B1) and (B3) can be com
bined to yield the general matrixvector equation (1),
with D/Dt deﬁned in equation (B2) and
u =
p
v
1
v
2
v
3
,s =
q
f
1
f
2
f
3
,A=
κ 0 0 0
0 ρ 0 0
0 0 ρ 0
0 0 0 ρ
and B =
b
p
∗ 0 0 0
0 b
v
∗ 0 0
0 0 b
v
∗ 0
0 0 0 b
v
∗
.(B5)
Matrices D
x
,N
x
and K are the same as in Appendix A.
The symmetry properties described by equations (7),(9)
and (13) are easily conﬁrmed.Finally,note that in the
frequency domain formulation,the temporal convolution
kernels b
p
(x,t)∗ and b
v
(x,t)∗ in matrix B are replaced
by complex frequencydependent functions
ˆ
b
p
(x,ω) and
ˆ
b
v
(x,ω),respectively.
APPENDIX C:MOMENTUM TRANSPORT
The nonlinear equation of motion for a viscous ﬂuid
reads [10,11]
ρ
Dv
i
Dt
−
∂τ
ij
∂x
j
= f
i
,(C1)
where v
i
is the particle velocity,τ
ij
the stress tensor,ρ
the mass density and f
i
the volume density of external
force.The stress tensor is symmetric,i.e.,τ
ij
= τ
ji
.
Stoke’s stressstrainrate relation reads
−τ
ij
+η
ijkl
∂v
k
∂x
l
= pδ
ij
,(C2)
5
where η
ijkl
is the anisotropic viscosity tensor and p the
hydrostatic pressure.The viscosity tensor obeys the fol
lowing symmetry relation
η
ijkl
= η
jikl
= η
ijlk
= η
klij
.(C3)
For isotropic ﬂuids the viscosity tensor reads
η
ijkl
=
¡
ζ −
2
3
η
¢
δ
ij
δ
kl
+η
¡
δ
ik
δ
jl
+δ
il
δ
jk
¢
,(C4)
where η is the dynamic viscosity coeﬃcient and ζ a vol
ume viscosity coeﬃcient (or second viscosity).From en
ergy considerations it follows that η is positive and from
entropy considerations that ζ is positive [12].Eliminating
τ
ij
fromequations (C1) and (C2) yields the NavierStokes
equation,according to
ρ
Dv
i
Dt
−
∂
∂x
i
³
¡
ζ −
2
3
η
¢
∂v
k
∂x
k
´
−
∂
∂x
j
"
η
³
∂v
i
∂x
j
+
∂v
j
∂x
i
´
#
= f
i
−
∂p
∂x
i
,(C5)
or,for constant η and ζ,
ρ
Dv
Dt
−
¡
ζ +
1
3
η
¢
∇(∇¢ v) −η∇
2
v +∇p = f,(C6)
where
v =
v
1
v
2
v
3
and f =
f
1
f
2
f
3
.(C7)
On the other hand,equations (C1) and (C2) can be com
bined to yield the general matrixvector equation (1).To
this end we ﬁrst rewrite these equations as
ρ
Dv
Dt
−D
α
τ
α
= f,(C8)
−τ
α
+h
αβ
D
β
v = pδ
α
,(C9)
(lower case Greek subscripts take on the values 1 and 2),
where
τ
1
=
τ
11
τ
22
τ
33
,τ
2
=
τ
23
τ
31
τ
12
,(C10)
δ
1
=
1
1
1
,δ
2
=
0
0
0
,(C11)
h
11
=
η
1111
η
1122
η
1133
η
1122
η
2222
η
2233
η
1133
η
2233
η
3333
,h
12
=
η
1123
η
1131
η
1112
η
2223
η
2231
η
2212
η
3323
η
3331
η
3312
(C12)
h
21
= h
T
12
,h
22
=
η
2323
η
2331
η
2312
η
2331
η
3131
η
3112
η
2312
η
3112
η
1212
,(C13)
and
D
1
=
∂
∂x
1
0 0
0
∂
∂x
2
0
0 0
∂
∂x
3
,D
2
=
0
∂
∂x
3
∂
∂x
2
∂
∂x
3
0
∂
∂x
1
∂
∂x
2
∂
∂x
1
0
.
(C14)
Hence,we obtain
¯
A
Du
Dt
+
¯
Bu +CD
x
u =
¯
s,(C15)
with
u =
v
−τ
1
−τ
2
,¯s =
f
pδ
1
0
,(C16)
¯
A=
ρI O O
O O O
O O O
,
¯
B =
O O O
O I O
O O I
,(C17)
C =
I O O
O h
11
h
12
O h
21
h
22
,D
x
=
O D
1
D
2
D
1
O O
D
2
O O
,(C18)
I being the 3 × 3 identity matrix and O the 3 × 3 null
matrix.For the isotropic situation we have
h
11
=
ζ +
4
3
η ζ −
2
3
η ζ −
2
3
η
ζ −
2
3
η ζ +
4
3
η ζ −
2
3
η
ζ −
2
3
η ζ −
2
3
η ζ +
4
3
η
,h
22
=
η 0 0
0 η 0
0 0 η
,
(C19)
and h
12
= h
21
= O.The determinant of Cequals 12ζη
5
.
Since η > 0 and ζ > 0 the inverse of C exists for the
isotropic situation.Whether the inverse of C exists in
general for anisotropic ﬂuids remains to be investigated.
Multiplication of all terms in equation (C15) by the in
verse of C and linearization of the term Du/Dt yields
A
∂u
∂t
+Bu +D
x
u = s,(C20)
with
A= C
−1
¯
A=
¯
A,B = C
−1
¯
B and s = C
−1
¯s.
(C21)
Matrices N
x
and K,appearing in the modiﬁed diver
gence theorems (6) and (8),read
N
x
=
O N
1
N
2
N
1
O O
N
2
O O
,(C22)
N
1
=
n
1
0 0
0 n
2
0
0 0 n
3
,N
2
=
0 n
3
n
2
n
3
0 n
1
n
2
n
1
0
(C23)
6
and
K= diag(1,−1,−1),with1 = (1,1,1).(C24)
Based on the structure of the matrices
¯
A,
¯
B,C,D
x
and
K we ﬁnd that the symmetry properties described by
equations (7),(9) and (13) are obeyed.
APPENDIX D:ELASTODYNAMIC WAVE
PROPAGATION IN A SOLID
The linearized equation of motion in solids reads [13–
15]
ρ
∂v
i
∂t
−
∂τ
ij
∂x
j
= f
i
,(D1)
where v
i
is the particle velocity associated to the elasto
dynamic wave motion,τ
ij
the stress tensor,ρ the mass
density of the medium in equilibrium and f
i
the volume
density of external force.The stress tensor is symmetric,
i.e.,τ
ij
= τ
ji
.
Hooke’s linearized stressstrain relation reads
−
∂τ
ij
∂t
+c
ijkl
∂v
k
∂x
l
= c
ijkl
h
kl
,(D2)
where h
kl
is the external deformation rate,with h
kl
=
h
lk
,and c
ijkl
is the anisotropic stiﬀness tensor.The stiﬀ
ness tensor obeys the following symmetry relation [16]
c
ijkl
= c
jikl
= c
ijlk
= c
klij
.(D3)
For isotropic solids this tensor reads
c
ijkl
= λδ
ij
δ
kl
+¹(δ
ik
δ
jl
+δ
il
δ
jk
),(D4)
where λ and ¹ are the Lam´e parameters.Eliminating τ
ij
from equations (D1) and (D2),using equation (D4) and
taking h
kl
= 0,yields the elastodynamic wave equation,
according to
ρ
∂
2
v
i
∂t
2
−
∂
∂x
i
³
λ
∂v
k
∂x
k
´
−
∂
∂x
j
"
¹
³
∂v
i
∂x
j
+
∂v
j
∂x
i
´
#
=
∂f
i
∂t
,
(D5)
or,for constant λ and ¹,
ρ
∂
2
v
∂t
2
−(λ +2¹)∇(∇¢ v) +¹∇×∇×v =
∂f
∂t
,(D6)
with
v =
v
1
v
2
v
3
and f =
f
1
f
2
f
3
.(D7)
On the other hand,equations (D1) and (D2) can be com
bined to yield the general matrixvector equation (1).
To this end we ﬁrst rewrite these equations (for the
anisotropic situation) as
ρ
∂v
∂t
−D
α
τ
α
= f (D8)
and
−
∂τ
α
∂t
+c
αβ
D
β
v = c
αβ
h
β
(D9)
(lower case Greek subscripts take on the values 1 and 2),
where
τ
1
=
τ
11
τ
22
τ
33
,τ
2
=
τ
23
τ
31
τ
12
,(D10)
h
1
=
h
11
h
22
h
33
,h
2
=
2h
23
2h
31
2h
12
,(D11)
c
11
=
c
1111
c
1122
c
1133
c
1122
c
2222
c
2233
c
1133
c
2233
c
3333
,c
12
=
c
1123
c
1131
c
1112
c
2223
c
2231
c
2212
c
3323
c
3331
c
3312
(D12)
c
21
= c
T
12
,c
22
=
c
2323
c
2331
c
2312
c
2331
c
3131
c
3112
c
2312
c
3112
c
1212
,(D13)
and
D
1
=
∂
∂x
1
0 0
0
∂
∂x
2
0
0 0
∂
∂x
3
,D
2
=
0
∂
∂x
3
∂
∂x
2
∂
∂x
3
0
∂
∂x
1
∂
∂x
2
∂
∂x
1
0
.
(D14)
Hence,we obtain
¯
A
∂u
∂t
+CD
x
u =
¯
s,(D15)
with
u =
v
−τ
1
−τ
2
,¯s =
f
c
1β
h
β
c
2β
h
β
,
¯
A=
ρI O O
O I O
O O I
,
(D16)
C =
I O O
O c
11
c
12
O c
21
c
22
,D
x
=
O D
1
D
2
D
1
O O
D
2
O O
,(D17)
I being the 3 × 3 identity matrix and O the 3 × 3 null
matrix.Note that C is a symmetric realvalued matrix.
From energy considerations it follows that it is positive
deﬁnite [15],hence its inverse exists.Multiplying all
terms in equation (D15) by the inverse of C yields
A
∂u
∂t
+D
x
u = s,(D18)
7
with
A= C
−1
¯
A and s = C
−1
¯s,(D19)
or
A=
ρI O O
O s
11
2s
12
O 2s
21
4s
22
,s =
f
h
1
h
2
,(D20)
with
s
11
=
s
1111
s
1122
s
1133
s
1122
s
2222
s
2233
s
1133
s
2233
s
3333
,s
12
=
s
1123
s
1131
s
1112
s
2223
s
2231
s
2212
s
3323
s
3331
s
3312
(D21)
s
21
= s
T
12
,s
22
=
s
2323
s
2331
s
2312
s
2331
s
3131
s
3112
s
2312
s
3112
s
1212
,(D22)
where the s
ijkl
are the elements of the compliance tensor,
with s
ijkl
= s
jikl
= s
ijlk
= s
klij
.Note that
c
ijkl
s
klmn
= s
ijkl
c
klmn
=
1
2
(δ
im
δ
jn
+δ
in
δ
jm
).(D23)
For an isotropic solid we have
c
11
=
λ +2¹ λ λ
λ λ +2¹ λ
λ λ λ +2¹
,c
22
=
¹ 0 0
0 ¹ 0
0 0 ¹
,
(D24)
c
21
= c
T
12
= O,(D25)
and
s
11
=
Λ+2M Λ Λ
Λ Λ+2M Λ
Λ Λ Λ+2M
,s
22
=
M 0 0
0 M 0
0 0 M
,
(D26)
s
21
= s
T
12
= O,(D27)
with
Λ = −
λ
(3λ +2¹)2¹
,M =
1
4¹
.(D28)
Matrices N
x
and K,appearing in the modiﬁed diver
gence theorems (6) and (8),read
N
x
=
O N
1
N
2
N
1
O O
N
2
O O
,(D29)
N
1
=
n
1
0 0
0 n
2
0
0 0 n
3
,N
2
=
0 n
3
n
2
n
3
0 n
1
n
2
n
1
0
(D30)
and
K= diag(1,−1,−1),with1 = (1,1,1).(D31)
Based on the structure of the matrices
¯
A,C,D
x
and
K we ﬁnd that the symmetry properties described by
equations (7),(9) and (13) are obeyed.
APPENDIX E:ELECTROMAGNETIC
DIFFUSION AND WAVE PROPAGATION
The Maxwell equations for electromagnetic phenomena
read [15,17]
ǫ
∂E
∂t
+σE−∇×H = −J
e
,(E1)
¹
∂H
∂t
+∇×E = −J
m
,(E2)
where E and H are the electric and magnetic ﬁeld
strengths,ǫ,¹ and σ are the permittivity,permeabil
ity and conductivity,respectively and,ﬁnally,J
e
and J
m
are source functions in terms of the external electric and
magnetic current densities.Eliminating the magnetic
ﬁeld strength Hfrom equations (E1) and (E2) yields the
electromagnetic wave equation,according toσ
∂E
∂t
+ǫ
∂
2
E
∂t
2
+∇×
³
1
¹
∇×E
´
= −
∂J
e
∂t
−∇×
³
1
¹
J
m
´
.
(E3)
Note that ﬁrst and second order time derivatives occur in
this equation,which means that this equation accounts
for diﬀusion as well as wave propagation.On the other
hand,equations (E1) and (E2) can be combined to
A
∂u
∂t
+Bu +D
x
u = s,(E4)
with
u =
µ
E
H
¶
,s =
µ
−J
e
−J
m
¶
,A=
µ
ǫI O
O ¹I
¶
,B =
µ
σI O
O O
¶
,
(E5)
D
x
=
µ
O D
T
0
D
0
O
¶
,D
0
=
0 −
∂
∂x
3
∂
∂x
2
∂
∂x
3
0 −
∂
∂x
1
−
∂
∂x
2
∂
∂x
1
0
,(E6)
with I being the 3×3 identity matrix and Othe 3×3 null
matrix.Matrices N
x
and K,appearing in the modiﬁed
divergence theorems (6) and (8),read
N
x
=
µ
O N
T
0
N
0
O
¶
,N
0
=
0 −n
3
n
2
n
3
0 −n
1
−n
2
n
1
0
,(E7)
K= diag(−1,1),1 = (1,1,1).(E8)
The symmetry properties described by equations (7),(9)
and (13) are easily conﬁrmed.
APPENDIX F:COUPLED ELASTODYNAMIC
AND ELECTROMAGNETIC WAVE
PROPAGATION IN POROUS SOLIDS
We brieﬂy review the theory for elastodynamic waves
coupled to electromagnetic ﬁelds in a dissipative inhomo
geneous anisotropic ﬂuidsaturated porous solid [6,18].
8
The linearized equations of motion read in the frequency
domain (using the vector notation introduced in Ap
pendix D)
jωρ
b
ˆv
s
+jωρ
f
ˆw−D
α
ˆτ
bα
=
ˆ
f
b
,(F1)
jωρ
f
ˆv
s
+η
ˆ
k
−1
³
ˆw−
ˆ
L
ˆ
E
´
+∇ˆp =
ˆ
f
f
,(F2)
with ˆw = φ(ˆv
f
− ˆv
s
).Here ˆv
s
and ˆv
f
are the averaged
solid and ﬂuid particle velocities associated to the wave
motion,ˆwis the ﬁltration velocity,φ the porosity,ˆτ
bα
the
averaged bulk stress (organized as in equation (D10)),ˆp
the averaged ﬂuid pressure and
ˆ
E the averaged electric
ﬁeld strength.Matrix D
α
contains the spatial diﬀerential
operators and is deﬁned in equation (D14).The source
functions
ˆ
f
b
and
ˆ
f
f
are the volume densities of exter
nal force on the bulk and on the ﬂuid,respectively.The
constitutive parameters ρ
b
and ρ
f
are the anisotropic
bulk and ﬂuid mass densities,respectively [19].In the
following we assume that these tensors are symmetric,
according to ρ
b
= (ρ
b
)
T
and ρ
f
= (ρ
f
)
T
,which is for
example the case when the anisotropy is the result of
parallel ﬁne layering at a scale much smaller than the
wavelength.The complex frequencydependent tensor
ˆ
k
is the dynamic permeability tensor of the porous mate
rial,with
ˆ
k =
ˆ
k
T
,and η is the ﬂuid viscosity param
eter.Finally,the complex frequencydependent tensor
ˆ
L accounts for the coupling between the elastodynamic
and electromagnetic waves.In the following we will as
sume that this tensor is symmetric as well,according to
ˆ
L =
ˆ
L
T
(Pride and Haartsen [6] discuss the conditions
for this symmetry).
The linearized stressstrain relations read
−jωˆτ
bα
+c
αβ
D
β
ˆv
s
+d
α
∇¢ ˆw = 0,(F3)
jωˆp +d
Tα
D
α
ˆv
s
+M∇¢ ˆw = 0,(F4)
with
d
1
=
d
11
d
22
d
33
,d
2
=
d
23
d
31
d
12
,(F5)
c
αβ
deﬁned in equations (D12) and (D13) and 0 a 3 ×1
null vector.M,d
ij
and c
ijkl
are the realvalued stiﬀness
parameters of the ﬂuidsaturated porous solid.
Maxwell’s electromagnetic ﬁeld equations read
jωǫ
ˆ
E+
ˆ
J −∇×
ˆ
H = −
ˆ
J
e
,(F6)
jω¹
ˆ
H+∇×
ˆ
E = −
ˆ
J
m
,(F7)
where
ˆ
H is the averaged magnetic ﬁeld strength,
ˆ
J the
averaged induced electric current density,ǫ and ¹ are the
anisotropic permittivity and permeability,with ǫ = ǫ
T
and ¹ = ¹
T
,and J
e
and J
m
are source functions in terms
of the external electric and magnetic current densities.
The induced electric current density is coupled to the
elastodynamic wave motion,according to
ˆ
J = ˆσ
ˆ
E−
ˆ
L
h
∇ˆp +jωρ
f
ˆv
s
−
ˆ
f
f
i
,(F8)
where ˆσ is the complex frequency dependent conductiv
ity,with ˆσ = ˆσ
T
.Substituting the constitutive relation
(F8) into the Maxwell equation (F6),and adding
ˆ
L times
equation (F2) to equation (F6) in order to compensate
for the term −
ˆ
L[∇ˆp +jωρ
f
ˆv
s
−
ˆ
f
f
],yields
jωǫ
ˆ
E+
³
ˆσ−η
ˆ
L
ˆ
k
−1
ˆ
L
´
ˆ
E+η
ˆ
L
ˆ
k
−1
ˆw−∇×
ˆ
H= −
ˆ
J
e
.(F9)
Equations (F9) and (F7),together with equations (F1),
(F2),(F3) and (F4) can be combined to yield
jω
¯
Aˆu +
¯
Bˆu +CD
x
ˆu =
ˆ¯s,(F10)
where
ˆu =
ˆ
u
1
ˆu
2
ˆu
3
,
ˆ¯s =
ˆ¯s
1
ˆ¯s
2
ˆ¯s
3
,
¯
A=
¯
A
11
O O
O
¯
A
22
¯
A
23
O
¯
A
T
23
¯
A
33
,
¯
B =
¯
B
11
O
¯
B
13
O O O
−
¯
B
T
13
O
¯
B
33
,(F11)
C =
I O O
O C
22
C
23
O C
T
23
C
33
and D
x
=
D
11
O O
O D
22
O
O O D
33
,
(F12)
where I and O are identity and null matrices of appro
priate size and
ˆu
1
=
µ
ˆ
E
ˆ
H
¶
,ˆu
2
=
ˆv
s
−
ˆ
τ
b
1
−ˆτ
b
2
,ˆu
3
=
µ
ˆw
ˆp
¶
,
ˆ¯s
1
=
µ
−
ˆ
J
e
−
ˆ
J
m
¶
,
ˆ
¯s
2
=
ˆ
f
b
0
0
,
ˆ
¯s
3
=
µ
ˆ
f
f
0
¶
,(F13)
¯
A
11
=
µ
ǫ O
O ¹
¶
,
¯
A
22
=
ρ
b
O O
O I O
O O I
,
¯
A
23
=
ρ
f
0
O 0
O 0
,
¯
A
33
=
µ
O 0
0
T
1
¶
,(F14)
¯
B
11
=
µ
(ˆσ −η
ˆ
L
ˆ
k
−1
ˆ
L) O
O O
¶
,
¯
B
13
=
µ
η
ˆ
L
ˆ
k
−1
0
O 0
¶
,
¯
B
33
=
µ
η
ˆ
k
−1
0
0
T
0
¶
,(F15)
9
C
22
=
I O O
O c
11
c
12
O c
21
c
22
,C
23
=
O 0
O d
1
O d
2
,
C
33
=
µ
I 0
0
T
M
¶
,(F16)
D
11
=
µ
O D
T
0
D
0
O
¶
,D
0
=
0 −
∂
∂x
3
∂
∂x
2
∂
∂x
3
0 −
∂
∂x
1
−
∂
∂x
2
∂
∂x
1
0
,
D
33
=
µ
O ∇
∇
T
0
¶
(F17)
and D
22
equal to D
x
in equation (D17).Note that C in
equation (F12) is a symmetric realvalued matrix.From
energy considerations it follows that it is positive deﬁnite
[6,15],hence its inverse exists.Multiplying all terms in
equation (F10) by the inverse of C ﬁnally yields
jωAˆu +Bˆu +D
x
ˆu =ˆs,(F18)
with A= C
−1
¯
A,B = C
−1
¯
B =
¯
B and ˆs = C
−1
ˆ¯s =
ˆ
¯s.Matrices N
x
and K,appearing in the modiﬁed diver
gence theorems (6) and (8),read
N
x
=
N
11
O O
O N
22
O
O O N
33
,N
11
=
µ
O N
T
0
N
0
O
¶
,
N
0
=
0 −n
3
n
2
n
3
0 −n
1
−n
2
n
1
0
,N
33
=
µ
O n
n
T
0
¶
,(F19)
K= diag(−1,1,1,−1,−1,1,−1),(F20)
and N
22
equal to N
x
in equation (D29).Based on the
structure of the matrices
¯
A,
¯
B,C,D
x
and K as well as
the symmetry relations discussed above,we ﬁnd that the
symmetry properties described by equations (7),(9) and
(13) are obeyed.
Finally,note that when the coupling tensor
ˆ
L is zero,
the matrix
¯
B
13
vanishes and hence equation (F10) de
couples into the electromagnetic wave equation for the
wave vector ˆu
1
and Biot’s poroelastic wave equation for
the wave vector (ˆu
T
2
,ˆu
T
3
)
T
,[20].For a nonporous solid
the matrices
¯
A
23
and C
23
vanish as well,so Biot’s wave
equation reduces to the elastodynamic wave equation for
the wave vector
ˆ
u
2
.Obviously the symmetry proper
ties described by equations (7),(9) and (13) are obeyed
for the matrices appearing in the electromagnetic wave
equation,Biot’s poroelastic wave equation and the elas
todynamic wave equation,respectively.
[1] A.T.de Hoop and H.J.Stam,Wave Motion 10,479
(1988).
[2] N.N.Bojarski,J.Acoust.Soc.Am.74,281 (1983).
[3] J.T.Fokkema and P.M.van den Berg,Seismic applica
tions of acoustic reciprocity (Elsevier,Amsterdam,1993).
[4] L.M.Lyamshev,Doklady Akademii Nauk 138,575
(1961).
[5] J.F.Allard,B.Brouard,and D.Lafarge,Wave Motion
17,329 (1993).
[6] S.R.Pride and M.W.Haartsen,Journal of the Acous
tical Society of America 100,1301 (1996).
[7] B.P.Belinskiy,On some general mathematical proper
ties of the system:Elastic plateacoustic medium (World
Scientiﬁc Publishing Co,2001),pp.193–218.
[8] M.Fisher and K.J.Langenberg,IEEE Trans.Ant.Prop.
32,1080 (1984).
[9] O.A.Godin,Physical Review Letters 97,054301 (2006).
[10] K.K.Kuo,Principles of combustion (John Wiley &Sons,
New York,1986).
[11] J.R.Welty,C.E.Wicks,and R.E.Wilson,Fundamen
tals of momentum,heat,and mass transfer (John Wiley
& Sons,New York,1976).
[12] L.D.Landau and E.M.Lifshitz,Fluid Mechanics (Perg
amon Press,1987).
[13] J.D.Achenbach,Wave propagation in elastic solids
(North Holland,1973).
[14] K.Aki and P.G.Richards,Quantitative seismology,Vol.
I (W.H.Freeman and Company,San Fransisco,1980).
[15] A.T.de Hoop,Handbook of radiation and scattering of
waves (Academic Press,London,1995).
[16] F.A.Dahlen and J.Tromp,Theoretical global seismology
(Princeton University Press,1998).
[17] L.D.Landau and E.M.Lifshitz,Electrodynamics of con
tinuous media (Pergamon Press,1960).
[18] S.Pride,Phys.Rev.B 50,15678 (1994).
[19] M.Schoenberg and P.N.Sen,Journal of the Acoustical
Society of America 73,61 (1983).
[20] M.A.Biot,Journal of the Acoustical Society of America
28,168 (1956).
Enter the password to open this PDF file:
File name:

File size:

Title:

Author:

Subject:

Keywords:

Creation Date:

Modification Date:

Creator:

PDF Producer:

PDF Version:

Page Count:

Preparing document for printing…
0%
Comments 0
Log in to post a comment