Earthquake response of concrete arch dams: a plastic-damage approach

EARTHQUAKE ENGINEERING & STRUCTURAL DYNAMICS Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 Published online 22 July 2013 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/eqe.2317

Earthquake response of concrete arch dams: a plastic–damage approach Omid Omidi*,† and Vahid Lotfi Department of Civil and Environmental Engineering, Amirkabir University of Technology, Tehran, Iran

SUMMARY There are several alternatives to evaluate seismic damage-cracking behavior of concrete arch dams, among which damage theory is the most popular. A more recent option introduced for this purpose is plastic– damage (PD) approach. In this study, a special finite element program coded in 3-D space is developed on the basis of a well-established PD model successfully applied to gravity dams in 2-D plane stress state. The model originally proposed by Lee and Fenves in 1998 relies on isotropic damaged elasticity in combination with isotropic tensile and compressive plasticity to capture inelastic behaviors of concrete in cyclic or dynamic loadings. The present implementation is based on the rate-dependent version of the model, including large crack opening/closing possibilities. Moreover, with utilizing the Hilber–Hughes–Taylor time integration scheme, an incremental–iterative solution strategy is detailed for the coupled dam–reservoir equations while the damage–dependent damping stress is included. The program is initially validated, and then, it is employed for the main analyses of the Koyna gravity dam in a 3-D modeling as well as a typical concrete arch dam. The former is a major verification for the further examination on the arch dam. The application of the PD model to an arch dam is more challenging because the governing stress condition is multiaxial, causing shear damage to become more important than uniaxial states dominated in gravity dams. In fact, the softening and strength loss in compression for the damaged regions under multiaxial cyclic loadings affect its seismic safety. Copyright © 2013 John Wiley & Sons, Ltd. Received 1 December 2011; Revised 27 April 2013; Accepted 30 April 2013 KEY WORDS:

plastic–damage; arch dam; Koyna dam; dam–reservoir interaction; seismic analysis

1. INTRODUCTION Seismic safety of concrete dams as one of the main infrastructures needed for flood control or water supply is still an engineering concern, because it is increasingly evident that the seismic design concepts used at the time of construction of most existing dams were simplistic and inadequate. Furthermore, the population at risk located downstream of dams continues to expand. Therefore, the collapse of dams would cause a socioeconomic tragedy. The growing concern of the seismic safety of such critical structures has caused great interest for reevaluating existing dams by using nonlinear models that are able to predict crack initiation and its propagation through dam body. Damage and failure analyses of arch dams are very complicated because they need to be considered as 3-D systems with semi-unbounded foundation rock and reservoir domains. The dynamic analysis of arch dams should consider the following factors [1]: dam–water interaction, wave absorption at the reservoir boundary, water compressibility, dam–foundation interaction, spatially varying ground motion around the canyon, vertical contraction and peripheral joints opening and slippage, and, finally, a possible cracking in mass concrete monoliths. *Correspondence to: Omid Omidi, Department of Civil and Environmental Engineering, Amirkabir University of Technology, Tehran, Iran. † E-mail: [email protected] Copyright © 2013 John Wiley & Sons, Ltd.

2130

O. OMIDI AND V. LOTFI

Constitutive theory of concrete materials has been one of the main themes of research for some decades. Plasticity, continuum damage mechanics, and combined plasticity and damage mechanics are the most common theories that have been used for the description of concrete constitutive behavior. By the term damage-based, a class of smeared crack models such as fixed and rotating crack models is also included. Besides, a robust and efficient algorithm as a crucial issue of numerical implementation needs to be employed when a crack growth analysis in concrete dams is carried out. Plasticity theory has been widely utilized in concrete modeling. Damage mechanics as a relatively simple approach simulates concrete fracture by a degradation of stiffness due to the development of micro-cracks. Many isotropic and anisotropic damage models have been developed for concrete. However, because of the composite nature of concrete, its complicated behavior especially in cyclic loadings cannot be satisfactorily reflected in the usual constitutive theories of materials such as pure plasticity and pure damage mechanics theories [2]. In fact, the plasticity theory fails to address the degradation of the material stiffness due to micro-cracking, and the continuum damage mechanics also fails to describe the irreversible deformations and the inelastic volumetric expansion in compression. Therefore, because both micro-cracking and irreversible deformations are contributing to the nonlinear response of concrete, a constitutive model should address equally the two physically distinct modes of irreversible changes. Several effective plastic– damage (PD) models have been developed for concrete [2–13]. Replacing the isotropic elastic properties with orthotropic constitutive relations, Vargas-Loli and Fenves used fixed smeared crack models to represent the fracture zone in concrete dams [14]. They employed the model in seismic analysis of the Pine Flat dam in which the resulting crack profiles were diffused while numerical instability was also reported because of sudden energy release in small fractured elements. El-Aidi and Hall developed a fixed smeared crack model and seismically applied it to concrete gravity dams [15]. However, arbitrariness of the shear retention factor and a buildup of shear stresses along the fixed crack surface have been reported as disadvantages of the fixed smeared crack approach. Bhatacharjee and Leger employed a rotating smeared crack model as an improved smeared crack approach in the seismic analysis of gravity dams [16]. The model addresses softening by the secant stiffness in the direction of maximum principal strain. They defined the viscous damping on the basis of the tangent stiffness to include the energy dissipation. Espandar and Lotfi compared the application results of a non-orthogonal smeared crack model with an elastic perfectly plastic model to concrete arch dams, and it was concluded that the nonorthogonal smeared crack approach can redistribute the state of stresses and produces a more realistic profile of stresses in the dam [17–19]. Mirzabozorg and Ghaemian proposed a smeared crack to model static and dynamic behaviors of mass concrete in three-dimensional space [20]. The developed model examined in the seismic fracture analysis of the Koyna dam simulates the tensile fracture on mass concrete and contains pre-softening behavior, softening initiation, fracture energy conservation, and strain-rate effects under dynamic loads. Several researchers have used continuum damage mechanics to simulate mass concrete behavior in dams. Formulating the damage evolution, loading/unloading conditions, and crack opening/closing rules on the basis of the total strain, Ghrib and Tinawi proposed an orthotropic damage model [21, 22]. Cervera et al. developed an isotropic damage model that accounts for the different behaviors of concrete in tension and compression, each with its own damage surface and evolution law [23, 24]. Valliappan et al. evaluated the seismic response of gravity and arch dams by damage mechanics approach [25, 26]. Gunn proposed a 3-D damage model for static analysis of arch dams [27, 28]. Ardakanian et al. successfully employed the damage-based model proposed by Gunn for dynamic evaluation of arch dams [29]. Although much effort has been devoted to the development of coupled damage-plasticity constitutive laws for concrete in the recent years, there are a few studies in which a PD model is employed to assess the earthquake-induced damage of concrete dams [30, 31]. Lee and Fenves proposed a PD model for seismic analysis of concrete dams [32]. They applied the PD model, which had been implemented in 2-D plane stress state, to the Koyna gravity dam as a widely used benchmark problem. The PD model proposed by Lee and Fenves is developed herein in threedimensional formulation and implemented in a special finite element program, SNACS [33], for the seismic application to a typical arch dam. In this respect, the present study can be considered to be the first for this class of PD models to study the nonlinear seismic behavior of concrete arch dams. Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

PLASTIC–DAMAGE ANALYSIS OF ARCH DAMS

2131

Dam–reservoir interaction affecting the response of concrete dams has been widely investigated in the context of linear analyses. However, there are a few investigations carried out on the nonlinear analysis of concrete dams, including dam–water interaction. Utilizing a staggered solution scheme for the coupled dam–reservoir equations, Ghaemian and Ghobarah evaluated the seismic cracking of concrete gravity dams by a smeared crack model. The smeared crack model was found to be influenced by the dam–water interaction solution in comparison with the added mass approach commonly used [34, 35]. Considering the dam–water interaction, Calayir and Karaton examined the application of the smeared crack model proposed by Bhatacharjee and Leger to the Koyna dam [36]. The rest of the paper is outlined as follows. Section 2 describes the main features of the ratedependent PD model in the usual cracking state along with introducing the yield surface and the flow rule employed. The theoretical issues involved in the large cracking concept to properly address the cyclic behavior of concrete subjected to large tensile strains are also highlighted in this section. By using single-element tests, some basic validations for the 3-D implementation are discussed in Section 3. The incremental–iterative solution strategy implemented for the time integration of the coupled dam–reservoir equations in a nonlinear context is derived in Section 4. The Koyna gravity dam is then analyzed in a 3-D modeling, and the resulting crack pattern is compared with the available experiment. Finally, by applying the extended model to the seismic analysis of an arch dam, it is demonstrated that the PD model is eminently suitable for damage and failure analyses of arch dams.

2. MAIN FEATURES OF PLASTIC–DAMAGE MODEL The PD model implemented herein for three-dimensional analysis of concrete arch dams was established by Lubliner and coauthors in 1989, which was mainly appropriate for analyses with monotonic loadings [3]. The model was later modified by Lee and Fenves in 1998 to simulate cracking and crushing of concrete under cyclic or dynamic loadings [2]. Their proposed model is based on isotropic damaged elasticity in combination with isotropic multi-hardening plasticity formulated in the effective (undamaged) stress space. Although the damage part of the model is isotropic, two distinct damage variables are defined, which correspond to the two main failure mechanisms as tensile cracking and compressive crushing. Furthermore, two hardening variables connected to these failure mechanisms control the evolution of the yield surface [2]. By following Lubliner et al. [3], these hardening variables are referred to as the PD variables herein. 2.1. Constitutive law and damage evolution equation An overview of the rate-independent model as the backbone of the rate-dependent extension is initially presented. The basic constitutive equations of the PD model in tensor notation are s ¼ ð1  D Þ s ; « ¼ «e þ «p ; «_ p ¼ l_ rs Φ ; k_ ¼ l_ Hðs; kÞ ;

s ¼ E0 : «e D ¼ Dðs; kÞ

(1)

where «, «e and «p are the total, elastic, and plastic strains, respectively; s and s are the stress and the effective stress tensors, respectively; E0 is the initial (undamaged) elastic stiffness tensor; l is the plastic consistency parameter; Φ is the flow potential function; D is the stiffness degradation variable; and the vector H is the plastic modulus derived considering plastic dissipation. Moreover, k is the PD vector containing the normalized PD variables in tension and compression (i.e., k = {kt,kc}T), playing the role of hardening variables herein. The evolution equations of kt and kc are initially constructed by considering uniaxial loading conditions, and subsequently, they are extended to multiaxial states [2]. By using exponential functions, the uniaxial stress is assumed to be related to the scalar plastic strain symbolized by ep. In this condition, the PD variables are represented as Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2132

O. OMIDI AND V. LOTFI

kℵ ¼

1 gℵ

Z

Z

ep

sℵ dep ; 0

gℵ ¼

1

sℵ dep

ℵ 2 ft; cg

(2)

0

where gℵ called the specific fracture energy is the fracture energy normalized by the characteristic length as the localization zone size, lℵ (i.e., gℵ = Gℵ/lℵ). In a finite element analysis, this is employed to ensure mesh objectivity. Different choices are available for lℵ. Here, it is included as the side of an equivalent cube having the same volume as the tributary volume at each Gauss point of a 3-D solid element. It is clear that 0 ≤ kℵ ≤ 1 and k_ ℵ ≥0. The PD model utilizes the yield function of Lubliner et al. with the modifications proposed by Lee and Fenves to give different evolutions of strength in tension and compression [2, 3]. It takes the following form in terms of effective stresses: F ðs; kÞ ¼ f ðs; kÞ  cc ðkÞ    i 1 hpffiffiffiffiffiffiffiffi ^max ^max  g s f ðs; kÞ ¼ 3 J 2 þ a I 1 þ bðkÞ s 1a

(3)

^max is the maximum principal effective stress. The Macauley bracket h  i is defined by in which s hxi = (x + |x|)/2. The parameters a, b and g have the following definitions: a¼

ðf b0 =f c0 Þ  1 ; 2 ðf b0 =f c0 Þ  1

bðkÞ ¼

cc ðkc Þ ð1  aÞ  ð1 þ aÞ; ct ðkt Þ

g¼

3ð1  roct Þ 2 r oct  1

(4)

where fb0/fc0 is the ratio of the initial yield strengths under biaxial and uniaxial compression; ct and cc denote the effective tensile and compressive cohesion (positive values utilized here), respectively; roct pffiffiffiffiffiffiffiffiffiffiffiffi is the ratio of the octahedral shear effective stress (i.e., toct ¼ 2J 2 =3) on the tensile meridian to that on the compressive meridian (i.e., roct ¼ ðtoct ÞTM =ðtoct ÞCM ) at the initial yield condition for any given I 1 such that 0.5 ≤ roct ≤ 1.0. The value of roct = 2/3 suggested by Lubliner et al. leads to g = 3 [3]. Containing the parameter g, which appears only in triaxial compression, the yield function better predicts the concrete behavior in compression under confinement. F and l_ obey the loading/ _ ðs; kÞ ¼ 0. Furthermore, a Drucker–Prager hyperbolic _ unloading conditions: l≥0; F ðs; kÞ≤0; lF function is employed here as the plastic potential function: Φ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b2H þ 2 J 2 þ ap I 1 ;

bH ¼ e0 ap f t0

(5)

where ap is the dilatancy parameter, ft0 is the maximum uniaxial tensile strength, and e0 is the eccentricity parameter. It is noted that the original 3-D formulation of Lee and Fenves utilizes the linear function (i.e., Equation (5) when e0 = 0), which causes severe numerical difficulties in returnmapping process because of its apex’s singularity. Such a modification, which is inevitable for the 3-D implementation, requires a part of the model to be reformulated and makes an extra iteration necessary within the local iteration to compute the plastic consistency parameter as the main step of the stress return algorithm [37]. SNACS program utilizes the spectral return-mapping algorithm efficiently derived for this class of PD models [33, 38]. Although the damaged elastic stiffness remains isotropic, the model is accurately capable of capturing the tensile damage and compressive damage as the two major damage phenomena. Based on Lee-Fenves PD model [2, 32], a multidimensional degradation is defined by interpolating between the two main damage variables:       ^ D t ðk t Þ ^ k ¼ 1  ð1  D c ðk c ÞÞ 1  s s D s;

(6)

where the damage variables in tension and compression denoted by Dt and Dc, respectively, are explicit functions of the PD variables in tension and compression [2]. They can have values ranging from zero, corresponding to the undamaged state, up to one, which represents the fully damaged state (i.e., 0 ≤ Dℵ ≤ 1). Moreover, s is the recovery parameter such that 0 ≤ s ≤ 1 and used to include Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

PLASTIC–DAMAGE ANALYSIS OF ARCH DAMS

2133

the elastic stiffness recovery during elastic unloading process from tension to compression. Setting a minimum value, s0, for s, it may be written as [2] ! ! 3   3   X X         s ^ ¼ s0 þ ð1  s0 Þ r s ^ ; ^ ¼ ^ ≤1 ^i = ^i  s s r s s (7) ; 0≤r s i¼1

i¼1

^i ’s (i = 1, 2, 3) are the principal effective stress components and r is a multiaxial stress weight where s factor ranging from zero when all principal stresses are negative to 1 when they are positive. The stiffness is recovered upon crack closure as the load changes from tension to compression. On the other hand, it is not recovered as the load changes from compression to tension once crushing micro-cracks have developed. 2.2. Large cracking modification In a cyclic response, tensile behavior is close to pure damage, whereas response in compression is closer to pure plasticity, that is, an approximately secant path in tension and an elastic path in compression are followed in an unloading state. In fact, the crack opening/closure mechanism becomes similar to a discrete crack after a large damage is sustained in tension region, that is, kt > kcr where kcr is an empirical value usually near to unity [32]. At such a tensile damage level, the evolution of plastic strain caused by tensile damage is stopped. The modified evolution relation proposed by Lee and Fenves is utilized herein to treat large crack opening/closing and reopening e , the process in the continuum context [32]. By employing an intermediate effective stress, s evolution equation of plastic strain is modified as e ¼ E 0 : ð«  e «p Þ s «_ ; «_ ¼ ð1  rÞ e p

p

e j F ðs e ; kÞ≤0g fs  ^ e e «_ ¼ l_ rs r¼r s ~ Φ; 2

p

(8)

where e «p is referred to as the intermediate plastic strain, and r is the weight function causing the large crack modification to be applied for the states, which have positive principal stresses. To make the effective stress on the basis of the plastic strain admissible in the stress space, one has to employ a new damage variable denoted by Dcr at the crack damage corrector, making the evaluated effective stress return back onto the yield surface. It is determined by the plastic consistency condition for a continued loading ( D_ cr > 0 ) such that F ðð1  Dcr Þ s; kÞ ¼ 0 , from which Dcr ¼ 1  cc ðkÞ=f ðs; kÞ [32]. The stiffness degradation variable is redefined considering large cracking as D ¼ 1  ð1  Dc Þð1  s Dt Þð1  s Dcr Þ

(9)

With utilizing the spectral return-mapping algorithm proposed in [38], the numerical implementation including the algorithmic tangent operator and the stress update algorithm modified to consider large cracking possibilities has been detailed by Omidi and Lotfi in [39]. 2.3. Visco-plastic regularization Some of the convergence difficulties in material models addressing softening and degradation in behavior are treated by using a visco-plastic regularization of constitutive equations. Besides, it will improve mesh objectivity in structural simulations [32]. Both plastic strain and degradation variable are regularized herein by adding viscosity with the Duvaut–Lions regularization. The visco-plastic strain rate tensor, «_ vp , and the viscous degradation variable, Dv, for the visco-plastic system are defined, respectively, by 1 «_ vp ¼ ð«p  «vp Þ ; m

1 v D_ ¼ ðD  Dv Þ m

(10)

where m, which is called the viscosity parameter, shows the relaxation time of the visco-plastic model; «p and D are the plastic strain and degradation variable computed in the inviscid backbone model. The stress– strain relation for the visco-plastic model is given as [32] Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2134

O. OMIDI AND V. LOTFI

s ¼ ð1  Dv Þ E0 : ð«  «vp Þ

(11)

2.4. Visco-elastic-damage damping stress Several approaches are available to include material damping effects in a transient analysis, among which the most common one is based on the constant damping mechanism. This kind of damping utilizing the linear form of visco-elastic damping stress introduces artificial damping forces as cracks open with large relative velocity, and the damping forces then restrains the crack opening by transferring artificial stresses across the crack. This artificial resistance results in unreal diffused cracks in the cases where cracking should be highly localized [40, 41]. Three alternatives have been employed to solve this problem: (i) setting the damping term to zero upon cracking of an element [15]; (ii) using the tangent stiffness matrix rather than the elastic one [16]; and (iii) utilizing the degraded stiffness in computation of the visco-elastic damping stress [32]. By following Lee and Fenves, a nonlinear form of visco-elastic damping, which includes the degradation damage variable, D, computed within the rate-independent backbone model, is employed here as xð«; «_ Þ ¼ bR ð1  DÞE0 : «_

(12)

where bR is the coefficient of viscous damping and calibrated to provide a damping ratio at the fundamental period. This kind of damping stresses introducing damage-dependent damping mechanism represents stiffness proportional term for damping of damaged material.

3. BASIC VALIDATIONS OF PLASTIC–DAMAGE MODEL By using SNACS program [33], several single-element tests initially validating the 3-D extended model are discussed herein. In all cases, a = 0.12 and ap = 0.2 are used as the model parameters. 3.1. Cyclic loading example under large tensile strain To examine how the large cracking phenomenon is being captured in the model, one subjects a single element to full cyclic uniaxial loadings under a large tensile strain. The following material properties are utilized: E0 = 31.7 GPa, n = 0.18, f 0 t = 3.48 MPa, f 0 c = 27.6 MPa, Gt = 12.3 N/m, Gc = 1750 N/m and lt = lc = 25.5 mm. This test exemplifies the ability of the model to simulate stiffness recovery when the status changes from tension to compression and vice versa. The result with the loading paths A to G is illustrated in Figure 1(a). During the tensile unloading and subsequent compressive loading (path A–B or D–E–F),

Figure 1. Crack opening/closing simulation: (a) stress–strain curve and (b) visco-plastic effects. Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2135

PLASTIC–DAMAGE ANALYSIS OF ARCH DAMS

the stiffness is properly recovered. Large crack opening/closing is simulated during path C–D–E. The crack closes at point E, and the elastic stiffness is recovered. It yields in compression again to reach point F. Subsequent reopening of the crack is noticed along path F–G during unloading. In Figure 1(b), the full cyclic loading response of the rate-dependent model is compared with the rate-independent result. The strain rate of 0.0002 sec-1 is imposed through the full cycles. The differences in the responses emerge only at the softening and hardening regions. The stiffness recovery is well represented during the unloading from tension to compression in the analysis including rate dependency. 3.2. Concrete under 3-D compression The g term in the yield function is adjusted to have better prediction in 3-D applications loaded in compression under confinement. The results of the model are compared for this type of loading between the cases of g = 0 and g = 3. The material properties are the same as those utilized for the cyclic loading test. The numerical results under three sets of confining stresses are presented in Figure 2. The confining stresses are increased proportional to the imposed strain in the third direction for all these cases. It is clearly observed that the enhancements of strength and ductility due to confinement are satisfactorily captured.

Figure 2. s1  e1 curves under three different confinements when (a) g = 0 and (b) g = 3.

3.3. Shear panel test With imposing a special in-plane boundary condition as illustrated in Figure 3(a), a 3-D single-element mesh with unit dimensions is used to study multiaxial cyclic behavior of a concrete panel. This example was adopted from [38]. Such a boundary condition imposed to the element does not provide a pure shear state. The out-of-the-plane boundary condition is considered as plane stress. The material properties are E0 = 30.0 GPa, n = 0.2, f 0 t = 3.3 MPa, f 0 c = 30.0 MPa, Gt = 250 N/m, Gc = 25000 N/m and lt = lc = 0.5 m. The results are compared in Figure 3(a) for three different sets of displacements as 20, 40, and 80 steps per one cycle of loading to assess the accuracy of the solution. Excellent agreement exists between the three analyses through all stages. This indicates good accuracy of the solution over a large range of strain increments. Figure 3(b) illustrates the effects of rate dependency on the results by considering different viscosity ratios. The imposed shear strain rate is 0.0004 sec-1. To examine how well the stiffness degrades during the cycles, one can subject the panel to several cyclic loadings. The hysteresis curves computed by the rate-independent model are given in Figure 4(a). The softening and strength loss of the panel are precisely observed under multiaxial cyclic loading. Furthermore, Figure 4(b) in which the time represents loading sequence illustrates the evolutions of tensile and compressive damages. This example clearly shows that multiaxial stress states may cause tensile and compressive damages to occur simultaneously, which is usually interpreted as shear damage although compressive stresses do not reach the compressive strength limit [38, 42]. Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2136

O. OMIDI AND V. LOTFI

Figure 3. Shear stress versus shear strain for the shear panel under cyclic loading: (a) with different size steps (20, 40, and 80 steps) and (b) effects of rate dependency.

Figure 4. Shear panel test under several loading loops: (a) shear stress vs. shear strain (b) evolutions of tensile and compressive damage variables.

4. TIME INTEGRATION SCHEME WITH DAM–RESERVOIR INTERACTION The Lagrangian–Eulerian formulation, in which displacement and pressure are unknown variables for dam and reservoir, respectively, is utilized to construct the governing coupled equations for dam–reservoir system [43]. Herein, the fundamental equations of equilibrium are briefly outlined first for the dam and the reservoir separately. Afterward, the handling of the coupled governing equations is presented by the Hilber–Hughes–Taylor (HHT) time integration scheme in an incremental–iterative solution strategy. 4.1. Coupled dam–reservoir equations in a nonlinear context The equations of motion may be written in vector-matrix notation as follows for the discretized dam under specified free-field ground acceleration: € þ Pðu; u_ Þ ¼ Rst þ QT p MU

;

€ ¼u € þ J ag U

(13)

€ represent the nodal vectors of the displacement, where M is the mass matrix of dam body; u, u_ and u € is the vector of absolute nodal velocity, and acceleration relative to the ground motion, respectively; U Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2137

PLASTIC–DAMAGE ANALYSIS OF ARCH DAMS

acceleration; P denotes the restoring force vector; Q is often referred to as interaction or coupling matrix; p represents the vector of nodal hydrodynamic pressures and Rst is the static loads vector due to weight of the dam and hydrostatic pressure applied to upstream face of the dam; ag denotes the vector of ground accelerations; and J is a matrix with each of the three rows equal to a 3  3 identity matrix. It is noted that in the current implementation, the element’s restoring force vector, Pe, is given by the integral of the total stress including damping stress over the element volume as Z

P ¼ e

Ve

BT ðs þ xÞ dVe

(14)

where B is the strain–displacement transformation matrix. In this condition, the stiffness and damping matrices of an element (i.e., Ke = @ Pe/@ ue and Ce ¼ @Pe =@ u_ e ) are derived as Z

Ke ¼

Z Ve

BT Eats B dVe ;

Ce ¼ b R

Ve

BT ð1  DÞE0 B dVe

(15)

where Eats = @ (s + x)/@ « denotes the modified algorithmic tangent stiffness with damping. The viscoelastic damping stress is included in Eats to keep the convergence rate rather unaffected by rapid changes of damping forces [32]. The tangent operator, which is consistent with the stress updating algorithm, is derived from the linearized constitutive relations around the given state. Considering large cracking possibilities and visco-plastic regularization, the derivation details of the algorithmic tangent stiffness have been given in [39]. Furthermore, according to El-Aidi and Hall [15], the mass-proportional term for the damping matrix has been omitted, because it would provide some artificial numerical stability during time-marching process. Furthermore, one can apply Galerkin method to the wave equation and impose related boundary conditions to obtain the following relation for the discretized reservoir domain [43]: € € þ L p_ þ H p ¼ r Q U Gp

(16)

where G, L and H are called generalized mass, damping, and stiffness matrices, respectively. To construct the coupled equations set of the system consistently, the governing equations for the dam and reservoir are restated as the following, respectively: € þ Fðu; u; _ pÞ ¼ Rst ; MU

_ pÞ ¼ Pðu; u_ Þ  QT p Fðu; u;

0 € ; € þ F ðp; p_ Þ ¼ r Q U Gp

(17)

0

F ðp; p_ Þ ¼ L p_ þ H p

(18)

Now, one may set up the following coupled equations:

M

0

  € u

rQ G

€ p

þ

F F

0



¼

R R

 (19)

0

€ as total displacement, velocity where R = Rst  M J ag and R0 =  r Q J ag. With considering u, u_ and u and acceleration vectors with the following definitions: u¼

 u p

; u_ ¼



 € u_ u €¼ ; u € p_ p

(20)

and denoting M as the total mass matrix, Equation (19) is written in a more compact form: Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2138

O. OMIDI AND V. LOTFI



€ þ F ¼ R; Mu

M¼

M

0

rQ G



; F¼

F F

0



; R¼

R R

 (21)

0

4.2. Iterative solution strategy with the Hilber–Hughes–Taylor time integration scheme Numerical damping associated with the time-stepping schemes used for integrating second-order systems of equations over time stabilizes the numerical integration by damping out the unwanted high-frequency modes. For the Newmark method, numerical damping also affects the lower modes and reduces the accuracy of integration scheme from second order to first order. For the HHT method, numerical damping affects only the higher modes and always maintains second-order accuracy. Larger numerical damping values are usually necessary for problems involving rigid body rotational motion and dynamic contact/impact [32]. Utilizing the HHT integration scheme, the discrete form of the governing equilibrium equation, Equation (21), is written at the time step n + 1 as   €nþ1 þ F unþ’ ; u_ nþ’  Rnþ’ ¼ 0 Mu

(22)

where ’ = 1 + a and a are the free parameter controlling the amount of the numerical dissipation, which is increased by decreasing a. The recommended range for a is [1/3, 0] to achieve unconditional stability, second-order accuracy, and good high-frequency dissipative characteristics. Furthermore, on the basis of the generalized midpoint scheme, the interpolated quantity is defined as ðÞnþ’ ¼ ð1  ’ÞðÞn þ ’ ðÞnþ1

(23)

It is noted that this interpolation scheme is only used for the equations of motion (i.e., unþ’ and u_ nþ’) and the applied force, Rnþ’ . Because employing the generalized midpoint scheme in the stress computation does not have robust numerical performance [44], the stresses needed for the restoring force, Pnþ’, are integrated on the basis of the shifted backward-Euler scheme [38], in which the internal variables are directly computed at tnþ’ rather than tnþ1. The shifted scheme becomeso the n

original method in [tnþ’  1,tnþ’], for which the given variable set «nþ’ ; «pnþ’1 ; knþ’1 is

utilized to find the updated variable set snþ’ ; «pnþ’ ; knþ’ : Because Equation (22) is nonlinear with respect to u, a global iteration scheme is employed to update it by defining the residual of Ψ at the time step n + 1 as €nþ1 þ Fnþ’  ð1  ’Þ Rn  ’ Rnþ1 Ψnþ1 ¼ Μ u

(24)

By applying the Taylor’s expansion to this residual and eliminating the higher-order terms, the  ðiÞ ðiÞ equation @Ψ=@u nþ1 du ¼ Ψnþ1 is established for the ith iteration in which @Ψ=@u is the Jacobian of the residual Ψ with respect to u. The vector of du is subsequently utilized to update ðiþ1Þ ði Þ unþ1 (i.e., unþ1 ¼ unþ1 þ du). From Equation (24), the Jacobian is derived as

@Ψ @u



€nþ1 @u ¼M þ @unþ1 nþ1

Copyright © 2013 John Wiley & Sons, Ltd.

@Fnþ’ @ u_ nþ’ @ u_ nþ1 @Fnþ’ @unþ’ þ @ u_ nþ’ @ u_ nþ1 @unþ1 @unþ’ @unþ1

! (25)

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2139

PLASTIC–DAMAGE ANALYSIS OF ARCH DAMS

€ and u_ are defined as Furthermore, according to the Newmark formulation, u €nþ1 ¼ a0 ðunþ1  un Þ  a2 u_ n  a3 u €n ; u

€n u_ nþ1 ¼ a1 ðunþ1  un Þ  a4 u_ n  a5 u

(26)

in which, a0 through a5 are the coefficients of the Newmark approach. By using the above relations along with Equation (23), the Jacobian matrix yields 

   @Ψ=@u nþ1 ¼ a0 M þ ’ a1 Cnþ’ þ Knþ’

(27)

where C and K are called total damping and stiffness matrices, respectively as !

 C 0 K QT C¼ ; K¼ 0 L 0 H

(28)

€ and u_ are updated at the end of the iteration by After obtaining du, u €ðiÞ þ a0 du ; €ðiþ1Þ ¼ u u nþ1 nþ1

ðiþ1Þ ðiÞ u_ nþ1 ¼ u_ nþ1 þ a1 du

(29)

Moreover, the initial values for these vectors at the beginning of a time step would be as €ð0Þ ¼ a2 u_ n  a3 u €n ; u nþ1

ð0Þ €n ; u_ nþ1 ¼ a4 u_ n  a5 u

ð 0Þ

unþ1 ¼ un

(30)

Finally, the HHT time integration algorithm is summarized in Table I. It is worthwhile to mention that under certain circumstances, it is possible to utilize the pseudo-symmetric technique relying on symmetric matrices for efficiency purposes [33, 43, 45]. Table I. The HHT time integration algorithm with dam–reservoir interaction. €ð0Þ ¼ a2 u_ n  a3 u €n ; u_ ð0Þ ¼ a4 u_ n  a5 u €n and uð0Þ ¼ un : 0. i ¼ 0; agnþ1 ; u nþ1 nþ1 nþ1

 Rnþ1 0 g g st 1. Rnþ1 ¼ R  M J anþ1 , R nþ1 ¼ r Q J anþ1 and Rnþ1 ¼ . 0 R nþ1 ðiÞ ð iÞ 2. unþ’ ¼ ð1  ’Þ un þ ’ unþ1 . iÞ 3. «nþ’ ¼ B ðue Þðnþ’ , then obtain snþ’ and Dnþ’ on the basis of the plastic–damage model. ðiÞ ð iÞ 4. u_ nþ’ ¼ ð1  ’Þ u_ n þ ’ u_ nþ1 .

  iÞ 5. «_ nþ’ ¼ B ðu_ e Þðnþ’ , then compute xnþ’ ¼ bR 1  Dnþ’ E0 «_ nþ’ .  XZ ðiÞ 6. Pnþ’ ¼ BT snþ’ þ xnþ’ dVe . 0 e 7.

ðiÞ Fnþ’

¼

V nel ðiÞ Pnþ’  QT

ð iÞ pnþ’ ,

0 ð iÞ

F nþ’ ¼ L

ð iÞ p_ nþ’

þH

ðiÞ pnþ’

ðiÞ €ðiÞ þ FðiÞ  ð1  ’Þ Rn  ’ Rnþ1 . 8. Ψnþ1 ¼ Μ u nþ1 nþ’    ðiÞ  9. If Ψnþ1 ≤TolG , then exit.    ð iÞ ðiÞ ðiÞ 10. @Ψ=@u nþ1 ¼ a0 M þ ’ a1 Cnþ’ þ Knþ’ .

and

ðiÞ Fnþ’

¼@

ðiÞ

Fnþ’ 0 ðiÞ

F nþ’

1 A.

 ðiÞ ðiÞ 11. Solve @Ψ=@u nþ1 du ¼ Ψnþ1 for du. €ðiþ1Þ ¼ u €ðiÞ þ a0 du; u_ ðiþ1Þ ¼ u_ ðiÞ þ a1 du and uðiþ1Þ ¼ uðiÞ þ du. 12. u nþ1 nþ1 nþ1 nþ1 nþ1 nþ1 13. i = i + 1, then go to Step 2.

Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2140

O. OMIDI AND V. LOTFI

4.3. Verification of the dam–water interaction modeling By utilizing the derived time integration algorithm implemented in SNACS program, a linear dynamic analysis of the Pine Flat dam is carried out, and the response is compared with the corresponding results reported by Chopra et al. [46]. It is clear that in a linear analysis, a nonlinear iterative procedure can be utilized. However, the system’s equations converge in the first iteration of each step. The 3-D analysis is performed in a state of plane stress over the tallest monolith of the dam. The idealization of 20-node isoparametric elements utilized to discretize solid and fluid domains is depicted in Figure 5(a). The truncation length of the reservoir was checked to be adequate for the considered horizontal excitation depicted in Figure 5(b). Rayleigh damping coefficients are determined such that damping would be equivalent to 5% of critical damping for frequencies close to the first and third vibration modes. The concrete properties are E0 = 22.75 GPa, n = 0.2 and unit weight of 24.8 kN/m3. Moreover, the bottom of the reservoir is assumed to be completely reflective, and pressure wave velocity of water is 1440 m/s. Prior to the seismic loading, the dam is subjected to gravity and hydrostatic loads. The crest stream displacement history is compared with the corresponding result of Chopra et al. [46] in Figure 6(a). It should be mentioned that the employed approach in the work of Chopra et al. is based on infinite elements for the reservoir discretization in a frequency domain analysis. The envelopes of maximum tensile stresses are also compared in Figure 6(b) between the two mentioned analyses. The comparison reveals that the accuracy of the present implementation is well enough for the solution of the coupled governing equations although it is being carried out in time domain.

Figure 5. Pine Flat dam: (a) mesh and (b) the S69E component of Taft earthquake.

Figure 6. SNACS’s results compared with the responses obtained by Chopra et al. [46]: (a) time history of crest stream displacement and (b) envelope of maximum tensile stresses. Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2141

PLASTIC–DAMAGE ANALYSIS OF ARCH DAMS

5. FRACTURE RESPONSE OF THE KOYNA DAM IN 3-D MODELING This section focuses on the seismic fracture response of the Koyna dam, which has been broadly examined by many investigators as the main validation of their proposed concrete models [16, 20, 32, 36, 47]. The Koyna dam is a 103-m high concrete gravity dam located in India and constructed in 1963. It was subjected to the 1967 earthquake of magnitude 6.5 and was severely damaged [48]. Its scaled model was also experimentally tested in 1988 [49]. 5.1. Model description and loading The simulation is based on the 3-D modeling of a slice of the tallest monolith depicted in Figure 7(a). The cracked configuration experimentally observed by a shaking table test is also shown in this figure [49]. The eight-node isoparametric solid and fluid elements having two integration points in each direction are utilized to discretize the dam and reservoir domains, respectively. As illustrated in Figure 7(b), two meshes are employed in the simulation to examine mesh objectivity. The material properties for the PD model are E0 = 30.0 GPa, n = 0.2, ft0 = 2.9 MPa, fc0 =  24.1 MPa, Gt = 200 N/m, Gc = 20000 N/m and m = 0.15 Δt. Moreover, Δt = 0.005 sec and a =  0.3 are selected as the analysis parameters. The thickness of the slice is 1 m, and the dam body is assumed to be in a state of plane stress. The length of the reservoir simulated in the model is 2.5 times of the reservoir height, whereas the reservoir bottom is assumed to be partially reflective with the wave reflection coefficient equal to 0.85. The unit weights of concrete and water are assumed as 25.8 and 9.81 kN/m3, respectively. The stiffness proportional damping factor is computed such that 3% damping is being captured for the fundamental vibration period. At the first step of loading, the dam is subjected to gravity load, and then, the hydrostatic pressure is applied to the upstream face. The static loads are applied at negative range of time in 10 increments each before the earthquake excitation is employed. Afterward, 10 sec of the horizontal and vertical components of the 1967 earthquake records [32] are applied to the base of the dam and reservoir starting at time 0. The time integration is performed for 15 sec of which the last 5 sec is free vibration. 5.2. Plastic–damage analysis The analysis commences with applying static loads while the dam remains elastic with no damage at the end of this step. Damage to the dam initiates during the earthquake loading. Figure 8(a) shows the stream and vertical displacements at the left corner of the crest relative to the ground motion, whereas positive horizontal displacement is in the downstream direction. Furthermore, the contour plots of tensile damage, compressive damage, large-crack damage, and degradation variables at the end of the analysis are illustrated in Figure 8(b). These are denoted as Dt, Dc, Dcr and D, respectively. The resulting crack profile is divided into three major parts. The rigid foundation causes the stress concentration to occur near the base, and consequently, a crack propagates into the dam along the lowest layer of the elements. As observed in Figure 8(b), a crucial damage has

Figure 7. Koyna gravity dam: (a) geometry (unit, meter) and (b) coarse and fine meshes. Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2142

O. OMIDI AND V. LOTFI

Figure 8. Seismic results of Koyna dam with the coarse mesh: (a) time histories of the crest stream and vertical displacements and (b) damage and degradation variables at the end of the analysis.

occurred near the downstream change of slope where the localized band of damaged elements is formed. Curving down due to the rocking motion of the upper block, this downstream crack spreads toward the upstream direction. Furthermore, a horizontal crack propagating toward the curved part of the downstream crack is formed during a major deformation. As noticed, the damage has very well localized on just one layer of the elements along the upstream face. The compressive damage variable, Dc, depicted in Figure 8(b) shows that a few number of the damaged elements experience shear damage. Contour of the large-crack damage variable, Dcr, shown in Figure 8(b) reveals that most of the damaged elements also undergo large cracking phenomenon. It is worthwhile to mention that all cracks have occurred in the body by the time 4.410 sec, which corresponds to the second major excursions of the crest in the upstream direction. Oscillating back and forth during the remainder of the analysis including the free vibration, the upper portion of the dam remains stable while the upstream and downstream cracks alternatively close and open because of rigid body rotation of the upper block. It is simply observed from the presented displacement histories. The distribution of the stiffness degradation variable, D, in Figure 8(b) indicates that most of the cracks close under compressive stresses at the end of the analysis. The well-localized damage pattern predicted by SNACS is consistent with those numerically captured by other investigators [20, 47] as well as the cracked zone obtained in the scaled experimental test [49]. To examine capability of the implemented model on mesh objectivity concern, the authors also present results of the fine mesh introduced in Figure 7(b) herein. The stream displacements at the dam crest for the fine mesh are compared with the corresponding results of the coarse mesh in Figure 9(a). This demonstrates reasonably objective results for the two meshes. Although 3-D models could have relatively higher level of mesh sensitivity because of the effect of possible cracking in the out-of-plane direction, the characteristic length along with the visco-plastic regularization provides a good mesh objective result in this model. Besides, the tensile damage pattern depicted in Figure 9(b) for the fine mesh at the end of the analysis is fairly similar to the damaged zone of the coarse mesh, although as expected, the crack bandwidth is smaller. The inclined part of the cracking has mostly two element widths, which is inevitable when the alignment of the elements is not coincident with the crack path. Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

PLASTIC–DAMAGE ANALYSIS OF ARCH DAMS

2143

Figure 9. Seismic response of Koyna dam with the fine mesh: (a) crest displacement history of the fine mesh compared with the coarse mesh’s result and (b) tensile damage, Dt, at the end of the analysis.

6. APPLICATION ON CONCRETE ARCH DAMS In this section, the seismic nonlinear behavior of Shahid Rajaee arch dam analyzed by SNACS is discussed. The dam is 130-m high, with the crest length of 420 m, and it is located in the north of Iran in the seismically active foothills of Alborz Mountains, near the city of Sari [17–19]. Two analysis cases are presented herein: a linear one mainly used for comparative purposes and a nonlinear case by the PD model. 6.1. Mesh, model parameters and loading The dam body and its reservoir are discretized by the isoparametric 20-node elements having three integration points in each direction. A relatively fine mesh employed in the analyses and two key Gauss points monitored are depicted in Figure 10. The reservoir length considered in the simulation is 2.75 times of the reservoir height. The finite element mesh consists of 480 solid elements and 714 fluid elements for the dam body and the reservoir domain, respectively. As it is clear in this case, with knowing that cross-canyon excitation is neglected, simulation might be performed on one half of the dam–reservoir system by applying symmetry condition in the midplane. However, this was disregarded, and the whole domain is discretized. Therefore, by neglecting the cross-canyon excitation, symmetric results are expected throughout the analyses, and this fact is used as an additional check for the accuracy of the obtained results. Dam–reservoir interaction is included in the analyses based on the technique discussed in Section 4. Meanwhile, the foundation is taken as rigid. This latter idealization was decided to reduce computational efforts. Of course, it is evident that the assumption of rigid foundation would cause large tensile stresses near the boundaries for the linear analysis. However, for the PD analysis, these high tensile stresses are expected to release, and this can be considered as a challenging test for the nonlinear model. Material properties for the PD model are E0 = 30.0 GPa, n = 0.18, ft0 = 3.0 MPa, fc0 =  25.0 MPa, Gt = 300 N/m, Gc = 45000 N/m, m = 0.15 Δt and unit weight of 24.0 kN/m3. Moreover, the analysis parameters are selected as Δt = 0.005 sec and a =  0.2. In the analyses carried out, the Rayleigh stiffness proportional damping is applied and the damping coefficient, bR, is determined such that equivalent damping for the fundamental frequency of vibration would be 8% of the critical

Figure 10. Finite element mesh of the Shahid Rajaee arch dam and the two Gauss points monitored. Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2144

O. OMIDI AND V. LOTFI

damping. The reservoir bottom is partially reflective with the wave reflection coefficient equal to 0.9. The dead weight and hydrostatic pressures at upstream face are applied as the static loading at negative time. Maximum water depth is 122 m. The earthquake excitations applied at time 0 include two components of the Friuli–Tolmezzo earthquake (as mentioned, the cross canyon is neglected to have a symmetric condition) whose records are normalized on the basis of the frequency content for maximum design earthquake condition with the peak ground acceleration of 0.56 g [18].

6.2. Analysis results For the linear case, it is observed that maximum tensile stresses for the spillway and abutment regions reach to maximum values of 8.06 and 24.84 MPa, respectively (Figure 11). These high tensile stresses are predicted to be limited to tensile strength of concrete in the PD analysis. The maximum tensile stresses of these zones occur in the arch direction and perpendicular to the abutments as expected. The maximum compressive stresses of this case are 14.43 and 11.65 MPa for the spillway and abutment regions, respectively (Figure 12). Although they are less than the compressive strength, it is shown in the nonlinear analysis that compressive damages would occur because of shear damages. In Figure 13(a) and (b), displacement histories of the mid-crest point for the PD case are compared with those of the linear case in the stream and vertical directions, respectively. The analyses did not expose large displacements in the dam under the severe earthquake excitations considered, giving an initial sign of a safe design of the dam on the basis of the adopted theories and assumptions. Comparison between the results of the PD model and those of the linear case shows a drift in the crest displacements in upward and upstream directions. This is the main feature of displacement histories in the PD approach, which is mostly due to extensive cracking of downstream face of the dam beneath the spillway block. When most of these cracks are open, they create a cumulative plastic strain as the major source of the drift. The envelope of maximum tensile stresses throughout the PD analysis is depicted in Figure 14. In the linear analysis, it was noticed that very high tensile stresses were computed at the base of the dam because of the rigid foundation assumption. In addition, high tensile stresses occurred in the spillway region in arch direction, which in reality are released with the opening of contraction joints. It is observed that the PD model has completely bounded the amount of tensile stresses to a value close to the tensile strength of concrete. In this regard, slight overstressing of tensile stresses (about 1.7%) is still noticed in the dam body for this case, because of three-dimensional stress states’ effects. Furthermore, the tensile damages obtained for different regions of the dam body at the end of the

Figure 11. Envelope of maximum tensile stresses (MPa) for the linear analysis.

Figure 12. Envelope of maximum compressive stresses (MPa) for the linear analysis. Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2145

PLASTIC–DAMAGE ANALYSIS OF ARCH DAMS

Figure 13. Comparison of displacements at mid-crest point between the linear (LN) and plastic–damage (PD) models in (a) stream direction and (b) vertical direction.

Figure 14. Envelope of maximum tensile stresses (MPa) for the plastic–damage analysis.

analysis are illustrated in Figure 15. The regions undergoing tensile damage are well coincident with the high tensile zones depicted in Figure 14 as the envelope of maximum tensile stresses. Time histories of the major stress components at the two key points being monitored are compared between the LN and PD cases in Figure 16. As presented in Figure 16(a), the linear response shows very high values for sz at the base of the central cantilever, whereas it is released during the static loading in the PD case. The arch stress component as illustrated in Figure 16(b) is limited to the tensile strength during the seismic loading while linear analysis estimates a value of 8.06 MPa. The envelope of maximum compressive stresses is captured in Figure 17. In comparison with the linear case, the maximum compression has increased about 26.1% and reached to 18.19 MPa in the upstream face near the spillway region. Because the high compressive stress values are still below the considered compressive strength, it is not expected to have compressive damages due to crushing. Although the maximum value of the stress components is less than that of the compressive strength, it cannot be concluded that no compressive damage is occurred in the dam body. This is mainly due to the governing equation for damage evolution in the implemented PD model. In fact, because all three principal components of stress contribute in the damage evolution, it is expected

Figure 15. Tensile damage, Dt, at the end of the plastic–damage analysis. Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2146

O. OMIDI AND V. LOTFI

Figure 16. Comparison of stresses in the index points at the upstream face of the dam:(a) sz at the base (cantilever stress) and (b) sx at the crest (arch stress). LN, linear; PD, plastic–damage.

that the compressive damage evolves when the tensile damage is occurring and at least one of the principal stresses is negative. As indicated in Section 3.3, the state of damage when both tensile and compressive damages occur simultaneously is commonly interpreted as shear damage. This phenomenon was demonstrated in the shear panel test previously discussed. The distribution of compressive damage variable at the end of the analysis is depicted in Figure 18. This figure together with the tensile damage pattern, which actually shows the regions that both damage variables are arising simultaneously, illustrates the extent of shear damage occurring in these zones of the dam body because of multiaxial stress states governed in such mass concrete structures. The compressive damage occurred at the zones in which tensile damage is also observed highlights that most part of the damaged area in arch dams sustains shear damages, which affects the safety margin with respect to the cases in which the tensile damage occurs alone. Besides, the large-crack damage variable, Dcr, computed for the points undergoing large cracking is presented in Figure 19. As observed, an extensive part of the damaged regions sustain large cracking. This emphasizes the necessity of such a modification in simulation of cyclic loaded concrete subjected to large tensile strains while this class of PD models is utilized.

Figure 17. Envelope of maximum compressive stresses (MPa) for the plastic–damage analysis.

Figure 18. Compressive damage, Dc, at the end of the plastic–damage analysis. Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2147

PLASTIC–DAMAGE ANALYSIS OF ARCH DAMS

Figure 19. Large-crack damage variable, Dcr, at the end of the plastic–damage analysis.

7. CONCLUSIONS A special finite element program called SNACS is developed in three-dimensional space on the basis of a well-established PD model, which was successfully applied to gravity dams in 2-D plane stress state, to investigate seismic nonlinear behavior of concrete arch dams. The model was initially proposed by Lubliner and coauthors in 1989 and later on modified by Lee and Fenves in 1998 to make it applicable in cyclic loadings such as earthquake. The model accounts for the different behaviors in tension and compression. Tensile and compressive damages are independently represented by the concept of the fracture-energy-based multi-hardening variables. Effects of damage on elastic stiffness and its recovery during crack opening/closing are well simulated by a simple scalar degradation model. Although the damaged elastic stiffness is isotropic, the model offers the different evolutions of the tensile and compressive strengths and induced directional damage by plastic strains. In this paper, an overview of the rate-dependent PD model was presented first, and then, the coupled dam–reservoir equations were integrated by the HHT scheme within an incremental–iterative solution procedure. By utilizing the damage-dependent damping mechanism, the varying restoring force vector is computed in which the damping stress is consistently included. The developed program was initially examined by several single-element tests. A full uniaxial cyclic test was carried out to emphasize the necessity of large crack modification to properly simulate the cyclic behavior of plain concrete subject to large tensile strains. Furthermore, to show the performance of the model in reproducing the softening and strength loss of concrete under multiaxial cyclic loading and the corresponding stiffness recovery during unloading and reloading, the shear panel example was studied. This test well demonstrates that the model is capable of capturing shear damage represented by the two distinct damages in tension and compression occurring simultaneously. Afterward, SNACS was verified by investigating the seismic fracture response of the Koyna gravity dam as a well-known benchmark problem and by comparing the resulting crack pattern with the corresponding experimental data. By employing two meshes for the Koyna dam, mesh objectivity of the current implementation of the model was well examined. Finally, by utilizing the 3-D extended PD model, the Shahid Rajaee concrete arch dam was analyzed. The dynamic instability was not observed despite extensive damages occurring in the dam body. The conducted analysis included the dynamic effects of the reservoir and also the wave absorption at the reservoir boundaries. It is confirmed that the PD model gives reasonable crack profiles, which is expected from the contour of the maximum tensile stresses obtained from linear dynamic analysis. The results also reveal that the concrete arch dams can suffer significant damage during a strong earthquake and still remain stable. The results are compared with the solutions from their corresponding linear analysis, and it is shown that the structural responses of dams vary considerably because of damage evolution being captured by the PD model. The prominent characteristic of the crest displacement histories in the PD analysis is a drift in upward and upstream directions, which is mostly due to a cumulative plastic strain caused by extensive cracking of downstream face of the dam beneath the spillway block when most of these cracks are open. Because multiaxial stress state leads to shear damages represented with two distinct damages occurring simultaneously, it is important to carefully evaluate the damaged regions observed in an earthquake damage analysis of an arch dam because the concrete may in certain cases partially lose its strength in compression while it has fully damaged in tension. Although it is less important for gravity dams, shear Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

2148

O. OMIDI AND V. LOTFI

damages causing the softening and strength loss in compression for the damaged regions under multiaxial loadings affect the seismic safety of an arch dam and need to be carefully captured. Future work includes seismic damage-cracking analysis of a concrete arch dam while preexisting joints are also simulated using discrete crack elements. Under those circumstances, it will be possible to evaluate the seismic safety of concrete arch dams in a more realistic manner. The comprehensive simulation including joint nonlinearity may need a mesh refinement to predict any sort of collapse mechanism. ACKNOWLEDGEMENTS

The supports of Iran’s National Elites Foundation and Ministry of Science, Research, and Technology of Iran are greatly appreciated. The writers also thank Professor Somasundaram Valliappan for insightful discussions during the first author’s visit to the University of New South Wales, Australia. REFERENCES 1. Chopra AK. Earthquake analysis of arch dams: factors to be considered. The 14th World Conference on Earthquake Engineering, October 12–17, 2008, Beijing, China. 2. Lee J, Fenves GL. A plastic–damage model for cyclic loading of concrete structures. Journal of Engineering Mechanics, ASCE 1998; 124:892–900. 3. Lubliner J, Oliver J, Oller S, Onate E. A plastic–damage model for concrete. International Journal of Solids and Structures 1989; 25:299–326. 4. Yazdani S, Schreyer HL. Combined plasticity and damage mechanics model for plain concrete. Journal of Engineering Mechanics, ASCE 1990; 116:1435–1450. 5. Abu-Lebdeh TH, Voyiadjis GZ. Plasticity–damage model for concrete under cyclic multiaxial loading. Journal of Engineering Mechanics, ASCE 1993; 119(7):1465–1484. 6. Meschke G, Lackner R, Mang HA. An anisotropic elastoplastic–damage model for plain concrete. International Journal for Numerical Methods in Engineering 1998; 42:703–727. 7. Gatuingt F, Pijaudier-Cabot G. Coupled damage and plasticity modeling in transient dynamic analysis of concrete. International Journal of Numerical and Analytical Methods in Geomechanics 2002; 26: 1–24. 8. Jason L, Pijaudier-Cabot G, Huerta A, Crouch R, Ghavamian S. An elastic plastic damage formulation for concrete: application to elementary tests and comparison with an isotropic damage model. Computer Methods in Applied Mechanics and Engineering 2006; 195(52):7077–7092. 9. Grassl P, Jirasek M. Damage–plastic model for concrete failure. International Journal of Solids and Structures 2006; 43:7166–7196. 10. Cicekli U, Voyiadjis GZ, Abu Al-Rub RK. A plasticity and anisotropic damage model for plain concrete. International Journal of Plasticity 2007; 23:1874–1900. 11. Cervenka J, Papanikolaou VK. Three dimensional combined fracture–plastic material model for concrete. International Journal of Plasticity 2007; 24:2192–2220. 12. Nguyen GD, Korsunsky AM. Development of an approach to constitutive modeling of concrete: isotropic damage coupled with plasticity. International Journal of Solids and Structures 2008; 45:5483–5501. 13. Abu Al-Rub RK, Voyiadjis GZ. Gradient-enhanced coupled plasticity-anisotropic damage model for concrete fracture: computational aspects and applications. International Journal of Damage Mechanics 2009; 18:115–154. 14. Vargas-Loli L, Fenves GL. Effects of concrete cracking on the earthquake response of gravity dams. Earthquake Engineering and Structural Dynamics 1989; 18:575–592. 15. El-Aidi B, Hall JF. Nonlinear earthquake response of concrete gravity dams, part 1: modeling; part 2: behavior. Earthquake Engineering and Structural Dynamics 1989; 18:837–865. 16. Bhattacharjee SS, Leger P. Seismic cracking and energy dissipation in concrete gravity dams. Earthquake Engineering and Structural Dynamics 1993; 22:991–1007. 17. Espandar R, Lotfi V. Application of the fixed smeared crack model in earthquake analysis of arch dams. Dam Engineering 2000; X(4):219–248. 18. Espandar R, Lotfi V. Comparison of non-orthogonal smeared crack and plasticity models for dynamic analysis of concrete arch dams. Computers and Structures 2003; 81:1461–1474. 19. Espandar R, Lotfi V, Razaqpur G. Influence of effective parameters of non-orthogonal smeared crack approach in seismic response of concrete arch dams. Canadian Journal of Civil Engineering 2004; 30:890–901. 20. Mirzabozorg H, Ghaemian M. Nonlinear behavior of mass concrete in three-dimensional problems using a smeared crack approach. Earthquake Engineering and Structural Dynamics 2005; 34:247–269. 21. Ghrib F, Tinawi R. An application of damage mechanics for seismic analysis of concrete gravity dams. Earthquake Engineering and Structural Dynamics 1995; 24:157–173. 22. Ghrib F, Tinawi R. Nonlinear behavior of concrete dams using damage mechanics. Journal of Engineering Mechanics, ASCE 1995; 121(4):513–526. 23. Cervera M, Oliver J, Faria R. Seismic evaluation of concrete dams via continuum damage models. Earthquake Engineering and Structural Dynamics 1995; 24:1225–1245. Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe

PLASTIC–DAMAGE ANALYSIS OF ARCH DAMS

2149

24. Cervera M, Oliver J, Manzoli O. A rate-dependent isotropic damage model for the seismic analysis of concrete dams. Earthquake Engineering and Structural Dynamics 1996; 25:987–1010. 25. Valliappan S, Yazdchi M, Khalili N. Earthquake analysis of gravity dams based on damage mechanics concept. International Journal for Numerical and Analytical Methods in Geomechanics 1996; 20:725–751. 26. Valliappan S, Yazdchi M, Khalili N. Seismic analysis of arch dams — a continuum damage mechanics approach. International Journal for Numerical Methods in Engineering 1999; 45:1695–1724. 27. Gunn RM. Nonlinear design and safety analysis of arch dams using damage mechanics, part 1: formulation. Hydropowers and Dams 2001; 2:67–74. 28. Gunn RM. Nonlinear design and safety analysis of arch dams using damage mechanics, part 2: applications. Hydropowers and Dams 2001; 3:72–80. 29. Ardakanian R, Ghaemian M, Mirzabozorg H. Nonlinear behavior of mass concrete in 3-D problems using damage mechanics approach. European Earthquake Engineering 2006; 2:89–65. 30. Alembagheri M, Ghaemian M. Seismic assessment of concrete gravity dams using capacity estimation and damage indexes. Earthquake Engineering and Structural Dynamics 2013; 42:123–144. 31. Alembagheri M, Ghaemian M. Damage assessment of a concrete arch dam through nonlinear incremental dynamic analysis. Soil Dynamics and Earthquake Engineering 2013; 44:127–137. 32. Lee J, Fenves GL. A plastic–damage concrete model for earthquake analysis of dams. Earthquake Engineering and Structural Dynamics 1998; 27:937–956. 33. Omidi O SNACS: a program for seismic nonlinear analysis of concrete structures. Department of Civil and Environmental Engineering, Amirkabir University of Technology, Tehran, Iran, 2010. 34. Ghaemian M, Ghobarah A. Staggered solution schemes for dam–reservoir interaction. Journal of Fluids and Structures 1998; 12:933–948. 35. Ghaemian M, Ghobarah A. Nonlinear seismic response of concrete gravity dams with dam–reservoir interaction. Engineering Structures 1999; 21:306–315. 36. Calayir Y, Karaton M. Seismic fracture analysis of concrete gravity dams including dam–reservoir interaction. Computers and Structures 2005; 83:1595–1606. 37. Omidi O, Lotfi V. Finite element analysis of concrete structures using plastic–damage model in 3-D implementation. International Journal of Civil Engineering 2010; 8(3):187–203. 38. Lee J, Fenves GL. A return-mapping algorithm for plastic–damage models: 3-D and plane stress formulation. International Journal for Numerical Methods in Engineering 2001; 50:487–506. 39. Omidi O, Lotfi V. Continuum large cracking in a rate-dependent plastic–damage model for cyclic-loaded concrete structures. International Journal for Numerical and Analytical Methods in Geomechanics 2013; 37:1363–1390. 40. Horii H, Chen SC. Computational fracture analysis of concrete gravity dams by crack-embedded elements: toward an engineering evaluation of seismic safety. Engineering Fracture Mechanics 2003; 70:1029–1045. 41. Omidi O, Valliappan S, Lotfi V. Seismic cracking of concrete gravity dams by plastic–damage model using different damping mechanisms. Finite Elements in Analysis and Design 2013; 63:80–97. 42. Wu JY, Li J, Faria R. An energy release rate-based plastic–damage model for concrete. International Journal of Solids and Structures 2006; 43:583–612. 43. Lotfi V. Seismic analysis of concrete dams using the pseudo-symmetric technique. Dam Engineering 2002; XIII(2): 119–145. 44. Simo JC. Algorithm for static and dynamics multiplicative plasticity that preserve the classical return-mapping schemes of the infinitesimal theory. Computer Methods in Applied Mechanics and Engineering 1991; 99:61–112. 45. Omidi O, Lotfi V. Application of pseudo-symmetric technique for seismic analysis of concrete arch dams. Proceedings of 4th International Conference on Fluid–Structure Interaction, the New Forest, UK, 2007; 153–162. 46. Chopra AK, Chakrabarti P, Gupta S. Earthquake response of concrete gravity dams including hydrodynamic and foundation interaction effects. Report No. UCB/EERC 80/01 1980, Earthquake Engineering Research Center, University of California, Berkeley, USA. 47. Guanglun W, Pekau OA, Chuhan Z, Shaomin W. Seismic fracture analysis of concrete gravity dams based on nonlinear fracture mechanics. Engineering Fracture Mechanics 2000; 65:67–87. 48. Chopra AK, Chakrabarti P. The Koyna earthquake and damage to Koyna dam. Bulletin of the Seismological Society of America 1973; 63:381–397. 49. Hall JF. The dynamic and earthquake behavior of concrete dams: review of experimental behavior and observational evidence. Soil Dynamics and Earthquake Engineering 1988; 7(2):57–121.

Copyright © 2013 John Wiley & Sons, Ltd.

Earthquake Engng Struct. Dyn. 2013; 42:2129–2149 DOI: 10.1002/eqe