- Email: [email protected]
Comparison of non-orthogonal smeared crack and plasticity models for dynamic analysis of concrete arch dams
Computers and Structures 81 (2003) 1461–1474 www.elsevier.com/locate/compstruc
Comparison of non-orthogonal smeared crack and plasticity models for dynamic analysis of concrete arch dams Radin Espandar a
a,*
, Vahid Lotfi
b
Department of Earth Sciences, Shahid Beheshti University, P.O. Box 13145-165, Tehran, Iran b Department of Civil and Environmental Engineering, Amirkabir University, Tehran, Iran Received 19 April 2002; accepted 20 January 2003
Abstract A special three-dimensional finite element program is developed to study the seismic response of concrete arch dams. Two types of continuum mechanics non-linear techniques are incorporated in the program, i.e. non-orthogonal smeared crack and elasto-plastic models. The formulation is presented initially and accuracy of the models and their correct implementations are verified by analyzing a notched beam and comparing its results with available experimental data. Following the satisfactory results of the beam analysis, the program is applied to a practical problem, viz. the non-linear behavior of the 130 m high Shahid Rajaee arch dam in Iran, subjected to the Friuli–Tolmezzo earthquake. Two cases of non-orthogonal smeared crack model are implemented for the dam, namely, ordinary and fictitious approaches. Meanwhile, an elastic–perfectly plastic model for mass concrete is also applied. The results are compared with those obtained from a linear elastic analysis under the same excitation. Based on the results, it is concluded that the nonorthogonal smeared crack approach can redistribute the state of stresses and produces a more realistic profile of stresses in the dam. Further, the elasto-plastic model could reduce significantly the high amount of overstressing noticed in the linear analysis. However, the elasto-plastic model is not perceptibly activated in the compressive state of stresses. It is also observed that a drift in the crest displacements is the prominent consequence of cracking or plasticity of the dam. Ó 2003 Elsevier Science Ltd. All rights reserved. Keywords: Non-linear dynamic analysis; Concrete arch dam; Smeared crack approach; Plasticity; Concrete modeling; Seismic load
1. Introduction There are several approaches to model the complicated stress–strain behavior of concrete under multiaxial stress states. In this regard, the theories based on plasticity and fracture mechanics are commonly applied in most of the engineering analysis. However, only few research works have been performed to study the dynamic behavior of concrete arch dams with concrete modeling. Kuo [1] suggested interface smeared crack
* Corresponding author. Tel.: +98-911-226-0861; fax: +9821-831-0043. E-mail addresses: [email protected] (R. Espandar), vahlotfi@aut.ac.ir (V. Lotfi).
approach to model the contraction joints and applied this technique for dynamic analysis of arch dams. Dungar [2] used a bounding surface model for static analysis of an arch dam. Cervera et al. [3] utilized a continuum damage model to analyze Talvacchia arch dam in Italy. In another research work, Noruziaan [4] applied bounding surface and orthogonal smeared crack models for compressive and tensile regions of stresses, respectively. However, it was reported that numerical problems occurred in some cases. The most recent related research work belongs to Hall [5], which a simple smeared crack model is utilized for modeling the contraction and construction joints. However, the diagonal cracks are neglected in that work. In the continuum models, such as elasto-plastic and smeared crack models, the dam body is assumed as a
0045-7949/03/$ - see front matter Ó 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0045-7949(03)00083-X
1462
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
monolith and the non-linear behavior is limited to the mass concrete of the dam body. The behavior of dam is investigated by studying the displacement histories and envelopes of principal stresses. The mentioned parameters are studied qualitatively, due to lack of any standard. In this manner, displacements of the dam are reasonable when they are not so large to endanger the general stability of the structure. Also, overstressing in tensile and compressive stresses, i.e. stresses higher than strengths, is used as a criterion for evaluating effectiveness of models in correct redistribution of stresses throughout the dam body. Furthermore, reasonable conclusions can be deduced by comparison among different models.
2. Description of analysis procedure The earthquake response of concrete arch dams is affected by dam–water and dam–foundation interactions as well as material non-linearities and joint movements. Although the interaction effects can be rigorously included in a frequency-domain analysis for a linear model, a time domain analysis is necessary for nonlinear models. In compliance with the main purpose and scope of the present work, i.e. examining the suitability of the various assumptions in concrete modeling on the seismic response of arch dams, certain simplifying assumptions are made with respect to the other parameters which influence the response in practice. In the present study, dam–water interaction effects are modeled using the conservative and computationally efficient modified WestergaardÕs method [1]. In this approach, consistent added mass matrices are introduced which can be easily combined with the mass of the dam body. Meanwhile, the foundation is assumed to be rigid, which means dam–foundation interaction is neglected. The equations of motion of the discretized dam for a specified free-field ground acceleration are € g þ ½CfU_ g þ fP ðfU gÞg ¼ fFst g fFeq g ½MfU
ð1Þ
where ½M is the mass matrix including the added mass for the impounded water, ½C, the damping matrix, fP g, the vector of stiffness forces, fFst g, the static load vector, and fFeq g, the earthquake load vector due to uniform free-field ground accelerations in x, y and z directions. In addition, fU g is the vector of unknown nodal displacements relative to the free-field input ground motions and superdots indicate time derivatives. Solution of Eq. (1) is carried out in time domain using the NewmarkÕs scheme as explained elsewhere [6].
3. Constitutive models of concrete 3.1. Non-orthogonal smeared crack model The NOSC is the first model adopted in this study. The model captures crack initiation and employs a bilinear softening tensile stress–strain relation based on fracture energy release equivalence criterion similar to the works of de Borst and Nauta [7] and Rots [8]. However, the formulation followed by them is more appropriate for considering a single crack modeling at any given location, i.e. uni-directional smeared crack model. Thus they were obliged to adopt some simplifications to extend the UDSC to a NOSC model [8]. In this regard, Channakeshava [9] improved the formulation to a computationally efficient form, which is the basis of the formulation in this study. However, the method is also modified to consider a threshold angle in the concept of initiation of new cracks. Therefore, the implemented formulation is rigorous, relatively efficient and it could form numerous cracks (up to 10) in a sampling point. The use of non-zero shear retention factor for fixed cracks as well as the cyclic nature of seismic loads implies that the axis of major principal stress rotates after crack formation. This means that magnitude of inclined principal tensile stress may significantly increase, even though the fixed crack correctly shows softening for the stress normal to the crack surface. Therefore, it must be decided whether or not a subsequent crack is allowed to initiate. In the formulation presented herein, this decision is taken based on the value of maximum tensile stress and the angle between new crack and former crack(s) directions. In this regard, a new crack initiates when: (1) the maximum tensile stress at a sampling point exceeds the tensile strength of material, and (2) the inclination angle between the maximum principal tensile stress direction and the normal vector of existing crack(s) exceeds an arbitrary threshold angle a. Based on the latter assumption, it is obvious that the value of r1 may violate the tensile strength condition and possible overstressing may occur. Further concepts of the model are summarized as follows. In the pre-cracked phase, for simplicity, a linear elastic stress–strain relation for concrete is assumed. A fundamental feature of the implemented crack model is decomposition of the total strain increment vector of cracked concrete fDeg into an uncracked concrete strain increment vector fDeco g and a crack strain increment cr vector fDecr g (Fig. 1). In this figure, rcr nn and enn are crack normal stress and crack normal strain, respectively. In the NOSC approach, the crack strain increment vector would represent all of the formed cracks in a sampling point. If ncr denotes the number of formed cracks in a sampling point, decomposition of the total strain increment vector could be written as follows:
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
1463
Fig. 1. Decomposition of cracked concrete strains.
fDeg ¼ fDeco g þ fDecr g ¼ fDeco g þ
ncr X fDecr gi
ð2Þ
i¼1
The relation between stress increment vector fDrg and total strain increment vector fDeg in a sampling point of cracked concrete can be written as [10]: ( ) 1 ncr h X cr 1 T i co 1 Te i ½D i Te i fDeg fDrg ¼ ½D þ i¼1
ð3Þ where ½Te i is a 6 3 matrix, transforming increments of the ith crack strain vector from the local to global system of coordinates. This is obtained by selecting appropriate columns of usual transformation matrix ½Te i . Meanwhile, the properties of ½Dcr i , which is referred to crack tangent operator, constitute the core of the smeared crack formulation. Another word, this matrix controls the stiffness degradation and crack opening/ closing states for the ith crack at each sampling point. In three-dimensional cases, this is a 3 3 matrix, which distinguishes between modes I and II of failure. Normally, there is no interaction assumed between these two modes of failure, thus the explicit coupling among normal and shear components in the crack traction-strain expression is neglected [8]. As a result, the crack tangent operator reduces to the diagonal form as follows: 2 3 0 Dc 0 ½Dcr i ¼ 4 0 Gcs 0 5 ð4Þ 0 0 Gct i
flects the shape of the softening curve. For a bilinear softening relation, two values of a2 are necessary to determine the slopes of the two lines or branches. The second branch starts at f^t ¼ a1 ft (Fig. 2), in which a1 is a user defined quantity. Other parameters in this figure are ft and Dsc , which are defined as reduced tensile strength and stiffness in normal crack direction on the unloadingreloading case (equal to infinity in the elastic unloading case), respectively. It must be noted, however, the values of a2 for the two branches should be specified such that the total area underneath the softening curve equals the fracture energy density of the material. For an ideal bilinear behavior, an end branch is added to the main softening curve in order to avoid zero stiffness in the crack normal direction, which is problematic in the calculation of ½Dcr 1 in relation (3). It should be noted i that for very high values of fracture energy; only the shape of softening curve is slightly changed. In this regard, the main branch of the curve would be approximately parallel to the horizontal axis. This status represents a fictitious case, which is referred to NOSC-U in this study. For smeared crack, fracture is distributed over a crack bandwidth l , which is related to the particular finite element dimensions and configuration.
where Dc is stiffness in the n-direction (normal to the crack plane), and Gcs , Gct are crack interface shear stiffness components along the crack plane (s t). These terms are defined below, respectively. Dc ¼
a2 ft2 l
2 Gf
ð5Þ
The parameter ft denotes the concrete tensile strength or crack initiation stress, and Gf is fracture energy release rate, which is defined as the amount of energy required for forming a unit area of a mode I crack. These parameters are material constants. The parameter a2 re-
Fig. 2. Bilinear idealization of strain softening branch.
1464
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
Consequently, the fracture energy should be released over this width in order to ensure objectivity with respect to mesh refinement [11]. Different alternatives are proposed for defining l . Here, the crack bandwidth (characteristic length) is assumed equal to side length of an equivalent cube with the same volume as the tributary volume at the cracked sampling point of a 20-noded isoparametric solid element [12]. The crack interface shear stiffness components, Gcs and Gct in Eq. (4) are expressed in terms of the uncracked concrete shear modulus G and so-called shear retention constants [13]: Gcs ¼
Gct ¼
ð6aÞ
bt G 1 bt
ð6bÞ
where bs and bt are the shear retention factors (constants) in the s and t directions, respectively. In this study, these factors are assumed to be equal, i.e. bs ¼ bt ¼ b. In general the value of b varies from 0 to 0.5. A retention of high shear stiffness (b close to 0.5) for crack interface can cause extensive cracking in certain cases while zero retention results in numerical instabilities in solving the governing equations numerically [14]. 3.2. Elasto-plastic model Considering the elasto-plastic behavior for solid elements, the tangent stiffness matrix of each element can be defined as: Z ½Kt ¼ ½BT ½Dep ½B dv ð7Þ v
In this equation, ½B and ½Dep are the strain–displacement and elasto-plastic rigidity matrices, respectively. The latter matrix depends on the yield surface and the assumptions being adopted. In general, the yield surface can be written as: ð8Þ
where frg is the stress vector and j is the hardening parameter, which controls expansion of the yield surface and utilized in defining subsequent yield surfaces. In the case of the associate flow rule, elasto-plastic rigidity matrix would be as follows: ! ½DfagfagT ½D ½Dep ¼ ½D ð9Þ A þ fagT ½Dfag In this relation, ½D is the linear elastic rigidity matrix and the other parameters are defined as follows: A¼
oF dj oj dk
oF ofrg
ð10aÞ
ð10bÞ
In this study, Mohr–Coulomb criterion is used to define the yield surface. Meanwhile, perfect plasticity idealization is adopted, i.e. A ¼ 0. The equation of this yield surface may be written as: pffiffiffiffi 1 p F ðfrgÞ ¼ I1 sin / þ J2 sin h þ 3 3 1 p sin / c cos / þ pffiffiffi cos h þ 3 3 ¼0
bs G 1 bs
F ðfrg; jÞ ¼ f ðfrgÞ KðjÞ ¼ 0
fag ¼
ð11Þ
The parameters c, / are cohesion and angle of internal friction, respectively. Furthermore, I1 , J2 and h are first invariant of stress tensor, second invariant of deviatoric stress tensor and angle of similarity, respectively. The latter parameter is being assumed as zero on the projection of r1 axis on the deviatoric plane. Thus, the angle of similarity varies between 0 and p/3. It is also worthwhile to mention that increments of effective plastic strain are calculated based on plastic works and effective stresses. The effective plastic strain can be used as a criterion for determining the amount of plastic behavior (damage) in the elasto-plastic model. Detailed description related to the method of calculation of each parameter and other aspects can be found elsewhere [6,10,15].
4. Verification of proposed models The proposed models are implemented in a special finite element program, called SNAP [16]. The program solves the equations of motion incrementally using either load or displacement control. To verify the mentioned models and their accuracy, a single edge notched beam with 102 102 cross section and 788 mm length subjected to three point static bending is analyzed. The notch-to-depth ratio for the beam is 0.5. The obtained results are compared with the corresponding experimental and numerical results published by Malvar and Fourney [17]. The concrete properties are taken from their work, which are obtained from ancillary tests on appropriate test specimens. The reported properties are modulus of elasticity E ¼ 21:7 GPa, tensile strength at 28 days ft ¼ 3:1 MPa, compressive strength at 28 days, fc0 ¼ 29 MPa, fracture energy Gf ¼ 0:0763 N/mm, and Poisson ratio m ¼ 0:2. In the present analysis, 20-node isoparametric solid elements with full Gaussian integration, i.e. 3 3 3, are used. Only one quarter of the beam is discretized due to double symmetry. The three-dimensional finite element mesh for the quarter beam comprises 152 elements as shown in Fig. 3a. The incremental analysis was per-
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
1465
Fig. 3. Notch beam verification example.
formed using displacement control. In addition, it is theoretically known that no shear stresses should be formed in the central section of the symmetric beam. Therefore, a very low value of shear retention factor (b ¼ 0:0001) is adopted for the NOSC model in this example. Meanwhile, crack bandwidth is assumed to be equal to 20 mm, which is 2.1 times the maximum aggregate size, i.e. 9.5 mm, consistent with the mentioned reference. Furthermore, The amounts of cohesion and internal angle of friction for elasto-plastic model can be easily calculated from the uniaxial tensile and compressive strengths of concrete. Herein, these would be c ¼ 4:74 MPa and / ¼ 53:79°. Malvar and Fourney [17] considered two cases in their analyses, case LE-CS with 500 linear elastic and 20 concrete (inelastic) elements and case C-CS composed of inelastic concrete elements only. They used CornelissenÕs non-linear softening model [18] and the 2 2 2 Gaussian integration scheme to compensate for the excessive stiffness of their elements. In Fig. 3b, the load-mid span deflection diagrams of the beam obtained by the ELPL and NOSC models are compared with the corresponding experimental curve and with the analytical results of [17]. In the NOSC model, the peak load predicted by SNAP is close to the experimental value (1.62% difference), and the total energy needed for fracturing is approximately the same. On the other hand in the ELPL model, although the maximum tensile principal stresses are limited to the tensile strength of concrete, the amount of force is raised with increasing the mid-span vertical displacement. This behavior is basically occurred based on the lack of appropriate mechanism for modeling the softening phenomenon in the elasto-plastic models. Meanwhile, the predicted initial stiffness of beam by SNAP is slightly lower than the experimental value. Compared to the Malvar and FourneyÕs model, the present analysis gives better results for NOSC model in this problem. Additional verifying examples can be
found elsewhere [10] which are not discussed here due to space limitations.
5. Analysis of Shahid Rajaee arch dam Consider the idealized symmetric model of Shahid Rajaee concrete arch dam as shown in Fig. 4. The dam is 130 m high with a crest length of 420 m. It is located in the north of Iran, in the seismically active foothills of Alborz Mountains near the city of Sari. The damÕs seismic response is analyzed by using a linear and three non-linear models described in the following. Let LIN designates the linear analysis case, used mainly for comparative purposes, NOSC-O and NOSCU the non-linear analysis cases using the non-orthogonal smeared crack model with ordinary and unlimited amounts of fracture energy, respectively. Meanwhile, ELPL is the non-linear analysis case corresponding to the elasto-plastic model. In all these cases, the same finite element discretization is utilized, consisting of 487 nodes and 76 isoparametric 20-noded solid elements in two layers through the thickness of the dam. The earthquake excitations include two components of the Friuli–Tolmezzo earthquake (the cross-canyon component is neglected to preserve symmetric conditions) whose records are normalized based on the frequency content for maximum design earthquake (MDE) condition with the peak ground acceleration of 0.56 g (Fig. 5). The Rayleigh damping coefficients were determined such that equivalent damping for frequencies close to the first and sixth modes of vibration would be 10% of the critical damping. The time integration is performed for 12 s, which covers the main response of the structure. As mentioned in Section 2, the foundation is assumed rigid to keep the computation time realistic. Although this assumption would highly influence the boundary stresses in the linear case, it is believed to be less important for the non-linear cases. In addition, the
1466
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
Fig. 4. Idealized model of Shahid Rajaee concrete arch dam.
Fig. 5. Accelerograms of applied earthquake (scale factor for MDE ¼ 9:81).
assumption of rigid foundation has less effect on the stresses in the vicinity of the crest spillway, where the major failure is expected to occur. As discussed earlier, the dam–water interaction effects are included by the modified WestergaardÕs method. Although different concrete mixes were used in the dam construction, the material properties are assumed uniform. For this study, the basic material properties are: E ¼ 30 GPa, q ¼ 2400 kg/m3 and m ¼ 0:18. Meanwhile, for the NOSC cases the concrete tensile strength is assumed as ft ¼ 3:0 MPa, which is more representative of tensile strength of mass concrete of dam body. The fracture energy release rate is assumed as Gf ¼ 2400 N/m and Gf ¼ 2400 106 N/m for the NOSC-O and NOSC-U cases, respectively. It should be mentioned that high discrepancy has been reported in the amount of fracture energy of concrete (Gf ) based on type of aggregates, cements and additives. In this regard, RILEM [19] shows a variation in the value of Gf from 70 to 200 N/m. Meanwhile, The results of a study performed on the concrete of existing dams in Switzerland show the
amount of Gf for dam concrete is two to three times higher than ordinary concrete [20]. Furthermore, in practical studies, researchers have often utilized relatively higher values of Gf . For instance, Vargas-Loli and Fenves [21] used Gf ¼ 2200 N/m for analysis of Pine Flat dam and Feltrin et al. [22] assumed Gf ¼ 1375 N/m. It must be also mentioned that it is not possible to lower the amount of fracture energy release rate considerably, due to the applied mesh sizes. Otherwise, the amount of crack tangent operator will be increased and if the value exceeds the elastic modulus of material, correct softening behavior is not achievable. Therefore, to study the effect of this parameter on the behavior of the dam, it is decided to raise this value significantly instead. The threshold angle between consecutive cracks in the NOSC cases is assumed to be equal to 15° and shear retention factors for all of the smeared crack models are considered constant and equal to b ¼ 0:1. A bilinear softening curve with constant parameters, a1 ¼ 0:01 and a2 ¼ 1 for the first branch, and a2 ¼ 0:0001 for the second branch is used. In addition, the cohesion and angle of
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
internal friction of concrete are determined based on the uniaxial compressive and tensile strengths for the ELPL case. These parameters are c ¼ 4:74 MPa and / ¼ 54:9°. The analysis commenced by applying the static loads. These include the dead weight of the dam body and the hydrostatic water pressures corresponding to 122 m height of reservoir water. These loads were incrementally increased in time until they reached their full values. Three time steps were initially considered, i.e. 0.01, 0.005 and 0.0025 s, and it was noticed that the results were almost the same for Dt ¼ 0:005 and 0.0025 s. Therefore, the former value was selected for time integration [10]. In this respect, the dead load is applied in 20 increments and the hydrostatic pressures thereafter in 20 increments at negative range of time. At time zero, the actual non-linear dynamic analysis is started with static displacements and stresses being applied as initial conditions.
1467
Table 1 Maximum displacements at dam crest (mm) Case
Right quarter point
Center point
U
V
W
V
W
LIN NOSC-O NOSC-U ELPL
)28.5 )26.2 )27.6 )27.3
64.2 64.0 64.2 63.5
5.4 7.5 4.9 5.0
117.1 )120.9 110.4 114.0
)13.1 64.0 17.3 17.0
5.1. Effect of fracture energy As mentioned earlier, it was not possible to decrease the amount of fracture energy significantly, based on the selected mesh size and amount of tensile strength of concrete. Therefore, to study the influence of fracture energy on the behavior the dam, it was decided to increase the amount of this parameter. Meanwhile, a fictitious case for smeared crack approach can be obtained by selecting a very large value for Gf . The basic physical characteristic of such a definition is making a constant maximum tensile principal stress, which is equivalent to the tensile strength of material, for all values of crack normal strain. This arrangement is almost similar to the elasto-perfect plastic behavior in tensile region and elastic behavior on compressive region. In this regard, the elastic unloading option is also utilized for the NOSC cases. To study the displacements of the dam, two nodes are selected on the dam crest. These are located at the center and at approximately the right quarter point. The maximum displacements through time for each of the x, y, and z directions at these nodes are summarized in Table 1, i.e. U , V , and W , respectively. The analysis did not reveal large displacements in the structure under the severe earthquake excitations considered, giving an initial indication of a safe design of the dam based on the adopted theories and assumptions. In Figs. 6 and 7, displacement histories of mid-crest point for the NOSC cases are compared with the linear case, in the stream and vertical directions, respectively. Comparison among the results of the smeared crack models and those of the linear model shows a drift in the crest displacements in upward and upstream directions. This is the main characteristic of displacement histories in smeared crack approach, which is due to extensive cracking of downstream face of the dam beneath the
Fig. 6. Comparison of displacements in the stream direction at mid-crest point among the linear and smeared crack models.
Fig. 7. Comparison of displacements in the vertical direction at mid-crest point among the linear and smeared crack models.
spillway block. When most of these cracks are in open state, they produce a cumulative crack strain, which is the major source of the drift. However, the results of the NOSC-U case are settled between the two other cases and the amounts of drift in stream and vertical directions are significantly less than the NOSC-O case. Stiffer behavior in the NOSC-U case in comparison with the NOSC-O case is the major cause of such response. Since
1468
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
in the NOSC-U case the tensile strength of the material is almost constant for all of the values of crack normal strain, the load bearing capacity of the material is almost the same during the analysis. Therefore, in this case the material is stiffer than the NOSC-O case and it shows lower tendency toward opening of cracks and increasing the crack strain. The envelope of maximum principal stresses throughout the analysis is displayed in Figs. 8–15, while a summary of the results is obtained and tabulated in Table 2. In the linear analysis, very high tensile stresses are calculated at the base of the dam (Fig. 8), due to the rigid foundation assumption. In addition, high tensile stresses occur in the spillway region, which in reality are released with the opening of contraction joints. A maximum compressive stress of )15.25 MPa is obtained in the upstream face near the spillway region (Fig. 12). The NOSC model can completely restrict the amount of tensile stresses. In this regard, no overstressing in tensile stresses is observed in the NOSC-O case, which means that tensile stresses are limited to the tensile strength of concrete. However, slight overstressing of tensile stresses (less than 6%) are still noticed in the dam body for the NOSC-U case, due to threshold angle condition for initiation of subsequent crack(s). As mentioned earlier, a new crack is initiated at a sampling point when not only the maximum principal stress exceeds the tensile strength but also the angle between the new crack and the existing cracks is more than the threshold angle.
Comparison among the principal (tensile and compressive) stress results shows that the general patterns of contours are very similar in both above-mentioned NOSC cases. The extent of softening obtained from different regions of the smeared crack models is shown in plots of reduced tensile strength at the end of analysis (Figs. 16– 18). These figures together with the crack patterns (Figs. 19 and 20), which actually show the normal vector to each crack plane, illustrate the extent of damage occurring in the dam body. It is noticed that the extension of cracking in both NOSC cases, i.e. NOSC-O and NOSC-U, is almost the same. However, the cracked zone in the NOSC-U case is more extended in upstream face of spillway region than the NOSC-O case, conversely for the base region of the dam, in the downstream face the extension of damaged zone is higher for the NOSC-O case. In addition, the minimum tensile strength of material is not changed through the analysis for the NOSC-U case (Fig. 17), due to its special form of softening curve. In this case, for all of the values of crack strain, the amount of crack normal stress is constant and equal to the initial tensile strength of concrete.
5.2. Elasto-plastic model In this section, the results of elasto-plastic model (ELPL case) is studied and compared with the linear
Fig. 8. Envelope of maximum tensile stresses (MPa) for the LIN.
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
1469
Fig. 9. Envelope of maximum tensile stresses (MPa) for the NOSC-O.
Fig. 10. Envelope of maximum tensile stresses (MPa) for the NOSC-U.
model (LIN case) and a fictitious state of NOSC model (NOSC-U case). With this arrangement, it is possible to study the elasto-plastic behavior of the material and compare it with linear elastic and a special form of smeared crack approach.
As mentioned earlier, the central and the right quarter points of the dam crest are selected to monitor the displacement history of the dam. The maximum displacements through time at these nodes are summarized in Table 1.
1470
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
Fig. 11. Envelope of maximum tensile stresses (MPa) for the ELPL.
Fig. 12. Envelope of maximum compressive stresses (MPa) for the LIN (U/S view).
Fig. 13. Envelope of maximum compressive stresses (MPa) for the NOSC-O (U/S view).
It is observed that large displacements are not occurring in the structure under the severe earthquake in the ELPL case similar to the other cases. In Figs. 21 and
22, displacement histories of mid-crest point for the NOSC-U case are compared with the ELPL case, in the stream and vertical directions, respectively. Based on
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
1471
Fig. 14. Envelope of maximum compressive stresses (MPa) for the NOSC-U.
Fig. 15. Envelope of maximum compressive stresses (MPa) for the ELPL.
the results, the displacements and amount of drift in the ELPL case are very close to the fictitious state of smeared crack approach, i.e. NOSC-U. The envelope of maximum principal stresses throughout the analysis are obtained and summarized in
Table 2 while the principal stress envelopes are displayed in Figs. 11 and 15. The linear and fictitious state of the NOSC cases are discussed in the former section and the elasto-plastic model is compared with the mentioned models in this section.
1472
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
Table 2 Maximum tensile and compressive stresses (MPa) Case
LIN NOSC-O NOSC-U ELPL
Maximum tensile stress (r1 )
Maximum compressive stress (r3 )
Spillway and central region
Base and abutments
Spillway and central region
Base and abutments
U/S
D/S
U/S
D/S
U/S
D/S
U/S
D/S
8.93 3.00 3.17 3.06
8.93 3.00 3.17 3.06
15.76 3.00 3.17 3.06
2.09 3.00 1.60 1.53
)15.25 )15.54 )15.26 )15.46
)13.58 )13.85 )15.26 )15.46
)6.91 )15.54 )13.60 )13.77
)11.92 )12.17 )11.93 )12.09
Fig. 16. Reduced tensile strength (MPa) at the end of analysis for the NOSC-O.
Fig. 17. Reduced tensile strength (MPa) at the end of analysis for the NOSC-U.
Fig. 18. Effective plastic strain at the end of analysis for the ELPL.
Fig. 19. Crack pattern at the end of analysis for the NOSC-O.
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
1473
Fig. 20. Crack pattern at the end of analysis for the NOSC-U.
Fig. 21. Comparison of displacements in the stream direction at mid-crest point between fictitious state of smeared crack and elasto-plastic models.
stress distribution and maximum value of tensile stresses are significantly varied in comparison with linear case, as the amounts of tensile stresses are limited in the spillway and base regions to a maximum value close to tensile strength of the material. Meanwhile, the stress distribution is similar to the NOSC-U case with maximum difference of less than 4%. For the ELPL case, it is noticed that the distribution of maximum compressive stresses is similar to the linear analysis in some regions of the dam. Moreover, the stress results are very close to the NOSC-U case through the dam body. This demonstrates that performance of the elasto-plastic model in compressive regions is similar to the NOSC-U case, and state of stresses has not reached the plastic surface in the compressive region. Therefore for compressive stresses, the conduct of the dam in all analyses is very close to linear elastic behavior under the presumed earthquake excitations. The contours of effective plastic strain are illustrated in Fig. 18, to indicate the extension of the plastic zones through the dam body. As it is depicted in this figure, the highest values of effective plastic strain are produced in the upstream face of the dam base, and afterwards in the spillway region, of course in lower magnitudes. 6. Conclusions
Fig. 22. Comparison of displacements in the vertical direction at mid-crest point between fictitious state of smeared crack and elasto-plastic models.
Based on the obtained results, the ELPL model can appropriately restrict the amount of tensile stresses. In this regard, the maximum tensile stresses through dam body are limited to r1 ¼ 3:06 MPa, which is very close to the tensile strength of mass concrete. In addition,
A computer program is developed for the non-linear dynamic analysis of concrete arch dams based on the fixed smeared crack approach and elasto-plastic modeling of the material. The program is initially verified by analyzing a notched concrete beam and comparing the results with the corresponding experimental data. Subsequently, an idealized symmetric model of Shahid Rajaee arch dam is considered and four cases are analyzed: a linear case (LIN), a non-orthogonal smeared crack with ordinary fracture energy (NOSC-O), a nonorthogonal smeared crack with a very high value of fracture energy (NOSC-U), and an elasto-plastic model (ELPL). The main conclusions of the analyses are: A drift in upward and upstream directions forms the prominent characteristic of histories of the crest displacements in the smeared crack and elasto-plastic models.
1474
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
The non-orthogonal smeared crack and elasto-plastic models can overcome the overstressing problem encountered in the analysis by linear model and predict realistic distribution and levels of stresses in the dam body. The stress results for each time step are changed when the amount of fracture energy is significantly increased. Correspondingly, in the results of envelopes of the principal stresses, some differences between the amounts of stresses can be seen for the different zones of dam body. However, general shape of the contours of these envelopes, indicate no major change between ordinary and fictitious smeared crack approaches. In this regard, the displacement histories of the monitoring points on the dam crest show a considerable differences and results of the fictitious model are inclined toward the linear model. The comparison between the ordinary form of smeared crack and elasto-plastic models reveals major differences in the results. However, comparison between the fictitious form of smeared crack and elasto-plastic models shows very similar results. In this regard, both envelopes of principal stresses, as well as time history of displacements at dam crest is very much similar. This indicates that the major part of plastic behavior is occurring in the tensile region. It must be finally emphasized that knowing the deficiencies of usual elasto-plastic models in tensile region in comparison with fracture mechanics based models, it is recommended to use fracture based models for dynamic analysis of concrete arch dams.
Acknowledgement This research investigation was supported by grant no. 31305301 from the National Research Council of Iran to the Amirkabir University, Tehran, Iran. The support is gratefully acknowledged.
References [1] Kuo JSH. On the nonlinear dynamic response of arch dams to earthquakes––I. Fluid–structure interaction: Added-mass computations for incompressible fluid. II. Joint opening nonlinear mechanism: Interface smeared crack model. PhD thesis. University of California, Berkeley, 1982. [2] Dungar R. A viscoplastic bounding surface model for concrete under compressive and tensile conditions and its application to arch dam analysis. Proceedings of Second International Conference on Constitutive Laws for Engineering Materials. Tucson/ AZ. 1987. p. 849–56.
[3] Cervera M, Oliver J, Faria R. Seismic evaluation of concrete dams via continuum damage models. Earthquake Eng Struct Dyn 1995;24:1225–45. [4] Noruziaan B. Nonlinear seismic analysis of concrete arch dams. PhD thesis. Department of Civil and Environmental Engineering. Carleton University. Ottawa, Canada, 1995. [5] Hall JF. Efficient non-linear seismic analysis of arch dams. Earthquake Eng Struct Dyn 1998;27:1425–44. [6] Owen DRJ, Hinton E. Finite element in plasticity. Swansea: Pineridge Press Ltd; 1980. [7] de Borst R, Nauta P. Non-orthogonal cracks in a smeared finite element model. Eng Comput 1985;2(3):35–46. [8] Rots JG. Computational modelling of concrete fracture. Doctorate dissertation. TH Delft. The Netherlands, 1988. [9] Channakeshava C. Elasto-plastic-cracking analysis of plain and reinforced concrete structures. PhD dissertation. Indian Institute of Science, Bangalore, India, 1987. [10] Espandar R. A study on nonlinear dynamic response of concrete arch dams. PhD dissertation. Department of Civil Engineering, Amirkabir University of Technology, Tehran, Iran, 2001. [11] Bazant ZP, Oh BH. Crack band theory for fracture of concrete. Mater Struct RILEM 1983;16(93):155– 77. [12] Cervera M. Nonlinear analysis of reinforced concrete structures using three dimensional and shell finite element models. PhD dissertation. Department of Civil Engineering. University College of Swansea, Wales, UK, 1986. [13] Rots JG, Blaauwendraad J. Crack models for concrete–– Discrete or smeared? Fixed, multi-directional or rotating? Heron 1989;34(1). [14] Barzegar F, Maddipudi S. Three-dimensional modeling of concrete structures. J Struct Eng ASCE 1997;123(10):1339– 56. [15] Zienkiewicz OC, Taylor RL. The finite element method. 4th ed. Maidenhead, Berkshire: McGraw-Hill; 1989. [16] Espandar R. SNAP: Seismic nonlinear analysis program. Tehran, Iran: Amirkabir University of Technology; 2000. [17] Malvar LJ, Fourney ME. A three dimensional application of the smeared crack approach. Eng Fract Mech 1990;35(1/ 2/3):251–60. [18] Cornelissen HAW, Hordijk DA, Reinhardt HW. Experiments and theory for the application of fracture mechanics to normal and lightweight concrete. Proceedings of International Conference on Fracture Mechanics of Concrete, Lausanne, Switzerland, 1985. [19] RILEM, Technical Committee 90 FMA. Fracture mechanics of concrete structures. Chapman and Hall Ltd, 1989. [20] Bruhwiler E, Wittmann FH. Failure of dam concrete subjected to seismic load conditions. Eng Fract Mech 1990;35(1/2/3):565–71. [21] Vargas-Loli L, Fenves GL. Effects of concrete cracking on the earthquake response of gravity dams. Earthquake Eng Struct Dyn 1989;18:575–92. [22] Feltrin G, Galli M, Bachmann H. Influence of cracking on the earthquake response of concrete gravity dam with reservoir. Proceedings of the 10th World Conference on Earthquake Engineering, vol. 8, Madrid, Spain, 1992. p. 4627–32.
Comparison of non-orthogonal smeared crack and plasticity models for dynamic analysis of concrete arch dams Radin Espandar a
a,*
, Vahid Lotfi
b
Department of Earth Sciences, Shahid Beheshti University, P.O. Box 13145-165, Tehran, Iran b Department of Civil and Environmental Engineering, Amirkabir University, Tehran, Iran Received 19 April 2002; accepted 20 January 2003
Abstract A special three-dimensional finite element program is developed to study the seismic response of concrete arch dams. Two types of continuum mechanics non-linear techniques are incorporated in the program, i.e. non-orthogonal smeared crack and elasto-plastic models. The formulation is presented initially and accuracy of the models and their correct implementations are verified by analyzing a notched beam and comparing its results with available experimental data. Following the satisfactory results of the beam analysis, the program is applied to a practical problem, viz. the non-linear behavior of the 130 m high Shahid Rajaee arch dam in Iran, subjected to the Friuli–Tolmezzo earthquake. Two cases of non-orthogonal smeared crack model are implemented for the dam, namely, ordinary and fictitious approaches. Meanwhile, an elastic–perfectly plastic model for mass concrete is also applied. The results are compared with those obtained from a linear elastic analysis under the same excitation. Based on the results, it is concluded that the nonorthogonal smeared crack approach can redistribute the state of stresses and produces a more realistic profile of stresses in the dam. Further, the elasto-plastic model could reduce significantly the high amount of overstressing noticed in the linear analysis. However, the elasto-plastic model is not perceptibly activated in the compressive state of stresses. It is also observed that a drift in the crest displacements is the prominent consequence of cracking or plasticity of the dam. Ó 2003 Elsevier Science Ltd. All rights reserved. Keywords: Non-linear dynamic analysis; Concrete arch dam; Smeared crack approach; Plasticity; Concrete modeling; Seismic load
1. Introduction There are several approaches to model the complicated stress–strain behavior of concrete under multiaxial stress states. In this regard, the theories based on plasticity and fracture mechanics are commonly applied in most of the engineering analysis. However, only few research works have been performed to study the dynamic behavior of concrete arch dams with concrete modeling. Kuo [1] suggested interface smeared crack
* Corresponding author. Tel.: +98-911-226-0861; fax: +9821-831-0043. E-mail addresses: [email protected] (R. Espandar), vahlotfi@aut.ac.ir (V. Lotfi).
approach to model the contraction joints and applied this technique for dynamic analysis of arch dams. Dungar [2] used a bounding surface model for static analysis of an arch dam. Cervera et al. [3] utilized a continuum damage model to analyze Talvacchia arch dam in Italy. In another research work, Noruziaan [4] applied bounding surface and orthogonal smeared crack models for compressive and tensile regions of stresses, respectively. However, it was reported that numerical problems occurred in some cases. The most recent related research work belongs to Hall [5], which a simple smeared crack model is utilized for modeling the contraction and construction joints. However, the diagonal cracks are neglected in that work. In the continuum models, such as elasto-plastic and smeared crack models, the dam body is assumed as a
0045-7949/03/$ - see front matter Ó 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0045-7949(03)00083-X
1462
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
monolith and the non-linear behavior is limited to the mass concrete of the dam body. The behavior of dam is investigated by studying the displacement histories and envelopes of principal stresses. The mentioned parameters are studied qualitatively, due to lack of any standard. In this manner, displacements of the dam are reasonable when they are not so large to endanger the general stability of the structure. Also, overstressing in tensile and compressive stresses, i.e. stresses higher than strengths, is used as a criterion for evaluating effectiveness of models in correct redistribution of stresses throughout the dam body. Furthermore, reasonable conclusions can be deduced by comparison among different models.
2. Description of analysis procedure The earthquake response of concrete arch dams is affected by dam–water and dam–foundation interactions as well as material non-linearities and joint movements. Although the interaction effects can be rigorously included in a frequency-domain analysis for a linear model, a time domain analysis is necessary for nonlinear models. In compliance with the main purpose and scope of the present work, i.e. examining the suitability of the various assumptions in concrete modeling on the seismic response of arch dams, certain simplifying assumptions are made with respect to the other parameters which influence the response in practice. In the present study, dam–water interaction effects are modeled using the conservative and computationally efficient modified WestergaardÕs method [1]. In this approach, consistent added mass matrices are introduced which can be easily combined with the mass of the dam body. Meanwhile, the foundation is assumed to be rigid, which means dam–foundation interaction is neglected. The equations of motion of the discretized dam for a specified free-field ground acceleration are € g þ ½CfU_ g þ fP ðfU gÞg ¼ fFst g fFeq g ½MfU
ð1Þ
where ½M is the mass matrix including the added mass for the impounded water, ½C, the damping matrix, fP g, the vector of stiffness forces, fFst g, the static load vector, and fFeq g, the earthquake load vector due to uniform free-field ground accelerations in x, y and z directions. In addition, fU g is the vector of unknown nodal displacements relative to the free-field input ground motions and superdots indicate time derivatives. Solution of Eq. (1) is carried out in time domain using the NewmarkÕs scheme as explained elsewhere [6].
3. Constitutive models of concrete 3.1. Non-orthogonal smeared crack model The NOSC is the first model adopted in this study. The model captures crack initiation and employs a bilinear softening tensile stress–strain relation based on fracture energy release equivalence criterion similar to the works of de Borst and Nauta [7] and Rots [8]. However, the formulation followed by them is more appropriate for considering a single crack modeling at any given location, i.e. uni-directional smeared crack model. Thus they were obliged to adopt some simplifications to extend the UDSC to a NOSC model [8]. In this regard, Channakeshava [9] improved the formulation to a computationally efficient form, which is the basis of the formulation in this study. However, the method is also modified to consider a threshold angle in the concept of initiation of new cracks. Therefore, the implemented formulation is rigorous, relatively efficient and it could form numerous cracks (up to 10) in a sampling point. The use of non-zero shear retention factor for fixed cracks as well as the cyclic nature of seismic loads implies that the axis of major principal stress rotates after crack formation. This means that magnitude of inclined principal tensile stress may significantly increase, even though the fixed crack correctly shows softening for the stress normal to the crack surface. Therefore, it must be decided whether or not a subsequent crack is allowed to initiate. In the formulation presented herein, this decision is taken based on the value of maximum tensile stress and the angle between new crack and former crack(s) directions. In this regard, a new crack initiates when: (1) the maximum tensile stress at a sampling point exceeds the tensile strength of material, and (2) the inclination angle between the maximum principal tensile stress direction and the normal vector of existing crack(s) exceeds an arbitrary threshold angle a. Based on the latter assumption, it is obvious that the value of r1 may violate the tensile strength condition and possible overstressing may occur. Further concepts of the model are summarized as follows. In the pre-cracked phase, for simplicity, a linear elastic stress–strain relation for concrete is assumed. A fundamental feature of the implemented crack model is decomposition of the total strain increment vector of cracked concrete fDeg into an uncracked concrete strain increment vector fDeco g and a crack strain increment cr vector fDecr g (Fig. 1). In this figure, rcr nn and enn are crack normal stress and crack normal strain, respectively. In the NOSC approach, the crack strain increment vector would represent all of the formed cracks in a sampling point. If ncr denotes the number of formed cracks in a sampling point, decomposition of the total strain increment vector could be written as follows:
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
1463
Fig. 1. Decomposition of cracked concrete strains.
fDeg ¼ fDeco g þ fDecr g ¼ fDeco g þ
ncr X fDecr gi
ð2Þ
i¼1
The relation between stress increment vector fDrg and total strain increment vector fDeg in a sampling point of cracked concrete can be written as [10]: ( ) 1 ncr h X cr 1 T i co 1 Te i ½D i Te i fDeg fDrg ¼ ½D þ i¼1
ð3Þ where ½Te i is a 6 3 matrix, transforming increments of the ith crack strain vector from the local to global system of coordinates. This is obtained by selecting appropriate columns of usual transformation matrix ½Te i . Meanwhile, the properties of ½Dcr i , which is referred to crack tangent operator, constitute the core of the smeared crack formulation. Another word, this matrix controls the stiffness degradation and crack opening/ closing states for the ith crack at each sampling point. In three-dimensional cases, this is a 3 3 matrix, which distinguishes between modes I and II of failure. Normally, there is no interaction assumed between these two modes of failure, thus the explicit coupling among normal and shear components in the crack traction-strain expression is neglected [8]. As a result, the crack tangent operator reduces to the diagonal form as follows: 2 3 0 Dc 0 ½Dcr i ¼ 4 0 Gcs 0 5 ð4Þ 0 0 Gct i
flects the shape of the softening curve. For a bilinear softening relation, two values of a2 are necessary to determine the slopes of the two lines or branches. The second branch starts at f^t ¼ a1 ft (Fig. 2), in which a1 is a user defined quantity. Other parameters in this figure are ft and Dsc , which are defined as reduced tensile strength and stiffness in normal crack direction on the unloadingreloading case (equal to infinity in the elastic unloading case), respectively. It must be noted, however, the values of a2 for the two branches should be specified such that the total area underneath the softening curve equals the fracture energy density of the material. For an ideal bilinear behavior, an end branch is added to the main softening curve in order to avoid zero stiffness in the crack normal direction, which is problematic in the calculation of ½Dcr 1 in relation (3). It should be noted i that for very high values of fracture energy; only the shape of softening curve is slightly changed. In this regard, the main branch of the curve would be approximately parallel to the horizontal axis. This status represents a fictitious case, which is referred to NOSC-U in this study. For smeared crack, fracture is distributed over a crack bandwidth l , which is related to the particular finite element dimensions and configuration.
where Dc is stiffness in the n-direction (normal to the crack plane), and Gcs , Gct are crack interface shear stiffness components along the crack plane (s t). These terms are defined below, respectively. Dc ¼
a2 ft2 l
2 Gf
ð5Þ
The parameter ft denotes the concrete tensile strength or crack initiation stress, and Gf is fracture energy release rate, which is defined as the amount of energy required for forming a unit area of a mode I crack. These parameters are material constants. The parameter a2 re-
Fig. 2. Bilinear idealization of strain softening branch.
1464
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
Consequently, the fracture energy should be released over this width in order to ensure objectivity with respect to mesh refinement [11]. Different alternatives are proposed for defining l . Here, the crack bandwidth (characteristic length) is assumed equal to side length of an equivalent cube with the same volume as the tributary volume at the cracked sampling point of a 20-noded isoparametric solid element [12]. The crack interface shear stiffness components, Gcs and Gct in Eq. (4) are expressed in terms of the uncracked concrete shear modulus G and so-called shear retention constants [13]: Gcs ¼
Gct ¼
ð6aÞ
bt G 1 bt
ð6bÞ
where bs and bt are the shear retention factors (constants) in the s and t directions, respectively. In this study, these factors are assumed to be equal, i.e. bs ¼ bt ¼ b. In general the value of b varies from 0 to 0.5. A retention of high shear stiffness (b close to 0.5) for crack interface can cause extensive cracking in certain cases while zero retention results in numerical instabilities in solving the governing equations numerically [14]. 3.2. Elasto-plastic model Considering the elasto-plastic behavior for solid elements, the tangent stiffness matrix of each element can be defined as: Z ½Kt ¼ ½BT ½Dep ½B dv ð7Þ v
In this equation, ½B and ½Dep are the strain–displacement and elasto-plastic rigidity matrices, respectively. The latter matrix depends on the yield surface and the assumptions being adopted. In general, the yield surface can be written as: ð8Þ
where frg is the stress vector and j is the hardening parameter, which controls expansion of the yield surface and utilized in defining subsequent yield surfaces. In the case of the associate flow rule, elasto-plastic rigidity matrix would be as follows: ! ½DfagfagT ½D ½Dep ¼ ½D ð9Þ A þ fagT ½Dfag In this relation, ½D is the linear elastic rigidity matrix and the other parameters are defined as follows: A¼
oF dj oj dk
oF ofrg
ð10aÞ
ð10bÞ
In this study, Mohr–Coulomb criterion is used to define the yield surface. Meanwhile, perfect plasticity idealization is adopted, i.e. A ¼ 0. The equation of this yield surface may be written as: pffiffiffiffi 1 p F ðfrgÞ ¼ I1 sin / þ J2 sin h þ 3 3 1 p sin / c cos / þ pffiffiffi cos h þ 3 3 ¼0
bs G 1 bs
F ðfrg; jÞ ¼ f ðfrgÞ KðjÞ ¼ 0
fag ¼
ð11Þ
The parameters c, / are cohesion and angle of internal friction, respectively. Furthermore, I1 , J2 and h are first invariant of stress tensor, second invariant of deviatoric stress tensor and angle of similarity, respectively. The latter parameter is being assumed as zero on the projection of r1 axis on the deviatoric plane. Thus, the angle of similarity varies between 0 and p/3. It is also worthwhile to mention that increments of effective plastic strain are calculated based on plastic works and effective stresses. The effective plastic strain can be used as a criterion for determining the amount of plastic behavior (damage) in the elasto-plastic model. Detailed description related to the method of calculation of each parameter and other aspects can be found elsewhere [6,10,15].
4. Verification of proposed models The proposed models are implemented in a special finite element program, called SNAP [16]. The program solves the equations of motion incrementally using either load or displacement control. To verify the mentioned models and their accuracy, a single edge notched beam with 102 102 cross section and 788 mm length subjected to three point static bending is analyzed. The notch-to-depth ratio for the beam is 0.5. The obtained results are compared with the corresponding experimental and numerical results published by Malvar and Fourney [17]. The concrete properties are taken from their work, which are obtained from ancillary tests on appropriate test specimens. The reported properties are modulus of elasticity E ¼ 21:7 GPa, tensile strength at 28 days ft ¼ 3:1 MPa, compressive strength at 28 days, fc0 ¼ 29 MPa, fracture energy Gf ¼ 0:0763 N/mm, and Poisson ratio m ¼ 0:2. In the present analysis, 20-node isoparametric solid elements with full Gaussian integration, i.e. 3 3 3, are used. Only one quarter of the beam is discretized due to double symmetry. The three-dimensional finite element mesh for the quarter beam comprises 152 elements as shown in Fig. 3a. The incremental analysis was per-
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
1465
Fig. 3. Notch beam verification example.
formed using displacement control. In addition, it is theoretically known that no shear stresses should be formed in the central section of the symmetric beam. Therefore, a very low value of shear retention factor (b ¼ 0:0001) is adopted for the NOSC model in this example. Meanwhile, crack bandwidth is assumed to be equal to 20 mm, which is 2.1 times the maximum aggregate size, i.e. 9.5 mm, consistent with the mentioned reference. Furthermore, The amounts of cohesion and internal angle of friction for elasto-plastic model can be easily calculated from the uniaxial tensile and compressive strengths of concrete. Herein, these would be c ¼ 4:74 MPa and / ¼ 53:79°. Malvar and Fourney [17] considered two cases in their analyses, case LE-CS with 500 linear elastic and 20 concrete (inelastic) elements and case C-CS composed of inelastic concrete elements only. They used CornelissenÕs non-linear softening model [18] and the 2 2 2 Gaussian integration scheme to compensate for the excessive stiffness of their elements. In Fig. 3b, the load-mid span deflection diagrams of the beam obtained by the ELPL and NOSC models are compared with the corresponding experimental curve and with the analytical results of [17]. In the NOSC model, the peak load predicted by SNAP is close to the experimental value (1.62% difference), and the total energy needed for fracturing is approximately the same. On the other hand in the ELPL model, although the maximum tensile principal stresses are limited to the tensile strength of concrete, the amount of force is raised with increasing the mid-span vertical displacement. This behavior is basically occurred based on the lack of appropriate mechanism for modeling the softening phenomenon in the elasto-plastic models. Meanwhile, the predicted initial stiffness of beam by SNAP is slightly lower than the experimental value. Compared to the Malvar and FourneyÕs model, the present analysis gives better results for NOSC model in this problem. Additional verifying examples can be
found elsewhere [10] which are not discussed here due to space limitations.
5. Analysis of Shahid Rajaee arch dam Consider the idealized symmetric model of Shahid Rajaee concrete arch dam as shown in Fig. 4. The dam is 130 m high with a crest length of 420 m. It is located in the north of Iran, in the seismically active foothills of Alborz Mountains near the city of Sari. The damÕs seismic response is analyzed by using a linear and three non-linear models described in the following. Let LIN designates the linear analysis case, used mainly for comparative purposes, NOSC-O and NOSCU the non-linear analysis cases using the non-orthogonal smeared crack model with ordinary and unlimited amounts of fracture energy, respectively. Meanwhile, ELPL is the non-linear analysis case corresponding to the elasto-plastic model. In all these cases, the same finite element discretization is utilized, consisting of 487 nodes and 76 isoparametric 20-noded solid elements in two layers through the thickness of the dam. The earthquake excitations include two components of the Friuli–Tolmezzo earthquake (the cross-canyon component is neglected to preserve symmetric conditions) whose records are normalized based on the frequency content for maximum design earthquake (MDE) condition with the peak ground acceleration of 0.56 g (Fig. 5). The Rayleigh damping coefficients were determined such that equivalent damping for frequencies close to the first and sixth modes of vibration would be 10% of the critical damping. The time integration is performed for 12 s, which covers the main response of the structure. As mentioned in Section 2, the foundation is assumed rigid to keep the computation time realistic. Although this assumption would highly influence the boundary stresses in the linear case, it is believed to be less important for the non-linear cases. In addition, the
1466
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
Fig. 4. Idealized model of Shahid Rajaee concrete arch dam.
Fig. 5. Accelerograms of applied earthquake (scale factor for MDE ¼ 9:81).
assumption of rigid foundation has less effect on the stresses in the vicinity of the crest spillway, where the major failure is expected to occur. As discussed earlier, the dam–water interaction effects are included by the modified WestergaardÕs method. Although different concrete mixes were used in the dam construction, the material properties are assumed uniform. For this study, the basic material properties are: E ¼ 30 GPa, q ¼ 2400 kg/m3 and m ¼ 0:18. Meanwhile, for the NOSC cases the concrete tensile strength is assumed as ft ¼ 3:0 MPa, which is more representative of tensile strength of mass concrete of dam body. The fracture energy release rate is assumed as Gf ¼ 2400 N/m and Gf ¼ 2400 106 N/m for the NOSC-O and NOSC-U cases, respectively. It should be mentioned that high discrepancy has been reported in the amount of fracture energy of concrete (Gf ) based on type of aggregates, cements and additives. In this regard, RILEM [19] shows a variation in the value of Gf from 70 to 200 N/m. Meanwhile, The results of a study performed on the concrete of existing dams in Switzerland show the
amount of Gf for dam concrete is two to three times higher than ordinary concrete [20]. Furthermore, in practical studies, researchers have often utilized relatively higher values of Gf . For instance, Vargas-Loli and Fenves [21] used Gf ¼ 2200 N/m for analysis of Pine Flat dam and Feltrin et al. [22] assumed Gf ¼ 1375 N/m. It must be also mentioned that it is not possible to lower the amount of fracture energy release rate considerably, due to the applied mesh sizes. Otherwise, the amount of crack tangent operator will be increased and if the value exceeds the elastic modulus of material, correct softening behavior is not achievable. Therefore, to study the effect of this parameter on the behavior of the dam, it is decided to raise this value significantly instead. The threshold angle between consecutive cracks in the NOSC cases is assumed to be equal to 15° and shear retention factors for all of the smeared crack models are considered constant and equal to b ¼ 0:1. A bilinear softening curve with constant parameters, a1 ¼ 0:01 and a2 ¼ 1 for the first branch, and a2 ¼ 0:0001 for the second branch is used. In addition, the cohesion and angle of
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
internal friction of concrete are determined based on the uniaxial compressive and tensile strengths for the ELPL case. These parameters are c ¼ 4:74 MPa and / ¼ 54:9°. The analysis commenced by applying the static loads. These include the dead weight of the dam body and the hydrostatic water pressures corresponding to 122 m height of reservoir water. These loads were incrementally increased in time until they reached their full values. Three time steps were initially considered, i.e. 0.01, 0.005 and 0.0025 s, and it was noticed that the results were almost the same for Dt ¼ 0:005 and 0.0025 s. Therefore, the former value was selected for time integration [10]. In this respect, the dead load is applied in 20 increments and the hydrostatic pressures thereafter in 20 increments at negative range of time. At time zero, the actual non-linear dynamic analysis is started with static displacements and stresses being applied as initial conditions.
1467
Table 1 Maximum displacements at dam crest (mm) Case
Right quarter point
Center point
U
V
W
V
W
LIN NOSC-O NOSC-U ELPL
)28.5 )26.2 )27.6 )27.3
64.2 64.0 64.2 63.5
5.4 7.5 4.9 5.0
117.1 )120.9 110.4 114.0
)13.1 64.0 17.3 17.0
5.1. Effect of fracture energy As mentioned earlier, it was not possible to decrease the amount of fracture energy significantly, based on the selected mesh size and amount of tensile strength of concrete. Therefore, to study the influence of fracture energy on the behavior the dam, it was decided to increase the amount of this parameter. Meanwhile, a fictitious case for smeared crack approach can be obtained by selecting a very large value for Gf . The basic physical characteristic of such a definition is making a constant maximum tensile principal stress, which is equivalent to the tensile strength of material, for all values of crack normal strain. This arrangement is almost similar to the elasto-perfect plastic behavior in tensile region and elastic behavior on compressive region. In this regard, the elastic unloading option is also utilized for the NOSC cases. To study the displacements of the dam, two nodes are selected on the dam crest. These are located at the center and at approximately the right quarter point. The maximum displacements through time for each of the x, y, and z directions at these nodes are summarized in Table 1, i.e. U , V , and W , respectively. The analysis did not reveal large displacements in the structure under the severe earthquake excitations considered, giving an initial indication of a safe design of the dam based on the adopted theories and assumptions. In Figs. 6 and 7, displacement histories of mid-crest point for the NOSC cases are compared with the linear case, in the stream and vertical directions, respectively. Comparison among the results of the smeared crack models and those of the linear model shows a drift in the crest displacements in upward and upstream directions. This is the main characteristic of displacement histories in smeared crack approach, which is due to extensive cracking of downstream face of the dam beneath the
Fig. 6. Comparison of displacements in the stream direction at mid-crest point among the linear and smeared crack models.
Fig. 7. Comparison of displacements in the vertical direction at mid-crest point among the linear and smeared crack models.
spillway block. When most of these cracks are in open state, they produce a cumulative crack strain, which is the major source of the drift. However, the results of the NOSC-U case are settled between the two other cases and the amounts of drift in stream and vertical directions are significantly less than the NOSC-O case. Stiffer behavior in the NOSC-U case in comparison with the NOSC-O case is the major cause of such response. Since
1468
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
in the NOSC-U case the tensile strength of the material is almost constant for all of the values of crack normal strain, the load bearing capacity of the material is almost the same during the analysis. Therefore, in this case the material is stiffer than the NOSC-O case and it shows lower tendency toward opening of cracks and increasing the crack strain. The envelope of maximum principal stresses throughout the analysis is displayed in Figs. 8–15, while a summary of the results is obtained and tabulated in Table 2. In the linear analysis, very high tensile stresses are calculated at the base of the dam (Fig. 8), due to the rigid foundation assumption. In addition, high tensile stresses occur in the spillway region, which in reality are released with the opening of contraction joints. A maximum compressive stress of )15.25 MPa is obtained in the upstream face near the spillway region (Fig. 12). The NOSC model can completely restrict the amount of tensile stresses. In this regard, no overstressing in tensile stresses is observed in the NOSC-O case, which means that tensile stresses are limited to the tensile strength of concrete. However, slight overstressing of tensile stresses (less than 6%) are still noticed in the dam body for the NOSC-U case, due to threshold angle condition for initiation of subsequent crack(s). As mentioned earlier, a new crack is initiated at a sampling point when not only the maximum principal stress exceeds the tensile strength but also the angle between the new crack and the existing cracks is more than the threshold angle.
Comparison among the principal (tensile and compressive) stress results shows that the general patterns of contours are very similar in both above-mentioned NOSC cases. The extent of softening obtained from different regions of the smeared crack models is shown in plots of reduced tensile strength at the end of analysis (Figs. 16– 18). These figures together with the crack patterns (Figs. 19 and 20), which actually show the normal vector to each crack plane, illustrate the extent of damage occurring in the dam body. It is noticed that the extension of cracking in both NOSC cases, i.e. NOSC-O and NOSC-U, is almost the same. However, the cracked zone in the NOSC-U case is more extended in upstream face of spillway region than the NOSC-O case, conversely for the base region of the dam, in the downstream face the extension of damaged zone is higher for the NOSC-O case. In addition, the minimum tensile strength of material is not changed through the analysis for the NOSC-U case (Fig. 17), due to its special form of softening curve. In this case, for all of the values of crack strain, the amount of crack normal stress is constant and equal to the initial tensile strength of concrete.
5.2. Elasto-plastic model In this section, the results of elasto-plastic model (ELPL case) is studied and compared with the linear
Fig. 8. Envelope of maximum tensile stresses (MPa) for the LIN.
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
1469
Fig. 9. Envelope of maximum tensile stresses (MPa) for the NOSC-O.
Fig. 10. Envelope of maximum tensile stresses (MPa) for the NOSC-U.
model (LIN case) and a fictitious state of NOSC model (NOSC-U case). With this arrangement, it is possible to study the elasto-plastic behavior of the material and compare it with linear elastic and a special form of smeared crack approach.
As mentioned earlier, the central and the right quarter points of the dam crest are selected to monitor the displacement history of the dam. The maximum displacements through time at these nodes are summarized in Table 1.
1470
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
Fig. 11. Envelope of maximum tensile stresses (MPa) for the ELPL.
Fig. 12. Envelope of maximum compressive stresses (MPa) for the LIN (U/S view).
Fig. 13. Envelope of maximum compressive stresses (MPa) for the NOSC-O (U/S view).
It is observed that large displacements are not occurring in the structure under the severe earthquake in the ELPL case similar to the other cases. In Figs. 21 and
22, displacement histories of mid-crest point for the NOSC-U case are compared with the ELPL case, in the stream and vertical directions, respectively. Based on
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
1471
Fig. 14. Envelope of maximum compressive stresses (MPa) for the NOSC-U.
Fig. 15. Envelope of maximum compressive stresses (MPa) for the ELPL.
the results, the displacements and amount of drift in the ELPL case are very close to the fictitious state of smeared crack approach, i.e. NOSC-U. The envelope of maximum principal stresses throughout the analysis are obtained and summarized in
Table 2 while the principal stress envelopes are displayed in Figs. 11 and 15. The linear and fictitious state of the NOSC cases are discussed in the former section and the elasto-plastic model is compared with the mentioned models in this section.
1472
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
Table 2 Maximum tensile and compressive stresses (MPa) Case
LIN NOSC-O NOSC-U ELPL
Maximum tensile stress (r1 )
Maximum compressive stress (r3 )
Spillway and central region
Base and abutments
Spillway and central region
Base and abutments
U/S
D/S
U/S
D/S
U/S
D/S
U/S
D/S
8.93 3.00 3.17 3.06
8.93 3.00 3.17 3.06
15.76 3.00 3.17 3.06
2.09 3.00 1.60 1.53
)15.25 )15.54 )15.26 )15.46
)13.58 )13.85 )15.26 )15.46
)6.91 )15.54 )13.60 )13.77
)11.92 )12.17 )11.93 )12.09
Fig. 16. Reduced tensile strength (MPa) at the end of analysis for the NOSC-O.
Fig. 17. Reduced tensile strength (MPa) at the end of analysis for the NOSC-U.
Fig. 18. Effective plastic strain at the end of analysis for the ELPL.
Fig. 19. Crack pattern at the end of analysis for the NOSC-O.
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
1473
Fig. 20. Crack pattern at the end of analysis for the NOSC-U.
Fig. 21. Comparison of displacements in the stream direction at mid-crest point between fictitious state of smeared crack and elasto-plastic models.
stress distribution and maximum value of tensile stresses are significantly varied in comparison with linear case, as the amounts of tensile stresses are limited in the spillway and base regions to a maximum value close to tensile strength of the material. Meanwhile, the stress distribution is similar to the NOSC-U case with maximum difference of less than 4%. For the ELPL case, it is noticed that the distribution of maximum compressive stresses is similar to the linear analysis in some regions of the dam. Moreover, the stress results are very close to the NOSC-U case through the dam body. This demonstrates that performance of the elasto-plastic model in compressive regions is similar to the NOSC-U case, and state of stresses has not reached the plastic surface in the compressive region. Therefore for compressive stresses, the conduct of the dam in all analyses is very close to linear elastic behavior under the presumed earthquake excitations. The contours of effective plastic strain are illustrated in Fig. 18, to indicate the extension of the plastic zones through the dam body. As it is depicted in this figure, the highest values of effective plastic strain are produced in the upstream face of the dam base, and afterwards in the spillway region, of course in lower magnitudes. 6. Conclusions
Fig. 22. Comparison of displacements in the vertical direction at mid-crest point between fictitious state of smeared crack and elasto-plastic models.
Based on the obtained results, the ELPL model can appropriately restrict the amount of tensile stresses. In this regard, the maximum tensile stresses through dam body are limited to r1 ¼ 3:06 MPa, which is very close to the tensile strength of mass concrete. In addition,
A computer program is developed for the non-linear dynamic analysis of concrete arch dams based on the fixed smeared crack approach and elasto-plastic modeling of the material. The program is initially verified by analyzing a notched concrete beam and comparing the results with the corresponding experimental data. Subsequently, an idealized symmetric model of Shahid Rajaee arch dam is considered and four cases are analyzed: a linear case (LIN), a non-orthogonal smeared crack with ordinary fracture energy (NOSC-O), a nonorthogonal smeared crack with a very high value of fracture energy (NOSC-U), and an elasto-plastic model (ELPL). The main conclusions of the analyses are: A drift in upward and upstream directions forms the prominent characteristic of histories of the crest displacements in the smeared crack and elasto-plastic models.
1474
R. Espandar, V. Lotfi / Computers and Structures 81 (2003) 1461–1474
The non-orthogonal smeared crack and elasto-plastic models can overcome the overstressing problem encountered in the analysis by linear model and predict realistic distribution and levels of stresses in the dam body. The stress results for each time step are changed when the amount of fracture energy is significantly increased. Correspondingly, in the results of envelopes of the principal stresses, some differences between the amounts of stresses can be seen for the different zones of dam body. However, general shape of the contours of these envelopes, indicate no major change between ordinary and fictitious smeared crack approaches. In this regard, the displacement histories of the monitoring points on the dam crest show a considerable differences and results of the fictitious model are inclined toward the linear model. The comparison between the ordinary form of smeared crack and elasto-plastic models reveals major differences in the results. However, comparison between the fictitious form of smeared crack and elasto-plastic models shows very similar results. In this regard, both envelopes of principal stresses, as well as time history of displacements at dam crest is very much similar. This indicates that the major part of plastic behavior is occurring in the tensile region. It must be finally emphasized that knowing the deficiencies of usual elasto-plastic models in tensile region in comparison with fracture mechanics based models, it is recommended to use fracture based models for dynamic analysis of concrete arch dams.
Acknowledgement This research investigation was supported by grant no. 31305301 from the National Research Council of Iran to the Amirkabir University, Tehran, Iran. The support is gratefully acknowledged.
References [1] Kuo JSH. On the nonlinear dynamic response of arch dams to earthquakes––I. Fluid–structure interaction: Added-mass computations for incompressible fluid. II. Joint opening nonlinear mechanism: Interface smeared crack model. PhD thesis. University of California, Berkeley, 1982. [2] Dungar R. A viscoplastic bounding surface model for concrete under compressive and tensile conditions and its application to arch dam analysis. Proceedings of Second International Conference on Constitutive Laws for Engineering Materials. Tucson/ AZ. 1987. p. 849–56.
[3] Cervera M, Oliver J, Faria R. Seismic evaluation of concrete dams via continuum damage models. Earthquake Eng Struct Dyn 1995;24:1225–45. [4] Noruziaan B. Nonlinear seismic analysis of concrete arch dams. PhD thesis. Department of Civil and Environmental Engineering. Carleton University. Ottawa, Canada, 1995. [5] Hall JF. Efficient non-linear seismic analysis of arch dams. Earthquake Eng Struct Dyn 1998;27:1425–44. [6] Owen DRJ, Hinton E. Finite element in plasticity. Swansea: Pineridge Press Ltd; 1980. [7] de Borst R, Nauta P. Non-orthogonal cracks in a smeared finite element model. Eng Comput 1985;2(3):35–46. [8] Rots JG. Computational modelling of concrete fracture. Doctorate dissertation. TH Delft. The Netherlands, 1988. [9] Channakeshava C. Elasto-plastic-cracking analysis of plain and reinforced concrete structures. PhD dissertation. Indian Institute of Science, Bangalore, India, 1987. [10] Espandar R. A study on nonlinear dynamic response of concrete arch dams. PhD dissertation. Department of Civil Engineering, Amirkabir University of Technology, Tehran, Iran, 2001. [11] Bazant ZP, Oh BH. Crack band theory for fracture of concrete. Mater Struct RILEM 1983;16(93):155– 77. [12] Cervera M. Nonlinear analysis of reinforced concrete structures using three dimensional and shell finite element models. PhD dissertation. Department of Civil Engineering. University College of Swansea, Wales, UK, 1986. [13] Rots JG, Blaauwendraad J. Crack models for concrete–– Discrete or smeared? Fixed, multi-directional or rotating? Heron 1989;34(1). [14] Barzegar F, Maddipudi S. Three-dimensional modeling of concrete structures. J Struct Eng ASCE 1997;123(10):1339– 56. [15] Zienkiewicz OC, Taylor RL. The finite element method. 4th ed. Maidenhead, Berkshire: McGraw-Hill; 1989. [16] Espandar R. SNAP: Seismic nonlinear analysis program. Tehran, Iran: Amirkabir University of Technology; 2000. [17] Malvar LJ, Fourney ME. A three dimensional application of the smeared crack approach. Eng Fract Mech 1990;35(1/ 2/3):251–60. [18] Cornelissen HAW, Hordijk DA, Reinhardt HW. Experiments and theory for the application of fracture mechanics to normal and lightweight concrete. Proceedings of International Conference on Fracture Mechanics of Concrete, Lausanne, Switzerland, 1985. [19] RILEM, Technical Committee 90 FMA. Fracture mechanics of concrete structures. Chapman and Hall Ltd, 1989. [20] Bruhwiler E, Wittmann FH. Failure of dam concrete subjected to seismic load conditions. Eng Fract Mech 1990;35(1/2/3):565–71. [21] Vargas-Loli L, Fenves GL. Effects of concrete cracking on the earthquake response of gravity dams. Earthquake Eng Struct Dyn 1989;18:575–92. [22] Feltrin G, Galli M, Bachmann H. Influence of cracking on the earthquake response of concrete gravity dam with reservoir. Proceedings of the 10th World Conference on Earthquake Engineering, vol. 8, Madrid, Spain, 1992. p. 4627–32.