Extended Finite Element Method for Analyzing Hydraulic Fracturing of Rock Cracks under Compression

1. Introduction

Hydraulic fracturing refers to the expansion and extension of fractures existing within rock masses under the action of high-pressure water flow or other liquids [1]. Hydraulic fracturing was first applied in the petroleum industry to increase the output of oil wells, and was later used in fields such as in-situ stress measurement, nuclear waste storage and geothermal development [2]. Meanwhile, hydraulic fracturing can also have serious consequences for engineering projects. For example, rock slopes may slide under the action of groundwater, high-pressure water conveyance structures and reservoir dams may crack, and water gushing may occur during tunnel construction. Therefore, the problem of hydraulic fracturing in rock masses has become an issue that urgently needs to be solved at present.

Numerical simulation methods provide important tools for studying the mechanism of hydraulic fracturing in rocks. Among them, the Extended Finite Element Method (XFEM) is a relatively effective numerical simulation method for analyzing discontinuous problems [3]. When this method is used to analyze fracture problems, its computational grid is independent of the physical boundaries or the internal geometry of the structure. It overcomes the need for high-density meshing in the crack tip region, and there is no need to re-mesh when simulating crack propagation. Therefore, it can conveniently analyze the problem of hydraulic fracturing in rocks. Sheng Mao et al. [4] used the Extended Finite Element Method to simulate the propagation of a single hydraulic crack under the action of a constant water pressure. Shi Luyang et al. [5] introduced the uncoupled model of hydraulic fracturing into the Extended Finite Element Method to simulate the propagation of hydraulic fractures and natural fractures. Zhang et al. [6] combined the phase field method with the Extended Finite Element Method to simulate the hydraulic fracturing process of natural fractures in shale formations. Lecampion et al. [7] used the Extended Finite Element Method based on the uncoupled model to simulate the crack width and pressure distribution in hydraulic cracks. At present, there are relatively few studies on the mechanism of hydraulic fracturing of rock compression cracks using the Extended Finite Element Method.

In this paper, a numerical model for solving the problems of hydraulic fracturing and frictional contact of rock cracks using the extended finite element method is established. The computational model is applied to the analysis of hydraulic fracturing of specimens with central cracks and gravity dams with initial cracks under the action of axial pressure, so as to study the influence of axial compressive stress and initial crack length on the critical water pressure of hydraulic fracturing of specimens as well as the influence of hydraulic fracturing on the crack stability of gravity dams.

2. The Hydraulic Fracturing Model and Contact Model of the Extended Finite Element Method

2.1. The Approximate Displacement Function of the Extended Finite Element Method

The extended finite element method is based on the partition of unity concept. Additional functions that reflect the singularity at the crack tip and the discontinuity of the crack surface are added to the traditional finite element approximate displacement function formula. Its approximate displacement function is as follows [8]:

$$u={\displaystyle \sum _{i\in \Omega }{N}_i}\left(x\right)\left[u_i+\underset{i\in {\Omega }_{\Gamma }}{\underbrace{H\left(x\right)a_i}}+\underset{i\in {\Omega }_{\Lambda }}{\underbrace{{\displaystyle \sum _{l=1}^4{F}_l\left(x\right)b_i^{\left(l\right)}}}}\right]$$

(1)

where $u_i$is the nodal degrees of freedom of a conventional element; $a_i$is the additional degrees of freedom of a crack-penetrating element; $b_i^{\left(l\right)}$is the additional degrees of freedom of a crack-tip element; $\Omega$is the set of all nodes within the computational domain; ${\Omega }_{\Gamma }$ is the set of nodes of the crack-penetrating element; ${\Omega }_{\Lambda }$ is the set of nodes of the crack-tip element; $N_i$is the shape function of the node.

$H(x)$is a jump function that reflects the displacement discontinuity of the crack-penetrating element:

$$H(x)=\left\{\begin{array}{l}+1,\left(x-x^{*}\right)\cdot n>0\\ -1,\mathrm{otherwise}\end{array}\right.$$

(2)

$F_l\left(x\right)$is an additional function of the crack-tip element that reflects the singularity at the crack tip [9]:

$$F_l\left(x\right)=\sqrt{r}\left[\sin \left({\displaystyle \frac{\theta }{2}}\right)cos\left({\displaystyle \frac{\theta }{2}}\right)\right.\left.\sin \left(\theta \right)cos\left({\displaystyle \frac{\theta }{2}}\right)sin\left(\theta \right)\sin \left({\displaystyle \frac{\theta }{2}}\right)\right]$$

(3)

where$\left(r,\theta \right)$are the local polar coordinates at the crack tip (see Figure 1).

For the elements containing cracks, the relative displacement of the same point on the crack surfaces can be obtained from Equation (1):

$$w=u^{+}-u^{-}=2{\displaystyle \sum _{i\in {\Omega }_{\Gamma }}{N}_i{a}_i}+2\sqrt{r}{\displaystyle \sum _{i\in {\Omega }_{\Lambda }}{N}_i{b}_i^{\left(1\right)}}$$

(4)

where $w$represents the relative displacement of the same point on the crack surfaces.

2.2. Discrete Equations for Hydraulic Fracturing

After the displacement mode is constructed, the extended finite element governing equations for the hydraulic fracturing problem are derived through the principle of virtual work:

$${\displaystyle {\int }_{\Omega }\sigma :\delta \epsilon \mathrm{d}\Omega }={\displaystyle {\int }_{\Omega }b}\cdot \delta u\mathrm{d}\Omega +{\displaystyle {\int }_{{\Gamma }_t}t}\cdot \delta u\mathrm{d}S+{\displaystyle {\int }_{{\Gamma }_c}p}\cdot \delta w\mathrm{d}S$$

(5)

where $p$ is the water pressure on the crack surface ${\Gamma }_c$; $t$ is the external force on the boundary${\Gamma }_t$; $b$is the body force on the computational domain $\Omega$;$\sigma$, $\delta u$ and$\delta \epsilon$are the cauchy stress tensor, virtual displacement and virtual strain respectively.

Substituting the extended finite element approximate displacement expressions (1) and (4) into the extended finite element governing equation (5), the discrete equations for hydraulic fracturing can be obtained:

$$KU=f$$

(6)

The global stiffness matrix$K$is assembled from the element stiffness matrices$K_{ij}^{rs}$

$$K_{ij}^{rs}={{\displaystyle {\int }_{{\Omega }^e}\left(B_i^r\right)}}^{T}D{B}_j^{s}d\Omega \left(r,s=u,a,b\right)$$

(7)

the functional expressions of the strain matrices$B_i^a$, $B_i^b$, $B_i^u$ can be found in Reference [10].

$U$is the node displacement vector:

$$U={\left\{u_i{a}_i{b}_i^{(1)}{b}_i^{(2)}{b}_i^{(3)}{b}_i^{(4)}\right\}}^T$$

(8)

$f_i$is the equivalent nodal load vector of the body force$b$, the external force$t$and the water pressure$p$, and it can be expressed as:

$$f_i={\left\{f_i^u{f}_i^a{f}_i^{b1}{f}_i^{b2}{f}_i^{b3}{f}_i^{b4}\right\}}^T$$

(9)

$$\left.\begin{array}{l}{f}_i^u={\displaystyle {\int }_{{\Gamma }_t}{N}_i}td{\Gamma }_t+{\displaystyle {\int }_{{\Omega }^e}{N}_i}bd\Omega \\ f_i^a={\displaystyle {\int }_{{\Gamma }_t}{N}_i}Htd{\Gamma }_t+{\displaystyle {\int }_{{\Omega }^e}{N}_i}Hbd\Omega +2{\displaystyle {\int }_{{\Gamma }_c}n\cdot N_i}pd\Gamma \\ f_i^{bl}={\displaystyle {\int }_{{\Gamma }_t}{N}_i}{F}_{l}td\Gamma +{\displaystyle {\int }_{{\Omega }^e}{N}_i}{F}_{l}bd\Omega +2{\displaystyle {\int }_{{\Gamma }_c}n\sqrt{r}\cdot N_i}pd\Gamma \end{array}\right\}$$

(10)

where $l=1,2,3,4$.

2.3. Contact Model

The extended finite element governing equations for the frictional contact problem can be obtained from the principle of virtual displacement as follows:

$${\displaystyle {\int }_{\Omega }\sigma :\delta \epsilon d\Omega }={\displaystyle {\int }_{\Omega }b}\cdot \delta ud\Omega +{\displaystyle {\int }_{s_{\sigma }}t}\cdot \delta udS+{\displaystyle {\int }_{{\Gamma }_c}p}\cdot \delta w\mathrm{d}S+{\displaystyle {\int }_{\Gamma }\overline{p}}\cdot \delta wd\Gamma$$

(11)

where $\overline{p}$ is the contact force vector on the crack surface${\Gamma }_c$(in the global coordinate system).

It can be obtained from the governing equation (11), formula (5) and the discrete equation (6):

$$KU=f+\overline{P}$$

(12)

where $\overline{P}$ is the equivalent nodal load vector of the contact force vector$\overline{p}$on the crack surface.

The above formula can be written in the following form:

$$\overline{d}=\overline{R}\overline{p}+\overline{q}$$

(13)

where $\overline{R}=D{K}^{-1}C$, $\overline{q}=D{K}^{-1}f$, $\overline{q}$ are the relative displacements generated by the external load$f$on the crack surface.

with

$$\left.\begin{array}{l}D=\left[02N_j\left(x\right)2\sqrt{r}{N}_k\left(x\right)000\right]\\ C={\left[0N_i\left(x\right)\sqrt{r}{N}_i\left(x\right)000\right]}^T\end{array}\right\}$$

(14)

Equation (13) can be expressed in the following form:

$$d=Rp+q$$

(15)

where $R=T^T\overline{R}T$, $q=T^T\overline{q}$, $T$ are the global transformation matrices.

The above formula is the relational expression between the relative displacements and forces on the contact surface.

When dealing with contact problems using the conventional finite element method, it is required that the crack surface be consistent with the element boundary and that there are nodes on the crack surface. However, when handling contact problems with the extended finite element method, the crack surface can penetrate through the interior of the element and there are no nodes on the crack surface. In this paper, the contact relationship on the crack surface is reflected by using the relationship between the displacements and forces at the Gaussian points on the crack surface (see Figure 2).

The Gaussian points are contact point pairs. Considering the -th contact point pair, it should satisfy the geometric compatibility condition and the friction condition (see Reference [11]). These contact conditions can be written in the form of the following non-smooth equations of the nonlinear complementary type:

$$h_1^i=\min \left\{d_n^i,p_n^i\right\}=0$$

(16)

$$\begin{array}{l}{h}_2^i=d_{\tau }^i+\min \left\{0,\mu max\left\{0,p_n^i-d_n^i\right\}+p_{\tau }^i-d_{\tau }^i\right\}\\ +\max \left\{0,-\mu max\left\{0,p_n^i-d_n^i\right\}+p_{\tau }^i-d_{\tau }^i\right\}=0\end{array}$$

(17)

where $\mu$ is the Coulomb friction coefficient.

It can be seen from Equation (15) that the vector$d$ is a function of the vector$p$. Therefore, Equations (16) and (17) can be written as non-smooth equations of the nonlinear complementary type with$p$as the variable:

$$\left.\begin{array}{l}{h}_1^i\left(p\right)=0\\ h_2^i\left(p\right)=0\end{array}\right\}$$

(18)

The above formula is a non-smooth equation system model of the nonlinear complementary type for the plane friction contact problem. In this paper, the non-smooth damped Newton method [11] is adopted to directly solve the above non-smooth equation system.

3. Stress Intensity Factor Calculation and Crack Propagation Criterion

3.1. Calculation of Stress Intensity Factor

The stress intensity factor at the crack tip is calculated by the interaction integral method. The expression of the interaction integral $I^{\left(1,2\right)}$is as follows [12]:

$$I^{\left(1,2\right)}={\displaystyle {\int }_A\left({\sigma }_{ij}^{\left(1\right)}\frac{\partial u_i^{\left(2\right)}}{\partial x_1}+{\sigma }_{ij}^{\left(2\right)}\frac{\partial u_i^{\left(1\right)}}{\partial x_1}-W^{\left(1,2\right)}{\delta }_{1j}\right)}{\displaystyle \frac{\partial Q}{\partial x_j}}dA-{\displaystyle {\int }_{{\Gamma }_c}\left(p_j^{\left(1\right)}\frac{\partial u_i^{\left(2\right)}}{\partial x_1}+p_j^{\left(2\right)}\frac{\partial u_i^{\left(1\right)}}{\partial x_1}\right)}Qd\Gamma$$

(19)

where $u_i$, ${\sigma }_{ij}$and${\epsilon }_{ij}$ are the displacement vector, stress tensor and strain tensor respectively;$A$ is the integral region around the crack tip;${\delta }_{1j}$ is the Kronecker delta;$Q$ is the weight function;$p_j$ is the water pressure on the crack surface; the superscript 1 represents the real state and takes the numerical solution of the extended finite element method; the superscript 2 represents the auxiliary state and takes the crack tip asymptotic field.

$W^{\left(1,2\right)}={\displaystyle \frac{1}{2}}\left({\sigma }_{ij}^{\left(1\right)}{\epsilon }_{ij}^{\left(2\right)}+{\sigma }_{ij}^{\left(2\right)}{\epsilon }_{ij}^{\left(1\right)}\right)$is the interaction strain energy density of states 1 and 2, and its expression is as follows:

$$W^{\left(1,2\right)}={\displaystyle \frac{1}{2}}\left({\sigma }_{ij}^{\left(1\right)}{\epsilon }_{ij}^{\left(2\right)}+{\sigma }_{ij}^{\left(2\right)}{\epsilon }_{ij}^{\left(1\right)}\right)$$

(20)

There is the following relationship between the stress intensity factor and the interaction integral [12]:

$$I^{\left(1,2\right)}={\displaystyle \frac{2}{E^{*}}}\left(K_{{\rm I}{\rm I}}^{\left(1\right)}{K}_{{\rm I}{\rm I}}^{\left(2\right)}+K_{{\rm I}}^{\left(1\right)}{K}_{{\rm I}}^{\left(2\right)}\right)$$

(21)

for plane strain,$E^{*}=E/1-{\nu }^2$; for plane stress,$E^{*}=E$.

Selecting state 2 as the asymptotic fields of mode I and mode II, the stress intensity factors of mode I and mode II for state 1 can be obtained: $K_{{\rm I}}={\displaystyle \raisebox{1ex}{$E^{\ast }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{I}^{\left(1,mode{\rm I}\right)}$, $K_{{\rm I}{\rm I}}={\displaystyle \raisebox{1ex}{$E^{\ast }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{I}^{\left(1,mode{\rm I}{\rm I}\right)}$.

3.2. Crack Propagation Criterion

The crack propagation criterion adopts the maximum circumferential stress criterion to obtain the crack propagation angle${\theta }_c$:

$${\theta }_c=2\arctan {\displaystyle \frac{1}{4}}\left(K_{{\rm I}}/K_{{\rm I}{\rm I}}\pm \sqrt{{\left(K_{{\rm I}}/K_{{\rm I}{\rm I}}\right)}^2+8}\right)$$

(22)

The equivalent stress intensity factor$K_{{\rm I}}^{\mathrm{eq}}$ is determined by the following formula:

$$K_{{\rm I}}^{\mathrm{eq}}={\displaystyle \frac{1}{2}}\cos \left({\displaystyle \frac{{\theta }_c}{2}}\right)\left[K_{{\rm I}}\left(1+\cos {\theta }_c\right)-3K_{{\rm I}{\rm I}}\sin {\theta }_c\right]$$

(23)

4. Numerical Examples

4.1. A Plate with a Single-Edge Crack

The rectangular plate with the length of a single-edge crack being$c=2m$ is shown in Figure 3. The width of the flat plate is$w=5m$, the height is$h=5m$ , and the uniformly distributed water pressure$p=1MPa$ acts on the crack surface. The elastic modulus of the material of the rectangular plate is$E=12\mathrm{GPa}$, and the Poisson's ratio is$\nu =0.2$.

Figure 3. A rectangular plate with edge crack.

Figure 3. A rectangular plate with edge crack.

Figure 4. Mesh generation.

Figure 4. Mesh generation.

According to the superposition principle [13], the exact solution of the stress intensity factor when the crack surface is subjected to the action of uniformly distributed water pressure is

$$K_{\mathsf{{\rm I}}}^{\mathrm{exac}}=\left[1.12-0.231\left({\displaystyle \raisebox{1ex}{$c$}\!\left/ \!\raisebox{-1ex}{$w$}\right.}\right)+10.55{\left({\displaystyle \raisebox{1ex}{$c$}\!\left/ \!\raisebox{-1ex}{$w$}\right.}\right)}^2\right.\left.-21.72{\left({\displaystyle \raisebox{1ex}{$c$}\!\left/ \!\raisebox{-1ex}{$w$}\right.}\right)}^3+30.39{\left({\displaystyle \raisebox{1ex}{$c$}\!\left/ \!\raisebox{-1ex}{$w$}\right.}\right)}^4\right]\sigma \sqrt{\mathsf{\pi }c}$$

(24)

Let the normalized stress intensity factor

$$K_{{\rm I}}^n={\displaystyle \frac{K_{{\rm I}}^{XFEM}}{K_{{\rm I}}^{exac}}}$$

(25)

where $K_{{\rm I}}^{XFEM}$is the stress intensity factor calculated by the method proposed in this paper.

Table 1. Normalized SIF values for different elements.

	the number of mesh elements
	1225	2025	4225	5625	9025	13225
$K_{{\rm I}}^n$	0.9011	0.9292	0.9627	0.9734	0.9894	1.0008

Table 1. Normalized SIF values for different elements.

	the number of mesh elements
	1225	2025	4225	5625	9025	13225
$K_{{\rm I}}^n$	0.9011	0.9292	0.9627	0.9734	0.9894	1.0008

The influence of the number of mesh elements on the value of the normalized stress intensity factor is shown in Figure 5. It can be seen from the figure that high computational accuracy can be achieved even when the number of elements is relatively small. As the number of elements increases, the computational error gradually decreases. When the number of mesh elements is around 10,000, the influence of the number of mesh elements on computational accuracy can be neglected.

Figure 6 shows the stress intensity factors under the action of different water pressures on the crack surface when the number of mesh elements is 13,225 (the mesh generation is shown in Figure 4). It can be seen from the figure that the values of the stress intensity factors calculated by the method in this paper are consistent with the exact solutions, indicating that the method in this paper is accurate in calculating the stress intensity factors when the crack surface is subjected to the action of uniformly distributed water pressures.

4.2. A Plate with a Through Crack

Figure 7 shows a rectangular plate with a through crack. The friction coefficient of the crack surface is$\mu =0.3$, the bottom is under fixed constraints. Both the length and width of the flat plate are 10 meters. Its top is subjected to the tangential uniformly distributed shear force$\tau =25KPa$ and the normal uniformly distributed pressure$\sigma =100KPa$, The elastic modulus is$E=200GPa$, and the Poisson's ratio is$\nu =0.3$.

Figure 8a shows the extended finite element mesh, which adopts four-node isoparametric elements. The computational model is divided into 100 elements in total. The numerical example has undergone 5 iterations of calculation and met the convergence requirements, with a tolerance of 9.1104e^-9. Figure 9 shows the displacement contour plot. It can be seen from the figure that under the action of compressive-shear stress, the horizontal and vertical displacements of the upper and lower parts of the crack surface are continuous, indicating that the contact algorithm in this paper is effective. The stress contour plot is shown in Figure 10. As can be seen from the figure, the normal normal stress on the crack surface is continuous, while the tangential normal stress is discontinuous, which conforms to the conditions of the compressive stress contact surface.

The contact algorithm in this paper is compared with the finite element penalty function method to verify the accuracy of the contact algorithm in this paper. The finite element mesh generation is shown in Figure 8b, and the computational model is divided into 400 elements in total. The contact stresses on the crack surface calculated by the finite element penalty function method and the contact algorithm in this paper are shown in Figure 11. It can be seen from the figure that the calculation results of the contact algorithm in this paper and the finite element penalty function method are quite consistent, indicating that the contact algorithm in this paper is accurate. Table 2 presents the normal contact forces, normal relative displacements as well as tangential contact forces and tangential relative displacements of Gaussian points 1 to 20 on the crack surface.

4.3. A Plate with a Central Crack

Figure 12 shows a rock specimen with a central crack. The length of the central crack is 2 meters, the height of the specimen is 15 meters, and the width is 10 meters. The top of the specimen is subjected to an external axial pressure$\sigma =1MPa$, and the crack surface is acted upon by a uniformly distributed water pressure$p=3MPa$. The elastic modulus of the rock specimen is$E=10\mathrm{GPa}$, the Poisson's ratio is$\nu =0.3$, the fracture toughness is$K_{{\rm I}c}=0.6M\mathrm{N}·{\mathrm{m}}^{-3/2}$, and the friction coefficient of the crack surface is$\mu =0.3$.

Considering the structural symmetry of the specimen, the right half of the specimen is taken. The mesh generation is shown in Figure 13. The computational model is divided into 12,675 elements in total. The numerical example has undergone 4 iterations of calculation and met the convergence requirements, with a tolerance of 1.1586e^-10. Figure 14 shows the stress contour plot. It can be seen from the figure that under the combined action of the external axial pressure and the internal water pressure, an obvious stress concentration phenomenon appears at the crack tip of the rock specimen with a central crack. The stress intensity factors under different crack lengths are shown in Figure 15, indicating that under the combined action of the external axial pressure and the internal water pressure, the stress intensity factor at the crack tip of the rock specimen increases as the crack length increases.

The water pressure inside the crack starts from 0 MPa and gradually increases with an increment of 0.1 Mpa. When the stress intensity factor at the crack tip is greater than the fracture toughness, the water pressure inside the crack at this time is taken as the critical water pressure. Figure 16 shows the crack propagation path under the action of the critical water pressure of 3 Mpa. The number of crack propagation steps is 8, and the crack propagation step length is taken as 0.5 m. The crack propagation criterion is the maximum circumferential stress criterion. It can be seen from the figure that the crack starts to propagate from the tip and extends along the horizontal direction. The numerical calculation results are the same as the experimental results in Reference [14].

Figure 17 shows the critical water pressures under different crack lengths. It can be seen from the figure that the critical water pressure decreases as the crack length increases. The critical water pressures under different external axial pressures are shown in Figure 18. It can be seen from the figure that as the external axial pressure increases, the critical water pressure also increases. The calculation results are consistent with the experimental results in Reference [14]. The axial pressure has an inhibitory effect on crack initiation, and the critical water pressure will increase.

4.4. A Gravity Dam with an Initial Crack

There is an initial crack with a length of 10 meters at the foundation of a concrete gravity dam. The angle between this crack and the horizontal section of the dam heel is 30^°. The height of the dam is 192 meters, the width at the bottom of the dam is 159.8 meters, and the width at the top of the dam is 10 meters. The scope of the dam foundation extends 192 meters upstream, downstream and in the depth of the dam foundation respectively, as shown in Figure 19. The upstream and downstream boundaries of the dam foundation are under normal constraints, and the bottom of the dam foundation is under fixed constraints. The computational parameters related to the materials of the gravity dam are shown in Table 3. The quadrilateral isoparametric elements are adopted for numerical analysis. The mesh generation is shown in Figure 20, with the number of elements being 5920 and the number of nodes being 6109.

The main loads considered in the calculation include the self-weight of the dam, the upstream hydrostatic pressure and the water pressure inside the crack. For the hydrostatic pressure, the method of gravity overload is adopted. Only the water gravity is increased while the water level remains at 192 meters. The overload coefficient is the ratio before and after the increase of the water pressure. In this example, the overload coefficient is taken as 2.5.

The maximum circumferential stress criterion is adopted as the crack propagation criterion, and the crack propagation step length is taken as 3 meters. There are a total of 12 crack propagation steps. Figure 21 shows the stress contour plot of the gravity dam obtained by the extended finite element calculation. It can be seen from the figure that stress concentration phenomena occur in both the dam heel and the crack tip areas, and the stress value in the dam heel area is significantly lower than that in the crack tip area.

stress; (b) ${\sigma }_{yy}$stress; (c) ${\sigma }_{xy}$stress.

The crack propagation paths with and without considering hydraulic fracturing are shown in Figure 22. The crack propagation angle in the case of considering hydraulic fracturing is smaller than that in the case of not considering hydraulic fracturing.

The relationship curve between the stress intensity factor and the crack length is shown in Figure 23. The mode-I stress intensity factor in the case of considering hydraulic fracturing has fluctuations, while the mode-I stress intensity factor in the case of not considering crack hydraulic fracturing is relatively flat. As the crack continues to propagate, the mode-I stress intensity factor in the case of considering hydraulic fracturing is greater than that in the case of not considering hydraulic fracturing, indicating that hydraulic fracturing will increase the mode-I stress intensity factor at the crack tip, thereby reducing the stability of the crack at the foundation of the gravity dam.

5. Conclusions

This paper takes advantage of the extended finite element method in solving problems on moving discontinuous surfaces. By considering the virtual work principle of the water pressure distribution on the crack surface and the frictional contact of the crack surface, the governing equations for analyzing hydraulic fracturing and frictional contact problems using the extended finite element method are derived. A numerical model of the extended finite element method for solving the problems of hydraulic fracturing and frictional contact of rock cracks is established. The calculation model is applied to the hydraulic fracturing analysis of the specimens with central cracks and the gravity dams with initial cracks under the action of axial pressure. The following conclusions are obtained:

(1)

The values of the stress intensity factors calculated by the method in this paper are consistent with the exact solutions, indicating that the method in this paper is accurate in calculating the stress intensity factors when the crack surface is subjected to the action of uniformly distributed water pressures. Moreover, as the number of elements increases, the computational error gradually decreases.

(2)

The contact algorithm in this paper can effectively prevent the crack surfaces from interpenetrating. Its results are consistent with the calculation results of the finite element penalty function method, indicating that the contact algorithm in this paper is accurate.

(3)

Under the combined action of the external axial pressure and the internal water pressure, the stress intensity factor at the crack tip of the rock specimen increases as the crack length increases. The critical water pressure decreases as the crack length increases; the critical water pressure increases as the external axial pressure increases. The axial pressure has an inhibitory effect on crack initiation, and the critical water pressure will increase.

(4)

Hydraulic fracturing will increase the mode-I stress intensity factor at the crack tip and reduce the stability of the cracks at the foundation of the gravity dam.