# Numerical Simulation of Hydraulic Fracturing in Earth and Rockfill Dam Using Extended Finite Element Method.

1. Introduction

In recent years, with the needs of hydropower construction, hydropower energy development is increasing rapidly in China. Therefore, many high earth and rockfill dams have been designed or under construction, such as the Changhe (240 m in height) in Sichuan province, Shiziping (136 m in height) in Hebei province, and Nuozhadu (261.5 m in height) in Yunnan province. In addition, a number of super high (up to 300 m in height) earth and rockfill dam have been built, such as Lianghekou (295 m in height) in Shanxi province and Shuangjiangkou (322 m in height) in Guizhou province. With the rapid development of earth and rockfill dam, many accidents have occurred. Although no devastating disaster is caused by the accidents, safety operation of the dam is seriously affected. Consequently, many researchers have gradually paid attention to the safety design of high earth and rockfill dams.

Hydraulic fracturing is one of the most important factors affecting the safety of earth and rockfill dam. A series of studies have been carried out to simulate hydraulic fracturing. Ng and Small  presented a method of predicting hydraulic fracturing by using special joint elements that allow fluid flow and fracture to be modeled. Secchi and Schrefler  proposed a cohesive fracture model to simulate crack propagation of a concrete dam, which needs mesh updating of the crack tip. Barani and Khoei  presented a 3D cohesive crack propagation method in saturated porous media. By using double-nodded zero-thickness cohesive interface elements, the fracture behavior was simulated. Wang et al.  thought that the mechanism of hydraulic fracturing can be explained by fracture mechanics, and the crack propagation may follow any of mode I, mode II, and mixed mode I-II. Liu et al.  presented a stabilized extended finite element formulation to simulate the hydraulic fracturing process in an elastoplastic medium. The fracture propagation process was governed by a cohesive fracture model. Li et al.  used the smeared cracking theory to simulate the process of hydraulic fracturing in a certain earth and rockfill dams. Tang et al. [7-9] gave different mechanical parameters to different material units by using a certain statistical distribution. Thus, the crack propagation can be simulated. In general, the existing simulation methods require specific means, like grid refinement and estimate the position and direction of the fracture or set cohesive element, which result in the large calculating quantity and complex simulation process.

In 1999, Belytschko and Black  and Moes et al.  presented the extended finite element method (XFEM). The unique advantages of simulating crack propagation in homogenous material such as rock and concrete materials have been demonstrated . By using the XFEM, the crack is allowed to pass through the element; thus, cracks with complex shapes can be calculated in regular mesh without remeshing. Because of these advantages, the XFEM has been considered as a more accurate and efficient method in many areas, such as dynamic crack propagation, shear zone evolution, and multiphase flow. Yet, no researchers have used the XFEM to simulate crack problem in earth and rockfill dam engineering. In theory, it is possible to simulate hydraulic fracturing in earth and rockfill dam by using the XFEM.

In this paper, the XFEM is used to simulate the hydraulic fracturing in an actual high earth and rockfill dam. The possibility of hydraulic fracturing occurrence is analyzed, and the critical crack length is obtained when hydraulic fracturing occurs. Then, the crack propagation path and length is analyzed by inserting different length initial cracks at different elevation. Furthermore, based on the results, the possibility of crack penetration upstream and downstream of the core is analyzed. Related research achievements are of great significance for the prevention and control of hydraulic fracturing in actual earth and rockfill dam construction.

2. The Mechanism of Hydraulic Fracturing

A number of researchers did some meaningful works toward the mechanism of hydraulic fracturing in the core of earth and rockfill dam . The "water wedging" action has been approved by many engineers, which means that the core is not completely homogeneous and there exists different planes of weakness during the construction of the dam. Thus, the "permeable weak surfaces/bodies" occur at the upstream surface of the core as can be seen by ABC area in Figure 1. The permeability of these weak areas is higher than that of the soil around these weakens which allows water to enter into these weak areas rapidly as the reservoir water level rises. Consequently, there exists a directional water pressure on the AC surface and a downward water pressure on the BC surface as shown in Figure 1. If the water pressure is large enough and the core wall has poor anticrack ability, cracks may occur at the point C. Furthermore, the entry of high-pressure water will cause crack propagation, which will result in hydraulic fracturing.

3. Simulation Method of XFEM

3.1. Displacement Mode. Based on the partition of the unity method, the XFEM displacement formula can be derived by introducing an additional function which reflects the local characteristic of the crack on the basis of the conventional FEM displacement mode.

The displacement formula can be expressed as

[mathematical expression not reproducible], (1)

where [u.sub.I] is the nodal displacement vector, [[alpha].sub.I] is the nodal enriched degree of freedom vector for crack penetration elements, H(x) is the jump function which represents the gap between the crack surfaces, [b.sup.i.sub.I] is the nodal enriched degree of freedom vector for crack tip elements, and [[psi].sub.i] (x) is the crack tip functions which represent the crack tip singularity. The jump function H (x) for a general crack is defined as

[mathematical expression not reproducible]. (2)

The crack tip enrichment function for an isotropic elastic material is

[mathematical expression not reproducible]. (3)

3.2. The Application of Water Pressure. The hydrostatic pressure in the crack is applied on the crack surface as a facial pressure. As the crack expands and propagates continuously, the upstream pressure water penetrates into the crack tip slowly. Therefore, the hydrostatic pressure needs to be updated with the propagation of the crack in real time. In this paper, to simplify the problem, it is assumed that the water pressure is evenly distributed on the surface of the crack and the application of water pressure is carried out for each crack segment. According to the hydrostatic water pressure at different elevations, the water pressure is applied on the nodes with additional degrees of freedom (DOFs) in the crack segments.

3.3. Critical Water Pressure Inducing Hydraulic Fracturing Investigators carried out a great deal of works on hydraulic fracturing criterion. An elastic-plastic solution to fracture initiation pressure is derived based on the Mohr-Coulomb shear failure criterion and the theory of cavity expansion, which can be given by  the following equation:

[mathematical expression not reproducible], (4)

where [eta] = [([r.sub.a]/[r.sub.c]).sup.2] and M = 3[eta] + (2/[eta]) - 3[eta] sin [phi] + 7 sin [phi] - 1.

Actually, critical water pressure [p.sub.f] is the minimum value of [u.sub.a]. Based on the theory of expansion of a hollow cylinder in finite length, the radius of elastic-plastic zone boundary is obtained as follows :

[mathematical expression not reproducible], (5)

where Y = (2c cos [phi])/ (1 - sin [phi]), [alpha] = (1 + sin [phi])/(1 - sin [phi]), c is the cohesion, [phi] is the internal friction angle, [r.sub.a] and [r.sub.b] are internal and external radii of the hollow cylinder, respectively, [r.sub.c] is the radius of the plastic zone, [p.sub.0] is the radial stress of the external boundary, and p is the water pressure applying on the internal surface of cavity.

3.4. Damage Evolution Criterion. The damage evolution is determined by the definition of damage variable D. The expression of the normal and tangential stress components affected by the damage can be defined as

[mathematical expression not reproducible], (6)

where D is the damage variable, which represents the average damage value between the two cracks and varies between 0 and 1, [t.sub.n], [t.sub.s], and [t.sub.t] are the normal stress component and two tangential stress components, respectively, [T.sub.n] is the normal stress component of the unit, and [T.sub.s] and [T.sub.t] are the first and second tangential stress components of the unit, respectively.

4. Simulation of Hydraulic Fracturing in Earth and Rockfill Dam

4.1. Mesh, Constitutive Model, and Parameters. To obtain the length of crack propagation accurately, the local mesh refinement method is conducted at different elevation of the initial crack. The average length of the element is about 0.5 m after refined. Figure 2 shows the schematic diagram of the mesh with an initial crack at 1/5 height of the core wall.

The bottom of the model is constrained in three directions, while the sides of the model are constrained in the horizontal direction. The vertical direction is free to allow the settlement of the model. To simulate the construction of the dam, the load applied on the model is divided into ten stages.

Quadrilateral plane strain elements are used in the mesh, and a small number of triangular elements are used for transition. In addition, the reduced integration method is used in the XFEM enrichment area, and the convergence is controlled by the hourglass. Differently, the full integration method is used in other areas, which is similar to the conventional finite element method.

The Duncan-Chang E-v nonlinear constitutive model is used in this paper. The constitutive parameters are shown in Table 1. The fracture energy is 4.8 N/m and the tensile strength is 47 kPa of the core material obtained from uniaxial tensile test.

4.2. Research Framework and Scheme. Firstly, the possibility of hydraulic fracturing occurrence without initial crack in the core is analyzed by comparing the water pressure on the upstream with the calculated values [p.sub.f] using (4). The criteria can be implemented in the XFEM platform automatically though compiling subroutine UMAT. Then, apply "permeable weak surfaces" (assumed as initial cracks) of different lengths at different elevations in the core. Similarly, (4) is used to determine the occurrence of hydraulic fracturing. If the calculated critical water pressure is less than the water pressure at the crack tip, then the new crack is determined to appear. The propagation of the initial crack is determined by whether the accumulated energy of the crack propagation is more than the fracture energy of the core material. Thus, the critical crack length under the full reservoir water level can be obtained by establishing the relationship between the initial crack length and the reservoir water level when the hydraulic fracturing occurs. Finally, according to the initial cracks of different lengths at different elevation, the crack propagation path and length are studied by increasing the critical water level to the full reservoir water level gradually.

The calculation scheme is shown in Table 2, in which the reservoir water level is applied in two stages. If the water level is less than the maximum reservoir water level (697m), then the water pressure on the crack surface is equal to the pressure on the surface of the core at the same elevation. If the water level is greater than the reservoir water level, the water pressure on the surface of the core wall is the maximum reservoir water level while the water level on the crack surface is applied according to the actual water level. The main purpose of this method is to research the critical length of initial crack under full water level.

4.3. Analysis of Calculated Results

4.3.1. The Possibility of Hydraulic Fracturing Occurrence. Figure 3 shows the relationship between the initial fracturing pressure [p.sub.f] (calculated from formula (4)) and the hydrostatic pressure at the elevation where the initial crack locates. It can be seen that the initial fracturing pressure calculated at each elevation is greater than the hydrostatic pressure at its corresponding elevation. The result indicates that hydraulic fracturing will not occur in the earth and rockfill dam in this paper without permeable weak surfaces.

As the reservoir water level rises, the pressure water enters into the crack; thus, "water wedge effect" appears . The water pressure on the crack surface increases with the increasing water level. Consequently, hydraulic fracturing is likely to occur. The larger the length of the initial crack and the water pressure is, the stronger the water wedge effect is. In other words, there exists a critical length of the initial crack for an earth and rockfill dam under the normal water level.

To obtain the critical length of the initial crack, the relationship between the length of initial crack and the water level required for the propagation of the initial crack is established in Figure 4.

It can be seen that under the condition of full reservoir water level, the critical length of the initial crack required for hydraulic fracturing is 5.3 m, 5.9 m, and 7.9 m at the elevation of 507 m, 577 m, and 649 m, respectively. The results illustrate that the minimum length of the initial crack required for hydraulic fracturing is 5.3 m. Since the water pressure on the bottom is the largest, the possibility of hydraulic fracturing occurrence at the bottom of the core wall is greater than that of the upper part.

4.3.2. Analysis of the Crack Propagation. In the dam engineering, it is meaningful to analyze the path and length of the crack propagation when hydraulic fracturing happens under extreme conditions. Figures 5-7 show the nephogram of crack propagation with different lengths of initial crack at different elevation. XFEMSTATUS is a variable to describe crack propagation and changes between 0 and 1 and 1 represents that the crack is completely open. Considering the ratio of the crack to the whole dam is small, local areas of nephograms are given in order to show the propagation of the cracks clearly.

It is obvious that, under the constant driving of water pressure, the crack propagates to the inner part of core continuously. The larger the length of initial crack is, the longer the crack propagation will be. The direction of crack propagation changes from the horizontal direction to the upstream slowly. The reason is that the water pressure on the lower part of the core is larger than that of the upper part, which leads the core has a tendency to deflect.

With the increase of initial crack elevation, the length of crack propagation tends to decrease. For the reason that, the water pressure is decreasing as the elevation grows, when the initial crack is located at 4/5 height of the dam (elevation of 649 m), the propagation length is merely 5 m with a 12 m length initial crack.

The relation curves between the length of crack propagation and reservoir water level at 507 m elevation are shown in Figure 8. The reservoir water level is increased to the full reservoir level (697 m) gradually. Obviously, the crack propagation length and the reservoir water level have a linear relationship. The longer the length of initial crack, the greater the degree of crack propagation is. When the initial crack length is 10 m, the crack propagation reaches a maximum value, 14.7 m, which accounts for 15% of the width of the core wall. Therefore, the crack propagation is not enough to penetrate the core wall.

Figures 9 and 10 show the relation curves of the length of crack propagation versus the reservoir water level at the initial crack elevation of 577 m and 649 m. The reservoir water level required for the crack propagation is relatively close to the full reservoir water level, and the change of water pressure is relatively small. Therefore, the crack propagation length is longer than that in the lower part of the core wall. However, it is noteworthy that the width of the core wall decreases with the increasing elevation, which is merely 28.6 m at the elevation of 649 m. The critical length of crack is 7.9 m, which accounts for 27.6% of the core width.

As a consequent, although the occurrence probability of hydraulic fracturing at the bottom of the core wall is larger, the possibility of hydraulic fracturing to penetrate the upstream and downstream is extremely high when the upper part of the core wall reaches the critical crack length.

5. Conclusions

On the basis of the research on the possibility of hydraulic fracturing occurrence in the earth and rockfill dam, the hydraulic fracturing behavior is simulated by using the extended finite element method.

(1) For the calculation model in this paper, hydraulic fracturing will not occur without the permeable weak surface (initial crack), and the critical initial crack length required for hydraulic fracturing is 5.3 m.

(2) The propagation length of hydraulic fracturing decreases with the increase of elevation, and the average propagation length decreases from 9.4 m to 3.4 m.

(3) After the water level gradually increased to full reservoir level, there is no crack penetration through the core wall of all calculation schemes. The direction of crack propagation has a certain angle with the horizontal plane toward the downstream.

(4) Considering the up-narrow and down-wide types of the core wall, although the occurrence probability of hydraulic fracturing at the bottom of the core wall is larger, the possibility of hydraulic fracturing to penetrate the core is extremely high when the upper part of the core wall reaches the critical crack length.

https://doi.org/10.1155/2018/1782686

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors acknowledge the financial support from the National Key Research (2017YFC0404806) and Development Program of China (51779152).

References

 A. K. L. Ng and J. C. Small, "A case study of hydraulic fracturing using finite element methods," Canadian Geotechnical Journal, vol. 36, no. 5, pp. 861-875, 1999.

 S. Secchi and B. A. Schrefler, "A method for 3-D hydraulic fracturing simulation," International Journal of Fracture, vol. 178, no. 1-2, pp. 245-258, 2012.

 O. R. Barani and A. R. Khoei, "3D modeling of cohesive crack growth in partially saturated porous media: a parametric study," Engineering Fracture Mechanics, vol. 124-125, pp. 272-286, 2014.

 J. J. Wang, J. G. Zhu, H. Mroueh, and C. F. Chiu, "Hydraulic fracturing of rock-fill dam," International Journal of Multiphysics, vol. 1, no. 2, pp. 199-219, 2016.

 F. Liu, P. Gordon, H. Meier et al., "A stabilized extended finite element framework for hydraulic fracturing simulations," International Journal for Numerical and Analytical Methods in Geomechanics, vol. 41, no. 5, pp. 654-681, 2017.

 Q. M. Li, B. Y. Zhang, Y. Z. Yu et al., "Numerical simulation of the process of hydraulic fracturing in earth and rockfill dams," Chinese Journal of Geotechnical Engineering, vol. 29, no. 2, pp. 212-217, 2007.

 C. A. Tang, Z. Z. Liang, and T. Xu, "Numerical modelling of the influence of heterogeneity on 3D polygonal fracture in layered materials," Failure Mechanics Letters, vol. 2, pp. 1-2, 2005.

 C. A. Tang, Y. B. Zhang, Z. Z. Liang et al., "Fracture spacing in layered materials and pattern transition from parallel to polygonal fractures," Physical Review E, vol. 73, no. 5, p. 056120, 2006.

 C. A. Tang, Z. Z. Liang, Y. B. Zhang et al., "Fracture spacing in layered materials: a new explanation based on two-dimensional failure process modeling," American Journal of Science, vol. 308, no. 1, pp. 49-72, 2008.

 T. Belytschko and T. Black, "Elastic crack growth in finite elements with minimal remeshing," International Journal for Numerical Methods in Engineering, vol. 45, no. 5, pp. 601-620, 1999.

 N. Moes, J. Dolbow, and T. Belytschko, "A finite element method for crack growth without remeshing," International Journal for Numerical Methods in Engineering, vol. 46, no. 1, pp. 131-150, 1999.

 E. Gordeliy and A. Peirce, "Enrichment strategies and convergence properties of the XFEM for hydraulic fracture problems," Computer Methods in Applied Mechanics and Engineering, vol. 283, pp. 474-502, 2015.

 J. J. Wang, J. G. Zhu, C. F. Chiu et al., "Experimental study on fracture toughness and tensile strength of a clay," Engineering Geology, vol. 94, no. 1, pp. 65-75, 2007.

 J. G. Zhu, E. Y. Ji, Y. F. Wen et al., "Elastic-plastic solution and experimental study on critical water pressure inducing hydraulic fracturing in soil," Journal of Central South University, vol. 22, no. 11, pp. 4347-4354, 2015.

 H. S. Yu, Cavity Expansion Methods in Geomechanics, Kluwer Academic Publishers, AH Dordrecht, Netherlands, 2000.

 J. J. Wang and Y. X. Liu, "Hydraulic fracturing in a cubic soil specimen," Soil Mechanics and Foundation Engineering, vol. 47, no. 4, pp. 136-142, 2010.

Enyue Ji, (1,2) Zhongzhi Fu [ID], (1,2) Shengshui Chen, (1,2) Jungao Zhu, (3) and Zhizhou Geng (1)

(1) Geotechnical Engineering Department, Nanjing Hydraulic Research Institute, Nanjing 210024, China

(2) Key Laboratory of Failure Mechanism and Safety Control Techniques of Earth-Rock Dam of the Ministry of Water Resources, Nanjing Hydraulic Research Institute, Nanjing 210029, China

(3) Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering, Hohai University, Nanjing 210098, China

Correspondence should be addressed to Zhongzhi Fu; fu_zhongzhi@yahoo.com

Received 29 December 2017; Accepted 1 March 2018; Published 1 April 2018

Caption: FIGURE 1: Sketch diagram of hydraulic fracturing in the earth and rockfill dam.

Caption: FIGURE 2: Schematic diagram of mesh generation (1/5 dam height).

Caption: FIGURE 3: Relation curves of the initial fracturing pressure versus elevation.

Caption: FIGURE 4: Relation curves of the reservoir water level versus the length of initial crack.

Caption: Figure 5: Crack propagation nephogram at 507m elevation. (a) Initial crack with a length of 6 m. (b) Initial crack with a length of 8 m. (c) Initial crack with a length of 10 m.

Caption: FIGURE 6: Crack propagation nephogram at 577 m elevation. (a) Initial crack with a length of 6 m. (b) Initial crack with a length of 8 m. (c) Initial crack with a length of 10 m.

Caption: FIGURE 7: Crack propagation nephogram at 649 m elevation. (a) Initial crack with a length of 8 m. (b) Initial crack with a length of 10 m. (c) Initial crack with a length of 12 m.

Caption: FIGURE 8: Relation curves of reservoir water level versus length of crack propagation at an elevation of 507 m.

Caption: FIGURE 9: Relation curves of reservoir water level versus length of crack propagation at elevation of 577 m.

Caption: FIGURE 10: Relation curves of reservoir water level versus length of crack propagation at elevation of 649 m.
```TABLE 1: Constitutive parameters of the materials in the earth and
rockfill dam.

Number      Material type          [rho]        [[phi].sub.0]
(g/[cm.sup.3])    ([degrees])

1               Core                2.07            29.2
2             Rockfill              2.37            52.9
3            Overburden             1.45             48
4         Grout compartment         1.4              38

Number   [DELTA][phi]    c (kPa)       Duncan-Chang E-v model
([degrees])                [R.sub.f]     K       n

1             --           96         0.81       420     0.41
2            8.6            0         0.79       1352    0.38
3             7            --         0.80       1100    0.34
4             --            5         0.78       1500    0.23

Number   Duncan-Chang E-v model
G        F      D

1         0.43    0.088   1.5
2         0.30    0.13    8.6
3         0.33    0.010    4
4         0.34    0.010    4

TABLE 2: Calculation schemes of hydraulic fracturing.

Number   Elevation      Width of the     Length of    Reservoir water
of initial     core wall at      initial        level (m)
crack (m)        certain        crack (m)
elevation (m)

1           507             97.4           2-12      575-697, 697-1000
2           577             62.3           2-12      630-697, 697-930
3           649             28.6           2-12      685-697, 697-910
```