Seismic response analysis of soil–structure interactive system using a coupled three-dimensional FE–IE method

Nuclear Engineering and Design 240 (2010) 1949–1966

Contents lists available at ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

Seismic response analysis of soil–structure interactive system using a coupled three-dimensional FE–IE method Jeong-Soo Ryu a , Choon-Gyo Seo b,∗ , Jae-Min Kim c , Chung-Bang Yun d a

Research Reactor Design & Engineering Division, Korea Atomic Energy Research Institute, Daejeon, Republic of Korea Steel Technology Department, Samsung Engineering Co. Ltd., Seoul, Republic of Korea Department of Civil & Environmental Engineering, Chonnam National University, Yeosu, Republic of Korea d Department of Civil & Environmental Engineering, KAIST, Daejeon, Republic of Korea b c

a r t i c l e

i n f o

Article history: Received 10 November 2009 Received in revised form 14 March 2010 Accepted 17 March 2010

a b s t r a c t This paper proposes a slightly new three-dimensional radial-shaped dynamic infinite elements fully coupled to finite elements for an analysis of soil–structure interaction system in a horizontally layered medium. We then deal with a seismic analysis technique for a three-dimensional soil–structure interactive system, based on the coupled finite–infinite method in frequency domain. The dynamic infinite elements are simulated for the unbounded domain with wave functions propagating multi-generated wave components. The accuracy of the dynamic infinite element and effectiveness of the seismic analysis technique may be demonstrated through a typical compliance analysis of square surface footing, an L-shaped mat concrete footing on layered soil medium and two kinds of practical seismic analysis tests. The practical analyses are (1) a site response analysis of the well-known Hualien site excited by all travelling wave components (primary, shear, Rayleigh waves) and (2) a generation of a floor response spectrum of a nuclear power plant. The obtained dynamic results show good agreement compared with the measured response data and numerical values of other soil–structure interaction analysis package. © 2010 Elsevier B.V. All rights reserved.

1. Introduction Many attempts have been made to apply artificial boundaries directly coupled to the finite element method to solve the unbounded diffusion and radiation problems such as elasto-statics, hydro-dynamics, elasto-dynamics, heat transfer, acoustics, electromagnetic field problems, etc. Artificial boundaries have been proposed to prevent a spacious reflection at the truncated boundary by artificially simulating the radiation and diffusion phenomena toward the infinite domain. Following much pioneering research over the past four decades, many kinds of modeling techniques are available for viscous boundary (White et al., 1969), transmitting boundary (Tassoulas, 1983), hybrid modeling method (Tzong and Penzien, 1986), boundary element method (Chen and Penzien, 1986), boundary solution method (Luco, 1974; Wong and Luco, 1985), scaled boundary finite element (Song and Wolf, 1998), infinite element method (Ungless, 1973; Bettess, 1977; Astley, 1983) and so forth. In this paper, we employ the coupled finite element–infinite element (FE–IE) method to analyze the soil–structure interaction

∗ Corresponding author. Tel.: +82 2 2148 2218; fax: +82 2 3458 4060. E-mail address: [email protected] (C.-G. Seo). 0029-5493/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2010.03.028

(SSI) system. The infinite element method is an attractive technique since its concept and mathematical formulations are similar to those of the FEM except for the infinite extent of the element region and shape functions. IEM has been successfully extended to various fields and there is a remarkable literature over 300 papers on the topic of infinite elements. These are mentioned in the below reference list, for instance, Chow and Smith (1981) and Rajapakse and Karasudhi (1985) for elasto-statics; Schrefler and Simoni (1987) for consolidation analysis; Zhao and Valliappan (1993a,b,c) for seepage flow; Zienkiewicz and Bettess (1983), and Astley et al. (1994) for the Helmholtz equation including hydro-dynamics and acoustics. The shape functions of IE for approximating field variables were usually derived using analytical solutions in the exterior region. This may be the reason why lots of IE formulations are available and still under development. When civil engineers analyze a structure supported or surrounded by soil medium subjected to dynamic loads, it is necessary to consider the flexibility of the soil (kinematic interaction) and the energy radiation into outer far-field region (inertial interaction). The physical interaction phenomena have to be considered in the elasto-dynamic field often called soil–structure interaction. However, limited work has been done concerning the dynamic IEs related with the SSI analysis field for a homogeneous or a layered unbounded domain. Medina and Penzien (1982) carried out the

1950

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

Several researchers also developed more effective infinite elements for the layered medium (Yun et al., 1995, 2000; Yang and Yun, 1992; Choi et al., 2001; Park et al., 1992). Yun et al. proposed a systematic procedure to obtain the shape functions of the decayfunction typed infinite elements in a multi-layered soil system in frequency or time domain, and resolved the erroneous behavior of infinite elements in high frequency range. An efficient integration scheme was also proposed for calculating the element matrices involving arbitrary number of multiple wave components with less integral points. The non-linear seismic behavior of the soil, particularly the strain-dependent shear modulus and damping ratio, was also analyzed using the dynamic IE method and equivalent linearization techniques (Choi et al., 2001).

Fig. 1. 3D soil–structure interactive system by 3D finite and infinite elements.

Fig. 2. Three kinds of infinite domains and geometric mappings.

leading work in this field, and then Rajapakse and Karasudhi (1986) suggested an infinite element which is capable of propagating multiple waves in a homogeneous half-space. Yang et al. proposed the mapped frequency-independent infinite element with single wave component in a homogeneous semi-infinite domain, and they recently analyzed the special problem such as a wave propagating field along the train moving load using the 2.5 D finite–infinite element approach (Yang et al., 1996; Yang and Hung, 2001, 2008). Table 1 Geometric mapping functions of each infinite element. x

y



HIE

(1 + )

CIE

Lj (, )xj

(1 + )

Lj ()xj j=1

Notes: N is number of nodes

N 

(1 + )

Lj (, )zj j=1

Lj (, )yj j=1

N 

(1 + )

Lj (, )yj j=1

N 

Lj (, )xj j=1

N 

N

j=1

N 

VIE

z



N

Lj ()yj

−zj (1 + )

−zj (1 + )

j=1

 = 8(HIE, VIE), = 3(CIE), when 20-node FEs are used = 4(HIE, VIE), = 2(CIE), when 8-node FEs are used

Fig. 3. Typical shape functions for an 8-node horizontal infinite element (Real parts only).

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

1951

Table 2 Displacement and shape functions for each mode. Nodal mode

Edge mode

Face mode







N/2

N

HIE

NH

Lj f1 d j1 j=1

NV N/2  

N−1 NH  

j=1

m e (2j−1)m

L(2j)

m=2

Lj f1 h1 d j11 j=1

j=1

L(2j−1) m

m=2

j=1 NH 

L(2j−1) f1 j=1 n (; ω)

m f (2j)m

m=2

NV N−1 NH  

L(2j−1) l h1 e(2j−1)m1 +

N−1 NV  

Notes: m (, , ; ω) = fm (, , ; ω) − f1 (, , ; ω),

m=2

NV N/2  

L(2j−1) j=1

N 

CIE

j=1

Internal mode

L(2j) m f (2j)m

m=2

Lj h1 d j1 j=1

NH

L(2j−1) m e(2j−1)m j=1

N 

VIE

N/2

n e (2j−1)1n

n=2

m=2

n f (2j−1)mn +

n=2

L2 m h1 f 2m1 +

m=2

NV NH  

L2 m m=2

n i 2mn

n=2

NV 

L2 f1

n f 21n

n=2

= gn (; ω) − g1 (; ω).

Meanwhile, the studies concerning the three-dimensional (3D) IE method in Cartesian coordinate system are rare, only attempted by a few researchers (Zhao and Valliappan, 1993a,b,c; Seo et al., 2007). In addition, the applicable research into a seismic response analysis with the 3D dynamic IE does not exist. The main reason is that it is very complicated and more trouble-some to construct the analytical model of 3D FE–IE system than that of two-dimensional (2D) or axisymmetric. The 3D simulation of the SSI system also includes the difficult and non-economic numerical model which needs an excellent computing system. Consequently most researchers and engineers got used to attempt the methods which the 3D SSI system was replaced with the 2D or axisymmetric simplification. Our purpose is to propose a slightly new radial-shaped dynamic IE of a decay-function type and a seismic analysis technique considering 3D SSI. We then intend to show that the proposed methodologies have a reliable enough performance and wide application for real SSI problems. The dynamic IE formulation procedure and a set of wave propagation functions superior to those in previous studies are introduced in this paper, we also discuss generalities concerning a modeling manner for an infinite soil region. The performance of the proposed dynamic IE will be shown in order to get the dynamic response of a typical compliance analysis of a surface square footing and a non-symmetric L-shaped mat concrete footing, compared with other analytical and numerical results. Seismic analysis technique using the 3D FE–IE method may be demonstrated through two practical analysis tests. The applica-

tion examples for seismic analysis are a site response analysis of well-known Hualien multi-layered site excited by all seismic wave components (primary, shear, Rayleigh waves) and a generation of a floor response spectrum at the top of an inner structure in typical nuclear power plant. The obtained dynamic response shows good agreement when compared with reference values given by other researchers and a SSI analysis package. 2. Three-dimensional dynamic infinite element 2.1. Modeling strategy and geometric mapping functions The modeling strategy for an unbounded three-dimensional domain makes use of a technique to improve another methodology defined by the previous study (Seo et al., 2007) with the geometric restriction of the truncated interface between the near-field and the far-field zones. Fig. 1 shows the geometry of the proposed infinite elements with the arbitrary or non-symmetric shaped interface. They are horizontal, corner and vertical infinite elements (HIE, CIE, VIE). The decay-function type 3D infinite element generally necessitates the following integration form for both the finite and infinite intervals:



I=

P(x, y, z)W (x, y, z; k) d˝

(1)

˝

where W (x, y, z; k) is a wave function with a complex-valued wavenumber k (=ω/c) tracing the attenuation and oscillation of dis-

Fig. 4. The effect of the IE refinements (x–y plane).

1952

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

Fig. 5. 3D seismic forces on the interface by the exterior rigid boundary method.

placement in 3D infinite space; P(x,y,z) is a polynomial of finite order and denotes the straight edge or flat plane on an arbitrarily shaped interface. It should be mapped into the 1D boundary line or 2D boundary plane expressed by Lagrangian polynomials L(n), L(n, ) in natural coordinate system. The precise relationship of geometric mappings of the 3D IE is illustrated and defined as in Fig. 2 and Table 1; where xj , yj and zj are the global coordinates at node j; N is the number of nodes for each infinite element.

2.2. Displacement fields and 3D IE formulation in frequency domain The general harmonic solution in an ideal unbounded region may be well-known and derived as resolving the Navier’s governing equation related to the special boundary conditions. In numerical practice, the boundary condition at infinite boundary  ∞ called the Sommerfeld radiation is more valid if the near field zone extends a sufficient distance from the excitation source and wave front is

Fig. 6. Flow chart of the seismic response analysis with SSI.

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

spherical or cylindrical rather than planar by imposing the stringent requirement (Eringen and Suhubi, 1975). The present study induces that linear combinations of planar wave functions propagating outward the finite domain can also ensure the regularity condition on infinity. The asymptotic behavior of wave functions is regarded as complex-valued exponential functions approximated from an analytical solution, and will be presented. The time-harmonic displacement solution in an unbounded domain can be established as U(x; ω)eiωt . The amplitude U(x; ω) is the complex displacement field in a Cartesian coordinate system x = x, y, zT . The displacement field in the infinite element region may be expressed as a combination of multiple wave components in frequency domain, and then its discretized general form is rewritten as (Yang and Yun, 1992):

U(x; ω) =

 N NW   l=1

=

 Li (x)pil (ω)

Wl (x; ω)

i=1

NW N  

Nil (x; ω)pil (ω) = N p p

(2)

i=1 l=1

1953

infinite element domain can be expressed as tensor forms: U(x; ω) =

NH N  

Njm (x; ω)pjm (ω)

in ˝HIE

Njn (x; ω)pjn (ω)

in ˝VIE

j=1 m=1

U(x; ω) =

N NV  

(3)

j=1 n=1

U(x; ω) =

NH NV N   

Njmn (x; ω)pjmn (ω)

in ˝CIE

j=1 m=1 n=1

where pjm (ω), pjn (ω) and pjmn (ω) are the generalized coordinates associated with the shape functions Njm , Njn and Njmn ; NH and NV are the respective number of the horizontal and vertical wave components. The equations above mean that displacement fields may be expressed as the simultaneous propagation of multi-wave components. To construct the dynamic stiffness matrices compatible to displacement field of the finite elements in the near field, it is required to transform and express that of each IE in terms of a nodal (d), edge (e), face (f) and internal (i) displacement in the physical coordinate (q). U(x; ω) = N q (x; ω)q(ω) = N q q

where NW is the number of total available wave components; pil is the unknown vector or generalized coordinates associated with node i on interface edge/face and lth wave component; Wl (x; ω) is a wave function representing lth wave component; The shape function Nil (x; ω) corresponds to the parameter pil . It should be noted that the nodeless parameter characterized as decay-function type IE is generated at l ≥ 2. The displacement fields for each 3D



(4)





T

where N q (x; ω) = N d , N e , N f , N i , q = d T , eT , f T , iT and Nd , Ne , Nf and Ni are the shape function matrices related to each displacement mode. Example shape functions for an 8-node horizontal IE coupled with a 20-node FE are illustrated in Fig. 3, and displacements and shape functions for each mode are expressed in Table 2. The relationship of the two coordinates p(ω) and q(ω) can

Fig. 7. Rigid square footing on layered half-space and two kinds of 1/4 FE models.

1954

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

Fig. 8. Compliance functions for a rigid square footing on a layered half-space (h/B = 2.0).

Fig. 9. L-shaped mat concrete footing on three layers ( = 0.3).

Fig. 10. Three kinds of FE–IE meshes concerning the radial IEs (h = 2B).

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

1955

Fig. 11. Frequency displacements of the L-shape footing.

be written as: p(ω) = T pq q(ω), N q (x; ω) = N p (x; ω) T pq

formed into q coordinates as: (5)

where Tpq is the constant transformation matrix that can be derived from Eqs. (2) and (4). The example for the transformation matrix of corner infinite element is shown in Appendix A. The dynamic matrix of SSI system may be obtained using the elementary matrices transformed by Tpq , to apply the effective integration scheme with less integral points (Yun et al., 1995). And then the element matrices and dynamic stiffness of infinite element are also trans-

Kqq (ω) = TTpq Kpp (ω) Tpq and Mqq (ω) = TTpq Mpp (ω) Tpq S¯ qq (ω) = (1 + i2ˇd )Kqq (ω) − ω2 Mqq (ω)

K(ω) and M(ω) are the frequency dependent stiffness and mass matrices of a infinite element; S¯ (ω) is the dynamic stiffness matrix or impedance function with hysteretic damping effect; ˇd is the hysteretic damping ratio; The integrations in the finite direction can be performed using the Gauss–Legendre quadrature, while those in the infinite direction can be performed using the Gauss–Laguerre quadrature. Finally, the discretized algebraic equation of motion for forced vibration analysis of SSI system is rearranged corresponding to the index of degree of freedom, and can be set down:



Fig. 12. A layered free field at the Hualien site and the IE-FE model for site response analysis.

(6)

S nn (ω)

S ne (ω)

S en (ω)

S ee (ω) + S¯ ee (ω)



U n (ω) U e (ω)

 =

F n (ω) 0

 (7)

where S(ω) is the dynamic stiffness matrix of the finite elements; F(ω) is the external force vector applied within the near field; subscript n stands for the degrees of freedom (DOF) of the structure and the near field soil region; e denotes the DOF along the interface  I.

1956

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

2.3. Wave function The shape functions N(x;ω) for each infinite element in Eq. (3) can be expressed in natural coordinates as: Njm (, , ; ω) = Lj (, )fm (, , ; ω)

for HIE

Njn (, , ; ω) = Lj (n, )gn (; ω)

for VIE

Njmn (, , ; ω) = Lj ()fm (, , ; ω)gn (; ω)

for CIE

(8)

In the above equations, fm (·) and gn (·) are wave propagation functions based on the complex-valued exponential functions satisfying the radiation condition at infinity. These wave functions include the multi-axes values and multiple wave components as follows: fm (, , ; ω) ∈

gn (; ω) ∈

 e



e−{ˇ+iks ()ro (,)} , e−{ˇ+ikp ()ro (,)} , S {e−{˛+ika ro (,)} }N a=1 (h)

−{ˇ+iks }

(h)

,e

−{ˇ+ikp }



S , {e− sa  , e− pa  }N a=1

(9)

 N

s } and in which the sequence examples are fm ∈ {f1 , f2 , {fa+2 }a=1 Ns gn ∈ {g1 , g2 , {ga+2 , ga+3 }a=1 }, while the number of surface wave components (NS ) are selected as 1 for a homogeneous soil medium

and 2 or 3 for layered media in later numerical analyses; and then the number of wave components are determined as NH = 2 + NS and NV = 2 + 2·NS . It can be assured that the wave functions have a unit value at (0,0,0) and the oscillating characteristic with geometric decay in the infinite directions. Variable parameters in these wave functions are calculated as below: ks,p () =

1 up 1 up low low (k + ks,p ) + (ks,p − ks,p ) 2 s,p 2

⎧⎛ ⎞2 ⎛ ⎞2 ⎫1/2 ⎪ ⎪ N N ⎨  ⎬  ⎝ ⎠ ⎝ ⎠ ro (, ) = Lj (, )xj + Lj (, )yj ⎪ ⎪ ⎩ j=1 ⎭ j=1

(10)

kS () and kP () are the body (shear and primary) wavenumbers varied to -axis in a horizontally multi-layered medium. This means a linear variation of body wave numbers in the HIE between an (h) (h) upper soil layer and a lower soil layer. Also kS and kp are the body NS wave numbers in an underlying half-space. {ka }a=1 are NS -multiple surface wave numbers with a dispersive characteristic in a multi2

(h) 2



(h) 2

layered soil medium. sa = ka − (ks ) and pa = ka2 − (kp ) are the coefficients related to the surface wave in the underly-

Fig. 13. Transfer functions H jc (ω/2 ) for each wave component at six points.

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

ing half-space in z-direction. ro (, ) denotes the radial distance to node j on any infinite element from origin as shown in Fig. 1. The ˛ and ˇ are the geometric damping coefficients taken to be frequency-independent from the location of an individual IE. ˛ for the surface wave has been chosen as 0.25 and ˇ for the body waves has been taken as 0.75 referring to the results of numerical parametric studies (Seo et al., 2007). Fig. 4 shows the typical propagating phenomena for the surface waves due to a vertical harmonic excitation at a circular disk on an elastic half-space using the radial-shaped IEs. The IE refinement exposes the interpolating effect as real wave fronts, superposing planar waves between individual infinite elements.

1957

incident to the near field soil region: f

f

f

pe (ω) = S¯ ee (ω)U e (ω) − At e (ω) f

(11)

f

where Ue (ω) and te is the displacement and traction vector along  I for the free-field soil medium subjected to the earthquake excif tation and the traction te may be separately defined according to the three kinds of seismic components shown in Fig. 5; A is the constant matrix to transform the free-field traction into the nodal forces along the interface between two fields. This matrix may be defined as (4 × 4) squared form:

 

A = Aij

⎡ "⎤ 4 NG 4 NG     ! ! =⎣ Li (m , n )Lj (m , n ) !Jmn ! Wm Wn ⎦ i=1 j=1

m=1 n=1

(12) 3. Earthquake analysis & equivalent earthquake force The exterior rigid boundary method is applicable when estimating the equivalent earthquake force along the interface in this study (Zhao and Valliappan, 1993a,b,c). Shown in Fig. 5, the force vector f pe could be calculated using the result of free-field analysis and the dynamic stiffness of the far-field region. This scheme makes the f earthquake forces by calculating the reaction force (−pe ) against the zero motion of the interface (Ue (ω) = 0), while the unit con¨ 0 (ω) = 1 is subjected to free-field system. Therefore trol motion X the earthquake inputs are regarded as traveling waves vertically

where NG is the number of Gauss integral points; (m , n ) is the

!nodal ! value of the integral point; W denotes the weighing factor; !Jmn ! Wm Wn = nAe is the projected quadrilateral area of the basis plane of a horizontal or a vertical infinite elements corresponding to the incident seismic waves in Fig. 5. Eq. (7) subjected to the unit seismic excitation in frequency domain can be recomposed as below (Wolf, 1985):



S nn (ω)

S ne (ω)

S en (ω)

S ee (ω) + S¯ ee (ω)

Fig. 14. Site responses for SV-wave excitation (GL: 0.0 m).



U n (ω) U e (ω)



 =

0 f

pe (ω)



(13)

1958

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

Eq. (13) may be used to calculate the absolute displacement U(ω) corresponding to the unit control motion in the overall frequency domain, so then the displacement vector must be equal to the transfer or frequency response function H jc (ω) at any point j on ¯ j (ω) for a structure or a near field zone. The interest displacement U ¨ c (ω) by transforming x¨ c (t) can be obtained the real control motion X as:

 ¨ 0 (ω) = H jc (ω), U j (ω) = H jc (ω)X

∞

¨ c (ω) = X

x¨ (t)e−iωt dt

−∞

¨ c (ω) ¯ j (ω) = H jc (ω)X U

(14)

The transfer function H jc (ω) has to be interpolated at sufficient ¨ c (ω) may be frequency points using its smooth features, while X very irregular. The time history functions can be finally calculated using the inverse Fourier transformation formula as following: uj (t) =

1 2



∞

¯ j (ω)eiωt dω, U −∞

¨ j (t) = − u

1 2



∞

¯ j (ω)eiωt dω ω2 U −∞

(15)

And then, the floor response spectrum can be easily estimated using the above signal response. The seismic analysis flow with SSI is tabulated in Fig. 6, however, the non-linear seismic behavior of the near or far field will be not dealt with in this study.

4. Verification examples 4.1. Square surface footing on layered soil The 3D FE–IE formulation may be verified in determining the dynamic responses of a typical example including the massless rigid square footing on two layered soil depicted in Fig. 7(a). The surface footing and the near field soil region are discretized with non-conforming brick elements and the remaining far-field soil region is modeled by proposed 3D infinite elements as shown in Fig. 7(b) and (c) refined by ellipsoidal shapes. As presenting the soil–surface footing system on layered soil, the ratio of a layer depth (h) to a half width (B) is 2.0; the ratios of the unit densities (2 /1 ) and shear velocities (cS2 /cS1 ) between two layers are 1.13 and 1.25, respectively. The dimensional ratios (D/B) are determined as 2.0–2.5 by referring to the parametric studies by other researchers (Yang and Yun, 1992; Seo et al., 2007). The dynamic response functions of rigid, surface and massless footings in the SSI system may be normalized by a non-dimensionless frequency ao (= ωB/cS1 ). The non-dimensionless compliance functions are defined as below (Wolf, 1985).

CHH (a0 ) =

G1 B x , H

CHM (a0 ) =

CVV (a0 ) =

G1 B2 x M

Fig. 15. Site responses for SV-wave excitation (GL: −12.15 m).

G1 B z , V

CMM (a0 ) =

G1 B3 y , M (16)

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

where G1 (= cS2 ) is the shear modulus of the first soil layer; 1

cS1 is the shear wave velocity of the first layer; and are the complex-valued displacement and rotation; CHH and CVV are the dimensionless horizontal and vertical compliances and CMM and CHM are the dimensionless rocking and coupling compliances. Fig. 8 shows the compliance functions of a rigid square footing for a layered half-space, compared with those of previous studies (Seo et al., 2007), Wong and Luco (1985) and the conventional SSI package SASSI (Lysmer et al., 1988). The horizontal, vertical and coupling compliances of the present results also show similar tendencies when compared with other values, while the rocking compliance of Luco has a different level. This phenomenon occurring remarkably in the rocking compliance has been shown for the analysis of other typical case such as a square footing on homogeneous medium. In Chow’s paper, it was pointed out that Luco’s results were overestimated by considering the flexible coefficient, and the static compliance (a0 = 0) for a square footing on homogeneous medium determined by Chow’s method agrees more closely with the exact static solutions. 4.2. L-shaped mat concrete footing Second numerical tests are carried out for the frequency response solutions of an L-shaped mat concrete footing using the various FE–IE meshes with the proposed IEs. The work is aimed

1959

at identifying whether a stable result is obtained regardless of the model shapes. Fig. 9 shows the geometric definition of the Lshaped foundation on three layered soil medium, subjected two translational harmonic loads at the centroid of the footing plane. The SSI system is modeled by three kinds of full scaled meshes as shown in Fig. 10. The exact dynamic solution of the arbitrary footing as like L-shaped does not exist, and therefore the present results are compared with SASSI values. Fig. 11(a) and (b) shows the horizontal and vertical frequency displacement responses. It can be found that the three sets of results for each direction are almost identical. The frequency responses for the unfavorable case such as the concave model shown in Fig. 10(c) are not different to those of the two remainder cases, and the relative difference against the mean value of the other two cases is about 0.1%.

5. Applicant seismic examples 5.1. Site response analysis The site response analysis (free-field analysis) under an earthquake excitation may become the general method to verify the earthquake force input and the modeling methodology for the farfield zone. The finite and infinite element mesh and the interest points for response evaluations are shown in Fig. 12. The material

Fig. 16. Site responses for SH-wave excitation (GL: 0.0 m).

1960

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

Fig. 17. Site responses for SH-wave excitation (GL: −12.15 m).

Table 3 Material properties for site response analysis. Soil layer

Sand 1 Sand 2 Gravel 1 Gravel 2 (half-space)

Properties Shear wave velocity (m/s)

Density (Mg/m3 )

Hysteretic damping (%)

Poison ratio

Depth (m)

133 231 317 476

1.69 1.93 2.42 2.42

2.0 2.0 2.0 2.0

0.38 0.48 0.47 0.47

2.00 3.15 7.00 –

Table 4 Relative errors for the results of the site response analyses. Control motion

A15NS (SV) A15EW (SH) A15UD (P)

Location: Relative errors of each point (%) A1

A2

B1

B2

C1

C2

0.17 0.20 0.18

0.17 0.21 0.18

0.25 0.45 0.29

0.26 0.45 0.30

0.73 0.85 0.80

0.73 0.85 0.81

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

properties of the free-field system are identified as the four layers at the Hualien site, Taiwan shown in Table 3. The input control motions measured at the site on January 20, 1994 are presented in Fig. 14(a), Fig. 16(a) and Fig. 18(a), which are the vertically traveling P, SV and SH waves subjected at the surface ground level (Ohsaki Research Institute, 1994). The transfer functions and time history accelerations at the interest points of the free-field system were obtained using the proposed seismic methodology. Fig. 13 shows the trans! analysis ! fer functions !H jc (f )! of the six points corresponding to the unit control motions. It can be seen that the transfer functions of SH and SV-wave must be identical each other, not concerning the inplane or out-plane inputs in horizontally layered system. It may be observed that the results may have a few errors in the high frequency region compared with the linear SHAKE analysis (Schnabel et al., 1991). One can overcome the problem to refine more FE–IE meshes. Using Eqs. (14) and (15), the acceleration time histories at the surface (A1, A2) and bottom (C1, C2) are depicted in Figs. 14–19. To ¨ investigate the accuracy of the results, the estimated results u(t) at the interest points were compared with the control motion x¨ (t) or ¨ SH (t) from the SHAKE the results u analysis. % % The % relative % errors in ¨ % / %u ¨ SH (t) − u(t) ¨ SH (t)% × 100, and Table 4 are defined as e = %u can be observed to have very few differences with each other. It

1961

is thereby seen that the present methodology gives accurate solutions. 5.2. Generation of floor response spectrum of nuclear power plant Floor response spectra are the important end products of seismic response study, as these are directly used in seismic qualification of a subsystem or a secondary system such as equipment, pipes and their supports for a nuclear power plant. There have been many studies on the generation of the floor response spectra for the seismic analysis/design of the nuclear power plant (Sackman and Kelly, 1980; Elaidi and Eissa, 1998). A recent tendency has been work aiming toward a re-evaluation of the seismic estimation for an existing nuclear power plant structures with consideration of SSI effects on the world. The present example is typical and same thing is included in SASSI user’s manual illustrated in Fig. 20. The schematic view of the stick model shows that the upper containment structure is modeled by the vertical beams with lumped nodal masses, while the floor slab perimeter is rigidly linked on the base concrete footing. The sectional and soil properties are listed in Tables 5 and 6, respectively. Two kinds of FE models for the SASSI analysis by the sub-structure method and FE–IE analysis by the direct method are also shown in Fig. 20(b) and (c). The dynamic infinite elements

Fig. 18. Site responses for P-wave excitation (GL: 0.0 m).

1962

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

Fig. 19. Site responses for P-wave excitation (GL: −12.15 m).

Table 5 Structural properties of a nuclear power plant. Node no.

Nodal mass (ton)

Node connectivity

Section area (m2 )

Shear area (m2 )

Second moment of area (m4 )

Inner structure 1 2 3 4 5 6 7

372.0 553.4 3873.7 1705.5 2853.1 1138.5 1260.1

2–1 3–2 4–3 5–4 6–5 7–6 Base to 7

16.7 72.5 161.7 182.1 205.3 237.8 185.8

6.5 33.4 55.7 67.8 135.6 144.9 122.6

34.5 1726.2 7767.9 11220.3 10357.2 10357.2 9494.1

9–8 10–9 11–10 12–11 13–12 14–13 15–14 16–15 17–16 18–17 Base to 18

92.0 92.0 92.0 92.0 130.1 130.1 130.1 130.1 130.1 130.1 130.1

46.5 46.5 46.5 46.5 65.0 65.0 65.0 65.0 65.0 65.0 65.0

1726.2 6904.8 12946.5 16398.9 24166.7 24166.7 24166.7 24166.7 24166.7 24166.7 24166.7

Containment building 8 86.2 9 961.6 10 1120.4 11 1369.9 12 2091.1 13 1905.1 14 1905.1 15 1905.1 16 1905.1 17 1905.1 18 2086.6

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

1963

Fig. 20. SSI analysis models of a nuclear power plant.

for far-field zone are coupled to finite elements at the position of two times of the radius of the circular footing in Fig. 20(c). A dynamic behavior of a soil–structure interactive system may be basically estimated to calculate the scattering motion of a massless rigid footing in the frequency domain. The impedance analysis is executed, and then compared with other results as like those of SASSI and Luco. Fig. 21 shows the real parts of the horizontal and rocking impedances of the circular footing. Substantially, the direct method using the FE–IE method of this study is not necessary for the impedance analysis, however, the comparative study is observed to be in good agreement between the two solutions. In order to generate a floor response spectrum for the horizontal motions at the top (Node 1) of inner structure, the seismic response

analysis of SSI system are calculated using the proposed seismic methodology in Fig. 6. The control motion is the Elcentro-NS earthquake scaled down the peak value of 0.1 g shown in Fig. 22, which is travelling for 8.8 s and equivalent to the OBE level. Using Eq. (14), the time history acceleration of the interest point is obtained in Fig. 23(a) compared with the control motion. One can observe that the response is amplified to be about twice the peak of the control motion. The floor response spectra curve at the point are also shown in Fig. 23(b), while the single objective damping ratio is selected as 2%. The present result is compared with those of SASSI and FASS analysis (Bechtel Power Corporation, 1974), and also shows good agreements on the interaction analysis of the soil–structure system with each other.

1964

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

Fig. 21. Impedance functions of the rigid disk on the site (real value).

6. Conclusions

Fig. 22. El Centro NS earthquake (0–8.8 s, t = 0.005).

This paper presented the radial-shaped 3D dynamic infinite element for the SSI analysis in a horizontally layered soil medium. The shape functions of the infinite elements were grouped, and composed of the wave functions containing multiple wave components and multi-axes valuables. To verify the effectiveness of the proposed FE–IE method, the compliance functions were calculated for a rigid square footing on layered half-space. Subsequently their practical appropriateness showing the frequency responses of an L-shaped mat concrete footing was demonstrated, and the present results are found to be in good agreement with those obtained by other researchers. The seismic response analysis using the 3D coupled FE–IE method was also executed for two practical examples as the site response analysis and the generation of floor response spectra. The equivalent earthquake forces are calculated and applied along the interfaces between the near- and far-field zones. The site response analysis of a layered soil medium excited by all travelling wave components was carried out, and the site responses at the points of interest were found to be identical to the target SHAKE solutions. The floor response spectrum at the top of the inner structure of the nuclear power plant was generated and compared with the results obtained from conventional SSI packages. The present results had good agreement on the interaction problem.

Acknowledgements

Fig. 23. Dynamic responses at the top of the inner structure with SSI (damping ratio = 2%).

Table 6 Soil properties under the nuclear power plant. Shear velocity (m/s) Unit mass (Mg/m3 ) Damping ratio Poisson ratio Depth (m)

609.6 2.08 0.05 1/3 –

The study was conducted under a MEST (Ministry of Education, Science and Technology) project of Korean Government. The authors would like to thank the Smart Infra-Structure Technology Center (SISTeC) at KAIST, and also the Research Reactor Design & Engineering Division at KAERI. Their financial supports are greatly acknowledged.

Appendix A. An example of the displacement modes and transform matrix in CIE The three-dimensional displacement fields expressed by q(ω) are entirely composed of four types of modes as mentioned in the main text. Total displacement field of CIE by p(ω) generalized coordinates may be expressed as a combination of multi-wave

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

components as: U(x; ω) =

⎧ NH NV ⎨ N    l=1 m=1

⎩

where

⎫ ⎬ Lj (x)pjlm (ω)

⎭

j=1

Flm (x; ω)

(A-1)

where N (=3), NH (≥ 3) and NV (≥ 4) are the respective number of nodes, horizontal and vertical wave components; a generalized parameter pj11 can be obtained, and then substituting it into Eq. (A-5) gives: U(x; ω) =

⎧ N & ⎨ ⎩

+

+

+

j=1

⎧ NH N ⎨  ⎩

1 F11 (x) F11 (xj )

⎭

&

F1m (xj )

F1m (x) −

Pj (x)

F11 (xj )

m=2 j=1

& Fl −

Pj (x)

l=2 m=2 j=1

F11 (x)

F11 (xj )

&

Flm (xj ) F11 (xj )

⎫ ⎬

'

Fl1 (xj )

Fl1 (x) −

Pj (x)

⎧ NH NV N ⎨  ⎩

pj11 (ω)

l=2 j=1

⎧ NV N ⎨  ⎩

⎫ ⎬

'

Pj (x)

⎭ ⎫ ⎬

' F11 (x)

pj1m (ω)

'

⎫ ⎬

F11 (x)

⎭

pjlm (ω)

⎭

In the above equation, it is noted that the first term is represented in terms of the nodal displacements, and the remainders are in terms of the nodeless variables pjlm (l ≥ 2, m ≥ 2); it is added with the nodeless variables, in that the second and third terms mean the edge modes (j = 1 and j = 2) and the face modes (j = 3); also last term becomes the face modes (j = 1 and j = 2) and internal modes (j = 3). Therefore, the displacement fields expressed by the shape functions N in q coordinates become to be equal to the third block in Table 2. The relationship between two generalized vectors p and q can be obtained. pj11 =

dj −

NH 

j=1

ejl1 −

pj1m =

NH 

 NV

ej1m −

j=1 m=2

pjlm =

NH

"

1×(NH ×NV )

,

⎡ −I · · · −I −I · · · −I ⎤ ⎢ I ⎥ ⎢ ⎥ ⎢ ⎥ .. ⎢ ⎥ . ⎢ ⎥ 0 ⎢ ⎥ ⎢ ⎥ I ⎢ ⎥ I2 = ⎢ ⎥ ⎢ ⎥ I ⎢ ⎥ ⎢ ⎥ .. ⎢ ⎥ . ⎢ ⎥ ⎢ ⎥ ⎣ 0 I ⎦ 0

⎡

I

···

,

(NH ×NV )×(NH +NV −2)

I

⎤

⎢ −I ⎥ −I ⎢ ⎥ ⎢ ⎥ . . .. . ⎢ ⎥ . ··· ⎢ ⎥ ⎢ ⎥ ⎢ −I −I ⎥ ⎢ ⎥ ⎢ −I · · · −I ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ I ⎥ ⎢ ⎥ ⎢ ⎥ .. ⎢ ⎥ . I3 = ⎢ 0 ⎥ ⎢ ⎥ ⎢ ⎥ I ⎢ ⎥ ⎢ ⎥ .. ⎢ ⎥ . ⎢ ⎥ ⎢ ⎥ ⎢ −I · · · −I ⎥ ⎢ ⎥ ⎢ ⎥ I ⎢ ⎥ ⎢ ⎥ .. ⎢ ⎥ . ⎣ ⎦ 0 I

(A-5)

(NH ×NV )×{(NH −1)×(NV −1)}

in which, I is a 3 × 3 identity matrix.

Appendix B. An alternative approach for shape function groups of IE

f jlm

NH 

T

f jlm

"

(A-3) Another method to construct the shape functions for the radial 3D infinite elements will be introduced. The shape function group may be derived from Eq. (9) in the main text, while the tensor notations are different to previous cases. The shape functions for each infinite element can be expressed in local coordinates as:

f jlm

l=2

 2

NH NV  

0 ··· 0

l=2 m=2

"

m=2

j=1 l=2 2

ej1m +

m=2

l=2

NH 2  

pjl1 =

ejl1 −

NV 



I1 = I

pjl1 (ω)

(A-2)

N 

1965

NV

f jlm

j=1 l=2 m=2

One can find that Eq. (A-3) includes some rules related to the transformation. The transformation matrix Tpq is the constant square matrix of degrees N × NH × NV independent to frequencies, as follows:

⎡

I1

0

0

I2

0

0

I3

0

0

T pq = ⎣ 0

I1

0

0

I2

0

0

I3

0

0

0

I1

0

0

I2

0

0

I3

⎢

⎤ ⎥ ⎦

(N×NH ×NV )×(N×NH ×NV )

(A-4)

Njlm (, , ; ω) = Lj (, ) fl (, , ; ω)gm (, , ; ω)

HIE

Njn (, , ; ω) = Lj (, )hn (; ω)

VIE

Njlmn (, , ; ω) =

(h) (h) Lj ()fl (, , ; ω)gm (, , ; ω)hn (; ω)

CIE (B-1)

where fl , gm and hn are wave functions for l, m = 1–NH and n = 1–NF ; these are represented by the following:

1966

J.-S. Ryu et al. / Nuclear Engineering and Design 240 (2010) 1949–1966

fl (, , ; ω) ∈ (h)

fl



(, , ; ω) ∈

gm (, , ; ω) ∈

S e−{ˇx (,)+iks ()·Cx (,)} , e−{ˇx (,)+ikp ()·Cx (,)} , {e−{˛x (,)+ika ·Cx (,)} }N a=1





(h) gm (, , ; ω) ∈

hn (; ω) ∈

(h)

e−{ˇx (,)+iks

·Cx (,)}

(h)

,e

−{ˇx (,)+ikp ·Cx (,)}

S , {e−{˛x (,)+ika ·Cx (,)} }N a=1



e



(h)

−{ˇy (,)+iks ·Cy (,)}

(h)

e−{ˇ+iks

}

(h)

,e

−{ˇ+ikp }

(h)

,e

−{ˇy (,)+ikp ·Cy (,)}

S , {e− sa  , e− pa  }N a=1

! ! !x(, )! ro (, )

=

+ *

N

!* ! ! N L (,)x ! ! j=1 j j! , 2 +*

L (,)xj j=1 j

y (, ) =



S e−{ˇy (,)+iks ()·Cy (,)} , e−{ˇy (,)+ikp ()Cy (,)} , {e−{˛y (,)+ika ·Cy (,)} }N a=1



S , {e−{˛y (,)+ika ·Cy (,)} }N a=1

in which each wave function contains two variable parameters (decay and position parameters: ˛X and ˇX , CX ) dependent to scalars X (, ) considered as direction cosines of normal for wave propagations in Cartesian coordinates. The X (, ) is defined in respective directions as:

x (, ) =



+

N

,2 1/2 ,

L (,)yj j=1 j

! ! !y(, )!

(B-3)

ro (, )

So, the parameters are easily calculated as ˛X (, ) = ˛ · X (, ) = ˛ · ˇX (, ) = ˇ · X (, ) = ˇ ·

! ! !x(, )! ro (, )

,

! ! !x(, )! ro (, )

,

CX (, ) = x(, ) · X (, ) = x(, ) ·

! ! !x(, )! ro (, )

(B-4)

The formulations with the above alternative approach may be more difficult and complex than those of the main approach in the previous section. However, one can obtain the same shape function plots and compliance solutions using infinite elements with the wave functions as Eq. (B-2). References Astley, R.J., 1983. Wave envelope and infinite elements for acoustic radiation. Int. J. Numer. Methods Flu. 3, 507–526. Astley, R.J., Macaulay, G.J., Coyette, J.P., 1994. Mapped wave envelope elements for acoustic radiation and scattering. J. Sound Vib. 170, 97–118. Bechtel Power Corporation, 1974. Computer Program CE933 (FASS): Fourier Analysis of Soil–Structure System. Bechtel Power Corporation, San Francisco, CA. Bettess, P., 1977. Infinite element. Int. J. Numer. Methods Eng. 11, 54–64. Chen, C.H., Penzien, J., 1986. Dynamic modeling of axisymmetric foundation. Earthquake Eng. Struct. Dynam. 14, 823–840. Choi, J.S., Yun, C.B., Kim, J.M., 2001. Earthquake response analysis of the Hualien soil–structure interaction system on based updated soil properties using forced vibration test data. Earthquake Eng. Struct. Dynam. 30, 1–26. Chow, Y.K., 1986. Simplified analysis of dynamic response of rigid foundations with arbitrary geometries. Earthquake Eng. Struct. Dynam. 14, 643–653. Chow, Y.K., Smith, I.M., 1981. Static and periodic infinite solid elements. Int. J. Numer. Methods Eng. 17, 503–526. Elaidi, B.M., Eissa, M.A., 1998. Soil–structure interaction in fuel handling building. Nucl. Eng. Des. 181, 145–156. Eringen, A.C., Suhubi, E.S., 1975. Elastodynamic. Academic-Press. Luco, J.E., 1974. Impedance functions for a rigid foundation on a layered medium. Nucl. Eng. Des. 31, 204–217.





(B-2)

Lysmer, J., et al., 1988. SASSI: A System for Analysis of Soil–Structure Interaction. User’s Manual. University of California, Berkeley. Medina, F., Penzien, J., 1982. Infinite elements for elastodynamics. Earthquake Eng. Struct. Dynam. 10, 699–709. Ohsaki Research Institute, 1994. Blind prediction analysis of 1/20/94 earthquake. In: Hualien LSST Meeting, Ohsaki Research Institute, Palo Alto, CA. Park, W.S., Yun, C.B., Pyun, C.K., 1992. Infinite elements for 3-dimensional wavestructure interaction problems. Eng. Struct. 14, 335–346. Rajapakse, R.K.N.D., Karasudhi, P., 1985. Elasto-static infinite elements for layered half space. J. Eng. Mech. Div. ASCE 111, 1144–1158. Rajapakse, R.K.N.D., Karasudhi, P., 1986. An efficient elastodynamic infinite element. Int. J. Solid Struct. 22, 643–657. Sackman, J.L., Kelly, J.M., 1980. Equipment response spectra for nuclear power plant systems. Nucl. Eng. Des. 57, 277–294. Schnabel, P.B., Lysmer, J., Seed, H.B., 1991. SHAKE91—A Computer Program for Earthquake Response Analysis of Horizontally Layered Sites. Earthquake Engineering Research Center. UCB, USA. Schrefler, B.A., Simoni, L., 1987. Non-isothermal consolidation of unbounded porous media using mapped infinite elements. Commun. Numer. Methods Eng. 3, 445–452. Seo, C.G., Yun, C.B., Kim, J.M., 2007. Three-dimensional frequency dependent infinite elements for soil–structure-interaction. Eng. Struct. 29, 3106–3120. Song, C.M., Wolf, J.P., 1998. The scaled boundary finite-element method: analytical solution in frequency domain. Comput. Methods Appl. Mech. Eng. 164, 249–264. Tassoulas, J.L., 1983. Elements for the numerical analysis of wave motion in layered media. Int. J. Numer. Methods Eng. 19, 1005–1032. Tzong, T.J., Penzien, J., 1986. Hybrid modeling of a single-layer halfspace system in soil–structure. Earthquake Eng. Struct. Dynam. 14, 517–530. Ungless, R.F., 1973. An infinite element. M,S,A. Thesis, University of British Columbia. White, W., Valliappan, S., Lee, I.K., 1969. Unified boundary for finite dynamic model. J. Eng. Mech. Div. ASCE 103, 160–174. Wolf, J.P., 1985. Dynamic Soil–Structure Interaction. Prentice-Hall. Wong, H.L., Luco, J.E., 1985. Tables of impedance functions for square foundations on layered media. Soil Dynam. Earthquake Eng. 4, 64–81. Yang, Y.B., Hung, H.H., 2001. A 2.5D finite/infinite element approach for modeling Visco-elastic bodies subjected to moving loads. Int. J. Numer. Methods Eng. 51, 1317–1336. Yang, Y.B., Hung, H.H., 2008. Soil vibration caused by underground moving trains using the 2.5D finite/infinite element approach. J. Geotechnol. Geoenviron. Eng. ASCE 134, 1633–1644. Yang, S.C., Yun, C.B., 1992. Axisymmetric infinite elements for soil–structure interaction analysis. Eng. Struct. 14, 361–370. Yang, Y.B., Kuo, S.R., Hung, H.H., 1996. Frequency-independent infinite elements for analyzing semi-infinite problems. Int. J. Numer. Methods Eng. 39, 3553–3569. Yun, C.B., Kim, J.M., Hyun, C.H., 1995. Axisymmetric elastodynamic infinite elements for multi-layered half-space. Int. J. Numer. Methods Eng. 38, 3723–3743. Yun, C.B., Kim, D.K., Kim, J.M., 2000. Analytical frequency-dependent infinite elements for soil–structure interaction analysis in two-dimensional medium. Eng. Struct. 22, 258–271. Zhao, C., Valliappan, S., 1993a. Transient infinite elements for seepage problems in infinite media. Int. J. Numer. Methods Eng. 17, 323–341. Zhao, C., Valliappan, S., 1993b. A dynamic infinite element for three-dimensional infinite-domain wave problems. Int. J. Numer. Methods Eng. 36, 2567–2580. Zhao, C., Valliappan, S., 1993c. An efficient wave input procedure for infinite media. Commun. Numer. Methods Eng. 9, 407–415. Zienkiewicz, O.C., Bettess, P., 1983. A novel boundary infinite element. Int. J. Numer. Methods Eng. 19, 393–404.