- Email: [email protected]
Earthquake response of concrete gravity dams including dam–foundation interface nonlinearities
Engineering Structures 30 (2008) 3065–3073
Contents lists available at ScienceDirect
Engineering Structures journal homepage: www.elsevier.com/locate/engstruct
Earthquake response of concrete gravity dams including dam–foundation interface nonlinearities Hamidreza Arabshahi ∗ , Vahid Lotfi 1 Department of Civil and Environmental Engineering, Amirkabir University of Technology, Tehran, Iran
article
info
Article history: Received 17 June 2007 Received in revised form 12 March 2008 Accepted 6 April 2008 Available online 27 May 2008 Keywords: Nonlinear dynamic analysis Interface element Sliding Uplifting Concrete gravity dam Plasticity
a b s t r a c t Various possibilities of natural isolation occurring along the dam–foundation interface during an earthquake are studied in order to evaluate the seismic response of an existing gravity dam. A plasticity -based formulation is used in the local stress space of interface elements to model sliding as well as partial opening along the dam base. The resistance to sliding is governed by two Mohr–Coulomb envelope surfaces representing the peak and residual strength of the joints. Effects of base deformation, sliding parameters and foundation flexibility are also investigated. Apart from sliding, other modes of failure of the dam are also identified. The results show that sliding generally reduces the maximum amount of tensile principal stresses in the dam body; however, the amount of reduction is generally not large enough to prevent cracking, especially at the upper parts of the dam. © 2008 Elsevier Ltd. All rights reserved.
1. Introduction Traditional approaches toward evaluating the sliding safety of gravity dams were based upon limit equilibrium approaches by assuming the dam to be a rigid body sliding along its base or lift joints under the action of equivalent static loads. Although it was accepted that providing a high safety factor against sliding be a requirement for the safety design of concrete gravity dams, recent results show that sliding can be used as a defensive design measure against earthquake, if it is properly controlled and an adequate safety factor against the critical stability limit be provided. This conclusion was based upon the results of analyses showing that a large reduction in tensile stresses [1] as well as base shear force [2] will occur when the dam is able to slide on its foundation. A more elaborated model used for the dam–water–foundation interaction by Chaves and Fenves [3], however, shows a different result indicating that the sliding displacement is generally not large enough to substantially reduce the stresses in the dam especially for one on flexible foundation rock.
∗ Corresponding address: Department of Civil and Environmental Engineering, Amirkabir University of Technology (Tehran Polytechnic), No. 424, Hafez Ave., Tehran, Iran. E-mail addresses: [email protected] (H. Arabshahi), [email protected] (V. Lotfi). 1 Address: Department of Civil and Environmental Engineering, Amirkabir University of Technology (Tehran Polytechnic), No. 424, Hafez Ave., Tehran, Iran. Tel.: +98 21 6414013; fax: +98 21 6414213. 0141-0296/$ – see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.engstruct.2008.04.018
These discrepancies in results lie mainly in the foundation model and the dam base condition used. Chaves and Fenves [3] showed that by taking account of the dam–foundation interaction in the frequency domain, a large amount of energy dissipation (as much as 80% of the total input energy) could occur through the foundation rock, thus lowering the response values including sliding as well as base shear force. A rigid dam base condition without any partial uplift was also considered in their analysis. These two factors cause sliding to have a less pronounced effect on response values. In contrast, the massless foundation model used by Leger and Katsouli [2] (which could not account for energy dissipation through radiation damping) as well as a flexible dam base condition resulted in an increase in the amount of sliding displacement, and thus cause sliding to play the primary role in lowering the response values. It should be mentioned that force vibration testing of concrete dams typically shows a small amount of damping and this amount should include all radiation damping as well [4]. So it seems that the model used by Chaves and Fenves [3] gives an upper bound solution for amount of radiation damping. Other investigators have also studied the sliding response of gravity dams and provided analytical procedures [5] as well as formulas [6] to estimate the seismic-induced sliding displacement at the base of concrete gravity dams. In all above mentioned studies as well as this paper a predefined surface (dam–foundation interface) was chosen for safety evaluation of the dam. This assumption may not be correct [7], due to the presence of substantial shearing at the tip of the joint which may favor the rotation of principal stresses in such
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
3066
obvious from Eq. (1), [Dl ] plays a key role on characterizing the joint’s behavior. Initially before the joint cracks, this matrix takes on penalty values to enforce the continuity condition, i.e.: KS 0 0 l e [D ] = [D ] = 0 KT 0 (3) 0 0 KN where KS , KT , KN are stiffness coefficients for two tangential and normal directions, respectively. In what follows a brief review on the four constitutive relationships used in the local stress space of the interface element will be presented: 2.1. Linear case (NCRIT = 0)
Fig. 1. The 16-node zero thickness interface element.
a way that the crack path dip into the surrounding materials [8]. It will be shown later that the current model used can also well predict this phenomenon. However, this assumption can be acceptable if it is assumed that the foundation rock is sufficiently stronger than the dam–foundation joint [9]. In the current study, several cases of natural isolation mechanisms are considered to obtain insight into the nonlinear seismic response of gravity dams with emphasis on the dynamic stress distribution in the dam body. Effects of various joints parameters are investigated and possible modes of failure are identified. In the upcoming sections, concepts of the constitutive relationships used for the interface element, modeling, loading and analysis parameters will be explained. Subsequently, the cases analyzed and the results of the analyses will be presented. Analyses are performed by a nonlinear finite element program developed by the authors. 2. Constitutive relationships for interface nonlinearities Formulation of constitutive relationships is mostly divided into elasto-plastic and fracture mechanics approaches where the former is often associated to the modeling of joints by means of contact/gap/interface elements and the latter is appropriate for evaluation of the formation and propagation of discontinuities. Based on the method used, the constitutive behavior of interfaces is governed by various parameters such as tensile strength, cohesion, friction angle (in elasto-plastic approach), fracture toughness and fracture energy (in linear and nonlinear fracture mechanics models). Framework of plasticity is used in the present study. In what follows, an outline of the joint element used (Fig. 1) as well as four different constitutive relationships formulated to model various separation mechanisms will be presented. The final form of the stiffness matrix of a zero thickness joint element can be written as [10]: Z (1) [K ] = [G]T [B]T [T ]T [Dl ][T ][B]dA [G] A
where [T ], [B] and [G] are matrices of rotation (between local and global coordinates of the joint element), strain interpolation and transformation of displacements to strains, respectively. Matrix [Dl ] relates the local stresses to the local strains of the interface element, i.e.: dτ1 du l dτ2 = [D ] dv (2) dσ dw where u, v, w are local strains in two tangential and normal directions and τ1 , τ2 , σ are their corresponding stresses. As it is
For the sake of comparison and showing the fundamental behavior of the dam, a linear model of the dam is analyzed in which relatively large values are set for the stiffness coefficients of the interface elements in order to prevent opening and sliding of the joints (joint strength is assumed to be infinity). Although interface elements are included in this model, the results should be the same as the model without these elements at the base. This occurs when the stiffness coefficients of the joint elements are chosen properly [11]. 2.2. Opening without sliding model (NCRIT = 1) In this case, joints at the base may open and close gradually during the analysis but cannot slide. This can be achieved by assigning constant penalty values to the two tangential stiffness coefficients of the joints which do not change throughout the analysis. However, as the tensile stress at an integration point reaches the tensile strength of the joint (ft ), it opens up and normal stress, tensile strength and normal stiffness coefficient of the joint will be set equal to zero. This model is appropriate for joints with shear keys of considerable height (Fig. 2(a)). 2.3. Opening with simple sliding model (NCRIT = 2) This is a simple model to account for sliding. The stiffness coefficients are set initially equal to penalty values but as tensile crack occurs at an integration point, perfect tensile and shear softening is assumed and the two tangential stiffness coefficients, tensile and shear stresses, tensile strength and the two shear strains of the joint will be set equal to zero (Fig. 2(b)). Later on, as the joint closes, penalty values would be activated for both normal and tangential directions. Care must be paid when the joint closes as only part of the shear strains would cause shear stresses, which can be obtained by multiplying the shear strain increments by the following factor: r1 =
(w)mi++11 (1w)im++11
(4)
where (w)im++11 is the normal strain at the end of iteration (i + 1) of i +1 load step (m + 1) and (1w)m +1 is the normal strain increment at that iteration. 2.4. Opening with sliding model (NCRIT = 3) The framework of plasticity is established in the local stress space of the interface element to model sliding as well as joint opening. Use of plasticity in the local stress space of interface element was also proposed by previous researchers [12–14]. The model is based on Mohr–Coulomb criterion and consists of two truncated cones representing an initial failure surface and a subsequent yield surface in order to account for the peak and
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
3067
Fig. 3. Various paths taken by the stress point in the local stress space of the interface element (NCRIT = 3); (a) the stress point has crossed the tension cut-off or the failure surface in tension region, (b) the stress point has crossed the failure surface and normal strain has become positive, (c) the stress point has crossed the failure surface in the compression region, (d) cracking and subsequently closing of the joint in one iteration.
at a positive normal strain (point A) and subsequently closes in that iteration. Only part of the strain vector can cause stress in this case which can be obtained by multiplying the strain increment vector by the factor obtained in Eq. (4).
Fig. 2. Various constitutive relationships used for the joint element; (a) joint opening model (NCRIT = 1), (b) opening with simple sliding model (NCRIT = 2), (c) opening with sliding model (NCRIT = 3).
residual strength of the joint (Fig. 2(c)). Also a tension cut-off is imposed on both surfaces to model opening of the joint. The tension cut-off is not considered as a yield surface but a criterion for restricting the tensile stress in the joint. Fig. 3 schematically represents some important features of the model used. A brief discussion is as follows: When the stress point reaches the failure surface in tension region or crosses the tension cutoff, it will be transferred to origin and local stresses as well as two tangential shear strains would become nil (Fig. 3(a)). Having crossed the failure surface, if the normal strain becomes positive during the analysis, the stress point will be transferred to the origin (Fig. 3(b)). If the stress point reaches the failure surface in compression region, it will be dropped onto the yield surface (Fig. 3(c)) and the two components of shear stresses will be modified by the following factor: n=
cres − µσ c − µσ
(5)
where c and cres are the amount of cohesion and its residual value, respectively, µ is the coefficient of friction and σ is the normal stress (compression is assumed to be negative). This factor is actually the ratio of shear stress on the yield surface to that of failure surface for a constant value of normal compressive stress. Other situations such as the one shown in Fig. 3(d) can also occur during the analysis. It shows a case in which the joint cracks
2.4.1. Elasto-plastic rigidity matrix formulation In this section the general procedure for obtaining the elastoplastic rigidity matrix of the interface elements will be explained. This matrix can be used in constructing the tangent stiffness of the interface element when the stress point is on the yield surface. The general equation of both failure and yield surfaces has the following form [12]: q F = τ12 + τ22 + µσ − c = 0 (6) which for the yield surface, the amount of cohesion (c) must be replaced by its residual value (cres ). The above equation represents a cone in the 3D local stress space of the interface element with (σ) as its axis. Following the standard procedure in plasticity and based on an associated flow rule, the elasto-plastic rigidity matrix can be formulated as:
[Dep ] = [De ] −
[De ] {a} {a}T [De ] (A + {a}T [De ] {a})
(7)
in which, A is equal to zero as no hardening is considered for the model and the vector {a} contains derivatives of the yield function with respect to local stresses, i.e.:
{a} =
∂F ∂F ∂F ∂τ1 ∂τ2 ∂σ
T
= {β1 β2 µ}T
(8)
where β1 and β2 are defined as:
τ1 β1 = q 2 τ1 + τ12
(9a)
τ2 β2 = q . 2 τ1 + τ12
(9b)
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
3068
By introducing Eqs. (3) and (8) into Eq. (7), the final form of the elasto-plastic rigidity matrix [Dep ] can be written as [14]: 1 KS 2 1 − β1 β2 KS KT − µβ1 KS KN KS 1 − H β1 H H 1 1 KT 2 (10) − β1 β2 KS KT − µβ2 KT KN KT 1 − β2 H H H 1 1 KN 2 − µβ1 KS KN − µβ2 KT KN KN 1 − µ H
H
H
in which H takes the following form: H = {a}T [De ]{a} =
β21 KS + β22 KT + µ2 KN .
Fig. 4. A case in which the loading function L is lower than zero, while part of the behavior is still elasto-plastic.
(11)
It is noted that shearing is accompanied by dilatancy in the above formulation due to off diagonal terms in matrix (10). Dilatancy has important consequences in confined problems, while it is anticipated to be less important for unconfined sliding mechanism such as the one considered in this study. 2.4.2. Calculating elasto-plastic stress The complete procedure for stress calculation is not presented here, only the procedure for calculating elasto-plastic stress will be reviewed. In the following formulation it is assumed that the stress point is on the yield surface. All the formulations are presented for the (n + 1)th step of sub-increment loading. Stress–strain relationship in plastic loading can be written as [15]:
{1σ˜ }im++11 = [De ] n+1 {1ε˜ }im++11 −[De ]{a}1λ (12) where {1ε˜ } is the strain sub-increment (obtained by dividing n+1
that portion of strain vector causing plastic loading into m parts) and {1σ˜ } is its corresponding stress sub-increment vector. 1λ is generally defined as:
1λ =
L H
=
{a}T [De ]{1ε}
(13)
H
in which L is the loading function. By introducing Eqs. (3) and (8) into Eq. (13) and make use of strain sub-increment vector {1ε˜ } instead of {1ε}in the above equation, 1λ can be written as:
1λ =
[β1
n+1
(1τ˜ 1e )im++11 +β2
n+1
(1τ˜ 2e )im++11 +µ H
n+1
(1σ˜ e )im++11 ]
(14)
where 1τ˜ 1e , 1τ˜ 2e , 1σ˜ e are the vector {1ε˜ } corresponding elastic stress sub-increments. 1λ is always greater than or equal to zero, so if it becomes negative during the analysis due to numerical error, it will be set equal to zero. Now considering Eq. (12), the elasto-plastic stress subincrements can be written as: n+1
(1τ˜ 1ep )im++11 = n+1 (1τ˜ 1e )im++11 −1λ KS β1
(15a)
n+1
(1τ˜ 2ep )im++11 = n+1 (1τ˜ 2e )im++11 −1λ KT β2
(15b)
n+1
(1σ˜ )
ep i+1 m+1
=
n+1
(1σ˜ )
e i +1 m+1
−1λ KN µ.
(15c)
Form the above equations one can observe that if 1λ = 0, the stress increment vector would become totally elastic. Obtaining the elasto-plastic stress increment, the stress vector can now be updated. An additional refinement is introduced at the end of each subincrement to scale the stress point to the yield surface. The scale factor used is: n=
cres − µσ
τ
(16)
where τ is the length of the resultant of the two shear stress components τ1 and τ2 . The two components of shear stress can now be modified as: updated
τ1 = nτ1
(17a)
updated
τ2 = nτ2 .
(17b)
Fig. 5. Finite element model of Pine Flat gravity dam.
2.4.3. A note on the effective stress and loading function As only three stress components are obtained for the interface element, general definition of the effective stress cannot be used here. However, a similar concept can be used in this case if the effective stress takes the following form: q σe = τ12 + τ22 + µσ. (18) By comparing the effective stress with cohesion or its residual value (cres ), the state of stress point can be determined at any instant. Also care must be paid in using the loading function to find the loading state of the stress point. As an example, consider Fig. 4 which shows a stress point on the yield surface subjected to an elastic stress increment vector. The loading function is negative (so as usually assumed unloading has occurred), but the stress point is outside the yield surface at the end of the loading step. This usually happens near the origin where the radius of the cross-section of the yield surface is small. Thus, apart from the loading function, the state of stress point at the end of loading step must also be checked. 3. Dynamic analysis 3.1. Modeling The procedure discussed is used to compute the earthquake response of Pine Flat gravity dam. It is a 561 m long concrete gravity dam near Fresno, California. The tallest of the thirty-seven monoliths is 122 m high, and each monolith is 12–15 m wide. The finite element model of the highest monolith of the dam with its foundation (for cases with flexible foundation) is shown in Fig. 5. The model consists of 40 and 48, 20-node isoparametric solid elements for the dam and foundation, respectively, and 4, 16node zero- thickness interface elements. A section of the dam with 10 m thickness is used for the analyses. A plane stress behavior is assumed for the dam concrete with a modulus of elasticity (E) of 22.75 GPa, a Poisson’s ratio (ν) of 0.2 and a unit weight (γc ) of
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
3069
Table 1 Vibration frequencies of Pine Flat gravity dam Mode
Vibration frequency (rad/s)
1 2 3 4 5 6 7 8 9 10
(Ef /Ed ) = ∞
(Ef /Ed ) = 1
(Ef /Ed ) = 0.25
15.48 33.24 53.71 57.11 82.59 106.23 113.95 122.44 129.18 136.95
12.66 27.04 34.32 45.70 70.39 84.90 97.28 105.30 109.70 121.25
8.91 19.36 21.68 36.97 62.59 76.99 92.41 95.00 96.57 114.61
Table 2 Cases analyzed
Model
NCRIT
Direction of the constraint applied at the base nodes
ft (MPa)
µ
c (MPa)
P0 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10
0 1 2 3 3 3 3 3 3 3 3
– – – z x, z – – – – – –
– 0.5 0.5 0 0 0.5 0 0 0.5 0.5 0.5
– – – 1 1 1 1 1 1.2 1 1
– – – 0 0 1 0.01 1 1 1 1
Ef Ed
Fig. 6. Envelope of maximum tensile principal stresses (MPa) for case P0.
∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ 1 0.25
24.8 kN/m3 . The foundation is assumed to be in a state of plane strain with a Poisson’s ratio of 0.33. The dam–reservoir interaction is accounted for by the Westergaard added mass approach and a finite massless foundation is used to include the dam–foundation interaction. The behavior of the dam, water and foundation system is assumed to be elastic and the nonlinear behavior is restricted to the dam–foundation interface zone. Table 1 lists the vibration frequencies of the dam–reservoir–foundation system for different ratios of the dam to foundation flexibility. 3.2. Loading and solution procedure The applied loads are divided into static and dynamic loads. Static loadings include gravity load which is applied first in one increment and hydrostatic and uplift pressure which are applied simultaneously in 10 increments. A triangular distribution is assumed for the uplift pressure which is zero at toe and increases up to 40% of the hydrostatic pressure at the heel of the dam. Static loadings are applied at negative range of time (before the earthquake loads). Also, 20 sec of the Taft record, S69E horizontal component scaled to a peak acceleration of 0.4g is used for dynamic loading which is applied at time zero. Analyses are performed in time domain using a time step of .001 sec. Also the Newmark method with a numerical damping equal to 20% [13] is used for the integration of equation of motion (β = 0.36, γ = 0.7). A Rayleigh damping matrix is used in the analyses and the coefficients are determined such that the damping ratios for the first and the fifth modes of vibration become 5% of the critical damping. 4. Discussion of results The cases analyzed referred to as P0 to P10 are given in Table 2 and the main results of the analyses are summarized in Table 3. Several cases are considered to obtain insight into nonlinear behavior of the dam. For all cases, the envelope of principal stresses
Fig. 7. Finite element model of Pine Flat gravity dam and its reservoir.
includes dynamic as well as static effects. Also it is assumed that the residual cohesion of the joint be one-tenth of its initial value. 4.1. Linear response of the dam (NCRIT = 0) For the first case, a linear model of the dam is analyzed. The envelope of maximum tensile stresses for this case is shown in Fig. 6. Large tensile stresses are developed near the base as a result of rigid foundation assumption and near the discontinuity in the slope of the downstream face of the dam due to stress concentration and narrow width of this area. The maximum amount of principal compression is well below the compressive strength of concrete and occurs at the downstream face of the crest. 4.1.1. Influence of water compressibility In order to estimate the effects of water compressibility on the response of the dam, the dam–reservoir interaction is taken into account with the use of 20-node fluid finite elements as shown in Fig. 7. The water is taken as compressible, inviscid fluid with the weight density of 10.0 kN/m3 and pressure wave velocity of 1440 m/s. The length of the reservoir discretized is 200 m, and Sommerfeld boundary condition is applied at the upstream boundary while, the reservoir bottom condition is assumed completely reflective. The envelope of maximum tensile stress is shown in Fig. 8. Compared with case P0, it is observed that the maximum tensile stress is approximately the same for both cases and has occurred at the heel of the dam. However, the maximum tensile stress near the crest is about 5.47 MPa for this case and increases to 6.32 MPa for case P0 (16% increase). Based on the above results, the added mass approach was considered to be used for the rest of the analyses. The justification of using added mass approach for the dam–water interaction has also been verified in previous studies [16].
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
3070 Table 3 Maximum responses of Pine Flat gravity dam Model
Crest horizontal displacement (mm)
Joint opening (mm)
Base sliding (mm)
Static and dynamic tensile principal stress (MPa)a
U/S
D/S
U/S
D/S
U/S
D/S
U/S
D/S
U/S
D/S
– – 24.6 110 116 156 177 156 127 193 199
– –
7.16 0.78 2.41 4.22 4.66 0.5 0.41 0.42 0.35 0.31 0.85
1.88 2.06 4.16 0.64 1.21 3.04 1.67 2.6 2.18 2.09 1.9
6.32 6.17 7.24 6.55 6.61 7.37 6.64 7.32 6.9 6.24 4.05
6.32 9.22 8.14 6.55 6.61 7.37 7.47 7.32 8.85 5.55 5.19
6.74 6.04 8.86 5.59 5.04 6.56 5.72 6.21 7.06 8.73 9.67
3.72 7.9 13.4 3.28 3.75 10.26 9.98 10.75 14.01 10.38 13.36
6.49 8.65 8.36 7.23 7.33 7.38 7.94 7.4 8.82 5.71 6.0
8.33 7.69 9.39 8.13 8.24 9.46 8.92 9.48 8.82 8.52 6.0
Base
P0 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10
107 105 102 86.5 86.7 92.2 93.7 92.1 96.2 106 140
– 44.3 40 – – 28.9 28 29 39.3 30.3 39.6
16 109 116 151 171 150 122 184 181
Static and dynamic compressive principal stress (MPa)a
Crest
Base
Crest
a Stresses are obtained at Gauss points, copied to the nearest node and an average value is subsequently obtained.
Fig. 8. Envelope of maximum tensile principal stresses (MPa)-water compressibility is included in the model.
Fig. 9. Envelope of maximum tensile principal stresses (MPa) for case P1.
4.2. Response of the dam to joint opening (NCRIT = 1) In case P1, tensile stresses at the base were released as they reach the tensile strength of the joints (only uplifting effect is considered). As can be seen from Fig. 9, a large amount of reduction in tensile stresses has occurred near the base of the dam but at the same time the maximum of tensile principal stresses has increased at the top by nearly 3 MPa (i.e. 46%) with respect to linear model. The time history of crest displacement is compared to that of model P0 in Fig. 10. It is observed that uplifting has reduced the response value over a wide range of frequencies, but it has also amplified the crest horizontal displacement in the part of strong shaking causing the maximum displacement to be nearly the same as that of model P0. Time history of joint opening at the heel is plotted in Fig. 11. It is noted that the crack first opened when the hydrostatic pressure was applied.
Fig. 10. Comparison of crest horizontal displacement for cases P0 and P1.
4.3. Response of the dam considering simple sliding model (NCRIT = 2) Fig. 11. Time history of joint opening at the heel for case P1.
A simple sliding model in along with joint opening is considered in case P2. The maximum value of tensile stresses has reduced with respect to model P1 (Fig. 12), but it is still above the linear case. Also it is noted that considerable tensile stresses are built up at the base of the dam. The directions of these principal stresses are nearly horizontal. This issue will be discussed further later on. The maximum amount of joint opening has reduced at the heel although the dam slides only 16 mm in this case (Table 3).
4.4. Response of the dam to sliding (NCRIT = 3) In cases P3 and P4, effects of sliding without joint opening on the response of the dam were investigated. In order to suppress joint opening at the dam base, the vertical degrees of freedom of the base nodes were constrained together, causing the base to be always in compression for these cases. The horizontal degrees of freedom were also constrained in case P4 to cause a rigid
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
3071
Fig. 14a. Comparison of crest horizontal displacement for cases P0 and P4.
Fig. 12. Envelope of maximum tensile principal stresses (MPa) for case P2.
Fig. 14b. Comparison of base sliding displacement for case P4 and that of reference 3.
Fig. 13. Envelope of maximum tensile principal stresses (MPa) for case P4.
dam base. The amount of cohesion is assumed to be nil for these two cases as the dam does not slide for c = 1 MPa. Envelope of maximum principal tensile stresses for case P4 is illustrated in Fig. 13. Maximum tensile stresses at the base have reduced by about 35% in comparison to linear model but they are still large enough to initiate cracking at the heel. In contrast, no such reduction in tensile stresses has occurred at the top part of the dam (Table 3). So it is observed that, in spite of sliding at the base, upper as well as lower parts of the dam are still vulnerable to substantial cracking. Time history of crest displacement for this case is compared to that of model P0 in Fig. 14a. It is noted that, in contrast to uplifting, sliding causes no amplification and outof-phaseness in the response value. Time history of base sliding displacement for this case is compared with that of reference 3 in Fig. 14b. The dam slides 116 mm which compares with 124 mm obtained by Chaves and Fenves [3]. The two curves are in good agreement; however, there are slight differences at sliding events which are mainly due to the assumption used in modeling (such as the effect of water compressibility). Base horizontal deformation in case P3 is also seen to have no substantial effect on the response values with respect to case P4 (Table 3). 4.5. Response of the dam to joint opening and sliding (NCRIT = 3) In case P5, the effects of sliding as well as uplifting on the response of the dam were investigated. The results (Fig. 15a) show that sliding has somehow reduced the enhancing effects of joint opening on the tensile principal stresses at the top of the dam (Table 3). However, the amount of sliding is still not large enough to lower these stresses with respect to the linear model (P0).
Fig. 15a. Envelope of maximum tensile principal stresses (MPa) for case P5.
Envelope of maximum joint opening at the base for case P5 is compared to that of model P1 in Fig. 15b. It is noted that sliding has reduced the maximum amount of joint opening at the heel while increasing its value at the toe. The envelope of maximum compressive stresses for case P5 is also shown in Fig. 15c. The peak principal compression has occurred at the toe. Although this amount of stress is well below the compressive strength of concrete, it may introduce some potential failure zone (i.e. crushing of concrete material at the toe) for stronger earthquakes. Time history of base sliding displacement is plotted in Fig. 16. The amounts of sliding at the heel and toe of the dam are not the same due to the dam base flexibility. It is noted that partial uplift has increased the amount of sliding by 30% in comparison to case P4. 4.6. Effects of joint parameters on response values In cases P6, P7 and P8 effects of joint parameters on response values are investigated (Table 2). The tensile strength of the joint is chosen to be nil in case P6 as it cannot be greater than (c/ tan φ) (Fig. 2(c)). The results for cases P6 and P7 show no substantial changes in response values with respect to case P5 for this range of variation in joint parameters. The reason is that most of the integration points crack during the initial phase of loading and thus these parameters will be set to their residual values for the rest of the analysis which are nearly the same as those of model P5. Fig. 17 shows an exaggerated deformation of the dam at time 8 sec of the
3072
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
Fig. 15b. Comparison of envelope of joint opening at base for cases P1 and P5. Fig. 17. Dam deformation at time 8 sec of analysis for case P7.
Fig. 15c. Envelope of maximum compressive principal stresses (MPa) for case P5. Fig. 18. Envelope of maximum tensile principal stresses (MPa) for case P10.
Fig. 16. Base sliding displacement of U/S and D/S nodes at base for case P5. Fig. 19. Base sliding displacement of U/S and D/S nodes at base for case P10.
analysis for case P7. Joint opening and sliding are clearly seen in this figure. Increasing the coefficient of friction in case P8 to 1.2 results in lower sliding displacement and higher amounts of tensile and compressive stresses (Table 3). 4.7. Effects of foundation rock flexibility The effect of foundation rock flexibility is investigated in cases P9 and P10. Table 3 shows that the maximum amount of sliding displacement has increased in these cases. Also the maximum amount of tensile stresses has reduced, but the amount of reduction is not so much even for the extreme case of Ef /Ed = 0.25 (Fig. 18). The time history of sliding displacement for case P10 (Fig. 19) shows a sudden downstream/upstream sliding event at nearly time 7 sec of the analysis. Studying the graph carefully, it is noted that this is due to the foundation flexibility as no such behavior is observed for the sliding displacement of the node at the toe. Fig. 20 shows an enlarged view of the tensile stress pattern at the dam–foundation interface for case P10. Like the previous
mentioned models with rigid foundation, it is observed that tensile stresses are built up at the base on a plane nearly perpendicular to that of the crack, suggesting that the crack path may tilt into the dam or foundation rock causing a curved surface. No such behavior is observed when uplifting is not allowed at the base. The amount of these stresses may seem to be low not causing cracking, but a finer mesh can catch higher stresses at this area. Similar results have been obtained previously especially for cracks at higher elevation [9]. It was concluded that these stresses are due to the high shear stresses transferred by the contact zone during maximum opening of the crack at the base. 5. Conclusions The nonlinear response of a concrete gravity dam was investigated with emphasis on dynamic stress distribution in the dam body. Three constitutive relationships were used in the local stress space of the interface elements to account for various isolation mechanisms (i.e. uplifting and sliding) occurring
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
Fig. 20. Envelope of maximum tensile principal stresses (MPa) at the dam–foundation interface for case P10.
along the dam–foundation interface. Simple practical models were also used to account for dam–water–foundation interactions. The following conclusions were drawn from results of the analyses. Results from linear analysis of the dam show high tensile stresses at the heel as well as the top of the dam. When uplifting is included in the model, tensile stresses at the base drastically reduce, however the maximum of principal tension considerably increases at the upper parts. Thus, uplifting can both reduce and amplify the response values. When sliding without any uplifting is forced to occur along the base, the tensile stresses at the heel reduce but no such reduction is observed for principal tension at the top of the dam. Upper as well as lower parts of the dam are still susceptible to substantial cracking in this case. In another case, sliding along with joint opening were allowed to occur. The results showed that sliding generally reduces the enhancing effects of uplifting on the tensile principal stresses at the top of the dam, but these stresses are still well above the tensile strength of concrete and thus can cause cracking at that zone. Also horizontal tensile stresses are developed at the base near the crack tip on a plane nearly perpendicular to that of the crack. Such stress distributions suggest that the crack path will not propagate along the predefined dam–foundation interface and may dip into the dam or foundation rock, respectively. No signs of overturning of the dam were observed for the cases considered as the maximum amount of joint opening at the base was limited to 44.3 mm. Joint tensile strength as well as cohesion do not seem to have substantial effects on the dam response values for the range considered in the analyses when uplifting as well as sliding were allowed to occur along the base but they can have substantial effects for cases with a rigid dam base assumption. The coefficient of friction, however, was found to be the most influential factor which can affect both tensile and compressive stresses, respectively.
3073
The massless foundation model used cannot well predict the behavior of the system as it generally overestimates the amount of sliding at the base and thus is not suitable if the sliding stability of the dam is of concern. In addition it underestimates the amount of reduction in tensile stresses in the dam body with respect to more detailed foundation models. Uplifting at the base of the dam cannot be prevented during an earthquake due to the low tensile strength of the dam–foundation interface zone. However, the sliding displacement can be controlled by providing shear keys at the base of the dam. Based on the above results, it is observed that the dam behaves much better if it is allowed to slide moderately on its foundation as sliding generally reduces the amount of tensile stresses in the dam body. However, it should be mentioned that such reduction in tensile stresses is generally not large enough to prevent cracking, especially at the upper parts of the dam. References [1] Hall JF, Dowling MJ, El-Aidi B. Defensive design of concrete gravity dams. Report no. EERL 91-02. Pasadena (CA): California Institute of Technology; 1991. 115 pp. [2] Leger P, Katsouli M. Seismic stability of concrete gravity dams. Earthquake Eng Struct Dynam 1989;18(6):889–902. [3] Chaves JW, Fenves GL. Earthquake analysis and response of concrete gravity dams including base sliding. Report no. UCB/EERC-93/07. Berkeley (CA): University of California at Berkeley; 1993. 184 pp. [4] Hall JF. Efficient non-linear seismic analysis of arch dams. Earthquake Eng Struct Dynam 1998;27(12):1425–44. [5] Chopra AK, Zhang L. Earthquake-induced base sliding of concrete gravity dams. J Struct Eng 1991;117(12):3698–719. [6] Danay A, Adeghe LN. Seismic induced slip of concrete gravity dams. J Struct Eng 1993;119(1):108–29. [7] Chandra Kishen JM. Recent developments in safety assessment of concrete gravity dams. Current Sci 2005;89(4):650–5. [8] ICOLD EUROPEAN CLUB. Sliding safety of existing gravity dams-final report. 2004. [9] El-Aidi B. Nonlinear earthquake response of concrete gravity dam systems. Report no. EERL 88-02. Pasadena (CA): California Institute of Technology, 1988, 170 pp. [10] Lotfi V. Application of discrete crack in nonlinear dynamic analysis of Shahid Rajaee arch dam. Internat J Eng 2001;14(1):9–23. [11] Lotfi V, Espandar R. An investigation of joints behavior in seismic response of arch dams. Electron J Struct Eng 2002;1(1):17–31. [12] Ghaboussi J, Wilson EL, Isenberg J. Finite element for rock joints and interfaces. J Soil Mech Foundation Division 1973. [13] Hohberg JM. A joint element for the nonlinear dynamic analysis of arch dams. Report no. 186. Switzerland: Institute of Structural Engineering, ETH Zurich. 1992. [14] Lau DT, Noruziaan B, Razaqpur AG. Modeling of contraction joints and shear sliding effects on earthquake response of arch dams. Earthquake Eng Struct Dynam 1998;27:1013–29. [15] Chen WF, Han DJ. Plasticity for structural engineers. NY: Springer-Verlag; 1988. [16] Asteris PG, Tzamtzis AD. Nonlinear seismic response analysis of realistic gravity dam–reservoir systems. Internat J Nonlinear Sci Numer Simul 2003; (4):329–38.
Contents lists available at ScienceDirect
Engineering Structures journal homepage: www.elsevier.com/locate/engstruct
Earthquake response of concrete gravity dams including dam–foundation interface nonlinearities Hamidreza Arabshahi ∗ , Vahid Lotfi 1 Department of Civil and Environmental Engineering, Amirkabir University of Technology, Tehran, Iran
article
info
Article history: Received 17 June 2007 Received in revised form 12 March 2008 Accepted 6 April 2008 Available online 27 May 2008 Keywords: Nonlinear dynamic analysis Interface element Sliding Uplifting Concrete gravity dam Plasticity
a b s t r a c t Various possibilities of natural isolation occurring along the dam–foundation interface during an earthquake are studied in order to evaluate the seismic response of an existing gravity dam. A plasticity -based formulation is used in the local stress space of interface elements to model sliding as well as partial opening along the dam base. The resistance to sliding is governed by two Mohr–Coulomb envelope surfaces representing the peak and residual strength of the joints. Effects of base deformation, sliding parameters and foundation flexibility are also investigated. Apart from sliding, other modes of failure of the dam are also identified. The results show that sliding generally reduces the maximum amount of tensile principal stresses in the dam body; however, the amount of reduction is generally not large enough to prevent cracking, especially at the upper parts of the dam. © 2008 Elsevier Ltd. All rights reserved.
1. Introduction Traditional approaches toward evaluating the sliding safety of gravity dams were based upon limit equilibrium approaches by assuming the dam to be a rigid body sliding along its base or lift joints under the action of equivalent static loads. Although it was accepted that providing a high safety factor against sliding be a requirement for the safety design of concrete gravity dams, recent results show that sliding can be used as a defensive design measure against earthquake, if it is properly controlled and an adequate safety factor against the critical stability limit be provided. This conclusion was based upon the results of analyses showing that a large reduction in tensile stresses [1] as well as base shear force [2] will occur when the dam is able to slide on its foundation. A more elaborated model used for the dam–water–foundation interaction by Chaves and Fenves [3], however, shows a different result indicating that the sliding displacement is generally not large enough to substantially reduce the stresses in the dam especially for one on flexible foundation rock.
∗ Corresponding address: Department of Civil and Environmental Engineering, Amirkabir University of Technology (Tehran Polytechnic), No. 424, Hafez Ave., Tehran, Iran. E-mail addresses: [email protected] (H. Arabshahi), [email protected] (V. Lotfi). 1 Address: Department of Civil and Environmental Engineering, Amirkabir University of Technology (Tehran Polytechnic), No. 424, Hafez Ave., Tehran, Iran. Tel.: +98 21 6414013; fax: +98 21 6414213. 0141-0296/$ – see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.engstruct.2008.04.018
These discrepancies in results lie mainly in the foundation model and the dam base condition used. Chaves and Fenves [3] showed that by taking account of the dam–foundation interaction in the frequency domain, a large amount of energy dissipation (as much as 80% of the total input energy) could occur through the foundation rock, thus lowering the response values including sliding as well as base shear force. A rigid dam base condition without any partial uplift was also considered in their analysis. These two factors cause sliding to have a less pronounced effect on response values. In contrast, the massless foundation model used by Leger and Katsouli [2] (which could not account for energy dissipation through radiation damping) as well as a flexible dam base condition resulted in an increase in the amount of sliding displacement, and thus cause sliding to play the primary role in lowering the response values. It should be mentioned that force vibration testing of concrete dams typically shows a small amount of damping and this amount should include all radiation damping as well [4]. So it seems that the model used by Chaves and Fenves [3] gives an upper bound solution for amount of radiation damping. Other investigators have also studied the sliding response of gravity dams and provided analytical procedures [5] as well as formulas [6] to estimate the seismic-induced sliding displacement at the base of concrete gravity dams. In all above mentioned studies as well as this paper a predefined surface (dam–foundation interface) was chosen for safety evaluation of the dam. This assumption may not be correct [7], due to the presence of substantial shearing at the tip of the joint which may favor the rotation of principal stresses in such
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
3066
obvious from Eq. (1), [Dl ] plays a key role on characterizing the joint’s behavior. Initially before the joint cracks, this matrix takes on penalty values to enforce the continuity condition, i.e.: KS 0 0 l e [D ] = [D ] = 0 KT 0 (3) 0 0 KN where KS , KT , KN are stiffness coefficients for two tangential and normal directions, respectively. In what follows a brief review on the four constitutive relationships used in the local stress space of the interface element will be presented: 2.1. Linear case (NCRIT = 0)
Fig. 1. The 16-node zero thickness interface element.
a way that the crack path dip into the surrounding materials [8]. It will be shown later that the current model used can also well predict this phenomenon. However, this assumption can be acceptable if it is assumed that the foundation rock is sufficiently stronger than the dam–foundation joint [9]. In the current study, several cases of natural isolation mechanisms are considered to obtain insight into the nonlinear seismic response of gravity dams with emphasis on the dynamic stress distribution in the dam body. Effects of various joints parameters are investigated and possible modes of failure are identified. In the upcoming sections, concepts of the constitutive relationships used for the interface element, modeling, loading and analysis parameters will be explained. Subsequently, the cases analyzed and the results of the analyses will be presented. Analyses are performed by a nonlinear finite element program developed by the authors. 2. Constitutive relationships for interface nonlinearities Formulation of constitutive relationships is mostly divided into elasto-plastic and fracture mechanics approaches where the former is often associated to the modeling of joints by means of contact/gap/interface elements and the latter is appropriate for evaluation of the formation and propagation of discontinuities. Based on the method used, the constitutive behavior of interfaces is governed by various parameters such as tensile strength, cohesion, friction angle (in elasto-plastic approach), fracture toughness and fracture energy (in linear and nonlinear fracture mechanics models). Framework of plasticity is used in the present study. In what follows, an outline of the joint element used (Fig. 1) as well as four different constitutive relationships formulated to model various separation mechanisms will be presented. The final form of the stiffness matrix of a zero thickness joint element can be written as [10]: Z (1) [K ] = [G]T [B]T [T ]T [Dl ][T ][B]dA [G] A
where [T ], [B] and [G] are matrices of rotation (between local and global coordinates of the joint element), strain interpolation and transformation of displacements to strains, respectively. Matrix [Dl ] relates the local stresses to the local strains of the interface element, i.e.: dτ1 du l dτ2 = [D ] dv (2) dσ dw where u, v, w are local strains in two tangential and normal directions and τ1 , τ2 , σ are their corresponding stresses. As it is
For the sake of comparison and showing the fundamental behavior of the dam, a linear model of the dam is analyzed in which relatively large values are set for the stiffness coefficients of the interface elements in order to prevent opening and sliding of the joints (joint strength is assumed to be infinity). Although interface elements are included in this model, the results should be the same as the model without these elements at the base. This occurs when the stiffness coefficients of the joint elements are chosen properly [11]. 2.2. Opening without sliding model (NCRIT = 1) In this case, joints at the base may open and close gradually during the analysis but cannot slide. This can be achieved by assigning constant penalty values to the two tangential stiffness coefficients of the joints which do not change throughout the analysis. However, as the tensile stress at an integration point reaches the tensile strength of the joint (ft ), it opens up and normal stress, tensile strength and normal stiffness coefficient of the joint will be set equal to zero. This model is appropriate for joints with shear keys of considerable height (Fig. 2(a)). 2.3. Opening with simple sliding model (NCRIT = 2) This is a simple model to account for sliding. The stiffness coefficients are set initially equal to penalty values but as tensile crack occurs at an integration point, perfect tensile and shear softening is assumed and the two tangential stiffness coefficients, tensile and shear stresses, tensile strength and the two shear strains of the joint will be set equal to zero (Fig. 2(b)). Later on, as the joint closes, penalty values would be activated for both normal and tangential directions. Care must be paid when the joint closes as only part of the shear strains would cause shear stresses, which can be obtained by multiplying the shear strain increments by the following factor: r1 =
(w)mi++11 (1w)im++11
(4)
where (w)im++11 is the normal strain at the end of iteration (i + 1) of i +1 load step (m + 1) and (1w)m +1 is the normal strain increment at that iteration. 2.4. Opening with sliding model (NCRIT = 3) The framework of plasticity is established in the local stress space of the interface element to model sliding as well as joint opening. Use of plasticity in the local stress space of interface element was also proposed by previous researchers [12–14]. The model is based on Mohr–Coulomb criterion and consists of two truncated cones representing an initial failure surface and a subsequent yield surface in order to account for the peak and
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
3067
Fig. 3. Various paths taken by the stress point in the local stress space of the interface element (NCRIT = 3); (a) the stress point has crossed the tension cut-off or the failure surface in tension region, (b) the stress point has crossed the failure surface and normal strain has become positive, (c) the stress point has crossed the failure surface in the compression region, (d) cracking and subsequently closing of the joint in one iteration.
at a positive normal strain (point A) and subsequently closes in that iteration. Only part of the strain vector can cause stress in this case which can be obtained by multiplying the strain increment vector by the factor obtained in Eq. (4).
Fig. 2. Various constitutive relationships used for the joint element; (a) joint opening model (NCRIT = 1), (b) opening with simple sliding model (NCRIT = 2), (c) opening with sliding model (NCRIT = 3).
residual strength of the joint (Fig. 2(c)). Also a tension cut-off is imposed on both surfaces to model opening of the joint. The tension cut-off is not considered as a yield surface but a criterion for restricting the tensile stress in the joint. Fig. 3 schematically represents some important features of the model used. A brief discussion is as follows: When the stress point reaches the failure surface in tension region or crosses the tension cutoff, it will be transferred to origin and local stresses as well as two tangential shear strains would become nil (Fig. 3(a)). Having crossed the failure surface, if the normal strain becomes positive during the analysis, the stress point will be transferred to the origin (Fig. 3(b)). If the stress point reaches the failure surface in compression region, it will be dropped onto the yield surface (Fig. 3(c)) and the two components of shear stresses will be modified by the following factor: n=
cres − µσ c − µσ
(5)
where c and cres are the amount of cohesion and its residual value, respectively, µ is the coefficient of friction and σ is the normal stress (compression is assumed to be negative). This factor is actually the ratio of shear stress on the yield surface to that of failure surface for a constant value of normal compressive stress. Other situations such as the one shown in Fig. 3(d) can also occur during the analysis. It shows a case in which the joint cracks
2.4.1. Elasto-plastic rigidity matrix formulation In this section the general procedure for obtaining the elastoplastic rigidity matrix of the interface elements will be explained. This matrix can be used in constructing the tangent stiffness of the interface element when the stress point is on the yield surface. The general equation of both failure and yield surfaces has the following form [12]: q F = τ12 + τ22 + µσ − c = 0 (6) which for the yield surface, the amount of cohesion (c) must be replaced by its residual value (cres ). The above equation represents a cone in the 3D local stress space of the interface element with (σ) as its axis. Following the standard procedure in plasticity and based on an associated flow rule, the elasto-plastic rigidity matrix can be formulated as:
[Dep ] = [De ] −
[De ] {a} {a}T [De ] (A + {a}T [De ] {a})
(7)
in which, A is equal to zero as no hardening is considered for the model and the vector {a} contains derivatives of the yield function with respect to local stresses, i.e.:
{a} =
∂F ∂F ∂F ∂τ1 ∂τ2 ∂σ
T
= {β1 β2 µ}T
(8)
where β1 and β2 are defined as:
τ1 β1 = q 2 τ1 + τ12
(9a)
τ2 β2 = q . 2 τ1 + τ12
(9b)
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
3068
By introducing Eqs. (3) and (8) into Eq. (7), the final form of the elasto-plastic rigidity matrix [Dep ] can be written as [14]: 1 KS 2 1 − β1 β2 KS KT − µβ1 KS KN KS 1 − H β1 H H 1 1 KT 2 (10) − β1 β2 KS KT − µβ2 KT KN KT 1 − β2 H H H 1 1 KN 2 − µβ1 KS KN − µβ2 KT KN KN 1 − µ H
H
H
in which H takes the following form: H = {a}T [De ]{a} =
β21 KS + β22 KT + µ2 KN .
Fig. 4. A case in which the loading function L is lower than zero, while part of the behavior is still elasto-plastic.
(11)
It is noted that shearing is accompanied by dilatancy in the above formulation due to off diagonal terms in matrix (10). Dilatancy has important consequences in confined problems, while it is anticipated to be less important for unconfined sliding mechanism such as the one considered in this study. 2.4.2. Calculating elasto-plastic stress The complete procedure for stress calculation is not presented here, only the procedure for calculating elasto-plastic stress will be reviewed. In the following formulation it is assumed that the stress point is on the yield surface. All the formulations are presented for the (n + 1)th step of sub-increment loading. Stress–strain relationship in plastic loading can be written as [15]:
{1σ˜ }im++11 = [De ] n+1 {1ε˜ }im++11 −[De ]{a}1λ (12) where {1ε˜ } is the strain sub-increment (obtained by dividing n+1
that portion of strain vector causing plastic loading into m parts) and {1σ˜ } is its corresponding stress sub-increment vector. 1λ is generally defined as:
1λ =
L H
=
{a}T [De ]{1ε}
(13)
H
in which L is the loading function. By introducing Eqs. (3) and (8) into Eq. (13) and make use of strain sub-increment vector {1ε˜ } instead of {1ε}in the above equation, 1λ can be written as:
1λ =
[β1
n+1
(1τ˜ 1e )im++11 +β2
n+1
(1τ˜ 2e )im++11 +µ H
n+1
(1σ˜ e )im++11 ]
(14)
where 1τ˜ 1e , 1τ˜ 2e , 1σ˜ e are the vector {1ε˜ } corresponding elastic stress sub-increments. 1λ is always greater than or equal to zero, so if it becomes negative during the analysis due to numerical error, it will be set equal to zero. Now considering Eq. (12), the elasto-plastic stress subincrements can be written as: n+1
(1τ˜ 1ep )im++11 = n+1 (1τ˜ 1e )im++11 −1λ KS β1
(15a)
n+1
(1τ˜ 2ep )im++11 = n+1 (1τ˜ 2e )im++11 −1λ KT β2
(15b)
n+1
(1σ˜ )
ep i+1 m+1
=
n+1
(1σ˜ )
e i +1 m+1
−1λ KN µ.
(15c)
Form the above equations one can observe that if 1λ = 0, the stress increment vector would become totally elastic. Obtaining the elasto-plastic stress increment, the stress vector can now be updated. An additional refinement is introduced at the end of each subincrement to scale the stress point to the yield surface. The scale factor used is: n=
cres − µσ
τ
(16)
where τ is the length of the resultant of the two shear stress components τ1 and τ2 . The two components of shear stress can now be modified as: updated
τ1 = nτ1
(17a)
updated
τ2 = nτ2 .
(17b)
Fig. 5. Finite element model of Pine Flat gravity dam.
2.4.3. A note on the effective stress and loading function As only three stress components are obtained for the interface element, general definition of the effective stress cannot be used here. However, a similar concept can be used in this case if the effective stress takes the following form: q σe = τ12 + τ22 + µσ. (18) By comparing the effective stress with cohesion or its residual value (cres ), the state of stress point can be determined at any instant. Also care must be paid in using the loading function to find the loading state of the stress point. As an example, consider Fig. 4 which shows a stress point on the yield surface subjected to an elastic stress increment vector. The loading function is negative (so as usually assumed unloading has occurred), but the stress point is outside the yield surface at the end of the loading step. This usually happens near the origin where the radius of the cross-section of the yield surface is small. Thus, apart from the loading function, the state of stress point at the end of loading step must also be checked. 3. Dynamic analysis 3.1. Modeling The procedure discussed is used to compute the earthquake response of Pine Flat gravity dam. It is a 561 m long concrete gravity dam near Fresno, California. The tallest of the thirty-seven monoliths is 122 m high, and each monolith is 12–15 m wide. The finite element model of the highest monolith of the dam with its foundation (for cases with flexible foundation) is shown in Fig. 5. The model consists of 40 and 48, 20-node isoparametric solid elements for the dam and foundation, respectively, and 4, 16node zero- thickness interface elements. A section of the dam with 10 m thickness is used for the analyses. A plane stress behavior is assumed for the dam concrete with a modulus of elasticity (E) of 22.75 GPa, a Poisson’s ratio (ν) of 0.2 and a unit weight (γc ) of
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
3069
Table 1 Vibration frequencies of Pine Flat gravity dam Mode
Vibration frequency (rad/s)
1 2 3 4 5 6 7 8 9 10
(Ef /Ed ) = ∞
(Ef /Ed ) = 1
(Ef /Ed ) = 0.25
15.48 33.24 53.71 57.11 82.59 106.23 113.95 122.44 129.18 136.95
12.66 27.04 34.32 45.70 70.39 84.90 97.28 105.30 109.70 121.25
8.91 19.36 21.68 36.97 62.59 76.99 92.41 95.00 96.57 114.61
Table 2 Cases analyzed
Model
NCRIT
Direction of the constraint applied at the base nodes
ft (MPa)
µ
c (MPa)
P0 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10
0 1 2 3 3 3 3 3 3 3 3
– – – z x, z – – – – – –
– 0.5 0.5 0 0 0.5 0 0 0.5 0.5 0.5
– – – 1 1 1 1 1 1.2 1 1
– – – 0 0 1 0.01 1 1 1 1
Ef Ed
Fig. 6. Envelope of maximum tensile principal stresses (MPa) for case P0.
∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ 1 0.25
24.8 kN/m3 . The foundation is assumed to be in a state of plane strain with a Poisson’s ratio of 0.33. The dam–reservoir interaction is accounted for by the Westergaard added mass approach and a finite massless foundation is used to include the dam–foundation interaction. The behavior of the dam, water and foundation system is assumed to be elastic and the nonlinear behavior is restricted to the dam–foundation interface zone. Table 1 lists the vibration frequencies of the dam–reservoir–foundation system for different ratios of the dam to foundation flexibility. 3.2. Loading and solution procedure The applied loads are divided into static and dynamic loads. Static loadings include gravity load which is applied first in one increment and hydrostatic and uplift pressure which are applied simultaneously in 10 increments. A triangular distribution is assumed for the uplift pressure which is zero at toe and increases up to 40% of the hydrostatic pressure at the heel of the dam. Static loadings are applied at negative range of time (before the earthquake loads). Also, 20 sec of the Taft record, S69E horizontal component scaled to a peak acceleration of 0.4g is used for dynamic loading which is applied at time zero. Analyses are performed in time domain using a time step of .001 sec. Also the Newmark method with a numerical damping equal to 20% [13] is used for the integration of equation of motion (β = 0.36, γ = 0.7). A Rayleigh damping matrix is used in the analyses and the coefficients are determined such that the damping ratios for the first and the fifth modes of vibration become 5% of the critical damping. 4. Discussion of results The cases analyzed referred to as P0 to P10 are given in Table 2 and the main results of the analyses are summarized in Table 3. Several cases are considered to obtain insight into nonlinear behavior of the dam. For all cases, the envelope of principal stresses
Fig. 7. Finite element model of Pine Flat gravity dam and its reservoir.
includes dynamic as well as static effects. Also it is assumed that the residual cohesion of the joint be one-tenth of its initial value. 4.1. Linear response of the dam (NCRIT = 0) For the first case, a linear model of the dam is analyzed. The envelope of maximum tensile stresses for this case is shown in Fig. 6. Large tensile stresses are developed near the base as a result of rigid foundation assumption and near the discontinuity in the slope of the downstream face of the dam due to stress concentration and narrow width of this area. The maximum amount of principal compression is well below the compressive strength of concrete and occurs at the downstream face of the crest. 4.1.1. Influence of water compressibility In order to estimate the effects of water compressibility on the response of the dam, the dam–reservoir interaction is taken into account with the use of 20-node fluid finite elements as shown in Fig. 7. The water is taken as compressible, inviscid fluid with the weight density of 10.0 kN/m3 and pressure wave velocity of 1440 m/s. The length of the reservoir discretized is 200 m, and Sommerfeld boundary condition is applied at the upstream boundary while, the reservoir bottom condition is assumed completely reflective. The envelope of maximum tensile stress is shown in Fig. 8. Compared with case P0, it is observed that the maximum tensile stress is approximately the same for both cases and has occurred at the heel of the dam. However, the maximum tensile stress near the crest is about 5.47 MPa for this case and increases to 6.32 MPa for case P0 (16% increase). Based on the above results, the added mass approach was considered to be used for the rest of the analyses. The justification of using added mass approach for the dam–water interaction has also been verified in previous studies [16].
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
3070 Table 3 Maximum responses of Pine Flat gravity dam Model
Crest horizontal displacement (mm)
Joint opening (mm)
Base sliding (mm)
Static and dynamic tensile principal stress (MPa)a
U/S
D/S
U/S
D/S
U/S
D/S
U/S
D/S
U/S
D/S
– – 24.6 110 116 156 177 156 127 193 199
– –
7.16 0.78 2.41 4.22 4.66 0.5 0.41 0.42 0.35 0.31 0.85
1.88 2.06 4.16 0.64 1.21 3.04 1.67 2.6 2.18 2.09 1.9
6.32 6.17 7.24 6.55 6.61 7.37 6.64 7.32 6.9 6.24 4.05
6.32 9.22 8.14 6.55 6.61 7.37 7.47 7.32 8.85 5.55 5.19
6.74 6.04 8.86 5.59 5.04 6.56 5.72 6.21 7.06 8.73 9.67
3.72 7.9 13.4 3.28 3.75 10.26 9.98 10.75 14.01 10.38 13.36
6.49 8.65 8.36 7.23 7.33 7.38 7.94 7.4 8.82 5.71 6.0
8.33 7.69 9.39 8.13 8.24 9.46 8.92 9.48 8.82 8.52 6.0
Base
P0 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10
107 105 102 86.5 86.7 92.2 93.7 92.1 96.2 106 140
– 44.3 40 – – 28.9 28 29 39.3 30.3 39.6
16 109 116 151 171 150 122 184 181
Static and dynamic compressive principal stress (MPa)a
Crest
Base
Crest
a Stresses are obtained at Gauss points, copied to the nearest node and an average value is subsequently obtained.
Fig. 8. Envelope of maximum tensile principal stresses (MPa)-water compressibility is included in the model.
Fig. 9. Envelope of maximum tensile principal stresses (MPa) for case P1.
4.2. Response of the dam to joint opening (NCRIT = 1) In case P1, tensile stresses at the base were released as they reach the tensile strength of the joints (only uplifting effect is considered). As can be seen from Fig. 9, a large amount of reduction in tensile stresses has occurred near the base of the dam but at the same time the maximum of tensile principal stresses has increased at the top by nearly 3 MPa (i.e. 46%) with respect to linear model. The time history of crest displacement is compared to that of model P0 in Fig. 10. It is observed that uplifting has reduced the response value over a wide range of frequencies, but it has also amplified the crest horizontal displacement in the part of strong shaking causing the maximum displacement to be nearly the same as that of model P0. Time history of joint opening at the heel is plotted in Fig. 11. It is noted that the crack first opened when the hydrostatic pressure was applied.
Fig. 10. Comparison of crest horizontal displacement for cases P0 and P1.
4.3. Response of the dam considering simple sliding model (NCRIT = 2) Fig. 11. Time history of joint opening at the heel for case P1.
A simple sliding model in along with joint opening is considered in case P2. The maximum value of tensile stresses has reduced with respect to model P1 (Fig. 12), but it is still above the linear case. Also it is noted that considerable tensile stresses are built up at the base of the dam. The directions of these principal stresses are nearly horizontal. This issue will be discussed further later on. The maximum amount of joint opening has reduced at the heel although the dam slides only 16 mm in this case (Table 3).
4.4. Response of the dam to sliding (NCRIT = 3) In cases P3 and P4, effects of sliding without joint opening on the response of the dam were investigated. In order to suppress joint opening at the dam base, the vertical degrees of freedom of the base nodes were constrained together, causing the base to be always in compression for these cases. The horizontal degrees of freedom were also constrained in case P4 to cause a rigid
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
3071
Fig. 14a. Comparison of crest horizontal displacement for cases P0 and P4.
Fig. 12. Envelope of maximum tensile principal stresses (MPa) for case P2.
Fig. 14b. Comparison of base sliding displacement for case P4 and that of reference 3.
Fig. 13. Envelope of maximum tensile principal stresses (MPa) for case P4.
dam base. The amount of cohesion is assumed to be nil for these two cases as the dam does not slide for c = 1 MPa. Envelope of maximum principal tensile stresses for case P4 is illustrated in Fig. 13. Maximum tensile stresses at the base have reduced by about 35% in comparison to linear model but they are still large enough to initiate cracking at the heel. In contrast, no such reduction in tensile stresses has occurred at the top part of the dam (Table 3). So it is observed that, in spite of sliding at the base, upper as well as lower parts of the dam are still vulnerable to substantial cracking. Time history of crest displacement for this case is compared to that of model P0 in Fig. 14a. It is noted that, in contrast to uplifting, sliding causes no amplification and outof-phaseness in the response value. Time history of base sliding displacement for this case is compared with that of reference 3 in Fig. 14b. The dam slides 116 mm which compares with 124 mm obtained by Chaves and Fenves [3]. The two curves are in good agreement; however, there are slight differences at sliding events which are mainly due to the assumption used in modeling (such as the effect of water compressibility). Base horizontal deformation in case P3 is also seen to have no substantial effect on the response values with respect to case P4 (Table 3). 4.5. Response of the dam to joint opening and sliding (NCRIT = 3) In case P5, the effects of sliding as well as uplifting on the response of the dam were investigated. The results (Fig. 15a) show that sliding has somehow reduced the enhancing effects of joint opening on the tensile principal stresses at the top of the dam (Table 3). However, the amount of sliding is still not large enough to lower these stresses with respect to the linear model (P0).
Fig. 15a. Envelope of maximum tensile principal stresses (MPa) for case P5.
Envelope of maximum joint opening at the base for case P5 is compared to that of model P1 in Fig. 15b. It is noted that sliding has reduced the maximum amount of joint opening at the heel while increasing its value at the toe. The envelope of maximum compressive stresses for case P5 is also shown in Fig. 15c. The peak principal compression has occurred at the toe. Although this amount of stress is well below the compressive strength of concrete, it may introduce some potential failure zone (i.e. crushing of concrete material at the toe) for stronger earthquakes. Time history of base sliding displacement is plotted in Fig. 16. The amounts of sliding at the heel and toe of the dam are not the same due to the dam base flexibility. It is noted that partial uplift has increased the amount of sliding by 30% in comparison to case P4. 4.6. Effects of joint parameters on response values In cases P6, P7 and P8 effects of joint parameters on response values are investigated (Table 2). The tensile strength of the joint is chosen to be nil in case P6 as it cannot be greater than (c/ tan φ) (Fig. 2(c)). The results for cases P6 and P7 show no substantial changes in response values with respect to case P5 for this range of variation in joint parameters. The reason is that most of the integration points crack during the initial phase of loading and thus these parameters will be set to their residual values for the rest of the analysis which are nearly the same as those of model P5. Fig. 17 shows an exaggerated deformation of the dam at time 8 sec of the
3072
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
Fig. 15b. Comparison of envelope of joint opening at base for cases P1 and P5. Fig. 17. Dam deformation at time 8 sec of analysis for case P7.
Fig. 15c. Envelope of maximum compressive principal stresses (MPa) for case P5. Fig. 18. Envelope of maximum tensile principal stresses (MPa) for case P10.
Fig. 16. Base sliding displacement of U/S and D/S nodes at base for case P5. Fig. 19. Base sliding displacement of U/S and D/S nodes at base for case P10.
analysis for case P7. Joint opening and sliding are clearly seen in this figure. Increasing the coefficient of friction in case P8 to 1.2 results in lower sliding displacement and higher amounts of tensile and compressive stresses (Table 3). 4.7. Effects of foundation rock flexibility The effect of foundation rock flexibility is investigated in cases P9 and P10. Table 3 shows that the maximum amount of sliding displacement has increased in these cases. Also the maximum amount of tensile stresses has reduced, but the amount of reduction is not so much even for the extreme case of Ef /Ed = 0.25 (Fig. 18). The time history of sliding displacement for case P10 (Fig. 19) shows a sudden downstream/upstream sliding event at nearly time 7 sec of the analysis. Studying the graph carefully, it is noted that this is due to the foundation flexibility as no such behavior is observed for the sliding displacement of the node at the toe. Fig. 20 shows an enlarged view of the tensile stress pattern at the dam–foundation interface for case P10. Like the previous
mentioned models with rigid foundation, it is observed that tensile stresses are built up at the base on a plane nearly perpendicular to that of the crack, suggesting that the crack path may tilt into the dam or foundation rock causing a curved surface. No such behavior is observed when uplifting is not allowed at the base. The amount of these stresses may seem to be low not causing cracking, but a finer mesh can catch higher stresses at this area. Similar results have been obtained previously especially for cracks at higher elevation [9]. It was concluded that these stresses are due to the high shear stresses transferred by the contact zone during maximum opening of the crack at the base. 5. Conclusions The nonlinear response of a concrete gravity dam was investigated with emphasis on dynamic stress distribution in the dam body. Three constitutive relationships were used in the local stress space of the interface elements to account for various isolation mechanisms (i.e. uplifting and sliding) occurring
H. Arabshahi, V. Lotfi / Engineering Structures 30 (2008) 3065–3073
Fig. 20. Envelope of maximum tensile principal stresses (MPa) at the dam–foundation interface for case P10.
along the dam–foundation interface. Simple practical models were also used to account for dam–water–foundation interactions. The following conclusions were drawn from results of the analyses. Results from linear analysis of the dam show high tensile stresses at the heel as well as the top of the dam. When uplifting is included in the model, tensile stresses at the base drastically reduce, however the maximum of principal tension considerably increases at the upper parts. Thus, uplifting can both reduce and amplify the response values. When sliding without any uplifting is forced to occur along the base, the tensile stresses at the heel reduce but no such reduction is observed for principal tension at the top of the dam. Upper as well as lower parts of the dam are still susceptible to substantial cracking in this case. In another case, sliding along with joint opening were allowed to occur. The results showed that sliding generally reduces the enhancing effects of uplifting on the tensile principal stresses at the top of the dam, but these stresses are still well above the tensile strength of concrete and thus can cause cracking at that zone. Also horizontal tensile stresses are developed at the base near the crack tip on a plane nearly perpendicular to that of the crack. Such stress distributions suggest that the crack path will not propagate along the predefined dam–foundation interface and may dip into the dam or foundation rock, respectively. No signs of overturning of the dam were observed for the cases considered as the maximum amount of joint opening at the base was limited to 44.3 mm. Joint tensile strength as well as cohesion do not seem to have substantial effects on the dam response values for the range considered in the analyses when uplifting as well as sliding were allowed to occur along the base but they can have substantial effects for cases with a rigid dam base assumption. The coefficient of friction, however, was found to be the most influential factor which can affect both tensile and compressive stresses, respectively.
3073
The massless foundation model used cannot well predict the behavior of the system as it generally overestimates the amount of sliding at the base and thus is not suitable if the sliding stability of the dam is of concern. In addition it underestimates the amount of reduction in tensile stresses in the dam body with respect to more detailed foundation models. Uplifting at the base of the dam cannot be prevented during an earthquake due to the low tensile strength of the dam–foundation interface zone. However, the sliding displacement can be controlled by providing shear keys at the base of the dam. Based on the above results, it is observed that the dam behaves much better if it is allowed to slide moderately on its foundation as sliding generally reduces the amount of tensile stresses in the dam body. However, it should be mentioned that such reduction in tensile stresses is generally not large enough to prevent cracking, especially at the upper parts of the dam. References [1] Hall JF, Dowling MJ, El-Aidi B. Defensive design of concrete gravity dams. Report no. EERL 91-02. Pasadena (CA): California Institute of Technology; 1991. 115 pp. [2] Leger P, Katsouli M. Seismic stability of concrete gravity dams. Earthquake Eng Struct Dynam 1989;18(6):889–902. [3] Chaves JW, Fenves GL. Earthquake analysis and response of concrete gravity dams including base sliding. Report no. UCB/EERC-93/07. Berkeley (CA): University of California at Berkeley; 1993. 184 pp. [4] Hall JF. Efficient non-linear seismic analysis of arch dams. Earthquake Eng Struct Dynam 1998;27(12):1425–44. [5] Chopra AK, Zhang L. Earthquake-induced base sliding of concrete gravity dams. J Struct Eng 1991;117(12):3698–719. [6] Danay A, Adeghe LN. Seismic induced slip of concrete gravity dams. J Struct Eng 1993;119(1):108–29. [7] Chandra Kishen JM. Recent developments in safety assessment of concrete gravity dams. Current Sci 2005;89(4):650–5. [8] ICOLD EUROPEAN CLUB. Sliding safety of existing gravity dams-final report. 2004. [9] El-Aidi B. Nonlinear earthquake response of concrete gravity dam systems. Report no. EERL 88-02. Pasadena (CA): California Institute of Technology, 1988, 170 pp. [10] Lotfi V. Application of discrete crack in nonlinear dynamic analysis of Shahid Rajaee arch dam. Internat J Eng 2001;14(1):9–23. [11] Lotfi V, Espandar R. An investigation of joints behavior in seismic response of arch dams. Electron J Struct Eng 2002;1(1):17–31. [12] Ghaboussi J, Wilson EL, Isenberg J. Finite element for rock joints and interfaces. J Soil Mech Foundation Division 1973. [13] Hohberg JM. A joint element for the nonlinear dynamic analysis of arch dams. Report no. 186. Switzerland: Institute of Structural Engineering, ETH Zurich. 1992. [14] Lau DT, Noruziaan B, Razaqpur AG. Modeling of contraction joints and shear sliding effects on earthquake response of arch dams. Earthquake Eng Struct Dynam 1998;27:1013–29. [15] Chen WF, Han DJ. Plasticity for structural engineers. NY: Springer-Verlag; 1988. [16] Asteris PG, Tzamtzis AD. Nonlinear seismic response analysis of realistic gravity dam–reservoir systems. Internat J Nonlinear Sci Numer Simul 2003; (4):329–38.