Propagation of partially coherent non-stationary random waves in a viscoelastic layered half-space

ARTICLE IN PRESS

Soil Dynamics and Earthquake Engineering 28 (2008) 305–320 www.elsevier.com/locate/soildyn

Propagation of partially coherent non-stationary random waves in a viscoelastic layered half-space Q. Gaoa,b, J.H. Lina, W.X. Zhonga, W.P. Howsonb,, F.W. Williamsb a

Department of Engineering Mechanics, State Key Laboratory of Structural Analysis of Industrial Equipment, Dalian University of Technology, Dalian 116023, PR China b Cardiff School of Engineering, Cardiff University, Cardiff CF24 3AA, Wales, UK Received 7 December 2005; received in revised form 18 June 2007; accepted 19 June 2007

Abstract This paper presents a highly accurate method based on the precise integration method (PIM) and on the pseudo excitation method (PEM). The method computes the propagation behaviour of partially coherent non-stationary random waves in a viscoelastic, transversely isotropic solid, which consists of a multi-layered soil resting on a homogeneous semi-infinite space. The excitation source is a local rupture between two layers, which causes a partially coherent non-stationary random field. The analysis of non-stationary random wave propagation is transformed into that for deterministic waves by using PEM. The resulting governing equations in the frequencywavenumber domain are linear ordinary differential equations, which are solved very precisely by using PIM. The evolutionary power spectral densities of the ground level responses are investigated and some typical earthquake phenomena are explained. r 2007 Elsevier Ltd. All rights reserved. Keywords: Layered half-space; Precise integration method; Pseudo excitation method; Non-stationary random waves; Wave propagation

1. Introduction Wave propagation in layered media has been studied extensively in many fields, with geophysics and earthquake engineering being prominent, as follows. The theory of wave propagation in elastic media has been studied for a very long time [1–7] and latterly progressively more complicated and realistic models have been developed, e.g. to cover wave propagation for anisotropic layered solids, viscoelastic layered media, porous material, inhomogeneous solids, random material and random loads. Relevant detail of some of this work is now presented. The theory and phenomena of wave propagation in anisotropic stratified media has been studied extensively [8–13] and transient wave propagation in viscoelastic media has also been investigated [14]. A method for simulating the propagation of elastic waves in stratified transversely Corresponding author. Tel.: +44 2920874263; fax: +44 2920874597.

E-mail addresses: [email protected] (Q. Gao), [email protected] (J.H. Lin), zwoffi[email protected] (W.X. Zhong), [email protected] (W.P. Howson), [email protected] (F.W. Williams). 0267-7261/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.soildyn.2007.06.003

isotropic inhomogeneous media has been proposed by Zhang and Li [15] in which a generalized Thomson–Haskell matrix method is used. The theory of wave propagation in porous materials was initiated by Biot [16,17], since then other papers have been published [18–24]. The random properties of waves propagating in layered materials can arise from several sources [25] and have been studied for the cases when they are caused by irregular interfaces [26–28], by random material properties [29–34] and by stationary random waves propagating in isotropic and anisotropic layered solids [35,36]. However, to date there is no published work about the propagation of partially coherent non-stationary random waves. This paper covers non-stationary waves, which are propagating in a transversely isotropic, viscoelastic solid, which consists of a layered soil resting on a homogeneous semi-infinite space. The non-stationary random excitations are assumed to be partially coherent and to be caused by discontinuity of the displacements and/or stresses at the bottom of a specified layer, so that they can be written as the product of stationary random excitations and a slowly varying modulation function [37]. The pseudo excitation

ARTICLE IN PRESS 306

Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

method (PEM) is used to solve this problem because it is an efficient and accurate algorithm for linear random vibration of systems, which has been used successfully for random analyses of many structures, e.g. buildings, bridges, dams and vehicles [38–41]. However, PEM has had to be generalized in this paper in order to deal with non-stationary random wave propagation. This simplified the problem to a remarkable extent, because the nonstationary random response analyses are transformed into much easier deterministic transient dynamic analyses [42–49]. The resulting deterministic wave motion equations are then transformed into first order ordinary differential equations (ODEs) in the frequency–wavenumber domain by using Fourier transformation. In this paper, the ODEs are solved numerically by using the precise integration method (PIM) to achieve very good accuracy, great economy and stability [50,51]. The validity and high accuracy of the proposed method are demonstrated by comparisons with analytical solutions, and some typical earthquake phenomena are also explained. 2. The governing equations for non-stationary random waves In this paper, the materials are assumed to be transversely isotropic, so that the parameters of the soil in the horizontal plane are independent of direction. Therefore, see Fig. 1, the horizontal direction of propagation of the plane waves can be taken as being the x axis, with the z axis pointing downwards such that z ¼ 0 at the free surface, i.e. at ground level. Hence all variables are independent of the coordinate y. As can be seen: the rth layer (r ¼ 1, 2, y l) is bounded by the horizontal planes z ¼ zr1 and z ¼ zr; the (l+1)th layer is the homogeneous semi-infinite space and; the non-stationary random source is located at the bottom of layer s (0ospl). The first part of the theory below is in the space–time domain but the final part is in the frequency–wavenumber domain. Therefore the space–time domain is indicated by putting a tilde above symbols. Hence let u~ m ðx; z; tÞ; ðm ¼ x; y; zÞ be the displacements along the inertia coordinates (x, y, z) and t~ mn ðx; z; tÞ; ðm; n ¼ x; y; zÞ be the components of the stress tensor. Then the equations of wave motion can

be written as q~txx q~txz rq2 u~ x þ ¼ ; qx qz qt2 q~tzx q~tzz rq2 u~ z þ ¼ , qx qz qt2

q~tyx q~tyz rq2 u~ y þ ¼ , qx qz qt2 ð1Þ

in which t is time and r is the density, which may differ from layer to layer. Using ~mn ðm; n ¼ x; y; zÞ to denote components of the strain tensor, the strain–displacement relationship is (2) ~mn ¼ 12ðu~ m;n þ u~ n;m Þ. Now let the viscoelastic isotropic stress-strain relationships be k1 X

k

2 X dk dk ð~tmn Þ ¼ qk k ðldmn y~ þ 2G~mn Þ, k dt dt k¼0 k¼0 ( 1; m ¼ n; y~ ¼ ~ xx þ ~yy þ ~ zz ; dmn ¼ 0; man;

pk

ð3Þ

where l and G are the Lame´ constants, which may differ from layer to layer and; pk and qk are the viscoelastic material constants. Hence: p0 ¼ 1, p16¼0 and q16¼0 for a Maxwell fluid; p0 ¼ 1, q06¼0 and q16¼0 for a Kelvin solid and; p0 ¼ 1, q06¼0, p16¼0 and q16¼0 with q14p1  q0 corresponds to a three-parameter solid [52]. The free boundary conditions at ground level (i.e. at z ¼ z0) are t~ mz ðx; 0; tÞ ¼ 0;

ðm ¼ x; y; zÞ,

(4)

Usually, in numerical analysis, the original sources are substituted by the ‘equivalent force distribution’, which is confined to a rather limited region with small explosions. In this way, we can represent such a source via a singular displacement and stress distribution at the bottom of the sth layer (0ospl), i.e. at z ¼ zs [5,7], then the continuity conditions at the interfaces are u~ m and t~ mz ðm ¼ x; y; zÞ are continuous at zr ðr ¼ 1; 2; . . . ; s  1; s þ 1; . . . ; lÞ

ð5Þ

and in addition ~ m ðx; z u~ m ðx; zþ ¯ m ðx; tÞ s ; tÞ ¼ u s ; tÞ þ u þ  t~ mz ðx; zs ; tÞ ¼ t~ mz ðx; zs ; tÞ þ t¯ mz ðx; tÞ ðm ¼ x; y; zÞ,

ð6Þ

where u¯ m ðx; tÞ and t¯ mz ðx; tÞ are the discontinuity functions  of the displacements and stresses, and zþ s and zs represent the lower and upper faces of the interface zs. In this paper it is assumed that u¯ m ðx; tÞ and t¯ mz ðx; tÞ are random fields, which are non-stationary in the time domain and have their attenuation functions given in the spatial domain. The initial conditions are qu~ m ðx; z; 0Þ ¼ 0 ðm ¼ x; y; zÞ. (7) qt Radiative conditions in the homogeneous semi-infinite space (i.e. layer l+1) are also required and they are presented in Section 4. u~ m ðx; z; 0Þ ¼

Fig. 1. Layered half-space.

ARTICLE IN PRESS Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

3. Pseudo excitation method for non-stationary random waves PEM is an efficient and accurate algorithm for solving linear random vibration problems and has been successfully used in the random analyses of many structures [38–41]. Its key idea is the introduction of pseudo deterministic excitations to transform the stationary or non-stationary random vibration analyses into much simpler harmonic or transient dynamic ones [42–49]. For brevity, the discontinuity functions in the source conditions (6) are written as the vector fðx; tÞ ¼ f¯ux ðx; tÞ; u¯ y ðx; tÞ; u¯ z ðx; tÞ; t¯ xz ðx; tÞ; t¯ yz ðx; tÞ; t¯ zz ðx; tÞgT ,

307

In which E[#] represents taking the expected value of # and the superscript H denotes Hermitian transpose. The relationship between the correlation function matrix ¯ ¯ RðtÞ and the PSD matrix SðoÞ for the stationary random vector fs(t) is Z þ1 ¯ ¯ RðtÞ ¼ E½f s ðt1 Þf H SðoÞ exp½ioðt1  t2 Þdo. ðt Þ ¼ 2 s 1

(12) Then Eq. (11) can be written as Z þ1 H ¯ Rðx; 0; t; tÞ ¼ Jðx; o; tÞSðoÞJ ðx; o; tÞ do

(13)

1

in which

(8)

Z tZ

þ1

Gðx; t  tÞAðx; tÞ expðiotÞ dx dt.

in which f(x, t) is a non-stationary random process vector in the time domain, and the superscript T represents the transpose of a matrix or vector. The non-stationary random process can usually be written as the product of a stationary random process and a deterministic slowly varying modulation function [37]. If its six non-stationary random discontinuity functions are assumed to have the same modulation function, Eq. (8) gives

Jðx; o; tÞ ¼

~ tÞf s ðtÞ, fðx; tÞ ¼ Aðx;

in which J(x, o, t) is a 3  6 matrix, such that its element Jij(x, o, t) represents the ith ground level displacement component caused when f(x, t) is null except that its jth element is assumed to be A(x, t)exp(iot). Computing S(x, 0, o, t) from Eq. (15) directly is very time consuming, which is why PEM [42–49] is used instead. The PSD matrix of the stationary component fs(t) of the non-stationary random f(t) is a known Hermitian matrix ¯ SðoÞ, of which the eigenvalues are generally non-negative real numbers. Now assume that the number of non-zero ¯ eigenvalues of SðoÞ is a¯ and let la and Ha ða ¼ 1; 2; . . . ; a¯ Þ be the eigenvalues and normalized eigenvectors, i.e.

(9)

~ tÞ is a deterministic function and fs(t) is the where Aðx; corresponding stationary random process. Since the elements of vector fs(t) are mutually partially coherent, ¯ the rank a¯ of the PSD matrix SðoÞ of fs(t) should satisfy 1o¯ap6. The aim of the present paper is to compute the non-stationary PSD matrix of random ground level ~ tÞ and SðoÞ ¯ response when Aðx; are given for the nonstationary source vector f(x, t). Assume that G(x, t) is a 3  6 Green function matrix such that its element Gij(x, t) represents the ith ground level displacement component caused when f(x, t) is null except that it has a unit impulse as its jth element. Since the governing equations are a linear system with constant coefficients, the random responses induced by source conditions (9) can be obtained by using the Green function matrix G(x, t) as Z t Z þ1 q~ ðx; 0; tÞ ¼ Gðx; t  tÞfðx; tÞ dx dt 0 1 Z t Z þ1 ¼ Gðx; t  tÞAðx; tÞf s ðtÞ dx dt, 0

1

q~ ðx; z; tÞ ¼ fux ðx; z; tÞ; uy ðx; z; tÞ; uz ðx; z; tÞgT .

ð10Þ

Then the correlation function matrix of the random ground level displacement responses can be computed by Rðx; 0; t1 ; t2 Þ ¼ E½~qðx; 0; t1 Þ~qH ðx; 0; t2 Þ Z t1 Z þ1 Z t2 Z þ1 Gðx; t1  t1 ÞAðx1 ; t1 Þ ¼ 0

1

0

1

H E½f s ðt1 Þf H s ðt2 ÞAðx2 ; t2 ÞG ðx; t2  t2 Þ dx1 dt1 dx2 dt2 .

0

(14) According to random theory [53,54], the PSD matrix S(x, 0, o, t) of the ground level random displacement responses q~ ðx; 0; tÞ is the integrand of the right-hand side of Eq. (13), i.e. H ¯ Sðx; 0; o; tÞ ¼ Jðx; o; tÞSðoÞJ ðx; o; tÞ,

¯ SðoÞH HH a ¼ l a Ha ; a Hb ¼ dab ; ¯ Then SðoÞ can be written as ¯ SðoÞ ¼

a¯ X

a; b ¼ 1; 2; . . . ; a¯ .

la H a H H a.

(15)

(16)

(17)

a¼1

Now let the following pseudo excitations be constructed at z ¼ zs [50,55] pffiffiffiffiffi ~ tÞ expðiotÞ la Ha ; a ¼ 1; 2; . . . ; a¯ . (18) f~ a ðx; tÞ ¼ Aðx; Then according to Eqs. (10) and (14), the pseudo displacement responses of the ground level (z ¼ 0) excited by the pseudo excitations f~ a ðx; tÞ are pffiffiffiffiffi (19) q~ a ðx; 0; tÞ ¼ Jðx; o; tÞ la Ha ; a ¼ 1; 2; . . . ; a¯ . Thus a¯ X a¼1

ð11Þ

1

q~ a q~ H a

¼ Jðx; o; tÞ

a¯ X

! la H a H H a

a¼1 H ¯ ðx; o; tÞ.  JH ðx; o; tÞ ¼ Jðx; o; tÞSðoÞJ

ð20Þ

ARTICLE IN PRESS Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

308

Comparison of Eqs. (15) and (20) gives Sðx; 0; o; tÞ ¼

a¯ X

q~ a q~ H a.

(21)

a¼1

It is seen from the above that q~ a is the deterministic response excited by the deterministic excitation f~ a ðx; tÞ. Therefore, the key idea of PEM is to transform the nonstationary random wave problem into a series of deterministic transient wave problems. Thus the main steps when applying PEM to non-stationary random wave problems are: (1) ¯ decompose the PSD matrix SðoÞ of the stationary random components of the non-stationary random excitations according to Eqs. (16) and (17); (2) construct pseudo excitations f~ a ðx; tÞ by using Eq. (18); (3) compute the pseudo responses q~ a ðx; 0; tÞ due to f~ a ðx; tÞ and; (4) compute the PSD matrix S(x, 0, o, t) of the non-stationary random responses of the ground level displacement by using Eq. (21). In the following sections, the pseudo responses will be computed by using PIM. 4. The state equations in frequency–wavenumber domain The Fourier transformation and the inverse Fourier transformation are defined as, respectively, Z þ1 Z þ1 ~ z; tÞ exp½iðkx gðk; gðx; dt, ¯ ¼ ¯ ¯ z; oÞ ¯  otÞdx 1 1 Z þ1 Z þ1 1 ~ z; tÞ ¼ 2 gðk; k¯ do, gðx; ¯ exp½iðkx ¯ ¯ ¯ z; oÞ ¯  otÞd 4p 1 1 ð22Þ where g~ represents any displacement, strain or stress, e.g. u~ x ; ~xy ; s~ xy , etc.; o ¯ is angular frequency and; k¯ is wavenumber. Through the Fourier transformation, the time–spatial coordinates t and x are transformed into the frequency–wavenumber coordinates o ¯ and k, ¯ so that the coordinate system (x, z, t) becomes ðk; ¯ ¯ z; oÞ. Performing the Fourier transformation for Eqs. (1)–(3), the governing equation for isotropic solids leads to the following two first order linear ODEs: " # ( ) qja Aj Dj 0 vja ¼ Hj vja ; Hj ¼ ; vja ¼ ; ðj ¼ 1; 2Þ pja Bj Cj (23) 0

in which ð#Þ ¼ qð#Þ=qz, and q2a ¼ fuxa ; uza gT ;

q1a ¼ uya ;

p1a ¼ tyza ;

p2a ¼ ftxza ; tzza gT ,

(24) A1 ¼ C1 ¼ 0;

B1 ¼ mG  ro ¯ 2; "

A2 ¼

CH 2

¼ ik¯

" 4GðlþGÞ B2 ¼ mk¯

2

lþ2G

0

0

1

l lþ2G

0 0

0

#

D1 ¼

1 , mG 2 1

;

#  ro ¯ 2 I2 ,

1 G D2 ¼ 4 m 0

0 1 lþ2G

Pk 2 ðioÞ ¯ k qk , (26) m ¼ Pk¼0 k1 ¯ k pk k¼0 ðioÞ pffiffiffiffiffiffiffi where i ¼ 1; I2 is the (2  2) unit matrix; q1a and q2a are the displacement vectors in the frequency-wavenumber domain; p1a and p2a are the dual stress vectors in the frequency–wavenumber domain and; m is the damping factor of the viscoelastic material. The ODEs with j ¼ 1 correspond to the SH wave and the ODE with j ¼ 2 corresponds to the P-SV wave. Subscript a represents the ath pseudo excitation. Note that Hj in Eq. (23) is independent of a. Transforming the boundary and interface conditions of Eqs. (4) and (5), as well as the pseudo source conditions of Eq. (18), into the frequency–wavenumber domain, gives: the free surface condition at ground level (i.e. at z ¼ z0 ¼ 0) pja ðk; ¯ oÞ ¼ 0; ¯ 0; o; the interface continuity condition qja ðk; ¯ oÞ; ¯ z; o;

pja ðk; ¯ oÞ are continuous at ¯ z; o;

zr ðr ¼ 1; 2; . . . ; s  1; s þ 1; . . . lÞ;

ð28Þ

and the source conditions qja ðk; ¯ oÞ ¼ qja ðk; ¯ oÞ þ Aðk; ¯  oÞsja , ¯ zþ ¯ z ¯ o s ; o; s ; o; pja ðk; ¯ oÞ ¼ pja ðk; ¯ oÞ þ Aðk; ¯  oÞtja , ¯ zþ ¯ z ¯ o s ; o; s ; o;

ð29Þ

in which pffiffiffiffiffi la fHa1 ; Ha3 gT ; pffiffiffiffiffi t1a ðk; ¯ oÞ ¼ la fHa4 ; Ha6 gT ; ¯ o;

s1a ðk; ¯ oÞ ¼ ¯ o;

pffiffiffiffiffi la Ha2 , pffiffiffiffiffi t2a ðk; ¯ oÞ ¼ la Ha5 , ¯ o; s2a ðk; ¯ oÞ ¼ ¯ o;

ð30Þ where Hai ði ¼ 1; 2; . . . ; 6Þ represents the ith component of Ha. In the frequency and wavenumber domain, s1a ðk; ¯ oÞ and t1a ðk; ¯ oÞ are the pseudo excitations for ¯ o; ¯ o; SH non-stationary random waves and s2a ðk; ¯ oÞ ¯ o; and t2a ðk; ¯ oÞ are the pseudo excitations for P–SV ¯ o; non-stationary random waves. They are functions of k; ¯ and o, where k¯ and o ¯ are the wavenumber and ¯ o frequency which resulted from the Fourier transformation of the space–time coordinates x and t. Here o represents the frequency, i.e. the abscissa of the PSD curves of the random excitations. The wave motion equations should also satisfy the radiative conditions in homogeneous semi-infinite space, which are independent of the source conditions, and so are independent of a. Therefore, the state equation for the homogeneous semi-infinite space can be written as ¯ j v¯ j v¯ 0j ¼ H

3

(27)

ðj ¼ 1; 2Þ.

(31)

¯ 1 is a 2  2 matrix, whereas for For SH waves j ¼ 1 and H ¯j ¯ 2 is a 4  4 matrix. Let T¯ j and K P-SV waves j ¼ 2 and H (of which Appendix A gives details) be the eigenvector and ¯ j and let eigenvalue matrices of H

5,

ð25Þ

¯ 1 b¯ j ¼ T ¯j . j v

(32)

ARTICLE IN PRESS Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

Then vectors b¯ j must satisfy 0 b¯ j

¯ j b¯ j , ¼K

(33)

which has the solution ¯ j ðz  zl Þb¯ j ðzl Þ; b¯ j ðzÞ ¼ exp½K

zXzl

(34)

¯ j and T¯ j ordered as in Appendix A. For SH waves, with K b¯ 1 ðzl Þ is a two element vector with its first and second components representing, respectively, waves travelling upwards and downwards. Similarly, for P-SV waves, b¯ 2 ðzl Þ is a 4 element vector and its first two components correspond to waves travelling upwards while its last two components correspond to waves travelling downwards. Furthermore, according to Eq. (32), the state vector at z ¼ zþ l can be written as " #( ) ¯ U;j ðzl Þ TUU;j TUD;j b vj ðzþ , (35) l Þ ¼ T TDD;j b¯ D;j ðzl Þ DU;j in which TUU;j ; TUD;j ; TDU;j and TDD;j are obtained by ¯ j . The radiative conditions require that no partitioning T upwards travelling waves exist. Hence b¯ U;1 ðzl Þ ¼ 0 for SH waves and b¯ U;2 ðzl Þ ¼ 0 for P-SV waves, so that, by using Eq. (35), ¯ qj ðzþ l Þ ¼ TUD;j bD;j ;

¯ pj ðzþ l Þ ¼ TDD;j bD;j .

(36)

All governing equations and boundary conditions for non-stationary random waves in the frequency and wavenumber domain have now been obtained. Hence if qa ðk; ¯ o; zÞ is the solution of the state Eq. (23) with ¯ o; boundary conditions (27), (28) and (36) and pseudo excitations (29), then the inverse Fourier transformation can be used to obtain the pseudo response at ground level as Z þ1 Z þ1 1 q~ a ðx; t; o; 0Þ ¼ q ðk; ¯ o; 0Þ ¯ o; 4p2 1 1 a  exp½iðkx k¯ do. ð37Þ ¯ ¯ ¯  otÞd Because the non-stationary random source is transformed into deterministic pseudo excitation and the modulation function is zero when the coordinate tends to infinity, the Dirichlet condition is always satisfied and thus the Fourier transform of Eq. (37) always exists. However, in general, an analytical solution to the inverse Fourier transformation of Eq. (37) is not available. Therefore, instead, discrete Fourier transformation is used to solve Eq. (37). If the integration ranges used for the wavenumber and for frequency are, respectively, k¯ 0 to k¯ f and o ¯ 0 to o ¯f, and if the corresponding integration steps are Dk¯ and Do, ¯ Eq. (37) can be written as q~ a ðx; t; o; 0Þ 

¯ f =Do o k¯ f =Dk¯ X¯ DkD ¯ X ¯ o q ðnDk; ¯ o; 0Þ ¯ mDo; 4p2 n¼k¯ =Dk¯ m¼o¯ =Do¯ a 0

0

 expðinDkxÞ ¯ ¯ expðimDotÞ.

ð38Þ

Hence, the PSDs of the non-stationary random ground level displacement responses can be computed from Eq. (21).

309

For isotropic solids, Kennett [7] proposed an algorithm for solving the wave motion ODEs. It decomposes the wavefield into upgoing and downgoing parts and uses the reflection and transmission matrices at the interface. However, the application of this method to anisotropic solids is not easy, because it needs the eigenvalues of the ODEs, which are very difficult to extract. Therefore PIM is used below to solve the ODEs, because it does not require the eigenvalues. The wave motion equations for P-SV and SH waves have similar mathematical forms and so can be solved in the same way when using PIM, as follows, in which the subscripts j and a are omitted. 5. The precise integration method for the state equations PIM is an accurate method for solving ODEs with twopoint boundary value conditions or initial value conditions [50,51] and can be used for anisotropic solids [8,35]. To solve Eq. (31), an interval [za, zb] is selected within an arbitrary layer and it is assumed that the displacements and stresses are continuous in this interval. Let qa be the displacement vector at za and pb be the force vector at zb. For linear systems, it is known [50,51] that qa ¼ Fðzb ; za Þqb  Gðzb ; za Þpa ;

pb ¼ Qðzb ; za Þqb þ Eðzb ; za Þpa ,

(39) in which F, Q, G, and E are complex matrices to be determined. Consider two adjacent intervals [za, zb] and [zb, zc] and assume that the displacements and stresses are continuous in the intervals [za, zb] and [zb, zc]. Then from Eq. (39) qa ¼ F 1 qb  G 1 pa ;

pb ¼ Q1 qb þ E1 pa ,

(40)

qb ¼ F 2 qc  G 2 pb ;

pc ¼ Q2 qc þ E2 pb ,

(41)

The intervals [za, zb] and [zb, zc] can be merged into a single interval [za, zc], for which qa ¼ F12 qc  G12 pa ;

pc ¼ Q12 qc þ E12 pa

for

½zb ; zc , (42)

Comparing Eqs. (40) and (41) with Eq. (42) gives [50,51] F12 ¼ F1 ðI þ G2 Q1 Þ1 F2 ;

Q12 ¼ Q2 þ E2 Q1 ðI þ G2 Q1 Þ1 F2 ,

E12 ¼ E2 ðI þ Q1 G2 Þ1 E1 ;

G12 ¼ G1 þ F1 G2 ðI þ Q1 G2 Þ1 E1 .

ð43Þ The interval matrices F, Q, G and E for each layer can be computed by using the dual Eq. (23) and PIM, see Refs. [8,35,50,51] for details, so that the pseudo responses can be computed as follows. First, combining the interval matrices of layer 1 with those of layer s using Eq. (43) gives Q0s ; G0s ; F0s and E0s . Secondly, combining the interval matrices of layer s+1 with those of layer l in similar fashion gives Qsl ; Gsl ; Fsl and Esl . However, intervals [z0, zs] and [zs, zl] cannot be combined directly because there are source conditions at interface zs. These are  qþ s ¼ qs þ s;

 pþ s ¼ ps þ t.

(44)

ARTICLE IN PRESS Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

310

Now the equations for intervals [z0, zs] and [zs, zl] are qþ 0

þ ¼ F0s q s  G0s p0 ;

þ þ qþ s ¼ Fsl ql  Gsl ps ;

 þ p s ¼ Q0s qs þ E0s p0 , þ þ p l ¼ Qsl ql þ Esl ps ,

(45) (46)

so that substituting Eq. (44) into Eqs. (45) and (46), then combining intervals [z0, zs] and [zs, zl] gives þ þ ¯ 0l ; qþ 0 ¼ F0l ql  G0l p0 þ q

þ þ pþ ¯ 0l l ¼ Q0l ql þ E0l p0 þ p

(47) Here F0l ; G0l ; Q0l and E0l are obtained by combining Q0s ; G0s ; F0s and E0s with Qsl ; Gsl ; Fsl and Esl by using Eq. (43), and q¯ 0l ¼  F0s ½s þ ðI þ Gsl Q0s Þ1 Gsl ðt  Q0l sÞ, p¯ 0l ¼ Esl ðI þ Q0s Gsl Þ1 ðt  Q0s sÞ.

ð48Þ

According to the free boundary condition of Eq. (27) at z0 and the radiative conditions of Eq. (36) for the homogeneous semi-infinite space, we have 1 qþ ¯ 0l þ q¯ 0l . 0 ¼ F0l TUD ðTDD  Q0l TUD Þ p

(49)

This enables the pseudo ground level response to be obtained by using Eq. (38), after which the ground level displacement PSD matrix can be obtained by using Eq. (21). 6. Numerical examples In the following examples, the viscoelastic parameters used were p0 ¼ 1.0, p1 ¼ 0.05, q0 ¼ 1 and q1 ¼ 0.1, while the integral limits and step sizes used in the discrete inverse Fourier transformation were k¯ 0 ¼ 0:005, k¯ f ¼ 0:005, Dk¯ ¼ 2  106 m1 , o o and ¯ 0 ¼ 10:0, ¯ f ¼ 10:0 Do ¯ ¼ 0:01 s1 . In addition, it was assumed that the modulation function of the non-stationary random source ~ tÞ ¼ g^ x ðxÞg~ t ðtÞ and that the non-stationary (9) was Aðx; random source can be rewritten as fðx; tÞ ¼ g^ x ðxÞg~ t ðtÞf s ðtÞ,

(50)

where: f(x, t) is the non-stationary random source vector; fs(t) is the corresponding stationary random vector; g~ t ðtÞ is the modulation function of the non-stationary random source vector and; g^ x ðxÞ is the distribution function of the source intensity in the space domain. Fig. 2 shows the g^ x ðxÞ and g~ t ðtÞ used, which are 8 0:2t; 0pto5; > > >   < 2 1; 5pto35; x ; g~ t ðtÞ ¼ g^ x ðxÞ ¼ exp  6 > 0:2t þ 8; 35pto40; 4  10 > > : 0; 40pt: (51) Example 1. consists of non-stationary random SH waves for a single layered soil, which overlies a homogeneous semi-infinite space. The shear modulus and density are: G 1 ¼ 0:35  1011 N=m2 and r1 ¼ 3:0  103 kg=m3 for the layer, which has thickness L ¼ 1:9  104 m; and G 2 ¼ 0:15  1011 N=m2 and r2 ¼ 3:0  103 kg=m3 for the semi-

Fig. 2. The profile for functions g~ t ðtÞ and g^ x ðxÞ.

infinite space. The non-stationary random source is a rupture along the y direction at z ¼ L, i.e. between the layer and the semi-infinite space, and the discontinuous functions are displacement v¯ and/or stress t¯ yz , which are assumed to be non-stationary random processes. Three alternative PSDs of the stationary random components of v¯ and t¯ yz are considered, see Cases 1–3 in the caption of Fig. 3, which gives the corresponding autoPSDs of the random ground level response of the y direction displacements at the observation locations x ¼ 0.0, 10.0, 20.0 and 30.0 km. It can be seen that the random response PSDs are all zero for a certain initial time interval, the size of which depends on the distance from the observation position to the source. The results for Case 1 show that if the discontinuity of the source consists solely of a displacement discontinuity, the random response PSDs of ground level displacement approach steady state after a time, which is long if x is large or the related frequency is low. In contrast, the results for Case 2 show that if the only discontinuity of the source is a stress one, the random response PSDs of ground level displacement no longer approach a steady state. This is because, for the specified parameters, the stress PSDs were distributed mainly in a very low frequency range, for which decay is very slow. Finally, the results of Fig. 3 for Case 3 show that when the source involves both displacement and stress discontinuities, the response possesses a combination of the characteristics exhibited by Cases 1 and 2. The reasons for these phenomena are both due to the viscoelastic stress–strain relationship used, in which the damping is bigger for higher frequencies, and in addition, even if possessing equal damping, the higher frequency components will decay more quickly. Therefore, for a given time interval the components with higher frequency will attenuate more quickly. Fig. 3 also shows that the random response PSDs at the earthquake epicentre (x ¼ z ¼ 0) affect the widest range in the frequency domain and that

ARTICLE IN PRESS Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

the effect range concentrates at the lower frequencies when the observation location is distant from the earthquake epicentre.

311

For this simple problem, one of the main aims is to check the accuracy of the results obtained by the proposed method by comparing them with an analytical formula. If

Fig. 3. Non-stationary random response PSDs of v at locations x ¼ 0, 10, 20 and 30 km for the following three excitations. Case 1: S¯ v¯ v¯ ¼ 1:0 m2 s and S¯ v¯ t¯ yz ¼ S¯ t¯ yz v¯ ¼ S¯ t¯ yz t¯ yz ¼ 0; Case 2: S¯ v¯ v¯ ¼ S¯ v¯ t¯ yz ¼ S¯ t¯ yz v¯ ¼ 0 and S¯ t¯ yz t¯ yz ¼ 1:0 N2 s=m4 ; Case 3: S¯ v¯ v¯ ¼ 1:0 m2 s, S¯ v¯ t¯ yz ¼ conjðS¯ t¯ yz v¯ Þ ¼ ð3:0 þ 4:0iÞ  105 Ns=m and S¯ t¯ yz t¯ yz ¼ 1:0  1012 N2 s=m4 . (a) Case 1, x ¼ 0; (b) Case 1, x ¼ 10 km; (c) Case 1, x ¼ 20 km; (d) Case 1, x ¼ 30 km; (e) Case 2, x ¼ 0; (f) Case 2, x ¼ 10 km; (g) Case 2, x ¼ 20 km; (h) Case 2, x ¼ 30 km; (i) Case 3, x ¼ 0; (j) Case 3, x ¼ 10 km; (k) Case 3, x ¼ 20 km and (l) Case 3, x ¼ 30 km.

ARTICLE IN PRESS Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

312

Fig. 3. (Continued)

the two pseudo excitations of the source are taken as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S v¯ t¯ yz S t¯ yz v¯ S v¯ t¯ yz , s1;2 ¼ S v¯ t¯ yz S t¯ yz v¯ þ ðl1;2  Sv¯ v¯ Þ2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Sv¯ t¯ yz S t¯ yz v¯ ðl1;2  S v¯ v¯ Þ, t1;2 ¼ Sv¯ t¯ yz S t¯ yz v¯ þ ðl1;2  S v¯ v¯ Þ2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðS v¯ v¯ þ S t¯ yz t¯ yz Þ  ðSv¯ v¯  St¯ yz t¯ yz Þ2 þ 4Sv¯ t¯ yz St¯ yz v¯ l1;2 ¼ ð52Þ 2 then the analytical solutions of the pseudo ground level displacements are (see Appendix B)

For the third source condition, the variances of ground level displacement at observation locations x ¼ 0.00, 10.0, 20.0 and 30.0 km were computed from the analytical formula of Eq. (52) and by using the proposed method, see Fig. 4. Clearly, the results given by the proposed method agree excellently with the analytical ones. In fact the accuracy is much higher than can be seen, because the relative errors are all of the order 1010. Fig. 4 also shows that the variances are zero initially within an initial time interval, which becomes longer as the distance between the observation location and the earthquake epicentre increases.

v~1;2 ðx; o; t; 0Þ Z Z 1 þ1 þ1 gx ðkÞg ¯  oÞðmG 2 b2 s1;2 þ t1;2 Þ expðb1 LÞ exp½iðkx ¯ ¯ t ðo ¯  otÞ dk¯ do, ¼ 2 ¯ 2p 1 1 m½ðG 2 b2 þ G 1 b1 Þ þ expð2b1 LÞðG 2 b2  G 1 b1 Þ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 o r ro ¯ ¯2 1 b1 ¼ k¯ 2  ; b2 ¼ k¯ 2  2 , mG 1 mG2 In which gx and gt are the Fourier transformations of g^ x and g~ t , b1 and b2 are the complex roots with positive real parts and m is given by Eq. (26). The variance of the non-stationary random ground level displacement response is s2v ðx; tÞ ¼

Z

1

S vv ðx; o; t; 0Þ do. 0

(54)

ð53Þ

Example 2. consists of three layers and a homogeneous semi-infinite space, for which the material properties are listed in Table 1 and were taken from the first four layers of the Gutenberg model [5]. The source is the rupture in the solid and the discontinuous displacements and/or stresses are regarded as non-stationary random processes. The modulation functions of the source were those given by Eq. (51) and the PSD matrix S¯ of the

ARTICLE IN PRESS Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

313

response auto-PSDs at different observation locations are all zero for an initial time interval, which increases for greater distances x from the epicentre to the observation location. Even when the auto-PSDs have reached their maximum values, they still need some time to settle down and this time becomes longer as the distance from the observation location to the epicentre increases, or for lower frequencies, which attenuate more slowly. When the observation location is distant from the earthquake epicentre, e.g. x ¼ 40.0 km, the response auto-PSDs concentrate in the lower frequency domain, which explains why structures with longer natural periods are more easily damaged when they are further from the earthquake epicentre. Fig. 4. Comparison of variances of the analytical formula and PIM. Key: __________ analytical J proposed method.

Table 1 Soil parameters (d ¼ 1, 2, 4, or 6) Layer

l (1010 N/m2)

G (1010 N/m2)

r (103kg/ m3)

Thickness(km)

1 2 3 4

3.3 4.4 8.0 8.2

3.5 4.3 7.2 7.0

2.74 3.00 3.32 3.34

9.5  d 9.5  d 6.0  d Semi-infinite

corresponding stationary random components was assumed to be 2

1:0

0:2  0:3i

0:2 þ 0:2i

3

6 7 0:3  0:1i 7 S11 ¼ 6 4 0:2 þ 0:3i 1:0 5, 0:2  0:2i 0:3 þ 0:1i 1:0 2 3 1:0 0:1  0:2i 0:1 þ 0:1i 6 7 12 0:3  0:1i 7 S22 ¼ 6 4 0:1 þ 0:2i 1:0 5  10 , 0:1  0:1i 0:3 þ 0:1i 1:0 2 3 3:0  2:0i 1:0 þ 2:0i 2:0  3:0i 6 7 5 3:0 þ 4:0i 2:0 þ 2:0i 7 S12 ¼ 6 4 2:0  3:0i 5  10 , 1:0 þ 1:0i 3:0  1:0i 1:0  2:0i " # S11 S12 S¯ ¼ . SH S22 12

ð55Þ

The random ground level displacement response autoPSDs were computed for the four different sets of layer thicknesses given by the alternative values d given in the caption of Table 1, i.e. d ¼ 1, 2, 4, and 6. The source was always at the bottom of the second layer. Figs. 5–7 give the auto-PSDs at observation locations x ¼ 0 or 40 km, of the random ground level displacements along the x, y and z directions, respectively. As for Example 1, the

Example 3. consists of observing the effect of the coherency of the random excitations to the random ground level responses when d ¼ 2 for the system of layers described in Example 2. The source was still at the bottom of the second layer but the PSD matrix S¯ of the stationary random excitations was modified to 2

1:0

6 S11 ¼ 6 4 gð0:2 þ 0:3iÞ gð0:2  0:2iÞ 2 1:0 6 S22 ¼ 6 4 gð0:1 þ 0:2iÞ gð0:1  0:1iÞ 2 3:0  2:0i 6 S12 ¼ g6 4 2:0  3:0i 1:0 þ 1:0i " # S11 S12 ¯S ¼ . SH S22 12

gð0:2  0:3iÞ gð0:2 þ 0:2iÞ 1:0

3

7 gð0:3  0:1iÞ 7 5,

gð0:3 þ 0:1iÞ 1:0

gð0:1  0:2iÞ gð0:1 þ 0:1iÞ 1:0

3

7 12 gð0:3  0:1iÞ 7 5  10 ,

gð0:3 þ 0:1iÞ 1:0

1:0 þ 2:0i 2:0  3:0i

3

7 5 3:0 þ 4:0i 2:0 þ 2:0i 7 5  10 , 3:0  1:0i 1:0  2:0i ð56Þ

Here the coefficient g represents the coherency of the six components of the non-stationary sources. Eq. (56) shows that if g ¼ 0, the PSD matrix S¯ becomes a diagonal matrix, which means that the six components of the non-stationary sources are independent to each other. A larger g means a stronger coherency between the six excitation components. It should be noted that to date, there are only a limited number of seismic records available that can be used to determine the value of g. Therefore, in this paper we assume that g is independent of frequency, from which it follows that the coherency of the six components of the non-stationary sources is also independent of frequency. g ¼ 1 gives Eq. (55), i.e. it corresponds to Example 2. The ground level displacement PSDs were computed for g ¼ 0, 0.2 and 1.0. Fig. 8 gives the auto-PSDs at g ¼ 0 of the displacement along the y direction at observation location x ¼ 0, which are clearly very close to those

ARTICLE IN PRESS 314

Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

Fig. 5. Non-stationary random response auto-PSDs of x direction ground level displacement at x ¼ 0 and 40 km, for d (see Table 1) ¼ 1, 2, 4 and 6. (a) d ¼ 1, x ¼ 0; (b) d ¼ 1, x ¼ 40 km; (c) d ¼ 2, x ¼ 0; (d) d ¼ 2, x ¼ 40 km; (e) d ¼ 4, x ¼ 0; (f) d ¼ 4, x ¼ 40 km; (g) d ¼ 6, x ¼ 0 and (h) d ¼ 6, x ¼ 40 km.

for g ¼ 1.0, which were given above, as Fig. 6(c). The auto-PSDs for g ¼ 0.2 were also very close to those of Fig. 8 and therefore are not shown. The corresponding numerical results for the auto-PSDs of displacement responses along the x and z directions are not presented, but

they also showed that the cross-PSDs of such excitations have insignificant effect on the response auto-PSDs. Fig. 9 gives the norms of the random ground level response cross-PSDs Suv and Suw between the displacement along the x direction and the displacements along the y or z

ARTICLE IN PRESS Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

315

Fig. 6. Non-stationary random response auto-PSDs of y direction ground level displacement at x ¼ 0 and 40 km, for d (see Table 1) ¼ 1, 2, 4 and 6. (a) d ¼ 1, x ¼ 0; (b) d ¼ 1, x ¼ 40 km; (c) d ¼ 2, x ¼ 0; (d) d ¼ 2, x ¼ 40 km; (e) d ¼ 4, x ¼ 0; (f) d ¼ 4, x ¼ 40 km; (g) d ¼ 6, x ¼ 0 and (h) d ¼ 6, x ¼ 40 km.

directions, observed at location x ¼ 0 and for g ¼ 0, 0.2 and 1.0. When g ¼ 0, the cross-PSD Suv is zero, because the displacements along the x and y directions are decoupled for an isotropic solid. However, the cross-PSD Suv between displacements along the x and z directions is non-zero, as

these displacements are coupled, and thus form the P-SV waves. Finally, Fig. 9 shows that the cross-PSDs are dependent upon g to a remarkable extent, such that as the excitation cross-PSDs increase the response cross-PSDs exhibit substantial increases.

ARTICLE IN PRESS 316

Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

Fig. 7. Non-stationary random response auto-PSDs of z direction ground level displacement at x ¼ 0 and 40 km, for d (see Table 1) ¼ 1, 2, 4 and 6. (a) d ¼ 1, x ¼ 0; (b) d ¼ 1, x ¼ 40 km; (c) d ¼ 2, x ¼ 0; (d) d ¼ 2, x ¼ 40 km; (e) d ¼ 4, x ¼ 0; (f) d ¼ 4, x ¼ 40 km; (g) d ¼ 6, x ¼ 0 and (h) d ¼ 6, x ¼ 40 km.

7. Conclusions This paper shows that PEM in combination with the PIM is an effective and highly accurate CQC (complete quadratic combination) [42–49] method for analyzing the

propagation of partially coherent non-stationary random waves in a layered half-space. PEM is used to transform the non-stationary random wave propagation problems into deterministic transient wave response problems, while PIM is used to solve the state equation of wave

ARTICLE IN PRESS Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

motion in the frequency–wavenumber domain. An example with non-stationary random SH waves traveling in a layered half-space, which consists of a single layer and a

Fig. 8. Non-stationary random response auto-PSDs of y direction ground level displacement when g ¼ 0.

317

homogeneous semi-infinite space, has been used to verify the validity and high accuracy of the proposed method. This simple example also shows that the non-stationary effect is more significant when the observation location is horizontally very far from the earthquake epicentre. Non-stationary SH and P-SV waves traveling in a more general multi-layered half-space are investigated in the second and third examples. The results show that for the auto-PSDs and cross-PSDs of the random ground level displacement responses at different observation locations: (1) the non-stationary effect is more evident when the source is deep; (2) at observation locations which are further from the earthquake epicentre, the nonstationary effect is more significant and the response PSDs concentrate in the lower frequency region; (3) the crossPSDs of the source have little effect on the auto-PSDs of the random ground level response, but have a significant effect on the cross-PSDs of the random ground level response.

Fig. 9. Non-stationary random response cross-PSDs of x and y, or x and z, direction ground level displacements at x ¼ 0, for g ¼ 0, 0.2, 1.0. (a) g ¼ 0; (b) g ¼ 0; (c) g ¼ 0.2; (d) g ¼ 0.2; (e) g ¼ 1.0 and (f) g ¼ 1.0.

ARTICLE IN PRESS Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

318

2

Acknowledgements The authors are grateful for the support of: the doctoral research fund of the Chinese Ministry of Education (No. 20040141020); the Natural Science Foundation of China (No. 10472023) and; the Cardiff Advanced Chinese Engineering Centre.

i

6 6 61 6 T2 ¼ 6 6 ¯ 6 i2kG 4 ¯ 2kG

¯

¯

lþ3G i 2kð ¯ GÞ ¯ lþ

¯

¯

i

lþ3G i 2kð ¯ GÞ ¯ lþ

¯ G ¯ lþ3  2kð ¯ GÞ ¯ lþ

1

¯ G ¯ lþ3 ¯ GÞ ¯ 2kðlþ

¯ iG

¯ i2kG

¯ iG

¯ G

¯ 2kG

¯ G

3 7 7 7 7 7. 7 7 5

ðA:6Þ

(3) When o6¼0 ¯j Appendix A. The eigenvector and eigenvalue matrices of H ¯ are the Lame´ constants and r¯ is In this section, l¯ and G the density of the semi-infinite layer l+1. For the SH wave, cases (1) and (2) below must be considered, with " # 1 0 ¯ mG ¯1 ¼ H . (A.1) ¯  ro k2 mG ¯ 2 0

1 T1 ¼ ¯ aG

K1 ¼ diagða; aÞ;

 1 ; ¯ aG

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ro ¯ 2 , a ¼ k2  ¯ mG

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ro ¯ 2 , a2 ¼ k 2  ¯ mG

K2 ¼ diagða1 ; a2 ; a1 ; a2 Þ 2 ik

a1

6 6 61 6 T2 ¼ 6 6 i2kmG ¯ 6 6 4 2 ¯

¯ 2Þ ð2k mGro a1

¯  ro (1) When k2 mG ¯ 2 a0 

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ro ¯ 2 a1 ¼ k2  ; ¯ mðl¯ þ 2GÞ 1

 aik

1

 aik

1

ik a2

¯ ro ¯ 2Þ ð2k2 mG a2

¯ i2kmG

¯ i2kmG

 ð2k

1

2

2 mG ¯ ro ¯ 2Þ

a1

3

7 7 7 7 7 2 2 ¯ ¯ Þ 7, ð2k mGro  7 a2 7 5 ¯ i2kmG

ðA:7Þ in which a1 and a2 take the roots with the positive real part.

(A.2) in which a takes the root with the positive real part. ¯  ro (2) When k2 mG ¯ 2¼0     0 0 0 1 K1 ¼ ; T1 ¼ . (A.3) ¯ 0 1 0 G For the P-SV wave, cases (1)–(3) below must be considered, with 3 2 1 0 ik 0 ¯ mG 7 6 1 7 6 ik ¯ l¯ 0 0 ¯ ¯ ¯ 6 lþ2G mðlþ2GÞ 7 7. H2 ¼ 6   2 7 6 2 4l¯ Gþ4 l¯ 7 2 6 k m ¯¯ ¯G¯  ro 0 0 ik ¯ ¯ G ¯ 5 4 lþ2 lþ2G 0 ro ¯ 2 ik 0

(1) When k ¼ 0 and o ¼ 0 2 2 3 0 0 0 0 0 60 60 0 0 07 6 6 7 K2 ¼ 6 7; T2 ¼ 6 ¯ 4G 41 0 0 05 0

1 0

0

0

(2) When ka0 and o ¼ 0 3 2 k 0 0 0 7 6 61 k 0 07 7 6 K2 ¼ 6 7, 6 0 0 k 0 7 5 4 0

0

1

k

1 0

0 0 ¯ l¯ þ 2G

0 17 7 7. 0 05

The wave motion equations for the finite viscoelatic layer are   2  q q v~1 ðx; z; tÞ q2 v~1 ðx; z; tÞ G 1 q0 þ q1 þ qt qx2 qz2   2 q q v~1 ðx; z; tÞ , ¼ r1 p 0 þ p 1 qt qt2     q q q~v1 ðx; z; tÞ , p0 þ p1 t~ 1yz ðx; z; tÞ ¼ G 1 q0 þ q1 qt qt qz ðB:1Þ

(A.4)

and the wave motion equations for the semi-infinite space are

(A.5)

  2  q q v~2 ðx; z; tÞ q2 v~2 ðx; z; tÞ G 2 q0 þ q1 þ qt qx2 qz2   2 q q v~2 ðx; z; tÞ ¼ r2 p 0 þ p 1 , qt qt2     q q q~v2 ðx; z; tÞ , p0 þ p1 t~ 2yz ðx; z; tÞ ¼ G 2 q0 þ q1 qt qt qz

3

0

Appendix B. Derivation of analytical solution for Example 1

0 0

ðB:2Þ If the PSD matrix for the stationary components of v¯ and t¯ yz is " # S v¯ v¯ S v¯ t¯ yz ¯ SðoÞ ¼ (B.3) S t¯ yz v¯ S t¯ yz t¯ yz

ARTICLE IN PRESS Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

319

infinite space, so C3 in Eq. (B.12) must be zero which gives

and if the following definitions are used: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S v¯ t¯ yz S t¯ yz v¯ S v¯ t¯ yz , s1;2 ¼ S v¯ t¯ yz S t¯ yz v¯ þ ðl1;2  Sv¯ v¯ Þ2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S v¯ t¯ yz S t¯ yz v¯ ðl1;2  S v¯ v¯ Þ, t1;2 ¼ Sv¯ t¯ yz S t¯ yz v¯ þ ðl1;2  S v¯ v¯ Þ2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðSv¯ v¯ þ St¯ yz t¯ yz Þ  ðS v¯ v¯  S t¯ yz t¯ yz Þ2 þ 4Sv¯ t¯ yz S t¯ yz v¯ , l1;2 ¼ 2

v2 ðk; ¯ oÞ ¼ C 2 expðb2 zÞ. ¯ z; o;

(B.13)

Substituting Eqs. (B.11) and (B.13) into Eq. (B.10) gives " # ½expðb1 LÞ þ expðb1 LÞ expðb2 LÞ

ðB:4Þ then the pseudo excitations at z ¼ L are

mG1 b1 ½expðb1 zÞ  expðb1 zÞ mG 2 b2 expðb2 zÞ ( ) ( ) s1;2 C1  ¼ gx ðkÞg . ðB:14Þ ¯  oÞ ¯ t ðo t1;2 C2 The solution of Eq. (B.14) is

v~2 ðx; L; tÞ ¼ v~1 ðx; L; tÞ þ g^ x ðxÞg~ t ðtÞ expðiotÞs1;2 , t~ 2yz ðx; L; tÞ ¼ t~ 1yz ðx; L; tÞ þ g^ x ðxÞg~ t ðtÞ expðiotÞt1;2 .

ðB:5Þ

C1 ¼

gx ðkÞg ¯  oÞðmG 2 b2 s1;2 þ t1;2 Þ expðb1 LÞ ¯ t ðo . m½ðG 2 b2 þ G1 b1 Þ þ expð2b1 LÞðG 2 b2  G 1 b1 Þ (B.15)

The free boundary condition at z ¼ 0 is t~ 1yz ðx; 0; tÞ ¼ 0

(B.6)

By using Fourier transformation, the wave motion equations, boundary conditions and pseudo excitations can be transformed into the frequency–wavenumber domain as   ro ¯2 v001 ðk; (B.7) v1 ðk; ¯ oÞ ¼ k¯ 2  1 ¯ oÞ ¯ z; o; ¯ z; o; mG 1   r2 o ¯2 00 2 v2 ðk; v2 ðk; ¯ oÞ ¼ k¯  ¯ oÞ, ¯ z; o; ¯ z; o; mG 2

(B.8)

t1yz ðk; ¯ oÞ ¼ 0, ¯ z; o;

(B.9)

v2 ðk; ¯ oÞ ¼ v1 ðk; ¯ oÞ þ gx ðkÞg ¯  oÞs1;2 , ¯ L; o; ¯ L; o; ¯ t ðo mG 2 v02 ðk; ¯ oÞ ¼ mG 1 v01 ðk; ¯ oÞ þ gx ðkÞg ¯  oÞt1;2 . ¯ L; o; ¯ L; o; ¯ t ðo ðB:10Þ The solution to Eq. (B.7), subject to the free boundary condition (B.9), is sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r o ¯2 v1 ðk; ¯ oÞ ¼ C 1 ½expðb1 zÞ þ expðb1 zÞ; b1 ¼ k¯ 2  1 ¯ z; o; mG 1 (B.11) and the solution to Eq. (B.8) is v2 ðk; ¯ oÞ ¼ C 2 expðb2 zÞ þ C 3 expðb2 zÞ; ¯ z; o;

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r o ¯2 b2 ¼ k¯ 2  2 . mG 2

(B.12) Here b1 and b2 are complex and the root taken is the one with a positive real part. The radiative boundary condition requires that there are no up-going waves in the semi-

Then the pseudo response at ground level in the frequency–wavenumber domain is v1 ðk; ¯ oÞ ¼ 2C 1 ¯ 0; o;

(B.16)

from which the pseudo response at ground level in the time-space domain is given by Eq. (53). References [1] Timoshenko SP, Goodier JN. Theory of elasticity. New York: McGraw-Hill; 1951. [2] Ewing WM, Jardetzky WS, Press F. Elastic waves in layered media. New York: McGraw-Hill; 1957. [3] Achenback JD. Wave propagation in elastic solids. Amsterdam: North-Holland; 1973. [4] Graff KF. Wave motion in elastic solids. Oxford: Clarendon Press; 1975. [5] Aki K, Richards PG. Quantitative seismology. San Francisco: Freeman; 1980. [6] Brekhovskikh LM. Waves in layered media. New York: Academic Press; 1980. [7] Kennett BLN. Seismic wave propagation in stratified media. Cambridge: Cambridge University Press; 1983. [8] Caviglia G, Morro A. Riccati equations for wave propagation in planarly stratified solids. Eur J Mech A/Solids 2000;19(4):721–41. [9] Caviglia G, Morro A. Wave propagation in multilayered anisotropic solids. Int J Eng Sci 2000;38(8):847–63. [10] Caviglia G, Morro A. Reflection and transmission in anisotropic dissipative multilayers. Eur J Mech A/Solids 2002;21(6):1055–67. [11] Verma KL. On the propagation of waves in layered anisotropic media in generalized thermoelasticity. Int J Eng Sci 2002;40(18): 2077–96. [12] Gulyayev VI, Lugovyy PZ, Ivanchenko GM. Discontinuous wave fronts propagation in anisotropic layered media. Int J Solids Struct 2003;40(1):237–47. [13] Gao Q, Zhong WX, Howson WP. A precise method for solving wave propagation problems in layered anisotropic media. Wave Motion 2004;40(3):191–207. [14] Alshaikh IABU, Turhan D, Mengi Y. Two-dimensional transient wave propagation in viscoelastic layered media. J Sound Vibr 2001; 244(5):837–58. [15] Zhang JF, Li YM. Numerical simulation of elastic wave propagation in inhomogeneous media. Wave Motion 1997;25(2):109–25.

ARTICLE IN PRESS 320

Q. Gao et al. / Soil Dynamics and Earthquake Engineering 28 (2008) 305–320

[16] Biot MA. Theory of propagation of elastic waves in a fluid-saturated porous solid. I: low frequency range. J Acoust Soc Am 1956;28(2): 168–78. [17] Biot MA. Theory of propagation of elastic waves in a fluid-saturated porous solid. II: higher-frequency range. J Acoust Soc Am 1956; 28(2):179–91. [18] Degrande G, Roeck GD, Broeck PVD, Smeulders D. Wave propagation in layered dry, saturated and unsaturated poroelastic media. Int J Solids Struct 1998;35(34-35):4753–78. [19] Yang J. Importance of flow condition on seismic waves at a saturated porous solid boundary. J Sound Vibr 1999;221(3):391–413. [20] Yang J. Influence of water saturation on horizontal and vertical motion at a porous soil interface induced by incident P wave. Soil Dyn Earthquake Eng 2000;19(8):575–81. [21] Nelson JT. Prediction of ground vibration from trains using seismic reflectivity methods for a porous soil. J Sound Vibr 2000;231(3):727–37. [22] Li XK, Zhang JB, Zhang HW. Instability of wave propagation in saturated poroelastoplastic media. Int J Numer Analyt Methods Geomech 2002;26(6):563–78. [23] Zhou XL, Wang JH, Lu JF. Transient foundation solution of saturated soil to impulsive concentrated loading. Soil Dyn Earthquake Eng 2002;22(4):273–81. [24] Vashishth AK, Khurana P. Inhomogeneous waves in anisotropic porous layer overlying solid bedrock. J Sound Vibr 2002;258(4):577–94. [25] Mamolis GD. Stochastic soil dynamics. Soil Dyn Earthquake Eng 2002;22(1):3–15. [26] Fishman KL, Ahmad S. Seismic response for alluvial valleys subjected to SH-waves, P-waves and SV-waves. Soil Dyn Earthquake Eng 1995;14(4):249–58. [27] Vai R, Castillo-Covarrubias JM, Sanchez-Sesma FJ, Komatitsch D, Vilotte JP. Elastic wave propagation in an irregularly layered medium. Soil Dyn Earthquake Eng 1999;18(1):11–8. [28] Heymsfield E. Two-dimensional scattering of SH waves in a soil layer underlain with sloping bedrock. Soil Dyn Earthquake Eng 2000;19(7):489–500. [29] Kiyono J, Toki K, Sato T, Mizutani H. Simulation of stochastic waves in a non-homogenous random field. Soil Dyn Earthquake Eng 1995;14(5):387–96. [30] Manolis GD, Show RP. Harmonic wave propagation through viscoelastic heterogeneous media exhibiting mild stochasticity: I. Fundamental solution. Soil Dyn Earthquake Eng 1996;15(2):119–27. [31] Manolis GD, Show RP. Harmonic wave propagation through viscoelastic heterogeneous media exhibiting mild stochasticity: II. Applications. Soil Dyn Earthquake Eng 1996;15(2):129–39. [32] Parra JO, Ababou R, Sablik MJ, Hackert CL. A stochastic wavefield solution of the acoustic wave equation based on the random Fourier–Stieltjes increments. J Appl Geophys 1999;42(2):81–97. [33] Parra JO, Zook BJ. Stochastic wave field solution of the 2D elastic wave equation based on the random Fourier–Stieltjes increments. J Appl Geophys 2001;48(1):43–63. [34] Wang S, Hong H. Effects of random variations of soil properties on site amplification of seismic ground motions. Soil Dyn Earthquake Eng 2002;22(7):551–64.

[35] Zhong WX, Lin JH, Gao Q. The precise computation for wave propagation in stratified materials. Int J Numer Methods Eng 2004;60(1):11–25. [36] Gao Q, Lin JH. Stationary random waves propagation in 3D viscoelastic stratified solid. Appl Math Mech 2005;26(6):785–96. [37] Priestley MB. Power spectral analysis of nonstationary random processes. J Sound Vibr 1967;6(1):86–97. [38] Li QS, Zhang YH, Wu JR, Lin JH. Seismic random vibration analysis of tall buildings. Eng Struct 2004;26(12):1767–78. [39] Fan LC, Wang JJ, Chen W. Response characteristics of long-span cable-stayed bridges under non-uniform seismic action. Chinese J Comput Mech 2001;18(3):358–63. [40] Liu TY, Liu GT. Seismic random response analysis of arch dams. Eng Mech 2000;17(6):20–5. [41] Zhao YQ, Zhang GY, Guo KH. Handling safety simulation of driver-vehicle closed-loop system with evolutionary random road input. Veh System Dyn 2000;33(3):169–81. [42] Lin JH. A fast CQC algorithm of PSD matrices for random seismic responses. Comput Struct 1992;44(3):683–7. [43] Lin JH, Williams FW, Zhang WS. A new approach to multiphaseexcitation stochastic seismic response. Microcomput Civil Eng 1993;8(4):283–90. [44] Lin JH, Zhang WS, Williams FW. Pseudo-excitation algorithm for nonstationary random seismic responses. Eng Struct 1994;16(4): 270–6. [45] Lin JH, Shen WP, Williams FW. A high precision direct integration scheme for non-stationary random seismic responses of nonclassically damped structures. Struct Eng Mech 1995;3(3):215–28. [46] Lin JH, Song GZ, Sun Y, Williams FW. Non-stationary random seismic response of non-uniform beams. Soil Dyn Earthquake Eng 1995;14(4):301–6. [47] Lin JH, Fan Y, Bennett PN, Williams FW. Propagation of stationary random waves along substructural chains. J Sound Vibr 1995;180(5): 757–67. [48] Lin JH, Fan Y, Williams FW. Propagation of non-stationary random waves along substructural chains. J Sound Vibr 1995;187(4): 585–93. [49] Lin JH, Wang XC, Williams FW. Asynchronous parallel computing of structural non-stationary random seismic responses. Int J Numer Methods Eng 1997;40(12):2133–48. [50] Zhong WX. Duality system in applied mechanics and optimal control. Boston: Kluwer Academic Pub.; 2004. [51] Zhong WX. On precise integration method. J Comput Appl Math 2004;163(1):59–78. [52] Shames IH, Cozzarelli FA. Elastic and inelastic stress analysis. New Jersey: Prentice-Hall; 1992. [53] Preumont A. Random vibration and spectral analysis. London: Kluwer; 1994. [54] Lin JH, Zhang YH. Pseudo excitation method of random vibration. Beijing: Science Press; 2004 [In Chinese]. [55] Lin J, Zhang YH. Seismic random vibration of long-span structures. In: de Silva CW, editor. Vibration and shock handbook. New York: CRC Press; 2005 [Chapter 30].