{'fraction_non_alphanumeric': 0.02367908724589125, 'fraction_numerical': 0.011273599855118395, 'mean_word_length': 5.053165250753631, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear seismic response of arch dams considering joint opening effects and boundary conditions of discontinuous foundation

ARTICLE INFO Article history: Received: January 2021. Revised: February 2021. Accepted: March 2021.


Introduction
Stability of an arch dam under ground motion is generally dependent on the strength of its abutments. Actually, even high safety margins for unexpected ground motions do not guarantee the stability of dam if it is established on uncertain abutments. In addition, due to the complex nature of rock foundation including non-homogeneous materials, existing joints sets and faults and propagation of seismic waves from far or near field, as well as errors due to the simplified analytical simulation, final judgment about the performance of dam will be challenging. Collapse of Malpasset Dam in France in 1959 is an obvious example of failure due to lack of foundation strength. During the past years, extensive research in various fields related to the analysis and design of concrete dams has been done, but the need for more accurate modeling of a dam in a coupling system considering the mass, flexibility and nonhomogeneity of discontinuous rock foundation seems crucial. In the present study, a numerical program for nonlinear dynamic analysis of concrete arch dams is developed in FORTRAN. For this purpose, Karun 4 Dam is considered as a case study.
The main features of this study are considering the effect of contraction joints and choosing a proper boundary condition for far-end which has a direct effect on the accuracy and precision of analytical results. In the analysis and design process of an arch dam, it is necessary to model the following features: 1) Dam reservoir interaction and distribution of hydrodynamic pressure, 2) Reservoir foundation interaction and effects of reservoir bottom sediments, 3) Dam foundation interaction and role of non-homogeneous and discontinuities in bed rock, 4) Nonuniform input of the free-field motions, and 5) Nonlinear behavior of quasi-brittle material of concrete and jointed rock and contact in contraction and peripheral joints of dam body [1,3,4]. Fundamentals and analytical methods of all the above mentioned are outside the scope of this article and just a brief review of the main issues related to research are presented here. Simple and primary models in earthquake analysis of dams are the added mass approach of Westergaard for fluidstructure interaction. Westergaard's analytical solution neglected dam flexibility and water compressibility. As a result, several researchers developed advanced numerical methods based on the finite elements, boundary elements, or Downloaded from nmce.kntu.ac.ir at 13:20 +0330 on Thursday October 14th 2021 both of them to model dynamic dam-reservoir interactions. Two common finite element approaches in fluid domain are Eulerian-and Lagrangian-based formulations [1,2]. The former approach is known as pressure-or potential-based formulations where fluid pressure or velocity potentials are selected as state variables. Thus, special contact algorithms are required at the fluid-structure interface. The Lagrangian approach in fluid domain is an extension of the solid finite element formulation with nodal displacements as degrees of freedom. Therefore, fluid domain is formulated using the same shape functions as structural elements and compatibility at the fluidstructure interface is automatically met. In this case, fluid elements are characterized by a volumetric elastic modulus equal to the fluid bulk modulus (or fluid compressibility), with a negligible shear resistance and Poisson constant being nearly 0.5 to simulate fluid flow realistically. Subsequently, many studies were carried out to improve the boundary condition at the far-end of reservoir in dynamic analysis of coupling system. A radiating condition (such as, Sommerfeld-1949, sharan-1985and Küçükarslan-2004) can be applied at the truncated boundary of reservoir. However, the effect of radiating condition on the solution can generally be negligible if the reservoir length is taken to be three or more times the dam height. The next important point is the interaction between impounded water and rock foundation. The partial absorption of pressure waves at sediment layers of reservoir bottom and lateral sides may significantly affect the magnitude of hydrodynamic forces owing to the response of dams due to ground motion [6]. In the current study, the Lagrangian approach is used for modeling the fluid and sediment domains. The interface elements with low shear stiffness are modeled for common surfaces of fluid and solid elements. Arch dam foundation interaction is related to the bedrock's flexibility, changes of physical properties of rock foundation, existence of joints and faults in the rock, geometry of dam body, water leakage, uplift pressure, etc. Different models were used for foundation modeling (such as mass less/massed and rigid/flexible foundation) in order to determine the seismic response of concrete arch dams. It has been proven that in a nonlinear dynamic analysis, it is crucial to include dam foundation interaction, the foundation's mass, flexibility and radiation damping. In order to model the highly complex behavior of jointed rock masses, the strength and deformability of jointed rock masses should be expressed as a function of joint orientation, joint size, and joint frequency. Moreover, it is not possible to represent each and every joint individually in a constitutive model. Therefore, it is necessary to use a simple technique such as the equivalent continuum method which can capture reasonably the behavior of jointed rock mass.
The finite element method applied in the present study recognizes that the foundation rock should act both as: 1) nonlinear solid element for modeling the jointed rock as an equivalent continuum whose properties represent material properties of the jointed rock, and 2) nonlinear interface element used to account for surface roughness of discontinuities. Therefore, this paper deals with the mass, flexibility and non-homogeneity of foundation rock, in addition to the discontinuity due to the master joints and faults. The shape and size of foundation model must be properly selected. Using the finite element procedure, a spherical and cylindrical system in the lower and upper half model are employed, respectively. A right way to determine the size of the foundation model is based on the ratio of deformation modulus of foundation (Ef) to the elastic modulus of concrete dam (Ec). For a flexible rock foundation with Ef /Ec less than ½, the foundation model should be extended at least twice the dam height in all directions [7]. The definition of suitable boundary conditions related to its surrounded domain is another important part of the modeling procedure. For the present study, governing equations and related boundary conditions consisting of free surface water and far-end truncated boundaries of the reservoir and rock foundation have been selected. By neglecting the effects of gravity waves, a zero-pressure boundary condition is prescribed at the horizontal free surface water (negligible surface tension). This assumption of no surface wave is a common assumption in analysis of concrete dams, particularly for deep reservoirs. For the truncated boundary condition of reservoir and rock foundation, interface elements and boundary elements have been used. Using these boundary conditions does not prevent the sliding in faults during seismic loading which can provide non-uniform excitation. In order to avoid modeling the complete absorption of wave propagation in the foundation boundaries, displacement history is used instead of acceleration history as seismic input in dynamic analysis. In evaluating the seismic performance of mega structures such as concrete dams, it is evident that the manner in which ground motions excite the dam-reservoir-foundation system is of major importance. Variations in the ground motion arise mainly from three sources, e.g.; the wave passage effect, the incoherence effect and the site response effect [8,9]. In this paper, the wave incoherence effects due to the different mass and environmental conditions of rock blocks were considered as non-uniform input sources. In addition, the nonlinear material behaviors of concrete dam and bedrock as well as the nonlinear effects of contraction joints of dam and discontinuities in the dam foundation are recognized by the suitable behavioral models which will be described later.

Fundamentals of analytical program and numerical modeling
In this study, the dam reservoir foundation system is modeled by assemblage of solid and interface elements. The isoparametric 8-node solid elements with 2x2x2 Gauss integration are used for all domains including dam body, reservoir, sediment and rock abutments. Output values of stresses, strains, invariants, principal values, etc. are the simple numerical average of Gaussian point values in center point. Also, eight-node interface elements are used in common surfaces of domains, such as the contraction and perimetral joints in dam body and discontinuities of rock abutments. The elements' formulations support geometric and material nonlinear analysis. Assuming that nonlinearities are limited to concrete dam, rock blocks, contraction joints of dam body and rock discontinuities, the elements' stiffness need to be updated at each iteration. Interface elements are placed between continuum (solid) elements. Notably, in the programing, the capabilities of several open source programs that are interested in the analysis of concrete dams such as ADAP-88 [10], EAGD-SLIDE[11], EACD-3D-96[12] was investigated and useful subroutines and subprogram of the finite elements were used [13,14]. The governing equations of motion for 3D nonlinear dynamic analysis of coupling system (subjected to earthquake loads) were discretized by Newmark's method. By adopting very large time increments, static nonlinear analysis can be accomplished as a dynamic analysis. The discretized nonlinear dynamic equation of motion is given by (K. J. Bathe 1996, Zeinkeiwicz et al. 1972): Where φ and c are the internal friction angle and cohesion, respectively. The 3D failure surface of the Mohr-Coulomb criterion can be expressed in terms of stress invariants "): Where the lode angle is: For the tension cut-off yield function, we have: With tension strength ′ . This criterion can be fully described by the following equation: In every time step, the program will check dam and rock foundation solid elements with Modified Mohr-Coulomb criterion. Figure 1 shows the failure envelope under these combined criteria. Different models have been developed to represent the contact between two surfaces. The cohesive law can be expressed such that the local traction ( ) across the interface is taken as a function of displacement jump (δ) across the cohesive surfaces. Defining a free energy density per unit undeformed area (Φ) such that traction acting in the interface is given by (Needleman and Ortiz-1999): Where δ c denotes the value of δ at peak traction ( max = σ c ) and we have: When the contact surfaces undergo a combination of shear deformation and normal compression, the effective separation δ is computed only from the shear components whereas, under normal compression the cohesive material behaves as a linear spring. Weighting coefficient β defines the ratio between the shear and normal critical tractions. For more details, see the reference [15]. The value of interface stiffness will depend on the roughness of contact surfaces; as well as the relevant properties of filling material and moisture. For an initially closed interface, normal stiffness Kn and tangential stiffness KS are set to have a high value. These values can be estimated from the lowest Young's modulus and shear modulus of the adhesive domain around the contact, according to following relations: Where, m i is a factor that controls contact properties, usually between 0.01 and 100 (only in penetration), E and G are the smaller elastic and shear modulus when considering contact between two different materials, l e is the characteristic thickness of the adjacent solid elements perpendicular to the interface, and A e is the interface element surface. In order to verify the accuracy and validity of the finite element modeling and developed computer code, the tallest monolith with unit width of well-known Pine Flat Dam, a concrete gravity dam in California, which is 122 m high, is selected. A water depth of 116 m is considered as the full reservoir condition. Properties of applied material in modeling are: Ec=22.75 GPa, v=0.2 and ρ=25 kN/m 3 . For nonlinear material analysis, the tensile strength of the concrete is taken to be 3.36 MPa which is 12% of its compressive strength.
The dynamic tensile strength shall be equivalent to the direct tensile strength multiplied by a factor of 1.50 (Raphael 1984, Cannon 1991). The analysis results consist of the weight of the dam, the static pressure of the impounded water and the earthquake excitation of horizontal-x component of Taft-1952 Lincoln California ground motion (S69E Component) with scaling to a PGA=0.4g. Proportional damping in the dam provides a critical damping ratio of 5% in the fundamental vibration mode of the dam. Figure 3 shows the comparison between the results of crest displacement for the first 15 seconds excitation obtained from the present developed program and commercially available ANSYS  program in case of linear material behavior. Also, figure 4 shows the analysis results obtained from developed program considering the effects of the material nonlinearity and base sliding. After verification of the developed program, in the next section, the highest concrete dam in Iran, (Karun-4) is considered for case study and investigation of the influence of more accurate modeling of foundation in the dynamic response of dam.  Table 1. These plates create six large blocks, as shown in figure 7.

Karun
The fluid body of reservoir has been modeled using Lagrangian approach by modified elastic elements to solve hydrostatic and hydrodynamic pressures on dam and a maximum reduction factor of 0.006 is applied to the shear modulus to approximate simulation of inviscid flow. The water in the reservoir has a constant depth of 155 m, mass density=10 kN/m 3 , bulk modulus=2131 MPa and Poisson's ratio = 0.495 for nearly incompressible nature. For the upstream/downstream sediments with assumption of about 60/30 meters depth: mass density ρs =15 kN/m 3 , bulk modulus Bs = ρs, and cs 2 = 38400 MPa, where the sound speed profile is estimated from physical sediment properties using Biot theory and assuming cs = 1600 m/s. In all contact surfaces of fluid-solid, the sliding layers are used to enforce the slip condition in order to decouple the tangential displacement components. The developed finite element mesh of reservoir and sediments is shown in figure 8.   Interface elements are used for modeling the rock discontinuities, vertical joints between cantilevers and the peripheral of dam in connection with canyon rock, as well as contacts between the reservoir and surrounding domains with negligible shear stiffness. At the truncated boundaries of reservoir and rock foundation, the appropriate methods such as Boundary Element and interface element are available in the developed numerical program which can be applied on the boundaries. The reference [16] gives the integrated description of BEM. In the case of dynamic analysis, using the interface elements can provide high analysis efficiency and present good estimation. Therefore, the interface elements have been used in far-end of the infinite domain to present modeling of the analytical system (called "Moving B.C." later). The properties of several interface elements are presented in Table 2. The coupled system model includes 11764 nodes and 9348 finite elements. The element meshes of interface elements (without discontinuity elements of foundation) are shown in Figure 9.
The loads applied to the model consist of static and dynamic loading. Static loads are dead weight and hydrostatic pressure at normal level of water and sediment weight. The effects of temperature, tail water load and uplift are neglected. However, in a complete safety evaluation analysis these loads should be taken into account. Knowledge of the in-situ stress field in rock foundation is a fundamental parameter of the dam analysis. The in-situ stress in rock mass is simply equal to the weight of the overlying material and discontinuities control the magnitude and direction of this stress field. In this research, first the static load of discontinuous rock weight is applied for investigation of in-situ stresses. For this loading case, the dam body should remain free of stress owing to canyon deformation. To overcome this problem, the numerical program has the ability to change the material properties such as Young's modulus and Poisson's ratio in loading steps. The analytical results indicate a high efficiency in correctly applying the in-situ stresses, despite the complexity and sensitivity of loading pattern. The loads applied to the model consist of static and dynamic loading. Static loads are dead weight and hydrostatic pressure at normal level of water and sediment weight. The effects of temperature, tail water load and uplift are neglected. However, in a complete safety evaluation analysis, these loads should be taken into account. Knowledge of the in-situ stress field in rock foundation is a fundamental parameter of the dam analysis. The in-situ stress in rock mass is simply equal to the weight of the overlying material and discontinuities that control the magnitude and direction of this stress field. In this research, first the static load of discontinues rock weight is applied for investigation of in-situ stresses. For this loading case, the   Figure 10.
In order to present the effect of foundation interaction on the seismic response of dam-reservoir system, several cases of massive foundation are chosen considering geometric, material and contact nonlinearity as follows:  C1: Continuous rock foundation-Fixed B.C (without interface elements between the rock blocks and on the truncated boundary);  C2: Discontinuous rock foundation-Moving B.C (with interface elements between the rock blocks and on the truncated boundary);  C3: Condition "C2" with applying reduction factor of 10% for deformation modulus and allowable bearing capacity of rock blocks RB1, RB2, RB3 and RB4 (demonstrated in Figure 7) materials; and  C4: Condition "C2" with applying reduction factor of 10% for deformation modulus and allowable bearing capacity of rock blocks RB5 and RB6 (demonstrated in Figure 7) materials. Figures 11, 12 & 13 show the crest displacement of the crown cantilever in upstream-downstream direction. As can be seen, modeling of foundation with more details plays a crucial role in coupling system analysis. Also, using the interface elements with appropriate characteristics on the far-end boundary and major fault zones of the foundation changed the seismic response of dam significantly. It should be noted that this boundary condition and modeling discontinuity in bedrock are critical for an actual response of the dam rather than applying non-homogenous material of foundation, as compared in figures 12 and 13. Generally, comparison of the results show that accurate modeling of the dam (C2 Case) in interaction with the rock abutments greatly affects the behavior of the dam and its effect cannot be ignored. To evaluate the results, Table 3 shows the comparison of the maximum amount of dam crown displacement. Figure 14 shows a comparison of the computed dynamic response of the sliding (relative displacement) between the two crown monoliths at the crest and mid-height levels and on the upper level for the first joint near the left abutment for the C2 case. In addition, sliding between the upper levels of adjacent monoliths has a much greater effect on the seismic behavior of the dam. To prove the significant effect of dam body joints on the dynamic response, according to Figure 15, a cracking scenario is considered and the Natural frequency of the system without contraction joints is evaluated. Figure  16 shows the frequency difference between the first and second modes under the influence of a crack in the dam body. It should be noted that crack modeling carried out as a reduction in the mechanical properties of the crack zone elements. The contour of maximum principal stress obtained from nonlinear analysis for cases of C1 and C2 with a deformed scale of 20 at time 20 seconds are shown in Figure  17. In this figure, the discontinuity of the stress contour due to the existence of contraction joints is clearly shown.

Conclusions
In this paper, the earlier analytical procedure to evaluate the seismic response of arch dams, considering the various effects of dam-foundation interaction in time domain has been investigated by implementing the effects of inertia and flexibility of foundation rock in analysis. For this purpose, the far-end boundary condition and major discontinuity of the foundation are modeled by interface elements. The results obtained from this study show that applying the displacement time history on the model with material non-homogeneity in foundation is an important factor in seismic response of arch dam. However, including the interface elements on the far-end boundary of the foundation and foundation discontinuities affected significantly the response of dam compared with applying non-homogenous material of foundation. In the case of Karun 4 arch dam under the load combination with earthquake, the maximum joint sliding can reach 9.5 cm, in which the faults in rock foundation not only have large influences on the opening of joints of the dam body but can also result in larger openings and serious damage.