A high-order time-domain transmitting boundary for cylindrical wave propagation problems in unbounded saturated poroelastic media

Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

Contents lists available at SciVerse ScienceDirect

Soil Dynamics and Earthquake Engineering journal homepage: www.elsevier.com/locate/soildyn

A high-order time-domain transmitting boundary for cylindrical wave propagation problems in unbounded saturated poroelastic media Peng Li, Er-xiang Song n Department of Civil Engineering, Tsinghua University, Beijing 100084, China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 16 August 2012 Accepted 16 January 2013 Available online 1 March 2013

Based on the u–p formulation of Biot equation with an assumption of zero permeability coefficient, a high-order transmitting boundary is derived for cylindrical elastic wave propagation in infinite saturated porous media. By this transmitting boundary the total stresses on the truncated boundaries of a numerical model, such as a finite element model, are replaced by a set of spring, dashpot and mass elements, with some additionally introduced auxiliary degrees of freedom. The transmitting boundaries are incorporated into the DIANA SWANDYNE II program and an unconditionally stable implicit time integration algorithm is adopted. Despite the assumption made in the derivation of the transmitting boundary, numerical examples show that it can provide highly accurate results for cylindrical elastic wave propagation problems in infinite saturated porous medium in case the u–p formulation is applicable. Although the direct applications of the proposed transmitting boundary to general two dimensional wave problems in infinite saturated porous media are not highly accurate, acceptable accuracy can still be achieved by placing the transmitting boundary at relatively large distance from the wave source. & 2013 Elsevier Ltd. All rights reserved.

1. Introduction The transient wave propagation in elastic porous media which characterizes the behavior of saturated soil occurs in many fields of engineering and physics, such as earthquake engineering and soil dynamic analysis, where the coupled interaction of the soil skeleton and the pore fluid should be considered. The theory of elastic wave propagation in saturated porous media was first established by Biot [1,2]. Since then, many researchers [3–8] have studied this theory and developed relevant finite element models. The finite element method is widely used in solving wave propagation problems for its versatility and reliability. However, the numerical simulation is directly applicable only to a finite domain. When the spatial domain of the problem is unbounded, a proper artificial boundary designed to avoid energy reflection of outgoing waves is to be imposed on the restraining boundary. Such artificial boundary condition is also named in many literatures as absorbing, silent, non-reflecting, transmitting, radiating, transparent boundary conditions. The term ‘transmitting boundary’ will be used here. For one-phase media, a large number of transmitting boundaries have been developed, and great advances have been achieved

n

Corresponding author. Tel.: þ86 10 62773545; fax: þ86 10 62771132. E-mail address: [email protected] (E.-x. Song).

0267-7261/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.soildyn.2013.01.006

[9–13]. Along these works, the transmitting boundaries can be divided clearly into two categories: the global exact boundary and the local approximate boundary [14]. The former category is global in time and space, i.e., the current response at a point on transmitting boundary is related to the responses at all transmitting boundary points during all previous time steps. Several representative methods such as the boundary element method [15], the thin-layer method [16], the boundary integral methods based on Kirchhoff integral [17,18] and the Dirichlet-to-Neumann (DtN) boundary [19,20] belong to this category. The global exact boundary can give stable solution due to its exactness even imposed quite close to the wave source, but its ‘non-local in space’ leads to large amount of calculation and low efficiency. Moreover, such boundary is formulated in the frequency domain so it cannot be used directly with time-step integration techniques in solving nonlinear problems. The latter category is local in time and space, i.e., the current response at a point on the transmitting boundary is related to the responses at a few adjacent points during a few previous time steps. Several wellknown transmitting boundaries such as the viscous boundary [21,22], the viscous-spring boundary [23–25], the superposition boundary [26,27], the extrapolation boundary [28,29], the paraxial boundary [30,31] and the multi-direction boundary [32,33] all belong to this category. These transmitting boundaries are derived from the various symmetric cases of wave radiations which are spatially one-dimensional, such as plane, cylindrical

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

and spherical waves, therefore, these transmitting boundaries are local in space. Compared with the viscous boundary which is constructed based on a travelling wave of constant magnitude, the viscous-spring boundary has a higher precision because its assumption of scattered wave field takes the geometric attenuation [23] into account. Therefore, it can simulate not only the radiation damping of the infinite media but also the elasticity recovery capacity of the far field media outside the boundary with a good stability both under high frequency and low frequency seismic loads [34]. In addition, the infinite element method [35–38] also belongs to this category. The ‘decoupling in time and space’ makes the local approximate boundary easily implemented and can be used directly with time-step integration techniques in solving nonlinear problems with a low computational costs and storage. However, such boundary should be located sufficiently away from the wave source in order to obtain a solution of acceptable accuracy. In the last two decades the high-order accurate transmitting boundary has been developed, which is the more powerful method for wave propagation problem in unbounded media. For the high-order local boundaries such as higher order paraxial boundary [31], some specially defined auxiliary variables are introduced to develop the high-order accurate transmitting boundary that can be implemented up to an arbitrarily high order. A review on such high-order accurate transmitting boundaries can be found in [13] and the references therein. For the exact transmitting boundaries that are derived from analytical solutions and are always global in time, a temporal localization method is used, which consists of the convolution-to-differentiation treatment and the auxiliary variable realization [14]. For example, to construct the viscous–spring transmitting boundary, Deeks and Randolph [23] used the approximation expression for the wave front of a traveling cylindrical wave, resulting in a low-order approximation boundary for cylindrical wave radiation problem. Based on the exact solution of cylindrical elastic wave radiation problem, the dynamic-stiffness coefficient in frequency domain was usually approximated by rational function (Pade) approximation, which is an effective method due to its large convergence range and high convergence rate as the order increases [39]. Then through partial fraction expansion, Wolf [40–42] and Wu and Lee [43,44] established high-order accurate transmitting boundaries in time domain as (consistent or systematic, respectively) lumped-parameter models, which contain more spring, dashpot and mass elements and auxiliary degrees of freedom. Recently Du and Zhao [45] proposed a high-order accurate local transmitting boundary for cylindrical wave propagation problem, which is dynamically stable and easily implemented into the commercial finite element codes. For two-phase (saturated porous) media, the interaction between the two phases, i.e., the elastic skeleton and the compressible pore fluid, makes it more complicated to construct the transmitting boundary for dynamic analysis. It is well known from Biot theory that there are two dilatational waves (P1 and P2) and one shear wave in saturated porous media. Inertial and mechanical coupling of the two phases makes the wave speeds not equal to that of single phase [46], whereas viscous coupling makes wave propagation dispersive. In condition of low permeability and low frequency, the relative motion of the pore fluid with respect to the solid skeleton disappears. The second dilatational wave is completely damped out by viscous coupling, thus, the saturated porous media behaves as a one-phase media and has only one dilatational wave and shear wave. In contrast, in condition of high permeability and high frequency, no viscous coupling exists in saturated porous media and there are two dilatational waves and one shear wave. The body wave speeds of saturated porous media in these two extreme conditions are

49

independent of frequency. For intermediate value of permeability and frequency, the body wave speeds are dispersive. Dilatational wave with higher frequency has higher wave speed and attenuates more rapidly [47]. The transition from low to high viscous coupling occurred over quite a narrow range of permeability or frequency values for a given distance [48]. Since the 1990s some studies have been carried out on the transmitting boundary of saturated porous media based on the Biot theory [1] and different forms of governing equations. Modaressi and Benzenati [49,50] first extended the paraxial approximation approach to the case of two-phase media using u–p formulation. By using Fourier transforms to solve the equations, three body wave speeds have been derived. Then Akiyoshi et al. [51,52] proposed paraxial approaches with u–U, u–w and u–p formulations, respectively, for linear isotropic, transverse isotropic and anisotropic porous media. It is noticed that the zeroth-order paraxial boundary is equivalent to a viscous boundary. In all the studies above the wave equation is assumed to be harmonic. To obtain the intensities of the dampers the second dilatational wave was neglected. By using Fourier transforms to solve the u–w form dynamic equations of plane wave motion, Degrande and De Roeck [53,54] have decomposed the solution of the wave field into incident and reflected waves and zeroed the amplitude of reflected wave in order to establish a transmitting boundary condition which is frequency dependent. Only in the low-frequency limit, it reduces to a standard viscous boundary. Based on one-dimensional hypothesis, Gajo et al. [55] established a first-order differential equation system that satisfied the condition of outgoing wave, then deduced the matrix defining the viscous boundary both for the two extreme conditions of zero permeability and infinite permeability. The first condition corresponds to high viscous coupling between the solid skeleton and the pore fluid, and the results that can be called the ‘undrained’ boundaries are similar to those obtained by Modaressi and Akiyoshi. The second condition corresponds to no viscous coupling between solid skeleton and pore fluid and both the two dilatational waves are taken into account. In this condition the results can be called the ‘drained’ boundaries. For intermediate viscous coupling, numerical examples showed that the error made using the solution for infinite permeability is limited to 1% at the most. Assuming infinite permeability, Zerfa and Loret [56] derived the body wave speeds of the elastic and isotropic porous media based on the constitutive equations proposed by Bowen’s formulation. Then viscous transmitting boundary was developed with the u–w formulation. This viscous boundary regards the existence of the second dilatational wave and works correctly for lower permeability. The above viscous boundaries or paraxial boundaries in particular circumstances are frequency independent and can be used for transient analysis in the time domain. Similar as that in one-phase media, the viscous boundary for saturated porous media also has the stability problem when dealing with low frequency seismic loads. Based on cylindrical elastic wave radiation problem in saturated porous media, Liu and Song [57] proposed a viscous-spring transmitting boundary assuming zero permeability, the derivation method of which is the same as Deeks and Randolph’s method [23] in establishing the viscousspring transmitting boundary in one phase media. Meanwhile, Wang and Zhao [58] developed a viscous-spring transmitting boundary assuming infinite permeability. The second dilatational wave has been neglected in both of the above mentioned studies. In this paper, a high-order accurate local time-domain transmitting boundary for simulating the transient scalar wave propagation in plane-strain unbounded saturated porous media is presented from the cylindrical elastic wave radiation problem. Governing equations for the elastic and isotropic saturated porous media follow Zienkiewicz’s u–p formulation with an assumption

50

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

of zero permeability. The derivation starts with the construction of the dynamic stiffness coefficients for decoupled cylindrical P and SH waves which are temporally global. Then, a temporal localization method is used, which consists of the Pade approximation and the auxiliary variable realization. Finally, the spring, dashpot and mass parameters of the high-order local timedomain transmitting boundary are determined. The approach is implemented into the DIANA SWANDYNE II program which is a transient dynamic finite element analysis code for saturated porous media based on u–p formulation, and an unconditionally stable implicit time integration algorithm is adopted in solving transient scalar wave propagation problems in the system consisting of saturated soil near-field and high-order accurate local transmitting boundary. Numerical examples are given to demonstrate the accuracy of the proposed transmitting boundary in solving cylindrical elastic wave radiation problem in saturated porous media, and its potential applications to general twodimensional infinite wave problems are also discussed.

sij ¼ s0ij aDij p ¼ 2meij þ lekk Dij aDij p ½kf ðp,i þ rf bi Þ,i þ ae_ ii þ

1 p_ ¼ 0 Q

ð7Þ

ð8Þ

This so-called u–p formulation was introduced by Zienkiewicz [8] and is valid for most problems of earthquake analysis and frequencies lower than this. Yang [59] has given quantitative demonstration on the applicability of u–p formulation by comparing the body wave speeds derived from u–p formulation and u–w formulation referring to the dimensionless abscissa f/fc, where f is the frequency and fc is the characteristic frequency [3] that can be considered as a frequency scale of the material. It is noticed that when logðf =f c Þ o 0:5

ð9Þ

the wave speeds calculated from both formulations are almost the same. This restriction can be treated as the applicable condition of u–p equation. Here the u–p formulation is adopted to carry out the following solution procedures.

2. Governing equations In the 1950s, the governing equations of motion of saturated linear elastic porous media were first given by Biot [1,2]. From then on many subsequent derivations and modifications have been made. The generally used forms of equations for dynamic analysis of saturated soils are proposed by Zienkiewicz [6]. The first equation is overall equilibrium relation for the soil– fluid mixture:

sij,j þ rbi ¼ ru€ i þ rf w€ i

ð1Þ

The second equilibrium equation ensures momentum balance of the fluid:   1 € i þ k1 _ ð2Þ p,i þ rf bi ¼ rf u€ i þ w f wi n

3. Dynamic stiffness coefficients The cylindrical elastic wave is a particular form of threedimensional elastic wave, and can be decoupled into three independent cylindrical waves which propagate along the radial coordinate, i.e., the cylindrical P wave with the radial displacement ur (r, t), the cylindrical SV wave with the circumferential displacement uy(r, t) and the cylindrical SH wave with the axial displacement uz (r, t). Here, only P wave and SV wave equations are considered in order to derive the transmitting boundary for plane-strain wave problems in two dimensions. For cylindrical SV wave, momentum equilibrium equation of the micro body and shear stress strain relationship are described as @2 uy @ðtry  r 2 Þ ¼ @r @t 2

ð10Þ

  @uy uy  @r r

ð11Þ

The third equation is the effective stress and constitutive relationship (elastic) of saturated soils:

rr2

sij ¼ s0ij aDij p ¼ 2meij þ lekk Dij aDij p

try ¼ G

ð3Þ

The last equation is the one accounting for mass balance of the flow: wi,i ¼ aei,i þp=Q

ð4Þ

The above equations form the u–w–p formulation. Where ui and wi are the solid displacement and the pore fluid relative displacement to solid, respectively; sij is the total stress; b is the body force; p is the pore fluid pressure; n is the porosity of the soil; r, rs and rf are densities of the assembly, the solid and fluid phases, respectively, and by their definitions r ¼(1 n)rs þnrf; kf is the permeability which is related with Darcy permeability coefficient k by kf ¼krfg in which g is the gravity acceleration; l and m are Lame constants for the solid skeleton; a and Q are two introduced parameters, which can be expressed as

a ¼ 1K b =K s 1=Q ¼ n=K f þðanÞ=K s

ð5Þ

where Ks, Kf and Kb are the bulk modulus of the solid, fluid and the assembly, respectively. For medium speed phenomena it is reasonable to assume that the terms relating to the fluid acceleration can be neglected, keeping the absolute displacements of the solid skeleton and the pore water pressure as the primary unknowns (u–p formulation). The governing equations in this case are given as

sij,j þ rbi ru€ i ¼ 0

ð6Þ

where try is circumferential shear stress; G (equal to m) is shear modulus of the solid skeleton. By substituting Eq. (11) into Eq. (10), the wave equation for cylindrical SV wave propagation may be obtained as @2 uy 1 @uy 1 1 @2 u  2 uy ¼ 2 2y þ ð12Þ 2 r @r @r r V s @t pffiffiffiffiffiffiffiffiffi where V s ¼ G=r is the shear wave speed. The shear wave calculated from u–p formulation is not dispersive. For usual values of permeability, two dilatational waves exist and the wave propagation is dispersive. Since the transmitting boundary conditions are expected to be local in time, assumption of zero permeability coefficient is made here which is the same as Akiyoshi [51] and Modaressi [50]. In this extreme case no relative motion exists between the solid skeleton and the pore fluid and the ‘frozen mixture’ behaves like a one-phase media. The P2 wave vanishes and the P1 wave is not dispersive in this condition. It is known that for low value of permeability and frequency content, in which condition u–p formulation is valid, the variation of P1 wave speed is quite small and P2 wave is quickly damped out [52]. Thus, the simplification relating to zero permeability coefficient may be considered acceptable with the restriction of u–p formulation. For cylindrical P wave, the governing Eqs. (11) and (12) in which the term of body force is neglected, are written in

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

cylindrical coordinate system as 2

@ ur @sr 1 þ ðsr sy Þ ¼ r @r @t 2   ( sr ¼ s0r ap ¼ l þ2m @u@rr þ l urr ap   sy ¼ s0y ap ¼ l @u@rr þ l þ2m urr ap

r

ð13Þ

ð14Þ

With the permeability coefficient taken as zero, Eq. (8) is expressed as   @ur ur þ ð15Þ p ¼ aQ @r r Eq. (15) implies that the pore fluid pressure is proportional to the volume strain of the solid skeleton when permeability coefficient is zero. By substituting Eqs. (14) and (15) into Eq. (13), the wave equation for cylindrical P wave propagation may be obtained as @2 ur 1 @ur 1 1 @ 2 ur  ur ¼ 2 2 þ ð16Þ r @r r 2 @r 2 V p @t pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where V p ¼ ðl þ2G þ a2 Q Þ=r is the dilatational wave speed. Eqs. (12) and (16) represent the cylindrical elastic wave radiation problem of saturated porous media under plain strain conditions. The analytical solutions for these equations in frequency domain may be generated in a similar manner as presented by Novak [60] and Novak and Mitwally [61]. Using Fourier transform, the cylindrical SV wave Eq. (12) becomes ! @2 1@ 1 o2  uy ðr, oÞ ¼ 0 þ þ ð17Þ r @r r 2 @r 2 V 2S The outgoing wave solution of Eq. (17) can be written in form of the second kind of Hankel function of order one   or ð18Þ uy ðr, oÞ ¼ Hð2Þ 1 Vs Eq. (11) implies the circumferential interaction between the near field and the truncated far field infinite domain. By substituting Eq. (18) into the Fourier transform of Eq. (11), the stress– displacement relationship on the artificial boundary at r ¼R can be exactly described in frequency domain as f y ðR, oÞ ¼ Sy ðR, oÞuy ðR, oÞ

ð19Þ

where Sy(R, o) is circumferential dynamic stiffness coefficient which can be expressed as " # os Hð2Þ 0 ðos Þ Sy ðR, oÞ ¼ S0y 1 ð20Þ 2 Hð2Þ ðos Þ 1

and its standard form is   Sy R, oÞ ¼ S0y ky ðos Þ þ ios cy ðos Þ

ð21Þ

where os ¼ oR/Vs is the dimensionless S wave frequency; ð2Þ S0y ¼ 2G=R is circumferential static stiffness; Hð2Þ 0 and H1 are the second kind of Hankel functions of order zero and one, respectively; ky(os) and cy(os) are dimensionless circumferential spring and damping coefficients that are functions of the dimensionless frequency os, and can be expressed in terms of Bessel functions of the first and second kind. ky ðos Þ ¼ 1

c y ð os Þ ¼

1

os J0 ðos ÞJ1 ðos Þ þY 0 ðos ÞY 1 ðos Þ 2 J 21 ðos Þ þ Y 21 ðos Þ 1

pos J21 ðos Þ þY 21 ðos Þ

In a similar way, by using Fourier transform, the cylindrical P wave Eq. (16) becomes ! @2 1@ 1 o2  þ þ ð24Þ ur ðr, oÞ ¼ 0 r @r r 2 @r 2 V 2p with the outgoing wave solution written as   or ur ðr, oÞ ¼ Hð2Þ 1 Vp

ð25Þ

By substituting Eq. (15) into Eq. (14), the radial interaction between the near field and the truncated far field infinite domain is     @ur ur 2m ur  l þ2m þ a2 Q þ ð26Þ f r ¼ sr ðr, oÞ ¼ r @r r By substituting Eq. (25) into the Fourier transform of Eq. (26), the radial stress and displacement at the internal radius r ¼R is related by f r ðR, oÞ ¼ Sr ðR, oÞur ðR, oÞ

ð27Þ

where Sr(R, o) is radial dynamic stiffness coefficient with the expression: " # Z2 op Hð2Þ 0 0 ð op Þ Sr ðR, oÞ ¼ Sr 1 ð28Þ 2 Hð2Þ ðop Þ 1

and its standard form is   Sr R, oÞ ¼ S0r kr ðop Þ þiop cr ðop Þ

ð29Þ

where op ¼ oR/Vp is the dimensionless P wave frequency; S0r ¼ 2G=R is the radial static stiffness; Z is the ratio of P wave speed to S wave speed: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Vp l þ2m þ a2 Q Z¼ ¼ ð30Þ Vs m where kr(op) and cr(op) are, respectively, the dimensionless radial spring and damping coefficients that are functions of the dimensionless frequency op, and can be expressed in terms of Bessel functions of the first and second kind.   Z2 op J0 ðop ÞJ1 ðop Þ þ Y 0 ðop ÞY 1 ðop Þ kr op ¼ 1 2 J 21 ðop Þ þ Y 21 ðop Þ

ð31Þ

  Z2 c r op ¼

ð32Þ

1

pop J21 ðop Þ þY 21 ðop Þ

The low- and high-frequency limits of the dimensionless spring and damping coefficients of radial and circumferential dynamic stiffness coefficient are listed in Table 1. Eqs. (19) and (27) are similar in form as the analytical solution of cylindrical elastic wave radiation problem in one-phase media, which can be used as an exact transmitting boundary in frequency domain [60,61]. Nevertheless, such transmitting boundary is temporally global because it corresponds to a convolution after taking the inverse Fourier transform from the frequency to the time domain. Table 1 Low- and high-frequency limits of the dimensionless spring and damping coefficients in Eqs. (19) and (27).

ð22Þ

ð23Þ

51

kr cr ky cy

o-0

o-N

1 0 1 0

1  Z2/4 Z2/2 0.75 0.5

52

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

4. Formation of high-order spring–damper–mass transmitting boundary conditions

Table 2 Rational approximations of dynamic stiffness coefficients for cylindrical P and SV wave radiation [70].

In order to represent the convolution in the stress– displacement relationship and to make it applicable in the time domain, several temporal localization methods have been developed, the procedure of which all consists of the rational function approximation in the frequency domain and the auxiliary variable realization into the time domain. The rational function approximation (also called the Pade approximation) has the form of polynomial-fraction, and can provide high accuracy and convergence rate with a relatively low order [39]. Through performing a proper kind of fraction expansion, the rational function approximation is equivalently realized into time domain, forming a second-order system of ordinary differential equations with auxiliary variables introduced. 4.1. Rational function approximation of radial and circumferential dynamic stiffness coefficients in the frequency domain The exact circumferential dynamic stiffness coefficients Eq. (21) and radial dynamic stiffness coefficients Eq. (28) have similar formulas. By introducing an instrumental variable H(o0): Hðo0 Þ ¼ o0

Hð2Þ 0 ðo0 Þ Hð2Þ 1 ðo0 Þ

ð33Þ

with its rational function approximation: N ~ ðo0 Þ ¼ io0 p0 þp1 io0 þ    þ pN ðio0 Þ Hðo0 Þ  H q0 þq1 io0 þ    þqN ðio0 ÞN

ð34Þ

where o0 is a dimensionless frequency, which can be approximated by the rational functions pffiffiffiffiffiffiffiof dimensionless frequency ios or iop, respectively, where i ¼ 1 is the imaginary unit. In Eq. (34), N is the order of the proposed transmitting boundary; pj (j¼1,y, N) and qj (j ¼1,y, N) are the real parameters to be determined. From Table 1, these rational functions naturally satisfy the exactness of dimensionless spring coefficient at low-frequency limit, and letting pN ¼qN ¼1 guarantees the exactness of dimensionless damping coefficient at high-frequency limit in both conditions. The remaining real parameters are determined by fitting the instrumental variable H(o0). Obviously, the corresponding parameters in the rational function approximation of both radial and circumferential dynamic stiffness coefficients are identical. Such dimensionless parameters are applicable widely in cylindrical wave propagation problems because the results are independent of the material and geometric constants. Different parameter identification methods have been proposed such as the method based on linearization and iteration [62,63] and the method using hybrid genetic-simplex optimization algorithm [64]. Here, the method proposed in [64] is used which results in a stable and accurate rational approximation. Providing the dimensionless frequency range of interest from 0 to 10 and the number of discrete frequency points being 0.01, the resulting rational approximation parameters in dynamic stiffness coefficients for cylindrical P and SV wave radiation are identified and listed in Table 2 with N ¼1, 2 and 3. In practical computation of cylindrical elastic wave radiation problem in one phase media, very accurate rational approximation can be obtained by using the order N ¼3 for radial and circumferential dynamic stiffness coefficients [45]. 4.2. Auxiliary variable realization into the time domain Based on the rational function approximation expression, several high-order accurate local time-domain transmitting

p0 p1 p2 q0 q1 q2

N ¼1

N¼2

N¼3

1.8786949E  01 – – 6.816956E  01 – –

4.5051674E  02 9.4617453E  01 – 3.9884816E  01 1.4453013Eþ 00 –

9.5891440E  03 5.4431167E  01 1.7748170E þ00 1.7545927E  01 1.3081337E þ00 2.2746358E þ00

boundaries as mechanical model have been developed, such as the consistent lumped-parameter model proposed by Wolf [41] and the systematic lumped-parameter models proposed by Wu and Lee [43]. These two models are both derived through performing partial-fraction expansion on the polynomial fraction. Recently a novel local time-domain transmitting boundary as mechanical model deduced from continued fraction expansion of the polynomial fraction has been proposed by Zhao [45]. This model has a simple combining form, and contains mass on each auxiliary degree of freedom. Thus, the total dynamic system which couples the finite element and the high-order transmitting boundary is a purely second-order system of ordinary differential equations, which have equivalent stability to that of rational function approximation and can be coded and solved by an implicit time integration method. This new model is applied here to construct the high-order transmitting boundary for saturated porous media. It is briefly summarized as follows. The configuration of proposed mechanical model is shown in Fig. 1, which consists of the lumped spring, dashpot and mass elements Kj,0, Cj,l and Mj,l, and the auxiliary degrees of freedom uj,l(t). Equations of motion of the proposed radial and circumferential transmitting boundaries are ðC j,0 þ C j,1 Þu_ j ðtÞ þ K j,0 uj ðtÞ ¼ C j,1 u_ j,1 ðtÞ þ f j ðtÞ

ð35Þ

Mj,l u€ j,l ðtÞ þ ðC j,l þ C j,l þ 1 Þu_ j,l ðtÞ ¼ C j,l u_ j,l1 ðtÞ þ C j,l þ 1 u_ j,l þ 1 ðtÞ,l ¼ 1,:::,N1

ð36Þ Mj,N u€ j,N ðtÞ þ C j,N u_ j,N ðtÞ ¼ C j,N u_ j,N1 ðtÞ

ð37Þ

with uj,0 ¼uj, where the first subscript j ¼ y or r represents, respectively, the values for circumferential and radial direction; uj(t) and fj(t) are, respectively, the displacement and stress on cylindrical transmitting boundary; N is the order of the transmitting boundary; uj,l(t) (l¼1,y, N) are the displacements of auxiliary degrees of freedom; a dot over variable denotes the derivative to time; Kj,0, Cj,l (l ¼0,y,N) and Mj,l (l ¼1,y,N) are, respectively, the undetermined spring, dashpot and mass parameters. By introducing the dimensionless real parameters k0, cl (l ¼0,y,N) and ml (l ¼1,y,N) which are independent of the material and geometric constants, and satisfy K j,0 ¼ S0j k0

ð38Þ

C j,l ¼ rV j cl ðl ¼ 0,:::,NÞ

ð39Þ

Mj,l ¼ rRml ðl ¼ 0,:::,NÞ

ð40Þ

where R is the radius location of the transmitting boundary; S0j is the static stiffness for the radial (j¼ r) and circumferential (j¼ y) boundaries and is equal to 2G/R; Vj denotes the wave speeds for P wave (j ¼r) and S wave (j¼ y), respectively. Through applying Fourier transform Eqs. (35)–(37) are rewritten in the frequency domain. After some manipulations, the dynamic stiffness coefficient of the mechanical model can be

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

fj (t)

Table 3 Dimensionless parameters of two-dimensional in-plane normal and tangential transmitting boundary conditions [45].

k0 c0 c1 m1 c2 m2 c3 m3

uj (t)

Cj,1 Kj,0

53

Cj,0

Mj,1

uj,1(t)

Cj,2

Mj,2

N¼1

N¼2

N¼3

1.0000000E þ00 2.7562626E  01 7.2437374E  01 1.0627400E þ 00 – – – –

1.0000000E þ00 1.1295445E  01 8.8704555E  01 1.5764527E þ 00 2.7396532E  01 3.8650309E  01 – –

1.0000000E þ00 5.4651680E  02 9.4534832E  01 1.7880150E þ00 3.8929369E  01 5.1546393E  01 1.7193378E  01 3.9127731E  01

coefficients (see Eqs. (20) and (28)) can be expanded into the following continued fraction by using the polynomial division   Sj oj

uj,2(t)

  k2 1 ioj 1 þ ¼ S0j 1 þ g 1 þ ioj h1 þð1=g 2 þioj h2 þ    ð1=g N þ ioj hN ÞÞ 2

ð45Þ

Cj,l

Mj,l

uj,l (t)

Cj,N

Mj,N

uj,N (t)

Fig. 1. Mechanical model of high-order local time-domain transmitting boundary.

derived as the function of the dimensionless frequency:

   k2 ioj ðc0 þc1 þA1 Þ S oj ¼ S0j k0 þ 2 Al ¼

AN ¼

c2l cl þ cl þ 1 þ ml ioj þ Al þ 1 c2N cN þmN ioj

ðl ¼ 0,:::,N1Þ

ð41Þ

ð42Þ

ð43Þ

where Al ¼  cluj,l(o)/uj,l  1(o) (l ¼1,y,N  1) with uj,0(o)¼uj(o) are the auxiliary variables introduced; oj denotes the dimensionless frequency, j ¼P or S for the P or S wave, respectively, (op and os are the dimensionless frequency introduced in Section 3); the constant k ¼ Z or 1 for the radial or circumferential boundaries, respectively, (Z ¼Vp/Vs has been introduced in Section 3). Eliminating all auxiliary variables, Eqs. (41)–(43) may be written as the following continued fraction:

where gl and hl are the newly introduced parameters during derivation of the continued fraction. The detailed procedure is described in the Appendix. The dimensionless parameters of the high-order transmitting boundaries are obtained by comparing the continued fraction forms of both the dynamic stiffness transmitting boundary Eq. (44) and the exact dynamic stiffness coefficient Eq. (45). It is obvious that k0 ¼1 and c0 þc1 ¼1. For computing the cl (l¼1,y,N) and ml (l ¼1,y,N) from the gl (l ¼1,y,N) and hl (l ¼1,y,N), the following two-step recursive algorithm is provided by Zhao [14]. First, starting from cN/cN  1 ¼  1/(1þgN  1gN), the cl þ 1/cl is computed recursively for l from N  2 to 1 "   #1 cl þ 1 cl þ 2 1 ðl ¼ N2,:::,1Þ ¼  1 þg l g l þ 1 1þ ð46Þ cl cl þ 1 Then, starting from c2/c1, the cl and ml are computed recursively for l from 1 to N   c2 1þ ð47Þ , m1 ¼ h1 c21 c1 ¼ g 1 1 c1 cl ¼

cl cl1

cl1 ,

ml ¼ hl1 hl

ðl ¼ 2,:::,N Þ

ð48Þ

The dimensionless parameters for the circumferential and radial boundaries with N ¼1, 2 and 3 are listed in Table 3. The physical parameters Kj,0, Cj,l (l ¼0,y,N) and Mj,l (l ¼1,y,N) can be obtained by using Eqs. (38)–(40).

5. Finite element implementation The proposed high-order local time-domain transmitting boundary for saturated porous media is derived from u–p formulation, with an assumption of zero permeability coefficient.

" !#   c21 k2 ioj c0 þc1 þ S oj ¼ S0j k0 þ 2 c1 þ c2 þ m1 ioj þðc22 =c2 þ c3 þ m2 ioj þ& þðc2N =cN þ mN ioj ÞÞ

Following the way as in [14], the rational function approximations for the circumferential and radial dynamic stiffness

c2l ml1

ð44Þ

Here the finite element program DIANA-SWANDYNE II developed by Chan at the University of Birmingham is used and redeveloped.

54

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

media, equations of motion of the transmitting boundaries can be written in matrix form in local coordinate as follows: ( f~ n ¼ Kn u~ n Cn u_~ n Mn u€~ n ð52aÞ f~ ¼ K u~ C u~_ M u~€

Scattered wave source High-order transmitting boundary

R

s

Near field The truncated boundary

Fig. 2. Schematic plot of N-order radial and circumferential high-order transmitting boundary.

5.1. Finite element scheme of the governing equations DIANA-SWANDYNE II is a two-dimensional program which incorporates plane strain and axisymmetric analysis for static, consolidation and dynamic problems in geomechanics. The governing equation that is being solved is the fully coupled Biot equation with the u–p simplification. The finite element scheme for the overall system of saturated porous linear elastic media can be written in the terms of the variable set ½u,pT as #" #

" € # " " # " f # _ 0 0 u M 0 u K -Q u s ð49Þ þ ¼ þ fp QT S p 0 0 p€ 0 H p_ where ui ¼ Nu u and pi ¼ Np p with uandp representing the vectors of solid displacement and pore fluid pressure on element nodes, and Nu and Np representing the shape functions of solid phase and R R fluid phase; M ¼ O NTu rNu dO is the mass matrix; K ¼ O BTu DBu dO R is the stiffness matrix;Q ¼ O BTu maNp dO is the coupling matrix R S ¼ O NTp Q1 Np dO is the mass with mT ¼(1,1,1,0,0,0); T   R  compressibility matrix; H ¼ O rNp kf rNp dO is the  T permeability matrix with r ¼ ð@=@xÞ,ð@=@yÞ,ð@=@zÞ . The external force vectors are Z Z ^ G þ NT rbdO NTu td ð49aÞ fs ¼ u Gs

O

and fp ¼ 

Z Gq

NTp q^ n dG þ

Z O

ðrNp ÞT kf rf bdO

ð49bÞ

where t^ and q^ n are the total external force vector applied on the boundary and the influx across the boundary, respectively. The link between the successive displacement and pore pressure is € n and Dp_ provided as the following two-order scheme where Du n are as yet undetermined quantities. € € € u n þ 1 ¼ un þ Dun _ € _ € u n þ 1 ¼ un þ un Dt þ b1 Dun Dt 1 € n Dt 2 þ 1b Du _ n Dt þ u € n Dt 2 un þ 1 ¼ un þ u 2 2 2 p_ n þ 1 ¼ p_ n þ Dp_ n p ¼ p þ p_ Dt þ b Dp_ Dt nþ1

n

n

1

ð50Þ

n

Such recurrence schemes are unconditional stable [66] when 1 2

b2 Z b1 Z , b1 Z

1 2

ð51Þ

5.2. Finite element implementation of the transmitting boundary After applying the N-order circumferential and radial transmitting boundaries (shown in Fig. 2) in two-dimensional planestrain cylindrical wave propagation problem in saturated porous

s

s

s

s

s

s

u~ s ¼ ðu~ s , u~ s,1 , u~ s,2 ,. . ., where u~ n ¼ ðu~ n , u~ n,1 , u~ n,2 ,. . ., u~ n,N1 , u~ n,N ÞT , u~ s,N1 , u~ s,N ÞT are displacement vectors for radial boundary and circumferential boundary, respectively; f~ n ¼ ðs~ n ,0,0,. . .,0,0ÞT , f~ s ¼ ðt~ ,0,0,. . .,0,0ÞT in which s~ n and t~ are the radial and circumferential stresses at the truncated boundary that are imposed by the transmitting boundary; 2 3 K j,0 0 0    0 0 6 0 0 0  0 07 6 7 6 7  6 0 0 0  0 07 6 7 Kj ¼ 6 ð52bÞ ^ ^ & ^ ^7 6 ^ 7 6 7 4 0 0 0  0 05 0 2 6 6 6  6 6 Cj ¼ 6 6 6 6 4

0

0

0



0

C j,0 þ C j,1

C j,1

0



0

C j,1

C j,1 þ C j,2

C j,2



0

0

C j,2

C j,2 þ C j,3



0

^

^

^

&

^

0

0

0



C j,N1 þC j,N

0

0

0



C j,N

0

3

7 7 7 0 7 7 ^ 7 7 7 C j,N 7 5 C j,N 0

ð52cÞ 2

0 60 6 6  6 60 Mj ¼ 6 6^ 6 60 4 0

0 M j,1

0 0

 

0 0

0 ^

Mj,2



0

^

&

^

0

0



Mj,N1

0

0



0

3

0 0 7 7 7 0 7 7 ^ 7 7 7 0 7 5 M j,N

ð52dÞ

where j ¼n, s denotes radial and circumferential mechanical model, respectively. After coordinate transformation Eq. (52a) is rewritten in the global coordinate. For example, the stress t^ at the truncated boundary that is imposed by the transmitting boundary is given as follows: " # ( ) ( ) u_ x ux t^ T ^t ¼ x ¼ PT BP P DP ð53Þ uy u_ y t^y where {(ux/uy)} is the displacement of the node on truncated " # 0 C n,0 þ C n,1 boundary in global coordinate; B ¼ and 0 C s,0 þ C s,1 " # 0 K n,0 ; P is the transformation matrix with the expression D¼ 0 K s,0 P¼



cos a

sin a

sin a

cos a

 ð54Þ

where a is the angle between the radial direction in the local coordinate and the positive x axis in the global coordinate. Eq. (53) can be called a kind of mix boundary conditions as it describes a relationship between stress and displacement on the boundary. When applying transmitting boundary, the computational region contains three kinds of boundary conditions: the natural boundary conditions Gs and Gq which describe the total stress and the fluid flux, respectively; the compulsive boundary conditions Gu and Gp which describe the displacement and pore fluid pressure, respectively, and the mix boundary conditions Ga as mentioned above. To ensure the existence and uniqueness of solutions, the components of the total boundary should meet the

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

following requirements:

Gs [ Gp [ Gu [ Gq [ Ga ¼ G Gu \ Gs ¼ 0Gu \ Ga ¼ 0Gs \ Ga ¼ 0 Gp \ Gq ¼ 0Gp \ Ga ¼ 0Gq \ Ga ¼ 0

ð55Þ

Using the same shape function as that in the adjacent interior element, which method is called the consistent implementation [65], the motion equations of the transmitting boundary in global coordinate can be discretized. The corresponding results are substituted into Eq. (49), and after some mathematical derivations, the finite element scheme for the overall system of the near filed and the high-order transmitting boundary becomes 32 _ 3 2 32 € 3 2 u u C þC 0 0 M 0 0 6 76 6 76 u € 7 _ 7 6 7 u 0 C b 7þ4 0 5 4 0 Mb 0 5 6 b 4 5 4 b 5 T _p €p 0 S Q 0 0 0 2 6 þ4

KþK

0

0

0

0

0

3 2 3 u fs 76 u 7 6 0 7 5 0 54 b 5 ¼ 4 fp p H

-Q

32

ð56Þ

where ub is the displacement vector for both radial and circumferential auxiliary degree of freedom; Mb and Cb are the mass matrix and damping matrix for the auxiliary degrees of freedom in the global coordinate; C and Kare the modified damping matrix and stiffness matrix corresponding to the degrees of freedom at

Truncated boundary

Standard transmitting boundary element

55

the truncated boundary and can be derived from Eq. (53) as: R C ¼  Ga NTu PT BPNb dG R ð57Þ K ¼  Ga NTu PT DPNb dG where Nb is the shape function of the boundary element. Due to the simple combining form of the proposed high-order transmitting boundary, only the ‘standard transmitting boundary element’ needs to be developed which is a spring–damper–mass mechanical model with an auxiliary degree of freedom as shown in Fig. 3. The standard element has only solid degrees of freedom. Higher order transmitting boundary is realized easily by assembling this standard element. During dynamic finite element analysis of saturated soil, the mixed 8-4 noded isoparametric quadrilateral element is usually adopted, i.e., 8 solid nodes located at the four corners and four middle points of each side of the rectangular element with quadratic interpolation, while 4 fluid nodes located at the four corners of the rectangular element with linear interpolation. A 3  3 order Gauss quadrature rule is adopted for the numerical integration. The standard transmitting boundary element is simulated by two line elements connected in parallel, with three solid nodes and three gauss points in each line. For cylindrical wave propagation problem, the mechanical model of the high-order transmitting boundary is dynamically stable with all physical parameters positive. So the unconditionally stable implicit scheme proposed by Zienkiewicz [66] is still valid as long as the requirement of Eq. (51) is met. Compared with the explicit integration scheme [67,45] used in cylindrical wave propagation problems of unbounded one-phase media, the limitation on the time step size in this implicit scheme is only the required accuracy. Clearly in dynamic, earthquake problem in saturated soil small time steps will generally be used to follow the time characteristic of the input motion, while after the passage of earthquake or dynamic load, the remaining motion is caused by consolidation process which has a slower response allowing larger time steps to be used. The unconditional stable implicit scheme presented here has the special advantage in this aspect.

6. Numerical examples To illustrate the efficiency of the proposed transmitting boundary, several plane-strain transient dynamic problems in unbounded fluid-saturated porous media are investigated. First, a cylindrical elastic wave radiation problem is analyzed to demonstrate the effectiveness of high-order transmitting boundary in absorbing outgoing cylindrical P and SV waves. Second, a general

10

10

8

8

Amplitude(103)

Load(104N/m)

Fig. 3. Standard transmitting boundary element.

6 4 2 0 0.0

6 4 2

0.2

0.4

0.6

0.8

1.0

Time(s)

1.2

1.4

1.6

1.8

2.0

0

0

10

20

30

Frequency(Hz)

Fig. 4. Triangular impulse load: (a) time history and (b) amplitude of Fourier spectrum.

40

50

56

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

semi-infinite wave problem that contains the more complex outgoing wave fields with the box-shape transmitting boundary is analyzed. Comparisons are also made with the viscous–spring transmitting boundary proposed by Liu [57]. The material under consideration is medium-dense sand. In all cases, the material properties are homogeneous with no material damping and nonlinearity for simplicity, which are chosen as follows: shear modulus G¼79.6 MPa; poison ratio u¼0.4; Kb ¼371.47 MPa, Kf ¼1.848 GPa, Ks ¼ 36 GPa; void ratio e¼0.4; rs ¼2650 kg/m3, rf ¼1000 kg/m3. The wave speeds are: Vp ¼1466.4 m/s, Vs ¼191 m/s. The parameters in time recurrence scheme are chosen as that proposed by Dewoolkar [68]: b1 ¼0.6, b2 ¼0.605 and b1 ¼ 0:6. The unconditional stability is maintained and usually some algorithmic damping is introduced by using such values.

The cylindrical P wave radiation problem is first analyzed. The problem is the two-dimensional plane strain saturated porous elastic unbounded media with a circular hole of radius 2 m on which the uniform triangular impulse load as shown in Fig. 4 is imposed on the radial direction. The time step is 0.01 s in order to follow the characteristic of this dynamic load. As it is symmetric, only a quarter of the truncated domain is analyzed, which consists of 10  10 elements. The element size is 2 m on the direction of wave propagation, which satisfies the element size requirement suggested by Kuhlemeyer and Lysmer [69]: le r lmin =M

Load A

Symmetric boundary 20m

B Transmitting boundary

Fig. 5. Schematic plot of finite element model of motions of cylindrical P wave.

0.5

Displacement (10-5m)

Displacement (10-3m)

0.1

0.0

-0.1 Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-0.2

-0.3 0.0

ð58Þ

where le is the element size, lmin is the minimum wavelength, and M should be larger than 8 for two-dimensional problems. The surface on which dynamic load is applied is impermeable. The transmitting boundaries, are applied at R¼22 m in the radial direction. The proposed high-order transmitting boundary and the viscous–spring boundary are applied separately in order to compare their results. The reference solution for this model is obtained by using an extended mesh (400  10 elements) with a fixed boundary, in which no reflected wave from the truncated boundary travels back into the domain of interest within the analyzed interval. The permeability coefficient in this example is 1  10  13 m/s. The schematic diagram of finite element model and boundary conditions of the truncated domain is shown in Fig. 5. Figs. 6 and 7 show the time histories of the vertical displacements and pore fluid pressures, respectively, at the points A and B which are indicated in Fig. 5. Clearly, the solutions using the highorder transmitting boundary with N ¼1, 2 and 3 agree very well

20m Symmetric boundary

C

6.1. Cylindrical P wave radiation problem

0.1

0.2

0.3

0.4

0.0 -0.5 -1.0 -1.5 -2.0 0.0

0.5

Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

0.1

0.2

Time (s)

0.3

0.4

0.5

Time (s)

Fig. 6. Time histories of vertical displacement of cylindrical P wave radiation problem at (a) point A and (b) point B.

1.5 1.0

Pore fluid pressure (103Pa)

Pore fluid pressure (103Pa)

2.0

0.0

Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-2.0

-4.0 0.0

0.1

0.2

0.3

Time (s)

0.4

0.5 0.0 -0.5 Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-1.0 -1.5

0.5

-2.0 0.0

0.1

0.2

0.3

0.4

0.5

Time (s)

Fig. 7. Time histories of pore fluid pressure of cylindrical P wave radiation problem at (a) point A and (b) point B.

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

with the extended mesh solution, even at the truncated boundary, and is more accurate than that using the viscous-spring boundary. The results obtained by using high-order transmitting boundary with N¼1, 2 and 3 are close, which again demonstrates the high accuracy and convergence rate of rational function approximation. The second step in this analysis concerns the application scope of the proposed high-order transmitting boundary, as well as the viscous-spring boundary for saturated porous media based on the assumption of zero permeability coefficient in u–p formulation. In fact, for the weak permeability, nearly no relative motion between the fluid and the solid occurs and the saturated porous media behaves like an undrained material. When the actual permeability is moderate, the accuracy of the proposed transmitting boundaries are proved by considering three conditions with the same material characteristics but gradually increasing permeabilities: k¼1  10  7,

57

1  10  2 and 1 m/s, respectively. The vertical displacement and pore fluid pressure of point C for these three cases are shown in Figs. 8–10. From these results we see that the efficiency of the proposed transmitting boundaries depends on the permeability. Figs.8 and 9 show that when permeability is moderate, accurate results are provided by the high-order transmitting boundaries with N¼1, 2 and 3 and are much more accurate than the results obtained by the viscous-spring transmitting boundary. In dynamic analysis of saturated porous media it is important to simulate the pore fluid pressure accurately, and the results with the high-order transmitting boundaries agree very well with the reference solutions, while the accuracy of the results with the viscous-spring transmitting boundary reduces distinctly. Therefore, although the model was developed assuming zero permeability coefficient, it works correctly

0.5

1.5

0.0

1.0

Pore fluid pressure (103Pa)

Displacement (10-5m)

Fig. 8. Time histories of (a) vertical displacement and (b) pore fluid pressure at point C of cylindrical P wave radiation problem when permeability is 1  10  7 m/s.

-0.5 -1.0 Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-1.5 -2.0 -2.5 0.0

0.1

0.2

0.3 Time (s)

0.4

0.5 0.0 -0.5

-1.5 -2.0 0.0

0.5

Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-1.0

0.1

0.2

0.3 Time (s)

0.4

0.5

0.5

1.5

0.0

1.0

Pore fluid pressure (103Pa)

Displacement (10-5m)

Fig. 9. Time histories of (a) vertical displacement and (b) pore fluid pressure at point C of cylindrical P wave radiation problem when permeability is 1  10  2 m/s.

-0.5 -1.0 -1.5 Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-2.0 -2.5 -3.0 0.0

0.1

0.2

0.3 Time (s)

0.4

0.5

0.5 0.0 -0.5 Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-1.0 -1.5 -2.0 0.0

0.1

0.2 0.3 Time (s)

0.4

0.5

Fig. 10. Time histories of (a) vertical displacement and (b) pore fluid pressure at point C of cylindrical P wave radiation problem when permeability is 1 m/s.

58

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

for moderate permeabilities. However, for the highest permeability in this numerical example, Fig. 10 shows that the results provided by both transmitting boundaries become significantly worse. Fig. 4 indicates the characteristic frequency range of the applied load is from 0 to 10 Hz, so the maximum load frequency can be treated as f¼10 Hz. For the highest permeability in this numerical example, the characteristic frequency is fc ¼n/2prfkf ¼ng/2pk¼0.446 Hz, and the corresponding logarithm of the dimensionless frequency ratio is log (f/fc)¼1.35, which is beyond the applicability scope of the u–p formulation [59]. That is, the reference solution obtained from extended mash based on the u–p formulation is no longer reliable. As a result, the proposed transmitting boundary derived based on the u–p formulation is also invalid. Thus, this example demonstrates that for moderate permeability within the framework of geomechanics, the proposed transmitting boundaries which are derived based

6.2. Cylindrical SV wave radiation problem The cylindrical SV wave radiation problem is then analyzed. The problem is still the two-dimensional plane strain saturated porous elastic unbounded media with a circular hole of radius 2m, on which the uniform triangular impulse load as shown in Fig. 4 is imposed on the circumferential direction. The truncated domain consists of 40  10 elements, with the element size still 2 m on the radial direction which satisfies the element size requirement. The transmitting boundaries, namely the viscous-spring boundary or the high-order transmitting boundary, are located at R¼22 m on the circumferential direction. The reference solution for this model is obtained by using an extended mesh (300  40 elements) with a fixed boundary. The time step is still 0.01 s while the permeability coefficient in this example is 1  10  7 m/s. The schematic diagram of the finite element model with boundary conditions of the truncated domain is shown in Fig. 11. The horizontal displacements of point A and B as indicated in Fig. 11 are shown in Fig. 12. Clearly, the solutions using the high-order transmitting boundary with N¼1, 2 and 3 agree very well with the extended mesh solution and exhibit better absorption of outgoing cylindrical SV

C

D

B

on the assumption of zero permeability exhibit well absorption of outgoing cylindrical wave propagating in saturated porous media in the condition when u–p formulation is applicable. As noticed by Gajo et al. [48,59], the transition zone of the dynamic behaviors of saturated soils between the low frequency range and high frequency range is narrow, and for the low frequency conditions, the P1 wave speeds hardly changes and the P2 wave dies out quickly. This assumption is in agreement with the other authors [50,51] and the results indicate again that neglecting the second dilatational wave effect can still lead to acceptable solutions in usual situations such as geotechnical earthquake engineering.

A Transmitting boundary Fig. 11. Schematic plot of motions of cylindrical wave propagation problems.

0.5 0.3

Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

0.3

Displacement (10-4m)

Displacement (10-4m)

0.4

0.2 0.1 0.0 -0.1 -0.2 0.0

Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

0.2

0.1

0.0

-0.1 0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

0.0 0.2

0.4

0.6

Time (s)

0.8

1.0

1.2

1.4

1.6

1.8

2.0

Time (s)

0.1

0.1

0.0

0.0

Displacement (10-4m)

Displacement (10-4m)

Fig. 12. Time histories of horizontal displacement of cylindrical SV wave radiation problem at (a) point A and (b) point B.

-0.1 -0.2 -0.3 Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-0.4 -0.5 -0.6 0.0

0.1

0.2

0.3

Time (s)

0.4

0.5

-0.1 -0.2 Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-0.3 -0.4

0.6

-0.5 0.0

0.1

0.2

0.3

0.4

0.5

0.6

Time (s)

Fig. 13. Time histories of vertical displacement of combining motions of cylindrical P and SV waves radiation problem at (a) point C and (b) point D.

1.5

1.5

1.0

1.0

Pore fluid pressure (103Pa)

Pore fluid pressure (103Pa)

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

0.5 0.0 -0.5

Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-1.0 -1.5 -2.0 0.0

0.1

0.2

0.3

0.4

0.5

0.5 0.0 -0.5

Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-1.0 -1.5

0.6

59

-2.0 0.0

0.1

Time (s)

0.2

0.3

0.4

0.5

0.6

Time (s)

Fig. 14. Time histories of pore fluid pressure of combining motions of cylindrical P and SV waves radiation problem at (a) point C and (b) point D.

y Load

x

A B

C

20m Symmetric boundary

Transmitting boundary

Fig. 15. Schematic plot of two-dimensional plane-strain transient dynamic problem.

wave propagating in saturated porous media than the viscousspring transmitting boundary. 6.3. Cylindrical P and SV wave coupling radiation problem The coupling of cylindrical P wave and SV wave radiation problem is now analyzed. The finite element model is the same as used in Section 6.2. The uniform triangular impulse load shown in Fig. 4 is imposed on the direction which has an included angle of 301 with the radial direction. The transmitting boundaries, namely the viscous-spring boundary or the high-order transmitting boundary, are located at R¼22 m on both the circumferential and radial directions. The vertical displacements and pore fluid pressures of point C and D as indicated in Fig. 11 are shown in Figs. 13 and 14. Again, the high-order transmitting boundary with N ¼1, 2 and 3 give satisfied results which agree very well with the extended mesh solution, and exhibit high accuracy and convergence rate. 6.4. General two-dimensional plane-strain transient dynamic problem The general plain-strain transient dynamic problem in a semiinfinite fluid-saturated porous medium is finally analyzed. The schematic diagram of the finite element mesh and the corresponding

boundary conditions are shown in Fig. 15. The truncated problem domain consists of 10  10 elements with the element size of 2 m  2 m. The uniform triangular impulse load as shown in Fig. 4 is imposed on the surface within the range of X¼0  2 m. The time step is 0.01 s. Due to the symmetry of problem, only the left half of the domain is analyzed. The surface on which dynamic load is applied is impermeable. The transmitting boundaries, i.e., the viscous–spring boundary or the high-order transmitting boundary, are located at X¼Y¼  20 m as shown in Fig. 15. The reference solution for this model is obtained by using an extended mesh with a fixed boundary at X¼Y¼  820 m. The permeability coefficient in this example is 1  10  13 m/s. Figs. 16 and 17 show the time histories of the vertical displacements and pore fluid pressures at the points A and B, respectively, which are indicated in Fig. 15. Clearly, as the proposed transmitting boundaries are derived from cylindrical wave propagation problem of saturated porous media, the direct application of the high-order transmitting boundaries with N ¼1, 2 and 3 as well as the viscous-spring boundary to general twodimensional wave propagation problems are not so accurate, and no improvement of accuracy is achieved by increasing the order of the transmitting boundaries. This unexpected result accords with the conclusion in one phase media [45]: a transmitting boundary that is high-order accurate for the symmetric case of wave radiation is not necessarily high-order accurate for the general infinite wave problem. Fortunately, when the high-order transmitting boundaries with N ¼1, 2 and 3 as well as the viscous–spring boundary are set sufficiently away from the wave source, solutions with acceptable engineering accuracy may still be achieved, and the accuracy slightly increases as the transmitting boundary moves away from the wave source from X¼Y ¼ 100 m to  200 m, as shown in Figs. 18 and 19. Further modification of the parameters of the viscous-spring transmitting boundary for saturated porous media can be investigated, as already carried out for one phase media [65].

7. Conclusions A high-order accurate local time-domain transmitting boundary as mechanical model for simulating the transient scalar wave propagation in unbounded saturated porous media, based on the u–p formulation with an assumption of zero permeability coefficient, is derived from the cylindrical elastic wave radiation problem. The transmitting boundaries are incorporated into the DIANA SWANDYNE II program and an unconditionally stable implicit time integration algorithm is adopted. The good waveabsorbing capabilities of this high-order transmitting boundary are demonstrated by several different numerical examples. Despite the

60

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

3.5 Pore fluid pressure (103Pa)

D-0.1isplacement (10-3m)

0.1

0.0

0.1 Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-0.2

3.0 2.5 Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

2.0 1.5 1.0 0.5 0.0

-0.3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

-0.5 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 Time (s)

Time (s)

Fig. 16. Time histories of (a) vertical displacement and (b) pore fluid pressure at point A.

1.2

0.5

1.0

Pore fluid pressure (103Pa)

Displacement (10-4m)

0.0

-0.5

Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-1.0

-1.5

Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

0.8 0.6 0.4 0.2 0.0 -0.2

-2.0 0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

-0.4 0.0

2.0

0.2

0.4

0.6

0.8

Time (s)

1.0

1.2

1.4

1.6

1.8

2.0

Time (s)

0.5

0.8

0.0

0.6

Pore fluid pressure (103Pa)

Displacement (10-4m)

Fig. 17. Time histories of (a) vertical displacement and (b) pore fluid pressure at point B.

-0.5

-1.0

Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-1.5

-2.0 0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

0.4

0.2

0.0

-0.2 0.0

2.0

Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

Time (s)

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

Time (s)

Fig. 18. Time histories of (a) vertical displacement and (b) pore fluid pressure at point C with the transmitting boundary located at X¼ Y¼  100 m.

0.8

0.0 -0.5 -1.0 -1.5

Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

-2.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 Time (s)

Pore fluid pressure (103Pa)

Displacement (10-4m)

0.5

0.6 0.4

Reference solution Viscous-spring boundary High-order boundary (N=1) High-order boundary (N=2) High-order boundary (N=3)

0.2 0.0 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 Time (s)

Fig. 19. Time histories of (a) vertical displacement and (b) pore fluid pressure at point C with the transmitting boundary located at X¼ Y¼  200 m.

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

assumption made in the derivation of this transmitting boundary, results have shown that it can provide sufficiently accurate results in case the u–p formulation of the basic equations is applicable. Although numerical results show that the direct applications of the proposed transmitting boundary to general two dimensional wave problems in infinite saturated porous media are not so accurate, solutions of acceptable engineering accuracy may still be achieved by setting the transmitting boundary sufficiently away from the scatter. Meanwhile, further modification of the parameters of the viscous–spring transmitting boundary for saturated porous media can be investigated, which will encourage the use of the u–p finite element formulation for dynamic analysis of soil–structure interaction problems in the time domain, taking into account the ‘twophase’ behavior of saturated soils.

Acknowledgements The work reported in this paper was financially supported by the National Natural Science Foundation of China (Project No. 51038007) and the Autonomous Research Program of Tsinghua University (Project No. 20121087977).

Appendix. Continued fraction of the rational function approximation By extracting the highest-order term in its numerator, Eq. (34) can be written as N1 ~ o0 Þ pð1Þ þpð1Þ io0 þ    þ pð1Þ Hð Pð1Þ ðio0 Þ N1 ðio0 Þ ¼ 1 þ 0 ð1Þ 1 ð1Þ ¼ 1þ ð1Þ ð1Þ N io0 Q ðio0 Þ q þ q io0 þ    þ q ðio0 Þ 0

N

1

ðA1aÞ with ¼ qj qð1Þ j

pð1Þ ¼ pj qj j

j ¼ 0,:::,N,

j ¼ 0,:::,N1

ðA1bÞ

where the superscript denotes the times of the division operation performed. Starting from P(1)(io0)/Q(1)(io0), by removing the first and second highest order terms of the new rational function, the following division process can be performed: P ðlÞ ðio0 Þ ðlÞ

Q ðio0 Þ

¼

1 g l þ hl io0 þ ðP

ðl þ 1Þ

ðio0 Þ=Q ðl þ 1Þ ðio0 ÞÞ

l ¼ 1,:::,N

ðA2aÞ

with the new rational function

ðlÞ

Q ðio

o0 ÞNl o0 ÞNl þ 1

ðlÞ ðlÞ pðlÞ 0 þp1 i 0 þ    þ pNl ði ðlÞ ðlÞ ðlÞ q0 þq1 i 0 þ    þ qNl þ 1 ði 0Þ P ðN þ 1Þ ði 0 Þ l ¼ 1,:::,N ðN þ 1Þ ¼0 ði 0 Þ Q

P ðlÞ ðio0 Þ

o

¼

o

o o

ðA2bÞ

and with the new parameters hl ¼

ðlÞ qNl þ1 ðlÞ pN1

qðlj þ 1Þ ¼ plj

,g l ¼

ðlÞ qNl ðlÞ pN1



j ¼ 0,:::,Nl,

ðlÞ pNl1 qðlÞ Nl þ 1

pðlÞ2 N1

ðA2cÞ

pðlj þ 1Þ ¼ qðlÞ hj pðlÞ g j pðlÞ j j1 j

j ¼ 0,:::,Nl1

ðA2dÞ

pðlÞ 1

¼ 0. By performing the above division process recurwhere sively, the rational approximation Eq. (34) is finally expanded into the following continued fraction: ~ ðo0 Þ ¼ io0 þ H

io0 g 1 þ io0 h1 þ ð1=g 2 þ io0 h2 þ &ð1=g N þ io0 hN ÞÞ ðA3Þ

61

References [1] Biot MA. Theory of propagation of elastic waves in a water-saturated porous solid. Journal of the Acoustical Society of America 1956;28:168–78. [2] Biot MA. Theory of propagation of elastic waves in a water–saturated porous solid. Journal of the Acoustical Society of America 1956;28:179–91. [3] Ishihara K. Approximation forms of wave equations for water-saturated porous materials and related dynamic moduli. Soils and Foundations 1970;10:10–38. [4] Stoll RD. Acoustic waves in ocean sediments. Geophysics 1977;42:715–25. [5] Ghaboussi J, Wilson EL. Variational formulation of dynamics of fluidsaturated porous elastic solid. Journal of Engineering Mechanics Division, ASCE 1972;104(4):947–63. [6] Zienkiewicz OC, Shiomi T. Dynamic behavior of saturated porous media: the generalized Biot formulation and its numerical solution. International Journal for Numerical and Analytical Methods in Geomechanics 1984;8:71–96. [7] Prevost JH. Wave propagation in fluid-saturated porous media: an efficient finite element procedure. International Journal of Soil Dynamics and Earthquake Engineering 1985;4:183–202. [8] Simon BR, Zienkiewicz OC, Paul DK. Evaluation of u–w and u–p finite element methods for the dynamic response of saturated porous media using onedimensional models. International Journal for Numerical and Analytical Methods in Geomechanics 1986;10:461–82. [9] Kausel E. Local transmitting boundaries. Journal of Engineering Mechanics, ASCE 1988;114:1011–27. [10] Givoli D. Numerical methods for problems in infinite domains. Amsterdam: Elsevier; 1992. [11] Givoli D. Non-reflecting boundary conditions: a review. Journal of Computational Physics 1991;94:1–29. [12] Hagstrom T. Radiation boundary conditions for the numerical simulation of waves. Acta Numerica 1999;8:47–106. [13] Givoli D. High-order local non-reflecting boundary conditions: a review. Wave Motion 2004;39:319–26. [14] Zhao M, Du X. Explicit finite element artificial boundary scheme for transient scalar waves in two-dimensional unbounded waveguide. International Journal for Numerical Methods in Engineering 2011;87:1074–104. [15] Hall WS, Oliveto G. Boundary element methods for soil–structure interaction. Dordrecht: Kluwer Academic Publishers; 2003. [16] Kausel E. Thin-layer method: formulation in the time domain. International Journal for Numerical Methods in Engineering 1994;37:927–41. [17] Givoli D, Cohen D. Nonreflecting boundary conditions based on Kirchhofftype formulae. Journal of Computational Physics 1995;117:102–13. [18] Teng ZH. Exact boundary condition for time-dependent wave equation based on boundary integral. Journal of Computational Physics 2003;190:398–418. [19] Givoli D. Recent advances in the DtN FE method. Archives of Computational Methods in Engineering 1999;6:71–116. [20] Givoli D. Exact representations on artificial interfaces and applications in mechanics. Applied Mechanics Reviews 1999;52:333–49. [21] Lysmer J, Kuhlmeyer RL. Finite dynamic model for innite media. Journal of the Engineering Mechanics Division, ASCE 1969;95:859–77. [22] Kuhlmeyer RL, Lysmer J. Finite element method accuracy for wave propagation problems. Journal of the Soil Mechanics and Foundations Division, ASCE 1973;99:421–7. [23] Deeks AJ, Randolph MF. Axisymmetric time-domain transmitting boundaries. Journal of Engineering Mechanics, ASCE 1994;120:25–42. [24] Kellezi L. Local transmitting boundaries for transient elastic analysis. Soil Dynamics and Earthquake Engineering 2000;19:533–47. [25] Liu J, Du Y, Du X, Wang Z, Wu J. 3D viscous–spring artificial boundary in time domain. Earthquake Engineering and Engineering Vibration 2006;5(1):93–102. [26] Smith WD. A nonreflecting plane boundary for wave propagation problems. Journal of Computational Physics 1974;15:492–503. [27] Kunar RR, Marti J. Computational methods for infinite domain media– structure interaction. Journal of the Applied Mechanics Division (ASME) 1981;46:183–204. [28] Liao Z, Wong H. A transmitting boundary for the numerical simulation of elastic wave propagation. Soil Dynamics and Earthquake Engineering 1984;3:174–83. [29] Liao Z. Extrapolation nonreflecting boundary conditions. Wave Motion 1996;24:117–38. [30] Clayton J, Enquist B. Absorbing boundary conditions for acoustic and elastic wave equations. Bulletin of the Seismological Society of America 1977;67:1529–41. [31] Engquist B, Majda A. Radiation boundary conditions for acoustic and elastic wave calculations. Communications on Pure and Applied Mathematics 1979;32:313–57. [32] Higdon RL. Numerical absorbing boundary conditions for the wave equation. Mathematics of Computation 1987;49:65–90. [33] Higdon RL. Radiation boundary conditions for elastic wave propagation. SIAM Journal on Numerical Analysis 1990;27:831–70. [34] Liu J, Lu Y. A direct method for analysis of dynamic soil–structure interaction based on interface idea. In: Chuhan Zhang, Wolf JP, editors. Dynamic soil– structure interaction. International Academic Publishers; 1997. p. 258–73. [35] Zienkiewicz OC, Emson C, Bettess P. A novel boundary infinite element. International Journal for Numerical Methods in Engineering 1983;19:393–404.

62

P. Li, E.-x. Song / Soil Dynamics and Earthquake Engineering 48 (2013) 48–62

[36] Chang YM. The use of innite element. Computers and Geotechnics 1996;18: 65–70. [37] Astley RJ. Infinite elements for wave problems: a review of current formulations and an assessment of accuracy. International Journal for Numerical Methods in Engineering 2000;49:951–76. [38] Zhao C, Valliappan S. A dynamic infinite element for three-dimensional infinite-domain wave problems. International Journal for Numerical Methods in Engineering 1993;36:2567–80. [39] Baker GA, Graves-Morris P. Pade approximants. 2nd ed. Cambridge: Cambridge University Press; 1996. [40] Wolf JP. Foundation vibration analysis using simple physical models. Englewood Cliffs, NJ: Prentice-Hall; 1994. [41] Wolf JP. Consistent lumped-parameter models for unbounded soil: physical representation. Earthquake Engineering and Structural Dynamics 1991;20: 11–32. [42] Wolf JP, Paronesso A. Errata: consistent lumped-parameter models for unbounded soil. Earthquake Engineering and Structural Dynamics 1991;20:597–9. [43] Wu WH, Lee WH. Systematic lumped-parameter models for foundations based on polynomial-fraction approximation. Earthquake Engineering and Structural Dynamics 2002;31:1383–412. [44] Wu WH, Lee WH. Nested lumped-parameter models for foundation vibrations. Earthquake Engineering and Structural Dynamics 2004;33:1051–8. [45] Du X, Zhao M. A local time-domain transmitting boundary for simulating cylindrical elastic wave propagation in infinite media. Soil Dynamics and Earthquake Engineering 2010;30:937–46. [46] Li P, Song E. Compressional wave velocity and its physical nature in saturated soils with extreme permeability values. Chinese Journal of Rock and Soil Mechanics 2012;33(7):1979–85. [47] Garg SK, Nayfeh H, Good AJ. Compressional waves in fluid-saturated elastic porous media. Applied Physics 1974;45:1968–74. [48] Gajo A. Influence of viscous coupling in propagation of elastic-waves in saturated soil. Journal of Geotechnical Engineering, ASCE 1995;19:399–414. [49] Modaressi H, Benzenati I. Absorbing boundary element for dynamic analysis of two-phase media, In: Proceedings of the world conference on earthquake engineering; 1992. p. 1157–1161. [50] Modaressi H, Benzenati I. Paraxial approximation for poroelastic media. Soil Dynamics and Earthquake Engineering 1994;13:117–29. [51] Akiyoshi T, Fuchida K, Fang HL. Absorbing boundary conditions for dynamic analysis of fluid-saturated porous media. Soil Dynamics and Earthquake Engineering 1994;13:387–97. [52] Akiyoshi T, Sun X, Fuchida K. General absorbing conditions for dynamic analysis of fluid-saturated porous media. Soil Dynamics and Earthquake Engineering 1998;17:397–406. [53] Degrande G, De Roeck G. An absorbing boundary condition for wave propagation in saturated porous media. Part I: Formulation and efficiency evaluation. Soil Dynamics and Earthquake Engineering 1993;12:411–21. [54] Degrande G, De Roeck G. An absorbing boundary condition for wave propagation in saturated porous elastic media, Part II: Finite element formulation. Soil Dynamics and Earthquake Engineering 1993;12:423–32.

[55] Gajo A, Saetta A, Vitaliani R. Silent boundary conditions for wave propagation in saturated porous media. International Journal for Numerical and Analytical Methods in Geomechanics 1996;13:253–73. [56] Zohra Zerfa BenjaminLoret. A viscous boundary for transient analyses of saturated porous media. Earthquake Engineering and Structural Dynamics 2004;33:89–110. [57] Liu G, Song E. Visco-elastic transmitting boundary for numerical analysis of infinite saturated soil foundation. Chinese Journal of Geotechnical Engineering 2006;28(12):2128–33. [58] Wang Z, Zhao CA. Viscous-spring dynamical artificial boundary for saturated porous media. Chinese Journal of Theoretical and Applied Mechanics 2006;38(5):605–11. [59] Yang J. Dynamic analysis of saturated soil and its application in calculating dynamic impedance of pile. PhD thesis. Department of Civil Engineering, Tsinghua University, Beijing, China, 1996.(in Chinese). [60] Novak M. Vertical vibration of floating piles. Journal of the Engineering Mechanics Division, ASCE 1977;103(1):153–68. [61] Novak M, Mitwally H. Transmitting boundary for axisymmetrical dilation problems. Journal of the Engineering Mechanics Division, ASCE 1988;114(1): 181–7. [62] Alpert B, Greengard L, Hagstrom T. Nonreflecting boundary conditions for the time-dependent wave equation. Journal of Computational Physics 2002;180:270–96. [63] Birk C, Ruge P. Representation of radiation damping in a dam–reservoir interaction analysis based on a rational stiffness approximation. Computers and Structures 2007;85:1152–63. [64] Du X, Zhao M. Stability and identification for rational approximation of frequency response function of unbounded soil. Earthquake Engineering and Structural Dynamics 2010;39:165–86. [65] Liu J, Gu Y, Du Y. Consistent viscous-spring artificial boundaries and viscousspring boundary elements. Chinese Journal of Geotechnical Engineering 2006;28(9):1070–5. [66] Zienkiewicz OC, Chan AHC, Pastor M, Schrefler BA, Shiomi T. Computational geomechanics with special reference to earthquake engineering. New York: John Wiley and Sons; 1998. [67] Wang J, Zhang C, Du X. An explicit integration scheme for solving dynamic problems of solid and porous media. Journal of Earthquake Engineering 2008;12:293–311. [68] Dewoolkar MM. A study of seismic effects on centiliver-retaining walls with saturated backfill. PhD thesis. Department of Civil Engineering, University of Colorado, Boulder, USA; 1996. [69] Kuhlemeyer RL, Lysmer J. Finite element method accuracy for wave propagation problems. Journal of Soil Mechanics and Foundations Division, ASCE 1973;99(SM5):421–7. [70] Zhao M. Stress-type time-domain artificial boundary condition for finiteelement simulation of near-field wave motion and its engineering application. PhD thesis, Department of Bridge and Tunnel Engineering, Beijing University of Technology, Beijing, China, 2009. (in Chinese).