Seismic Performance Evaluation of Concrete Gravity Dams with Penetrated Cracks Considering Fluid– Structure Interaction

In this paper, a comprehensive study on the seismic behavior of fractured concrete gravity dams during ground shakings is carried out considering dam–reservoir interaction effects. To gain the seismic behavior of the whole system, finite and boundary elements are employed to model the liquid region and the cracked structure, respectively. Formulation and different computational aspects of the suggested staggered hybrid approach are thoroughly argued. A computer code was developed in order to discuss the presented hybrid BE–DE technique and comparisons are made between the obtained results and those reported in the literature. To gain this goal, several problems of seismic excitations in frequencyand timedomains are presented employing the proposed approach, showing that the present results agree well with the results from other numerical procedures. The cracked Koyna Dam is scrutinized, considering the dynamic interaction between dam and reservoir with focus on the nonlinear behavior due to its top profile crack. The developed numerical model is rigorously validated by extensive comparisons with available results in the literature in which the dam–reservoir interaction were simplified by added masses. It can be concluded that there is significant disparity between the overturning and sliding response schemes of the nonlinear analysis and those of added mass technique. Keywords-seismic behavior; concrete gravity dams; boundary elements; distinct elements; dam–reservoir interaction


INTRODUCTION
Concrete gravity dams behave differently from other structures because of their size and their interactions with the reservoir and foundation.In fact, concrete gravity dams might have well cracked at their foundation or at a specific height caused by seasonal temperature variations, concrete shrinkage or earlier earthquakes.Apart from those cracks with confined depth that exist on both the upstream and downstream faces and which will possibly develop subjected to static or dynamic conditions, some cracks may have expanded so much that they almost penetrate the monoliths.Consequently, these existing cracks weaken the stability of concrete gravity dams principally due to the nonlinear behavior of sliding and overturning, during strong earthquakes.Therefore, nonlinear response of cracked concrete gravity dams has been investigated by different techniques to study the propagation of cracks in this type of dams.A number of researches focused on the propagation of cracks in the dam, something that includes opening and closing of the cracks [1][2][3].
Several researches have been conducted based on experiments, focusing mainly on the process of crack initiation and propagation [4][5][6].On the other hand, little attempts have been made for the case in which cracks are expected to penetrate from upstream to downstream in non-overflow monoliths of the dam.In comparison with various researches concerning crack event and propagation in concrete gravity dams, little efforts have been done to demonstrate the seismic behavior of cracked concrete gravity dams.The overturning analysis of the Koyna Dam was first investigated by authors in [8] assuming no sliding of the top block, which was regarded as a rigid body, ignoring the interaction between the top and bottom blocks and adopting a horizontally penetrated crack.Their results indicated that under these assumptions the top block remained stable and overturning of the top profile would not happen under future earthquake of similar magnitude.However, the top block is comparatively large, and the interaction between the two blocks may considerably affect their behavior.Authors in [9] employed contact elements located at the dam-foundation interface in order to ascertain the seismic sliding and uplifting response of concrete gravity dam monoliths.Authors in [10] scrutinized the frictional base sliding of concrete gravity dams subjected to seismic ground motion, in which the whole dam was modeled as a singledegree-of-freedom system.Presented results indicate that the dam flexibility affects significantly the sliding displacement in the downstream direction.Based on the single-degree-offreedom system criterion, authors in [11] developed an empirical formula for estimating the seismic induced slip of concrete gravity dams at the rock interface with finite-element analysis.Employing finite element method (FEM), authors in [12] calculated the base sliding of the gravity dams during ground excitation.Results of their parametric study of typical gravity dams revealed how the earthquake-induced sliding is affected by the characteristics of the ground motion and dam system.By using FEM, authors in [13] investigated experimentally the friction model of concrete lift joint interfaces and applied the model to assess the static and dynamic behavior of concrete gravity dams.They simulated the contact conditions with gapfriction elements that are appropriate for the case of small displacement and concluded that the sliding displacement is reduced in multi-crack cases.Using the continuous mass of liquid as lumped masses in the horizontal direction for the reservoir effects, authors in [14] investigated some important phenomena such as the jounce of the upper block of the fractured Koyna dam in detail, by using the distinct element method (DEM).Afterwards, authors in [15] developed a rigid model with three-degrees-of-freedom considering sliding, rocking, rock-sliding, and drifting modes.They concluded that a large coefficient of friction does not surely prevent sliding, and rocking and drifting modes should be mainly considered in estimating the stability of cracked concrete gravity dams.Employing rigid block models, authors in [16] suggested a convenient and simple tool to supply a conservative estimate of upper cracked block residual sliding displacements.Authors in [17] inspected the stability of a gravity dam on jointed rock foundation and also the seismic stability of the upper part of the fractured Konya dam using discontinuous deformation analysis (DDA).Authors in [1] examined the seismic response of Pine Flat Dam with penetrated cracks using the finite element modeling technique but limited in scope to small deformation of the continuum.
Hydrodynamic interaction in Dam-reservoir system is a significant aspect in seismic response of concrete gravity dams that may be taken into account by various numerical methods [18,19].A number of previous researches on the dynamic behavior of concrete gravity dams with penetrated cracks were performed considering the simplified added mass technique for reservoir effects [14,16,20].In fact this technique, applying the results of rigid dam-incompressible fluid interaction, leads to approximate results in the seismic response of fractured concrete gravity dams.More authentic results may be achieved through the use of numerical methods, among which FEM and BEM (boundary element modeling) are most popular.To particularly perform the nonlinear dam analysis, due to the needs for iterative solutions and very short consecutive time steps, minor degrees of freedom for reservoir domain is a more proper way.Besides, in the FEM modeling of reservoir, a radiating condition (for instance, Sommerfeld condition) should be employed to the far-field truncated boundary of reservoir to consider radiation damping.Thanks to substantial characteristic of BEM that is more felicitous for boundless domains, the technique is able to consider the radiating damping with much less additional attempt compared with the FEM analysis of reservoir.
Consequently, the dynamic analysis of dam-reservoir system is achieved by employing the BE equations governing the reservoir region coupled with the dam FE equations of motion as a straightforward approach for observation of the whole system behavior [21,22].Using a hybrid DE-BE method, authors in [23] modeled the seismic behavior of the cracked Koyna dam focusing on the upper profile crack.
With the significant improvements of discontinuous media methods like DEM, it becomes not very difficult to investigate the stability of fractured dams in destructive earthquakes.In this paper, DEM is employed to assess the seismic response of the Koyna dam with penetrated cracks subjected to earthquakes.This analysis method allows relative motion between elements, and is specifically appropriate for problems in which the relative motion between blocks contributes considerably to the general structure deformation.As a result, reservoir effects are accurately modeled by the use of BEM, and then a coupled BEM-DEM algorithm is developed to analyze the seismic behavior of cracked gravity dams examining dam-reservoir interaction effect.

II. HYBRID FORMULATION
A coupled numerical method which combines the distinct elements for the fractured concrete gravity dam and the boundary elements for the reservoir domain is developed in the time domain.A sketch of dam-reservoir system is exposed in Figure 1.Cracked dam scheme on rigid foundation.
The following assumptions are considered in the formulation of governing equations:  Seismic response analyses are carried out using the horizontal acceleration component.
 The dynamic dam-foundation interaction is ignored.
 The dam-reservoir system is analyzed under twodimensional (2D) condition.
 The liquid in the reservoir is inviscid and compressible resulting in an irrotational flow with negligible amplitudes.
 The nonlinear analysis of the dam is employed for the existing cracked zone and the linear elastic analysis is adopted for the deformable blοcks.
 Uplift pressures are disregarded in the cracked zone.

A. BEM Formulation for the Reservoir Domain
A short explanation of the reservoir BEM follows.Emphasis is principally devoted to the important aspects of the method, subjected to some modifications compared to previous works [24,25].These modifications may be remarked as:

www.etasr.com Behshad and Shekari: Seismic Performance Evaluation of Concrete Gravity Dams with Penetrated …
 At the far away from the upstream face of the dam, the domain is left opened because the radiation boundary condition is automatically satisfied by employing the Hankel fundamental solution as the most exact Green's function of boundless domains.
 Dynamic boundary conditions of the free surface gravitational waves and the bottom partial absorption are considered.
For a compressible, irrotational and inviscid fluid, the governing equation of liquid motion is expressed by the wellknown Helmholtz differential equation ( 1): in which,  is the velocity potential, and c denotes the acoustic wave velocity in fluid.Considering harmonic ground motion, the velocity potential in the reservoir can be described in the frequency domain as φ(t) = φ(ω)exp(iωt), where ω represents the exciting frequency and φ(ω) the complex-valued frequency response function for the velocity potential.The boundary element method is used to solve the Helmholtz equation with the following boundary conditions (see Figure1): a) Dynamic and kinematic boundary condition of the reservoir free surface (S4): b) The boundary condition at the fluid-solid interface (S1): where n is the unit normal vector, a n is the normal acceleration on the interface and ρ is the mass density of the fluid.
c) The boundary condition at the foundation-reservoir interface (S2) for the reservoir bottom partial absorption: in which, β=(ρ b c b /ρc) is the acoustic impedance ratio of the foundation to reservoir water, ρ b and ρ represent the mass densities of foundation material and water respectively, and c b represents the compression wave velocity of the foundation.d) And eventually, for the reservoir up stream radiation condition (S3) of acoustic waves, no additional discretization is required as elucidated above.
The linearized Bernoulli equation for pressure represents the hydrodynamic pressure acting on the dam system.
Employing the constant boundary element formulation, the following equations may be derived to solve the motion of liquid in discretized form [25].
 Hφ Gq (6) where H and G stand for boundary element coefficient matrices for the potential vector φ and the flux vector q, respectively.
The components of H and G coefficient matrices can be obtained as followings: ( ) in which s i is the boundary of the ith constant element, and φ i * and q i * denote the values of the fundamental solution φ * and q * at the ith constant element respectively.The fundamental solution φ * and its normal derivative q * are: in which r is the distance between the points on the boundary, i.e. source and field, H 0 (1) (kr) and H 1 (kr) show the Hankel functions of the first kind of order zero and of order one respectively, and n represents the outward normal to the boundary.Discretized form of boundary element equations of the reservoir domain may be written and separated into the following expressions, according to three types of assumed boundary conditions: where the subscripts w, i and b indicate the boundary nodes on the reservoir free surface, dam-reservoir interface and reservoir bottom surfaces, respectively.Substituting (2) -( 4) into (9) the following expression is obtained: Inspection of (10) reveals that φ could be determined, since all sub-matrix components of D are known.

B. Distinct Element Method
The distinct (or discrete) element method (DEM) has been exhaustively illustrated in many papers.This method was primarily proposed in [20] for the investigation of jointed and cracked rock masses.In this section, a brief review of some substantial aspects, especially those that would be subject to modification for their usage to the coupled dam-reservoir system, are reviewed.The openings are considered as interface boundary conditions between distinct blocks, and put into the problem domain according to their locations and orientations.The domain in DEM is divided into polygonal blocks by joints and the blocks are supposed to be quite deformable.Deformable blocks are discretized into a mesh of finite difference triangular elements in the DEM program [26] for stress-strain analysis, whose formulation is similar to the socalled constant strain triangle (CST) finite element.Failed contacts will be effaced and new contacts will be recognized and established when the block contact relationship varies.Normal and shear springs and dashpots are placed between contact elements for simulating the interaction behavior of contact points of block elements (see Figure 2).Spring-dashpot contact model.
The normal spring is linear-elastic in compression.Zero tension in the normal direction of contacts is presumed in the normal springs for modeling penetrated cracks.When normal compressive forces attain zero, joints open and blocks separate at that point.As a result, the analogous contact will be deleted.The shear spring is specified as linear elastic-perfectly plastic and assumed to follow Coulomb's friction law.The stiffness damping forces regarding the normal and shear dashpots are linearly related to the corresponding normal and shear velocities of the blocks at contact.The linear coefficients may be supposed to be identical for both normal and shear dashpots.The increments of normal and tangential spring forces in each time step (Δt) at the interfaces of blocks are found for updating the normal and shear contact spring forces at time (t+Δt).To investigate the motion of each point in all the blocks Newton's second law is used which is resulted from the known forces acting on them.
Explicit time marching scheme is applied to solve the dynamic equilibrium equations, which by employing the central difference method the movement of the ith grid point can be described individually as and where m i is the lumped mass of the ith grid point and g is the ground acceleration.
The resultant forces F xi and F yi in ( 11) and ( 12) represent the summation of forces applied at the ith grid point including the resultant contact force if the grid point is on the boundary of the block, and the elastic forces owing to elastic deformation of the block.Dynamic response is achieved by using an explicit time marching scheme, which is equal to that used by the explicit finite difference method for continuum analysis.This usually requires a sufficiently small time step to overcome numerical instabilities principally caused by disturbances propagation between adjacent discrete elements during a single step [20,27].

C. Hybrid BEM-DEM Algorithm
In present section, the proposed hybrid BEM-DEM algorithm is described.In order to consider the dam-reservoir interaction, the influence of the reservoir when the dam is subjected to ground excitation should be taken into account.This may be achieved by the assessment of hydrodynamic pressure at any instant of time t, by solving the fluid domain applying (10).In this equation, q s depends on the horizontal (u i ) and vertical (v i ) displacement components of the dam upstream face under water which are calculated by ( 13) and ( 14).The foregoing equations are, in turn, updated employing ( 11) and (12) in which, the resultant forces F xi and F yi are functions of hydrodynamic pressures caused by ground excitation at the DEM grid points of the dam upstream face.To vanquish this problem, a staggering solution technique is used, in which at each time step (Δt), both the hydrodynamic pressure and the dam displacement components are iterated concurrently till a desired level of convergence is attained.For displacement iterations in the staggering scheme, in all nodes of the dam upstream face, the convergence criterion is prescribed by the following norm: 1) ( (15) in which j shows the iteration number and ε is a small convergence tolerance value.Thanks to intrinsic non-linear behavior of the whole system (i.e., hydrodynamic interaction problem coupled with cracked concrete dam) dynamic response, the most straightforward procedure may be considered as the time domain analysis of BEM-DEM formulation.Since the time domain analysis of BE formulation is more intricate in comparison with the frequency domain formulation, the latter has been selected to construct the coupled formulation in the present research [25].This frequency domain analysis (i.e., forward and inverse Fourier transformations) is carried out for very short consecutive time steps of Δt (e.g.0.1s) of ground motion time history to permit the algorithm to follow the correct path of a non-linear phenomenon.It is gainful remarking that, as the maximum sliding displacement occurring does not exceed a couple of centimeters, no considerable errors may be resulted in the boundary elements whose lengths are of the order of several meters.In fact, the present method is an Euler-Lagrangian one, in which an Eulerian and Lagrangian formulation are applied for the reservoir and the cracked concrete gravity dam respectively.In this method, the cracked dam is explained in a Lagrangian frame of reference, which deforms under the reservoir forces, while the Eulerian reservoir mesh does not need re-meshing since the sliding displacements are very small.

III. BENCHMARK EXAMPLES FOR ACCURACY VERIFICATION
A program has been developed to implement the proposed algorithm.Several numerical examples are scrutinized to determine the validity and accuracy of the proposed modeling technique.

A. Rigid Dam Example
In this section, a rigid dam example with infinite reservoir whose depth (h) has a constant value of 180m is excited by ramp acceleration, in which the reservoir foundation is assumed to be rigid.In the analysis, the effects of surface waves are neglected, as in the exact analytical solution [29], which is employed for comparison.The time variation of hydrodynamic pressure on the dam-reservoir interface is acquired by the present algorithm considering both constant and linear boundary elements, and is shown in Figure 3.Comparison of hydrodynamic pressure acting on the base of the dam.
In the figure, the dimensionless hydrodynamic pressure is given by / P C P ah

 
, where a indicates the maximum value of the ramp acceleration and h is the height of the reservoir.

B. Flexible Dam Example
A time history analysis for a flexible dam with infinite horizontal reservoir of 180 m in height, 15 m in width and a constant depth extending to infinity is carried out as the second verification example, using the 1940 El Centro earthquake as input ground motion.The reservoir and the dam bottoms are assumed to be rigid, the sound velocity of reservoir c is 1440 m/s and the water is assumed to be compressible and inviscid with a mass density of 1000 kg/m 3 .Figure 4 shows the time variation of hydrodynamic pressure at the base of the dam for the present study and the same FE-FE model in which, the reservoir is modeled by finite elements in the same method as the dam.One can judge by considering the figure, that the present algorithm gives good results in comparison with benchmark solution results.Variation of hydrodynamic pressure at the bottom of the dam under El Centro earthquake's ground acceleration.

IV. SEISMIC ANALYSIS OF FRACTURED KOYNA DAM
Koyna dam experienced a strong earthquake of 6.5 R magnitude on December 11, 1967.During the earthquake, the dam suffered severe damage encompassing horizontal cracks on both the upstream and downstream faces on a number of non-overflow monoliths.In this section, seismic analyses are carried out to observe the dynamic behavior of non-overflow monoliths of cracked Koyna dam (Figure 5) with a penetrated horizontal crack subjected to various ground motions.For this purpose, two real earthquake ground motions are selected: the 1967 Koyna earthquake in the stream direction with a PGA of 0.49g and the1994 Newhall earthquake with a PGA of 0.678g (Figure 6).It is observed that the peak value of bottom block's acceleration due to Koyna earthquake at the crack level is approximately 2.89 g which is about 5.9 times the value of the earthquake input, while the analogous acceleration of the top at the same level is a little smaller at 2.41 g.Furthermore, the shape of the acceleration at the top of the dam is mostly different from those at the crack level while the peak value is approximate.Similarly, due to Newhall earthquake, the corresponding accelerations of bottom and top blocks at the crack level are 4.35 g and 3.62 g respectively.According to above results, the ground motion has been amplified about 5.4 and 5.88 times under Koyna and Newhall earthquakes respectively.
Figures 10 and 11 plot the horizontal and vertical displacement histories of the crest of the cracked Koyna dam obtained by the present hybrid BE-DE method compared to the results of FE-DE technique, in which, the positive values of the horizontal and vertical displacement are in the downstream and upward directions respectively.It is observed that there is no damage during the moderately small amplitude shaking.The first large horizontal displacement of the dam in the downstream direction based on present algorithm is 2.65 cm at 4.02 s due to Koyna earthquake, while the corresponding value due to the other earthquake is 5.56 cm at 4.72 s. Figure 12 shows the results of horizontal sliding of the top block under foregoing ground motions obtained by the two numerical coupling models.The figure denotes that the residual displacements calculated by using the present hybrid BE-DE approach due to Koyna earthquake is 13.1mm, which is comparable to the residual displacement of FE-DE analysis (12 mm).Similarly the corresponding displacement due to Newhall earthquake is 20.3 mm for the present study and 18.4 mm for the FE-DE model.Although the present study results are in reasonable agreement with the results of FE-DE model, details of the two plots are somehow different mainly because damreservoir interaction has been evaluated using FE analysis.Time variation of upper block rotations of the cracked Koyna dam under both ground motions is illustrated in Figure 13.We observe that the upper block is safe against overturning, as the peak rotations of the top block are 0.0018 and 0.0026 under Koyna and Newhall earthquake records, respectively.Acceleration time histories at different positions under the1994 Newhall earthquake.

V. CONCLUSION
A numerical method has been developed to analyze the seismic response of cracked concrete gravity dams placed on rigid foundation subjected to the excitation of horizontal seismic waves including dam-reservoir interaction effects.This is achieved by coupling DEM, employed to model the cracked dam, and BEM, used to analyze the dynamic behavior of the reservoir.In the proposed method, the two different media of cracked dam and reservoir are solved independently and coupling effects are taken into account through the propounded iterative scheme where equilibrium condition at the damreservoir interface is satisfied.The proposed hybrid DEM-BEM approach elucidated in this paper presents the following chief advantages in comparison with other numerical techniques accessible for seismic analysis: 1. Due to the reduction of problem dimensions in BEM, this numerical method performs the dynamic analysis with very small size of matrices of reservoir domain and considerably less computational efforts.
2. Reservoir domain dimensions may be remarkably reduced as the far end of the infinite fluid domain is directly considered with no boundary elements.
3. DEM is applied to model the fractured concrete gravity dam, which is widely employed for the investigation of blocky media.
4. A satisfying level of convergence is simply attained through a few iterations.Several harmonic and transient excitation problems are examined employing the proposed method, and are compared with the results from other analytical solutions and numerical methods existing in the literature to appraise the accuracy and validity of the present method.
Result comparison shows that the present method works well.Furthermore, considering dam-reservoir interaction effects, seismic performance of the Koyna Dam with a horizontal crack at a height has been evaluated.It can be concluded that the appropriate modeling of the damreservoir interaction is necessary in the seismic behavior of fractured gravity dams.
Considering the above, this method seems to be well adapted to handle more complex discontinuous fluid-structure domains involving finite and semi-infinite subdomains.Of course, more parametric studies are proposed in order to calculate optimized time steps that may significantly improve the present algorithm.

Fig. 3 .
Fig. 3.Comparison of hydrodynamic pressure acting on the base of the dam.

Fig. 4 .
Fig. 4.Variation of hydrodynamic pressure at the bottom of the dam under El Centro earthquake's ground acceleration.

Fig. 5 .Fig. 6 .
Fig. 5.Details of fractured Koyna dam with a penetrated horizontal crack from upstream to downstream, at its top profile.

Fig. 8 .
Fig. 7.Hydrodynamic pressure time histories at the dam bottom for different conditions under various ground motions Fig. 9.Acceleration time histories at different positions under the1994 Newhall earthquake.