Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 52, 3, pp. 579-594, Warsaw 2014 EARTHQUAKE ANALYSIS OF ARCH DAMS INCLUDING THE EFFECTS OF FOUNDATION DISCONTINUITIES AND PROPER BOUNDARY CONDITIONS Adel Ferdousi, Ahmad R. Mostafa Gharabaghi Department of Civil Engineering, Sahand University of Technology, Tabriz, Iran e-mail: A ferdousi@iaut.ac.ir; mgharabaghi@sut.ac.ir Mohammad T. Ahmadi Department of Civil Engineering, Tarbiat Modarres University, Tehran, Iran e-mail: mahmadi@modares.ac.ir Mohammad R. Chenaghlou, Mehrdad Emami Tabrizi Department of Civil Engineering, Sahand University of Technology, Tabriz, Iran e-mail: mrchenaghlou@sut.ac.ir; m.emami@sut.ac.ir The interaction between dam foundation and rocks plays a significant role in the design and analysis of concrete arch dams. In the current study, comprehensive dynamic analysis of the dam-reservoir-foundation system is performed by developing a non-linear finite element program. The objective of this paper is to study the seismic behavior of a concrete arch dam taking into account the effects of foundation discontinuities and foundation inhomoge- neity on the stability. Propermodeling truncated boundary conditions at the far-end of the foundation and reservoir fluid domain as well as to correctly apply the in-situ stresses for the foundation rock represents field conditions in themodelmore realistically. In this paper, the nonlinear seismic response of the dam-reservoir-foundation system includes dam-canyon interaction, dam body contraction joint opening, discontinuities (possible sliding planes) of the foundation rock and the failure of the joined rock and concretematerials. The results for the Karun 4 dam as a case study revealed the substantial effect of modeling discontinuities and boundary conditions of the rock foundation under seismic excitation that needs to be considered in the design of new dams and can be applied to seismic safety evaluation of the existing dams. Keywords: concrete arch dam, non-homogeneous and discontinuous rock foundation, nonli- near dynamic analysis, boundary conditions 1. Introduction The seismicbehavior and safety evaluation of concrete archdamshasbeen the subject of extensi- ve research during the last years. The structural strength of aconcrete damundergroundmotion is highly dependent on the stability and strength of its abutments. In arch dam foundations, failure mechanisms are typically defined by rock discontinuities, the dam-abutment interface or strata with a lower strength, and choosing a high safetymarginin the design does not guarantee immunity. Generally, due to complex nature of rock foundation including non-homogeneousma- terials, the existing joint sets and faults and propagation of seismic waves from the far or near field as well as errors due to simplified analytical simulation and etc., the final judgment about the performance of the dam is difficult. Collapse of Malpasset Dam in France which failed in 1959 is an obvious example of lack of foundation strength. During thepastdecades, extensive research invariousfields related to theanalysis anddesign of concrete damshas been conducted, but theneed formore accuratemodeling of abutments in a 580 A. Ferdousi et al. coupling systemwith considering themass,flexibility andnon-homogeneity ofdiscontinuous rock foundation is still being felt. In the present study, a numerical program for nonlinear dynamic analysis of concrete arch dams is developed in FORTRAN. For this purpose, Karun 4 dam is considered as a case study. The main features of this study are related with considering the correct in-situ stress for discontinuous bedrock and choosing a proper boundary condition for far-end which has a direct effect on the accuracy and precision of analytical results. 2. A review of research areas and solving methods Themain factors that influence significantly the three-dimensional nonlinear dynamic analysis of archdamsare identified: 1)dam-reservoir interaction anddistributionof hydrodynamicpressure, 2) reservoir-foundation interaction and effects of reservoir bottom sediments, 3) dam-foundation interaction and role of non-homogeneous anddiscontinuities in bedrock, 4) boundary conditions, 5) nonuniform input of free-fieldmotions, and 6) nonlinear behavior of the quasi-brittlematerial of concrete and joined rock as well asdiscontinuities in the dammonoliths and peripheral joints of dam body. Fundamentals and analytical methods of all the abovementioned is outside the scope of this article and just a brief review of the main issues related to research are presented here. Simple and primary models in earthquake analysis of dams are the added mass approach of Westergaard for fluid-structure interaction. Westergaard’s analytical solution neglected dam flexibility and water compressibility, so several research works developed advanced numerical methods based on finite elements, boundary elements or both of them to model dynamic dam- reservoir interactions. Two common finite element approaches in the fluid domain are Eulerian- andLagrangian-based formulations (Bouaanani andLu, 2009). The former approach is knownas pressure- or potential-based formulations where fluid pressure or velocity potentials are selected as state variables. Thus, special contact algorithms are required at the fluid-structure interface. The Lagrangian approach in the fluid domain is an extension of the solid finite element formula- tion with nodal displacements as degrees of freedom. Therefore, the fluid domain is formulated using the same shape functions as structural elements and compatibility at the fluid-structure interface is automatically satisfied. In this case, fluid elements are characterized by a volumetric elasticmodulus equal to the fluidbulkmodulus (or fluid compressibility), with a negligible shear resistance and Poisson’s constant being nearly 0.5 to simulate fluid flow realistically. The next important point is the interaction between impoundedwater and rock foundation. The partial absorption of pressurewaves at sediment layers of reservoir bottom and lateral sides may significantly affect the magnitude of hydrodynamic forces as the response of dams to the ground motion. In general, assuming an absorptive reservoir boundary leads to a realistically less response for damswith impoundedwater due to groundmotions (Mirzabozorg et al., 2003). In the current study, the Lagrangian approach is used formodeling of the fluid and sediment domains. Also, the interface elements with low shear stiffness aremodeled for common surfaces of fluid and solid elements. The arch dam-foundation interaction is related to the bedrock flexibility, changes of physical properties of rock foundation, existence of joints and faults in the rock, geometry of the dam body,water leakage and uplift pressure, etc. Differentmodelswere used for foundationmodeling (such as mass less/massed and rigid/flexible foundation) in order to determine the seismic re- sponse of concrete arch dams. It has been proven that in a nonlinear dynamic analysis including the dam-foundation interaction, the foundationmass, flexibility and radiation dampingis impor- tant. Also, in order to model the highly complex behavior of joined rock masses, the strength and deformability of joined rock masses should be expressed as a function of joint orientation, joint size, and joint frequency.Moreover, it is not possible to represent each joint individually in Earthquake analysis of arch dams including the effects... 581 a constitutive model. Therefore, it is necessary to use a simple technique such as the equivalent continuum method which can capture reasonably the behavior of the joined rock mass. In the developed program, the finite element method recognizes that the foundation rock should act both as: 1) nonlinear solid element for modeling the joined rock as an equivalent continuum whose properties represent material properties of the joined rock, and 2) nonlinear interface element has been used to account for surface roughness of discontinuities. Therefore, this paper deals with the mass, flexibility and non-homogeneity of the foundation rock, in addition to the discontinuity due to the master joints and faults. Moreover, the shape and size of the foundation model must be properly selected. Using the finite element procedure, a spherical and cylindrical system in the lower and upper half model is employed respectively. A right way to determine the size of the foundation model is based on the ratio of deformation modulus of foundation Ef to the elastic modulus of the concrete dam Ec. For a flexible rock foundationwith Ef/Ec less than 1/2, the foundationmodel should extend at least twice of the damheight in all directions (Federal EnergyRegulatoryCommission Division of Dam Safety and Inspections,Washington DC, 1999). The definition of suitable boundary conditions related to its surrounded domain is another important part of the modeling procedure. For the present study, governing equations and re- lated boundary conditions consist of free surface water and far-end truncated boundaries of the reservoir and rock foundation. By neglecting the effects of gravity waves, a zero-pressure boun- dary condition is prescribed at the horizontal free surface of water (negligible surface tension). This assumption of no-surface wave is a common assumption in the analysis of concrete dams, particularly for deep reservoirs. To simulate the cut-off boundaries of the foundation, numerical models must include nor- mal and tangential energy absorption elements to model wave propagation problems. Several sophisticated waves transmitting boundaries have been developed bymany researchers (Khalili et al., 1997; Roesset and Ettouney, 1977; Bettes, 1977; Felippa and Park, 1980; Wolf et al., 1985-1996; Valiappan and Zhao, 1992; Zhao and Liu, 2003; Jahromi et al., 2007, 2009; Degroote et al., 2010). The viscous boundaries proposed by Lysmer and Kuhlemeyer (1969) consist of independent dashpots attached to the boundary have been applied in the rock foundation in all directions. The boundary element method has also been widely used (Zienkiewicz et al., 1977; Estorff and Kausel, 1989). However, it is not very easy to derive the fundamental solutions for many types of problems by the BEM. The reference (Beer et al., 2008) gives the integrated description of the BEM. Similarly, a radiating condition (such as Sommerfeld, 1949; Sharan, 1987; Maity and Bhat- tacharyya, 1999; Küçükarslan, 2004) can be applied to the truncated boundary of the reservoir in dynamic analysis of structures. However, the effect of the radiating condition on the solutions is generally negligible if the reservoir length is taken as three or more times of the dam height. In this paper, the truncated boundary condition for the reservoir and foundation domains includes the interface elements with normal and shear stiffness, and with the option to apply viscous boundaries. Using these boundary conditions, does not prevent the rock sliding along fault zones during seismic loading and it can provide the non-uniform excitation more realisti- cally. At the foundation boundaries, the displacement history is used instead of the acceleration history as the seismic input in dynamic analysis. In evaluating the seismic performance ofmega structures such as concrete dams, it is evident that themanner inwhichgroundmotions excite thedam-reservoir-foundation system is ofmajor importance. Variations in the ground motion arise mainly from three sources, e.g. the wave passage effect, the incoherence effect and the site response effect(Alves and Hall, 2006). In this paper, thewave incoherence effects due to differentmass andboundary conditions of rock blocks in the modeled domain were considered as non-uniform input sources at the dam foundation surface more realistically. 582 A. Ferdousi et al. Furthermore, peripheral and vertical construction joints of the dam body as well as discon- tinuities in the foundation may slip or open, and the concrete or joined bedrock materials may crack and crash during intense earthquake motions; thus, nonlinear dynamic analyses must be recognized by appropriate material models, which will be described later. 3. Numerical modeling and analysis procedure for the dam-reservoir-foundation system In this study, the dam-reservoir-foundation system is modeled by an assemblage of solid and interface elements, as shown in Fig. 1. The isoparametric 8-node solid elements with 2×2×2 Gauss integration are used for all domains including the dam body, reservoir, sediment and rock abutments. The output values of stresses, strains, invariants, principal values, etc., are the simple numerical average of Gaussian point values in the centerpoint. Also, eight-node interface elements are used in common surfaces of domains, such as the contraction and perimetral joints in the dam body and discontinuities of rock abutments (interface elements with normal and shear stiffness are set between contacts blocks of the foundation for simulating the behavior of major joints). The interface elements are placed between continuum (solid) elements. Figure 2 shows a detail of the 3D dam-reservoir contact including the interface elements. The elements formulations support geometric and material nonlinear analysis. Assuming that non-linearities are limited to the concrete dam, rock blocks, contraction joints of the dambodyand rock discon- tinuities, elements stiffness need to beupdated at each iteration.Notably, in the programing, the capabilities of several open source programs that are interested in the analysis of concrete dams such as ADAP-88 (Mojtahedi et al., 1992), EAGD-SLIDE (Chavez and Fenves, 1994), EACD- 3D-96 (Tan and Chopra, 1996) were investigated, and useful subroutines and subprograms of the finite elements were used (Smith and Griffiths, 2004; Bathe, 1996). Fig. 1. Detailed view of solid-interface elements Fig. 2. Interface element sketch at the reservoir-dam contact The governing equations of motion for 3D nonlinear dynamic analysis of the coupling sys- tem (subjected to earthquake loads) were discretized by the Newmarkmethod. The discretized nonlinear dynamic equation of motion was given byBathe (1996) and Zienkiewicz et al. (1972) Earthquake analysis of arch dams including the effects... 583 ( KT + γ β∆t CT + 1 β∆t2 M ) Ün+1 =Rn+1+C [ γ β∆t Un+ (γ β −1 ) U̇n+∆t ( γ 2β −1 ) Ün ] +M [ 1 β∆t2 Un+ 1 β∆t U̇n+ ( 1 2β −1 ) Ün ] (3.1) where M is themass matrix, C is the dampingmatrix, K is the stiffness matrix and R is the nodal vector of external forces. Ü, U̇ and U are the acceleration, velocity and displacement vectors, respectively, at the (n+1)-th time step.Thedampingmatrix is determinedbasedon the well-known Rayleigh damping (C = aK+ bM where the parameters a and b are pre-defined constants). Also, KT is the tangent stiffness matrix and the updated damping matrix CT reduces at the same time as the stiffness reduces. The full Newton-Raphson iteration scheme can be employed to resolve the nonlinearity. The parameters β and γ determine the stability and accuracy characteristics of the algorithm. The constant acceleration method is obtained when β = 1/4 and γ = 1/2. The method is implicit, unconditionally stable and second-order accurate. There are different finite element formulations to handle the fluid-structure interaction pro- blem that can be classified into the displacement formulation, the pressure formulation and velo- city potential formulation.Thefinite element formulation of the fluid is based on theLagrangian approach in which the fluid strains are calculated from the strain-displacement equations.The pressure volume relationship for a linear fluid is expressed by p=λεv (3.2) where p, λ and εv are pressure that is equal to the mean stress, the bulk modulus and the volumetric strains of the fluid, respectively. The volumetric strain can be expressed in terms of Cartesian displacement components as follows εv = ∂ux ∂x + ∂uy ∂y + ∂uz ∂z (3.3) where ux, uy and uz are displacement components related to the axes x, y and z, respectively. The major disadvantage of the Lagrangian approach is the generation of spurious circulation modes, due to the zero shearmodulus.This causes somenumerical problemswhich, in turn, lead to false, zero-energy deformation modes for the fluid elements and/or fluid-structure system. There are different methods for eliminating these zero energy modes that have been used by a large number of researchers (Chopra et al., 1969; Shugar and Katona 1975; Hamdi et al., 1978; Zienkiewicz and Bettess, 1978; Akkaş et al., 1979; Wilson and Khalvati 1983; Olson and Bathe, 1983; and Doǧangün et al., 1996; etc.). One approach is based on assuming that the shear modulus of the fluid is numerically very small. This approach has been adopted by this study. Another approach admits inviscid behavior and uses the implication that the fluid must be irrotational in nature. This behavior can be enforced by the use of a penalty function. Combinations of the Mohr-Coulomb yield function with a tension cut-off (i.e. the Modified Mohr-Coulomb model suggested by Paul (1961)) are used for both concrete and joined rock materials in the system modeling. In literature, the Mohr-Coulomb criterion (1882-1900) is widely used. Also, Rankine’s maximum stress criterion (Rankine criterion) crack model is used to simulate crack formation under tensile conditions. TheMohr-Coulomb criterion is based on Coulomb’s equation (1773). If σ11 ­σ22 ­σ33 are principal stresses, we can write this criterion as 1 2 (σ11−σ33)cosϕ+ [1 2 (σ11+σ33)+ σ11−σ33 2 sinϕ ] tanϕ−c=0 (3.4) where ϕ and c are the internal friction angle and cohesion, respectively. The 3D failu- re surface of the Mohr-Coulomb criterion can be expressed in terms of stress invariants 584 A. Ferdousi et al. (I1 = σii = σ11 + σ22 + σ33, J2 = sijsij/2, J3 = sijsjkski/3 with sij = σij − I1δij/3, δij – Kronecker delta) f1(I1,J2,J3)= I1 3 sinϕ+ √ J2 sin ( θ+ π 3 ) + √ J2 3 cos ( θ+ π 3 ) sinϕ− ccosϕ=0 (3.5) where theLode angle is θ= 1 3 cos−1 ( 3 √ 3 2 J3 √ J32 ) For the maximum-tension-stress yield function, we have (Rankine criterion) σ11 = f ′ t σ22 = f ′ t σ33 = f ′ t (3.6) with the tension strength f ′t. This criterion can be fully described by the following equation f2(I1,J2,J3)= 2 √ 3J2cosθ+ I1−3f ′t =0 0¬ θ¬ π 3 (3.7) At every time step, the programwill check dam and rock foundation solid elements with the Modified Mohr-Coulomb criterion. Figure 3 shows the failure envelope under these combined criteria. Fig. 3. ModifiedMohr-Coulombmaterial model for concrete and joined rock solid element Different models have been developed to represent the contact between two surfaces. The cohesive law can be expressed such that the local traction t across the interface is taken as a function of the displacement jump δ across the cohesive surfaces. Defining a free energy density per unit undeformed area Φ such that the traction acting at the interface is given by the exponential function (Needleman andOrtiz, 1991) t= ∂Φ ∂δ = eσc δ δc e − δ δc if δ= δmax and δ̇­ 0 (3.8) where δc denotes the value of δ at the peak traction (tmax =σc), and we have δ= ‖δ‖= √ δ2n+κ 2δ2s δn = δ ·n δs = ‖δ− δnn‖= √ δ2S1+ δ 2 S2 t= tn+ts = ∂Φ ∂δn n+ ∂Φ ∂δs ∂δs ∂δs = t δ (δnn+βκ 2 δs) Φ= eσcδc [ 1− ( 1+ δ δc ) e − δ δc ] (3.9) Earthquake analysis of arch dams including the effects... 585 When the contact surfaces undergo a combination of shear deformation andnormal compres- sion, the effective separation δ is computed only from the shear components.Also, undernormal compression, the cohesive material behaves as a linear spring. The weighting coefficient κ defi- nes the ratio between the shear and normal critical tractions. Formore details, see the reference (Ruiz et al., 2001). The value of interface stiffness will depend on the roughness of contact surfaces as well as the relevant properties of the filling material andmoisture. For an initially closed interface, the normal stiffness Kn and tangential stiffness Ks are set to have a high value. These values can be estimated from the lowest Young’s modulus and shear modulus of the adhesive domain around the contact, according to following relations Kn =m1 EAe le Ks =m2 GAe le (3.10) where mi is a factor that controls contact properties, usually between 0.01 and 100 (only in penetration), E and G are the smaller elastic and shearmodulus when considering the contact between two different materials, le is the characteristic thickness of the adjacent solid elements perpendicular to the interface, and Ae is the interface element surface. The far-end boundary conditions are employed as a viscous boundary (Mirzabozorg et al., 2003). Fig. 4. FEmodel of Pine Flat dam on rigid base andsubjected to Taft-1952 groundmotion In order to verify the accuracy and validity of the finite element modeling and developed computer code, the tallest monolith with unit width of well-known Pine Flat dam, a concrete gravity dam in California, which is 122m high, has been selected. A water depth of 116m is consideredas the full reservoir condition.Thegeometry andFEmodel ofPineFlat dammonolith with unit width is shown in Fig. 4. The properties of applied material in the modeling are: Ec =22.75GPa, ν =0.2 and ρ=25kN/m 3. For nonlinearmaterial analysis, the tensile strength of the concrete is taken to be 3.36MPa which is 12% of its compressive strength. The dynamic tensile strength shall be equivalent to the direct tensile strength multiplied by a factor of 1.50 (Raphael, 1984; Cannon, 1991). The analysis results consist of the weight of the dam, the static pressure of the impoundedwater and the earthquake excitation of thehorizontal x component of Taft-1952 Lincoln California groundmotion with scaling to PGA=0.4g. Proportional damping in the dam provides a critical damping ratio of 5% in the fundamental vibration mode of the dam. Figure 5 shows the comparison between the results of crest displacement for the first 15 seconds of excitation obtained from the present developmental program and commercially available ANSYS program in the case of a fixed base and the dam body with a linear elastic material. For transient structural analysis in ANSYS, a corresponding 8-node 3-D concrete 586 A. Ferdousi et al. element SOLID65 is selected. Also, Fig. 6 showsthe analysis resultsobtained from developed programwith considering the effects of the material nonlinearity and base sliding. Fig. 5. Comparison of displacement results of the Pine flat dam crest due to horizontal component of Taft-1952 with ANSYS program results Fig. 6. Comparison of horizontal displacement for the nonlinear behavior of materials and the possibility of slip at the base of dam After verification of the developed program, in the next section, the highest concrete dam in Iran (Karun-4) is considered for a case study and investigation of the influence ofmore accurate modeling of the foundation in the dynamic responseof dam. 4. Case study analysis results Karun-4damisadoublecurvaturearchdamontheKarunRiver in theprovinceofChaharmahal- e Bakhtiari, Iran. The main objective of this project is power generation and flood control. The whole crest length of the dam (440m) is divided by 20 contraction joints. Some geo- metric characteristics of the dam are: crest elevation=1032.0m, maximum height=230.0m, crest thickness=7.0m, base thickness=37.0m, maximum thickness=50.5m and concrete vo- lume=1675000m3. Figure 7a shows the dam structure with its appurtenance. The modulus of elasticity, Poisson’s ratio and unit weight of the concrete are taken as 23.6GPa, 0.2 and 24kN/m3, respectively. The tensile strength of the concrete is taken to be 2.75MPa, and dyna- micmagnification factors of 1.5, 1.3 and 1.25 are considered for themodulus of elasticity, tensile and compressive strengths, respectively. The damping ratio of the dam and the foundation equ- als 5%.Based on thegeotechnical investigations and studies, thefinal classification and estimated geomechanical parameters of the rock mass are mostly composed of: unit weight=25kN/m3, deformation modulus=11.0GPa, Poison’s ratio=0.25, friction angle=42◦, cohesion=0.5MPa and allowable bearing capacity=9-14 (→ 12)MPa. The nonlinearity in the finite element ana- lysis was incorporated in form of material nonlinearity of the equivalent rock with uniaxial Earthquake analysis of arch dams including the effects... 587 compressive and tensile strength of 12 and 1.2MPa, respectively. The FE model of the founda- tion extends 2.5 times of the dam height in all directions. Fig. 7. 3D finite elementmodel of: (a) Karun 4 arch dam and (b) the rock foundation divided into six blocks The main idea of the study is to investigate the effects of non-homogeneous characteristics of the rock foundation on the seismic performance of arch dams. Themodulus of deformation of the abutments and foundation is an important element in analyzing the performance of the dam since the flexibility of the foundation directly affects the stresses in the dam.On the other hand, for a discontinuous foundation, the effect of a large faulted zone on themodulus of deformation of the foundation must be taken into consideration. A large change in the modulus may result in formation of concentrated stresses in the concrete of the dam. For this purpose and regarding the geometry of discontinuities in each abutment, as primary analysis results (based on site investigations reported byMahab Ghodss Consulting Engineers), “F4-a & F6-a” and “MJ67-c, MJ28” are defined as critical discontinuities in the left and right abutment, respectively. The characteristics of the critical discontinuities are presented inTable 1.These platesmake six large blocks, as shown in Fig. 7b. Table 1.Geomechanical parameters of the critical discontinuities in the left abutment Geometrical Discontinuities specification F4-a F6-a MJ28New MJ67-c Dip direction 52◦ 1◦ 349◦ Us/: 15◦ Ds/:030-070 Dip 30◦ 41◦ 35◦ Us/: 35◦ Ds/: 30◦ Leakage condition Wet Wet Wet NA Geomechanical condition Rock fractured – calcium filling thickness 10-15cm, 2m displaced, planar-smooth Fractured zone, Fe gravel clay filling 10-30cm, 2-8m displaced, planar, rough, smooth Rock fractured, filling 2cm, planar, rough NA Thefluidbodyof the reservoir has beenmodeled using theLagrangian approach bymodified elastic elements to solve hydrostatic and hydrodynamic pressures acting on the dam and a maximum reduction factor of 0.006 applied to the shear modulus to approximate simulation of inviscid flow.Thewater in the reservoir has a constant depthof 155m,mass density=10kN/m3, bulk modulus=2131MPa and Poisson’s ratio=0.495 for nearly incompressible nature. For the upstream/downstream sediments with the assumptio n of about 60/30 meters depth, the mass 588 A. Ferdousi et al. density is ρs = 13.6kN/m 3 and the bulk modulus Bs = ρsc 2 s = 3312MPa, where the sound speed profile is estimated from physical sediment properties using Biot’s theory and assumes cs = 1560m/s. In all contact surfaces of the fluid-solid, the sliding layers are used to enforce the slip condition in order to decouple the tangential displacement components. The developed finite element mesh of the reservoir and sediments is shown in Fig. 8a. Fig. 8. (a) Finite element mesh of the reservoir and sediment elements, (b) schematic view of interface elements The interface elements are used formodeling the rock discontinuities, vertical joints between cantilevers and the peripheral of the dam in connectionwith the canyon rock, aswell as contacts between the reservoir and surrounding domains with negligible shear stiffness. At the truncated boundaries of the reservoir and rock foundation, the appropriatemethods such as the boundary element and interface element methods are available in the developed numerical programwhich can be applied at the boundaries. Bouaanani and Lu (2009) gave an integrated description of BEM. In the case of dynamic analysis, using the interface elements can provide a high analysis efficiency and give a good estimation. Therefore, the interface elements have been used in the far-end of the infinite domain for present modeling of the analytical system (called “Moving B.C.” later). The properties of several interface elements are presented in Table 2. The coupled systemmodel includes 11764 nodes and9348 finite elements. The elementmeshes of the interface elements (without discontinuity elements of the foundation) are shown in Fig. 8b. Table 2. Interface elements parameters Interface stiffness Position of contact surfaces Normal direction Tangential direction [N/mm3] [N/mm3] Contraction joints in the dam 3.0 ·109 1.5 ·109 Peripheral joints at the dam-foundation 4.0 ·109 2.0 ·109 Discontinuities in rock masses 1.0 ·109 0.8 ·109 Far-end boundaries of rock foundation 6.0 ·109 3.0 ·109 The loads applied to the model consist of static and dynamic loading. The static loads are dead weight and hydrostatic pressure at a normal level of water and sediment weight. The effects of temperature, tail water load and uplift are neglected, however in the complete safety evaluation these loads should be taken into account. Information about the in-situ stresses of the rock field is a fundamental parameter for the dam-foundation analysis and has a direct effect on the dynamic design of such a coupled system. The need for understanding of the in-situ stresses in the rock foundation has been recognized by geologists and engineers for a long time and several methods for their measurement have been proposed since the early 1930’s. Generally, the state of stress in the dam foundation is a result of Earthquake analysis of arch dams including the effects... 589 superposition of several stress fields, first and foremost being that due to the gravitational force. The in-situ stress in a rockmass is simply equal to theweight of the overlyingmaterial, therefore the discontinuities control the magnitude and direction of this stress field. Most often, a very variable in-situ stress field exists in a fractured rock mass. In this study, firstly the static load of discontinuous rock weight has been applied to investigate the in-situ stress. For this loading case, the dam body should remain free of stress because of canyon deformation. To overcome this problem, the numerical program has the ability to change the material properties such as Young’smodulusandPoisson’s ratio in loading steps.Therefore, in thepre-loading step,Young’s modulus of the dam body and the region of rock abutments near the dam (the radius of about 15m) gradually decreases, and Poisson’s ratio increases to 0.49. This causes the stresses of the dam body to be negligible during the loading step of the rock weight. Such a process demands a suitable selection of material properties and a sequence of loading steps which would include a sensitive loading pattern. In the next dummy load step, material properties gradually change to real values. The first 40 seconds of the three components of the Taft Lincoln School records of the 1952 Kern County, California earthquake are used as input ground motion in the dynamic loading. The peak ground acceleration of x,y (horizontal components), and z (vertical component) are 0.156g, 0.178g and 0.108g, respectively, and to develop seismic hazard study for Karun 4 dam site, the earthquake time histories are scaled to the maximum credible level at the middle height of the canyon (PGAhor = 0.49g, PGAver = 0.26g). A time step of 0.01s ischosen for the analysis. The displacement time histories for the three components of Taft earthquake are shown in Fig. 9. Fig. 9. Three displacement components of KernCounty, California earthquake of July 21, 1952 recorded at the Taft Lincoln School Tunnel Fig. 10. Maximum principal stress contours with deformed scaleof 20 at time 20 second of excitation for the cases of (a) C1, and (b) C2 (upstream views [Kgf/m2]) In order to present the effect of foundation interaction on the seismic response of the dam- reservoir system, several cases of massive foundations are chosenby considering geometric and material nonlinearity as follows (under Taft earthquake): • C1: Continuous rock foundation – fixed boundary condition (without interface elements between the rock blocks and on the truncated boundary), 590 A. Ferdousi et al. • C2: Discontinuous rock foundation –moving boundary condition (with interface elements between the rock blocks and on the truncated boundary), • C3: Condition C2 with applying a reduction factor of 10% for the deformation modulus and allowable bearing capacity of rock blocks RB1, RB2, RB3 andRB4 (demonstrated in Fig. 7b) materials, • C4: Condition C2 with applying a reduction factor of 10% for the deformation modulus and allowable bearing capacity of rock blocks RB5 and RB6 (demonstrated in Fig. 7b) materials. For another comparison, the first 40 seconds of the three components of the El Centro, Imperial Valley – 1940 records are used. Figure 11 shows the displacement time histories for the ElCentro earthquake. Two cases considering geometric andmaterial nonlinearity have been investigated as follows: • C5: Continuous rock foundation – fixed boundary condition, • C6: Discontinuous rock foundation – moving boundary condition. Fig. 11. Three displacement components of the El Centro, Imperial Valley – 1940 recorded In order to design an economical and safe concrete arch dam, the dynamic analysis procedure should consider the dam-reservoir interaction, realistic modeling of dam-foundation interaction, reservoir bottom sediment andboundary conditions, aswell as nonlinearity effects of contraction joints, discontinuities of the foundation and failure of the joined rock and concrete materials, as shown in this paper. In the current study, the influences of foundation discontinuities and inhomogeneity on the response of the dam is highlighted. Figures 12-15 show the crest displacement of the crown cantilever in the upstream-downstreamdirection for several cases. As can be seen, themodeling of the foundation withmore details have a very important role in the coupling system analysis. Also, using the interface elements with an appropriate characteristic on the far-end boundary and major fault zones of the foundation changed the seismic response of the dam significantly. The foundationmodel consists of blocks, divided into twomaterial regions in the upstream and downstream parts that have been chosen based on the effects of impounded water in the rock material characteristic and site investigation. As can be seen, the foundation flexibility with a reduction factor of 10% for the deformation modulus and allowable bearing capacity of the domain parts has significantly affected the dam response. It should be noted that this boundary condition and modeling discontinuity in the bedrock are critical for an actual response of the dam rather than application of a non-homogeneous material of the foundation, as compared in Figs. 13 and 14. A summary of the upstream-downstream and vertical displacements for all analyses are provided for comparison in Table 3. Earthquake analysis of arch dams including the effects... 591 Fig. 12. Comparison of upstream/downstream crest displacement of Karun 4 dam under Taft earthquake for C1 and C2 cases Fig. 13. Comparison of upstream/downstream crest displacement of Karun 4 dam under Taft earthquake for C1, C3 and C4 cases Fig. 14. Comparison of upstream/downstream crest displacement of Karun 4 dam under Taft earthquake for C2, C3 and C4 cases Fig. 15. Comparison of upstream/downstream crest displacement of Karun 4 dam under El Centro earthquake for C5 and C6 cases 592 A. Ferdousi et al. Table 3. Upstream-downstream and vertical displacement comparison of crest at crown canti- lever Cases index Seismic analysis results under Taft earthquake Max. upstream Time Max. downstream Time Max./min. Time [s] at displacement [m] [s] displacement [m] [s] vertical dis. [m] max/min dis. C1 0.681 32.40 0.529 31.60 0.174/-0.441 31.71/32.44 C2 0.567 17.19 0.601 16.14 0.178/-0.248 14.43/36.07 C3 0.573 38.60 0.514 16.30 0.212/-0.284 14.52/38.70 C4 0.603 17.34 0.561 39.39 0.164/-0.266 14.47/36.27 C5 0.517 22.58 0.674 17.67 0.115/-0.392 16.24/36.39 C6 0.592 38.25 0.458 39.22 0.131/-0.324 15.68/34.48 5. Conclusions In this paper, an early analytical procedure to evaluate the seismic response of an arch dam, consideringvarious effects of thedam-foundation interaction in the timedomainhasbeen investi- gated by implementing the effects of inertia andflexibility of the foundation rock in the analysis. For this purpose, the far-end boundary condition andmajor discontinuity of the foundation are modeled by interface elements. Asmentioned in the previous section, seven cases have been considered: estimation of nonli- nearmodels with flexibility, non-homogeneous and discontinuities effects of the foundation rock on the dam response. In each case, the time history response of crest displacement of the crown cantilever in the upstream-downstream direction have been obtained. The results obtained from this study show that application of the displacement time history to the model with material non-homogeneity in the foundation is an important factor in the seismic response of the arch dam.However, the including of the interface elements on the far-end boundary of the foundation and foundation discontinuities can significantly affect the response of the dam compared with the application of the non-homogeneousmaterial of the foundation. Acknowledgment The authors acknowledge the documentary support of Karun 4 Project fromMahabGhodss Consul- ting Engineers, Iran. References 1. AkkaşN.,AkayH.U.,Yilmaz Ç., 1979,Applicabilityofgeneral-purposefinite elementprograms in solid-fluid interaction problems,Computers and Structures, 10, 773-783 2. Alves S.W., Hall J.F., 2006,Generations of spatially non-uniform groundmotion for nonlinear analysis of a concrete arch dam,Earthquake Engineering and Structural Dynamics, 35, 1339-1357 3. Bathe K.J., 1996,Finite Element Procedures, Prentic Hall Prerss 4. BeerG., Smith I.,DuenserC., 2008,TheBoundary ElementMethod with Programming, Sprin- gerWien Prerss, USA 5. Bettes P., 1977, Infinite elements, International Journal NumericalMethods for Engineering, 11, 53-64 6. Bouaanani N., Lu F.Y., 2009, Assessment of potential-based fluid finite elements for seismic analysis of dam-reservoir systems,Computers and Structures, 87, 206-224 7. CannonR.W., 1991,Tensile strengthof roller compactedconcrete,U.S.ArmyCorps of Engineers, North Pacific Earthquake analysis of arch dams including the effects... 593 8. Chavez J.W., Fenves G.L., 1994, EAGD-SLIDE: A computer program for the earthquake ana- lysis of concrete gravity dams including base sliding, Department of Civil and Environmental Engineering, Report No. UCB/SEMM-94/02, University of California, Berkeley California 9. Chopra A.K., Wilson E.L., Farhoomand I., 1969, Earthquake analysis of reservoir-dam sys- tems,Proceedings of the Fourth World Conference on Earthquake Engineering, Santiago, Chile, 2, 1-10 10. Degroote J., Haelterman R., Annerel S., Bruggeman P., Vierendeels J., 2010, Per- formance of partitioned procedure in fluid structure interaction,Computers and Structures, 88, 7, 446-457 11. Doǧangün A., Durmus A., Ayvaz Y., 1996, Static and dynamic analysis of rectangular tanks by using the Lagrangian fluid finite element,Computers and Structures, 59, 547-552 12. Estorff O. Von, Kausel E., 1989, Coupling of boundary and finite elements for soil-structure interaction problems,Earthquake Engineering and Structural Dynamics, 18, 1065-1075 13. Federal Energy Regulatory Commission Division of Dam Safety and Inspections Washington, DC 20426, 1999, Engineering Guidelines for the Evaluation of Hydropower Projects. Chapter 11 – ArchDams 14. Felippa C.A., Park K.C., 1980, Staggered transient analysis procedures for coupledmechanical systems: formulation,Computer Methods in Applied Mechanics and Engineering, 24, 61-111 15. Hamdi M.A., Ousset Y., Verchery G., 1978, A displacement method for the analysis of vi- bration of coupled fluid-structure systems, International Journal for Numerical Methods in Engi- neering, 13, 139-150 16. Jahromi H.Z., Izzuddin B.A., Zdravkovic L., 2007, Partitioned analysis of nonlinear soil- structure interaction using iterative coupling, Interaction and Multiscale Mechanics, 1, 1, 33-51 17. Jahromi H.Z., Izzuddin B.A., Zdravkovic L., 2009, A domain decomposition approach for coupledmodeling nonlinear soilstructure interaction,ComputerMethods in AppliedMechanics and Engineering, 198, 33, 2738-2749 18. Khalili N., Valiappan S., Tabatabaie Yazdi J., Yazdchi M., 1997, 1D infinite element for dynamic problems in saturatedmedia,Communications in NumericalMethods in Engineering, 13, 727-738 19. KüçükarslanS., 2004,Dynamic analysis of dam-reservoir-foundationinteraction in timedomain, Computational Mechanics, 33, 274-281 20. Lysmer J., Kuhlemeyer R.L., 1969, Finite dynamic model for infinite media, Journal of Engi- neering Mechanics Division, ASCE, 95, EM4, 859-877 21. MaityD.,BhattacharyyaS.K., 1999,Timedomainanalysis of infinite reservoirbyfinit element method using a novel far boundary condition, International Journal of Finite Elements in Analysis and Design, 32, 85-96 22. MirzabozorgH.,GhaemianM.,KhalooA.R., 2003,Effect of reservoir bottom absorption on the seismic response of arch dams,Proceedings of the 4th International Conference of Earthquake Engineering and Seismology, Tehran, Iran 23. Mirzabozorg H., Ghannad M.A., Ghaemian M., 2003, Foundation interaction effect on the seismic response of Amir-Kabir arch dams, Proceedings of the 4th International Conference of Earthquake Engineering and Seismology, Tehran, Iran 24. Mojtahedi S., FenvesG.L., ReimerR.B., 1992,ADAP-88:A computer program for nonlinear earthquake analysis of concrete arch dams, Structural Engineering, Mechanics and Materials, De- partment of Civil Engineering, Report No. UCB/SEMM-92/11, University of California, Berkeley, California 25. Needleman A., Ortiz M., 1991, Effect on boundaries and interfaces on shear-band localization, International Journal of Solids and Structures, 28, 7, 859-877 594 A. Ferdousi et al. 26. OlsonL.G.,BatheK.J., 1983,A studyofdisplacement-basedfluidfinite elements for calculating frequencies of fluid and fluid-structure systems,Nuclear Engineering and Design, 76, 137-151 27. Ortiz M., Pandolfi A., 1999, finite-deformation irreversible cohesive elements for three- dimensional crack- propagation analysis, International Journal for Numerical Mechanics in En- gineering, 44, 1267-1282 28. Paul B., 1961,Modification of the Coulomb-Mohr theory of fracture, Journal of Applied Mecha- nics, 259-268 29. Raphael J.M., 1984, Tensile strength of concrete, Title No. 81-17, ACI Journal, March-April, 158-165 30. Roesset J.M., Ettouney M.M., 1977, Transmitting bouandaries: a comparison, International Journal for Numerical and Analytical Methods in Geomechanics, 1, 151-176 31. Ruiz G., Pandol A., Ortiz M., 2001, Three dimensional cohesive modeling of dynamic mixed- mode fracture, International Journal for Numerical Methods in Engineering, 52, 97-120 32. Sharan S.K., 1987, Time domain analysis of infinite fluid vibration, International Journal for Numerical Method in Engineering, 24, 945-958 33. Shugar T.A., Katona M.G., 1975, Development of finite element head injury model, Journal of the Engineering Mechanics Division, ASCE, 101, EM3, 11367, 223-239 34. Smith I.M.,GriffithsD.V., 2004,Programming the Finite ElementMethod, JohnWiley&Sons Press, USA 35. Sommerfeld A., 1949,Partial Differential Equations in Physics, Academic Press, NewYork 36. Tan H.C., Chopra A.K., 1996, EACD-3D-96: A computer program for three-dimensional ear- thquake analysis of concrete dams, Structural Engineering, Mechanics and Materials, Department of Civil and Environmental Engineering, Report No. UCB/SEMM-96/06, University of California, Berkeley California 37. Valiappan S., Zhao C., 1992, Dynamic response of concrete gravity dams including dam-water- foundation interaction, International Journal Numerical Methods for Engineering, 16, 79-99 38. WilsonE.L.,KhalvatiM., 1983,Finite elements for the dynamic analysis of fluid-solid systems, International Journal for Numerical Methods in Engineering, 19, 1657-1668 39. Wolf J.P., 1985,Dynamic Soil-Structure Interaction, Prentice Hall, EnglewoodCliffs, NJ 40. Wolf J.P, Darbre G.R., 1986, Nonlinear soil-structure interaction analysis based on the bo- undary element method in time domain with application to embedded foundation, Earthquake Engineering and Structural Dynamics, 13, 83-100 41. Wolf J.P., Obernhuber P., 1985, Nonlinear soil-structure interaction analysis using Green’s function of soil in the time domain,Earthquake Engineering and Structural Dynamics, 13, 213-223 42. Wolf J.P., Song C., 1996,Finite Element Modeling of Unbounded Media, Wiley, NewYork 43. Zhao C., Liu T., 2003, Non-reflecting artificial boundaries for transient scalar wave propagation in a two dimensional infinite homogenous layer, International Journal for Numerical Methods in Engineering, 58, 1435-1456 44. Zienkiewicz O.C., Bettes P., 1978, Fluid-structure dynamic interaction and wave forces an introduction to numerical treatment, International Journal for NumericalMethods in Engineering, 13, 1-16 45. Zienkiewicz O.C., Kelly D.W., Bettes P., 1977, The coupling of the finite elementmethods and boundary solution procedures, International Journal for Numerical Methods in Engineering, 11, 355-377 Manuscript received January 29, 2013; accepted for print December 30, 2013