A New Method and Application of Full 3D Numerical Simulation for Hydraulic Fracturing Horizontal Fracture

Abstract

The numerical simulation of hydraulic fracturing fracture propagation is the core content of hydraulic fracturing design and construction. Its simulation results directly affect the effect of fracturing, and can effectively guide the fracturing construction plan and reduce the construction risk. At present, two-dimensional or quasi-three-dimensional models are mainly used, but most of them are used to simulate the vertical fracture of hydraulic fracturing. There are errors in the application process. In this paper, a three-dimensional mathematical model, including an elastic rock mechanics equation and a material flow continuity equation, is established to simulate horizontal fracture propagation in shallow reservoirs. The emphasis of this paper is to propose a new method for solving equations. The basic idea of the iteration method has been proposed by previous scholars: Firstly, assuming that the initial pressure of each point in the fracture is uniform, the fracture height of each initial point can be obtained by using Equation (20). Using the initial height values, the pressure values at each point of continuous variation are calculated by Equation (16), and then the new fracture height values are obtained by Equation (20). Because of the equal initial pressure, this method leads to too many iterations in the later stages, which makes the calculation more complicated. In this paper, a new Picca iteration method is proposed. The iteration parameters are changed sequentially. Firstly, the distribution value of fracture height is assumed. Then, the pressure distribution value is calculated according to Equation (16). Then, the new distribution value of fracture height is obtained by bringing the obtained pressure distribution value into Equation (20). Then, the new distribution value of the fracture height is calculated according to Equation (16). The pressure distribution value completes an iteration process until the iteration satisfies the convergence condition. In addition, Sneddon’s model is introduced into the hypothesis of fracture height to obtain the maximum fracture height and assume that the initial fracture profile is a parabola. Finally, the proposed method can rapidly improve the convergence rate. Next, on the basis of investigating the solutions of previous equations, the Galerkin finite element method is used to solve the above equations. The new Picard iteration sequence method is applied to solve the height and pressure at different points in the fracture. By calculating the stress intensity factor, we can judge whether the fracture continues to extend or not, and then simulate the full three-dimensional horizontal fracture of the hydraulic fracturing expansion process. The infiltration process of three types of oil reservoirs in Daqing Changyuan oilfield is simulated. The results show that during the initial fracture stage, the radius and height of fractures increase rapidly, and the rate of increase slows down with the increase of construction time. The height and net pressure of each point in the fracture are unequal. The height and net pressure of the fracture in the wellbore reach the maximum, and gradually decrease to the front of the fracture. Compared with conventional fracturing, the fracturing-flooding percolation process has the characteristics of short fracture-making and large vertical percolation distance, which can greatly increase the swept volume of flooding fluid and thus enhance oil recovery. With the increase in the rock modulus of elasticity, the radius of fractures decreases and the height of fractures increases. With the increase in construction displacement, the radius of fractures hardly changes, the height of fractures increases, and the vertical infiltration distance of the fractures increases. It is suggested that the construction displacement should be 4.0 m³/min. In the range of fracturing fluid viscosity in the studied block, with the change of fracturing fluid viscosity, the change of fracture radius and height is not obvious. In order to further increase sweep volume, the fracturing fluid viscosity should be further reduced.

1. Introduction

Hydraulic fracturing technology has been used in oil and gas field development engineering since the 1950s, and has been widely used, having broad prospects. But before the 1980s, two-dimensional models were used in hydraulic fracturing design and calculation at China and abroad. The representative models include the PKN (Perkin & Kern) model, KGD (J. Geertsma & De Klerk) model [1,2]. However, these models are used to simulate vertical fractures. A large number of experiments and field studies [3,4,5,6,7,8] show that the height of fractures is not a fixed value, which limits the use of the model. At the same time, there is no two-dimensional classical calculation model for horizontal fractures. For large-scale fracturing process, it is necessary to use a three-dimensional model to simulate accurately. The three-dimensional model includes a quasi-three-dimensional model and full-three-dimensional model. The quasi-three-dimensional model assumes three-dimensional fracture extension and one-dimensional fluid flow [9,10]. (Fracturing fluid flows in one dimension along the fracture length inside the fracture. Fluid does not flow in the direction of fracture height.) While the full three-dimensional model considers three-dimensional fracture propagation and two-dimensional fluid flow in the fracture [11] and the filtration factor. Taking it into account, the fracturing fluid seeps into the formation perpendicular to the fracture wall, which is suitable for various reservoir conditions and simulates the physical process of hydraulic fracturing more truly. In recent years, many scholars have studied [12,13,14,15,16,17,18] the mechanism of fracture initiation and propagation, the filtration of fracturing fluid, and the change of temperature field. Chen [12] applied the theory of rock mechanics, seepage mechanics and fracture mechanics, adopted the maximum tensile stress criterion, used cohesive element of pore pressure to simulate the fracture, established the calculation model of hydraulic fracturing fracture propagation, and analyzed the influencing factors of fracture propagation. However, the assumption of the model takes into account the two-dimensional flow of liquid in both horizontal and vertical directions, which is actually a quasi-three-dimensional model. Zhang et al. [13] summarized the development history of hydraulic fracturing simulation, compared and analyzed the two-dimensional model with the quasi-three-dimensional model and the full-three-dimensional model, and realized that the full-three-dimensional model could predict and optimize the field hydraulic fracturing operation more truly and accurately. Shen et al. [14] proposed a new solution to the quasi-three-dimensional problem, which reduced the complexity of the numerical solution of the quasi-three-dimensional problem. Although the computation efficiency is improved, the fluid in the fracture is still a one-dimensional flow, which is quite different from the actual one. Xue et al. [15] carried out a three-dimensional finite element numerical simulation study on the hydraulic fracturing process. Presupposed fractures were simulated by bonding elements, and fracture initiation and opening were characterized by damage factors of elements. In order to make the fluid pressure dynamically track and load on the fracture surface with the fracture propagation, the fluid pressure on the fracture surface was given. The model is passed internally and the user subroutine is written to implement it. Qiao et al. [16] deduced the two-dimensional temperature field equation of fracturing fluid in hydraulic fracturing. The heat conduction, convection, and dissipation of fracturing fluid and the heat exchange between fracturing fluid and rock were taken into account. The equation and the integral equation describing the relationship between fracture width and pressure on the fracture surface in the three-dimensional fracture model of hydraulic fracturing and the flow equation of fracturing fluid were also considered. The criterion for fracture extension is the stress intensity factor. It can be seen that when considering the temperature effect, the two-dimensional model is still used to effectively simulate the fracture growth. Wu et al. [17] assumes that the fracture is an ellipse, describes the width equation, pressure drop equation, extension criterion, and continuity equation of the fracture, and solves the model by Hamming method. However it assumes that the fracture is ellipse, and the equation representing the width and pressure of the fracture is basically the empirical equation of two-dimensional hydraulic fracturing. The simulation results are quite different from the three-dimensional ones. Zhang et al. [18] established a new three-dimensional hydraulic fracturing fracture extension model. The main development of the model is to take into account the effects of stress between production layer, cap rock, and bottom layer and the changes of rock mechanics parameters (modulus of elasticity, Poisson’s ratio, and fracture toughness). It can simulate various stress distribution patterns and the extension before and after fracture penetration. However, the method adopted is effective differential method, which causes large errors. Although many scholars have [12,13,14,15,16,17,18] studied the propagation of fractures differently, most adopted two-dimensional or quasi-three-dimensional models, which have too many assumptions and are not suitable for field engineering application. Therefore, it is necessary to calculate and analyze the full three-dimensional equations and put forward a new solution to realize the simulation of fracture propagation in full three-dimensional hydraulic fracturing horizontal fracture.

In this paper, the elastic rock mechanics equation and fluid flow continuity equation of full three-dimensional hydraulic fracturing horizontal fracture are given. The matrix discretization of the equation is carried out by using the finite element method, and the calculation methods of double integral and area integral are defined. The triangular element is used for mesh generation. A new iteration method and fracture expansion judgement basis are given. Through the analysis of three types of reservoir fracturing, horizontal fracture height, fracture pressure, fracture radius, and the difference between fracturing percolation flooding technology and traditional fracturing technology are clarified. The simulation results can effectively guide engineering fracturing construction design.

2. Mathematical Physics Equation

The pseudo-three-dimensional model [19] assumes that the fracturing fluid flows one-dimensional along the fracture length inside the fracture. Fluid does not flow in the direction of fracture height, so the calculation model can be simplified to realize a simple calculation, as shown in Figure 1.

Full three-dimensional hydraulic fracturing propagation model [19] assumes that fracturing fluid flows in two-dimensional vector direction parallel to the fracture wall, and takes the filtration behavior of fracturing fluid into account. It is different from the pseudo three-dimensional model. In the pseudo three-dimensional model, the fracture length (R) is much larger than the fracture height (h, vertical fracture), the ratio of fracture length to height is more than 3.5–5 times [9,10] and it applies the one-dimensional flow assumption of fluid in the fracture. So the full three-dimensional hydraulic fracturing propagation model can be applied in various lithological conditions. The basic equations of full three-dimensional hydraulic fracturing expansion model include the fluid pressure and elastic deformation equation of fracture height acting on fracture wall, the fluid flow equation of fracturing fluid flow and hydraulic gradient in the fracture, namely the elastic rock mechanics equation and material flow continuity equation. The full three-dimensional horizontal fracture plan is shown in Figure 2.

2.1. Elastic Rock Mechanics Equation

Assuming that the stratum is an isotropic linear elastic body and the horizontal fracture is perpendicular to the vertical stress, without considering the effect of the horizontal stress on the fracture morphology, the surface integral method can be used to simplify the full three-dimensional mechanical problem in an infinite medium into two-dimensional problems in a finite area, and the surface integral method can be used to characterize the fractures. The relationship between normal stress and fracture height of wall [20]:

$$\Delta p(x,y)\equiv p(x,y)-\sigma (x,y)=E{\displaystyle {\int }_A[\frac{\partial w}{\partial x^{\prime }}\frac{\partial }{\partial x}(\frac{1}{R})+\frac{\partial w}{\partial y^{\prime }}\frac{\partial }{\partial y}(\frac{1}{R})]dA}$$

(1)

In the Equation (1), $\Delta p(x,y)$ is the net fluid pressure acting on fracture wall, $\mathrm{Pa}$; $p(x,y)$ is the liquid pressure of fracturing fluid in fracture, $\mathrm{Pa}$; $\sigma (x,y)$ is the minimum principal stress perpendicular to fracture wall, $\mathrm{Pa}$; $E$ is the equivalent elastic modulus of rock, $\mathrm{Pa}$; $w$ is the horizontal fracture height, $\mathrm{m}$; $R$ is the distance between the integral points $(x^{\prime },y^{\prime })$ of the integrand function and pressure action point $(x,y)$, $\mathrm{m}$; and $A$ is the horizontal fracture extent in the integral plane area, ${\mathrm{m}}^2$.

$$R={[{(x^{\prime }-x)}^2+{(y^{\prime }-y)}^2]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$$

(2)

$$E=\frac{G}{4\pi (1-\nu )}$$

(3)

In the Equation (3), $G$ is the rock shear modulus, $\mathrm{Pa}$; and $\nu$ is the rock Poisson’s ratio, decimal.

Because the integral term in the equation has the same integral point and pressure action point, the integral function is infinite at this time. Hongren Gu and Cheng H. Yew [21] introduced the potential function and transformed the term $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$R$}\right.$ differential of the integral function into the differential of the potential function. The equation after singularity reduction is as follows [20,21]:

$${\displaystyle {\int }_A[p(x,y)-\sigma (x,y)]V(x,y)}=E{\displaystyle {\int }_A{\displaystyle {\int }_A\frac{1}{R}}[\frac{\partial w}{\partial x^{\prime }}\frac{\partial V}{\partial x}+\frac{\partial w}{\partial y^{\prime }}\frac{\partial V}{\partial y}]d{x}^{\prime }d{y}^{\prime }dxdy}$$

(4)

2.2. Material Flow Continuity Equation

The fracturing fluid is an incompressible Newtonian fluid and has a laminar flow on two basically parallel fracture surfaces. We ignore the inertia effect and the volume force of fracturing fluid in the direction of horizontal fracture height (vertical fracture must be considered). The filtration rate of fracturing fluid is determined according to the comprehensive filtration coefficient and filtration time. The Navier–Stoke equation of fluid flow is used to determine the filtration rate of the fracturing fluid equation [22]:

$$\frac{\partial p(x,y)}{\partial x}=\mu \frac{\partial }{\partial z}(\frac{\partial v_x}{\partial z})$$

(5)

$$\frac{\partial p(x,y)}{\partial y}=\mu \frac{\partial }{\partial z}(\frac{\partial v_y}{\partial z})$$

(6)

In the Equations (5) and (6), $\mu$ is the viscosity of fracturing fluid, $\mathrm{Pa}·\mathrm{s}$; z is the fracture height direction, $\mathrm{m}$; and $v_x$,$v_y$ is the flow velocity along X and Y directions, $\mathrm{m}/\mathrm{s}$.

Equation (5) is integrated twice and the boundary condition on the fracture wall is that the velocity of flow is equal to 0 [22]:

$$v_x=-\frac{1}{2\mu }[{(\frac{w}{2})}^2-z^2]\frac{\partial p(x,y)}{\partial x}$$

(7)

$$v_x=-\frac{1}{2\mu }[{(\frac{w}{2})}^2-z^2]\frac{\partial p(x,y)}{\partial x}$$

(8)

The volume flow rate per unit fracture length is:

$$q={\displaystyle \sum _{-w/2}^{w/2}vdz}$$

(9)

The equilibrium conditions for the flow equation are as follows:

$$-\frac{\partial q_x}{\partial x}-\frac{\partial q_y}{\partial y}-q_l-q_I=\frac{\partial w}{\partial t}$$

(10)

In the Equation (10), $q_x$ and $q_y$ are the volume flow rate along X direction (Y direction) per unit length in Y direction (X direction), ${\mathrm{m}}^2/\mathrm{s}$; $q_I$ is the injection rate of fracturing fluid on unit length of fracture cross section, m/s; and $q_l$ is the volume filtration rate on unit area of fracture wall, $\mathrm{m}/\mathrm{s}$, its expression is as follows [23]:

$$q_l=\frac{2c_l}{\sqrt{t-\tau (x,y)}}$$

(11)

In the Equation (11), $c_l$ is the comprehensive fracturing fluid filtration coefficient, $\mathrm{m}/\sqrt{\mathrm{s}}$ and $\tau (x,y)$ is the time to start filtration at a point $(x,y)$ on the fracture wall, $\mathrm{s}$. By introducing Equations (7) and (8) into (9) we establish Equation (10), and the governing equations of fluid flow in fractures can be obtained:

$$\frac{\partial }{\partial x}(\frac{w^3}{12\mu }\frac{\partial p(x,y)}{\partial x})+\frac{\partial }{\partial y}(\frac{w^3}{12\mu }\frac{\partial p(x,y)}{\partial y})=q_I-\frac{\partial w}{\partial t}-\frac{2c_l}{\sqrt{t-\tau (x,y)}}$$

(12)

The condition equation for this equation has a unique solution:

$${\displaystyle {\int }_{\partial A_p}Qds}={\displaystyle {\int }_A{q}_l}dxdy+{\displaystyle {\int }_A\frac{\partial w}{\partial t}}dxdy$$

(13)

The boundary conditions are:

$$-\frac{w^3}{12\mu }\frac{\partial p(x,y)}{\partial n}=Q$$

(14)

The flow rates within the fractures and other points on the wall are 0.

3. Galerkin Finite Element Method and Equation Solution

3.1. Galerkin Finite Element Method

The finite element method (FEM) divides the fracture body into some small parts (called finite elements), which are connected to each other through the designated nodes. For the fracture problem, because we do not know the real change of fracture height and fluid pressure in the fracture, we assume that the change of field variables in the finite element method is represented by a simple approximate function (shape function), which can be determined by the value of field variables at the node. When the equilibrium equation of the whole variable is determined, the system of equations can be solved (usually in the form of matrix), that is, the node value of the field variable can be obtained, and then the field variable value of the whole unit set can be determined by the approximate function. The steps of finite element solution: (1) discretization of the structure or solution domains; (2) selection of appropriate interpolation mode; (3) element analysis; (4) overall synthesis; (5) introduction of constraints; (6) solution of equation; and (7) calculation of other parameters [24].

Galerkin method adopts the weak form corresponding to differential equation. Its principle is that by selecting finite polynomial trial functions (also called basis functions or shape functions), superimposing them, and requiring the results to satisfy the original equation in solving the weighted integral in the domain and on the boundary (the weight function is the trial function itself)—a set of easy-to-obtain can be obtained. The linear algebraic equation is solved, and the natural boundary condition can be automatically satisfied [25].

According to the above method, the Galerkin finite element method is used to discretize Equations (4) and (12), and the linear binary polynomial (shape function) model is used to interpolate the equation. The fracture height and the pressure in the fracture are in the same form, and the potential function is equal to the shape function, $V={\varphi }_n$.

$$w(x^{\prime },y^{\prime })={\displaystyle \sum _{j=1}^N{\varphi }_i}(x,y)w_j$$

(15)

Equation (4) forms the following matrix equation:

$$[A]\{w\}=\{f\}$$

(16)

In the Equation (16),

$$A_{ij}=E{\displaystyle {\int }_A{\displaystyle {\int }_A\frac{1}{R}}[\frac{\partial {\varphi }_i}{\partial x}\frac{\partial {\varphi }_j}{\partial x^{\prime }}+\frac{\partial {\varphi }_i}{\partial y}\frac{\partial {\varphi }_j}{\partial y^{\prime }}]d{x}^{\prime }d{y}^{\prime }dxdy}$$

(17)

$$f_i=-{\displaystyle {\int }_A[p(x,y)-\sigma (x,y)]{\varphi }_{i}dxdy}$$

(18)

$$p(x,y)={\displaystyle \sum _{i=1}^N{\varphi }_i}(x,y)p_i$$

(19)

Equation (12) forms the following matrix equation:

$$[K]\{p\}=-\{f_l\}-\{f_w\}+\{f_p\}$$

(20)

In the Equation (20),

$$K_{ij}={\displaystyle {\int }_A\frac{w^3}{12\mu }}[\frac{\partial {\varphi }_i}{\partial x}\frac{\partial {\varphi }_j}{\partial x}+\frac{\partial {\varphi }_i}{\partial y}\frac{\partial {\varphi }_j}{\partial y}]dxdy$$

(21)

$$f_{li}={\displaystyle {\int }_A\frac{2c_l}{\sqrt{t-\tau (x,y)}}{\varphi }_{i}dxdy}$$

(22)

$$f_{wi}={\displaystyle {\int }_A\frac{\partial w}{\partial t}{\varphi }_{i}dxdy}$$

(23)

$$f_{qi}={\displaystyle {\int }_{\partial A_p}{q}_I{\varphi }_{i}ds}$$

(24)

3.2. Three Node Triangular Isoparametric Element

Because triangular elements have strong adaptability to complex boundaries, it is easy to discretize a two-dimensional region into finite triangular elements. The original curve boundary is approximated by several straight lines on the boundary. As the number of elements increases, the fitting becomes more and more accurate [26]. Considering the existence of a wellbore, the fracture area is divided into two main parts—the fluid flow area and the fracture front end area. The mesh of fluid flow region is sparse and the fracture front area is encrypted, which can be used to simulate the front of the fracture in detail [27]. The three node triangular element dissection is shown in Figure 2.

The equation s of the shape functions and their derivatives for triangular isoparametric elements are as follows [21]:

$${\varphi }_1=1-\xi -\eta , {\varphi }_2=\xi , {\varphi }_3=\eta$$

(25)

$$\xi =\frac{1}{|J|}[(y_3-y_1)(x-x_1)-(x_3-x_1)(y-y_1)]$$

(26)

$$\eta =\frac{1}{|J|}[-(y_2-y_1)(x-x_1)+(x_2-x_1)(y-y_1)]$$

(27)

$$\frac{\partial {\varphi }_i}{\partial x}=\frac{1}{|J|}(y_i-y_k), \frac{\partial {\varphi }_i}{\partial y}=\frac{1}{|J|}(y_k-y_j),$$

(28)

in Equation (28), i = 1, 2, 3; j = 2, 3, 1; k = 3, 1, 2 (Indicating subscript rotation [19]).

$$|J|=(x_2-x_1)(y_2-y_1)-(x_3-x_1)(y_2-y_1)$$

(29)

3.3. The Solution of Integral Equation

For the double area integral Equation (17), if the defining mode of the shape function $\frac{\partial {\varphi }_i}{\partial x}\frac{\partial {\varphi }_j}{\partial x^{\prime }}+\frac{\partial {\varphi }_i}{\partial y}\frac{\partial {\varphi }_j}{\partial y^{\prime }}$ is constant, the integral equation is transformed into ${\displaystyle {\int }_A{\displaystyle {\int }_A\frac{1}{R}}d{x}^{\prime }d{y}^{\prime }dxdy}$. According to the theory of Ioakimidis N.I. [28], the internal integral of the integral equation can be expressed as

$${\displaystyle {\int }_A\frac{1}{R}}d{x}^{\prime }d{y}^{\prime }={\displaystyle \sum _1^{3}sign(J_k)}{\displaystyle {\int }_{\Delta k}\frac{1}{R}}d{x}^{\prime }d{y}^{\prime }=sign(J_1)h_{1}In|\frac{s}{s-a}|+sign(J_2)h_{2}In|\frac{s}{s-b}|+sign(J_3)h_{3}In|\frac{s}{s-c}|$$

(30)

In the Equation (30), $sign(J_k)$ is the symbolic operation on the Jacobian $J_k$ value of a triangle. If it is greater than 0, the value is 1, if it is less than 0, the value is 1. $s=(a+b+c)/2$, $\mathrm{m}$; $a$, $b$, $c$ are the lengths of three sides, respectively, $\mathrm{m}$; $h_1=2{[s(s-a)(s-b)(s-c)]}^{1/2}/a$; $h_2=2{[s(s-a)(s-b)(s-c)]}^{1/2}/b$; and $h_3=2{[s(s-a)(s-b)(s-c)]}^{1/2}/c$.

In this way, ${\displaystyle {\int }_A\frac{1}{R}}d{x}^{\prime }d{y}^{\prime }$ is transformed into a constant term, and the external integral is transformed into a local coordinate according to the three-node triangular isoparametric element, and then the specific numerical value is obtained.

For the area integral Equation (18), the corresponding shape function expression is introduced, the local coordinates of x and y are transformed, and the triangular isoparametric element integration is carried out.

For the area integral Equation (20), the control condition for the unique solution of the equation is injection increment = filtration increment + fracture volume increment. The discrete form is as follows:

$${\displaystyle \sum _1^n{f}_{qi}}={\displaystyle \sum _1^n{f}_{li}}+{\displaystyle \sum _1^n{f}_{wi}}$$

(31)

For the area integral Equation (21), from the defining mode of the shape function $\frac{\partial {\varphi }_i}{\partial x}\frac{\partial {\varphi }_j}{\partial x}+\frac{\partial {\varphi }_i}{\partial y}\frac{\partial {\varphi }_j}{\partial y}$ is constant, and the integral term is simplified as ${\displaystyle {\int }_A\frac{w^3}{12\mu }}dxdy$. Similarly, the local coordinates are transformed according to the three-node triangular isoparametric element, and the specific values are obtained.

For the area Equation (22), according to the similar theory, the local coordinate transformation of x and y is carried out. Because there is a time term in the integral equation, in order to realize the executability of the algorithm, take $t=\Delta t$, $\tau (x,y)=0$ if this process is expressed as the initial state the fracture does not open, the fracture opens for a certain distance in time $\Delta t$, at this time, the time of beginning filtration in each point of the fracture is 0, and the total time of filtration is $\Delta t$. In the next time interval, the fractures continue to extend. The cumulative filtration time of the original open joint is $2\Delta t$, the cumulative filtration time of the new open joint is $\Delta t$, and so on until the end of construction.

For the area Equation (23), the finite difference $\frac{\partial w}{\partial t}$ processing is performed:

$$\frac{\partial w}{\partial t}=\frac{w^n-w^{n-1}}{\Delta t}$$

(32)

In the equation, $w^n$ and $w^{n-1}$ are the height of horizontal fracture at $t$ and $t-\Delta t$ time, $m$ respectively. Assuming that the initial time of the fracture is not open, the fracture inhales liquid, and the time interval is $\Delta t$, then the initial time $\Delta t$, $w^n$ is the height of the internal points of the horizontal fracture at the initial time interval. After the initial interval $\Delta t$, the integral Equation (23) is transformed into

$$f_{wi}={\displaystyle {\int }_A\frac{w_i}{\Delta t}{\varphi }_{i}dxdy}$$

(33)

For the area integral Equation (24), when the wellbore is considered, the integral value is the area of the longitudinal cylinder formed by the fracture height at the perforation position of the wellbore multiplied by the injection rate of fracturing fluid per unit area.

3.4. Picard Iteration Method

The basic idea of the iteration method has been proposed by previous scholars [29]: Firstly, assuming that the initial pressure of each point in the fracture is uniform, the fracture height of each initial point can be obtained using Equation (20). Using the initial height values, the pressure values at each point of continuous variation are calculated by Equation (16), and then the new fracture height values are obtained by Equation (20). Because of the equal initial pressure, this method leads to too many iterations in the later stage, which makes the calculation more complicated.

In this paper, a new Picard iteration method is proposed. The iteration parameters are changed sequentially. First, the distribution value of fracture height $w_{k,i}^{(n)}$ is assumed. Then, the pressure distribution value $p_{k,i}^{(n)}$ is calculated according to Equation (16). Then, the new distribution value of fracture height $w_{k+1,i}^{(n)}$ is obtained by bringing the obtained pressure $p_{k,i}^{(n)}$ distribution value into Equation (20). Then, the new distribution value of pressure distribution value $p_{k+1,i}^{(n)}$ is calculated according to Equation (16). Thus it completes an iterative process. The convergence conditions of iteration are as follows,

$$\frac{{\displaystyle \sum _1^n|p_{k+1,i}^{(n)}-p_{k,i}^{(n)}|}}{{\displaystyle \sum _1^n|p_{k+1,i}^{(n)}|}}<\epsilon$$

(34)

In the Equation (34), $\epsilon$ is the specified tolerance, 1 × 10⁻⁴.

Hypothesis of initial fracture height [19]: According to Sneddon’s research [30] the relationship between fracture height and normal stress acting on joint avoidance, it is assumed that a horizontal fracture with a circle is opened in an infinite homogeneous single elastic medium. If the normal stress acting on the fracture wall is a function $\sigma (r)$ the height of the fracture is determined below [30]

$$w=\frac{8(1-v^2)R^{\prime }}{\pi E}{\displaystyle {\int }_{\rho }^1\frac{\mu d\mu }{\sqrt{{\mu }^2-{\rho }^2}}}{\displaystyle {\int }_0^1\frac{\lambda \sigma (\lambda \mu R^{\prime })}{\sqrt{1-{\lambda }^2}}}d\lambda$$

(35)

In the Equation (35), $w$ is the height of the fracture at the radius $r$, $\mathrm{m}$; $v$ is the Poisson’s ratio of rock, decimal; $R^{\prime }$ is the maximum radius of the fracture, $\mathrm{m}$; $E$ is the elastic modulus of rock, $\mathrm{Pa}$; $\rho$ is the ratio $\frac{r}{R}$, less than 1; $\lambda \mu$ is the integral variable representing the fraction of the radius of the fracture.

When r = 0, the maximum height of the fracture is obtained. The relationship between the maximum height of the fracture and the normal stress acting on the fracture wall is obtained by simplifying and integrating Equation (34), $\rho =0$ [30]:

$$w_{\max }=\frac{8(1-v^2)}{\pi E}{\displaystyle {\int }_0^{R^{\prime }}\sigma (r)}{\cos }^{-1}\frac{r}{R^{\prime }}dr$$

(36)

As long as the normal stress distribution function on the fracture wall is known, this value is equal at all points of the fracture in this study, the maximum height can be obtained.

Assuming that the longitudinal section of the initial fracture is parabola, the variation of the fracture along the radial height is determined by the following equation [19]:

$$w=w_{\max }{(1-\frac{r}{R^{\prime }})}^{1/2}$$

(37)

So far, through the new Picard iteration method and coupling of Equations (16) and (20), after the initial time interval $\Delta t$ and the convergence of iteration, the values of fluid pressure and fracture height at each point in the fracture can be determined. Next, according to the stress intensity factor at the leading edge of the fracture, we can judge whether the fracture continues to extend.

3.5. Fracture Extension Judgement

The fracture criterion of linear elastic fracture mechanics is used as the criterion for judging the fracture propagation, that is, the stress intensity factor $K_I$ at the extension point is kept approximately equal to the critical stress intensity factor $K_{IC}$, and the stress intensity factor at any point on the fracture boundary is as follows [31]:

$$K_I=\frac{G}{2(1-\nu )}{(\frac{2\pi }{a})}^{1/2}{w}_a(x,y)$$

(38)

In the Equation (38), $a$ is a small distance from the leading edge of the fracture, m and $w_a(x,y)$ is the height of the fracture away from the leading edge of the fracture, $\mathrm{m}$.

If the calculated stress intensity factor is larger than the critical stress intensity factor $K_{IC}$, according to Mastorjannis [32], the distance of each point in the front of the fracture extending forward is as follows:

$$\Delta d=\frac{K_I-K_{IC}}{K_{IC}+\frac{\sigma w_a(x,y)}{\sqrt{h}}}$$

(39)

In the Equation (39), $\sigma$ is the crustal stress value at the front of the fracture, $\mathrm{Pa}$ and $h$ is the depth of the fracture, $\mathrm{m}$.

If the calculated intensity factor at the small distance of the fracture front is less than the critical intensity factor, the fracture does not expand, that is to say $\Delta d=0$.

After determining the position of the fracture front in the next time interval, the triangular mesh is reconstructed for the new fracture area. The height and pressure of the fracture at the previous moment at the node are calculated by linear interpolation. Then a symmetric positive definite stiffness matrix is formed for the coupled Equation (16) and Equation (20) and solved using the Picca iteration method. In such a cycle, until the end of construction, we can not only simulate the hydraulic fracturing horizontal fracture full three-dimensional expansion process, determine the pressure and fracture height at each point of the fracture at different times, as well as the value of the fracture extension radius, but also can analyze the impact of different factors.

4. Applications

4.1. Example Analysis

Daqing Changyuan Oilfield is situated around the anticlinal structural belt in the central depression of Songliao Basin, which is the most favorable area for oil generation and storage. It is the largest Mesozoic and Cenozoic sedimentary basin in Northeast China. The main sedimentary rock series are Cretaceous and Jurassic of fault-depression deposits in large sags, with an area of about 24 × 10⁴ km².

The Saertu structure is a short-axis anticline in the northern part of Daqing Placanticline. Its axis is NE22 degrees, its short axis is 14–22 km, its north width is narrow, its long axis is 32 km, its closure area is 439 km², its closure height difference is 600 m, its east wing inclination angle is 2–5 degrees, and its west wing inclination angle is 3–14 degrees. The Sazhong Development Zone is located in the central part of the Saertu structure, with an area of 161.25 km². There are 82 top faults of Pu I formation in the whole area, most of which are distributed in the west, and most of the strike are NW and NNW, all of which are normal faults.

Fault No. 1, in the north, is located in the Central Development Zone of the Saertu structure, with three rows from north to north, three rows from south to middle, faults from west to 98 degrees, and transitional zone from east to east. The North-first fault is located from the top to the east wing of the Saertu anticline, and the experimental area is located at the top of the anticline. The structure is gentle. The dip angle of the stratum is 1–2 degrees. There are no faults in the east-third fault area of the North-first fault.

The Sazhong Development Zone was put into development in 1960. Five development wells have been deployed successively. The average oil–water interface of each reservoir is 1050–1060 m above sea level. The reservoir has a unified pressure system with inactive edge and low water and no interlayer water.

The three types of reservoirs mainly consist of delta front facies deposits, which cover many types of sand bodies. It mainly includes thick and stable outer front sheet sand, thin and stable outer front sheet sand, thin and unstable outer front sheet sand and outer sheet sand.

A production well (PRO-1) in the West Fault Block of North 1 Area of Changyuan Oilfield, Daqing Oilfield was selected as a case study, as shown in Figure 3. The effective thickness of perforated interval in this well is 10.4 m, which is divided into six layers for stratified mining. The thickness of each layer of sandstone is less than 2 m. It is a typical thin sandstone and mudstone interbedding.

Saertu formation water is of sodium bicarbonate type. The salinity is 5611 mg/L, the chloride ion content is 1870.8 mg/L. The formation temperature of Saertu is 42.4 degrees, the viscosity of crude oil is 8.2 mPa·s, and the original gasoline ratio is 46.6 m³/t. Crude oil belongs to wax type. Its density is 0.856 g/cm³, wax content is 29.61%, gum content is 16%, and solidification point is 28 °C. The specific gravity of natural gas is 0.667 g/cm³, in which the content of methane is 86.37%, the content of hydrocarbons above ethane is 12.71%, and the content of carbon dioxide is 0.92%.

Three types of reservoir fracturing in Daqing Changyuan Oilfield adopt a new type of fracturing fluid system. Its main components include surfactants and alkalis. This fracturing fluid system has the characteristics of low viscous and Newtonian fluid flow. At the same time, by entering the formation, it can effectively drive the formation crude oil [33]. The fracturing process does not add proppant, so that the fracturing fluid fully enters the formation. By means of making a horizontal fracture, a large number of such fracturing fluids are infiltrated into reservoirs, thus expanding the swept volume of flooding fluids and enhancing recovery. The parameters of reservoir elastic modulus, Poisson’s ratio, vertical stress, comprehensive filtration coefficient, construction displacement, viscosity of fracturing fluid, and filtration coefficient are set up comprehensively. Model were calculated parameters as shown in

Table 1. Calculation of model parameters.

Model Parameter	Numerical Value	Reference Range
Reservoir elastic modulus/Pa	12 × 109	10 × 109~40 × 109
Poisson’s ratio of interlayer	0.26	0.24~0.26
Vertical stress/Pa	23.2 × 106	22 × 106~24 × 106
Construction displacement/(m3/min)	4.5	4~5
Comprehensive filtration coefficient of conventional fracturing/($m/\sqrt{\min }$)	1 × 10‒4	0.01 × 10−4~1 × 10−4
Comprehensive filtration coefficient of pressure drive technology/($m/\sqrt{\min }$)	10 × 10−4	10 × 10−4~100 × 10−4
Viscosity of fracturing fluid/P·s	2 × 10−3	1 × 10−3~5 × 10−3
Viscosity of conventional fracturing fluids/P·s	25 × 10−3	10 × 10−3~50 × 10−3
Construction fluid volume/m3	450	100~1000
Construction time/s	6000	1500~15,000

.

The relationship curves of horizontal fracture radius and time and maximum fracture height and time were obtained by simulation. At the end of construction, the relationship curves of fracture height and radius, net pressure, and radius in the center longitudinal section of the wellbore were obtained. The comparison of fracture radius and longitudinal maximum percolation distance between conventional fracturing and fracturing flooding processes under different filtration coefficients was obtained, as shown in Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8.

As shown in Figure 4, in the initial fracture stage, the fracture opens rapidly and expands. With the increase of time, the increase of the fracture radius decreases. The main reason is that with the increase in fracture radius, constant displacement injects fluid, the filtration of fracturing fluid increases, and the volume used for fracturing becomes less and less.

As shown in Figure 5, in the initial stage of fracture, the rock opens rapidly, and the maximum height of the fracture reaches 3.4 mm. With the increase of time, the height of the fracture increases slowly, with a small increase from 3.4 mm to 5.3 mm.

As shown in Figure 6, the height of each point in the fracture is different, and the height of the fracture in the wellbore is the highest, and gradually decreases to the front of the fracture.

As shown in Figure 7, the net pressure of each point in the fracture is different. The net pressure of the fracture in the wellbore is the largest, and the fracture closes gradually near the front.

As shown in Figure 8, under the same reservoir conditions and construction parameters, compared with conventional fracturing, fracturing flooding technology has a larger filtration coefficient, which results in a larger filtration rate, smaller fracture length, and larger vertical liquid infiltration distance.

Because the simulated parameters are the actual well construction parameters, and the well has been subjected to microseismic monitoring, the monitoring results are similar to the simulated fracture radius. Microseismic data show that the fracture propagation process is basically circular, and the average fracture radius during the propagation process is close to 31.8 m, as shown in Figure 9 and

Table 2. Calculation of model parameters.

Fracturing Horizon	Fracture Network Length (m)	Fracture Network Width (m)	Orientation of Fracture Network (°)	Number of Microseismic Events (number)
PI111~112	32.246	31.214	133.870	23

. Therefore, the applicability and accuracy of the method are verified.

4.2. Factor Analysis

Horizontal fracture fracturing is a very complex multifield coupling problem. The actual fracture morphology is affected by various factors, including geological factors and engineering factors. The inherent properties of the original reservoir cannot be changed artificially when geological factors are involved. In order to achieve a certain design purpose, we must control it effectively by changing the construction plan. The horizontal fracture of vertical well fracturing described in the above section is the basic model, and the influence of parameter changes on the expansion characteristics of horizontal fracture is discussed.

4.2.1. Rock Elastic Modulus

The elastic modulus of rocks in the target block varies from 10 × 10⁹ Pa to 40 × 10⁹ Pa (Thin Sandstone Reservoir), and in this range the elastic modulus of rocks is changed for fracturing simulation calculation, the results are shown in Figure 10. Obviously, the elastic modulus of rock has obvious influence on the plane shape of fractures. With the increase of elastic modulus, the fracture radius increases significantly, while the maximum fracture width decreases linearly. Generally speaking, under the same conditions, reservoirs with large elastic modulus are prone to produce fractures with large plane range and small fracture height, while reservoirs with small elastic modulus are the opposite.

4.2.2. Construction Displacement

Construction displacement is an important measure and means, which is controlled by human factors. The displacement range of the construction block is from 0.0167 to 0.0833 ${\mathrm{m}}^3/\mathrm{s}$. The total injected liquid volume is controlled and thus unchanged. Different displacement injection time is simulated. The comparison curves of fracture height, radius, and longitudinal maximum percolation distance are obtained as shown in Figure 11.

From Figure 11, it can be seen that the influence of construction displacement on fracture radius is not significant, and the increase in displacement cannot greatly increase the fracture radius. The increase of displacement will increase the height of fractures. In addition, the increase in displacement can effectively increase the longitudinal maximum infiltration distance of the fracturing fluid, that is, the sweep range of fracturing fluid (flooding agent), which is quite beneficial to enhance oil recovery. But when the displacement exceeds 4.0 m³/min, the improvement effect is reduced. The construction capacity of the site optimization design is 4.0 m³/min.

4.2.3. Viscosity of Fracturing Fluid

The viscosity range of fracturing fluid in the study block is ${15\times 10}^{-3} \mathrm{Pa}·\mathrm{s}$, and this range is simulated in our research. The effect of fracturing fluid viscosity on fracture morphology and filtration is shown in Figure 12.

From the comparative study above, it can be seen that when the viscosity of traditional fracturing fluid is 12.5 times that of the fracturing fluid in the fracturing flooding process—the fracture radius of traditional fracturing fluid is larger and the infiltration distance of fluid is smaller. However, in the target block, the range of viscosity variation of fracturing fluid used is small, and viscosity has little effect on fracture radius and fracture height, but has a great influence on vertical filtration distance. In order to further enlarge the sweep volume of fracturing fluid, smaller viscosity should be chosen.

5. Results and Discussion

In this paper, the full three-dimensional mathematical and physical equation of horizontal fracture in hydraulic fracturing is given, and the equation is solved by the finite element method and a new iteration method. The fracturing and percolation process of three types of reservoirs in Daqing Chang-yuan Oilfield is simulated, which proves that this process can greatly increase the swept volume of flooding fluid.

In the initial stage of the traditional method, it is assumed that the pressure in the fracture is equal, which leads to too many iterations in the later stage and increases computational complexity. The iteration theory proposed in this paper reduces this complex type, that is, the selection of initial iteration value problem, so that the convergence and accuracy of the calculation can be improved rapidly.

The results of understanding the horizontal fracture propagation process include: in the initial stage, that the radius and height of the fracture increase sharply, and the increase slows down in the later stage; the height and net pressure of each point in the fracture are unequal, the difference is largest in the wellbore, and gradually decreases to the front of the fracture. When the elastic modulus is larger, the radius of fracture is longer and the height is smaller. With the increase of displacement, the radius of fracture does not change significantly, the height increases, and the maximum filtration distance increases. It is suggested that the displacement of construction should be 4.0 m³/min. The viscosity of the target block is small, and smaller viscosity should be chosen to further enlarge the sweep volume.

Our research work is aimed at the horizontal fracture of hydraulic fracturing, when the vertical stress is less than the minimum horizontal principal stress, it is easy to form horizontal fracture. In most cases, the vertical stress of the shallow reservoir is compounded above the law. The horizontal fracture propagation studied in this paper is suitable for shallow reservoir.

Because of the hypothesis of numerical simulation, only elastic modulus, Poisson’s ratio, shear modulus, and fracture toughness are considered. For different lithological reservoirs, these values are different. If these rock mechanics parameters of different lithological reservoirs are determined, simulation can be carried out. But the actual results are also affected by many other factors, such as reservoir heterogeneity, lithological mutation, reservoir bending, fluid solid coupling, etc. Therefore, this model is a relatively ideal model, and its applicability needs further modification of the parameters or requires a new object equation to be introduced. This paper does not carry out this work.