# Modeling of Hydraulic Fracture of Concrete Gravity Dams by Stress-Seepage-Damage Coupling Model

Authors: Sha Sha; Guoxin Zhang
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

High-pressure hydraulic fracture (HF) is an important part of the safety assessment of high concrete dams. A stress-seepage-damage coupling model based on the finite element method is presented and first applied in HF in concrete dams. The coupling model has the following characteristics: (1) the strain softening behavior of fracture process zone in concrete is considered; (2) the mesh-dependent hardening technique is adopted so that the fracture energy dissipation is not affected by the finite element mesh size; (3) four coupling processes during hydraulic fracture are considered. By the damage model, the crack propagation processes of a 1 : 40 scaled model dam and Koyna dam are simulated. The results are in agreement with experimental and other numerical results, indicating that the damage model can effectively predict the carrying capacity and the crack trajectory of concrete gravity dams. Subsequently, the crack propagation processes of Koyna dam using three notches of different initial lengths are simulated by the damage model and the coupling model. And the influence of HF on the crack propagation path and carrying capacity is studied. The results reveal that HF has a significant influence on the global response of the dam.


---

## Body

## 1. Introduction

Hydraulic fracture is a phenomenon of crack propagation after high-pressure water or other kinds of fluid entering into an existing crack. Hydraulic fracture is an important issue in the hydraulic engineering, the petroleum engineering, mining, and geotechnical industries. In the hydraulic engineering, mass concrete dams are likely to experience cracking on their upstream, downstream, and base surfaces due to the low tensile strength of concrete and the action of internal and external temperature changes, shrinkage of the concrete, differential settlement of the foundation, and other factors. Concrete gravity dam is a type of concrete structure that interacts with high-pressure water. With time, cracks are filled with water and penetrate deep into the dam under the action of high water pressure, resulting in reduced carrying capacity and safety of the dam. Therefore, for the safety of high or ultrahigh concrete dams, it is necessary to consider the influence of hydraulic fracturing [1].In hydraulic fracture of concrete and concrete gravity dams, there are some studies on it. Brühwiler and Saouma [2, 3] conducted hydraulic fracturing tests on concrete specimens of different gradations and examined the water pressure distribution in their cracks. They found that the hydrostatic pressure in the cracks was a function of the crack opening displacement, decreasing from the maximum pressure to zero along the fracture process zone. Slowik and Saouma [4] investigated the pressure distribution inside a crack with respect to time and the crack opening rate. The effect of sudden crack closure on the pressure distribution was also investigated. Using the fluid mass and momentum conservation theories, Li et al. [5] established a differential equation of the water pressure distribution in rock and concrete fractures caused by hydraulic fracturing and derived a formula for calculating the pressure at an arbitrary time during the crack propagation process. Bary et al. [6] presented an approach using mechanics of saturated porous media to model strongly coupled hydromechanical effects in concrete. Fang and Jin [7] studied the fracture process of concrete under the action of water pressure in fissure based on the extended finite element method (XFEM). Barpi and Valente [8] simulated hydraulic fracturing of the dam-foundation joint and analyzed the effect of fracture process zone on the path of crack formation. Dong and Ren [9] simulated the propagation of the crack at the gravity dam heel under uniform pressure by the XFEM. Wang et al. [10] studied hydraulic fracturing in concrete gravity dam considering fluid-structure interaction by the XFEM and finite volume method.Hydraulic fracture is a complex problem, involving four coupled processes: (1) the deformation of the surrounding medium induced by the water pressure on the fracture surface; (2) fluid flow within the fracture; (3) propagation of the fracture; (4) the leak-off of the fracturing fluid from the fracture into the surrounding medium [11]. However, among numerical studies of hydraulic fracturing in concrete gravity dams mentioned above, the fourth coupling process is not taken into account. Further, concrete is a quasibrittle material. Fracture in concrete is characterized by the existence of a nonlinear fracture process zone at the front of the real crack tip. Depending on how the fracture process zone of the concrete is assessed, crack analysis can be conducted using linear elastic fracture mechanics (LEFM) or nonlinear fracture mechanics (NLFM). LEFM may produce inaccurate results because of the neglect of the fracture process zone in the concrete. Theoretically, NLFM is more reasonable, being based on the fictitious crack model with the application of the strain softening law of the fracture process zone [12]. However, among the numerical studies of hydraulic fracturing mentioned above, LEFM is used mostly.There are many methods for the simulation of hydraulic fracture, such as phase-field method [13], peridynamics [14], cellular automata method [15], discrete element method (DEM) [16], numerical manifold method (NMM) [17], element free method [18], finite element method (FEM) [19], and boundary element method (BEM) [20]. For the excellent abilities in LEFM, the boundary element method and especially the displacement discontinuity method (DDM) have also been extensively applied to hydraulic fracture modeling. Higher order elements to DDM are introduced to improve the precision due to singularity variations near the crack tip [21]. In the DEM framework, cracks propagate along prescribed element boundaries when the stress intensity factor meets the criteria, and the crack opening is estimated by a Coulomb friction model [22]. FEM is the most mature and the most widely used. FEM is a method based on the continuum mechanics essentially. It must be improved to simulate the discontinuity. The improving methods are divided into two categories: the variable mesh method and the fixed mesh method. In the variable mesh method, the mesh needs updating as the crack tip advances. The cracked surface must be consistent with the edge of the element and a fine mesh or singular element has to be adopted at the crack tip, which could be computationally expensive. In the fixed mesh method, the mesh keeps invariant and the crack is simulated by modifying the interpolation or constitutive relation of the cracked element, such as the XFEM, the smeared crack model, and the continuum damage model. The use of the XFEM for hydraulic fracture problem can avoid remeshing. However, the XFEM also needs objective crack propagation criteria, which are usually based on quantities such as crack-tip stresses and stress intensity factors. So fine crack-tip meshes are necessary for accurate calculation of these quantities. This means that fine meshes are needed if cracks are unknown a priori, leading to high computational cost. The smeared crack model describes a cracked solid by an equivalent anisotropic continuum with degraded material properties in the direction normal to the crack orientation and no remeshing is needed [23]. In contrast, the fixed mesh method is more convenient.Among the numerical studies of hydraulic fracture in concrete gravity dams mentioned above, the XFEM is used mostly while the continuum damage model is less used. In this study, by regarding concrete as a saturated porous medium and employing the effective stress principle of porous media, a stress-seepage-damage coupling model based on FEM is developed and first applied in hydraulic fracture in concrete gravity dams. The coupling model has the following characteristics: (1) the constitutive law considers the strain softening characteristic of fracture process zone, damage-dependent pore-pressure-influence coefficient, stress-dependent permeability for the prepeak stage, and deformation-dependent permeability for the postpeak stage; (2) based on the principle of conservation of fracture energy, the damage model is combined with fracture mechanics to prevent the fracture energy dissipation from being affected by the finite element mesh size; (3) four coupling processes during hydraulic fracture are considered. By this coupling model, hydraulic fracture of Koyna dam is simulated. And the influences of hydraulic fracture on the crack trajectory and the dam bearing capacity are discussed.


## 2. Stress-Seepage-Damage Coupling Model

In this study, the dam concrete is assumed to be a saturated porous medium. In practice, dam concrete can hardly attain a saturated state owing to the small size of the pores and the consequent very low permeability in the intact condition. Except near cracks, the pore pressure of most of the zones of a concrete dam do not vary [24]. The consideration of the dam concrete as a saturated porous medium therefore constitutes a significant simplification. However, the purpose of this study is to investigate the impact of hydraulic fracture on the dam. Isotropic damage models are not appropriate for concrete because the crack trajectory follows principal compressive stresses which are perpendicular to principal tensile stresses promoting crack initiation and propagation. However, if anisotropic damage models are adopted, both theoretical derivation and numerical calculation are difficult. In contrast, isotropic damage models are easier to implement. When considering the effects of hydraulic fracturing, water flows along the crack, permeability coefficients in the directions parallel to the crack plane are different from that perpendicular to the crack plane. And the pore-pressure-influence coefficient in the direction perpendicular to the crack plane is also different from those parallel to the crack plane. Thus, in this study, an isotropic damage model and an anisotropic seepage model are used. The following nonlinear behaviors are defined in this study: (1) the stress-strain relationship of the damaged element, (2) the variation of the pore-pressure-influence coefficient caused by the damage, and (3) the variation of the permeability coefficients due to stress and the damage. In the calculations, it is assumed that the compressive stress-strain relationship is linear elastic because compressive stresses higher than the compressive strength generally do not occur in a gravity dam. The characteristics of the coupling model are discussed below.


### 2.1. Damage Model

At the front of a real concrete crack tip, there is a nonlinear fracture process zone where cohesive stresses can be transferred between the crack interfaces through aggregate interlock and interface friction. The existence of the fracture process zone causes the concrete to exhibit strain softening. The mechanical properties of the fracture process zone can be well simulated using the cohesive crack model proposed by Hillerborg et al. [12]. The majority of researchers have suggested the adoption of a linear stress-strain relationship for the ascending segment under uniaxial tension. However, different modes may occur in the descending segment, including single-line descending, segmented-line descending, and curved descending [25]. Irrespective of the adopted mode, the fracture energy of the stress-strain curve should always remain the same. In this study, the negative exponential equation developed by Jiang et al. [25] is adopted for the descending segment. The stress-strain relationship under uniaxial tension can therefore be expressed as(1)σ=E0εt,εt≤ε0ft0e-αεt-ε0,εt>ε0,where E0 is Young’s modulus of the intact material; ft0 is the tensile strength; εt is the uniaxial tensile strain; ε0 is the strain at cracking (=ft0/E0); and α is the softening coefficient for controlling the descending segment.In Figure1, the fracture energy GF is given by the area enclosed by the stress-crack width curve and the coordinate axes, as expressed by (2). As shown in Figure 2, the fracture energy per unit crack width gf is given by the area enclosed by the descending segment (represented by the line ε=ε0) and the horizontal axis ε, as expressed by (3).(2)GF=∫σdw,w=∫hdε,where w is the sum of the opening displacements of the microcracks in the fracture process zone and h is the width of the distribution area of the microcracks.(3)gf=∫ε0∞σdεt.Figure 1
Stress versus crack width curve and fracture energyGF.Figure 2
Fracture energy per unit crack widthgf.From (1),(4)gf=∫ε0∞σdεt=∫ε0∞ft0e-αεt-ε0dεt=ft0α.The relationship between the fracture energyGF and the fracture energy per unit crack width gf is as follows:(5)GF=gflt,where lt is a geometrical constant which is introduced as a measure of the characteristic length of the element. For a plane element, lt is calculated as the square root of the related area for each integration point. In the case of a solid element, it is calculated as the cube root of the related volume for each integration point.It can be determined from (4) and (5) that the softening coefficient of the descending segment should satisfy (6)α=ft0ltGF.According to the principle of equivalent strain proposed by Lemaitre [26], it can be assumed that the strain generated by the Cauchy stress σ exerted on the damaged material is equal to the strain generated by the effective stress σ~ exerted on the undamaged material. That is,(7)ε=σE~=σ~E0=σ1-dE0.The following is thus obtained:(8)σ=1-dE0ε,where E0 is Young’s modulus of the undamaged material; E~ is Young’s modulus of the damaged material; d is the damage variable, with d=0 corresponding to the undamaged state,d = 1 to the completely damaged state, and 0<d<1 to different degrees of damage.From the perspective of damage mechanics, the nonlinearity of the stress-strain relationships of rock and concrete is due to the formation and propagation of microcracks in the materials through load-induced continuous damage. The consequent brittleness is more obvious under tension. It is thus appropriate to use an elastic damage constitutive model to describe the mechanical properties. It has been confirmed that the results of an elastoplastic damage model do not significantly differ from those of an elastic damage model [27]. The elastic damage constitutive model is therefore adopted in this study.Based on (1) and (8), the equation of the damage evolution under uniaxial tension can be expressed as(9)d=0,εt≤ε01-ε0εt·e-αεt-ε0,εt>ε0.


### 2.2. Effective Stress Principle of Porous Media

The concept of the effective stress as applied to rock and concrete saturated by a single-phase fluid can be regarded as the development of Terzaghi’s effective stress principle for soil. Based on the effective stress principle, Biot first introduced a scalar parameter (known as Biot coefficient) to reflect the effect of the pore pressure on the effective stress [28–30]. In [29], the porosity bi is defined as the ratio of the effective pore area to the cross-sectional area Ai. The effective pore area is defined as the sum of all the micropore areas per unit length in the direction perpendicular to the cross-sectional area Ai. The porosity bi is also the ratio of the pore volume Vp to the bulk volume Vb. The elastic behavior of a porous dam concrete is generally assumed to be homogeneous and isotropic. An isotropic value of bi for which bx = by = bz = b0 is often used, where b0 is the initial Biot coefficient. Fauchet et al. [31] suggested b0 ≈ 0.2 for the elastic analysis of an arch dam. For soils, b0 ≈ 1 is assumed for stress calculations. However, for concrete, the parameter b0 is still lower than unity even in the state of imminent failure [24].Tensile stresses are assigned a positive sign. As shown in Figure3, the total stress exerted on the element surface can be divided into two parts: the external stressbp in equilibrium with the internal water pressurep and the effective stress. The average deformation of the element is related to the effective stress, which can be expressed as(10)σ′=Dε,σ′=σ-bp,where σ′ is the effective stress vector (σ′ = σx′,σy′,σz′,τxy′,τyz′,τzx′T); D is the elasticity matrix; ε=εx,εy,εz,γxy,γyz,γzxT; σ is the total stress vector;p is the pore water pressure; b is the pore-pressure-influence vector, defined as b=bx,by,bz,0,0,0 in three-dimensional analysis.Figure 3
Decomposition of stresses in saturated porous media.


### 2.3. Relationship between Pore-Pressure-Influence Coefficient and Damage

For the undamaged element,bx=by=bz is assumed. If the element is damaged, the influence coefficient perpendicular to the crack orientation would be different from that parallel to the crack orientation. Figure 4 qualitatively illustrates the effects of the cracking on the pore-pressure-influence coefficients, where 1, 2, and 3 represent the directions of the first, second, and third principal stresses, respectively. A local coordinate system is generated with the principal stress directions as the axes. b2 and b3 may be within the same range, namely, b0<b2=b3<b1 [32]. b2>b0 is justified because the diffusion of a crack in a finite area may increase the pore pressure action area paralleling to the crack direction [33]. Based on experimental observations, Bary et al. [6] proposed an anisotropic Biot tensor associated with the damage degree and the pore pressure in each direction. Owing to the lack of experimental data on the effects of cracking on the pore-pressure-influence coefficients parallel to the crack direction, it is assumed in this study that the coefficients remain isotropic after the concrete is damaged, implying b1=b2=b3. In the elastic state, b1=b2=b3=b0, where b0 is the initial Biot coefficient. In the fully damaged state, b1=b2=b3=1, and the evolution of the pore-pressure-influence coefficient with the damage can be expressed as(11)b1=b2=b3=b0+d1-b0.Figure 4
Effects of cracking on pore-pressure-influence coefficients and permeability.


### 2.4. Relationship between Permeability Coefficients and Stress/Damage

Concrete is composed of solid skeleton and microscopic pores. This structural characteristic may cause changes in the microscopic geometry and pore structure when the material is loaded or disturbed, resulting in variation of the porosity and permeability [34]. Porosity change essentially comprises two parts: (1) the change in pore volume caused by structural deformation and (2) changes of the pore structure and size due to the formation, propagation, and penetration of defects such as microcracks (i.e., the pore structure change caused by the material damage) [35]. Many studies have shown that the deterioration of concrete most significantly affects its permeability coefficients. Based on tests conducted on several groups of concrete samples, Picandet et al. [36] developed a coupled equation of the permeability coefficients of the material and its low-degree damage. Souley et al. [37] and Gawin et al. [38] also proposed equations of the permeability coefficients and damage variable, based on the results of experiments. Chatzigeorgiou et al. [39] used a discrete model to describe the coupling relationship between the permeability coefficients and the damage variable. Pijaudier-Cabot et al. [40] established the relationship between progressive damage (from diffused to localized damage distribution) and the permeability coefficients. Further, by introducing a seepage mutation coefficient, Zhao et al. [41] and Zhang and Meng [42] developed seepage-stress coupled equations of concrete during damage evolution under the action of complex stresses.When the element is in the elastic state (i.e.,d=0), the permeability coefficients are exponentially related to the effective stresses as follows [43]:(12)k11=k0ea0σ1k22=k0ea0σ2k33=k0ea0σ3,where k0 is the initial permeability coefficient; a0 is the coupling coefficient; σ1, σ2, and σ3 are the effective principal stresses, which are positive when tensile; and k11, k22, and k33 are the principal permeability coefficients in directions respectively corresponding to the principal stresses.When the first principal stress reaches the tensile strength, the element is damaged (i.e.,d>0) and a fracture is generated in it. Figure 4 qualitatively demonstrates the effect of cracking on the permeability of the element. The water in the concrete mainly flows along the fracture, with the permeability coefficients increasing in directions 2 and 3. Assuming that the fracture is planar and has parallel sides, the aperture of the fracture is given approximately by(13)u=ε1-εtlt,where εt and lt are as defined in Section 2.1. The so-called cubic law gives the flow rate between smooth parallel plates as(14)q=gu312υΔHl,where ΔH is the water head loss across the two ends; υ is the kinematic viscosity coefficient of water, and g is the acceleration of gravity. In (14), the hydraulic conductivity is given by the term gu2/12υ. Therefore, the hydraulic conductivity for a damaged element can be expressed as(15)k11=k0ea0σ1k22=k33=gu212υ=gε1-εtlt212υ.The permeability matrix[k′] in the local coordinate system contains k11, k22, and k33. By transforming this matrix into the global coordinate system through coordinate conversation, the permeability matrix in the global coordinate system [k] can be expressed as(16)k=βTk′β,where k=kxxkxykxzkyxkyykyzkzxkzykzz, k′=k11000k22000k33, and β=l1m1n1l2m2n2l3m3n3.  [β] is the transformation matrix between the two coordinate systems, and li,mi,ni(i=1,2,3) are the respective direction cosines of the principal stress σi.


### 2.5. Basic Differential Equations of Seepage Field

Assuming that the water is incompressible, according to Darcy’s law, the calculation of the steady seepage field with free surface (without internal sources) can be reduced to solving the quasiharmonic equations that satisfy the boundary conditions(17)∂∂xikij∂H∂xj=0HΓ1=ϕx,y,zk∂H∂ηΓ2=qx,y,zHΓ3=Zx,y,k∂H∂ηΓ3=0HΓ4=Zx,y,k∂H∂ηΓ4≤0,where H is the water head; kij is the permeability tensor; and Г1, Г2, Г3, and Г4, respectively, denote the boundaries of the water head, flow, free surface, and overflow.


## 2.1. Damage Model

At the front of a real concrete crack tip, there is a nonlinear fracture process zone where cohesive stresses can be transferred between the crack interfaces through aggregate interlock and interface friction. The existence of the fracture process zone causes the concrete to exhibit strain softening. The mechanical properties of the fracture process zone can be well simulated using the cohesive crack model proposed by Hillerborg et al. [12]. The majority of researchers have suggested the adoption of a linear stress-strain relationship for the ascending segment under uniaxial tension. However, different modes may occur in the descending segment, including single-line descending, segmented-line descending, and curved descending [25]. Irrespective of the adopted mode, the fracture energy of the stress-strain curve should always remain the same. In this study, the negative exponential equation developed by Jiang et al. [25] is adopted for the descending segment. The stress-strain relationship under uniaxial tension can therefore be expressed as(1)σ=E0εt,εt≤ε0ft0e-αεt-ε0,εt>ε0,where E0 is Young’s modulus of the intact material; ft0 is the tensile strength; εt is the uniaxial tensile strain; ε0 is the strain at cracking (=ft0/E0); and α is the softening coefficient for controlling the descending segment.In Figure1, the fracture energy GF is given by the area enclosed by the stress-crack width curve and the coordinate axes, as expressed by (2). As shown in Figure 2, the fracture energy per unit crack width gf is given by the area enclosed by the descending segment (represented by the line ε=ε0) and the horizontal axis ε, as expressed by (3).(2)GF=∫σdw,w=∫hdε,where w is the sum of the opening displacements of the microcracks in the fracture process zone and h is the width of the distribution area of the microcracks.(3)gf=∫ε0∞σdεt.Figure 1
Stress versus crack width curve and fracture energyGF.Figure 2
Fracture energy per unit crack widthgf.From (1),(4)gf=∫ε0∞σdεt=∫ε0∞ft0e-αεt-ε0dεt=ft0α.The relationship between the fracture energyGF and the fracture energy per unit crack width gf is as follows:(5)GF=gflt,where lt is a geometrical constant which is introduced as a measure of the characteristic length of the element. For a plane element, lt is calculated as the square root of the related area for each integration point. In the case of a solid element, it is calculated as the cube root of the related volume for each integration point.It can be determined from (4) and (5) that the softening coefficient of the descending segment should satisfy (6)α=ft0ltGF.According to the principle of equivalent strain proposed by Lemaitre [26], it can be assumed that the strain generated by the Cauchy stress σ exerted on the damaged material is equal to the strain generated by the effective stress σ~ exerted on the undamaged material. That is,(7)ε=σE~=σ~E0=σ1-dE0.The following is thus obtained:(8)σ=1-dE0ε,where E0 is Young’s modulus of the undamaged material; E~ is Young’s modulus of the damaged material; d is the damage variable, with d=0 corresponding to the undamaged state,d = 1 to the completely damaged state, and 0<d<1 to different degrees of damage.From the perspective of damage mechanics, the nonlinearity of the stress-strain relationships of rock and concrete is due to the formation and propagation of microcracks in the materials through load-induced continuous damage. The consequent brittleness is more obvious under tension. It is thus appropriate to use an elastic damage constitutive model to describe the mechanical properties. It has been confirmed that the results of an elastoplastic damage model do not significantly differ from those of an elastic damage model [27]. The elastic damage constitutive model is therefore adopted in this study.Based on (1) and (8), the equation of the damage evolution under uniaxial tension can be expressed as(9)d=0,εt≤ε01-ε0εt·e-αεt-ε0,εt>ε0.


## 2.2. Effective Stress Principle of Porous Media

The concept of the effective stress as applied to rock and concrete saturated by a single-phase fluid can be regarded as the development of Terzaghi’s effective stress principle for soil. Based on the effective stress principle, Biot first introduced a scalar parameter (known as Biot coefficient) to reflect the effect of the pore pressure on the effective stress [28–30]. In [29], the porosity bi is defined as the ratio of the effective pore area to the cross-sectional area Ai. The effective pore area is defined as the sum of all the micropore areas per unit length in the direction perpendicular to the cross-sectional area Ai. The porosity bi is also the ratio of the pore volume Vp to the bulk volume Vb. The elastic behavior of a porous dam concrete is generally assumed to be homogeneous and isotropic. An isotropic value of bi for which bx = by = bz = b0 is often used, where b0 is the initial Biot coefficient. Fauchet et al. [31] suggested b0 ≈ 0.2 for the elastic analysis of an arch dam. For soils, b0 ≈ 1 is assumed for stress calculations. However, for concrete, the parameter b0 is still lower than unity even in the state of imminent failure [24].Tensile stresses are assigned a positive sign. As shown in Figure3, the total stress exerted on the element surface can be divided into two parts: the external stressbp in equilibrium with the internal water pressurep and the effective stress. The average deformation of the element is related to the effective stress, which can be expressed as(10)σ′=Dε,σ′=σ-bp,where σ′ is the effective stress vector (σ′ = σx′,σy′,σz′,τxy′,τyz′,τzx′T); D is the elasticity matrix; ε=εx,εy,εz,γxy,γyz,γzxT; σ is the total stress vector;p is the pore water pressure; b is the pore-pressure-influence vector, defined as b=bx,by,bz,0,0,0 in three-dimensional analysis.Figure 3
Decomposition of stresses in saturated porous media.


## 2.3. Relationship between Pore-Pressure-Influence Coefficient and Damage

For the undamaged element,bx=by=bz is assumed. If the element is damaged, the influence coefficient perpendicular to the crack orientation would be different from that parallel to the crack orientation. Figure 4 qualitatively illustrates the effects of the cracking on the pore-pressure-influence coefficients, where 1, 2, and 3 represent the directions of the first, second, and third principal stresses, respectively. A local coordinate system is generated with the principal stress directions as the axes. b2 and b3 may be within the same range, namely, b0<b2=b3<b1 [32]. b2>b0 is justified because the diffusion of a crack in a finite area may increase the pore pressure action area paralleling to the crack direction [33]. Based on experimental observations, Bary et al. [6] proposed an anisotropic Biot tensor associated with the damage degree and the pore pressure in each direction. Owing to the lack of experimental data on the effects of cracking on the pore-pressure-influence coefficients parallel to the crack direction, it is assumed in this study that the coefficients remain isotropic after the concrete is damaged, implying b1=b2=b3. In the elastic state, b1=b2=b3=b0, where b0 is the initial Biot coefficient. In the fully damaged state, b1=b2=b3=1, and the evolution of the pore-pressure-influence coefficient with the damage can be expressed as(11)b1=b2=b3=b0+d1-b0.Figure 4
Effects of cracking on pore-pressure-influence coefficients and permeability.


## 2.4. Relationship between Permeability Coefficients and Stress/Damage

Concrete is composed of solid skeleton and microscopic pores. This structural characteristic may cause changes in the microscopic geometry and pore structure when the material is loaded or disturbed, resulting in variation of the porosity and permeability [34]. Porosity change essentially comprises two parts: (1) the change in pore volume caused by structural deformation and (2) changes of the pore structure and size due to the formation, propagation, and penetration of defects such as microcracks (i.e., the pore structure change caused by the material damage) [35]. Many studies have shown that the deterioration of concrete most significantly affects its permeability coefficients. Based on tests conducted on several groups of concrete samples, Picandet et al. [36] developed a coupled equation of the permeability coefficients of the material and its low-degree damage. Souley et al. [37] and Gawin et al. [38] also proposed equations of the permeability coefficients and damage variable, based on the results of experiments. Chatzigeorgiou et al. [39] used a discrete model to describe the coupling relationship between the permeability coefficients and the damage variable. Pijaudier-Cabot et al. [40] established the relationship between progressive damage (from diffused to localized damage distribution) and the permeability coefficients. Further, by introducing a seepage mutation coefficient, Zhao et al. [41] and Zhang and Meng [42] developed seepage-stress coupled equations of concrete during damage evolution under the action of complex stresses.When the element is in the elastic state (i.e.,d=0), the permeability coefficients are exponentially related to the effective stresses as follows [43]:(12)k11=k0ea0σ1k22=k0ea0σ2k33=k0ea0σ3,where k0 is the initial permeability coefficient; a0 is the coupling coefficient; σ1, σ2, and σ3 are the effective principal stresses, which are positive when tensile; and k11, k22, and k33 are the principal permeability coefficients in directions respectively corresponding to the principal stresses.When the first principal stress reaches the tensile strength, the element is damaged (i.e.,d>0) and a fracture is generated in it. Figure 4 qualitatively demonstrates the effect of cracking on the permeability of the element. The water in the concrete mainly flows along the fracture, with the permeability coefficients increasing in directions 2 and 3. Assuming that the fracture is planar and has parallel sides, the aperture of the fracture is given approximately by(13)u=ε1-εtlt,where εt and lt are as defined in Section 2.1. The so-called cubic law gives the flow rate between smooth parallel plates as(14)q=gu312υΔHl,where ΔH is the water head loss across the two ends; υ is the kinematic viscosity coefficient of water, and g is the acceleration of gravity. In (14), the hydraulic conductivity is given by the term gu2/12υ. Therefore, the hydraulic conductivity for a damaged element can be expressed as(15)k11=k0ea0σ1k22=k33=gu212υ=gε1-εtlt212υ.The permeability matrix[k′] in the local coordinate system contains k11, k22, and k33. By transforming this matrix into the global coordinate system through coordinate conversation, the permeability matrix in the global coordinate system [k] can be expressed as(16)k=βTk′β,where k=kxxkxykxzkyxkyykyzkzxkzykzz, k′=k11000k22000k33, and β=l1m1n1l2m2n2l3m3n3.  [β] is the transformation matrix between the two coordinate systems, and li,mi,ni(i=1,2,3) are the respective direction cosines of the principal stress σi.


## 2.5. Basic Differential Equations of Seepage Field

Assuming that the water is incompressible, according to Darcy’s law, the calculation of the steady seepage field with free surface (without internal sources) can be reduced to solving the quasiharmonic equations that satisfy the boundary conditions(17)∂∂xikij∂H∂xj=0HΓ1=ϕx,y,zk∂H∂ηΓ2=qx,y,zHΓ3=Zx,y,k∂H∂ηΓ3=0HΓ4=Zx,y,k∂H∂ηΓ4≤0,where H is the water head; kij is the permeability tensor; and Г1, Г2, Г3, and Г4, respectively, denote the boundaries of the water head, flow, free surface, and overflow.


## 3. Realization of Finite Element Method

In this study, the cracking problem is formulated in effective stresses, namely, total stresses minus the pore pressures. And the effective stresses are used to perform the computation. In the finite element analysis, the equilibrium between the external load and the internal load is as follows:(18)∭ΩBTσ′dV=F-Fp,where [B] is the stress-displacement transformation matrix, σ′ is the effective stress vector, and F and Fp are the external load and internal pore pressure load, respectively.(19)Fp=∭ΩBTbpdV.The seepage, stress, and damage fields mutually affect and are coupled to each other, and solving the problem using the fully coupled equations would require a significant amount of calculation. The weak coupled method is used to solve each equation independently. In the solution of the seepage field of theNth load step, the permeability coefficients modified by the stress field and the damage field of the (N-1)th load step are used to form the new permeability stiffness matrix, and the seepage field of theNth load step is obtained. The displacement field, stress field, and damage field of theNth load step are then obtained. Figure 5 shows the flowchart of the stress-seepage-damage coupling program.Figure 5
Flowchart of coupling program.During the calculation, the loads are applied by the incremental method and the load increment is set to be quite small. The total quantity method is used for the iteration in each load step. The damage degree could increase or decrease during an iteration process but could not reduce below that corresponding to the last load step. It is therefore necessary to preserve the damage degree of the last load step. The load of a current iteration step is the sum of the load increment of the current load step and the equivalent load increment of the unbalanced stress caused by the damage in the last iteration step. Further, the equivalent load increment of the unbalanced stress is the load increment due to the difference between the computational stress and the bearing stress. It should be noted that the load of the first iteration step is the sum of the external load increment of the current load step and the residual unbalanced force of the last load step. The equivalent load increment of the unbalanced stress and the residual unbalanced force can be calculated using (20) and (21), respectively.(20)Δfistepi=∭ΩBTσistep-1-σistepidV,istep≥1,i≥1(21)Δfistep=Δfistepj-Δfistepj-1,assuming  that  the  calculation  is  convergent  in  thejth  iteration  step,where Δfistepi is the equivalent load increment of the ith iteration step of the (i step)th load step, required to be added to the load increment of the (i+1)th iteration step if the calculation remained divergent in the ith iteration step; and Δfistep is the residual unbalanced force of the (i step)th load step, required to be added to the load increment of the first iteration step of the (i step + 1)th load step; that is,(22)Δfistep0=Δfistep-1,istep≥1.When the iteration converges, the residual unbalanced force is generally very small.The permeability stiffness matrixKsistep of each load step is derived from the permeability coefficients modified in accordance with the stress and damage states of the last load step. Wheni step = 1, Ks1 is derived from the initial permeability coefficients. The stiffness matrix of the first iteration step Kistep1 of each load step is derived from the elasticity coefficient modified based on the damage state of the last load step. Wheni step = 1, K11 is derived from the initial elasticity coefficient. The stiffness matrix of the other iteration steps Kistepi(i>1) in each load step is obtained from the elasticity coefficients modified in accordance with the damage state of the (i-1)th iteration step.Let us assume that I is the seepage module, II is the damage iteration module, and III is the coupling module. In this study, the computations are performed in a weakly coupled form, in a sequential manner by applying I-II-III without iteration. However, there are iterations in the damage iteration module. To guarantee the convergence of the results, a small load step should be applied. The smaller the load step is, the more accurate the results are. When the load step is small enough, the results of the weakly coupled algorithm can be regarded as those of the strongly coupled algorithm (applying I-II-III with iteration).The coupling program can also be used to investigate the crack propagation without consideration of the stress-seepage-damage coupling effect. By inputting the control parameters, the seepage and coupling modules can be skipped, and the damage iteration module directly calculated.


## 4. Case Study

### 4.1. Crack Propagation Analysis of Model Gravity Dam

Carpinteri et al. [44] once conducted fracture tests on a 1 : 40 scaled gravity dam subjected to equivalent hydraulic loads in the structural and material laboratory of University of Turin in Italy. The model was cast with conventional concrete mix and was 2.4 m high, 2.0 m wide at the bottom, and 0.3 m thick. A notch of length 0.2W (30 cm) was fabricated on its upstream face, at a quarter of its height, where W is the dam thickness at the location of the notch. The dimensions of the models are shown in Figure 6. The applied force was distributed as four concentrated loads with intensities as shown in Figure 6(a). The force was gradually increased until failure of the dam. The mechanical properties of the utilized concrete are listed in Table 1.Table 1
Material parameters of gravity dam models.
Young’s modulus,E (GPa)
Poisson’s ratio,ν
Density,ρ (kg/m3)
Tensile strength,ft  (MPa)
Fracture energy,Gf (N/m)
35.7
0.1
2400
3.6
184Figure 6
Dimensions and load distribution of gravity dam models (unit: mm).
(a)
(b)
(c)Various researchers have performed numerical simulations of the above tests using different models and approaches such as a combination of FEM and the cohesive crack model of Barpi and Valente [45], the XFEM of Du et al. [46]. In the present study, the experimental model is analyzed by the damage model, and the predicted responses are compared with the experimental results and the results of Barpi and Valente [45] and Du et al. [46].Figure7 shows the relationship between the total applied force and the CMOD, as obtained by the present damage model, the cohesive crack model of Barpi and Valente [45], and the XFEM of Du et al. [46], and the experimental results. As can be seen from Figure 7, the peak load predicted by the present damage model is in good agreement with those determined by the experiment and Barpi and Valente [45]; that predicted by the XFEM is relatively low. It can also be seen from the figure that the predictions of the present damage model indicate an initial stiffness, which agrees well with the experimental results. The results of Barpi and Valente [45] and Du et al. [46] indicated a softer prepeak response and larger CMOD at the peak load. There is a little difference between the postpeak response of this study and that of the experiment.Figure 7
Load versus CMOD.Figure8 shows the damage distribution determined in the present study. Figure 9 compares the numerically and experimentally determined crack propagation paths. The experimentally observed crack trajectory is initially horizontal and then turns toward the dam toe. The crack trajectories obtained by Barpi and Valente [45] and Du et al. [46] are initially downward and then begin to propagate with a nearly constant slope, representing significant deviation from the experimentally observed crack profile. The damage distribution determined in this study also agrees well with the experimentally observed trajectory, further substantiating the validity of the damage model for predicting the crack trajectory. The discrete crack model and XFEM thus do not appear to be significantly better than the present damage model for predicting crack trajectory in a case study.Figure 8
Damage distribution.Figure 9
Comparison of crack propagation paths.


### 4.2. Crack Propagation Analysis of Practical Gravity Dam (Koyna Dam)

In this section, the crack propagation processes of Koyna dam under overflow are simulated by the damage model and the coupling model, respectively. Koyna gravity dam was severely damaged near the dam neck by a strong earthquake [47] and has become a typical reference point for static and dynamic loading of a dam over the last few decades. Gioia et al. [48] used linear elastic fracture mechanics to investigate the fracture response of a dam with a preset notch on its upstream face and subjected to overflow, while Bhattacharjee and Léger [49] employed a smeared crack model for the same purpose. The nonlinear response of Koyna dam under overflow is investigated in the present study. To facilitate comparison of the results with those of previous studies, the employed geometry and configuration are the same as those used by Bhattacharjee and Léger [49]. The employed geometry and finite element mesh are shown in Figure 10. A prefabricated horizontal notch is set on the upstream face at the elevation of the downstream slope variation. The initial depth of the notch is here denoted bya, and the dam thickness at the location of the initial notch by l (l=19.3 mm). A plane stress state is assumed in the study. The considered loads include the self-weight of the dam, hydrostatic pressure of the full dam, and increment of the water pressure due to reservoir overflow. The bottom of the dam is rigidly supported. Table 2 gives the mechanical properties of the dam concrete.Table 2
Mechanical properties of Koyna dam concrete.
Young’s modulus,E (GPa)
Poisson’s ratio,ν
Density,ρ (kg/m3)
Tensile strength,ft (MPa)
Fracture energy,Gf (N/m)
25
0.2
2450
1.0
100Figure 10
Geometry and finite element mesh of Koyna dam (unit: mm).


#### 4.2.1. Crack Propagation Analysis without Consideration of Coupling Effect

(
1) Initial Notch Length
a
=
0.1
l. The relationships between the overflow and the horizontal displacement of the dam crest determined by the present damage model, the smeared crack model of Bhattacharjee and Léger [49], and Gioia et al. [48] using LEFM are shown in Figure 11. The initial displacement of the dam crest corresponds to the combined load of the self-weight and the full reservoir pressure without overflow. There is consistency among the displacements obtained by the three methods before the obvious occurrence of nonlinearity. In the nonlinear stage, the results of the present study agree well with those of Bhattacharjee and Léger [49]. The ultimate overflows determined by Bhattacharjee and Léger [49], in the present study, and by Gioia et al. [48] are approximately 10.0, 10.2, and 14.0 m, respectively. The structural resistances predicted by the smeared crack model and the damage model appear to be slightly lower than those determined by LEFM analysis.Figure 11
Overflow versus horizontal displacement of dam crest.The presently determined damage distribution for an overflow of 10.2 m is shown in Figure12. The distribution is compared with the results of Gioia et al. [48] and Bhattacharjee and Léger [49] in Figure 13. According to the present numerical simulation, the crack initially propagates horizontally. With increasing overflow, the crack gradually turns downward due to the increase in compressive stress downstream of the dam. The crack propagation path obtained by Gioia et al. [48] actually corresponds to an overflow of about 14 m, which explains its longer path in Figure 13. Nevertheless, its consistency with that of the present study is apparent. Both paths initially run horizontally over a very short distance and then turn downward. The crack growth path obtained by Bhattacharjee and Léger [49] corresponds to an overflow of about 10.0 m. It initially runs nearly horizontally before turning downward when the crack length reaches half the width of the dam neck. The damage zone is in good agreement with the results of Gioia et al. [48].Figure 12
Damage distribution (overflow = 10.2 m).Figure 13
Comparison of crack trajectories obtained by different methods.(
2) Different Initial Notch Lengths. To examine the effect of the initial notch length a on the carrying capacity and crack path, three initial lengths are, respectively, used for analysis, namely, a=0.1l, 0.2l, and 0.3l. The relationships between the overflow and the horizontal displacement of the dam crest for the three cases are shown in Figure 14. The corresponding damage distributions for an overflow of 10.2 m are shown in Figure 15.Figure 14
Overflow versus horizontal displacement of dam crest.Figure 15
Damage distributions for overflow of 10.2 m. (a)a= 0.1l, (b)a= 0.2l, and (c)a= 0.3l.
(a)
(b)
(c)It can be seen from Figure10 that (1) the critical overflow at which the dam is transformed from the linear state into the nonlinear state differs with the length of the initial notch, being 8.6, 7.4, and 7.0 m fora= 0.1l, 0.2l, and 0.3l, respectively, and (2) the ultimate overflows for the various values ofa are extremely close, being about 10.2 m. It can also be seen from Figure 15 that the cracks initially extend in the horizontal direction and then turn downward at nearly the same point and with almost the same deflection angle. The crack paths are also almost the same. The length of the initial notch thus only affects the critical overflow for transformation into the nonlinear state. The shorter the length, the higher the critical overflow. The notch length has no obvious effect on the ultimate carrying capacity of the dam and the crack path.


#### 4.2.2. Crack Propagation considering Coupling Effect

To investigate the impact of the coupling effect during crack growth on the crack propagation path and carrying capacity of the dam, nonlinear stress-seepage-damage coupling analyses are conducted on the dam using the three initial notch lengths in Section4.2.1(2). The initial permeability coefficient k0 is set as 5 × 10−9 m/s and the coupling coefficienta0 is set as 0.01. It is really true that the dam which is 100 m high is constructed in reality with drains. Due to lack of information about the drains of Koyna dam, the drains are not taken into account in this study.In this section, in order to select an appropriate load step, five different load steps are chosen, 1.0 m, 0.5 m, 0.1 m, 0.05 m, and 0.01 m, respectively. Taking Koyna dam with an initial notcha = 0.1l as an example, the ultimate overflows under the five load steps are 8.0 m, 7.5 m, 7.2 m, 7.0 m, and 6.98 m, respectively. Assuming the result of the load step 0.01 m as the precise solution, the errors of the load steps 1.0 m, 0.5 m, 0.1 m, and 0.05 m are 14.61%, 7.45%, 3.15%, and 0.29%, respectively. It can be seen that when the load step is taken as 0.05 m, the value of the load step has little influence on the solution. Thus, the result of the load step 0.05 m is convergent and the load step 0.05 m can be chosen in the calculation of this section.(
1) Initial Notch Length
a
=
0.1
l. The relationships between the overflow and the horizontal displacement for an initial notch lengthaof 0.1l with and without consideration of the coupling effect are shown in Figure 16. The damage distribution with consideration of the coupling effect is shown in Figure 17. It can be seen from Figure 16 that, before the occurrence of obvious nonlinearity, the results for the two cases are in agreement. When the coupling effect is taken into consideration, the critical overflow that transforms the dam from the linear state into the nonlinear state is 4.1 m, and the ultimate overflow is 7.0 m. This represents significant decreases in both parameters compared to when not considering the effect. Further, from a comparison of Figures 15(a) and 17, it can be seen that consideration of the coupling effect decreases the horizontal crack propagation distance, while increasing the deflection angle of the downward turn and the total length of the path.Figure 16
Overflow versus horizontal displacement of dam crest.Figure 17
Damage distribution for overflow of 7.0 m.(
2) Different Initial Notch Lengths. Figure 18 shows the damage distributions for the three initial notch lengths corresponding to an overflow of 4.2 m with consideration of the coupling effect. The relationships between the overflow and the horizontal displacement of the dam crest for the three notch lengths are shown in Figure 19. The critical and ultimate overflows with and without consideration of the coupling effect are given in Table 3.Table 3
Comparison of dam carrying capacities with and without consideration of coupling effect.
Initial crack length
Without consideration of coupling effect
With consideration of coupling effect
Percentage change
Critical load (m)
Carrying capacity (m)
Critical load (m)
Carrying capacity (m)
Critical load
Carrying capacity
0.1l
6.0
10.2
4.1
7.0
−31.7%
−31.4%
02l
5.0
10.2
2.0
5.3
−60.0%
−48.0%
0.3l
3.0
10.2
1.0
4.2
−66.7%
−58.8%Figure 18
Damage distributions for overflow of 4.2 m. (a)a= 0.1l, (b)a= 0.2l, and (c)a= 0.3l.
(a)
(b)
(c)Figure 19
Overflow versus horizontal displacement of dam crest.It can be seen from Figure18 that, for a given overflow, the crack only propagates over a short distance in the horizontal direction whena= 0.1l, while it initially propagates in the horizontal direction and then turns downward whena= 0.2l. Whena= 0.3l, the crack turns downward almost directly and the propagation distance is nearly twice that fora= 0.2l. With consideration of the coupling effect, the crack propagation distance increases with increasing initial notch length. It can also be seen from Figure 19 that, with consideration of the coupling effect, the critical and ultimate overflows decrease with increasing initial notch length. Further, Table 3 indicates that, with consideration of the coupling effect, the rates of reduction of the critical and ultimate overflows increase with increasing initial notch length.


## 4.1. Crack Propagation Analysis of Model Gravity Dam

Carpinteri et al. [44] once conducted fracture tests on a 1 : 40 scaled gravity dam subjected to equivalent hydraulic loads in the structural and material laboratory of University of Turin in Italy. The model was cast with conventional concrete mix and was 2.4 m high, 2.0 m wide at the bottom, and 0.3 m thick. A notch of length 0.2W (30 cm) was fabricated on its upstream face, at a quarter of its height, where W is the dam thickness at the location of the notch. The dimensions of the models are shown in Figure 6. The applied force was distributed as four concentrated loads with intensities as shown in Figure 6(a). The force was gradually increased until failure of the dam. The mechanical properties of the utilized concrete are listed in Table 1.Table 1
Material parameters of gravity dam models.
Young’s modulus,E (GPa)
Poisson’s ratio,ν
Density,ρ (kg/m3)
Tensile strength,ft  (MPa)
Fracture energy,Gf (N/m)
35.7
0.1
2400
3.6
184Figure 6
Dimensions and load distribution of gravity dam models (unit: mm).
(a)
(b)
(c)Various researchers have performed numerical simulations of the above tests using different models and approaches such as a combination of FEM and the cohesive crack model of Barpi and Valente [45], the XFEM of Du et al. [46]. In the present study, the experimental model is analyzed by the damage model, and the predicted responses are compared with the experimental results and the results of Barpi and Valente [45] and Du et al. [46].Figure7 shows the relationship between the total applied force and the CMOD, as obtained by the present damage model, the cohesive crack model of Barpi and Valente [45], and the XFEM of Du et al. [46], and the experimental results. As can be seen from Figure 7, the peak load predicted by the present damage model is in good agreement with those determined by the experiment and Barpi and Valente [45]; that predicted by the XFEM is relatively low. It can also be seen from the figure that the predictions of the present damage model indicate an initial stiffness, which agrees well with the experimental results. The results of Barpi and Valente [45] and Du et al. [46] indicated a softer prepeak response and larger CMOD at the peak load. There is a little difference between the postpeak response of this study and that of the experiment.Figure 7
Load versus CMOD.Figure8 shows the damage distribution determined in the present study. Figure 9 compares the numerically and experimentally determined crack propagation paths. The experimentally observed crack trajectory is initially horizontal and then turns toward the dam toe. The crack trajectories obtained by Barpi and Valente [45] and Du et al. [46] are initially downward and then begin to propagate with a nearly constant slope, representing significant deviation from the experimentally observed crack profile. The damage distribution determined in this study also agrees well with the experimentally observed trajectory, further substantiating the validity of the damage model for predicting the crack trajectory. The discrete crack model and XFEM thus do not appear to be significantly better than the present damage model for predicting crack trajectory in a case study.Figure 8
Damage distribution.Figure 9
Comparison of crack propagation paths.


## 4.2. Crack Propagation Analysis of Practical Gravity Dam (Koyna Dam)

In this section, the crack propagation processes of Koyna dam under overflow are simulated by the damage model and the coupling model, respectively. Koyna gravity dam was severely damaged near the dam neck by a strong earthquake [47] and has become a typical reference point for static and dynamic loading of a dam over the last few decades. Gioia et al. [48] used linear elastic fracture mechanics to investigate the fracture response of a dam with a preset notch on its upstream face and subjected to overflow, while Bhattacharjee and Léger [49] employed a smeared crack model for the same purpose. The nonlinear response of Koyna dam under overflow is investigated in the present study. To facilitate comparison of the results with those of previous studies, the employed geometry and configuration are the same as those used by Bhattacharjee and Léger [49]. The employed geometry and finite element mesh are shown in Figure 10. A prefabricated horizontal notch is set on the upstream face at the elevation of the downstream slope variation. The initial depth of the notch is here denoted bya, and the dam thickness at the location of the initial notch by l (l=19.3 mm). A plane stress state is assumed in the study. The considered loads include the self-weight of the dam, hydrostatic pressure of the full dam, and increment of the water pressure due to reservoir overflow. The bottom of the dam is rigidly supported. Table 2 gives the mechanical properties of the dam concrete.Table 2
Mechanical properties of Koyna dam concrete.
Young’s modulus,E (GPa)
Poisson’s ratio,ν
Density,ρ (kg/m3)
Tensile strength,ft (MPa)
Fracture energy,Gf (N/m)
25
0.2
2450
1.0
100Figure 10
Geometry and finite element mesh of Koyna dam (unit: mm).


### 4.2.1. Crack Propagation Analysis without Consideration of Coupling Effect

(
1) Initial Notch Length
a
=
0.1
l. The relationships between the overflow and the horizontal displacement of the dam crest determined by the present damage model, the smeared crack model of Bhattacharjee and Léger [49], and Gioia et al. [48] using LEFM are shown in Figure 11. The initial displacement of the dam crest corresponds to the combined load of the self-weight and the full reservoir pressure without overflow. There is consistency among the displacements obtained by the three methods before the obvious occurrence of nonlinearity. In the nonlinear stage, the results of the present study agree well with those of Bhattacharjee and Léger [49]. The ultimate overflows determined by Bhattacharjee and Léger [49], in the present study, and by Gioia et al. [48] are approximately 10.0, 10.2, and 14.0 m, respectively. The structural resistances predicted by the smeared crack model and the damage model appear to be slightly lower than those determined by LEFM analysis.Figure 11
Overflow versus horizontal displacement of dam crest.The presently determined damage distribution for an overflow of 10.2 m is shown in Figure12. The distribution is compared with the results of Gioia et al. [48] and Bhattacharjee and Léger [49] in Figure 13. According to the present numerical simulation, the crack initially propagates horizontally. With increasing overflow, the crack gradually turns downward due to the increase in compressive stress downstream of the dam. The crack propagation path obtained by Gioia et al. [48] actually corresponds to an overflow of about 14 m, which explains its longer path in Figure 13. Nevertheless, its consistency with that of the present study is apparent. Both paths initially run horizontally over a very short distance and then turn downward. The crack growth path obtained by Bhattacharjee and Léger [49] corresponds to an overflow of about 10.0 m. It initially runs nearly horizontally before turning downward when the crack length reaches half the width of the dam neck. The damage zone is in good agreement with the results of Gioia et al. [48].Figure 12
Damage distribution (overflow = 10.2 m).Figure 13
Comparison of crack trajectories obtained by different methods.(
2) Different Initial Notch Lengths. To examine the effect of the initial notch length a on the carrying capacity and crack path, three initial lengths are, respectively, used for analysis, namely, a=0.1l, 0.2l, and 0.3l. The relationships between the overflow and the horizontal displacement of the dam crest for the three cases are shown in Figure 14. The corresponding damage distributions for an overflow of 10.2 m are shown in Figure 15.Figure 14
Overflow versus horizontal displacement of dam crest.Figure 15
Damage distributions for overflow of 10.2 m. (a)a= 0.1l, (b)a= 0.2l, and (c)a= 0.3l.
(a)
(b)
(c)It can be seen from Figure10 that (1) the critical overflow at which the dam is transformed from the linear state into the nonlinear state differs with the length of the initial notch, being 8.6, 7.4, and 7.0 m fora= 0.1l, 0.2l, and 0.3l, respectively, and (2) the ultimate overflows for the various values ofa are extremely close, being about 10.2 m. It can also be seen from Figure 15 that the cracks initially extend in the horizontal direction and then turn downward at nearly the same point and with almost the same deflection angle. The crack paths are also almost the same. The length of the initial notch thus only affects the critical overflow for transformation into the nonlinear state. The shorter the length, the higher the critical overflow. The notch length has no obvious effect on the ultimate carrying capacity of the dam and the crack path.


### 4.2.2. Crack Propagation considering Coupling Effect

To investigate the impact of the coupling effect during crack growth on the crack propagation path and carrying capacity of the dam, nonlinear stress-seepage-damage coupling analyses are conducted on the dam using the three initial notch lengths in Section4.2.1(2). The initial permeability coefficient k0 is set as 5 × 10−9 m/s and the coupling coefficienta0 is set as 0.01. It is really true that the dam which is 100 m high is constructed in reality with drains. Due to lack of information about the drains of Koyna dam, the drains are not taken into account in this study.In this section, in order to select an appropriate load step, five different load steps are chosen, 1.0 m, 0.5 m, 0.1 m, 0.05 m, and 0.01 m, respectively. Taking Koyna dam with an initial notcha = 0.1l as an example, the ultimate overflows under the five load steps are 8.0 m, 7.5 m, 7.2 m, 7.0 m, and 6.98 m, respectively. Assuming the result of the load step 0.01 m as the precise solution, the errors of the load steps 1.0 m, 0.5 m, 0.1 m, and 0.05 m are 14.61%, 7.45%, 3.15%, and 0.29%, respectively. It can be seen that when the load step is taken as 0.05 m, the value of the load step has little influence on the solution. Thus, the result of the load step 0.05 m is convergent and the load step 0.05 m can be chosen in the calculation of this section.(
1) Initial Notch Length
a
=
0.1
l. The relationships between the overflow and the horizontal displacement for an initial notch lengthaof 0.1l with and without consideration of the coupling effect are shown in Figure 16. The damage distribution with consideration of the coupling effect is shown in Figure 17. It can be seen from Figure 16 that, before the occurrence of obvious nonlinearity, the results for the two cases are in agreement. When the coupling effect is taken into consideration, the critical overflow that transforms the dam from the linear state into the nonlinear state is 4.1 m, and the ultimate overflow is 7.0 m. This represents significant decreases in both parameters compared to when not considering the effect. Further, from a comparison of Figures 15(a) and 17, it can be seen that consideration of the coupling effect decreases the horizontal crack propagation distance, while increasing the deflection angle of the downward turn and the total length of the path.Figure 16
Overflow versus horizontal displacement of dam crest.Figure 17
Damage distribution for overflow of 7.0 m.(
2) Different Initial Notch Lengths. Figure 18 shows the damage distributions for the three initial notch lengths corresponding to an overflow of 4.2 m with consideration of the coupling effect. The relationships between the overflow and the horizontal displacement of the dam crest for the three notch lengths are shown in Figure 19. The critical and ultimate overflows with and without consideration of the coupling effect are given in Table 3.Table 3
Comparison of dam carrying capacities with and without consideration of coupling effect.
Initial crack length
Without consideration of coupling effect
With consideration of coupling effect
Percentage change
Critical load (m)
Carrying capacity (m)
Critical load (m)
Carrying capacity (m)
Critical load
Carrying capacity
0.1l
6.0
10.2
4.1
7.0
−31.7%
−31.4%
02l
5.0
10.2
2.0
5.3
−60.0%
−48.0%
0.3l
3.0
10.2
1.0
4.2
−66.7%
−58.8%Figure 18
Damage distributions for overflow of 4.2 m. (a)a= 0.1l, (b)a= 0.2l, and (c)a= 0.3l.
(a)
(b)
(c)Figure 19
Overflow versus horizontal displacement of dam crest.It can be seen from Figure18 that, for a given overflow, the crack only propagates over a short distance in the horizontal direction whena= 0.1l, while it initially propagates in the horizontal direction and then turns downward whena= 0.2l. Whena= 0.3l, the crack turns downward almost directly and the propagation distance is nearly twice that fora= 0.2l. With consideration of the coupling effect, the crack propagation distance increases with increasing initial notch length. It can also be seen from Figure 19 that, with consideration of the coupling effect, the critical and ultimate overflows decrease with increasing initial notch length. Further, Table 3 indicates that, with consideration of the coupling effect, the rates of reduction of the critical and ultimate overflows increase with increasing initial notch length.


## 4.2.1. Crack Propagation Analysis without Consideration of Coupling Effect

(
1) Initial Notch Length
a
=
0.1
l. The relationships between the overflow and the horizontal displacement of the dam crest determined by the present damage model, the smeared crack model of Bhattacharjee and Léger [49], and Gioia et al. [48] using LEFM are shown in Figure 11. The initial displacement of the dam crest corresponds to the combined load of the self-weight and the full reservoir pressure without overflow. There is consistency among the displacements obtained by the three methods before the obvious occurrence of nonlinearity. In the nonlinear stage, the results of the present study agree well with those of Bhattacharjee and Léger [49]. The ultimate overflows determined by Bhattacharjee and Léger [49], in the present study, and by Gioia et al. [48] are approximately 10.0, 10.2, and 14.0 m, respectively. The structural resistances predicted by the smeared crack model and the damage model appear to be slightly lower than those determined by LEFM analysis.Figure 11
Overflow versus horizontal displacement of dam crest.The presently determined damage distribution for an overflow of 10.2 m is shown in Figure12. The distribution is compared with the results of Gioia et al. [48] and Bhattacharjee and Léger [49] in Figure 13. According to the present numerical simulation, the crack initially propagates horizontally. With increasing overflow, the crack gradually turns downward due to the increase in compressive stress downstream of the dam. The crack propagation path obtained by Gioia et al. [48] actually corresponds to an overflow of about 14 m, which explains its longer path in Figure 13. Nevertheless, its consistency with that of the present study is apparent. Both paths initially run horizontally over a very short distance and then turn downward. The crack growth path obtained by Bhattacharjee and Léger [49] corresponds to an overflow of about 10.0 m. It initially runs nearly horizontally before turning downward when the crack length reaches half the width of the dam neck. The damage zone is in good agreement with the results of Gioia et al. [48].Figure 12
Damage distribution (overflow = 10.2 m).Figure 13
Comparison of crack trajectories obtained by different methods.(
2) Different Initial Notch Lengths. To examine the effect of the initial notch length a on the carrying capacity and crack path, three initial lengths are, respectively, used for analysis, namely, a=0.1l, 0.2l, and 0.3l. The relationships between the overflow and the horizontal displacement of the dam crest for the three cases are shown in Figure 14. The corresponding damage distributions for an overflow of 10.2 m are shown in Figure 15.Figure 14
Overflow versus horizontal displacement of dam crest.Figure 15
Damage distributions for overflow of 10.2 m. (a)a= 0.1l, (b)a= 0.2l, and (c)a= 0.3l.
(a)
(b)
(c)It can be seen from Figure10 that (1) the critical overflow at which the dam is transformed from the linear state into the nonlinear state differs with the length of the initial notch, being 8.6, 7.4, and 7.0 m fora= 0.1l, 0.2l, and 0.3l, respectively, and (2) the ultimate overflows for the various values ofa are extremely close, being about 10.2 m. It can also be seen from Figure 15 that the cracks initially extend in the horizontal direction and then turn downward at nearly the same point and with almost the same deflection angle. The crack paths are also almost the same. The length of the initial notch thus only affects the critical overflow for transformation into the nonlinear state. The shorter the length, the higher the critical overflow. The notch length has no obvious effect on the ultimate carrying capacity of the dam and the crack path.


## 4.2.2. Crack Propagation considering Coupling Effect

To investigate the impact of the coupling effect during crack growth on the crack propagation path and carrying capacity of the dam, nonlinear stress-seepage-damage coupling analyses are conducted on the dam using the three initial notch lengths in Section4.2.1(2). The initial permeability coefficient k0 is set as 5 × 10−9 m/s and the coupling coefficienta0 is set as 0.01. It is really true that the dam which is 100 m high is constructed in reality with drains. Due to lack of information about the drains of Koyna dam, the drains are not taken into account in this study.In this section, in order to select an appropriate load step, five different load steps are chosen, 1.0 m, 0.5 m, 0.1 m, 0.05 m, and 0.01 m, respectively. Taking Koyna dam with an initial notcha = 0.1l as an example, the ultimate overflows under the five load steps are 8.0 m, 7.5 m, 7.2 m, 7.0 m, and 6.98 m, respectively. Assuming the result of the load step 0.01 m as the precise solution, the errors of the load steps 1.0 m, 0.5 m, 0.1 m, and 0.05 m are 14.61%, 7.45%, 3.15%, and 0.29%, respectively. It can be seen that when the load step is taken as 0.05 m, the value of the load step has little influence on the solution. Thus, the result of the load step 0.05 m is convergent and the load step 0.05 m can be chosen in the calculation of this section.(
1) Initial Notch Length
a
=
0.1
l. The relationships between the overflow and the horizontal displacement for an initial notch lengthaof 0.1l with and without consideration of the coupling effect are shown in Figure 16. The damage distribution with consideration of the coupling effect is shown in Figure 17. It can be seen from Figure 16 that, before the occurrence of obvious nonlinearity, the results for the two cases are in agreement. When the coupling effect is taken into consideration, the critical overflow that transforms the dam from the linear state into the nonlinear state is 4.1 m, and the ultimate overflow is 7.0 m. This represents significant decreases in both parameters compared to when not considering the effect. Further, from a comparison of Figures 15(a) and 17, it can be seen that consideration of the coupling effect decreases the horizontal crack propagation distance, while increasing the deflection angle of the downward turn and the total length of the path.Figure 16
Overflow versus horizontal displacement of dam crest.Figure 17
Damage distribution for overflow of 7.0 m.(
2) Different Initial Notch Lengths. Figure 18 shows the damage distributions for the three initial notch lengths corresponding to an overflow of 4.2 m with consideration of the coupling effect. The relationships between the overflow and the horizontal displacement of the dam crest for the three notch lengths are shown in Figure 19. The critical and ultimate overflows with and without consideration of the coupling effect are given in Table 3.Table 3
Comparison of dam carrying capacities with and without consideration of coupling effect.
Initial crack length
Without consideration of coupling effect
With consideration of coupling effect
Percentage change
Critical load (m)
Carrying capacity (m)
Critical load (m)
Carrying capacity (m)
Critical load
Carrying capacity
0.1l
6.0
10.2
4.1
7.0
−31.7%
−31.4%
02l
5.0
10.2
2.0
5.3
−60.0%
−48.0%
0.3l
3.0
10.2
1.0
4.2
−66.7%
−58.8%Figure 18
Damage distributions for overflow of 4.2 m. (a)a= 0.1l, (b)a= 0.2l, and (c)a= 0.3l.
(a)
(b)
(c)Figure 19
Overflow versus horizontal displacement of dam crest.It can be seen from Figure18 that, for a given overflow, the crack only propagates over a short distance in the horizontal direction whena= 0.1l, while it initially propagates in the horizontal direction and then turns downward whena= 0.2l. Whena= 0.3l, the crack turns downward almost directly and the propagation distance is nearly twice that fora= 0.2l. With consideration of the coupling effect, the crack propagation distance increases with increasing initial notch length. It can also be seen from Figure 19 that, with consideration of the coupling effect, the critical and ultimate overflows decrease with increasing initial notch length. Further, Table 3 indicates that, with consideration of the coupling effect, the rates of reduction of the critical and ultimate overflows increase with increasing initial notch length.


## 5. Conclusions

In this study, a stress-seepage-damage coupling model based on FEM is developed and first applied in hydraulic fracture of concrete gravity dams. The coupling model has the following characteristics: (1) The constitutive law of this coupling model considers the strain softening characteristic of fracture process zone, damage-dependent pore-pressure-influence coefficient, stress-dependent permeability for the prepeak stage, and deformation-dependent permeability for the postpeak stage; (2) the mesh-dependent hardening technique is adopted so that the fracture energy dissipation is not affected by the finite element mesh size; (3) four coupling processes during hydraulic fracture are considered: (a) the deformation of the surrounding medium induced by the water pressure on the fracture surface; (b) fluid flow within the fracture; (c) propagation of the fracture; (d) the leak-off of the fracturing fluid from the fracture into the surrounding medium.The damage model is first used to simulate crack propagation of a 1 : 40 scaled model dam. The crack trajectory and structural response agreed well with the experimental results. Then, the damage model is used to simulate crack propagation of Koyna dam with an initial notch on the upstream face under overflow. The numerically determined crack trajectory and structural response of Koyna dam agree well with the previously reported numerical results. These results validate the present damage model for predicting the carrying capacity and crack trajectory in concrete gravity dams.Subsequently, the crack propagation processes of Koyna dam using three notches of different initial lengths are simulated by the damage model and the coupling model. By comparing the results obtained by the two models, it is found that hydraulic fracture produces the following influences: (1) decrease in the horizontal propagation distance of the crack path, increases in the deflection angle when the crack turns downward, and the total propagation distance; (2) significant decreases in the critical load and ultimate carrying capacity of the dam; (3) increase in the crack propagation distance at the same load with increasing initial notch length; (4) increases in the rates of decrease of the critical and ultimate loads with increasing initial notch length. This demonstrates that hydraulic fracture significantly affects the crack growth path, propagation distance, and ultimate carrying capacity of the dam.Hydraulic fracture is a complex problem and the effects of damage on pore-pressure-influence coefficients and permeability coefficients need to be further studied.



---

*Source: 8523213-2017-04-23.xml*