- Email: [email protected]
A discrete crack joint model for nonlinear dynamic analysis of concrete arch dam
Computers and Structures 79 (2001) 403±420
www.elsevier.com/locate/compstruc
A discrete crack joint model for nonlinear dynamic analysis of concrete arch dam M.T. Ahmadi a,*, M. Izadinia a, H. Bachmann b a
Department of Structural Engineering, School of Engineering, Tarbiat Modarres University, P.O. Box 14115-143, Tehran, Iran b IBK, Swiss Federal Institute of Technology (ETH), Hoenggerberg, Switzerland Received 30 November 1998; accepted 27 May 2000
Abstract Concrete arch dams are generally constructed of massive plain concrete with almost no tensile resistance. To control tensile forces due to concrete shrinkage, temperature variations and for construction facilitation, arch dams are built in cantilever monoliths separated by vertical contraction joints. Earlier studies show that the modeling of such joints has signi®cant in¯uence on the seismic safety evaluation of arch dams. This fact is due to the tensile and shear failures of joints causing a redistribution of internal forces during and after a big earthquake. In the present study, a nonlinear joint element model with a coupled shear±tensile behavior for realistic ®nite element analysis of dam±reservoir system is presented. Reservoir upstream radiation, and bottom partial absorption of acoustic waves, as well as water compressibility are considered. The model when applied to simpler cases solved by other workers shows good performances. However, it is much more useful to solve problems not considered so far, e.g., shear keys behavior, joint damages, etc. The model could be employed eectively and conveniently for earthquake safety evaluation of arch dams in highly active seismic regions. Ó 2000 Elsevier Science Ltd. All rights reserved. Keywords: Arch dam; Joint element; Contraction joint; Dynamic analysis; Finite elements; Shear failure; Discrete crack; Joint opening; Joint slip; Nonlinear analysis
1. Introduction Various methods have been proposed to model the discontinuous behavior of arch dam vertical joints. This phenomenon is treated by using constraint equations or by adopting special joint elements such as discrete springs [1]. However, such generalized nonlinear springs do not include all the mechanical properties of contraction joints, e.g., shear keys properties, initial tensile strength due to joint grouting, shear softening, etc. in numerical analysis. For example, shear and tensile failures of joints have an important coupling feature de-
* Corresponding author. Tel.: +98-911-230-0463; fax: +9821-800-6652. E-mail addresses: [email protected], [email protected] (M.T. Ahmadi).
pending on the amount of normal opening, shear key height and damages due to nonlinear deformations. Furthermore, considering the geometrical and physical features of joints, additional examinations will be needed to determine such springs stiness. Other workers have treated the joint/interface by using a quasicontinuum ®nite element of small thickness containing planes of weakness. A special joint ®nite element has been developed by Beer [2] and is successfully applied to rock joints studies and other geo-mechanical phenomena. Further studies are needed to extend such applications to dynamic analysis of jointed systems. Some de®ciencies have been observed in the nonlinear behavior, i.e., in both normal and tangential displacements of the joint model developed by Fenves [3]. He studied the nonlinear behavior of joints related only to the joint normal displacement. He assumed fully elastic behavior for joint response to tangential forces, without adequate
0045-7949/01/$ - see front matter Ó 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 0 ) 0 0 1 4 8 - 6
404
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
interaction with normal displacements (i.e., using shear keys with in®nite shear strength). Further shortcomings correspond to modeling the reservoir hydrodynamic forces by added mass approximate method. Hohberg [4] studied the nonlinear behavior of joints with relative normal and shear displacements extensively. He developed mathematical constitutive relations idealizing the behavior of concrete contraction joints. However, his model for joint softening behavior showed de®ciencies in convergence and diculties for parameters decision. Indeed it requires more descriptive examples, experiments and comparative studies with the observed damages to concrete contraction joints and shear keys. He also did not include the reservoir hydrodynamic forces as desired for large dams. Noruziaan [5] studied the nonlinear behavior due to joint failure and concrete cracking in the ®nite element analysis of ÔMorrow PointÕ concrete arch dam. He concluded that concrete nonlinear behavior have negligible eects on the damÕs overall response. This is doubtful in the case of strong earthquakes. He did not disclose the details of his time integration and numerical algorithm to deal with the proposed constitutive law of concrete joints. His proposed model should be tested by known simple examples of jointed systems. Besides, he also used the added mass concept for the reservoir effects. Discrete crack model for concrete joints taken by Lot® [6] did not describe all the physical features observed in concrete joints. He also applied the added mass concept for reservoir eects. Ahmadi provided a discrete crack ®nite element model for vertical and foundation joints of an arch dam [7]. He considered only tensile cracking due to persistent static loads and provided a lower bound solution for the failure load. He also studied the eect of dierent abutment ¯exibilities on joint opening due to weight and water loads. A recent release of program D I A N A [8] has the coupled shear± tensile failure modeling capability using CoulombÕs friction law, but apparently it does not consider the joint asperity height which controls the shear interlock. Recent developments have shown that in addition to reservoir hydrodynamic interaction, ¯uid compressibility could also considerably aect the earthquake responses of arch dams [9]. It is now well understood that the multivariable discrete nature of contraction joints of arch dam should be taken into consideration, and the present research deals with a realistic and novel handling of such a task. 2. Methodology 2.1. Constitutive relations for the adopted joint element Providing shear key members in arch dam upper contraction joints is the usual construction practice in
many countries. This approach is considered very effective for ensuring dam safety in seismic regions or sites with variable deformability in the foundation. With the present research introduced in this paper, eects of shear keys on earthquake response of the dam±reservoir system could be studied. Local discontinuities in solid displacement ®eld due to existing joints, cause the variational problem and minimization of total potential energy to be subjected to constraint equations for joint displacements. This constrained minimization problem could be analyzed by the Lagrange multipliers method or penalty function approach. However, even when assuming linear elastic solids and small deformations, this problem will remain strongly nonlinear because of the unknown state of boundary conditions in the interface [4]. Relative tangential and normal deformations at the joint surfaces cause internal resisting forces or stresses. In the elastic range, constitutive law for the interface surface is given by frg De fdg; 8 9 2 r kn0 0 > < n> = 6 ss 4 0 ks0 > : > ; 0 0 st
38 9 0 > = 7 0 5 us ; > : > ; ks0 ut
1
where rn , ss , st are the normal and the two tangential (shear) stresses, and v, us , ut are the normal and the tangential relative displacements. ks0 , kn0 are the penalty parameters and guarantee the no-slip and no-penetration conditions of the interface surfaces for service loads. These parameters have the dimension of force per unit volume and are also called the initial shear and normal joint stiness coecients. In Eq. (1), it is assumed that the coupling between normal and shear displacements is negligible, and the elastic modulus matrix De is diagonal. In a typical analysis of systems including joint elements, relative tangential and normal displacements in each quadrature point of joint surface are calculated by nodal displacements vector. Then, the resisting forces and stiness matrix can be determined by Eq. (1). Element geometry and local and global coordinates for the adopted 16-node joint element are presented in Fig. 1. It is possible to develop elasto-plastic modulus matrix in nonlinear stages of deformations by a yield function equation. This function in the three-dimensional (3-D) stress space ``rn , ss , st '' could be based on Mohr± Coulomb yield criteria and a bounding tensile strength value. It could be shown that the ¯ow rule for slip is now nonassociative. In this context, Hohberg [4] studied the softening behavior or the damage to dierent components of shear strength by multimechanism plasticity theory. However, actual or safe application of this model to dynamic analysis of arch dams requires extended studies. The elasto-plastic modulus matrix can
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
405
Fig. 1. Joint element with zero thickness: (a) element geometry, (b) element geometry, local and global coordinates, and (c) isoparametric coordinate system.
thus be written for isotropic joint element without softening behavior as 2 kn0 ks0 4 ÿlbs kn0 D ks0 lmkn0 ÿlb k ep
t n0
ÿmbs kn0 lmkn0 b2t ks0 ÿbs bt ks0
3 ÿmbt kn0 ÿbs bt ks0 5; lmkn0 b2s ks0 2
where l is the coecient of friction, m, a quantity equal or less than l in the equation of plastic potential function, bs cos c, bt sin c, where c arctan st =ss is the angle between the tangential stress components in the s and t perpendicular directions on the joint surface [4]. kn0 and ks0 correspond to undamaged joint properties. In Eq. (2), the interface stresses can be represented as functions of dierent variables in the following manner:
or
rn rn u; v; l; m; kn0 ; ks0 ;
3
frg Dfdg:
s s u; v; l; m; kn0 ; ks0 :
4
Path-dependent stiness functions ks , kn (later on called joint stiness coecients) are de®ned later. It is important to note that such presentation of a nonassociative elasto-plastic relationship does not overlook its coupled and nonsymmetric nature and is eectively serving similar to the original full matrix form of relationship. Thus, the nonsymmetric form of equations is avoided without losing consistency. Further, it is well understood that unlike concrete mass cracks, softening for a joint is quite abrupt, for both shear and tensile failure modes. After studying the contraction joints response of concrete arch dams to static and dynamic loads and the related experimental works done so far in this context, the following assumptions are considered for the joint constitutive law: 1. The joints have zero or small initial tensile strength. 2. Normal stiness coecient kn , and normal stress rn , will vanish after joint opening (or crack mode-I), but tangential stiness coecient ks will decrease in
Here, ÔuÕ stands for both us and ut as ÔsÕ stands for ss and st . Therefore, the form of matrix ÔDep Õ could now be modi®ed into an equivalent diagonalized form as the modulus matrix ÔDÕ still implying coupling between shear and normal forces adequately. Typical couplings between physical parameters rn , s, v, and u are shown in Fig. 2 for one-dimensional (1-D) joint laboratory observations [4]. Considering the nature of arch dam contraction joints and based on engineering intuition, the more in¯uential parameters including the joint initial tensile strength Ft , joint asperity height Dn , and joint cohesion c, are chosen and taken into consideration in a similar way as expressed above, and thus, rn kn u; v; Ft ; kn0 v;
5
s ks v; l; c; Dn ; rn ; ks0 u;
6
Fig. 2. Observed joint behavior in normal and tangential displacements where su is the ultimate shear strength, ss , the shear stress, m, the joint dilatancy and rn , the compressive stress.
7
406
3.
4.
5. 6.
7. 8.
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
steps with joint opening. At this stage, the joint undergoes elasto-plastic deformations. While the joint opening exceeds a speci®ed limit ÔDn Õ (characteristic asperity or shear key height), both the tangential stress and the tangential stiness will disappear. If shear stress for an intact joint is less than the shear strength (as provided by cohesion and friction according to CoulombÕs relation), then the joint will have elastic tangential displacements, otherwise after a shear failure (or crack mode-II) and loss of shear strength, joint will have a perfectly plastic state. The value of joint dilatation is assumed to be equal to zero. After the joint closes, friction will be activated against slippage. In this manner, the joint behaves elastically or in an elasto-plastic manner with coupling shear and normal displacements. The joint surfaces have isotropic properties. The joint shear failure will be followed by elimination of shear strength. However, upon re-loading, a residual shear strength as provided by a reduced friction coecient will arise.
9. Upon any type of joint failure, tensile strength will be lost, and shear stiness is reduced. The above elaborated joint model, namely the simpli®ed discrete crack joint (SDCJ) model for numerical computations is introduced in Fig. 3a±d with examples of dierent load paths and scenarios. In this ®gure, the following reduction factors (or softening parameters) are assigned to the joint after or during joint opening or joint slippage: r is the permanent reduction factor for friction coecient after shear failure, m, the permanent reduction factor for shear stiness after shear failure, and n, the temporary reduction factor for shear stiness after joint opening less than shear key height. These reduction factors or softening parameters along with other properties of the joint, such as tensile strength Ft , cohesion c, friction coecient l, asperity height Dn , have to be decided through experimental or design speci®cation examinations. With reference to assumptions 2, 3 and 6 for the joint constitutive law and Fig. 3d, the coupling between shear and normal displacements have been fully considered in static and dynamic nonlinear analyses.
Fig. 3. Relation between stresses and normal or tangential displacements. (a) Normal displacement: (1) Tensile fracture and (2) joint opening exceeding the shear key height. (b) Constitutive law for joint in shear: (0) loading, (1) shear fracture, (2) shear unloading, (3) new shear fracture, (4) new loading and (5) new shear fracture. (c) Subsequent yield surface for joint damage. (d) Shear and tensile fractures with coupling: (0) shear loading in tension, (1) reduction of shear stiness with partial joint opening, (2) shear unloading in compression, (3) complete joint opening, (4) shear loading in compression, (5) shear fracture, (6) shear loading, (7) complete slip in tension and (8) shear unloading.
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
The model, although seemingly simple to apply, is novel, new and practically similar to fuzzy approaches, which deals with quite complicated problems consistently. 2.2. Reservoir mathematical model The governing equation for ¯uid domain is the Helmholtz equation for hydrodynamic pressure: r2 p
1 p; C2
8
where p is the hydrodynamic pressure, and C, the acoustic wave velocity in water. The above equation implies small displacements of invicid compressible ¯uid with irrotational motion. Water compressibility has a signi®cant in¯uence on the ¯uid±structure interaction for a wide range of ratio of natural frequencies of structure to ¯uid domain, including the case of higher and stier dams [9]. Thus, for general applicability and completeness of the dam±reservoir formulation, one needs to include the reservoir water compressibility. Boundary conditions are expressed as op 1 p0 oz g
9
for the reservoir-free surface, op 1 ÿ p_ on bC
10
for the reservoir bottom partial absorption, op 1 ÿ p_ on C
11
for the reservoir upstream face radiation of acoustic waves, and op ÿqans on
12
for the interaction boundary between dam and reservoir. In the above equations, z is the vertical coordinate, b, the acoustic impedance ratio of rock to water, n, the vector perpendicular to the boundary, q, the mass density of water, g, the gravitational acceleration, and ans , the acceleration of dam upstream face in the normal direction. Here, we have assumed that the hydrodynamic waves satisfy the 1-D wave propagation equation (11), through the upstream reservoir near-®eld truncation surface. This boundary, sometimes known as the Sommerfeld or viscous boundary, performs well in time domain analysis when applied suciently far from the structure. It is applicable only if compressibility is included [10]. The above equations along with the governing equation for the structure would lead to a simultaneous
407
dierential equations set for the coupled dam±reservoir system. These equations are discretized by the ®nite element method in a standard way similar to that of Ref. [11]. To avoid prohibitively high number of nonsymmetric equations with large bandwidth, the staggering solution method [11] is employed. Here, the displacement and the pressure ®elds are solved alternatively in each time step to achieve ``inter-domain compatibility'' or convergence. 3. Computer implementations The SDCJ is established to study contraction joints in static and dynamic analysis of arch dams. In this model, two types of cracks (or failure modes) are de®ned at each Gauss point of joint elements, i.e., crack mode-I (or crack I) due to tensile failure, and crack mode-II (or crack II) due to shear failure. A state parameter (ITEN) describes the joint behavior in the following manner: ITEN 0 for intact joint, ITEN 1 for crack-I in compression, ITEN 2 for crack-I in tension with v P Dn , ITEN 3 for crack-I in tension with v < Dn , ITEN 4 for crack-II in compression with s P su (su , being the joint shear strength), ITEN 5 for crack-II in compression with s < su , and ®nally ITEN 6 for crack-II in tension. In each load step, the joint incremental stresses in local coordinate system will be calculated by the relation dr D dd, where dd is the relative displacement increment for the current iteration, and D is a variable consistently diagonalized matrix de®ned by Eqs. (5)±(7), as expressed in Fig. 4. It is similar to the elastic matrix De of Eq. (1), but the stiness coecients ks and kn of Eqs. (5) and (6) are determined by Table 1, regarding the current state parameter (ITEN) of the corresponding Gauss point. In this table, ks0 and kn0 are initial stiness coecients. According to the above descriptions, a ¯owchart summarizing the constitutive behavior of the joint is established as shown in Fig. 4 and implemented in a computer program. The joint quadrature order is proposed to be three. 4. Numerical results and discussion 4.1. Preliminary example The proposed constitutive law has been examined by some simple examples of jointed systems. 4.1.1. Cantilever beam with roller support Dynamic analysis of a cantilever beam with roller support is studied. A prismatic ¯exural beam with one ®xed and one moment resistant no-shear roller end
408
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Fig. 4. Flowchart of the SDCJ model: determination of new D or ITEN, based on current a, r and ITEN for each time step (or loading step), each iteration, in each element at each Gauss point (a: nodal displacement vector, r: stress vector of one Gauss point).
support is modeled by four 20-node brick elements and one joint element as visualized in Fig. 5. In order to represent correctly the roller support behavior by joint element, the following reasonable values are assumed. Joint initial stiness coecients ks0 0:0, kn0 9:81 1012 N/m3 , tensile strength Ft
9:81 1010 N/m2 , cohesion c 9:81 1010 N/m2 . By adopting these quantities, the joint element will behave like a roller but has an elastic normal deformation without separation or opening. This ¯exural beam is analyzed for a base acceleration of sine half-wave with an amplitude equal to 20 m/s2 in
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
roller support proving desirable similarity as shown in Fig. 6.
Table 1 Stiness coecients in matrix D Stiness coecient
ks kn
409
Joint state parameter (ITEN) Intact
Tensile damage
Shear damage
0
1
2
3
4
5
6
ks0 kn0
ks0 kn0
0 0
nks0 0
0 kn0
mks0 kn0
0 0
4.2. Earthquake response of Morrow Point arch dam
the transverse (vertical) direction for 0.2 s. Mass density of the beam material is equal to q 2400:0 kg/m3 and no damping is assumed. This problem is re-analyzed by modeling the ¯exural beam as a 1-D frame element in the y-direction with appropriate support nodes, and using classical theory available in current commercial ®nite element programs. Time history results of the two approaches have been drawn for displacement and bending moment at the
The response results presented in this section are for a well studied arch dam. It is the Morrow Point dam located on the Gunnison river in Colorado, USA. This dam is a 141.73 m high, approximately symmetric, single centered arch dam. Detailed description of the geometry of this dam is available in Ref. [12]. To limit the problem size, the dam and reservoir system is assumed to be symmetric about the y±z plane. Moreover, the dam-foundation interaction is neglected. The ground motion recorded at Taft Lincoln school during the Kern County, California earthquake of 21 July 1952 is selected as the free-®eld ground acceleration (Fig. 7).
Fig. 5. Cantilever beam, example 4.1.1.
Fig. 6. Results of dynamic analysis for cantilever beam.
410
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Fig. 7. Original ground motion of Taft Lincoln California Earthquake of 21 July 1952.
The ground motions acting in the stream (y), and the vertical (z) directions are de®ned as the S69E and the vertical components of the recorded ground motion respectively. The ®nite element idealization of the system have the following characteristics: modulus of elasticity for concrete E 2 1010 N/m2 , concrete PoissonÕs ratio mc 0:2, concrete unit mass qc 2400:0 kg/m3 , and internal viscous damping ratio 0.05. Ten 22-node brick elements and ten 16-node joint elements with a total number of nodal points equal to 187 (Fig. 8) are employed in the dam-body ®nite element model. Thus, three contraction joints are considered in a half-dam-body model. Reservoir domain includes 40 27-node elements, with a total number of pressure nodes equal to 495. Water level elevation for both hydrostatic and hydrodynamic pressure calculations is equal to the dam crest elevation (141.73 m). Acoustic wave velocity in water is C, 1440.0 m/s. The acoustic impedance ratio of rock to water is, b 3. After the literature review on contraction joints, properties of the joint elements adopted for contraction joint with shear keys, namely the ``original'' properties of joints, are as follows: joint initial shear stiness coecient ks0 1 109 N/m3 , initial normal stiness coecient kn0 2 109 N/m3 , coecient of friction l 0:9, cohesion c 1:5 106 N/m2 , tensile strength
Fig. 8. Dam±reservoir ®nite element model of Morrow Point arch dam.
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Ft 1 106 N/m2 , reduction factor for friction coef®cient due to shear failure r 0:9, reduction factor for shear stiness due to joint opening n 0:7, reduction factor for shear stiness due to shear failure m 0:2, shear key height of joint surface Dn 0:5 m. A time step Dt 0:005 s, with the full Newton± Raphson nonlinear solution algorithm is used. Therefore, the structural stiness matrix has to be updated in each iteration. The classical Newmark method with standard parameters (0.25 and 0.5) is adopted for the time integration of the system of equations of motion. An explicit Runge±Kutta stress increment computation is applied using half the value of strain increment to determine the D matrix from the ¯owchart of Fig. 4 in each iteration. This ensures the calculated total stresses not to fall beyond the margin of the yield surface. Convergence tolerance for nonlinear displacement iterations is based on the energy norms de®ned as jE ij=jE 1j 1:0 10ÿ12 . For pressure iterations in the staggering scheme, convergence is based on the pressure norm as jjDpjj=jjpjj 0:001. Maximum number of iterations for pressure is 8, and for displacement is 10. Integration order for all elements is 3 except for the dam thickness direction for which 2 is deemed as sucient. 4.2.1. Veri®cation under medium ground motion In order to verify the numerical results, Morrow Point dam±reservoir system was analyzed by program A N S Y S Ver.5.0 [13]. Its corresponding model consists of ten 20-node brick elements (SOLID95) for the dam, 59 2-node 1-D joint elements (CONTAC52) for contraction joints and 80 8-node elements (FLUID30) for the ¯uid
411
domain. None of the available versions of the above program (including recent ones) has an account on either cohesion or initial tensile strength for our purpose. Also, interaction between shear and tensile failures due to shear keys height could not be considered by dierent existing versions of this program. The joint elements CONTAC52 stiness values are obtained from the ks0 , and kn0 values described above, multiplied by the tributary area of each joint nodal point. Due to A N S Y S limitations, here the joint properties have to be modi®ed for both the programs as follows hereafter. This set of joint properties will be called ``modi®ed'' properties of joints. Tensile strength, cohesion and friction coecient of joints are assumed as Ft 0:0, c 0:0, and l 1 respectively, implying no shear failure except when joints open. Shear keys height is Dn 0:0, and reduction factors for joints failure are m r 1:0, and n 0 (no shear interlock after joint opening). The Taft ground motion records with the peak horizontal acceleration normalized to Ô0.2g' are taken as the free-®eld input excitation for the initial comparative studies of the two codes. Static analysis of the system under self-weight and the hydrostatic pressure was carried out to establish the initial condition for dynamic analysis. Sequentially, dynamic analysis of the system was performed by both authorsÕ and the A N S Y S programs. The history of displacement in the stream direction of nodal point 186 in the middle of dam crest and the history of hydrodynamic pressure at nodal point 31 in the mid-height of dam upstream face are shown in Figs. 9 and 10. The history of normal relative displacement v, normal stress
Fig. 9. History of displacement in the stream direction at nodal point 186 with maximum ground acceleration equal to 0.2g.
412
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Fig. 10. History of hydrodynamic pressure at nodal point 31 with maximum ground acceleration equal to 0.2g.
Fig. 11. History of normal relative displacement (or joint opening) for the Gauss point 6 of joint element 8 with maximum ground acceleration equal to 0.2g.
rn , and resultant tangential stress s, for the Gauss point 6 of joint element 8 below the crest, are shown in Figs. 11±13. For this level of acceleration, only three upper joint elements encountered shear or tensile failure.
Maximum opening for joint elements is still small and equal to 0.6 cm (at joint element 10 at time 8.405 s). The comparison denotes that the results of the two programs are close to each other. Most of the dierences
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
413
Fig. 12. History of normal stress for the Gauss point 6 of joint element 8 with maximum ground acceleration equal to 0.2g.
Fig. 13. History of resultant tangential stress for the Gauss point 6 of joint element 8 with maximum ground acceleration equal to 0.2g.
are due to local behavior of joints, i.e., joint-opening, joint normal and shear stresses. Such discrepancies are generally attributed to dierent ®nite element modeling characteristics as A N S Y S uses a 1-D one-to-one contact
elements rather than the 3-D surface contact elements used by the authors. The spurious hydrodynamic pressure obtained by A N S Y S in the ®rst seconds of motion is not appreciated (Fig. 10).
414
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Fig. 14. History of displacement in the stream direction at nodal point 186 with maximum ground acceleration equal to 0.5g.
Fig. 15. History of hydrodynamic pressure at nodal point 31 with maximum ground acceleration equal to 0.5g.
4.2.2. Veri®cation under strong ground motion Using the modi®ed joint properties again, but under 0.5g ground acceleration records, another comparative study is carried out between the authorsÕ algorithm and
that of A N S Y S . Its results are shown in Figs. 14±18. There is an overall agreement between the two programs. However, as seen in Figs. 17 and 18, according to the present model analysis, from about the time of 7.5 s,
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
415
Fig. 16. History of normal relative displacement (or joint opening) for the Gauss point 6 of joint element 8 with maximum ground acceleration equal to 0.5g.
Fig. 17. History of normal stress for the Gauss point 6 of joint element 8 with maximum ground acceleration equal to 0.5g.
joint element 8 undergoes complete tensile and shear release as a result of large shift and damage in the dam body while A N S Y S still shows small ¯uctuations of both
stresses beyond this time. This seems to be an outcome of the initial contact stress distribution at the interface which is rather dierent for the two models.
416
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Fig. 18. History of resultant tangential stress for the Gauss point 6 of joint element 8 with maximum ground acceleration equal to 0.5g.
Fig. 19. History of linear versus nonlinear displacements in the stream direction at nodal point 186 with maximum ground acceleration equal to 0.5g (according to the present model).
4.2.3. Nonlinear versus linear solution for original model under strong ground motion Further analysis is carried out with the original joint properties input data described above as the nonlinear
model, versus a fully linear (without joints) model. Other properties of dam and reservoir model are similar as before. This time, the peak horizontal acceleration of ground is again set to 0.5g to visualize the ability of the
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
417
Fig. 20. History of linear/nonlinear solution-based hydrodynamic pressures at nodal point 31 with maximum ground acceleration equal to 0.5g (according to the present model).
methodology to predict the in¯uence and damages due to very strong earthquakes. Figs. 19 and 20 illustrate the crest displacement responses and the corresponding hydrodynamic pressures on the upstream mid-height point of the dam in the cases of nonlinear (jointed) and the linear (monolithic) models of the dam body. This could well prove the signi®cance of inclusion of contraction joint failures in the dam response as the dam body is pronouncedly shifted upstream up to about 8 cm permanently. At the same time, maximum dynamic displacement is magni®ed by about 50% when joint failures are admitted. In this case, when the ground motion is terminated, all the joint elements other than elements 1 and 2 have encountered shear, or tensile, or a combination of both types of damages. Generally, joint opening happens for joint elements on the plane of symmetry while combined joint opening
and slippage happen for other joint elements. Maximum shear deformations (tangential displacements) for the joints which experienced shear failure according to the present nonlinear model, is shown in Table 2. It is evident that such values of shear deformations are signi®cant and could not be neglected during very strong ground motions in contrary to what was assumed in previous works [3]. Maximum dam crest displacement in the stream direction is equal to 34.4 cm for node 186, and maximum opening for joint elements is equal to 4.9 cm (joint element 10 on the plane of symmetry, at time 7.245 s) which is still less than the assumed shear key height. Maximum tensile and compressive principal stresses induced in the dam body are equal to 9.8 and 14.6 MPa, respectively. Maximum tensile principal stresses are more than the dynamic strength values of mass concrete, and thus the concrete mass material proves to behave in a nonlinear way under tension as a result of such a strong ground motion. It is interesting to note that compression is still well in the elastic range.
Table 2 Maximum shear displacement for joint elementsa
4.2.4. Joint properties parametric study Lastly, a parametric study is carried out to observe the sensitivity of dam responses to some joint properties under the latter load and with the adopted model. By changing a single parameter of joints (either joint initial tensile strength, Ft or joint cohesion, c) while other properties are constant, the stream-direction maximum
Joint element number Maximum shear displacement (cm) a
At Gauss points.
3
5
6
8
9
3.1
2.4
17.9
10.2
36.9
418
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Fig. 21. Sensitivity study of maxima of displacement and joint opening with respect to (a) joint tensile strength and (b) joint cohesion.
displacement of nodal point 186 (upper curve in both Fig. 21a and b) and the maximum opening of joint element 10 on the dam crest (lower curve in the same ®gures) are computed. It is interesting that for a wide
range of variation for both cohesion and joint tensile strength the general response is not much aected. However, the local stress distributions proved not to follow this insensitivity.
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
5. Conclusions and future works After reviewing previous researches on contraction joints behavior in dynamic response analysis of arch dams, a discrete crack joint model namely SDCJ, is proposed based on the most signi®cant mechanical features of arch dam joints. Full dam±reservoir interaction with upstream radiation and reservoir bottom partial absorption of pressure waves along with surface gravity wave boundaries are satis®ed. Water compressibility which has a major role in certain realistic cases of such systems is also included. It is well understood from the results presented in this paper that vertical contraction joints play an important role in response analysis of concrete arch dams in case of very strong earthquakes; a conclusion in agreement with other researchers ®ndings. Eciency of the proposed joint model is promising for earthquake safety evaluation of concrete arch dams. It is a powerful model which takes into account all the major factors and parameters affecting the joints of arch dams such as joint initial tensile strength, joint cohesion, joint friction coecient, shear key height, joint shear softening after damages, and most important, the coupling between joint opening and joint slippage similar to nonassociated plasticity theory. Besides, the model has physically meaningful parameters which are presumably more easier to determine. It also does not suer convergence and solution de®ciencies related to complicated mathematical nonassociated plasticity models proposed so far. Although the shear resistance is based on Mohr±Coulomb relationship, but the model is versatile enough to accommodate other laws as well. Finally, it is easy to implement in standard codes. By applying the new model to a large existing arch dam, it is observed that under very strong ground motion, vertical joints generally encounter shear and/or tensile failure, and cause redistribution of internal forces. For earthquakes in the stream or vertical direction and under symmetric loading, joints on the midplane generally have normal opening or tensile failure, and joints on the quarter section of arch span show shear failure. Joints opening magnitude is generally about a few centimeters and does not seem to exceed regular shear keys height at geometrically similar sites. However, sucient strength for shear keys is needed, otherwise their failure would endanger the integrity and stability of arch dams against severe earthquake. Careful implementation of appropriate shear keys for vertical joints causes increased strength and stability. It is claimed that shear keys supply the cohesion of joints eciently. They also provide an interlock shear-carrying system to help integrity of the jointed structure after joint opening. Unlike concrete compressive stresses, tensile stresses generally enter into the nonlinear range for high intensity earthquakes with the peak acceleration
419
of about 0.5g, and nonlinear material laws for concrete fracture might be necessary to consider. In order to increase the accuracy of computations, more precise boundary conditions for the reservoir boundaries and an inclusion of the dam±foundation interaction is desirable. Additional study together with laboratory experiments is necessary for determining the physical properties of contraction joints. The values of cohesion, friction coecient, and joint softening parameters should be determined regarding properties of concrete material and shear keys geometry. An extensive parametric study on the in¯uence of shear keys is the subject of future researches of the authors using the powerful methodology introduced here. Acknowledgements Parts of this research was made possible while the ®rst author had his sabbatical in the IBK, Swiss Federal Institute of Technology, ETH-Zurich. The supports of International Institute of Earthquake Engineering and Seismology, and Moshanir Consultant Engineers, Tehran, Iran is also greatly appreciated. References [1] Dowling MJ, Hall JF. Nonlinear seismic analysis of arch dams. J Engng Mech ASCE 1989;115(4):768±89. [2] Beer G. An isoparametric joint/interface element for ®nite element analysis. Int J Num Meth Engng 1985;21:585±600. [3] Fenves GL, Mojtahedi S, Reimer RB. Eect of contraction joints on earthquake response of an arch dam. J Struct Engng ASCE 1992;118(4):1039±55. [4] Hohberg JM. Seismic arch dam analysis with full joint nonlinearity. Proc Int Conf Dam Fracture, Denver, Colorado, 1991. p. 61±75. [5] Noruziaan B. Nonlinear seismic analysis of concrete arch dams, PhD Thesis, Department of Civil and Environmental Engineering, Carleton University, Ottawa, Canada, 1995. [6] Lot® V. Comparison of discrete crack and elasto-plastic models in nonlinear dynamic analysis of arch dams. Dam Engng 1996;VII(1):65±110. [7] Ahmadi MT, Razavi S. A three-dimensional joint opening analysis of an arch dam. Comput Struct 1992;44(1/2):187± 92. [8] TNO Building and Construction, Displacement Analyser D I A N A ver 7.2, January 2000, Netherlands. [9] Fok KL, Chopra AK. Water compressibility in earthquake response of arch dams. J Struct Engng ASCE 1987;113(5): 958±75. [10] Zienkiewicz OC, et al. The Sommerfeld radiation conditions on in®nite domains and its modelling in numerical procedures, IRIA Third Int Symp Comput Meth Appl Sci Engng, 1977. [11] Zienkiewicz OC, Chan AHC. Coupled problems and their numerical solution. In: Doltsinis IS, editor. Advances in
420
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
computational nonlinear mechanics, Springer, Berlin, 1988. p. 109±76 [chapter 3]. [12] US Department of the Interior, Bureau of Reclamation, Morrow Point Dam and Powerplant, Technical Record of
Design and Construction, US Government Printing Oce, Denver, Colorado, 1983. [13] Kohnke P. A N S Y S UserÕs Manual for Revision 5.0, Swanson Analysis Systems Inc., Houston, PA, USA, 1992.
www.elsevier.com/locate/compstruc
A discrete crack joint model for nonlinear dynamic analysis of concrete arch dam M.T. Ahmadi a,*, M. Izadinia a, H. Bachmann b a
Department of Structural Engineering, School of Engineering, Tarbiat Modarres University, P.O. Box 14115-143, Tehran, Iran b IBK, Swiss Federal Institute of Technology (ETH), Hoenggerberg, Switzerland Received 30 November 1998; accepted 27 May 2000
Abstract Concrete arch dams are generally constructed of massive plain concrete with almost no tensile resistance. To control tensile forces due to concrete shrinkage, temperature variations and for construction facilitation, arch dams are built in cantilever monoliths separated by vertical contraction joints. Earlier studies show that the modeling of such joints has signi®cant in¯uence on the seismic safety evaluation of arch dams. This fact is due to the tensile and shear failures of joints causing a redistribution of internal forces during and after a big earthquake. In the present study, a nonlinear joint element model with a coupled shear±tensile behavior for realistic ®nite element analysis of dam±reservoir system is presented. Reservoir upstream radiation, and bottom partial absorption of acoustic waves, as well as water compressibility are considered. The model when applied to simpler cases solved by other workers shows good performances. However, it is much more useful to solve problems not considered so far, e.g., shear keys behavior, joint damages, etc. The model could be employed eectively and conveniently for earthquake safety evaluation of arch dams in highly active seismic regions. Ó 2000 Elsevier Science Ltd. All rights reserved. Keywords: Arch dam; Joint element; Contraction joint; Dynamic analysis; Finite elements; Shear failure; Discrete crack; Joint opening; Joint slip; Nonlinear analysis
1. Introduction Various methods have been proposed to model the discontinuous behavior of arch dam vertical joints. This phenomenon is treated by using constraint equations or by adopting special joint elements such as discrete springs [1]. However, such generalized nonlinear springs do not include all the mechanical properties of contraction joints, e.g., shear keys properties, initial tensile strength due to joint grouting, shear softening, etc. in numerical analysis. For example, shear and tensile failures of joints have an important coupling feature de-
* Corresponding author. Tel.: +98-911-230-0463; fax: +9821-800-6652. E-mail addresses: [email protected], [email protected] (M.T. Ahmadi).
pending on the amount of normal opening, shear key height and damages due to nonlinear deformations. Furthermore, considering the geometrical and physical features of joints, additional examinations will be needed to determine such springs stiness. Other workers have treated the joint/interface by using a quasicontinuum ®nite element of small thickness containing planes of weakness. A special joint ®nite element has been developed by Beer [2] and is successfully applied to rock joints studies and other geo-mechanical phenomena. Further studies are needed to extend such applications to dynamic analysis of jointed systems. Some de®ciencies have been observed in the nonlinear behavior, i.e., in both normal and tangential displacements of the joint model developed by Fenves [3]. He studied the nonlinear behavior of joints related only to the joint normal displacement. He assumed fully elastic behavior for joint response to tangential forces, without adequate
0045-7949/01/$ - see front matter Ó 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 0 ) 0 0 1 4 8 - 6
404
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
interaction with normal displacements (i.e., using shear keys with in®nite shear strength). Further shortcomings correspond to modeling the reservoir hydrodynamic forces by added mass approximate method. Hohberg [4] studied the nonlinear behavior of joints with relative normal and shear displacements extensively. He developed mathematical constitutive relations idealizing the behavior of concrete contraction joints. However, his model for joint softening behavior showed de®ciencies in convergence and diculties for parameters decision. Indeed it requires more descriptive examples, experiments and comparative studies with the observed damages to concrete contraction joints and shear keys. He also did not include the reservoir hydrodynamic forces as desired for large dams. Noruziaan [5] studied the nonlinear behavior due to joint failure and concrete cracking in the ®nite element analysis of ÔMorrow PointÕ concrete arch dam. He concluded that concrete nonlinear behavior have negligible eects on the damÕs overall response. This is doubtful in the case of strong earthquakes. He did not disclose the details of his time integration and numerical algorithm to deal with the proposed constitutive law of concrete joints. His proposed model should be tested by known simple examples of jointed systems. Besides, he also used the added mass concept for the reservoir effects. Discrete crack model for concrete joints taken by Lot® [6] did not describe all the physical features observed in concrete joints. He also applied the added mass concept for reservoir eects. Ahmadi provided a discrete crack ®nite element model for vertical and foundation joints of an arch dam [7]. He considered only tensile cracking due to persistent static loads and provided a lower bound solution for the failure load. He also studied the eect of dierent abutment ¯exibilities on joint opening due to weight and water loads. A recent release of program D I A N A [8] has the coupled shear± tensile failure modeling capability using CoulombÕs friction law, but apparently it does not consider the joint asperity height which controls the shear interlock. Recent developments have shown that in addition to reservoir hydrodynamic interaction, ¯uid compressibility could also considerably aect the earthquake responses of arch dams [9]. It is now well understood that the multivariable discrete nature of contraction joints of arch dam should be taken into consideration, and the present research deals with a realistic and novel handling of such a task. 2. Methodology 2.1. Constitutive relations for the adopted joint element Providing shear key members in arch dam upper contraction joints is the usual construction practice in
many countries. This approach is considered very effective for ensuring dam safety in seismic regions or sites with variable deformability in the foundation. With the present research introduced in this paper, eects of shear keys on earthquake response of the dam±reservoir system could be studied. Local discontinuities in solid displacement ®eld due to existing joints, cause the variational problem and minimization of total potential energy to be subjected to constraint equations for joint displacements. This constrained minimization problem could be analyzed by the Lagrange multipliers method or penalty function approach. However, even when assuming linear elastic solids and small deformations, this problem will remain strongly nonlinear because of the unknown state of boundary conditions in the interface [4]. Relative tangential and normal deformations at the joint surfaces cause internal resisting forces or stresses. In the elastic range, constitutive law for the interface surface is given by frg De fdg; 8 9 2 r kn0 0 > < n> = 6 ss 4 0 ks0 > : > ; 0 0 st
38 9 0 >
1
where rn , ss , st are the normal and the two tangential (shear) stresses, and v, us , ut are the normal and the tangential relative displacements. ks0 , kn0 are the penalty parameters and guarantee the no-slip and no-penetration conditions of the interface surfaces for service loads. These parameters have the dimension of force per unit volume and are also called the initial shear and normal joint stiness coecients. In Eq. (1), it is assumed that the coupling between normal and shear displacements is negligible, and the elastic modulus matrix De is diagonal. In a typical analysis of systems including joint elements, relative tangential and normal displacements in each quadrature point of joint surface are calculated by nodal displacements vector. Then, the resisting forces and stiness matrix can be determined by Eq. (1). Element geometry and local and global coordinates for the adopted 16-node joint element are presented in Fig. 1. It is possible to develop elasto-plastic modulus matrix in nonlinear stages of deformations by a yield function equation. This function in the three-dimensional (3-D) stress space ``rn , ss , st '' could be based on Mohr± Coulomb yield criteria and a bounding tensile strength value. It could be shown that the ¯ow rule for slip is now nonassociative. In this context, Hohberg [4] studied the softening behavior or the damage to dierent components of shear strength by multimechanism plasticity theory. However, actual or safe application of this model to dynamic analysis of arch dams requires extended studies. The elasto-plastic modulus matrix can
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
405
Fig. 1. Joint element with zero thickness: (a) element geometry, (b) element geometry, local and global coordinates, and (c) isoparametric coordinate system.
thus be written for isotropic joint element without softening behavior as 2 kn0 ks0 4 ÿlbs kn0 D ks0 lmkn0 ÿlb k ep
t n0
ÿmbs kn0 lmkn0 b2t ks0 ÿbs bt ks0
3 ÿmbt kn0 ÿbs bt ks0 5; lmkn0 b2s ks0 2
where l is the coecient of friction, m, a quantity equal or less than l in the equation of plastic potential function, bs cos c, bt sin c, where c arctan st =ss is the angle between the tangential stress components in the s and t perpendicular directions on the joint surface [4]. kn0 and ks0 correspond to undamaged joint properties. In Eq. (2), the interface stresses can be represented as functions of dierent variables in the following manner:
or
rn rn u; v; l; m; kn0 ; ks0 ;
3
frg Dfdg:
s s u; v; l; m; kn0 ; ks0 :
4
Path-dependent stiness functions ks , kn (later on called joint stiness coecients) are de®ned later. It is important to note that such presentation of a nonassociative elasto-plastic relationship does not overlook its coupled and nonsymmetric nature and is eectively serving similar to the original full matrix form of relationship. Thus, the nonsymmetric form of equations is avoided without losing consistency. Further, it is well understood that unlike concrete mass cracks, softening for a joint is quite abrupt, for both shear and tensile failure modes. After studying the contraction joints response of concrete arch dams to static and dynamic loads and the related experimental works done so far in this context, the following assumptions are considered for the joint constitutive law: 1. The joints have zero or small initial tensile strength. 2. Normal stiness coecient kn , and normal stress rn , will vanish after joint opening (or crack mode-I), but tangential stiness coecient ks will decrease in
Here, ÔuÕ stands for both us and ut as ÔsÕ stands for ss and st . Therefore, the form of matrix ÔDep Õ could now be modi®ed into an equivalent diagonalized form as the modulus matrix ÔDÕ still implying coupling between shear and normal forces adequately. Typical couplings between physical parameters rn , s, v, and u are shown in Fig. 2 for one-dimensional (1-D) joint laboratory observations [4]. Considering the nature of arch dam contraction joints and based on engineering intuition, the more in¯uential parameters including the joint initial tensile strength Ft , joint asperity height Dn , and joint cohesion c, are chosen and taken into consideration in a similar way as expressed above, and thus, rn kn u; v; Ft ; kn0 v;
5
s ks v; l; c; Dn ; rn ; ks0 u;
6
Fig. 2. Observed joint behavior in normal and tangential displacements where su is the ultimate shear strength, ss , the shear stress, m, the joint dilatancy and rn , the compressive stress.
7
406
3.
4.
5. 6.
7. 8.
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
steps with joint opening. At this stage, the joint undergoes elasto-plastic deformations. While the joint opening exceeds a speci®ed limit ÔDn Õ (characteristic asperity or shear key height), both the tangential stress and the tangential stiness will disappear. If shear stress for an intact joint is less than the shear strength (as provided by cohesion and friction according to CoulombÕs relation), then the joint will have elastic tangential displacements, otherwise after a shear failure (or crack mode-II) and loss of shear strength, joint will have a perfectly plastic state. The value of joint dilatation is assumed to be equal to zero. After the joint closes, friction will be activated against slippage. In this manner, the joint behaves elastically or in an elasto-plastic manner with coupling shear and normal displacements. The joint surfaces have isotropic properties. The joint shear failure will be followed by elimination of shear strength. However, upon re-loading, a residual shear strength as provided by a reduced friction coecient will arise.
9. Upon any type of joint failure, tensile strength will be lost, and shear stiness is reduced. The above elaborated joint model, namely the simpli®ed discrete crack joint (SDCJ) model for numerical computations is introduced in Fig. 3a±d with examples of dierent load paths and scenarios. In this ®gure, the following reduction factors (or softening parameters) are assigned to the joint after or during joint opening or joint slippage: r is the permanent reduction factor for friction coecient after shear failure, m, the permanent reduction factor for shear stiness after shear failure, and n, the temporary reduction factor for shear stiness after joint opening less than shear key height. These reduction factors or softening parameters along with other properties of the joint, such as tensile strength Ft , cohesion c, friction coecient l, asperity height Dn , have to be decided through experimental or design speci®cation examinations. With reference to assumptions 2, 3 and 6 for the joint constitutive law and Fig. 3d, the coupling between shear and normal displacements have been fully considered in static and dynamic nonlinear analyses.
Fig. 3. Relation between stresses and normal or tangential displacements. (a) Normal displacement: (1) Tensile fracture and (2) joint opening exceeding the shear key height. (b) Constitutive law for joint in shear: (0) loading, (1) shear fracture, (2) shear unloading, (3) new shear fracture, (4) new loading and (5) new shear fracture. (c) Subsequent yield surface for joint damage. (d) Shear and tensile fractures with coupling: (0) shear loading in tension, (1) reduction of shear stiness with partial joint opening, (2) shear unloading in compression, (3) complete joint opening, (4) shear loading in compression, (5) shear fracture, (6) shear loading, (7) complete slip in tension and (8) shear unloading.
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
The model, although seemingly simple to apply, is novel, new and practically similar to fuzzy approaches, which deals with quite complicated problems consistently. 2.2. Reservoir mathematical model The governing equation for ¯uid domain is the Helmholtz equation for hydrodynamic pressure: r2 p
1 p; C2
8
where p is the hydrodynamic pressure, and C, the acoustic wave velocity in water. The above equation implies small displacements of invicid compressible ¯uid with irrotational motion. Water compressibility has a signi®cant in¯uence on the ¯uid±structure interaction for a wide range of ratio of natural frequencies of structure to ¯uid domain, including the case of higher and stier dams [9]. Thus, for general applicability and completeness of the dam±reservoir formulation, one needs to include the reservoir water compressibility. Boundary conditions are expressed as op 1 p0 oz g
9
for the reservoir-free surface, op 1 ÿ p_ on bC
10
for the reservoir bottom partial absorption, op 1 ÿ p_ on C
11
for the reservoir upstream face radiation of acoustic waves, and op ÿqans on
12
for the interaction boundary between dam and reservoir. In the above equations, z is the vertical coordinate, b, the acoustic impedance ratio of rock to water, n, the vector perpendicular to the boundary, q, the mass density of water, g, the gravitational acceleration, and ans , the acceleration of dam upstream face in the normal direction. Here, we have assumed that the hydrodynamic waves satisfy the 1-D wave propagation equation (11), through the upstream reservoir near-®eld truncation surface. This boundary, sometimes known as the Sommerfeld or viscous boundary, performs well in time domain analysis when applied suciently far from the structure. It is applicable only if compressibility is included [10]. The above equations along with the governing equation for the structure would lead to a simultaneous
407
dierential equations set for the coupled dam±reservoir system. These equations are discretized by the ®nite element method in a standard way similar to that of Ref. [11]. To avoid prohibitively high number of nonsymmetric equations with large bandwidth, the staggering solution method [11] is employed. Here, the displacement and the pressure ®elds are solved alternatively in each time step to achieve ``inter-domain compatibility'' or convergence. 3. Computer implementations The SDCJ is established to study contraction joints in static and dynamic analysis of arch dams. In this model, two types of cracks (or failure modes) are de®ned at each Gauss point of joint elements, i.e., crack mode-I (or crack I) due to tensile failure, and crack mode-II (or crack II) due to shear failure. A state parameter (ITEN) describes the joint behavior in the following manner: ITEN 0 for intact joint, ITEN 1 for crack-I in compression, ITEN 2 for crack-I in tension with v P Dn , ITEN 3 for crack-I in tension with v < Dn , ITEN 4 for crack-II in compression with s P su (su , being the joint shear strength), ITEN 5 for crack-II in compression with s < su , and ®nally ITEN 6 for crack-II in tension. In each load step, the joint incremental stresses in local coordinate system will be calculated by the relation dr D dd, where dd is the relative displacement increment for the current iteration, and D is a variable consistently diagonalized matrix de®ned by Eqs. (5)±(7), as expressed in Fig. 4. It is similar to the elastic matrix De of Eq. (1), but the stiness coecients ks and kn of Eqs. (5) and (6) are determined by Table 1, regarding the current state parameter (ITEN) of the corresponding Gauss point. In this table, ks0 and kn0 are initial stiness coecients. According to the above descriptions, a ¯owchart summarizing the constitutive behavior of the joint is established as shown in Fig. 4 and implemented in a computer program. The joint quadrature order is proposed to be three. 4. Numerical results and discussion 4.1. Preliminary example The proposed constitutive law has been examined by some simple examples of jointed systems. 4.1.1. Cantilever beam with roller support Dynamic analysis of a cantilever beam with roller support is studied. A prismatic ¯exural beam with one ®xed and one moment resistant no-shear roller end
408
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Fig. 4. Flowchart of the SDCJ model: determination of new D or ITEN, based on current a, r and ITEN for each time step (or loading step), each iteration, in each element at each Gauss point (a: nodal displacement vector, r: stress vector of one Gauss point).
support is modeled by four 20-node brick elements and one joint element as visualized in Fig. 5. In order to represent correctly the roller support behavior by joint element, the following reasonable values are assumed. Joint initial stiness coecients ks0 0:0, kn0 9:81 1012 N/m3 , tensile strength Ft
9:81 1010 N/m2 , cohesion c 9:81 1010 N/m2 . By adopting these quantities, the joint element will behave like a roller but has an elastic normal deformation without separation or opening. This ¯exural beam is analyzed for a base acceleration of sine half-wave with an amplitude equal to 20 m/s2 in
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
roller support proving desirable similarity as shown in Fig. 6.
Table 1 Stiness coecients in matrix D Stiness coecient
ks kn
409
Joint state parameter (ITEN) Intact
Tensile damage
Shear damage
0
1
2
3
4
5
6
ks0 kn0
ks0 kn0
0 0
nks0 0
0 kn0
mks0 kn0
0 0
4.2. Earthquake response of Morrow Point arch dam
the transverse (vertical) direction for 0.2 s. Mass density of the beam material is equal to q 2400:0 kg/m3 and no damping is assumed. This problem is re-analyzed by modeling the ¯exural beam as a 1-D frame element in the y-direction with appropriate support nodes, and using classical theory available in current commercial ®nite element programs. Time history results of the two approaches have been drawn for displacement and bending moment at the
The response results presented in this section are for a well studied arch dam. It is the Morrow Point dam located on the Gunnison river in Colorado, USA. This dam is a 141.73 m high, approximately symmetric, single centered arch dam. Detailed description of the geometry of this dam is available in Ref. [12]. To limit the problem size, the dam and reservoir system is assumed to be symmetric about the y±z plane. Moreover, the dam-foundation interaction is neglected. The ground motion recorded at Taft Lincoln school during the Kern County, California earthquake of 21 July 1952 is selected as the free-®eld ground acceleration (Fig. 7).
Fig. 5. Cantilever beam, example 4.1.1.
Fig. 6. Results of dynamic analysis for cantilever beam.
410
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Fig. 7. Original ground motion of Taft Lincoln California Earthquake of 21 July 1952.
The ground motions acting in the stream (y), and the vertical (z) directions are de®ned as the S69E and the vertical components of the recorded ground motion respectively. The ®nite element idealization of the system have the following characteristics: modulus of elasticity for concrete E 2 1010 N/m2 , concrete PoissonÕs ratio mc 0:2, concrete unit mass qc 2400:0 kg/m3 , and internal viscous damping ratio 0.05. Ten 22-node brick elements and ten 16-node joint elements with a total number of nodal points equal to 187 (Fig. 8) are employed in the dam-body ®nite element model. Thus, three contraction joints are considered in a half-dam-body model. Reservoir domain includes 40 27-node elements, with a total number of pressure nodes equal to 495. Water level elevation for both hydrostatic and hydrodynamic pressure calculations is equal to the dam crest elevation (141.73 m). Acoustic wave velocity in water is C, 1440.0 m/s. The acoustic impedance ratio of rock to water is, b 3. After the literature review on contraction joints, properties of the joint elements adopted for contraction joint with shear keys, namely the ``original'' properties of joints, are as follows: joint initial shear stiness coecient ks0 1 109 N/m3 , initial normal stiness coecient kn0 2 109 N/m3 , coecient of friction l 0:9, cohesion c 1:5 106 N/m2 , tensile strength
Fig. 8. Dam±reservoir ®nite element model of Morrow Point arch dam.
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Ft 1 106 N/m2 , reduction factor for friction coef®cient due to shear failure r 0:9, reduction factor for shear stiness due to joint opening n 0:7, reduction factor for shear stiness due to shear failure m 0:2, shear key height of joint surface Dn 0:5 m. A time step Dt 0:005 s, with the full Newton± Raphson nonlinear solution algorithm is used. Therefore, the structural stiness matrix has to be updated in each iteration. The classical Newmark method with standard parameters (0.25 and 0.5) is adopted for the time integration of the system of equations of motion. An explicit Runge±Kutta stress increment computation is applied using half the value of strain increment to determine the D matrix from the ¯owchart of Fig. 4 in each iteration. This ensures the calculated total stresses not to fall beyond the margin of the yield surface. Convergence tolerance for nonlinear displacement iterations is based on the energy norms de®ned as jE ij=jE 1j 1:0 10ÿ12 . For pressure iterations in the staggering scheme, convergence is based on the pressure norm as jjDpjj=jjpjj 0:001. Maximum number of iterations for pressure is 8, and for displacement is 10. Integration order for all elements is 3 except for the dam thickness direction for which 2 is deemed as sucient. 4.2.1. Veri®cation under medium ground motion In order to verify the numerical results, Morrow Point dam±reservoir system was analyzed by program A N S Y S Ver.5.0 [13]. Its corresponding model consists of ten 20-node brick elements (SOLID95) for the dam, 59 2-node 1-D joint elements (CONTAC52) for contraction joints and 80 8-node elements (FLUID30) for the ¯uid
411
domain. None of the available versions of the above program (including recent ones) has an account on either cohesion or initial tensile strength for our purpose. Also, interaction between shear and tensile failures due to shear keys height could not be considered by dierent existing versions of this program. The joint elements CONTAC52 stiness values are obtained from the ks0 , and kn0 values described above, multiplied by the tributary area of each joint nodal point. Due to A N S Y S limitations, here the joint properties have to be modi®ed for both the programs as follows hereafter. This set of joint properties will be called ``modi®ed'' properties of joints. Tensile strength, cohesion and friction coecient of joints are assumed as Ft 0:0, c 0:0, and l 1 respectively, implying no shear failure except when joints open. Shear keys height is Dn 0:0, and reduction factors for joints failure are m r 1:0, and n 0 (no shear interlock after joint opening). The Taft ground motion records with the peak horizontal acceleration normalized to Ô0.2g' are taken as the free-®eld input excitation for the initial comparative studies of the two codes. Static analysis of the system under self-weight and the hydrostatic pressure was carried out to establish the initial condition for dynamic analysis. Sequentially, dynamic analysis of the system was performed by both authorsÕ and the A N S Y S programs. The history of displacement in the stream direction of nodal point 186 in the middle of dam crest and the history of hydrodynamic pressure at nodal point 31 in the mid-height of dam upstream face are shown in Figs. 9 and 10. The history of normal relative displacement v, normal stress
Fig. 9. History of displacement in the stream direction at nodal point 186 with maximum ground acceleration equal to 0.2g.
412
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Fig. 10. History of hydrodynamic pressure at nodal point 31 with maximum ground acceleration equal to 0.2g.
Fig. 11. History of normal relative displacement (or joint opening) for the Gauss point 6 of joint element 8 with maximum ground acceleration equal to 0.2g.
rn , and resultant tangential stress s, for the Gauss point 6 of joint element 8 below the crest, are shown in Figs. 11±13. For this level of acceleration, only three upper joint elements encountered shear or tensile failure.
Maximum opening for joint elements is still small and equal to 0.6 cm (at joint element 10 at time 8.405 s). The comparison denotes that the results of the two programs are close to each other. Most of the dierences
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
413
Fig. 12. History of normal stress for the Gauss point 6 of joint element 8 with maximum ground acceleration equal to 0.2g.
Fig. 13. History of resultant tangential stress for the Gauss point 6 of joint element 8 with maximum ground acceleration equal to 0.2g.
are due to local behavior of joints, i.e., joint-opening, joint normal and shear stresses. Such discrepancies are generally attributed to dierent ®nite element modeling characteristics as A N S Y S uses a 1-D one-to-one contact
elements rather than the 3-D surface contact elements used by the authors. The spurious hydrodynamic pressure obtained by A N S Y S in the ®rst seconds of motion is not appreciated (Fig. 10).
414
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Fig. 14. History of displacement in the stream direction at nodal point 186 with maximum ground acceleration equal to 0.5g.
Fig. 15. History of hydrodynamic pressure at nodal point 31 with maximum ground acceleration equal to 0.5g.
4.2.2. Veri®cation under strong ground motion Using the modi®ed joint properties again, but under 0.5g ground acceleration records, another comparative study is carried out between the authorsÕ algorithm and
that of A N S Y S . Its results are shown in Figs. 14±18. There is an overall agreement between the two programs. However, as seen in Figs. 17 and 18, according to the present model analysis, from about the time of 7.5 s,
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
415
Fig. 16. History of normal relative displacement (or joint opening) for the Gauss point 6 of joint element 8 with maximum ground acceleration equal to 0.5g.
Fig. 17. History of normal stress for the Gauss point 6 of joint element 8 with maximum ground acceleration equal to 0.5g.
joint element 8 undergoes complete tensile and shear release as a result of large shift and damage in the dam body while A N S Y S still shows small ¯uctuations of both
stresses beyond this time. This seems to be an outcome of the initial contact stress distribution at the interface which is rather dierent for the two models.
416
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Fig. 18. History of resultant tangential stress for the Gauss point 6 of joint element 8 with maximum ground acceleration equal to 0.5g.
Fig. 19. History of linear versus nonlinear displacements in the stream direction at nodal point 186 with maximum ground acceleration equal to 0.5g (according to the present model).
4.2.3. Nonlinear versus linear solution for original model under strong ground motion Further analysis is carried out with the original joint properties input data described above as the nonlinear
model, versus a fully linear (without joints) model. Other properties of dam and reservoir model are similar as before. This time, the peak horizontal acceleration of ground is again set to 0.5g to visualize the ability of the
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
417
Fig. 20. History of linear/nonlinear solution-based hydrodynamic pressures at nodal point 31 with maximum ground acceleration equal to 0.5g (according to the present model).
methodology to predict the in¯uence and damages due to very strong earthquakes. Figs. 19 and 20 illustrate the crest displacement responses and the corresponding hydrodynamic pressures on the upstream mid-height point of the dam in the cases of nonlinear (jointed) and the linear (monolithic) models of the dam body. This could well prove the signi®cance of inclusion of contraction joint failures in the dam response as the dam body is pronouncedly shifted upstream up to about 8 cm permanently. At the same time, maximum dynamic displacement is magni®ed by about 50% when joint failures are admitted. In this case, when the ground motion is terminated, all the joint elements other than elements 1 and 2 have encountered shear, or tensile, or a combination of both types of damages. Generally, joint opening happens for joint elements on the plane of symmetry while combined joint opening
and slippage happen for other joint elements. Maximum shear deformations (tangential displacements) for the joints which experienced shear failure according to the present nonlinear model, is shown in Table 2. It is evident that such values of shear deformations are signi®cant and could not be neglected during very strong ground motions in contrary to what was assumed in previous works [3]. Maximum dam crest displacement in the stream direction is equal to 34.4 cm for node 186, and maximum opening for joint elements is equal to 4.9 cm (joint element 10 on the plane of symmetry, at time 7.245 s) which is still less than the assumed shear key height. Maximum tensile and compressive principal stresses induced in the dam body are equal to 9.8 and 14.6 MPa, respectively. Maximum tensile principal stresses are more than the dynamic strength values of mass concrete, and thus the concrete mass material proves to behave in a nonlinear way under tension as a result of such a strong ground motion. It is interesting to note that compression is still well in the elastic range.
Table 2 Maximum shear displacement for joint elementsa
4.2.4. Joint properties parametric study Lastly, a parametric study is carried out to observe the sensitivity of dam responses to some joint properties under the latter load and with the adopted model. By changing a single parameter of joints (either joint initial tensile strength, Ft or joint cohesion, c) while other properties are constant, the stream-direction maximum
Joint element number Maximum shear displacement (cm) a
At Gauss points.
3
5
6
8
9
3.1
2.4
17.9
10.2
36.9
418
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
Fig. 21. Sensitivity study of maxima of displacement and joint opening with respect to (a) joint tensile strength and (b) joint cohesion.
displacement of nodal point 186 (upper curve in both Fig. 21a and b) and the maximum opening of joint element 10 on the dam crest (lower curve in the same ®gures) are computed. It is interesting that for a wide
range of variation for both cohesion and joint tensile strength the general response is not much aected. However, the local stress distributions proved not to follow this insensitivity.
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
5. Conclusions and future works After reviewing previous researches on contraction joints behavior in dynamic response analysis of arch dams, a discrete crack joint model namely SDCJ, is proposed based on the most signi®cant mechanical features of arch dam joints. Full dam±reservoir interaction with upstream radiation and reservoir bottom partial absorption of pressure waves along with surface gravity wave boundaries are satis®ed. Water compressibility which has a major role in certain realistic cases of such systems is also included. It is well understood from the results presented in this paper that vertical contraction joints play an important role in response analysis of concrete arch dams in case of very strong earthquakes; a conclusion in agreement with other researchers ®ndings. Eciency of the proposed joint model is promising for earthquake safety evaluation of concrete arch dams. It is a powerful model which takes into account all the major factors and parameters affecting the joints of arch dams such as joint initial tensile strength, joint cohesion, joint friction coecient, shear key height, joint shear softening after damages, and most important, the coupling between joint opening and joint slippage similar to nonassociated plasticity theory. Besides, the model has physically meaningful parameters which are presumably more easier to determine. It also does not suer convergence and solution de®ciencies related to complicated mathematical nonassociated plasticity models proposed so far. Although the shear resistance is based on Mohr±Coulomb relationship, but the model is versatile enough to accommodate other laws as well. Finally, it is easy to implement in standard codes. By applying the new model to a large existing arch dam, it is observed that under very strong ground motion, vertical joints generally encounter shear and/or tensile failure, and cause redistribution of internal forces. For earthquakes in the stream or vertical direction and under symmetric loading, joints on the midplane generally have normal opening or tensile failure, and joints on the quarter section of arch span show shear failure. Joints opening magnitude is generally about a few centimeters and does not seem to exceed regular shear keys height at geometrically similar sites. However, sucient strength for shear keys is needed, otherwise their failure would endanger the integrity and stability of arch dams against severe earthquake. Careful implementation of appropriate shear keys for vertical joints causes increased strength and stability. It is claimed that shear keys supply the cohesion of joints eciently. They also provide an interlock shear-carrying system to help integrity of the jointed structure after joint opening. Unlike concrete compressive stresses, tensile stresses generally enter into the nonlinear range for high intensity earthquakes with the peak acceleration
419
of about 0.5g, and nonlinear material laws for concrete fracture might be necessary to consider. In order to increase the accuracy of computations, more precise boundary conditions for the reservoir boundaries and an inclusion of the dam±foundation interaction is desirable. Additional study together with laboratory experiments is necessary for determining the physical properties of contraction joints. The values of cohesion, friction coecient, and joint softening parameters should be determined regarding properties of concrete material and shear keys geometry. An extensive parametric study on the in¯uence of shear keys is the subject of future researches of the authors using the powerful methodology introduced here. Acknowledgements Parts of this research was made possible while the ®rst author had his sabbatical in the IBK, Swiss Federal Institute of Technology, ETH-Zurich. The supports of International Institute of Earthquake Engineering and Seismology, and Moshanir Consultant Engineers, Tehran, Iran is also greatly appreciated. References [1] Dowling MJ, Hall JF. Nonlinear seismic analysis of arch dams. J Engng Mech ASCE 1989;115(4):768±89. [2] Beer G. An isoparametric joint/interface element for ®nite element analysis. Int J Num Meth Engng 1985;21:585±600. [3] Fenves GL, Mojtahedi S, Reimer RB. Eect of contraction joints on earthquake response of an arch dam. J Struct Engng ASCE 1992;118(4):1039±55. [4] Hohberg JM. Seismic arch dam analysis with full joint nonlinearity. Proc Int Conf Dam Fracture, Denver, Colorado, 1991. p. 61±75. [5] Noruziaan B. Nonlinear seismic analysis of concrete arch dams, PhD Thesis, Department of Civil and Environmental Engineering, Carleton University, Ottawa, Canada, 1995. [6] Lot® V. Comparison of discrete crack and elasto-plastic models in nonlinear dynamic analysis of arch dams. Dam Engng 1996;VII(1):65±110. [7] Ahmadi MT, Razavi S. A three-dimensional joint opening analysis of an arch dam. Comput Struct 1992;44(1/2):187± 92. [8] TNO Building and Construction, Displacement Analyser D I A N A ver 7.2, January 2000, Netherlands. [9] Fok KL, Chopra AK. Water compressibility in earthquake response of arch dams. J Struct Engng ASCE 1987;113(5): 958±75. [10] Zienkiewicz OC, et al. The Sommerfeld radiation conditions on in®nite domains and its modelling in numerical procedures, IRIA Third Int Symp Comput Meth Appl Sci Engng, 1977. [11] Zienkiewicz OC, Chan AHC. Coupled problems and their numerical solution. In: Doltsinis IS, editor. Advances in
420
M.T. Ahmadi et al. / Computers and Structures 79 (2001) 403±420
computational nonlinear mechanics, Springer, Berlin, 1988. p. 109±76 [chapter 3]. [12] US Department of the Interior, Bureau of Reclamation, Morrow Point Dam and Powerplant, Technical Record of
Design and Construction, US Government Printing Oce, Denver, Colorado, 1983. [13] Kohnke P. A N S Y S UserÕs Manual for Revision 5.0, Swanson Analysis Systems Inc., Houston, PA, USA, 1992.