The Failure Mechanism of Concrete Gravity Dams considering Different Nonlinear Models under Strong Earthquakes

The paper focuses on the failure process and mechanism of the concrete gravity dam considering diﬀerent nonlinear models under strong earthquakes. By taking a typical monolith of a concrete gravity dam as a case study, a comparative analysis of the failure process and mechanism of the dam considering the plastic damage model and the dynamic contact model, respectively, is performed using the seismic overload method. Moreover, the ultimate seismic capacity of the dam is evaluated for both of the nonlinear models. It is found that the ultimate seismic capacity of the dam is slightly diﬀerent, but the failure process has signiﬁcant distinctions in each model. And, the damage model is recommended when the conditions permit.


Introduction
Seismic safety of high concrete dams is extremely important. Once a reservoir dam holding hundreds of millions of tons of lake water fails, it will cause unimaginable catastrophic consequences [1]. To ensure the seismic safety of dams, it is necessary to study the seismic response analysis and failure mechanism of high dams in the strong earthquake area further. e damage, cracking, or even collapse of dams could happen leading to the stress redistribution when the stress exceeds the allowable strength under strong earthquakes. erefore, the assumption of linear behavior is no longer applicable, and the nonlinear situation of the material should be considered to conform to the actual situation better [2].
e Standard for Seismic Design of Hydraulic Structures [3] also stipulates that, for gravity dams of a seismic fortification of class A with complex structures or complex geological conditions, the nonlinear influence of materials should be considered in the finite element analysis. In view that the concrete material is a kind of quasi-brittle material, the tensile damage occurs before the compressive damage generally. erefore, the tensile damage of dam concrete is the major concern ignoring the impact of compression damage in the nonlinear seismic analysis [4]. However, for high concrete dams in the strong earthquake area, the upstream water pressure is already huge. Once subjected to strong earthquakes, the static and dynamic comprehensive compressive stress of the dam at the dam toe often exceed the compressive strength of the dam concrete causing compression damage.
us, it is necessary to study and demonstrate the influence of concrete compression damage on the damage process, mechanism, and ultimate seismic capacity of the dams.
Besides, the tensile strength and shear strength of the RCC gravity dam with the rolled surface are usually lower than that of normal concrete, which may lead to the decrease of stability and safety of the dam [5]. e Code for Seismic Design of Hydraulic Structures [6] also stipulates that the stability of the rolled layer should be checked for RCC gravity dams. Moreover, according to the existing research studies, such as Koyna gravity dam [7], Xinfengjiang dam [8], and Sefid Rud dam [9], the static and dynamic synthesis maximum principal stress of the dam head are relatively large [10]. It can be concluded that the dam head is a weak part of gravity dams, and it is easy to form macrocracks in the upstream and downstream under the strong earthquakes [11]. Although no overall instability occurs under strong earthquakes, penetrating cracks in the upstream and downstream at the head of the dam appear [12]. erefore, the influence of weak layers and parts of the dams should be considered in the seismic analysis.
At present, a more commonly used model for simulating material nonlinearity is the plastic damage model proposed by Lubliner et al. [13] and improved and developed by Lee and Fenves [14,15]. e contact model can be used for gravity dams with deep sliding problems and weak parts with a local cracking inside the dam and the dam-foundation interface that may crack and slip during strong earthquakes [16]. Based on the existing research studies, the concrete plastic damage model [17,18] and the dynamic contact model [19,20] are used for seismic analysis of a typical nonoverflow monolith of a gravity dam in this paper, respectively. e impact of tensile and compressive damage on the seismic damage process and failure mechanism of the dam under strong earthquakes is analysed and discussed in detail. Moreover, the ultimate seismic capacity of the concrete dam for two models is compared and analyzed, which provides a reference for the engineering application.

Plastic Damage Model.
e concrete plastic damage model [21,22] and both the compression and tensile damage are considered in this paper. e model uses the yield criterion in plastic mechanics to judge whether the concrete is in the plastic state. If it enters the damaged state, the flow rule is used to evaluate the plastic strain (residual deformation), and the stiffness degradation of the material is calculated according to the damage evolution curve.

Yield Function.
When the stress is inside the yield surface, the concrete material is in an elastic state. When the stress is on the yield surface, the concrete material begins to enter the plastic state. e yield function can be expressed as follows: where σ is the effective stress tensor, α and β are dimensionless material constants, S is the effective stress deviator, defined as S � σ + pI, I is the second-order tensor, p is the effective hydrostatic pressure, q is the Mises equivalent effective stress, σ 1 is the maximum principal effective stress, σ c (ε pl c ) is the effective compressive cohesion stress, σ t (ε pl t ) is the effective tensile cohesion stress, and ε pl t and ε pl c are the tensile and compressive equivalent plastic strain; the typical experimental values of the ratio σ b0 /σ c0 for concrete are in the range from 1.10 to 1.16; yielding values of α are between 0.08 and 0.12, where σ b0 /σ c0 � 1.16 and α � 0.12.

Flow Rule.
e plastic damage model assumes the nonassociated potential flow: where ε pl is the plastic strain and λ is the plastic flow factor. e flow potential G is chosen for the Drucker-Prager hyperbolic function: where ψ is the dilation angle measured in the p − q plane at high confining pressure, ϵ is a parameter, referred to as the eccentricity, and σ t0 is the uniaxial tensile stress at failure. is paper refers to the value of parameters in the concrete gravity dam calculation example given in Reference [23], taking ϵ � 0 and ψ � 36.31 ∘ .

Damage and Stiffness Degradation.
When the concrete is subjected to the uniaxial tension and uniaxial compression, the plastic strain can be expressed as where d t and d c are the tensile damage factor and compression damage factor, E 0 is the initial (undamaged) modulus of the material, σ t and σ c are the tensile stress and compressive stress corresponding to the total strain, ε ck t is the cracking strain, and ε in c is the inelastic strain. e stress-strain relations under uniaxial tension and compression loading are, respectively, Under the uniaxial cyclic loading conditions, the stiffness degradation mechanisms of the material can be expressed as where E is the elastic modulus of the material, d is the stiffness degradation variable, s t and s c are functions of the stress state that are introduced to model stiffness recovery effects associated with stress reversals, and the weight factors w t and w c , which are assumed to be material properties, 2 Shock and Vibration control the recovery of the tensile and compressive stiffness upon load reversal; in this paper, w t � 0 and w c � 1. Under the action of multiaxial cyclic loading, the elastic stiffness degradation is assumed to be isotropic, which can be expressed by a single scalar variable d. e stress-strain relation of the viscoplastic model is given as where D el 0 is the initial (undamaged) elastic stiffness of the material, ε is the total strain, and ε pl is the plastic part of the strain.
e equivalent plastic strain rates are evaluated according to the expressions: where _ ε pl max and _ ε pl min are, respectively, the maximum and minimum eigenvalues of the plastic strain rate tensor _ ε pl and r(σ) is a stress weight factor that is equal to one if all principal stress σ i (i � 1, 2, 3) and are positive and equal to zero if they are negative.

Dynamic Contact Model.
In the dynamic contact model, the contact problem established by adding a supplementary equation constructed by the contact constraints to the dynamic equation is solved in the form of point-to-point contact. e explicit integration method of the central differential method combined with the Newmark constant average acceleration method is adopted, and the detailed procedures are shown in Reference [24,25]. e displacement expression of jth DOF of node i can be expressed as follows: where N t ij and τ t ij are the components of N i and τ i in the j direction at time t, respectively, K and C are the stiffness matrix and damping matrix of the structure, F, u, and _ u are the external load, displacement, and velocity column vectors, the subscript ij represents the component of the vector in the j-DOF direction of node i, m i is the concentrated mass on the node i, and Δt is the time step.
For the convenience of derivation, formula (13) is simplified into three parts: To get the normal contact force, taking the contact point pairs i and i ′ for instance, the normal contact condition including the initial tensile strength should meet the geometric condition of normal mutual inviolability: where n i is the normal vector of the contact surface at point i; (16) and (19), and N t i � − N t i′ can be obtained: When the normal contact between i and i ′ occurs, the tangential friction force should be calculated; let , and the tangential contact force can be obtained as follows: When the static friction force exceeds μ s N t i (μ s is the static friction coefficient), it means that i and i ′ are in the dynamic friction state: From the above calculation, u t+1 ij , Δu t+1 ij , and Δv t+1 ij can be obtained, and finally, u t+1 ij can be obtained.  Figure 1. e finite element model of the dam is shown in Figure 2. e dam foundation system is discretized into 14592 nodes and 14320 elements. e element size for the dam body is about 2 m. Material properties of the dam concrete and foundation rock are described in

Static and Dynamic Loads.
e static loads mainly include the self-weight, upstream and downstream hydrostatic pressure, and silting pressure. Westergaard's added mass without considering the compressibility of the reservoir water is used to consider the hydrodynamic interaction between the dam and the reservoir. e normal upstream water level is 2150 m, and the corresponding downstream water level is 2019.25 m. e elevation of silt is 2023.7 m, and its floating bulk density is 8 kN/m 3 with the friction angle of 12°. e horizontal and vertical peak ground accelerations of the design earthquake are both 0.4005 g, and the time histories of normalized acceleration are shown in Figure 3. Figures 4-7 show the damage evolution curve used in this paper, which is obtained by the concrete material tests. According to the CDP model parameters in Reference [27], the compression damage evolution is converted and adjusted appropriately by formula (6).

Discussion of the Results.
e damage evolution process of the concrete gravity dam is discussed in detail through the earthquake overload analysis [28]. e influence of the compression damage on the ultimate seismic capacity and damage mechanism of the gravity dam is mainly studied.

Damage Rupture Process Only considering the Tension Damage Model.
Only the tension damage of the dam concrete is considered, and the earthquake overload analysis is performed with increasing overload coefficients of 1.2, 1.4, 1.52, 1.53, and 2.0. e results are shown in Figure 8. e damage occurs at the dam heel first, and the macrocrack depth with a damage coefficient greater than 0.8 is about 18.5 m under the design earthquake. In view that the bedrock material is considered as linear elasticity, the damage at the dam heel is caused by stress concentration and the strong restraint of the foundation. erefore, the dam heel damage is not discussed in this study.
e damage occurs at the downstream face of the dam head for the overload coefficient of 1.2 and gradually expands to the upstream face with the increasing of the overload coefficient. Penetrating cracks (i.e., when a crack extends completely between the upstream and the downstream faces) appear on the upstream and downstream surfaces for the overload coefficient of 1.53. If the penetrating crack of the dam body is taken as the failure criteria, it is suggested that the ultimate seismic capacity of the concrete dam is 1.52 times the design earthquake. Figure 9 shows the number of damaged elements of the dam head under different overload coefficients, and Figures 10-14 show the time histories of tensile damage variables and relative displacements of point pairs i.e., upper and lower relative nodes of an element under different earthquake overload coefficients. According to these results, the damage begins to appear at the time of 6.47 s for the earthquake overload coefficient of 1.2. Meanwhile, relative displacements of the node pairs change rapidly, and the damage gradually increases. For the earthquake overload coefficient of 1.53, the damage of the dam head slope first appears at the time of 6.43 s. en, the cracking gradually spreads across the upstream surface, and the relative displacements of the node pairs gradually increase simultaneously with the deformation of the new damage element accumulating to the downstream element. e damage of the upstream surface begins to appear at the time of 16.10 s. e maximum tensile strain of the Gauss point of element 3400 is 4.36 × 10 − 5 at the time of 16.11 s, which exceeds the tensile strain corresponding to the ultimate stress. Afterward, the damage expands rapidly with the material entering the softening stage. At the time of 16.16 s, the penetrated cracking through the upstream and downstream is formed. Meanwhile, the relative displacement of the upstream node pairs at the dam head changes rapidly. Moreover, the crack has residual displacement after the earthquake. e damage process will accelerate with the increase of the overload coefficients, and the damage appears and has penetrated the dam head at the time of 3.91 s and 8.18 s, respectively, for the earthquake overload coefficients of 2.0. e relative displacement of the node pairs also has some increase with the expansion of the damage domain. Figure 15 shows the tension damage development process under different overload coefficients considering the tension and compression damage simultaneously. Figure 16 shows the compression damage development process under different overload coefficients.

Damage Rupture Process considering Both the Tension and Compression Damage Model.
It can be obtained that the compression damage near the dam toe expands slightly with the increasing of the earthquake overload coefficients ( Figure 16). e compression damage mainly occurs at the dam heel for small overload coefficients. As the overload coefficient increases, the compression damage occurs at the upstream and downstream slopes and gradually expands and the damage at the dam toe begins to appear. erefore, tension damage and compression damage occur simultaneously at the dam heel.
To explain this phenomenon, the results of the overload coefficient of 1.0 are taken as an example of a detailed discussion. e concrete strength grade at the dam heel is C 90 25, and the initial compressive yield stress is 14.8 MPa in the uniaxial state ( Figure 6). Figure 17 describes that the compression damage first appears at 5.82 s for the overload coefficient of 1.0. Table 2 describes the relevant variables at the integral point of dam heel for the overload coefficient of 1.0 at 5.82 s. e maximum compressive stress near the dam heel does not exceed the initial compressive yield stress (Table 1). Meanwhile, the maximum principal stress at the Gauss point of the element at the dam heel is 1.925 MPa, and the minimum principal stress is -0.839 MPa. Hence, the yield surface function F � 3.276 > 0 can be obtained from formula (1),which exceeds the yield surface. According to the flow rule, the plastic strain tensor can be calculated as 3.093 × 10 − 5 , − 1.338 × 10 − 5 T . Since it is in a biaxial tension    and compression state, the stress weight factor r(σ) is 0.696 by the formula, and the equivalent tension and compression plastic strains are 2.296 × 10 − 5 and 1.365 × 10 − 5 , respectively, by formulas (10) and (11). e tension and compression damage variables are obtained as 0.63228 and 0.00246, respectively, and the compression damage variable is not 0. us, with the increase of the overload coefficient, the tension damage and compression damage appear simultaneously at the upstream and downstream slope. Figure 18 shows the compressive stress-strain of C 90 25 concrete. Table 3 shows the related variable values of the element 1327 near the dam toe under different overload coefficients. It can be seen that the compression damage has occurred for the overload coefficient of 1.0, but the         maximum compressive strain and compressive equivalent plastic strain of the Gauss node are less than the values of the ultimate compressive state. erefore, the dam concrete is in the hardening stage, and the loading capacity is not decreased resulting in a small damage propagation.
From Figure 19, the damaged area is still slightly even for the earthquake overload coefficient of 3.00.
ough the compression damage may appear at the compressive stress concentration region under strong earthquakes, the compression constitutive of the concrete reflects the characteristics of ductile materials with strengthening properties, and the ultimate compressive strain is relatively large. As a result, serious tension damage of the dam has already developed when the compressive strain reaches the softening stage. So, the tension damage of the dam is dominant under strong earthquakes.
Moreover, the tension damage begins to appear at the downstream slope of the dam head for the overload coefficient of 1.2. Penetrating macrocracks through the upstream and downstream appear when the overload coefficient is 1.55. erefore, it is suggested that 1.54 times the design earthquake can be considered as the ultimate seismic capacity, which is similar to that only considering the tensile damage.

Seismic Analysis considering the Dynamic Contact Model
As shown in Figure 20, the joints are preset according to the damage penetrating position of the dam in the damage model. e number of joint pairs from upstream to downstream is 1-16, and its tensile strength is 1.6 MPa, which is consistent with the damage model. erefore, the dynamic contact model with the tensile and shear strength of joints is adopted for the seismic overload analysis of the concrete gravity dam. e overload coefficients are 1.00, 1.44, 1.47, 1.48, 1.49, and 1.50 times of the design earthquake, respectively.
Results show that the preset contact joints did not open until the overload coefficient reaches 1.49, indicating that no     contact joint causes tensile cracks and completely breaks after the earthquake.
To sum up, the earthquake overload coefficient of the cracking at the dam head in the damage model is lower than that in the contact model, but the failure process of the contact model is faster than that of the damage model. e main reasons lie in that the damage cracking of concrete in the damage model is mainly judged by the state of the maximum principal stress at the element Gauss point, and there is no restriction on the damage propagation direction. However, for the contact model, the contact status changes among the separation, adhesion, and contact states. e failure of the contact surface is mainly judged by the normal and tangential force, and the influence area of the node is also larger, so the earthquake overload coefficient of the cracking at the dam head of the damage model is lower than that of the contact model. Besides, in the damage model, the dam concrete will enter the softening stage when the maximum principal stress exceeds the tensile strength. In the softening stage, the damage develops and the stiffness gradually degenerates, while the bearing capacity will not be completely lost immediately. Moreover, the elements near the penetrating part of the dam head are also damaged, which plays a role in energy dissipation. However, for the contact model, once the cracks appear in the normal direction, the contact joints cannot bear the normal tensile  stress, and the tensile bearing capacity will be completely lost, so the damage and cracking process are faster than those of the damage model.

Comparative Analysis of the Ultimate Seismic Capacity.
If the penetrating crack of the dam is taken as the failure criteria, the ultimate seismic capacities of the damage model only considering tensile damage and considering tensile and compression damage are 1.52 and 1.54 times the design earthquake, respectively. For the contact model, it is suggested as 1.48 times the design earthquake. erefore, the difference in the ultimate seismic capacities of the concrete gravity dam obtained by the two different models is slight.
Compared with the contact model, the damage model is more consistent with the actual situation without presetting the contact joints. It only needs to consider the damage softening process of the material. However, the damage model requires a refined finite element mesh leading to a large number of computing efforts. e contact model is more suitable for RCC gravity dams with clear weak interfaces. e coarse finite element mesh is sufficient for the analysis, and computing efforts are reduced. Moreover, it can also give a reasonable result to evaluate the ultimate seismic capacity of concrete dams, while, for the dams without clear weak interfaces, the damage model is recommended when the calculation condition permits.

Conclusion
In this paper, seismic overload analysis is performed by taking a typical nonoverflow monolith of a concrete gravity dam considering the damage model and contact model. e tension and compression damage development process and failure mechanism of the concrete gravity dam under the strong earthquake including its ultimate seismic capacity are of comparative analysis. e main conclusions are as follows: For the concrete dam with the material damage model, with the increase of the overload coefficients, the tension and compression damage may occur simultaneously in multiaxial stress conditions, but the damage of the concrete gravity dam is mainly controlled by the tensile damage. A slight difference in the development pattern of the damage is found between the damage model considering only the tension damage and that considering the tension and compression damage.
Compared with the contact model, the earthquake overload coefficient of the cracking of the damage model is earlier, but the damage development process is slower. Based on the failure criteria of the penetrating crack of the dam, the ultimate seismic capacity of the gravity dam with the damage model considering only tensile damage is 1.52 times the design earthquake, and that, with the damage model considering the tension and compression damage, it is 1.54 times the design earthquake. However, the ultimate seismic capacity based on the contact model is 1.48 times the design earthquake. In summary, the ultimate seismic capacities of the gravity dam under the two models are similar, and the contact model is slightly lower.
It is suggested that the contact model should be used for the RCC gravity dam with clear weak interfaces and coarse meshes, while for the dams without clear weak interfaces, the damage model is recommended when the calculation conditions permit.

Data Availability
All data generated or analyzed during this study are included within this article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.