Soil Dynamics and Earthquake Engineering 115 (2018) 663–672

Contents lists available at ScienceDirect

Soil Dynamics and Earthquake Engineering

journal homepage: www.elsevier.com/locate/soildyn

Generalized stiffness matrix method for layered soils T

Eduardo Kausel
Massachusetts Institute of Technology, Cambridge, MA 02139, USA

ARTICLE INFO ABSTRACT

Keywords: This article provides an effective extension to the well-known Stiffness Matrix Method (SMM) for layered soils.
Green's functions The modified algorithm allows considering elastic media subjected to both 2.5-D and 3-D sources that are neither
Fundamental solutions plane nor axisymmetric. These sources can act on –or within– viscoelastic, horizontally layered media of either
Layered media finite or infinite depth. Although a direct approach in the context of a fully 3-D formulation in Cartesian co-
Wave propagation
ordinates in combination with a 2-D spatial Fourier transform is possible and could be used for the solution of the
Seismology
Elastodynamics
problem at hand, it will be shown that a substantial saving in computational effort can be achieved by for-
Soil dynamics mulating the evaluation of the Green's functions with the same cost as that expended for plane-strain or for
Seismic sources axisymmetric sources, provided the material properties are either isotropic or transversely isotropic. We de-
Moving loads monstrate herein that the flexibility functions in the frequency-wavenumber domain for the general 3-D problem
2.5-D problem can be inferred from the exact rigidity matrices for the 2-D (i.e. plane-strain) problems of SVP and SH waves, and
Thompson-Haskell problem only these need to be inverted via Gaussian reduction to obtain the requisite flexibility functions in the fre-
quency-wavenumber domain. This strategy relates to what in the technical literature is sometimes referred to as
the “inversion of the descent of dimensions”, a principle that states that the response for point loads can be inferred
from the response to line loads. It allows obtaining highly effective numerical solutions to the problem at hand.

1. Introduction interaction, and in the evaluation of pavements under impact loads,

among many other fields of application. A detailed, modern rendition of
In general, closed-form solutions for layered strata or layered half- the method is given in Kausel [3], which includes also the formulations
spaces subjected to dynamic sources do not exist. Indeed, even the for layered cylinders as well as layered spheres that have not been
homogeneous free-free plate—the so-called Mindlin plate—is ulti- published anywhere else. In the case of horizontally layered media, the
mately intractable by purely analytical means. Thus, other than formulation is available for both two-dimensional, plane strain sources
homogeneous elastic half-spaces for which a limited number of closed- (i.e. SVP and SH sources), and for three-dimensional sources expressible
form solutions exist, such problems must be solved with the aid of in cylindrical coordinates, such as disk loads, point loads, or double
numerical tools, which are needed to evaluate the Fourier or Hankel couples. In the ensuing, we extend the formulation to truly three-di-
transforms from the wavenumber domain back into the spatial domain. mensional sources with arbitrary, non-symmetric spatial distributions,
Among the available formulations, a very widely used and well-known and demonstrate its application by means of examples. For reasons of
scheme is the stiffness or impedance matrix method (SMM) of Kausel space, the presentation will be rather terse, and we shall assume the
and Röesset [1], which is based on a rigorous formulation for waves in reader to be already familiar with the original stiffness matrix method.
layered media. Since such matrix formulation is exact and does not The reader should notice that although the ensuing formulation can
entail any discretization of the medium, it follows that the layers can be readily be generalized to transversely isotropic media, that is not the
arbitrarily thick and the frequencies as high as necessary to capture the case for orthotropic and anisotropic media, for those material config-
temporal characteristics of the source. Thus, the exact rigidity matrices urations would require an altogether different formulation because of
for the continuum considered herein should not be confused with the the material dependence on the azimuth (i.e. the horizontal angle in
discrete dynamic stiffness matrices obtained with finite elements or cylindrical coordinates). In that case, a direct formulation in Cartesian
even with those in the semi-discrete Thin-Layer Method [2]. coordinates would be required. Indeed, one such general formulation
Since its inception in 1981, the stiffness matrix method (SMM) has has already been given by Rokhlin and Wang [4]. Also, an extension of
found wide use in elastodynamics, especially in geophysics and seis- the SMM to isotropic poroelastic media can be found in a 1995 paper by
mology, in non-destructive testing, in soil dynamics and soil-structure Rajapakse and Senjuntichai [5]. We mention in passing that the

E-mail address: [email protected].

https://doi.org/10.1016/j.soildyn.2018.09.003
Received 15 July 2018; Received in revised form 6 September 2018; Accepted 6 September 2018
Available online 11 October 2018
0267-7261/ © 2018 Elsevier Ltd. All rights reserved.
E. Kausel Soil Dynamics and Earthquake Engineering 115 (2018) 663–672

computational strategy to be described could just as well be ported to where is the azimuthal angle, i.e. it defines the direction of wave
other methods of solving waves in layered strata, such as the trans- propagation. Then
mission-reflection method, although that would probably require some n1 = ±kp, n2 = ± ks , n3 = ±ks (8)
considerable effort because of the symmetries thereby lost.
cos sin ± i s cos
= exp( kpz ), = exp( ksz ), = ± i s sin exp( ksz )
2. Layer stiffness matrices in three-dimensional space 1 sin 2 cos 3
ip 0 1

Consider a horizontally layered, viscoelastic medium subjected to (9)

dynamic sources acting at some arbitrary location. The system com- The proof that (6)–(9) is indeed the solution of the eigenvalue
prises a set of homogenous, isotropic layers of arbitrary finite thickness, problem can be found by substitution of these expressions into (5). We
which are (optionally) underlain by an infinitely deep, homogeneous define next the modified displacement vector ũ(k x , k y, z, ) in the fre-
half-space. In the absence of body loads, the wave equation within each quency-wavenumber space, to which we add a convenience factor
layer in 3-D space is of the form i= 1 so as to arrive at stiffness matrices which are symmetric.
From the preceding, it follows that
+ 2µ 0 0 2u
µ 0 0 2u
ux
u= 0 µ 0 + 0 + 2µ 0 cos
uy = sin
sin s cos cos sin s cos
x2 y2 ˜ =
u cos s sin sin cos s sin Ea exp i( t kx x ky y )
0 0 µ 0 0 µ i uz p 0 1 p 0 1

µ 0 0 2u 0 +µ 0 2u (10)
+ 0 µ 0 + +µ 0 0
z2 x y where
0 0 + 2µ 0 0 0
E = diag [exp( kpz ), exp( ksz ), exp( ksz ), exp(kpz ), exp(ksz ), exp(ksz )]
0 0 0 2u 0 0 +µ 2u
+ 0 0 +µ + 0 0 0 (11)
y z x z
0 +µ 0 +µ 0 0 (1) and a is a vector of arbitrary constants. In turn, this can be written
compactly as
where (z ), µ (z ), (z ) are, respectively, the Lamé constant, the shear
modulus and the mass density (which are constant within each layer), ux
u y = R { R1E 1 R2Ez } a1 exp i( t kx x ky y)
x , y are the two horizontal Cartesian coordinates, and z is the vertical z a2
coordinate, which we shall assume to be positive in the upward direc- i uz (12)
tion. Later on we shall also need the P and S wave velocities where the a1, a2 are again vectors of arbitrary constants, E = { E z 1
Ez } ,
= ( + 2µ )/ and = µ/ . Also, the displacement vector u is of and
the form
EZ = diag [exp(kpz ) exp(ksz ) exp(ksz )] (13)
u(x , y, z, t ) = {ux u y uz }T (2)
cos sin 0
where the superscript T denotes the transposed vector. Assume next a R= sin cos 0 = rotation matrix
plane wave solution of the form 0 0 1 (14)
u(x , y, z, t ) = exp i( t kx x k y y ) exp( nz ) (3)
1 0 s 1 0 s
where i = 1 , is the frequency, k x , k y are the two horizontal wa- R1 = 0 1 0 , R2 = 0 1 0
venumbers, n = i k z is the vertical wavenumber, and is a vector of p 0 1 p 0 1 (15)
constants. It follows that By taking appropriate spatial derivatives, we can also compute the
u u u u 2u stresses in horizontal planes associated with the above solution. After
= ik x u , = ik y u , = nu, =i u 2 = k x2 u , some relatively tedious algebra, we find
x y z t x
2u 2u 2u
xz
= k y2 u , = n2 u , = 2u
a1
y2 z2 t2 s˜ = yz = kµ R { Q1 E z 1 Q2 E z } a exp i( t kx x ky y)
2
(4) i z (16)
and by substitution into the wave Eq. (1), we obtain with

2 0 0 2 0 0 2 0 0 0 1 0
2I 0 2 0 k x2 2
0 2 0 ky + 0 2 0 n2 ( 2 2)
1 0 0 kx k y
0 0 2 0 0 2
0 0 2 0 0 0

0 0 1 0 0 0
+( 2 2) 2 2)
0 0 0 i kx n + ( 0 0 1 i ky n =0
1 0 0 0 1 0 (5)

For fixed values of the horizontal wavenumbers {C}k x , k y and fre-

2p 0 (1 + s 2 ) 2p 0 1 + s2
quency , the above expression constitutes a 3 × 3 eigenvalue problem
Q1 = 0 s 0 , Q2 = 0 s 0
in the vertical wavenumber, which has three distinct eigenvalues nj and
(1 + s 2) 0 2s 1 + s2 0 2s (17)
eigenvectors j for coupled P, SH, and SV waves, respectively. They are
as follows. Define: By considering an isolated layer as a free body in space subjected to
k= k x2 + k y2 , k x = k cos , k y = k sin external tractions at its upper and lower bounding planes, equilibrating
(6)
these tractions with the internal stresses at those elevations given by Eq.
2 2 (16), evaluating the displacements at those same elevations via Eq. (12)
kp = k2 , ks = k2 and eliminating the integration constants a1, a2 between these two sets
(7) of equations, one arrives at a matrix expression relating the external

664
E. Kausel Soil Dynamics and Earthquake Engineering 115 (2018) 663–672

tractions to the displacements at the two bounding horizons. The the elements of the flexibility submatrices depend on the pair of re-
symmetric matrix relating the external tractions to the displacement is ceiver-source indices m , n , i.e. fij = fijmn (k , ) , but do not depend ex-
the stiffness (or rigidity, or impedance) matrix that we wished to ob- plicitly on the Cartesian wavenumbers k x , k y . In general f13 f31, unless
tain. We can then proceed to assemble the stiffness matrix for the m = n , in which case they are indeed equal.
complete layered medium by overlapping appropriately the matrices Please notice that the flexibility elements Fmn can be obtained in-
for the individual layers. dependently of the relative values of k x , k y , i.e. they are identical to the
Examination of Eqs. (15, 17) reveals that R1, R2, Q1, Q2 are identical flexibility elements for the uncoupled plane-strain (2-D) problems in-
to the matrices used in the plane-strain and cylindrical problems, e.g. volving SVP and SH waves in k, . Hence, the response functions in the
Chapter 10 in [3]. Thus, we can already discern that the layer stiffness frequency-wavenumber domain can be obtained directly from the in-
matrices for the fully 3-D case can be written as version of the 2-D stiffness matrices for any layered system, which
1 brings about considerable savings in the computational effort.
R O K11 K12 R O 1 With the shorthand
K3D = , R = RT
O R K21 K22 O R
2D kx ky
C = cos = , S = sin = (23)
cos sin 0 k k
= sin cos 0 we can write Eq. (21) as
0 0 1 (18)
ux f11 0 f13 px
where the symmetric stiffness matrix at the middle is the same as the C S 0 C S 0
uy = S C 0 0 f22 0 S C 0 py
one used for the 2-D plane-strain and 3-D cylindrical problems. Hence,
i uz 0 0 1 f31 0 f33 0 0 1 i pz
the system matrix is exactly as for these two cases, except that each m mn n

triplet of rows and columns must be multiplied by R and RT , respec- f11 C 2 + f22 S 2 (f11 f22 ) SC f13 C px
tively. Hence, the global equation for the system of layers is of the form = (f11 f22 ) SC f11 S 2 + f22 C 2 f13 S py
p1 K11 K12 O O f31 C f32 S f31 S + f32 C f33 i pz
R O O O RTO O O u1 mn n (24)
p2 O R K21 K22 K23 O RT u2
=
O O K32 KN 1, N O Hence,
pN O O R O KN , N KNN O O RT uN
1 2D
ux f11 cos2 + f22 sin2 (f11 f22 )sin cos f13 cos px
(19) uy = (f11 f22 )sin cos f11 sin2 + f22 cos2 f13 sin py
i uz i pz
whose inverse is m f31 cos f31 sin f33
mn n

1
RTp1 (24′)
RT u1 K11 K12 O O
RT u2 K21 K22 K23 RTp2 Observe that by setting = 0 , then cos = 1, sin = 0 and we re-
= cover the classical plane-strain solution. At this point in time we now
O K32 KN 1, N
RTun O KN , N 1 KNN RT pN proceed to undo the convenience imaginary factor in which case Eq.
2D
(24) changes into
F11 F12 F1N RT p1
ux f11 cos2 + f22 sin2 (f11 f22 )sin cos i f13 cos px
F21 F22 RT p2 um = uy py
= = (f11 f22 )sin cos f11 sin2 + f22 cos2 i f13 sin = Gmn pn
FN 1, N uz m
i f31 cos i f31 sin f33 pz
n
mn
FN 1 FN , N FNN RTpN
1 2D (20) (25)
where the elements Kij of the stiffness matrix agree perfectly with those which provides the formal solution for the response functions to the
used in the two plane strain cases. general three-dimensional source problem. Since the coefficients f13 , f31
We observe that the dependence on k x , k y of the solution for the are usually distinct, this implies that the vertical displacement at the
layered system in the frequency-wavenumber domain is separable, interface m due to a horizontal source at interface n does not equal the
since the flexibility matrix is solely a function of the radial wavenumber horizontal displacement at interface m due to a vertical source at in-
k = k x2 + k y2 and of the frequency , while the explicit dependence on terface n . In other words, it is not enough to change source and receiver
the two horizontal wavenumbers comes only through the rotational directions, but also their locations for reciprocity to be satisfied.
matrices RT , which are common to all layers. This is also the reason why The displacements in the space domain can now be obtained —at
orthotropic and anisotropic media —except for transverse isotropy— least formally— from the two-dimensional, inverse Fourier transform
had to be disallowed.
( ) 1 2 + +
um ( x , y , ) = um (k x , k y, )e i kx x e i ky y dk
The displacements in the frequency-wavenumber domain observed 2 x dk y

at the m-th elevation caused by sources placed at the n-th elevation are
=( )
1 2 + +
Gmn pn e i kx x e i ky ydk
x dk y (26)
then 2

um (k x , k y, ) = R Fmn (k, ) RT pn (k x , k y, ) (21) in which the flexibility elements fij = fij (k , ) = fij ( + ) de- k x2 k y2 ,
pend only indirectly on both k x and k y , i.e. on the radial wavenumber
where each of the 3 × 3 flexibility matrices Fmn = Fmn (k, ) (= sub- k= k x2 + k y2 .
matrices of the inverse of the block-tridiagonal, symmetric dynamic
stiffness matrix in Eq. (20)) are of the form 3. Summary of procedure
f11 0 f13
In a nutshell, the flexibility or Green's functions Gmn are obtained as
Fmn = 0 f22 0
follows:
f31 0 f33 (22)
mn

in which the elements f12 = f21 = f23 = f32 = 0 are null because the in- • Find the Fourier transform of the sources, p (k , k ,
x x y ) , py (k x , k y , ) ,
plane (SV-P) degrees of freedom are uncoupled from the out-of-plane pz (k x , k y, ) . In most cases, this can be accomplished in closed form,
(SH) degrees of freedom. This is because the [1,2], [2,1], [2,3], and and is usually a straightforward task.
[3,2] elements in the R1, R2, Q1, Q2 matrices are all zero. Observe that • Assemble the global, symmetric SVP system stiffness matrix for each

665
E. Kausel Soil Dynamics and Earthquake Engineering 115 (2018) 663–672

k, of interest and find f11 , f13 , f31 , f33 at the desired elevations. kPz = ± kP2 k2 kSz = ± kS2 k2 (28)
• Assemble the global, symmetric SH system stiffness matrix for each
We define
k, of interest, and find f22 at the desired elevations.
• Obtain the displacements in the frequency-wavenumber space via
kp = ikPz = k2 kP2 ks = ikSz = k2 kS2 (29)
Eq. (26). The inversion of this equation into the space domain is
considered in the next section.
Hence, the propagation of P or S waves within a specific layer is
given by expressions of the form
4. Fourier inversion into the space-time domain AP exp[i( t kx )]exp( kp z ), AS exp[i( t kx )]exp( ks z ) (30)

In the preceding sections on the Stiffness Matrix Method (SMM) for where AS , AP are arbitrary amplitudes (they do not matter here) and
arbitrary, spatially distributed loads we learned how to calculate the the sign of the exponential depends on whether the wave moves up or
response functions for layered media when harmonic sources are ap- down. If we focus attention on the last factor and evaluate it at z = h
plied somewhere. Ultimately, we saw in Eq. (26) that this task requires (where z is measured from the interface being considered, and h is the
evaluation of improper integrals over the horizontal wavenumbers. In thickness of the layer in question), and ask ourselves what value the
most cases, such integrals are not amenable to closed form or analytical horizontal wave-number k > kS > kP must achieve so that a wave that
integration, so one must necessarily rely on numerical methods. emanates from one interface will have decayed to nearly zero upon
We mention in passing that such improper integrals are fraught with reaching the neighboring interface at a distance h , i.e. for what value of
difficulties that must adequately be dealt with, lest the results be highly k does the exponential term becomes negligibly small. Clearly this
inaccurate or even totally wrong. Although space considerations do not factor is given by
allow us to cover exhaustively all aspects of this problem, we can
exp ( ks z ) exp( ksh) = exp( h k 2 kS2 ) < = 10 m
(31)
nonetheless provide a summary of the important issues and suggest
some remedies. These issues are as follows: (for example, m=2, or m=3 will do). The solution is then

• The integrals are improper, i.e. their upper limit is infinitely large, h k2 kS2 > m ln 10 or k>
m ln 10 2
+ kS2
and so can also be the lower limit. In addition, the kernels of the h (32)
integrals decay relatively slowly with the wavenumber k , especially
For m= 2, 2 ln 10 = 4.61 this would imply for the
3
when the source and the receiver are at the same elevation. Thus, it 4.71 = 2
behooves to decide appropriately on the maximum wavenumber, current layer
say k max , and at the same time provide a strategy to circumvent 3 2
2
mere truncation so as to capture the contribution of the tail to the 2
k> +
integral above the upper limit k max . h
•
(33)
In addition, the layers strongly couple the various interfaces, but at
sufficiently high wavenumbers this coupling becomes negligible. While a similar limit exists for the P-wave term, it need not be
One must then decide on how large the needed truncation wave- considered because kS > kP . The actual maximum wavenumber used in
number k max should be. a computation depends on the shear wave velocities of the various
• The integrands are highly wavy, and exhibit sharp peaks in the vi- layers as well as their thicknesses, but for practical purposes, it is easiest
cinity of the wavenumbers associated with the normal modes at the to choose the maximum for all of the layers on the basis of the
given frequency. One then needs to be able to estimate a-priori how minimum shear wave velocity and the minimum thickness, i.e.
many such peaks can be expected. min( ), min(h ) , say
• The integrals must be discretized by an appropriate choice of the
wavenumber step. This step has to be small enough so as to properly 3 2

( )
2
2
model the sharp peaks in the integrands. min h
+ min
k max = max
2
We discuss briefly each of these issues in turn. min (34)

Observe that the maximum wavenumber is a function of the current

5. Maximum wavenumber for truncation and layer uncoupling frequency. One simple, conservative alternative is to use max to make
the choice of maximum wavenumber in Eq. (34).
We begin by showing that at sufficiently high wavenumber, the
flexibility functions (the integrands) behave quasi-statically (i.e. as if
5.1. Static asymptotic behavior: tail of integrals
= 0 ), and furthermore, that the layer interfaces decouple as if the
interface belonged to two infinitely deep half-spaces (an upper and a
At high wavenumbers not only do the impedance matrices decouple
lower half-space) with material properties equal to those of the two
into block-diagonal form, but they also converge to the static solution,
layers meeting at that interface. Thus, the stiffness (or impedance)
i.e. ks k , kp k . Moreover, at a specific interface, the static SVP and
matrices for each layer attain a block-diagonal structure, i.e. K12 O so
SH impedances of two dissimilar upper (U) and lower (L) half-spaces
that Klayer [K11, O ; O, K22]. The simplest way to demonstrate this is
with shear moduli µU , µL and Poison's ratios U , L are
by considering the dispersion equations for S and P waves for any of the
layers, while omitting for simplicity the layer sub-index. These equa- µL 1 aL2 µU 1 aU2
tions are then K SVP = 2k + , K SH = k (µL + µU )
1+ aL2 aL2 1 1+ aU2 aU2 1
2 2 (35)
2
k 2 + kPz = kP2 2
k 2 + kSz = kS2
(27) 1 2 L 1 2 U
aL2 = , aU2 =
2 2 L 2 2 U
where k is the radial wavenumber and we have added additional sub-
scripts P, S to the vertical wavenumbers so as to distinguish one from The elements of the inverses of these stiffnesses to which the diag-
the other. From here, we obtain onal blocks Gnn tend to are

666
E. Kausel Soil Dynamics and Earthquake Engineering 115 (2018) 663–672

1
f11 f13 wavenumbers in k , and place these functions into a look-up table for
FSVP = KSVP = , FSH = KSH1 = f22 later use. You should observe that determining those Green's func-
f31 f33 (36)
tions will entail solving a system of linear equations in the block-
with tridiagonal, symmetric matrix in Eq. (26) for a total of N times, and
not (2N ) 2 times, as might have been expected on account of the total
1 µL µU
f11 = f33 = + number of points in the two-dimensional wavenumber grid (i.e.
k 1 + aL2 1 + aU2 (37a) discretizations in k x , k y ). This is a very substantial saving in
1 computational cost, inasmuch as the Gaussian reduction for a multi-
f22 = (µ + µU ) layered soil can be computationally expensive, especially in the
k L (37b)
presence of many layers —as may be needed to simulate vertically
1 aL2 µL aU2 µU inhomogeneous soils whose properties vary gradually with depth.
•
f31 = f13 =
k 1+ aL2 1 + aU2 (37c) Retrieve from the lookup table appropriate functional values and
interpolate these appropriately to obtain the Green's functions at
µL µU
2
aL2 µL aU2 µU
2
any specific pair of horizontal wavenumbers k x , k y defined over a
=2 + uniform rectangular grid of points spaced at equal intervals
1+ aL2 1+ aU2 1+ aL2 1+ aU2 (37d)
k x , k y . The need for interpolation arises from the fact that
It follows that one can readily construct a “static” flexibility matrix k x2 + k y2 does not generally coincide with any point in the radial
(0)
Gnn = Gnn (k 0, = 0) with the elements in Eq. (36, 37), and that for wavenumber set k . Use the interpolated values in the rectangular
sufficiently large k such matrix will have a simple structure where the grid to evaluate the terms of Gmn in Eq. (38).
layers are fully decoupled. Also, Gm n (k 0, = 0) = O. Armed with
this knowledge, we now proceed to replace the integral in Eq. (26) with
• If m = n , subtract the contribution of the “static” terms defined by
Eq. (36). For this purpose, use the static flexibilities FSVP , FSH given
the following: previously to form Gnn (0)
, then subtract these from the diagonal blocks
of the flexibility matrices Gnn , using the appropriate material prop-
( ) 1 2 +kmax +kmax
(0) (0) i kx x e i ky y dk
um = um um = (Gmn Gmn ) pn e x dk y
2 kmax kmax erties for each pair of layers (which play here the roles of “upper”
(38) and “lower” half-spaces). Do this for each and every wavenumber in
The integrand in Eq. (38) is now virtually free from any tail, and so the grid. The net effect is that the tails will have vanished.
the truncation at k max has no deleterious effect on the results. Moreover, • Evaluate numerically the wavenumber integrals up to k max as spe-
the subtraction affects only the diagonal blocks of Gnn , that is, the lo- cified by Eq. (38). This yields um .
cations where the source and the receiver are at the same elevation. • Evaluate (just once) the static solution associated with Gnn (0)
, but this
This is because the off-diagonal blocks are virtually zero on account of time integrating to infinity, either numerically or perhaps even in
the reasons considered in the previous section. It should be noted, closed form. The latter would be possible in the case of regular load
however, that the subtraction does not make the zero-frequency term of shapes such as a rectangular load. This step merely requires the
Gnn to be zero. The reason is that Gnn (0)
(k 0) Gnn (k 0, = 0) , inverse Fourier transforms of Eq. (36, 37) multiplied by the load
because nn neglects the effect of layer thicknesses and the coupling to
(0) terms. Add this contribution at the appropriate locations.
•
G
the distant layers. Only for large k is the cancelation of these two terms Add the static solution to the results of the numerical inversion.
virtually perfect. But this poses no problems for evaluating (38). This
has also the advantage that the Fourier inversion of Gnn (0)
depends only However, if the last step were done by purely numerical means,
on the material properties of the two adjoining half-spaces and not on then it might be worth to discard altogether the subtraction
their thicknesses, so the inversion can in some cases be accomplished in procedure involving Gnn (0)
and try integrating the tail numerically
closed form. as follows: Do a clockwise contour integration along the path
After the dynamic component u = um um (0)
has been obtained, it k max (k max i ) k max and since no poles are enclosed, set
behooves to add the static part to recover the full solution, which, as we the result to zero; also, by Jordan's lemma, set also the integral “around
already mentioned, is needed only when m = n . The contribution of infinity” to be zero; finally, equate the desired tail integral with the
this integral is not only rather straightforward to evaluate, but in some integral along the vertical complex branch, and reverse the limits as
cases can even be done in closed form, including the case when the soil may be appropriate. For example, consider an integrand of the form
should be layered —and this is because of the layer decoupling at high
wavenumbers already mentioned. In essence, this requires the static
1
k max k
e ikx dk = {Ci(k max x ) +i 2
Si(k max x ) }
solution for a load applied at the interface of two half-spaces joined at kmax i e i(kmax i ki) x
the location of the load. For example, setting GU = 0 in Eq. (37) and = i dki
k max k max i ki
assuming pn = e (a unit vector), it can readily be shown that these (kmax x )
e
expressions lead to the classical Cerruti and Boussinesq solutions for = ie ikmax x d
point loads applied onto the surface of an elastic half-space. The pro-
=0 1 i (39)
cedure is then as follows:
where on the first line we added the exact result in terms of the integral

• At the desired frequency

sine and cosine functions (sinint and cosint functions in Matlab) so as
, define a dense set of wavenumbers
to have a yardstick against which the numerical results could be com-
k = [0, k , k max 2 ], i.e. up to the maximum wavenumber k max in
pared. Because the exponential term in the final numerical integral dies
Eq. (34) times 2 , to account for the fact that the length of the set of
out quickly, just a few complex points in the integral will suffice to
radial wavenumbers must be somewhat longer to allow for projec-
obtain rather good numerical approximations for the tail. Indeed, trying
tions along k x , k y . This set of wavenumbers must be dense enough
out a simple Matlab script for the above integral, we obtained results to
that the response functions fij (k, ) will change smoothly from one
step to the next. This will allow intermediate values to be obtained an accuracy of 99.96% integrating for k max x = 1 in the interval
0 6 at steps = 0.05. Of course, in an actual problem one would
reliably be interpolation, if need be. The wavenumber step k
needed for this purpose will be discussed later on. not integrate the function 1/ k , but the actual components of Gnn pn , and

• Compute the Green's functions [f11 , f22 , f13 , f31 , f33 ]mn at the desired just a few integration points should suffice to estimate the tail; the nice
exception is the case when pn is a constant, in which case the above
pair of source-receiver stations (m , n) for the dense set of N radial
integral can be used directly, inasmuch as in the tail Gnn ~1/ k . The

667
E. Kausel Soil Dynamics and Earthquake Engineering 115 (2018) 663–672

preceding assumes also that k max exceeds the wavenumber for the frequency response function for seismograms in layered media, but in
slowest Rayleigh waves, i.e. k max > 1.2 / min . Care must also be ex- such a way that this damping could be extracted later on, i.e. fully re-
ercised to use the correct complex values for all variables involved, and moved after the transformation into the time domain was completed.
also modify the integration direction as may be needed depending on An extension of this technique to other dynamic systems, which goes by
the sign of the exponential terms. Finally, if the integrand includes the name of Exponential Window Method (EWM), was proposed by
trigonometric functions, these may need to be expressed in terms of Kausel and Roësset [7] and by Hall and Beck [8]. In essence, this
exponentials and integrated along different contours, if need be. These method constitutes a numerical implementation of the Laplace Trans-
are all, however, details which are left to the readers to sort out. form. It entails using complex frequencies of the form c = i ,
where is a small positive, constant quantity, say = 2 , with
1

5.1.1. Wavenumber step being the frequency step used to convert back into the time domain.
The wavenumber step k is a rather delicate matter. It is influenced This allows considering soils without damping and virtually eliminates
by all of the following considerations: wraparound and sampling issues. For further details, see Kausel [9],
Chapter 7.
• A discrete wavenumber summation automatically implies a periodic
source with a spatial period of L = 2 / k . Hence, the step must be 6. Examples
small enough to achieve a good separation to those “neighboring
sources”. If one is not interested in the Green's functions in the fre- Example 1. 2.5-D line load acting on the surface of an elastic half-
quency domain themselves, but only care to use these as tools to space.
obtain the time response, then the separation need only be as large
as the distance traveled by the fastest waves from the neighboring Consider a line load of magnitude Pz and spatial frequency k 0 that
source to the farthest receiver (maximum range) up to the maximum points upward (i.e. in the z direction) and acts along the y -axis at the
time of interest. But if the actual objective is the Green's functions in free surface z = 0 of an elastic half-space, That half-space has a shear
the frequency domain, one generally needs a larger distance to modulus µ , P-wave and S-wave velocities , , and Poisson's ratio .
prevent contamination, in which case one relies on damping and The load is thus of the form
geometric attenuation to effectively eliminate the contribution of p0 = { 0 0 pz }T (40)
pz (x , y ) = Pz (x ) e i k0 y,
those neighboring sources.
• Depending on the method used to compute the wavenumber in- where k 0 is the prescribed wavenumber of the load in y . Its Fourier
transform is
tegrals, it is necessary to use a sufficient number of points to resolve
the oscillations of the trigonometric functions that multiply the + +
kernels. Hence, the higher the frequency or the farther the range p˜ 0 (k x , k y ) = Pz p0 (x , y )e i kx x e i k 0 ydx dy
rmax (the largest distance to the receivers), the more points are +
= Pz {0 0 1}T e i k 0 ye i ky ydx dy
needed.
• The step must also be small enough so that all poles (or nearly all of = 2 Pz {0 0 1}T (k 0 k y) (41)
the poles) are properly resolved. By and large, as the frequency in-
creases, so does the number of poles, and the last pole entering into which is the same as saying that k y = k 0 is a fixed, known parameter.
the computation usually influences the results the most. In practice, Now, the explicit expressions for the flexibility functions of the half-
we do not know ahead of time how many poles exist at any given space in cylindrical coordinates are as follows ([3], Table 10.1, p. 146)
frequency, but their number can be estimated by determining the 1 1
f11 f13 1 s 2 (1 s2) ps (1 + s 2)
frequencies at which the stratified soil has near resonances when the 1
FSVP = KSVP = = 2

waves move vertically, i.e. k = 0 . f31 f33 2kµ ps

1
(1 + s2)
1
p 2 (1 s 2)
•
2
An alternative to direct numerical integration might be to find the
(42a)
real (or nearly real) poles at some frequency via search techniques,
evaluate their residues by contour integration, subtract their con- FSH = KSH1 = {f22 } =
1
tribution to the integral in (38), and use numerical integration solely ksµ (42b)
for the remainder (the branch integrals).
1
= ps 4
(1 + s 2) 2 = Rayleigh function (42c)
One final observation: unless the response is needed at many equi-
distant receiver stations, it is generally not convenient to use the Fast 2 2

Fourier Transform to carry out the integrations over wavenumbers. This p= 1 , s= 1

k k (42d)
is because a direct numerical integration for the response at a few re-
ceiver stations is usually more economical than using the FFT, and From Eq. (21), the response in the frequency-wavenumber domain is
because that allows also using uneven spacing of the receivers. In ad- then
dition, it also allows for uneven wavenumber steps, say smaller in the
u˜ x f11 cos2 + f22 sin2 (f11 f22 )sin cos i f13 cos
vicinity of resonances, and fewer elsewhere (i.e. adaptive integration, 0
u˜ y = 2 Pz (f11 f22 )sin cos f11 sin2 + f22 cos2 i f13 sin 0
Filon's method, Clenshaw-Curtiss quadrature, etc.). 1
u˜ z i f31 cos i f31 sin f33

5.2. Inversion into the time domain via complex frequencies i f13 cos
= 2 Pz i f13 sin
By and large, the frequency response functions considered in this f33 (43a)
section, i.e. the kernels in the frequency-wavenumber integrals, tend to
where
exhibit sharp peaks at the systems’ cutoff frequencies, and the number
of these peaks grows steadily with frequency. When little or no material k= k x2 + k y2 , cos = k x /k , sin = k y/k (43b)
damping (i.e. attenuation) is present, these peaks are sharp and tend to
make the numerical inversion (via Fourier transform) into the time In the spatial domain, the displacements on the surface of the half-
domain rather problematic. To avoid these difficulties, complex fre- space due to a spatially harmonic line load follow from the inverse
quencies were used by Phinney [6] as a means to artificially flatten the Fourier transform

668
E. Kausel Soil Dynamics and Earthquake Engineering 115 (2018) 663–672

i f13 cos displacements due to a rectangular load of intensity q acting in the x

direction, measured at the lower left corner A of the rectangular area is
+
u 0 (x , y, 0) = Pz e iky y i f13 sin e i kx x dk
x
f33 given by
(44)
q q (1 2 )q
which can readily be evaluated numerically. uAx = {I1 + (I2 I1)}, uAy = I3, uAz = ( ) I4
µ µ µ (48a)
Observe that in this example k y is fixed; if that were the only value
of k y for which the computation is done, then it might be advantageous b a a
I1 (a, b) = a sinh 1 + b sinh 1 , I2 (a, b) = b sinh 1
to sample k x instead and avoid the interpolation when discretizing the a b b (48b)
radial wavenumber k . But if the computation were to be repeated for
several wavenumbers k y , then it behooves to heed the earlier re- b 1 a 2
I3 (a, b) = a + b a2 + b2 I4 (a , b) = a arctan + 2 b ln 1 +
commendation concerning the discretization along the radial wave- a b
number. (48c)
Example 2. Rectangular load on the surface of an elastic half-space Because of symmetry, the horizontal displacement is equal at all four
Consider next a harmonic uniform rectangular load px0 of size corners. Hence, to obtain the horizontal displacement at any other interior
2a × 2b with net horizontal force in the x direction, Px = 4ab px0 acting point underneath the load, simply divide the surface load into four rec-
on the free surface of an elastic half-space with shear modulus µ and tangular sub-loads that meet at a common corner, and then superpose the
Poisson's ratio . Its double Fourier transform into the wavenumber four corner displacements for each of those loads. For exterior points, do
space is likewise, but combine two large positive loads (encompassing the actual
load) with two smaller negative loads (which are external to the actual
sin k x a sin k y b load) to model appropriately the applied load. We leave these well-known
px (k x , k y , ) = Px
kx a ky b (45) details for the reader to sort out. Still, the reader is encouraged to also
evaluate the static solution by numerical integration as described in earlier
while py = pz = 0 . On the other hand, the Green's functions fij at the
sections and compare that with the closed-form solution given above. This
surface of a homogeneous half-space due to SVP and SH sources placed
would provide a nice check to the formulation.
at the surface are as in the previous example
In particular, the elements of the static flexibilities of the half-space Example 3. Moving bell load
inferred from Eq. (37) are
We consider next a load in the form of a Gaussian bell that moves
1 with constant speed V along the y axis on the surface of a 3-layer
f11(0) = f 33
(0)
=
µk (46a) medium. Although the applied load exhibits cylindrical symmetry, that
1 does not help us in this case because as the load moves, its center moves
(0)
f 22 = along and the problem ceases to be axisymmetric. Now, although a
µk (46b)
cylindrical off-center formulation might be possible by expressing the
1 2 cylindrical functions in terms of infinite summations via Graf's theorem,
f13(0) = f 31
(0)
=( )
µk (46c) convergence would be exceedingly slow and the solution would involve
many terms in the azimuthal direction, because the resulting dis-
From Eq. (38), the partial response (i.e. without the static contribution)
placement field does not have cylindrical symmetry. Thus, we resort to
observed at the surface is then
a 2-D Cartesian formulation. The bell load is then
ux
x 2 + y2
( ) r2
1 2 +kmax +kmax
uy = Px (G11 1 1 a2
q (r , a ) = e a2 = e
2 kmax kmax
a2 a2 (49)
uz surface

(0) sin k x a sin k y b where the a parameter is a radial measure of how fast the distributed
G11 ) e e i kx x e i ky ydk
x dk y load decays in r . Indeed, the distributed load will be negligible by the
kx a ky b (47a)
time that r > a , because exp( 2) = 0.517 × 10 4 . This load will have
where e = (1 0 0) is a unit vector and
T
a net resultant of 1 because
2
(f11 f11(0) ) [(f11 f11(0) ) (f22 (0)
f 22 )] sin2 q (r , a) r dr d = 2 e
2
d = 1, = r /a
(0) 0 0 0 (50)
(G11 G11 )e = [(f11 f11(0) ) (f22 (0)
f 22 )]sin cos
(0) Its double Fourier transform into the wavenumber domain is
i (f31 f 31 )cos
2 2
(47b) Q (k x , k y , a ) = a
1 +
e ( ax ) ei kx xdx
a
1 +
e ( ay ) ei ky ydy
If the soil has some damping (i.e. attenuation), then the Rayleigh 2 2
ky a
pole and the P and S branch points are complex, in which case the =e
kx a
2 e 2 (51)
flexibility functions have no singularities for real wavenumbers and can
thus be integrated numerically. However, f11 , f13 , f33 will exhibit a sharp which again has the form of a Gaussian bell, and thus decays fast with
peak in the neighborhood of k = / CR (where CR is the speed of the two wavenumbers. The displacement in the time domain for a
Rayleigh waves), and f22 will be strongly peaked in that very same moving bell load can straightforwardly be show to be given by
neighborhood, namely k = / (because CR ).
u (x , =y Vt )
At this point we can proceed to add the contribution of the static 2
2 ky a
part, namely the displacements on the surface produced by the static, 1 +
iky x
+
= exp e G (k x , k y ,
uniformly distributed rectangular load. The closed-form solution to this 2 2
static problem was first reported by Schleicher [10] and extensions kx a 2
= k y V )exp e ikx dk x dk y
thereto given by Steinbrenner [11], which are contained in mis- 2 (52)
cellaneous textbooks such as Harr [12] or Florin and Schiffman [13].
These solutions can also be obtained directly by integration of the where = y Vt is a coordinate centered on the load that moves with
Cerruti and Boussinesq solutions. It can thus be shown that the it at the same speed.

669
E. Kausel Soil Dynamics and Earthquake Engineering 115 (2018) 663–672

We define a grid of size 30 × 30 meters, but because of symmetry

with respect to y we only display the results for the right half
of that grid, i.e. X = [0: 15, 15: 15]. We chose the radius of the bell
load to be a = 0.25 m, the spatial discretization to be at steps
x = y = 0.10 m (i.e. 150 intervals in each direction), and a spatial
period of L = 16x max = 240 m, implying a wavenumber step of
k = 2 / L = 0.0262 m 1. Also, the integral can safely be truncated at
the upper limit k max = 2 /a = 25.13. This implies in turn that the wa-
venumber domain has N = k max / k = L/ a = 960 intervals (i.e. 961
points). The material properties for this problem are given in Table 1.
Observe that the second sandwiched layer is softer than the surface
layer above it.
Fig. 1 shows the phase velocity spectrum for this layered profile (SV-
P modes), obtained with the same SMM formulation in the context of a
search of the real zeros of the determinant of the global stiffness matrix.
This figure displays the phase velocity for the various wave modes, the
first of which exhibits a local minimum at a phase velocity of Vc = 154
m/s, which can be referred to as the critical velocity. At this frequency,
the group and phase velocities become equal, as can easily be demon-
strated: Fig. 1. Dispersion spectrum for the 3-layer soil used in Example 3: Phase ve-
locity [m/s] vs frequency [Hz].
dVph d 1 dk 1 dk 1 Vph
= = = 1 = 1 to be ticks in the intensity scale on the right are actually negative signs.
d d k k k2 d k k d k Vgr
Each change in shade represents an order or magnitude decrease in
=0 Vph = Vgr motion intensity, with the largest intensity being on the order of
10 7.8 = 1.5849 × 10 8 . This is comparable to the displacement under-
More generally,
neath a circular load of unit intensity and radius a = 0.25 underlain by
Vph dVph an elastic half-space with the properties of the surface layer.
=1 , At 100 m/s, the displacement pattern is rather similar to that of a
Vgr Vph d
static, non-moving load. By the time we reach 150 m/s (which is still
Thus, when the slope is negative as it is in all curves except in the below but near the critical speed Vc = 154 m/s), we observe a sub-
first segment of the first dispersion curve, the ratio is greater than 1, stantial increase in the area of strong motion as well as in the dis-
which means that the phase velocity exceeds the group velocity. placement pattern. At 160 m/s, which exceeds modestly the critical
Beyond the local minimum, however, the slope turns positive and the speed, the motion is already rather intense throughout, and patterns of
group velocity becomes faster than the phase velocity, and if so, it starts wave fronts appear which evoke the Doppler Effect.
controlling the speed at which energy propagates in the medium. That
is, it is the first dispersion curve that controls the resonance effects
observed in the soil when the load moves at some speed, and especially
so past the critical velocity.
Observe that at higher frequencies, the various higher modes in-
terlace without crossing at the Rayleigh wave velocity of the surficial
layer, namely Vph = 184 m/s.
Figs. 2–5 show the response at the surface of the soil within a co-
ordinate grid that moves together with the load. In each case, the
computation is carried out over the entire grid of points and 961 wa-
venumbers, for a total computational time of about 440 s (~7 min) for
each load speed, and using the Matlab program Moving_Bell_Load [14]
running on a HP desktop with an Intel i7 3.2 GHz processor. A com-
plimentary copy of the software used herein is available at the link
provided, see the reference section. The program includes functions for
the assembly of the rigidity matrices for the layers, which can be ar-
bitrary in number, and each can have arbitrary material properties. The
program also enforces the required complex symmetries and anti-
symmetries exhibited by the flexibility functions with respect to posi-
tive and negative values of the horizontal wavenumbers.
The load speeds chosen are V =100, 150, 160, 220 m/s. The motion
intensity is displayed in a logarithmic color scale, even if color is
available only in the online version. Please observe also that what seem

Table 1
Soil Properties.
Layer [m/s] [m/s] [kg/m3] Depth Damping

1 200 346 2000 0.25 2 0.005 Fig. 2. Vertical motion intensity log uz ; load speed V= 100 m/s, quasi-static.
2 141 244 2000 0.25 3 0.005
The color intensity corresponds to the (negative) exponential scale on the right.
3 200 346 2000 0.25 Inf 0.005
Note: Colored figures only available in the online version.

670
E. Kausel Soil Dynamics and Earthquake Engineering 115 (2018) 663–672

Fig. 5. Vertical motion intensity log uz , load speed V= 220 m/s, transonic (i.e.
Fig. 3. Vertical motion intensity log uz , load speed V= 150 m/s, sub-critical. super-Rayleigh in the top layer). Angle of Mach cone is sin = VR / V = 184/220 .

inclination, and indirectly confirming the accuracy of the computation.

Observe that the motion behind the Mach cone is very strong in com-
parison with the motions observed at lesser speeds. Also, ahead of that
Mach cone the motion is some two orders of magnitude smaller than the
motion behind the wave front, even if it does not vanish altogether.
That motion is elicited by P waves that move faster than the load.

7. Conclusions

This article provided the technical details on a straightforward ex-

tension of the well-known Stiffness Matrix Method (SMM) for layered
soils, which is known to be valid only for either plane strain or for
cylindrical coordinates. The extension herein allows for loads with ar-
bitrary spatial distribution, so it is represents a generalization of the
method to the full 3-D space. However, inasmuch as the formulation
takes full advantage of the isotropy (or transverse isotropy) of the
medium, to allow calculating initially the Green's functions in cylind-
rical coordinates that only in the last step are converted to fully
Cartesian coordinates, the method herein cannot be used for ortho-
tropic or anisotropic media.
The method is highly effective and could readily be extended further
to allow evaluations of the strains and stresses at arbitrary point in the
medium.

Acknowledgements

The material described herein is the result of curiosity-driven, basic

research by the writer. This research has not received any funding from
Fig. 4. Vertical motion intensity log uz , load speed V= 160 m/s, super-critical.
Government, Industrial or Private Institutions. The writer wishes to
thank the three reviewers for their rather constructive and precise
By the time that we reach 220 m/s, which is above the Rayleigh comments.
wave speed of the surface layer (VR = 184 m/s), we finally observe the
clear manifestation of the first Mach cone with an angle defined by the References
ratio of the Rayleigh wave velocity and the load speed. The angle of
that cone is = arcsinVR/ V = 55.8o . Fig. 5 includes a segmented line [1] Kausel E, Röesset JM. Stiffness matrices for layered soils. Bull Seism Soc Am
inclined at this angle, confirming the correctness of the Mach cone 1981;71:1743–61.

671
E. Kausel Soil Dynamics and Earthquake Engineering 115 (2018) 663–672

[2] Kausel E, Peek R. Dynamic loads in the interior of a layered stratum: an explicit [8] Hall JF, Beck JL. Linear system response by DFT: analysis of a recent modified
solution. Bull Seismol Soc Am 1982;72(5):1459–81. [See also the errata in BSSA 74, method. Earthq Eng Struct Dyn 1993;22:599–615.
1984, p. 1508]. [9] Kausel E. Advanced structural dynamics. Cambridge University Press; 2017.
[3] Kausel E. Fundamental solutions in elastodynamics: a compendium. Cambridge [10] Schleicher FA. Zur Theorie des Baugrundes. Der Bauing 1926;48:931–52.
University Press; 2006. [949–952].
[4] Rokhlin S, Wang L. Stable recursive algorithm for elastic wave propagation in [11] Steinbrenner W. Tafeln zur Setzungsberechnung. Die Str 1934;1:121–4.
layered anisotropic media: stiffness matrix method. J Acoust Soc Am [12] Harr ME. Foundations of theoretical soil mechanics. McGraw Hill; 1966.
2002;112:822–34. [13] Florin VA. Fundamentals of soil mechanics. In: Schiffman Robert L, editor. Volume
[5] Rajapakse RKND, Senjuntichai T. Dynamic response of a multi-layered poroelastic 1, Translation from the Russian. Washington, D.C: National Technical Information
medium. Earthq Eng Struct Dyn 1995;24:103–722. Service, U.S. Department of Commerce; 1972.
[6] Phinney RA. Theoretical calculation of the spectrum of first arrivals in a layered [14] Moving_Bell_Load: A computer program to evaluate the dynamic response of a
medium. J Geophys Res 1965;70:5107–23. layered elastic half-space when subjected to a bell-shaped load that moves over the
[7] Kausel E, Roësset JM. Frequency domain analysis of undamped systems. J Eng surface at constant velocity. 〈https://www.dropbox.com/s/qxbvte9syhczat8/
Mech, ASCE 1992;118(4):721–34. Moving_Bell_Load.m?dl=0〉.

672