Seismic fracture analysis of concrete gravity dams including dam–reservoir interaction

Computers and Structures 83 (2005) 1595–1606 www.elsevier.com/locate/compstruc

Seismic fracture analysis of concrete gravity dams including dam–reservoir interaction Yusuf Calayir *, Muhammet Karaton Fırat University, Faculty of Engineering, Department of Civil Engineering, 23279 Elazig, Turkey Received 28 January 2004; accepted 16 February 2005 Available online 13 April 2005

Abstract In this study, the seismic fracture response of concrete gravity dams is investigated with considering the effects of dam–reservoir interaction. A co-axial rotating crack model (CRCM), which includes the strain softening behavior, is selected for concrete material. The dynamic equilibrium equations of motion for the dam–reservoir system are solved by using the improved form of the HHT-a time integration algorithm. Two-dimensional seismic analyses of Koyna gravity dam are performed by using the 1967 Koyna earthquake records for numerical application. The effects of cracking on the seismic response of concrete gravity dams are discussed.
2005 Elsevier Ltd. All rights reserved. Keywords: Seismic fracture; Concrete gravity dam; Dam–reservoir interaction; Non-linear analysis

1. Introduction Many operating concrete dams present internal zones of extensive cracking or micro-cracking. The origin of cracking in concrete dams may be found in the phases of design, construction or in operation. The reasons for cracking may be many folds: thermal variations due to internal or external sources, seismic actions, irregular settlements of the foundation, internal chemical reactions during the setting, hardening and life of the concrete, etc. [1]. The cracking based on seismic actions in dams is generally more important and it may cause catastrophic consequences such as loss of life and property if the dam fails. Concrete dams are designed, in seismically active areas, to resist two levels of earthquakes:

*

Corresponding author. Fax: +96 424 241 5526. E-mail address: ycalayir@firat.edu.tr (Y. Calayir).

design base earthquake (DBE) as well as a maximum credible earthquake (MCE). The dam is expected to remain in the elastic range for the first level while a low annual probability of exceedance, the order of 10
3 or 10
4, is assumed for the second event [2,3]. Cracking and structural damage are expected under such conditions as long as the reservoir is contained. For concrete gravity dams designed according to current design criteria, the static and earthquake compressive stresses are generally much less than the compressive strength of the concrete. However, linear dynamic analyses of gravity dams show that the earthquake ground motion can produce tensile stresses that exceed the tensile strength of the mass concrete [4,5]. In such cases, a linear analysis is no longer valid because tensile cracks will form and propagate in the concrete, affecting the vibration properties and dynamic response of the dam. Even, water may fill inside the cracks on the dam upstream face and hence result in pressures in there, influencing crack

0045-7949/$ - see front matter
2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2005.02.003

1596

Y. Calayir, M. Karaton / Computers and Structures 83 (2005) 1595–1606

propagation as well as the behavior of dam. Consequently, tensile cracking is an important non-linear phenomenon to consider in the seismic response of concrete gravity dams. The effect of water pressures inside cracks is smaller than that of cracks only on the global modal properties and response of these dams [6]. Extensive research has been performed over the last 25 years into the fracture behavior of concrete. A tensile-strength-based crack propagation analysis is generally considered unreliable due to the mesh-dependent response prediction [7]. Non-linear behavior in the fracture process zone (FPZ), which is significantly large for the concrete normally used in dams, is neglected in the conventional linear elastic fracture mechanics (LEFM) models. Under very slowly applied loads and also under impulsive loads, concrete fracture behavior seems to be adequately predicted by the LEFM models [8,9]. In the intermediate range, from short-term static loading to seismic-induced strain rates, non-linear fracture mechanics (NLFM) models considering the strain softening behavior in the FPZ appear to be more appropriate. Continuum mechanics approaches to represent the tensile crack propagation are computationally very attractive for applications in complex structural analysis when the crack profiles are not known a priori. General drawbacks of the conventional smeared crack analysis model, such as the mesh-sensitive response prediction and the stress-locking in finite elements, can largely be overcome using the energy conserving co-axial rotating crack model (CRCM), which is an improved NLFM model [8,10]. Bhattacharjee and Le´ger [8] investigated the dynamic fracture response of Koyna dam by using the CRCM concrete model. In the analysis of this dam, the foundation was assumed to be rigid, and the dynamic interaction of dam–reservoir system was modeled by the Westergaard added mass technique. The response results of the dam indicated that cracks formed at the base and at the upper portion near the change in downstream slope for the Koyna earthquake records. Lee and Fenves [11] also obtained similar solutions for Koyna dam under the same boundary and loading conditions by using a plastic-damage concrete model which includes tension and compression strain softening effects. They studied in detail tension damage in the dam and noted that damage in compression is less important in comparison with damage in tension for dams. Guanglun et al. [12] used a numerical scheme based on non-linear crack band theory for tension cracks to investigate the seismic fracture behavior of Koyna dam under the Koyna earthquake accelerograms. This study does not take into account the dynamic interaction of dam–reservoir system, only including the static water pressures of reservoir. The analysis of dam with rigid foundation showed cracks at the upper portion of the dam near the change in

downstream slope. In another paper, Cervera et al. [13] developed an isotropic damage model for the seismic evaluation of concrete gravity and arch dams. This model accounts for the different behavior of concrete in tension and compression by splitting the stress tensor into tensile and compressive components, each with its own damage surface and evolution law. The seismic analyses of their selected gravity dam which resembles very closely Koyna dam, including dam–reservoir–foundation interactions, subjected to artificially generated earthquakes of different intensities indicated tension damage in the dam, but no compression damage occurred. The formed tension damage was seen at the upper portion of the dam. In the lights of these investigations, the tension cracks have only important influence on the seismic response of concrete gravity dams. Therefore, tension cracking is only considered in the present study. Dam–foundation and dam–reservoir interactions are two important aspects in the dynamic analysis of dams. The first one can be considered by simplified or rigorous methods. Solutions depending on the simplified methods are approximate. More realistic results may be obtained by using rigorous methods, which require finite element or boundary element modeling of foundation domain and result in increasing computation time and computer storage requirement. Many DOFs are not desired in the non-linear analysis of dams because of the necessities of an iterative solution and using very small time integration step. However, dam–foundation rock interaction effects in the dynamic response of dams are strongly depended on the elastic modulus ratio of foundation to dam. The solutions including flexible foundation approximates to those of the infinitely rigid foundation case as this ratio increases [14,15]. Therefore, the foundation can be modeled as being rigid if the foundation rock is much stiffer than the dam concrete, and also if the earthquake acceleration components are recorded near the dam foundation [12]. The latter aspect, dam– reservoir interaction, may be considered by three ways: added mass, Eulerian and Lagrangian approaches [16]. Added mass approach is the most simple and first technique and it employs the results of rigid dam–incompressible fluid interaction. More realistic solutions for dam–reservoir interaction can be obtained by using Eulerian and Lagrangian approaches. In the Eulerian approach, displacements and pressures are the variables in the structure and reservoir domains, respectively, while displacements are variables in the both domains for the Lagrangian approach. In the finite element analyses of dam–reservoir systems using these approaches, a radiating condition (for example, Sommerfeld radiation condition) may be applied to the truncated boundary of reservoir as one way for considering of radiating damping. However, the effect of the radiating condition on the solutions is generally negligible if the reservoir length is

Y. Calayir, M. Karaton / Computers and Structures 83 (2005) 1595–1606

1597

taken three or more times of its height [17]. Many researchers investigated the dynamic interaction of dam–reservoir systems taking the reservoir length as three times of its height [16–18]. In this paper, the seismic fracture response of the concrete gravity dams is investigated with considering the effects of dam–reservoir interaction. A CRCM which includes the strain softening behavior for tension is selected for the concrete material model. Lagrangian approach is used for the finite element modeling of dam–reservoir interaction problem. The reservoir domain is considered to be linearly elastic. For numerical application, two dimensional seismic analyses of Koyna gravity dam are performed by using the transverse and vertical acceleration components of 11 December 1967 Koyna earthquake records. A plane stress finite element idealization is adopted for the dam monoliths, which vibrate independently each other in the transverse plane because there was evidence of relative movement between adjacent monoliths of the dam during 11 December 1967 Koyna earthquake. The effects of cracking on seismic response of the concrete gravity dam are studied.

tion of the stress resistance, together with an increase in deformation. This material behavior is called strain softening [20]. Hence, the stress–strain relationship for the concrete like materials is divided into two regions: (1) from 0 up to ft and (2) a softening region. A widely used assumption in numerical modeling of mass concrete is to adopt a triangular stress–strain diagram for uniaxial loading. This gives a linear strain softening relationship. But various experimental evidences indicate that it is more realistic to assume a strain softening curve with a steep initial decline followed by an extended tail [19]. Thus, using an exponential curve for strain softening model is more appropriate than linear one. In present study, the behavior of dam concrete is represented by the co-axial rotating crack model (CRCM) including the exponential strain softening of concrete. According to Fig. 2(a), an exponential strain softening model yields

2. A non-linear smeared crack model for concrete material

rðeÞ ¼ 0:0;

As is well known, concrete and geo-materials eventually exhibit strain softening in tension, leading to a complete loss of strength. In these materials, the secant modulus decreases with increasing strain [19]. A typical stress–strain relationship for mass concrete from uniaxial tension test is shown in Fig. 1. At the beginning, a linear relationship between stress and strain exists for up to 60% of the maximum stress. Then, micro-cracks develop within the specimen, and they cause the non-linearity in the curve up to the tensile strength (ft). In the post-peak regime, more micro-cracks are developed in the weakest cross-section of the specimen (fracture process zone, FPZ) and coalescence of them causes a gradual reduc-

where ft is the tensile strength, e0 is the corresponding strain threshold, E0 is the modulus of elasticity and a is a dimensionless constant which can be determined from the following relation [19,21]:

rðeÞ ¼ E0 e;

a¼ e0

h

3 2E0 Gf lch ft2

Strees (MPa)

e0 < e < ecr ;

i P 0:0;
1

ð2Þ

where the characteristic length lch is a geometrical constant which is introduced as a measure of the length of the FPZ in a sample, and is equal to the length of the element normal to the cracks. Gf is the fracture energy per unit area (Fig. 2(a)).

σ ft lch

FPZ {

1.5 1.0

σ

0.5 εo

0.0 0.0

0.1

0.2

0.3

ð1bÞ ð1cÞ

e P ecr ;

2.5 2.0

ð1aÞ

e 6 e0 ;


 rðeÞ ¼ ft 2e
aðe
e0 Þ
e
2aðe
e0 Þ ;

0.5

0.4

0.6

0.7

0.8

0.9

1.0

-2

Strain x 10

Fig. 1. Typical stress–strain curve for mass concrete from simple tension test [20].

1598

Y. Calayir, M. Karaton / Computers and Structures 83 (2005) 1595–1606

(a)

(b) σ Softening initation

y

ft

s

Unloading Reloading

n

Eo

θ

Completely Fractured

En εo Reopening

x

εmax

(c)

ε

εcr

(d)

Fig. 2. Constitutive modeling for non-linear smeared crack analysis. (a) Idealized stress–strain curve for concrete; (b) biaxial stress effect in the strain softening initiation; (c) local axis system; and (d) closing and reopening of cracks.

In the relations given with Eq. (1), a maximum strain ecr has been adopted that may not be exceeded in strain softening, and is consistent with the study carried out by Bazant [22]. In this study, the value of ecr is calculated when its corresponding stress is equal to 0.02ft, which is a reasonable value. Then, ecr can be expressed as h pffiffiffiffiffiffii ln 1þ d1
d : ð3Þ ecr ¼ e0 þ a When d = 0.02, Eq. (3) gives ecr ¼ e0 þ

4:6 : a

ð4Þ

The stresses and strains at a material point are related as frg ¼ ½Dfeg;

ð5Þ

where [D] is the constitutive matrix, and {r} and {e} are the stress and strain vectors, respectively. If there are no cracks, the standard elastic plane stress–strain matrix for an isotropic material is used. Then, [D] is given by

2

1 m E0 6 m 1 ½D ¼ 4 1
m2 0 0

3 0 0 7 5;

ð6Þ

ð1
mÞ 2

where m is the PoissonÕs ratio. Softening initiation is determined using the following material index U0 ¼

ft e0 ft2 ¼ ; 2 2E

ð7Þ

which represents approximately the area under the uniaxial stress–strain curve up to the tensile strength. In the 2D finite element analyses, the strain softening is assumed to initiate when the tensile strain energy density (r1e1/2) is greater or equal to the material index U0. r1 and e1 are the maximum principal stress and strain, respectively. Hence, the biaxial effect in the strain softening initiation criterion may be given by r1 P

ft2 : Ee1

ð8Þ

This equation represents a biaxial failure envelope schematically given in Fig. 2(b). As seen in this figure, the

Y. Calayir, M. Karaton / Computers and Structures 83 (2005) 1595–1606

failure occurs for tension–tension and tension–compression stress states due to tension cracking. However, no failure forms in compression–compression stress state. After the softening initiation, a smeared band of micro-cracks is assumed to appear in the direction perpendicular to the principal tensile strain. The material reference axis system, referred to as the local axis system, is aligned with the principal strain directions [direction n-s in Fig. 2(c)]. The constitutive matrix in the local axis system may be written as 2 E0 g 3 E0 gm 0 1
gm2 1
gm2 6 E0 gm 7 En E0 0 7 ½Dns ¼ 6 ð9Þ 4 1
gm2 1
gm2 5; g ¼ E 0 ; E0 0 0 l 2ð1þmÞ where g (0 6 g 6 1) is the ratio between the softened elastic modulus (En) in the direction normal to a fracture plane and the initial isotropic elastic modulus (E0), and l is the shear resistance factor and is defined for the CRCM as follows [8]:   1 þ m gen
es l¼
gm ð0 6 l 6 1Þ; ð10Þ 1
gm2 en
es where en and es are the normal strain components in perpendicular and parallel directions to the fracture plane, respectively. [D]ns matrix obtained in the local axis system can be transformed to the global co-ordinate directions, [D]xy, using the standard transformation procedure. Under reversal loading cases, closing and reopening of cracks occur. In these cases, the tensile strain en in an element may alternatively increase or decrease. With increasing strain softening, the softened elastic modulus En (Fig. 2(d)) and hence the parameters g and l decrease gradually and may eventually reach zero values after complete fracture (en P ecr). The shear resistant factor l increases with the decreasing of en. The crack is assumed to be closed if the parameter l is greater than a threshold value lc, which may be taken as 0.95 [8]. In this case, the softened elastic modulus En is replaced by the initial elastic modulus E0. The residual strain upon closing of a crack is given by en = mes. Fig. 2(d) shows the closing/reopening behavior for a special case when es  0. The secant modulus En remains unchanged during unloading/reloading when the strain en is less than the previously attained maximum value emax. However, the parameter l changes during this process. In the completely fractured case, closing/reopening path goes on the axis e (Fig. 2(d)). While the values of g and l parameters change, the matrix [D]ns is updated for all cases mentioned above. In finite element implementation, finite element model of the dam is constructed using four-node isoparametric elements. 2 · 2 Gaussian integration rule is used to compute the element stiffness matrices. Strains are computed at each integration point and the average

1599

of the Gaussian point strains is taken as representative of the behavior of the element as a whole [23]. Principal strains computed from the average strain response of an element are used for determination of the material parameters. The parameter of lch, which is used in Eq. (2), is calculated from square root of related element area. The constitutive matrix, [D]xy, in the post-elastic state represent the average resistance of an element. The stresses at individual Gaussian points are computed from the respective total strains using the most recent [D]xy matrix. The element stiffness matrix is also updated using the same [D]xy. A finite element solution initially starts with the elastic stiffness matrix for the un-cracked elements. Among all candidate elements that would initiate softening at a particular iteration, the one with the highest tensile strain energy density (r1e1/2) is used first. The fracture propagation in the structure is achieved by using one new softening element per iteration.

3. Lagrangian formulation for dynamic interaction of fluid–structure systems In the Lagrangian approach, displacements are selected as the variables in both fluid and structure domains [16,24]. Fluid is assumed to be linearly elastic, inviscid and irrotational. Two-dimensional stress–strain relationships of the fluid undergoing small amplitude motion are given by      P b 0 ev ¼ ; ð11Þ Pz 0 az Wz where P is the pressure (tension has a positive sign), b is the bulk modulus of fluid and
v is the volumetric strain. Wz is rotation about axis z, and Pz and az are the stress and constraint parameter related with Wz, respectively. Including the small amplitude free-surface waves (sloshing waves) cause pressure at the free-surface of fluid. This pressure can be given by P ¼
cw ufn ;

ð12Þ

where cw is the weight density of fluid and ufn is the normal component of the free-surface displacement. The free surface stiffness for fluid is obtained from the discrete form of Eq. (12). Finite element equations of motion for the fluid system can be expressed as [16] ½M f faf g þ ½K f fuf g ¼ fF lf g

ð13Þ

where {af} and {uf} are the nodal point vectors of acceleration and displacement for the fluid system, respectively. [Kf] is the stiffness matrix including the free surface stiffness, [Mf] the mass matrix, and fF lf g the load vector for the fluid system. Reduced integration orders are utilized in the formation of the fluid element matrices [16]. In this study, one-point integration is applied to four-node plane element used in the finite model of fluid.

1600

Y. Calayir, M. Karaton / Computers and Structures 83 (2005) 1595–1606

For obtaining the coupled equations of the fluid– structure system, it is required to determine the interface condition. Because the fluid is assumed to be inviscid, only the displacement in normal direction to the interface is continuous at the interface of the system. This condition can be imposed by the penalty method [16,24]. Using the interface condition, the equations of motion of the coupled system may be given as ½M c fac g þ ½C s fvs g þ ½K f fuf g þ fF is g ¼ fF lc g;

ð14Þ

where [Mc] is the mass matrix of the coupled system and [Cs] is the damping matrix of the structure. {ac} and fF lc g are the relative acceleration and the load vectors for the coupled system, respectively. {vs} and fF is g are the relative velocity and the restoring force vectors for the structure system, respectively. The damping matrix [Cs] is defined to be stiffnessproportional as ½C s  ¼ bd ½K s ;

4. Numerical solution of dynamic equilibrium equations Dynamic equilibrium equations for fluid–structure systems are given by Eq. (14) in the preceding section. Solution of this equation is carried out in time domain using the improved form of HHT-a integration method presented by Miranda et al. [26]. The integration method retains the Newmark difference formulas: uiþ1 ¼ ui þ Dtvi þ 12Dt2 ½ð1
2bÞai þ 2baiþ1 ;

ð17aÞ

viþ1 ¼ vi þ Dt½ð1
cÞai þ caiþ1 ;

ð17bÞ

90.45 m.

ð15Þ

where bd is determined by specifying a desired damping ratio at a given frequency, and [Ks] is the structural stiffness matrix. A mass-proportional term is omitted because it would provide some artificial stability to the portion of the dam above a crack [25]. In the case of earthquake ground motion, the load vector fF lc g can be defined as ¼

fF stat c g


½M c fag g;

0.5 0.4

ð16Þ

where fF stat c g is the static load vector of coupled system and {ag} is the acceleration vector of ground motion for the same system.

1029

Fig. 4. Finite element model of reservoir.

14.80 m.

DETAIL A

0.2 0.1 0 -0.1 -0.2

28

-0.3

28

-0.4

1/6.54

-0.5

(a)

DETAIL A 28 28

0

1

2

3

4

5

6

7

8

Time (s) 0.4 Vertical Component

0.3

1/1.38 66.50 m. 1/24

Acceletration, g

103.00 m.

Transverse Component

0.3

Acceletration, g

fF lc g

271.35 m.

39.00 m.

0.2 0.1 0 -0.1 -0.2 -0.3 -0.4

913

(b)

0

1

2

3

4

5

6

7

8

Time (s)

70.20 m.

Fig. 3. Finite element model of Koyna dam.

Fig. 5. The Koyna accelerograms. (a) Transverse component and (b) vertical component.

Y. Calayir, M. Karaton / Computers and Structures 83 (2005) 1595–1606

the discrete equations in the time domain, Eq. (14), are modified using ½M c fac giþ1 þ ð1 þ aÞ½C s fvs giþ1
a½C s fvs gi þ ð1 þ aÞ½K f fuf giþ1
a½K f fuf gi þ ð1 þ aÞfF is giþ1
afF is gi ¼ ð1 þ aÞfF lc giþ1
afF lc gi ;

ð18Þ

where b and c are NewmarkÕs coefficients and a is a parameter controlling the numerical dissipation. To ensure second-order accuracy and unconditional stability, these parameters should be chosen such that
 ð19Þ a 2
13; 0 ; b ¼ 14ð1
aÞ2 ; c ¼ 12
a: In this study, a is selected as
0.10. For this selection, the values of b and c are being 0.605 and 0.60, respectively. The predictor-corrector technique is used in conjunction with the Newton–Raphson method to solve the non-linear dynamic equation of fluid–structure interaction problem. The improved form of HHT-a time integration algorithm is given in Ref. [26]. If the solution at the time step i is known, the predicted displacement and velocity vectors for the time step i + 1 can be calculated as f~ uc giþ1 ¼ fuc gi þ Dtfvc gi þ 12Dt2 ð1
2bÞfac gi ;

ð20aÞ

f~vc giþ1 ¼ fvc gi þ Dtð1
cÞfac gi ;

ð20bÞ

where subscript c indicates that the quantity is related to the coupled system. Substitution of the last two equations into Eqs. (17a) and (17b) leads the following relationships: fuc giþ1 ¼ f~ uc giþ1 þ bDt2 fac giþ1 ;

ð21aÞ

fvc giþ1 ¼ f~vc giþ1 þ Dtcfac giþ1 :

ð21bÞ

These quantities are substituted into Eq. (18) and then time marching algorithm is applied to Eq. (18) as given in Appendix A. 5. Response of a concrete gravity dam Koyna concrete gravity dam, which is extensively analyzed in several previous studies [2,8,11,12,27], in India, is selected for numerical application. This dam is one of a few concrete dams that have experienced a destructive earthquake. The earthquake of December 11, 1967, with maximum acceleration around 0.5g, caused significant structural damage to the dam, including horizontal cracks on the upstream and downstream faces of a number of non-overflow monoliths around the elevation at which the slope of the downstream face changes abruptly [28]. In this section, the non-linear dynamic fracture analysis of Koyna concrete gravity dam is performed by

-1 1

-2 2

2

-2 -4 -3

3

3

5 7

4 5

8

-8

-7 5

2

1

4

-6

-2

-7

4

-2

3

-6

2

3

-8

6

1

-4 -3

-4-3

-5

-5

2 1

-4 1

1 3

2

-4

5 4

(a)

1601

-3

-3

(b)

Fig. 6. Envelopes of maximum (a) and minimum (b) principal stresses of Koyna dam for linear dynamic analysis case.

1602

Y. Calayir, M. Karaton / Computers and Structures 83 (2005) 1595–1606

considering dam–reservoir interaction. Finite element model for the tallest section of the dam is shown in Fig. 3. Two nodal points and two elements, for which the time history graphs of the response quantities are plotted, are marked on the mesh. The foundation is modeled as being rigid since the foundation rock of the dam is basalt with an elastic modulus of 7 · 104 MPa which is about 2.25 times of that of the dam concrete, and also the earthquake acceleration components are recorded near the dam foundation [12]. The fluid domain (reservoir) is assumed to be linearly elastic and the length of reservoir is considered as three times of its height. Finite element model for the reservoir system is given in Fig. 4. The displacements of the nodal points in the normal direction at the bottom and the truncated boundary of the reservoir are assumed to be zero. The following data is considered for the fluid domain: the bulk modulus is 2070 MPa, the mass density is 1000 kg/m3, and the rotation constraint parameter is chosen as 100 times of the bulk modulus. The values of the material parameters for the dam concrete are selected as follows: The elastic modulus is 31027 MPa, the PoissonÕs ratio is 0.2, the mass density is 2643 kg/ m3, the tensile strength is 1.5 MPa and the fracture

energy is 150 N/m. Dynamic loading affects the concrete material parameters. Elastic modulus of concrete is generally considered to be less sensitive to strain rate than the tensile strength and fracture energy [8,29]. Because of strain rate effects, the tensile strength and the fracture energy are increased by 20% approximately, leading to the values of 1.8 MPa and 180 N/m, respectively. The variation of the elastic modulus due to dynamic loading is neglected. The damping in the dam is assumed as being proportional to the stiffness, and it provides a critical viscous damping ratio of 5% in the fundamental vibration mode of the dam alone with no cracking. The damping forces are then computed as proportional to the tangent stiffness matrix for the dam. Integration time step is selected as 0.001 s. For dynamic input, the transverse and vertical acceleration components of 11 December 1967 Koyna earthquake are selected, and assumed that they are acting in the stream and vertical directions of the dam–reservoir system, respectively. The records of this ground motion are given in Fig. 5. The linear and non-linear solutions obtained from the dynamic analysis of the dam–reservoir system are compared with each other. The cracking effects of the concrete material on

20

60 Non-linear Linear

20 0 -20 -40 -60

(a)

5 0 -5

-15

0

1

2

3

4

5

6

7

8

0

1

2

3

Time (s)

5

6

7

8

4

5

6

7

8

Non-linear Linear

15

10 5 0 -5 -10

10 5 0 -5

-15 -20 0

4

Time (s)

20 Non-linear Linear

Displacement (mm)

15

Displacement (mm)

10

-10

20

(b)

Non-linear Linear

15

Displacement (mm)

Displacement (mm)

40

1

2

3

4

5

6

7

8

Time (s)

Fig. 7. Time history graphs of the horizontal (a) and vertical (b) displacements of the nodal point 1029 of the dam.

-10

0

1

2

3

Time (s)

Fig. 8. Time history graphs of the horizontal (a) and vertical (b) displacements of the nodal point 28 of the dam.

Y. Calayir, M. Karaton / Computers and Structures 83 (2005) 1595–1606

1603

dynamic analysis case are given in Fig. 6. The numeric values marked on these envelopes are given in terms of MN/m2. The response results represent larger tensile stresses in the upper part of the dam, especially around the elevation at which the slope of the downstream face changes abruptly. These stresses, which exceed 5 MPa on the upstream face and 8 MPa on the downstream face, are approximately three to five times the tensile strength selected for the concrete. Tensile stresses around at the heel also exceed 6 MPa. Compressive stresses also take larger values in the upper part of the

the seismic response of the concrete gravity dam are studied. A computer program written in Fortran 90 language by the authors for the linear and non-linear dynamic analyses of dam–reservoir systems is used for the solutions. The static solutions of the dam–reservoir system due to its gravity loads are taken as initial conditions in the dynamic analyses of the system. For previously predicting the regions where cracking may occur depending on the level of earthquake ground motion, the envelopes of maximum and minimum principal stresses for linear

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 9. Cracking cases in the dam at six selected times. (a) t = 2.156 s; (b) t = 2.846 s; (c) t = 3.698 s; (d) t = 4.092 s; (e) t = 4.142 s; and (f) t = 4.536 s.

Y. Calayir, M. Karaton / Computers and Structures 83 (2005) 1595–1606

faces are not under tensile stress at the same time. The occurring crack profiles are compatible with the results of previous researches for the same dam [2,8,11,12,27]. The time history graphs of the maximum and minimum principal stresses occurred in center of the element 28 on the downstream face of the dam located at the slope change are given in Fig. 10. While the maximum principal stresses for the linear case take larger peak values, the maximum peak value of those for the non-linear case are about the tensile strength of the concrete. Fig. 10(a) confirms that the tensile strength is completely removed after cracking. The minimum principal stresses for the linear case also generally take larger peak values except that the maximum peak value for the non-linear solution is some larger than that of the linear case. The maximum peak values of the minimum principal stresses for the linear and non-linear cases are 11.89 and 12.74 MPa, respectively. Fig. 11 shows the time history graphs of the maximum and minimum principal stresses occurred in center of the element 913 located at the heel of the dam. The variations in these stresses are similar those of the element 28. Cracking in the dam causes redistribution of the stresses.

14

Maximum principal stress (MPa)

dam, especially around the elevation at which the slope of the downstream face changes abruptly. The maximum compressive stress in the monolith (Fig. 6(b)) is little larger than 12 MPa, which is smaller than the compressive strength of concrete used. In the regions where larger tensile stresses occur, cracking in the concrete is expected due to the effect of the Koyna earthquake accelerograms. Fig. 7 shows the time history graphs of the horizontal and vertical displacements of the nodal point 1029 located at the dam crest when cracking in the dam is allowed and prevented (linear response). Also, Fig. 8 illustrates those of the nodal point 28 located at the discontinuity in the slope of the downstream face. There is no cracking during the relatively small amplitude motion. Cracking in the dam is firstly observed at 1.862 s near the discontinuity in the slope of the downstream face. At this location, stresses are concentrated and the tensile stresses reach the largest value. The displacement time histories obtained from the linear and non-linear solutions separate with each other as the cracks propagate in the dam (Figs. 7 and 8). The differences in the displacement amplitudes, which are initially small, reach important levels (especially between t = 3.8 s and 7.0 s) as the cracks propagate deeper inside of the dam. Maximum peak value of the vertical displacements of the nodal point 28 at the time of 4.117 s corresponds to the maximum opening state of the crack at the downstream face. Crack propagation also changes the vibration period of the dam–reservoir system. Because of the cracking effect, especially elongation of the vibration period is quite evident between times 4.0 and 5.0 s for the crest displacements (Fig. 7). For the cracking model, the crest displacements are dominated by rigid-body rocking of the upper portion of the dam after the formation of cracks near the change in downstream slope. At all times the non-linear response analysis shows that the dam remains stable. The cracking cases in the dam at six selected times are shown in Fig. 9. The state of an open crack in a finite element is indicated by shading the entire element area. Closed or incomplete crack conditions are represented by dots in the element centers. Elements that have never softened are unmarked in the presentations. Due to the infinite rigidity of the foundation, a stress concentration induces a crack propagation at the base of the dam. The cracks at the downstream face are initially in a horizontal region (Fig. 9(a) and (b)) propagate deeper inside of the dam and a crack band occurs which curves down due to the compressive stresses resulting from rocking of the top block (Fig. 9(c)–(f)). At the same time, previously initiated cracks in the heel of the dam propagate from the heel to the toe. The softened and complete cracked elements in the crack bands can be closed and reopened at different times as seen in Fig. 9. The dam retains its overall stability because the upstream and downstream

Non-linear Linear

12 10 8 6 4 2 0 -2

(a)

0

1

2

3

4

5

6

(b)

7

8

Time (sec) 2

Minimum principal stress (MPa)

1604

Non-linear Linear

0 -2 -4 -6 -8 -10 -12 -14

0

1

2

3

4

5

6

7

8

Time (sec)

Fig. 10. Time history graphs of the maximum (a) and minimum (b) principal stresses occurred in center of the element 28.

Y. Calayir, M. Karaton / Computers and Structures 83 (2005) 1595–1606

Maximum principal stress (MPa)

7

(a)

Non-linear Linear

6 5 4 3 2 1 0 -1 0

1

2

3

4

5

6

7

8

Time (sec)

1605

ences, which are initially small, reach important levels while the cracks propagate deeper inside of the dam. Cracking also changes the vibration period of the dam–reservoir system. The dynamic instability was not seen in spite of much cracks occurred in the dam. The selected time marching algorithm is based on its ability to damp out high-frequency noise resulting from sudden changes in the stiffness due to closing of the cracks. The co-axial rotating crack model selected for the concrete may be used successfully in the earthquake response of concrete gravity dam–reservoir systems in conjunction with the improved form of the HHT-a algorithm for time marching.

Minimum principal stress (MPa)

1

(b)

Non-linear Linear

Appendix A. Time marching algorithm using fracture mechanics

0 -1

The algorithm applied to Eq. (18) is as follows:

-2

1. Compute integration parameters

-3

A1 ¼

-4 -5 0

1

2

3

4

5

6

7

8

Time (sec)

Fig. 11. Time history graphs of the maximum (a) and minimum (b) principal stresses occurred in center of the element 913.

6. Conclusions In this study, the seismic fracture response of concrete gravity dams is investigated with considering the effects of dam–reservoir interaction. A co-axial rotating crack model (CRCM) which includes the strain softening behavior is selected for the concrete. Lagrangian approach is used for the finite element modeling of the dam–reservoir interaction problem. Viscous damping in the dam is assumed to be proportional to the tangent stiffness of the dam. The dynamic equilibrium equations of motion for this problem are solved by using the improved form of the HHT-a time integration algorithm. Two dimensional seismic analysis of Koyna gravity dam are performed by using the 1967 Koyna earthquake records for numerical application. The effects of cracking on seismic response of the dam are discussed. Cracks tend to initiate near stress concentrations in the dam monolith, mainly at the base and near the changes in the slope of the faces. Cracking in the Koyna dam is firstly observed at the location where the slope in downstream face of the dam changes. Displacements for the linear and non-linear solutions become different with each other as cracks propagate in the dam. The differ-

1 ; bDt2

A2 ¼

c : bDt

ð22Þ

2. {uc}i, {vc}i and {ac}i is known; set iteration counter, n = 1. 3. Predict response at ti+1, uc giþ1 ; fuc gniþ1 ¼ f~

ð23aÞ

fvc gniþ1 ¼ f~vc giþ1 ;

ð23bÞ

fac gniþ1 ¼ f0g:

ð23cÞ

4. Compute the effective dynamic stiffness matrix and the vector of unbalanced forces b n ¼ A1 ½M c  þ ð1 þ aÞA2 ½C s n þ ð1 þ aÞ ½K iþ1 iþ1  ½K f  þ ð1 þ aÞ½K s niþ1 ;

ð24Þ

fDF gniþ1 ¼ ð1 þ aÞfF lc giþ1
afF lc gi
ð1 þ aÞ  fF s gniþ1 þ afF s gi
½M c fac gniþ1
ð1 þ aÞ½C s niþ1 fvs gniþ1 þ a½C s i fvs gi
ð1 þ aÞ½K f fuf gniþ1 þ a½K f fuf gi ; ½K s niþ1

ð25Þ

fF s gniþ1

and are the tangent stiffness mawhere trix and the resisting force vector of the structure system, respectively, which are assembled from the element contributions. 5. Solve for incremental displacements, b n fDuc gn ¼ fDF gn : ½K iþ1 iþ1 iþ1

ð26Þ

6. Update displacement, velocity and acceleration vectors, n n fuc gnþ1 iþ1 ¼ fuc giþ1 þ fDuc giþ1 ;

ð27aÞ

1606

Y. Calayir, M. Karaton / Computers and Structures 83 (2005) 1595–1606

 fvc gnþ1 vc giþ1 þ A2 fuc gniþ1
f~uc giþ1 ; iþ1 ¼ f~

ð27bÞ

 n fac gnþ1 uc giþ1 : iþ1 ¼ A1 fuc giþ1
f~

ð27cÞ

7. Check for convergence of iteration process fDF gn using a force tolerance (Tol). iþ1 P (a) If fDF gniþ1 = nm¼1 fDF gmiþ1 6 Tol, convergence is achieved. Set fuc giþ1 ¼ fuc gnþ1 iþ1 , nþ1 fvc giþ1 ¼ fvc gnþ1 and fa g ¼ fa g . c iþ1 c iþ1 iþ1 (b) If convergence is not achieved; set n = n + 1 and return step 4. 8. Set i = i + 1 and return step 2.

References [1] Cervera M, Oliver J, Herrero E, Onate E. A computational model for progressive cracking in large dams due to the swelling of concrete. Eng Fract Mech 1990;35(1–3): 573–85. [2] Ghrib F, Tinawi R. An application of damage mechanics for seismic analysis of concrete gravity dams. Earthquake Eng Struct Dynam 1995;24(2):157–73. [3] NCR. Earthquake engineering of concrete dams-design, performance, and research needs. National Academy Press; 1990. [4] Vargas-Loli LM, Fenves GL. Effects of concrete cracking on the earthquake response of gravity dams. Earthquake Eng Struct Dynam 1989;18(4):575–92. [5] Hall JF. Study of the earthquake response of Pine Flat dam. Earthquake Eng Struct Dynam 1986;14(2):281–95. [6] Tinawi R, Guizani L. Formulation of hydrodynamic pressures in cracks due to earthquakes in concrete dams. Earthquake Eng Struct Dynam 1994;23(7):699–715. [7] Bazant ZP, Cedolin L. Blunt crack band propagation in finite element analysis. J Eng Mech Div ASCE 1979; 105:297–315. [8] Bhattacharjee SS, Le´ger P. Seismic cracking and energy dissipation in concrete gravity dams. Earthquake Eng Struct Dynam 1993;22(11):991–1007. [9] Du J, Yon JH, Hawkins NM, Arakawa K, Kobayashi AS. Fracture process zone for concrete for dynamic loading. Mater J ACI 1992;89:252–8. [10] Rots JG. Smeared and discrete representations of localized fracture. Int J Fract 1991;51:45–59. [11] Lee J, Fenves GL. A plastic-damage concrete model for earthquake analysis of dams. Earthquake Eng Struct Dynam 1998;27(9):937–56. [12] Guanglun W, Pekau OA, Chuhan Z, Shaomin W. Seismic fracture analysis of concrete gravity dams based on nonlinear fracture mechanics. Eng Fract Mech 2000; 65(1):67–87. [13] Cervera M, Oliver J, Faria R. Seismic evaluation of concrete dams via continuum damage models. Earthquake

Eng Struct Dynam 1995;24(9):1225–45. [14] Priscu R, Popovici A, Stematiu D, Stere C. Earthquake engineering for large dams. Bucuresßti: Editura Academiei; 1985. [15] Chopra AK, Gupta S. Hydrodynamic and foundation interaction effects in frequency response functions for concrete gravity dams. Earthquake Eng Struct Dynam 1982;10(1):89–106. [16] Calayır Y, Dumanog˘lu AA, Bayraktar A. Earthquake analysis of gravity dam–reservoir systems using the Eulerian and Lagrangian approaches. Comput Struct 1996; 59(5):877–90. [17] Dumanog˘lu AA, Calayır Y, Karaton M. Earthquake analysis of concrete gravity dams including reservoir and foundation interaction effects by the Eulerian approach. _ Technical Report TDV/TR 047-82, Istanbul: Turkish Earthquake Foundation, 2003 [in Turkish]. [18] Greeves EJ. The modelling and analysis of linear and nonlinear fluid–structure systems with particular reference to concrete dams. PhD thesis, University of Bristol, UK, 1991. [19] Lubliner J, Oliver J, Oller S, Onate E. A plastic damage model for concrete. Int J Solids Struct 1989;25(3):299–326. [20] Bruhwiler E. Fracture of mass concrete under simulated seismic action. Dam Eng 1990;I(3):153–75. [21] Valliappan S, Yazdchi M, Khalili N. Seismic analysis of arch dams-A continuum damage mechanics approach. Int J Numer Meth Eng 1999;45(11):1695–724. [22] Bazant ZP. Snapback instability at crack ligament tearing and its implication for fracture micromechanics. Cement Concr Res 1987;17(6):951–67. [23] Hinton E, Campbell JS. Local and global smoothing of discontinuous finite element functions using a least squares method. Int J Numer Meth Eng 1974;8:461–80. [24] Olson LG, Bathe KJ. A study of displacement-based fluid finite elements for calculating frequencies of fluid and fluid-structure systems. Nucl Eng Des 1983;76:137–51. [25] El-Aidi B, Hall JF. Non-linear earthquake response of concrete gravity dams Part 1: Modelling. Earthquake Eng Struct Dynam 1989;18(6):837–51. [26] Miranda I, Ferencz RM, Hughes TJR. An improved implicit-explicit time integration method for structural dynamics. Earthquake Eng Struct Dynam 1989;18(5): 643–53. [27] Calayır Y, Karaton M. Non-linear dynamic cracking analysis of concrete gravity dams. In: Pala S, Ilki A, editors. ACE2002 Advances in civil engineering 5th inter_ national congress. Istanbul: ITU; 2002 [in Turkish]. [28] Chopra AK. Earthquake response analysis of concrete dams. In: Jansen RB, editor. Advanced dam engineering for design, construction, and rehabilitation. New York: Van Nostrand Reinhold; 1988. [29] Reinhardt HW. Loading rate, temperature, and humidity effects. In: Shah SP, Carpinteri A, editors. Fracture mechanics test methods for concrete. London: Chapman & Hall; 1991. RILEM Report No. 89-FMT.