Numerical Methods in Civil Engineering Push-over analysis of concrete gravity dams due to flood inflow

: This study numerically predicts the failure scenario of a cracked concrete gravity dam through a push-over nonlinear analysis. The mixed-mode Linear Elastic Fracture Mechanics (LEFM) was employed using the extended finite element method (XFEM). The dam base was considered to be fully fixed without any foundation effects. The hydrostatic pressure on the dam's upstream face was increased incrementally. The failure scenario type of the concrete gravity dams with an initial crack at an arbitrary level exhibited a ductile manner. A mixed-mode crack path with a positive mode II was observed in all initial crack levels. Based on the results, it can be claimed that provided water intrusion prevented, "cracked" concrete gravity dams may not be at noticeable risk of failure due to flood inflow.


Introduction
The push-over analysis (POA) is mainly used to perform nonlinear analysis of structures. In this method, the magnitude of the applied load or displacement is increased incrementally according to a predefined function up to the point of structural collapse. POA can be classified based on the type of the load vector as force-based, displacementbased, and energy-based, where the displacement is applied based on energy (area under the force-displacement curve) increments [1]. The load pattern can be selected according to the type of service or ultimate loads that the structure will experience. In concrete gravity dams, it can be determined based on the hydrodynamic or hydrostatic pressure distribution. Utilizing the static push-over analysis (SPOA), Azizan et al. studied a concrete gravity dam based on the crack profiles while assuming no sliding and rigid foundation. They used two types of lateral loads; inertial and hydrodynamic. The limit states proposed in this study are the yielding state, first crack initiation, and ultimate state, which is identified based on the crack pattern on the structure model [2]. Kamali and Vahdani proposed an energy-based nonlinear static analysis of concrete gravity dams by converting the force-displacement diagram into an energy capacity graph. The push-over capacity curve is calculated considering that, in each step, the work of the floor forces is equal to the structure's internal work and is demonstrated in terms of energy capacity. The procedure was applied to the Koyna Dam, and the results from the nonlinear dynamic procedure and SPOA were compared. The results show that the energybased analysis can be used effectively to estimate the seismic capacity of a concrete gravity dam under a given ground motion [3]. Alembagheri and Ghaemian proposed a new concept to determine the damage state in concrete gravity dams assuming no sliding plane and rigid foundation using nonlinear SPOA and incremental dynamic analysis (IDA) with the case study of the Pine Flat Dam. The dam's performance and its various limit states were determined utilizing the latter methods. It was concluded that SPOA and IDA could develop some indexes for seismic performance evaluation and damage assessment of concrete gravity dams [4]. Utilizing the smeared crack method and assuming a linear elastic compressive behavior of the homogeneous isotropic concrete, Bhattacharjee and Léger statically studied the cracking scheme in a concrete gravity dam with a rigid foundation. Initially, a full reservoir hydrostatic pressure was applied. Then incrementally, the water overflow pressure was increased [5]. Therein, the overflow was considered the net water height above the dam crest. Dias investigated tensile crack propagation in concrete gravity dams utilizing a procedure that inserts specific strain fields to enhance the performance of the underlying finite elements in fracture modeling. Representative numerical simulations verified the accuracy and robustness of the methodology [6]. Modeling moving discontinuities within the classical finite element method is quite cumbersome due to the necessity of the mesh to conform to discontinuity surfaces. Thus, the generalized finite element method (G-FEM) and the X-FEM have been developed to facilitate the modeling of arbitrary moving discontinuities through the partition of unity enrichment of finite elements (PUM), in which the main idea is to extend a classical approximate solution basis by a set of locally supported enrichment functions [7]. Floods are and will continue to be a critical problem in many countries. The possible mitigation of floods by dams and the risk to dams from floods are fundamental problems that need to be addressed [8]. Some attempts have been made to find the failure scenario of a concrete arch [9,10,11] or gravity dam [12,13]. However, there is no comprehensive study investigating the behavior of concrete gravity dams under flood loading. Thus, the POA of concrete gravity dams due to flood inflow was investigated in this study.

Governing Equations and Numerical Approaches
According to the above, the equations used in this study include the cracked and non-cracked parts of the structure. The details of these equations are explained below.

Displacement Components in LEFM-based XFEM
In the XFEM, there are two types of degrees of freedom from two different origins; geometrical and mechanical, respectively simulating the displacement jumps across the crack with an approximate function such as Heaviside and the stress singularity in the crack tip vicinity with some asymptotic functions. In other words, as shown in Figure 1, the enriched displacement should be either a split or a tip one. Thus, the XFEM-based displacement approximation can generally be written as in equations 1 and 2.
The XFEM-based displacement approximation, d, can thus be written mathematically as in equations 3 to 6, where N standard , N split , and N tip , are the standard, the split, and the tip shape functions, respectively, and Bi(x,y) corresponds to the functions representing the singularity effects in the crack tip region.
It should be noted that Bi functions are defined in the local polar coordinates firstly. Next, they are transferred to the local cartesian and then to the global cartesian coordinates. As shown in Figure 2, the crack tip is the origin of the local coordinate system (x1,x2), and α is the crack angle with respect to the global coordinate system.
In equation 9, D is the linear elastic stiffness matrix of the solid materials, and ε, which is related to displacement vector (d) in equations 1 and 2, denotes the strain tensor considering small strains and displacements. Equation 10 is obtained mathematically by converting the strong form of equation 9 to the weak form by employing the virtual work principle, satisfying the boundary conditions and definitions mentioned above.
The discrete form of the linear equations, equation 10, is obtained as equation 11.
. = K d f (11) In equation 11, K is the stiffness matrix, d, defined in equations 1 and 2, is the vector of nodal unknowns, and f is the force vector. Details of equation 11 in the XFEM are described precisely in the following sections.

Discrete form of the virtual work equation in XFEM
The constituting parts of equation 11 are explained herein. According to equations 1 and 2, d is rewritten as equation 12.
In equation 12, u, a, and b mean the standard, split, and singular degrees of freedom, respectively. The stiffness matrix, K, and the forces matrix, f, in equation 11, are defined in equations 13 and 14, respectively, where the indexes have been written according to the parts of the displacement field, defined in equation 12. It is worth mentioning that because of LEFM assumptions, fi b is entirely ignored in the governing equations [14].

Subdominant scheme in XFEM
As shown in Figure 4, based on the relative position of an arbitrary point of a split element on the crack path, the body is divided into two subdomains, Ω + and Ω -. So, a sign function is necessary to modify the displacement jump function of these elements.  1 0 Ω Ω In equation 16, x * is the closest point on the crack path to an arbitrary point x, and n is the normal vector to crack path in x * .

CRACK GROWTH
As shown in Figure 5, the geometric deformations of the two crack faces have three independent movement components relative to each other, and each crack movement can be represented by combining these three modes. In the first mode (I), the two crack faces move away in the y-direction (tensile or opening mode), in the second mode (II), the crack faces slide in the x-direction (shear or sliding mode), in the third mode (III), the crack faces slide on each other in the zdirection (tearing or twisting mode) [14].

STRESS INTENSITY FACTORS
The governing equations of stress, strain, and displacement fields near the crack tip in LEFM were detailed in the previous steps. As mentioned previously (in equations 5 and 6), there are some functions for considering the stress singularity effects in the crack tip region. In LEFM, these functions and some other parameters related to stress singularity in the crack tip region are described with the concept of stress intensity factors (SIFs). These parameters are used to evaluate the crack growth criteria, i.e., the path and rate of the crack growth and the stress distribution in the near-field of the crack tip. Thus, it is essential to precisely evaluate the SIFs for a finite element analysis based on LEFM. The boundary and domain integral methods are the most convenient and accurate methods proposed to compute the SIFs in mixed-mode conditions based on interaction energy integrals. Invoking the virtual work principle and applying the divergence theorem, the J contour integral is obtained [14]. The J-integral method was first introduced by Cherepanov [16] and Rice [17] to calculate the energy release rate of a crack, using a local crack tip coordinate system (x1, x2) in which the x1 axis is parallel to the crack faces. Li, Shih, and Needleman performed the computation of J-integral by transforming the contour integral to an equivalent area integral, described in equation 17 [18].
In equation 17, c is a weight function defined over the integration domain, and es is the strain energy density. As shown in Figure 6, c has a value of unity at an interior node and vanishes on an outer node of an element cut by a contour line [7]. For computing SIFs utilizing the J-integral method, an innovative method called interactive integral is used herein [7]. Consider two states of a cracked body: the present real state (1) by (ui (1) , εij (1) , σij (1) ) and an auxiliary state (2) by (ui (2) , εij (2) , σij (2) ); the J-integral for the sum of states (1) and (2) can be written according to equation 18.
where I (1,2) is called the interaction integral for states (1) and (2) as defined in equation 20, where W (1,2) is the interaction strain energy as defined in equation 21.
In LEFM, the energy release rate for 2-D mixed-mode crack problems can be generally defined based on the modes I and II SIFs, i.e., KI and KII, as in equation 22 [17]. KI and KII can be finally obtained according to equation 23 by assuming the auxiliary state (2) as a pure mode I and/or a pure mode II asymptotic field, through evaluating the displacement and stress fields at the crack tip area, and by substituting these variables into equation 19.

Crack growth and angle criteria
The fracture toughness for a brittle material in a pure mode I is denoted by KIC. For a general mixed-mode case, a criterion to determine the angle between the last and the new segments of the crack line and a function of the stress intensity factors are used for evaluating if the crack propagates or not. According to the maximum circumferential tensile stress theory, presented first by Erdogan and Sih [19] based on the state of stress near the crack tip, considering the expression for the tangential stress, σθθ, equation 24 is proposed. The solution of this equation results in the crack propagation angle, θc in equation 25, which can be expressed using the angle between the crack line and the crack growth direction, with the positive value defined in the anti-clockwise direction, in which maximum σθθ occurs [20].
The equivalent stress intensity factor, Keq, is defined in equation 26.

Verification
In order to evaluate the proposed numerical method, the Koyna Dam was selected as the prototype. It is a concrete gravity dam in India with 103 meters of height and has been used frequently as a benchmark in static analyses. In this work, it was loaded incrementally with hydrostatic pressure from the upstream reservoir. Its mechanical properties are shown in Table 1. In addition, the geometry of this dam and the initial crack position with an adopted crack length of 1.92 meters are shown in detail in Figure 8.   In continuous methods, the fracture energy of a material is a fraction of the stress-strain diagram area, as shown in Figure  9. In some discontinuous methods based on LEFM, this parameter is the same as Gc, related to KIc, when KI and KII are equal to KIc and zero, respatively (based on equation 22). The Gc in Table 1 is the fracture energy or the critical released energy per unit length of crack growth under the assumption of LEFM.
a) Bazant's proposed method [21] b) Sha's proposed method [22]  As shown in Figure 10, the predicted crack path for the Koyna Dam case obtained with the method used in this study and the one from the Bhattacharjee study are well-matched. The downward crack path indicates that positive mode II occurs.
According to the classical LEFM, a crack is considered with a straight line, α=π in Figure 2, and no body force and no traction in the singular region. Thus, the crack growth length must be larger than a minimum value so that the crack rotation does not affect the tip nodes. Therefore, considering a maximum element size of about 0.6m in the probable crack path area, a constant crack growth length of 1.9m was considered. Fig. 10: Crack path of the Koyna Dam under upstream hydrostatic load in this study and in Bhattacharjee's work [5] By raising the reservoir level, three stages in the dam crest horizontal displacement diagram were observed: the linear elastic, the plastic, and the hardening stages, as in Figure 11. As shown in this figure, the maximum overflow difference at the starting point of the crack growth between the two studies is about 1.5 and 2 meters. However, the whole level of the water over the crack should be taken into account to make a more accurate investigation. So, as shown in Figure  8, 36.5 meters must be added to the reported overflow in Figure 11. Therefore, the maximum difference reduces to 4%, which is negligible.

Numerical Example
Pine Flat Dam is a 561-meter-long concrete gravity dam. The highest monolith of this dam is 122 meters high, which was considered in this research. The geometric details of the numerical model used in this research are shown in Figure  12. The effect of the initial crack level on the behavior of a concrete gravity dam is investigated herein. Pine Flat Dam was analyzed statically with an initial crack length of 1.9 meters at the levels of 95, 85, 75, 65, 55, 45, 35, 25, and 15 meters above the foundation. The dam was subjected to hydrostatic pressure and selfweight by increasing the reservoir water head with respect to the dam crest level. The results of analyses conducted are shown in Figures 13 and 14 for cracks initiated from different elevations at the upstream face. It was found that hydrostatic pressure and gravity load amplify the behavior of modes I and II positive and mode I negative, respectively. As shown in these figures, the structure behaved in a ductile manner without general structural failure. The downward slope of the step-by-step crack path indicates a mixed-mode behavior, with the positive modes I and II, regardless of the initial crack levels.  As shown in Figure 15, in the Pine Flat concrete gravity dam, with increasing the initial crack level, the overflow level corresponding to the first step of the crack growth has an increasing trend when the initial crack level is 15 to 35 meters over the base. Conversely, if the initial crack level is more than 35 meters over the base, the overflow level corresponding to the first step of the crack growth has a decreasing trend.

Conclusions
Numerous numerical methods have been used to analyze the stability of the cracked concrete gravity dams. In this study, for the first time, a numerical method was used to investigate the failure scenario of a concrete gravity dam with an initial crack due to flood inflow. This method simulates the crack singularity and discontinuity, utilizing a LEFM-based XFEM formulation. Accordingly, Pine Flat Dam was investigated statically with an initial crack at low to high elevations.
The failure scenario type of the concrete gravity dams with a crack at an arbitrary level showed a ductile manner. A mixed-mode crack path was observed in all initial crack levels. Also, by studying the cracked Pine Flat Dam with an initial crack of 1.9 meters in length, it was seen that the corresponding overflow of the first step of crack growth shows a dual effect. Indeed, adequate resilience of the Pine Flat Dam against excessive overflows, in case of structural cracks at its lower parts was observed. On the contrary, the upper part of the dam showed considerably less resistance to such hazards. In summarizing the above results, it can be claimed that "cracked" concrete gravity dams may not be at noticeable risk of failure due to flood inflow. However it should be kept in mind that the above conclusion is valid only when probable other phenomena such as the effects of water intrusion into the crack opening are neglected.