1. Introduction

The real seismic response of a high arch dam to earthquake loadings is a very complicated problem and depends on various factors such as interaction of the dam with foundation rock and reservoir water, numerical modelling, selecting appropriate material properties, suitable boundary conditions for system and finally extracting eligible earthquake record for analysing dam-reservoir-foundation system. Therefore, an efficient method is required to properly assess the safety of an arch dam located in the region with significant seismicity. From this point of view, nonlinear seismic analysis of arch dams for earthquake ground motions should be considered based on a detailed analytical model. In addition, this model should be capable of representing nonlinear behaviour of mass concrete in dam body. Interaction between the dam and impounded water is one of key factors affecting the seismic response of arch dams during earthquake ground motion. When dam-reservoir system is subjected to a dynamic effect such as earthquake, hydrodynamic pressures in excess of hydrostatic pressures occur on the dam upstream face due to the vibration of the dam and water in the reservoir. These hydrodynamic pressures and the deformation of the dam interact with each other (Perumalswami and Kar, 1973;[1]). Therefore, the fluctuation of reservoir water level and its effects on the seismic response of arch dams must be considered in nonlinear analysis. The magnitude of hydrodynamic effects on an arch dam during earthquakes is based on the water level in the reservoir. The reservoir water level usually varies based on the operational regime and annual weathering conditions and sometimes reservoir level may be decreased to visit the upstream face or to repair probable damages or even because of drought and they can considerably affect the response of arch dams during earthquakes.

Importance of hydrodynamic effects on behaviour of dams subjected to earthquake ground motions has long been recognized by investigators such as Chopra 1967[2], Chwang 1978[3], Sharan 1985[4], Fok and Chopra 1986[5], Hung and Wang 1987[6], Hung and Chen 1990[7], Tsia and Lee 1990[8], Dominguez and Maeso 1993[9], Fahjan et al. 2003[10], Pani and Bhattacharyya 2007 and 2008[11,12], Mitra and Sinhamahapatra 2008[13], Najib and Fei 2009[14]. They have studied dam-reservoir system with different geometrical features of the dam and the reservoir. In Fok and Chopra studies, the complex frequency functions were used to define dam–reservoir interactions and the concrete material was assumed to be linearly elastic. They found that hydrodynamic pressures play an important role on dynamic response of arch dams. In these studies, the reservoir was generally considered as empty or full. Proulx et al. 2001[15] investigated experimentally and numerically variation of resonant frequencies of EMOSSON arch dam for four different water levels. Akkose et al. 2008[16] studied reservoir water level effects on nonlinear dynamic response of arch dams. In their investigation, concrete idealized as elasto-plastic using Drucker-Prager model and the reservoir was modelled by Lagrangian fluid finite elements. They found that the reservoir water level effects must be considered in nonlinear analysis of arch dams.

On the other hand, equilibrium equations for the system and kinematic relations for compatibility between strain and displacements should be verified for any nonlinear stress analysis problems. In addition, constitutive relations for stress and strain behaviour of the material should be determined. There are several approaches to simulate the stress–strain relationship of mass concrete under various stress states. In this regard, theories based on plasticity and fracture mechanics approaches have been developed and used in most of engineering analyses. In addition, several researches have been conducted to study dynamic behaviour of concrete arch dams. Kuo 1982[17] suggested interface smeared crack approach to model contraction joints and used this technique for dynamic analysis of arch dams. Dungar 1987[18] used a bounding surface model for static analysis of an arch dam. Cervera et al. 1995[19] utilized a continuum damage model to analyse Talvachia arch dam in Italy. In another research, Noruziaan 1995[20] applied bounding surface and orthogonal smeared crack models for compressive and tensile regions of stresses, respectively. However, numerical instability was reported in some cases (Noruziaan 1995). Hall 1998[21] proposed a simple smeared crack model to simulate contraction joints while diagonal cracking is neglected in his work. Espandar and Lotfi 2002[22] applied non-orthogonal smeared crack and elasto-plastic approaches to Karadj arch dam in Iran. Mirzabozorg et al. 2003[23] developed the staggered method to analyse the dam-reservoir interaction in 3D space. Also, Mirzabozorg et al. 2004[24] utilized damage mechanics approach to conduct seismic nonlinear analysis of concrete gravity dams in 2D space including dam-reservoir interaction effects. Mirzabozorg and Ghaemian 2005[25] developed a model based on smeared crack approach in 3D space. In their work, they analysed 3D models including dam-reservoir interaction effects and considered nonlinear behaviour of the structure, Ardakanian et al. 2006[26] considered nonlinear seismic behaviour of mass concrete in 3D space which is based on an anisotropic damage mechanics model, Mirzabozorg et al. 2007[27] investigated non- uniform cracking in smeared crack approach for 3D analysis of concrete dams, Mirzabozorg et al. 2010[28] studied nonlinear behaviour of concrete dams under non-uniform earthquake ground motion records.

At the present study, reservoir level fluctuation effects on nonlinear seismic response of a high arch dam including dam-reservoir-foundation interaction is investigated utilizing Eulerian approach in fluid domain and Lagrangian approach in solid domain. For this purpose the finite element model of dam-reservoir-foundation was developed. In this model dam body is modelled defining concrete materials property and considering smeared crack method (NL-SCM) and dam-reservoir-foundation interaction is taken into account (FSI-SSI). The heights of water in reservoir are considered in four different levels so that in the first 28.9%, in the second 52.1%, in the third 79.9% and in the forth 96.9% of the total height of the simulated dam is filled in order to study responses of the body.

2. Modeling Material Nonlinearity in Mass Concrete

In theories based on continuum crack models (CCM), the fracture is idealized to propagate as a blunt front smeared over an entire element or a certain band width of the element. After initiation of fracture process, determined by a suitable constitutive model, the pre-crack material stress-strain relation is replaced by an orthotropic relation with material reference axis system aligned with the fracture direction. The tension stiffness across the crack plane is either eliminated suddenly or a gradual stress-release criterion is applied. Thus, only the constitutive relation is updated with propagation of cracks and the finite element mesh is kept unchanged. The advantage of the models lies in its simplicity and cost effectiveness, although the physical nature of crack representation is questionable (Gunn 2001[29]).

One of the most important branches of CCM is nonlinear-smeared crack model (NL-SCM), which considers concrete nonlinear characteristics. These kinds of model must be capable of describing material constitutive behaviour before and after cracking. In the present model it is assumed that concrete material is initially (before cracking) isotropic and linear. The stress-strain matrix is defined by the following equation:

(1)

Figure 1. Typical uniaxial compressive and tensile stress-strain curve for concrete

where E is isotropic Young’s modulus for concrete and ν is Poisson’s ratio. The behaviour of concrete elements is assumed linear elastic until it reaches ultimate strength (Figure 1) and the elasticity modulus of concrete is considered as the average, E, instead of the linear actual, E_o. Stress increases linearly along with an increase in strain. In this step, each reloading of elements leads to elastic returning of strain. The tension failure of concrete is characterized by a gradual growth of cracks, which join together and finally disconnect larger parts of the structure. It is a usual assumption that forming cracks is a brittle process and the strength in the tension-loading direction abruptly goes to zero after big cracks or it can be simulated with gradually decreasing strength.

In concrete, cracking occurs when the principal tensile stress in any direction lies outside the failure surface. After cracking, the elastic modulus of the element is set to zero in the direction parallel to the principal tensile stress direction. In this model cracking is permitted in three orthogonal directions at each integration point. When cracking occurs at an integration point, the stress-strain relation is modified by defining a weak plane normal to the crack direction, which is unable to endure any tensile stresses. The presence of crack at an integration point and in special direction, represent through modification of stiffens matrix by exerting shear transfer coefficient in cracked plane. Also crushing occurs when all principal stresses are compressive and lie outside the failure surface; in this condition, the elastic modulus is set to zero in all directions.

Based on the fact that concrete has been cracked in one, two or three orthogonal directions, the stiffness matrix can be represented in following forms:

I) Concrete has been cracked in one direction and the crack is open:

(2)

where E_s is the secant modulus of elasticity, υ is Poisson’s ratio and β_open is open shear transfer coefficient and defined as the factor that represents shear strength reduction across the cracked face.

II) Concrete has been cracked in one direction and the crack is closed:

(3)

where β_close is closed shear transfer coefficient. The shear transfer coefficient, β, represents conditions of the crack face. The value of β ranges from 0.0 to 1.0, with 0.0 representing smooth crack (complete loss of shear transfer) and 1.0 representing a rough crack (no loss of shear transfer). The value of β used in many studies of reinforced concrete structures, however, varied between 0.05 and 0.25 (Bangash 1989,[30]). Convergence problems were encountered at low loads with β less than 0.2.

III) Concrete has been cracked in two directions and the cracks are open:

(4)

IV) Concrete has been cracked in two directions and both cracks are closed:

(5)

V) Concrete has been cracked in three directions and the cracks are open:

(6)

VI) Concrete has been cracked in three directions and all cracks are closed. In this situation equation 5 can be written again. It must be noted that all above stress-strain relations are written in coordinate system that is parallel to principal stress directions.

3. Finite Element Formulation of Reservoir

Hydrodynamic pressure distribution in reservoir is governed by pressure wave equation. Assuming that water is linearly compressible and neglecting viscosity, small- amplitude irrotational motion of water is governed by 3D wave equation given as (Ghaemian and Ghobarah 1998[31]):

(7)

where P is hydrodynamic pressure and c_o is velocity of pressure wave in water. For earthquake excitation, condition at boundaries of the dam-reservoir, the reservoir-foundation, the reservoir-far-end and the reservoir-free-surface as shown in Figure 2 are governed as presented in following sections.

Figure 2. Boundary conditions for dam-reservoir-foundation

3.1. Dam-Reservoir Boundary Condition

At the surface of fluid-structure, there must be no flow across the interface. This is based on the fact that face of concrete dams is impermeable. In the following equation, superscript “s” refers to the structure.

(8)

where a^s_n is normal acceleration of the dam body on the upstream face and n is normal vector on the interface of the dam-reservoir outwards the dam body and finally, ρ is the density of reservoir water.

3.2. Reservoir-Foundation Boundary Condition

If there is no energy absorption on reservoir bottom, the same boundary condition represented in equation 8 can be used. Reservoir bottom absorption effect is implemented as:

(9)

where q is admittance coefficient and relation between q and α (wave reflection coefficient at the reservoir bottom and sides) is expressed as:

(10)

3.3. Free-Surface Boundary Condition

In high dams, surface waves are negligible and hydrodynamic pressure on the free surface is set to be zero.

3.4. Reservoir-Far-End Boundary Condition

For modelling far-end truncated boundary, viscous boundary condition (called as Sommerfeld boundary condition) is utilized to absorb completely the outgoing pressure waves given as:

(11)

4. Finite Element Model of Dam- Reservoir-Foundation System

DEZ is a high double curvature arch dam with a peripheral joint separating the main body from its concrete saddle structure (PULVINO foundation). Its total height is 203m but the height of the dam above concrete plug (the simulated dam) is 194m. Also the height of main body (without PULVINO) is 186m. Thickness of the dam at the crest is 4.5m and its maximum thickness at the base is 21m.

4.1. Finite Element Model

The finite element model prepared for dam, foundation rock and reservoir is presented in Figure 3. The model consists of 792 8-node solid elements for modelling concrete dam and saddle and 3770 8-node solid elements for modelling foundation. The 8-node solid elements have three translational degrees of freedom (DOF) at each node. The reservoir is modelled using 3660 8-node fluid elements. Utilized Eulerian fluid elements have three translational DOFs and one pressure DOF in each node. It should be noted that translational DOFs are active only at nodes that are on the interface with solid elements.

Figure 3. Finite element model of the dam, foundation and the reservoir

4.2. Material Properties

Material properties for mass concrete and foundation are described in Table 1. Furthermore, reservoir water density is assumed 1000kg/m³, sound velocity is 1440m/s in water and wave reflection coefficient for reservoir around boundary is 0.8 conservatively.

Table 1. Material properties of the mass concrete and foundation rock Staticproperties Dynamicproperties MassConcrete Isotropic Elasticity 40GPa 46GPa Poisson’s Ratio 0.2 0.14 Density 2400kg/m3 2400kg/m3 Uniaxial Tensile Strength 3.4MPa 5.1MPa Compressive Strength 35MPa 36.5MPa FoundationRock Isotropic Elasticity (Saturated) 13GPa 13GPa Isotropic Elasticity (Unsaturated) 15GPa 15GPa Poisson’s Ratio 0.25 0.25

5. Loading History

Applied loads on the system are dam body self-weight, hydrostatic pressure in various water levels and earthquake load. It should be noted that thermal loads is not considered in the current study because of the main object of the study which is considering fluctuation effects of reservoir water level. The β-Newmark method is utilized to solve the coupled nonlinear problem of dam-reservoir-foundation model.

The system is excited at the foundation boundaries using TABAS earthquake records. Total time duration of earthquake is 48.9s. For specifying significant duration of excitation, Arias intensity on Husid diagram was used and then time interval between 10% and 80% of Arias intensity (totally 70%) was selected (Trifunac and Brady 1975[32]). With this method significant time duration is obtained to be 10s (as shown in Figure 4). Moreover structural damping is taken to be 7% of critical damping.

Figure 4. TABAS earthquake recorded at TABAS station (L: stream direction, T: cross section direction, V: vertical direction)

6. Numerical Results

Considering the great importance of natural vibration modes, 30 first periods of dam-reservoir-foundation system in four different water levels for monolithic dam body were extracted (as shown in Figure 5). It is found that the amount of period grows with arising reservoir water level which is because of increasing total mass of the system and differences are meaningful for ten first modes.

Figure 5. Comparison of dam periods for monolithic body in four different reservoir levels

Initial (static) stresses and displacements have strong effects on seismic response of arch dams. Therefore, static analysis of dam-reservoir-foundation rock system under self-weight considering stage construction and impounding is carried out to establish the initial condition for seismic analysis. Subsequently, nonlinear seismic analyses of the system are performed.

Water levels in the reservoir are considered as 56, 101, 155 and 188m that stand to level I, level II, level III and level IV to investigate reservoir level variation effects on nonlinear seismic response of the selected arch dam. Figure 6 shows different reservoir levels.

Figure 6. Reservoir levels considered in linear and nonlinear analyses

Figure 7 shows displacement time-history for half of the dam crest in maximum water level extracted from linear elastic model. As can be seen, maximum displacement belongs to the block number 9 which is one of the central blocks of the dam. Hence all the investigations on blocks limit to this one.

Figure 7. Displacement time-history for half of the dam crest in maximum water level

Table 2. Static (initial) displacement along the crest points in stream direction; NL-SCM Centre point (mm) Left quarter point (mm) Right quarter point (mm) Level (IV) 20.96 9.06 12.49 Level (III) 5.02 1.03 1.27 Level (II) -0.5 -0.62 -1.83 Level (I) -1.34 -0.96 -2.36

6.1. Study along Crest Nodes

Time-history of crest displacement at the central cantilever resulted from four conducted analyses corresponding to four water reservoir levels are shown in Figure 8. Table 2 represents initial (static) displacement along the crest in stream direction for different reservoir water level. Considering that the positive direction of y axis is toward downstream (DS), it becomes evident that in low levels of the reservoir, the main body of dam tends to move in upstream (US) direction.

Figure 8. Time-history of US/DS displacement at crest of the crown cantilever; NL-SCM; (a) Level IV, (b) Level III, (c) Level II, (I) Level I

Table 3 represents maximum displacement of the crest in US/DS direction during seismic excitation. As can be found, raising the reservoir level leads to increasing the crest displacement in the downstream direction and decreasing it in the upstream. However, this variation is negligible between reservoirs levels I and II.

Table 3. Maximum displacement of crest in stream direction; NL-SCM Centre point (mm) Upstream Downstream Level (IV) -19.1 77.2 Level (III) -36.2 45.4 Level (II) -45.8 39.5 Level (I) -44.8 36.5

Table 4 represents absolute maximum velocity and acceleration of the crest centre in the stream direction. It is found that arising reservoir level from I to IV leads to decreasing in velocities and accelerations. However, in high levels, these changes are negligible. The reason can be attributed to free movement of the dam body in low water levels.

Table 4. Absolute maximum velocity and acceleration at the crest centre; NL-SCM Velocity (m/s) Acceleration (m/s2) Level (IV) 1.093 24.0 Level (III) 1.064 27.54 Level (II) 1.087 31.75 Level (I) 1.157 40.64

6.2. Study along the Height of Crown Cantilever

Figure 9. Non-concurrent displacement envelope for upstream nodes of central cantilever; NL-SCM

Figure 9 represents non-concurrent displacement envelope for upstream nodes of the central cantilever. As shown, displacement in lower half of the central cantilever grows significantly with increasing of reservoir level. Displacements corresponding to levels I and II are approximately the same at upper half of the cantilever. In addition, maximum displacements are the same for levels I, II and III at the crest of central block, while the crest displacement corresponding to level IV of the reservoir water is about 65% more than the others.

Figure 10. Non-concurrent velocity envelope for upstream nodes of central cantilever; NL-SCM

Figure 10 represents non-concurrent velocity envelope along the height of the central cantilever. As shown rate of velocity variations in crown cantilever remains same for different levels of the reservoir water. In lower half of the cantilever maximum velocities are approximately the same in all considered cases. Although there is little difference in upper half, rate of growth is similar. A sudden increase can be seen in the height of 140m for all cases but at last at crest of the central block, velocities corresponding to low water levels are more than those obtained for high water levels.

Figure 11. Non-concurrent acceleration envelope for upstream nodes of central cantilever; NL-SCM

Figure 11 shows non-concurrent acceleration envelope corresponding to nodes considered for the two recent figures. As shown rate of acceleration variations along the crown cantilever are the same for all cases like as shown for velocity. Rate of acceleration increases slowly, up to 150m but a sudden growth occurs in this height for all cases so that accelerations in 194m are about triple of accelerations in 150m. Finally at the crest maximum acceleration is corresponding to minimum considered level of reservoir and vice versa.

6.3. Study of Pressure at the Base of the Crown Cantilever

Table 5 gives hydrostatic, absolute maximum value of hydrodynamic and total pressures at base of the main dam body for different reservoir levels. In addition, Figure 12 shows time-history of total pressure at base of the dam for levels I to IV. As shown absolute maximum of hydrodynamic pressure in level IV is more than the rest and in level I is less than the others. Between these levels hydrodynamic pressure varies according to previous maximum and minimum amounts. On the whole the total pressure grows with increasing reservoir level.

Table 5. Pressures at base level of crown cantilever; NL-SCM Hydrostatic (MPa) Hydrodynamic (MPa) Total (MPa) Level (IV) 1.647 0.357 2.005 Level (III) 1.324 0.300 1.624 Level (II) 0.794 0.326 1.120 Level (I) 0.353 0.253 0.606

Figure 12. Time-history of total pressure at base of the dam body; NL-SCM

6.4. Study of Stresses and Crack Profiles within the Dam Body

Figures 13 and 14 show non-concurrent envelopes of maximum and minimum principal stresses experienced by dam body nodes. It should be noted that these parameters must be limited to dynamic tensile and compressive strength of mass concrete. As can be seen, the extensions of high tensile stressed regions are more when reservoir water level is decreased. In addition, due to intensive arch action of the dam body in upper parts due to high curvature of horizontal arches, cantilevers act like as fixed support beams and extension of high tensile regions on the downstream face is expected as can be seen in Figure 13. On the other hand, due to arch action of the dam body, raising reservoir water level leads to increasing compressive stresses as shown in Figure 14. However, as expected from arch dams’ general behaviour, there is a large margin between experienced compressive stresses and dynamic compressive strength of mass concrete.

Cracked areas on the upstream and downstream faces at considered reservoir levels are shown in Figure 15. As it is predicted based on non-concurrent envelopes shown in figures 13 and 14, the crack profiles becomes critical by reservoir dewatering. Generally, cracked regions are located at upper parts on the upstream and downstream faces and some of them are in vicinity of PULVINO and also near the base.

7. Discussion on Results

As shown in Figure 16, comparing displacements resulted from NL-SCM with linear elastic (LE) model along the height of central block; rate of growth along with reservoir impounding is the same for both models. Displacements for levels III and IV coincident to each other in both models but by reservoir dewatering this correspondence losses so that for levels I and II, the upper half of the cantilever (dry region) experience larger displacements in NL-SCM in comparison with those obtained from LE. At the crest, percentage of displacement discrimination for levels I and II are about 11% and 4%, respectively. It can be shown that the same conclusion can be obtained comparing velocities and accelerations resulted from LE and NL-SCM models.

Figure 13. Non-concurrent envelope of maximum principal (S1) stresses for upstream and downstream faces of the dam body; NL-SCM, (a) Level IV, (b) Level III, (c) Level II, (I) Level I

Figure 14. Non-concurrent envelope of minimum principal (S3) stresses for upstream and downstream faces of the dam body; NL-SCM, (a) Level IV, (b) Level III, (c) Level II, (I) Level I

Figure 15. Cracked areas on upstream and downstream faces of the dam body; NL-SCM

Figure 16. Comparing non-concurrent displacement envelopes for upstream nodes of central cantilever

For discussion on stress distribution and crack profiles within the dam body, it should be reminded that in the current study the William-Warnke failure surface in principal stress state is utilized to identify state of stressed nodes. Based on obtained results dewatering the reservoir leads to extension of high tensile stress region on downstream face of the dam body, which is due to acting cantilevers like as semi-fixed supports beams. However, raising the reservoir level causes higher compressive stresses due to dome like action of the body. This conclusion is confirmed looking at extracted crack profiles from NL-SCM model.

8. Conclusions

In the present study, effect of reservoir water level fluctuation on nonlinear seismic behaviour of high concrete arch dams is considered. The reservoir is modelled as a compressible medium and the foundation rock is assumed to be mass-less. Mass concrete is assumed to have tensile cracking and compressive crushing ability. In addition, for identifying stress state of each Gaussian point within elements William-Warnke failure surface in stress space is utilized. For case study, DEZ double curvature arch dam in Iran is considered.

Based on extracted results from analysing NL-SCM models reservoir water fluctuation affects intensely on seismic response of high concrete arch dams in such a way that raising the reservoir water level leads to lowering extension of high tensile stressed regions and increasing compressive stresses within the dam body. In fact, based on the results, it can be confirmed that the worst case for high arch dams occurs in low water levels.