Numerical Investigation of Hydraulic Fracture Extension Based on the Meshless Method

The fracture propagation in hydraulic fracturing is described as a nonlinear problem dynamic boundary. Due to the limitation of mesh refinement, it is difficult to obtain the real crack propagation path using conventional numerical methods. Meshless methods (MMs) are an effective method to eliminate the dependence on the computational grid in the simulation of fracture propagation. In this paper, a hydraulic fracture propagation model is established based on the element-free Galerkin (EFG) method by introducing jump and branch enrichment functions. Based on the proposed method, three types of fracturing technology are investigated. The results reveal that the stress interference between fractures has an important impact on the propagation path. For the codirectional fracturing simultaneously, fractures propagate in a repel direction. However, the new fracture is attracted and eventually trapped by the adjacent fracture in the sequential fracturing case. For the opposite simultaneous fracturing in multiwells, two fractures with a certain lateral spacing will deflect toward each other. The effect of stress shadow should be used rationally in the optimization of construction parameters; for the single well multistage fracturing, the stage spacing should be out of stress inversion area, while for the simultaneous fracturing of multiple wells, stress inversion zones should be used to maximize communication between natural fractures. Overall, this study establishes a novel and effective approach of using MM to simulate the propagation of hydraulic fractures, which can serve as a useful reference for understanding the mechanism of hydraulic fracture propagation under various conditions.


Introduction
Hydraulic fracturing is a vital technology for the stimulation of unconventional tight reservoirs with low permeability and complex pore-throat system [1][2][3][4][5][6]. Complex fracture network is of great significance for economic exploitation of shale gas. The prediction of hydraulic fracture propagation paths plays an important role in fracturing design; consequently, it has received considerable research attention [7][8][9]. In recent years, except for physical experiments [10,11], various numerical methods have been adopted by researchers to investigate this problem, including finite element method (FEM), displacement discontinuity method (DDM) [12], extended finite element method (XFEM) [13], and phase field method (PFM) [14].
Hunsweck et al. [15] proposed a FEM-based algorithm to examine the effect of fluid lag during the fracturing process by considering the nonlinear coupling between the fluid pressure and fracture opening. A multiple fracture propagation 3D model based on the FEM established by Guo et al. [16], the inhibit effect on the fracture length and width produced by the stress interference effect between hydraulic fractures was presented. In addition, Shauer and Duarte [17] proposed an improved generalized finite element method (GFEM) to simulate the 3D hydraulic fracture propagation, and the results suggested that this algorithm could decrease the number of Newton iterations and increase the robustness. As one of the methods based on FEM, the cohesive zone method (CZM) has also been widely applied for the simulation of hydraulic fracture propagation [18][19][20]. Its main idea is to embed the cohesive zone element into the region between solid elements to simulate the fractures. The main disadvantage of this method is that the propagation path of the fracture is predefined, and the extension direction of the fracture cannot be changed according to the real-time stress state. Meanwhile, FEM depends on the mesh structure. Especially after the crack propagation, the mesh needs to be restructured, and its quality seriously affects the accuracy of the simulation results [21][22][23].
To better solve the problem of the discontinuous structure in the traditional FEM. Belytschko and Black [13] introduced an enrichment function that could describe the discontinuous structure in the shape function of the conventional finite element, so the requirement in the conventional FEM that the mesh must coincide with the boundary of the discontinuous structure was avoided, and the grids describing the crack could be independent of other grids. This is the core idea of XFEM [24]. However, in this method, additional equations are required to describe the fluid flow equations in fractures and reservoir filtration for simulating the hydraulic fracture propagation. Based on XFEM, Taleghani and Olson [25] introduced lubrication equation to describe the fluid flow in the fracture and examined the interaction between hydraulic fracture and natural fracture. Suo et al. [26] considered the effect of fluid leak-off and mixed-mode failure to investigate the interaction between hydraulic and natural fractures. Similarly, Zheng et al. [27] used XFEM to study the interaction between hydraulic and natural fractures, and the results showed that the approach angle and stress difference had a predominant effect on the reactivation of natural fractures as compared to the other parameters. In addition, it has been proposed to combine XFEM with CZM or other methods to simulate fracture propagation [28][29][30].
As one of the boundary element method (BEM), DDM is widely used to solve the fracture extension problem, especially for staged multicluster fracturing [31][32][33]. The advantage of this method is that only the hydraulic fractures need to be discretized. Based on this method, Zhou et al. [4] proposed a factor "H" to evaluate the relationship between the initial angle and deflection angle of fracture. It may be noted for a reservoir containing natural fractures, the path of hydraulic fractures can be disturbed by natural fractures, and the natural fractures may open and shear due to stress interference from hydraulic fractures before hydraulic fractures intersect with natural fractures [3]. Based on DDM method, Ren et al. [34] simulated the nonplanar propagation of hydraulic fractures and obtained the stimulated reservoir volume (SRV) of tensile or shear failure zone. Wu and Olson [35] suggested that the mechanical interaction between cracks was underestimated in the pseudo-3D DDM proposed by Olson [36], and a modified method called fully 3D DDM was proposed to overcome this problem. The results showed that the extra elements generated by the height of cracks could be eliminated. However, the filtration of liquid was not considered in this method, and it was only suitable for linear elastic and homogeneous materials. In addition, it was difficult to describe the variation of fluid pressure and its distribution at the intersection point.
The basic idea of PFM is to introduce a field order parameter or fracture field to describe the smooth transition between complete and incomplete states to solve the fracture problem. In other words, the sharp fractures are transformed into diffused cracks by introducing an auxiliary phase field variable [37]. The main advantage is that no additional rules are required to track the fracture interface, and the grid need not be restructured after fracture extension [38][39][40].
As pointed by Lecampion et al. [37], meshless methods (MMs) are effective for solving the problem of fracture propagation in elastic and porous-elastic [41]. The MMs are primarily divided into two categories: element-free Galerkin (EFG) and smoothed particle hydrodynamics (SPH) methods. The EFG method was originally developed by Belytschko et al. [42], and its main idea is to use moving least squares method for the approach functions. To solve the discontinuous boundary problem, Belytschko et al. [43] proposed to use a jump function for the displacement discontinuity along the crack faces and the Westergard's solution enrichment around the crack tip. The main advantage of this method is that it does not require any connectivity data between the nodes; in addition, the accuracy of the results can be ensured even if the node arrangement is irregular. However, there are only few studies on the simulation of hydraulic fracture propagation based on MMs. Oliaei et al. [44] simulated the propagation path of hydraulic fractures in homogeneous and heterogeneous saturated soils using the EFG method to verify the feasibility of this method. Furthermore, Wen et al. [45] used the 3D EFG method to simulate the fracture propagation in coalbed methane (CBM). In this model, the fluid pressure distribution in the fracture was considered, and the simulation results were consistent with the experimental results.
In this study, a hydraulic fracture propagation model based on the EFG method is established. Further, the evolution of deflection angle and the effect of spacing on the fracture propagation are investigated. Compared to the existing studies on the simulation of fracture propagation using MM, the morphology of multiple fractures occurring simultaneously or sequentially are presented.

Governing Equations of Hydraulic Fracturing Problem.
As shown in Figure 1, for the hydraulic fracture in the control region Ω, the fracture surface boundary, the external force boundary, and the displacement constraint boundary can be expressed as Γ c , Γ t , and Γ u , respectively. The model established in this paper is based on the following assumptions: (1) The solid model is isotropic in two-dimensional (2D) space (2) The rock deformation is considered as a small linear elastic deformation (3) The fluid in the fracture is considered as an incompressible Newtonian fluid, and the effect of fluid filtration is not taken into account 2 Geofluids Without considering the influence of the body force, the equilibrium differential equation can be expressed as follows: The boundary conditions of the equation can be written as According to the principle of virtual displacement, the equilibrium equation can be obtained as follows: whereδu and δv represent the virtual displacement in the x and y directions, respectively. Using the partial integral formula, we obtain ð Ω σ x δε x + τ xy δγ xy + σ y δε y dΩ = ð Γ σ x n x + τ xy n y À Á δu Â wheren x = cos ðn, xÞ andn y = cos ðn, yÞ.ε x , ε y , and γ xy represent the strain in the x direction, y direction, and the shear strain, respectively. On the force boundary of the elastomer, f x and f y represent the surface forces along the x and y directions, respectively. Combined with the elastic mechanics equation, equation (4) yields the following: ð 2.2. Fluid Flow within the Fractures. The flow of fracturing fluid within the hydraulic fracture is considered as a flow between parallel plates, which satisfies the following cubic law: where q is the volume flow rate along the crack direction, t is the fracturing time, w denotes the crack opening, p f represents the fluid pressure in the fracture, and μ is the viscosity of the fracturing fluid.
Assuming that the fluid has no filtration loss and is incompressible in the fracture, the material balance equation in the fracture can be expressed as Because the filtration of the crack is not considered, the total material balance equation is where L is the length of the crack, h represents the height of the crack, and Q 0 is the volume flow of injection.
In addition, at the tip of the fracture, the flow rate of the fracturing fluid is zero, and the crack width is also zero. Thus, the boundary conditions can be expressed as follows: The crack width, which is formed by the displacement of the crack wall surface (as shown in Figure 1), can be expressed as where n is the normal unit vector outside the crack surface, u + up represents the displacement of the upper surface of crack, and u − bot is the displacement of the bottom surface of crack. The displacement of the fracture surface must be calculated by MM. When the fluid and stress are coupled, the fluid  [43] proposed to use the partition of unity (PU) enhanced method to obtain a new displacement function for MM. The expression can be divided into two parts: continuous displacement and discontinuous displacement. The approximate displacement function can be expressed as follows: where I is the node set, i = 1,2,3, ⋯, n; n represents the total number of nodes; S represents the general node set; u iI represents the node displacement set; S c is the step strengthening point set; S f denotes the crack tip strengthening point set; N is the shape function; α iI and β iI j represent unknown vector coefficients. H is the Heaviside function, which is defined as Branch enrichment function B j ðxÞ indicates the displacement function of the crack tip. It can be written as follows: where r and θ represent the coordinates of the crack tip in the local polar coordinate system.

Shape Function.
The displacement and stress in a solid field can be calculated by establishing field nodes based on MM, where the mesh dependency of the traditional method can be eliminated. However, the shape function generated using the traditional method may not be suitable for MM as there is no connectivity between the nodes. Therefore, a new method is needed to generate the shape function for MM. To this end, a circular domain of influence is introduced, as shown in Figure 2 (yellow line). For this case, the radius of influence region is 1.7 times the maximum spacing between the field nodes. The shape function is generated by using the moving least squares method due to its high stability in EFG [46]. Considering the displacement variable function uðxÞ in a continuous domain, the moving least squares approximation at x can be expressed as follows: where pðxÞ is the basis function in 2D space coordinate system, and aðxÞ is an undetermined parameter. A monomial is usually adopted as the basis function.
A reasonable value of aðxÞ can make the approximate function as the optimal solution of uðxÞ function in one neighbourhood of the calculation node. To obtain this value, the weight function is defined on each node, and the error estimation function is introduced. Finally, the shape function can be expressed as follows: 2.3.3. Weight Function. The selection of weight function is extremely important for the generation of shape function. Usually, they are expressed using 3 rd order and 4 th order spline functions. Liu et al. [47] provided a general expression to construct any continuous-order weight function. To meet the corresponding conditions, the 4 th order spline function can be expressed as follows: where x − x i represents the distance between the nodes x and x i , and r w represents the radius of domain of influence. (5) can be written in vector form as follows:

Overall Stiffness Matrix. Equation
Equation (18) is the integral form generated by the where In contrast to the FEM, the expansion term brought by enriched degree of freedom is included in the B matrix, which is expressed as where B std and B enr are written as follows: Here, φ i is the branching function of crack tip. The specific formula is shown in equation (13).

Fracture Extension Criteria.
In this paper, the maximum circumferential stress theory is considered as the fracture propagation criteria. For the 2D mixed-mode fracture, the stress component at the crack tip is expressed in the polar coordinate system as follows [48]: The fracture propagation occurs when σ θθ max reaches the fracture toughness, and the deflection angle is expressed as follows [49]: where K I and K II represent the stress intensity factor (SIF) for the opening mode and shear mode, respectively; K IC is the fracture toughness, and θ 0 is the deflection angle.

Computational Flow
Based on the theory mentioned above, the computation flow is presented in Figure 3. The main calculation process is as follows: (i) The information on freedom of nodes generated (ii) As shown in Figure 4, the nodes are classified into three types. A circular region is formed around the crack tip with a radius of γ a . All the nodes in this region are classified as crack tip enhanced nodes. The nodes located on both sides of the crack within the minimum distance γ a from the crack are marked as step enhanced nodes. The remaining nodes are classified as normal nodes.
(iii) The sparse matrix form is adopted to store the overall stiffness matrix (iv) Similar to FEM, the displacement boundary is established by adopting the Lagrange multiplier method (v) Similar to the finite element, the displacement boundary was loading by adopting the Lagrange multiplier method (vi) The singular value decomposition (SVD) or QR decomposition method is employed to obtain the stress and displacement fields

Model Verification
In this section, the SIF and stress distribution during fracture propagation obtained using MM are compared with those reported in the literature to validate the feasibility of the proposed model. 6 Geofluids the SIF is independent of the length when the length is greater than three times the width, and it can be expressed as follows [48]: where a represents the length of crack, and L and b are the width and the length of plate, respectively; σ represents the tensile stress on the plate. Table 1 Table 2. Figures 6(a) and 6(b) show the induced stress distribution along and perpendicular to the fracture direction, respectively. The shear stress distribution is presented in Figure 6(c). To present a clear comparison of stress between the proposed and existing models, the stress along the wellbore is presented in Figure 7. It is obvious from Figures 6(c) and 6(d) that the induced stress components are consistent with the DDM results reported by Ren et al. [34]. As shown in Figure 7, dimensionless stress is adopted for a clear comparison. Similarly, σ xx , σ yy , and τ xy are in good agreement with the DDM results in the direction perpendicular to the crack.

Basic Parameters.
To verify the practicality and reliability of the proposed method, three typical hydraulic fracturing examples are established based on MM, including codirectional simultaneous fracturing, sequential fracturing, and opposite simultaneous fracturing, as shown in Figure 8. In the first case, two fractures extend simultaneously. In the second case, fractures propagate in sequence from bottom to top. In the third case, two fractures extend toward each other simultaneously. The law of multifracture interaction during the hydraulic fracturing process is investigated by using MM. The basic input parameters are shown in Table 2. The direction of maximum principal stress is parallel to the hydraulic fracture. Figure 9 shows the displacement distribution of two simultaneously propagating fractures. The stress interference between the fractures is weak during the initial period of fracturing, while it becomes strong later. The reason is that with the injection of liquid, the fluid pressure in the fracture gradually facilitates the stress interference between the fractures. In addition, the influence of cluster spacing and stress difference on crack propagation trajectory is also studied. The variation in the deflection angle of the crack tip is presented in Figures 10  and 11 under the cluster spacing of 10, 20, and 30 m. It is clear that the deflection angle gradually decreases with the increase in the crack spacing, while it increases with the increase in time. The reason is that with the increase of spacing, in the initial time, the stress interference has not spread to the adjacent fracture, and the intensity of the interference between fractures is relatively small. As the liquid continues injected, the fractures gradually propagate from each other.

Codirectional Simultaneous Fracturing.
The impact of stress difference (3, 5, and 7 MPa) on the deflection angle of crack is shown in Figures 12 and 13. It can be seen that the deflection angle of fracture decreases with the increase in stress difference. Obviously, during the early stage, the deflection angle increases rapidly with the increase in time but tends to decrease later until a plateau is reached. This is because the effect of stress interaction on the fractures becomes weak after the fracture is extended to a certain distance. There have two explanations for this  phenomenon. One is that the increases of fluid pressure required for the fracture deflect increase as the stress difference becomes larger. When the pressure required is higher than that along the original direction, the deflection tendency of the fracture would become weak according to the principle of least resistance; another reason is that the increase of stress difference makes it more difficult to form stress inversion in the area around the fracture.
Previous studies have suggested that the reservoirs with large stress difference generally exhibit poor fracability. However, from another perspective, in a reservoir with a higher horizontal stress difference, the fracture can extend to a long distance under uniform propagation. Therefore, the reservoir can shatter far away from wellbore. Other studies have shown that for the reservoirs with well-developed natural fractures, the natural fractures can easily intersect with hydraulic fractures under large stress difference. However, it is not beneficial for stimulating the natural fractures in the reservoirs and for maximizing the SRV of reservoirs. Therefore, the con-struction parameters must be optimized based on the development and distribution of natural fractures.

Sequential Fracturing.
Over the recent years, the average number of fracturing stages per well has increased continuously to boost the production from unconventional reservoirs [50]. However, single-stage fracturing is generally used due to the limitations of the pump pressure and injection rate. In addition, Wang [51] suggested that the fracture propagation behavior can change from complex to simple, and the number of main fractures decreases when the propagation distance from the wellbore increases. This leads to a stage that is fractured first, which causes stress interference on the later stage that affects the crack extension path. Therefore, it is important to accurately predict the extension of back pressure fractures in sequential fracturing. Figure 14(a) shows the displacement contour when the first fracture is initiated. Figure 14(b) presents the second fracture propagation after the first fracturing process is (d) Figure 6: Comparison of dimensionless stress contour between our model and literature [34]: (a, c) direction of maximum principle stress; (b, d) direction of minimum principle stress.

Geofluids
finished. Obviously, the second fracture is attracted by the first fracture and eventually coalesces with the first fracture. This phenomenon is quite different from that reported earlier in which the new fracture deviate from the first fracture. Earlier studies have suggested when the second fracture occurs, it is difficult to open it along the original direction due to the additional induced stress in the matrix around the initial fracture, and the new fracture extends away from the adjacent fracture. However, our results indicate that the phenomenon of fracture fusion may occur in one wellbore. This has not been demonstrated earlier. This is because there is a stress inversion region around the fracture, which implies that the new fracture extends along the new direction of maximum horizontal principal stress, causing it to be attracted by the previous fracture and intersect it finally. This explains why some of the fracturing stages have almost no contribution to production according to the field test [52]. Accordingly, the fracture propagation patterns with stage spacings of 50, 60, and 70 m are simulated. As shown in Figure 15, the tendency of fractures to merge becomes weak as the stage spacing increases. The fractures exhibit a trend of uniform propagation when the stage spacing is larger than 70 m. Further, according to the vector distribution of maximum principal stress, this distance is just outside the stress inversion region. Therefore, it can be concluded that in the current staged sequential fracturing technology, although tight stage spacing leads to a higher initial production, but it is not conducive to improve the long-term cumulative production.

Geofluids
In general, the stage spacing should be outside the stress inversion zone of the fractured fractures.

Opposite Simultaneous Fracturing.
It is generally believed that the shale gas wells exhibit a high production rate in a short period, but the stable production rate is lower during the entire production life. Therefore, the simultaneous fracturing of horizontal infill wells is usually recommended. For this case, the fracture paths under the lateral spacing of 0 and 15 m are shown in Figure 16. When the spacing is 0, the two fractures do not repel but merge into one fracture finally. When the spacing is 15 m, the two fractures attract each other, and this behavior is consistent with the report of Zhou and He [53]. Based on this observation, the cases with different lateral spacing are investigated, and the evolution of deflection angle as a function of time is shown in Figure 17. Initially, the fractures propagate with a small deflection angle because the stress interference region still exists even though the distance is larger, and the fractures repel slightly. This is similar to the case of simultaneous propagation of multiple fractures. However, the fractures gradually appear close to each other, and the tip starts to deflect when the influence regions of two fractures begin to overlap. Finally, they return to the direction of maximum horizontal principal stress. This is reflected in the evolution of deflection angle under the lateral spacing of 40 m.
The attraction between fractures can be attributed to the fact that stress inversion occurs in this region due to the superposition of induced stress of two cracks when the fractures begin to overlap. While implementing simultaneous fracturing of adjacent wells in the field, it is important to appropriately utilize the stress inversion area. When stress inversion occurs in the shale reservoir containing rich natural fractures, natural fractures that remain closed under the original in situ stress state may suffer shear failure even without proppant, and a new self-propping SRV is formed. The stress inversion area becomes smaller when the lateral spacing is small, while it is difficult to form the stress inversion area under extremely large spacing. Therefore, there is an optimal value of lateral spacing, which also depends on the development and distribution of natural fractures in the reservoir.

Conclusions
A 2D fracture propagation model was established based on MM. The validity of the proposed model was verified through the calculation of SIF and stress distribution. Based on this method, the nonplanar propagation of multiple fractures was investigated. The main results of the study are summarized as follows: (1) It was proved that the MM with coupled fluid flow is suitable for the simulation of hydraulic fracture propagation. The interdependencies between grids can be eliminated by using this method (2) For the case of codirection simultaneous propagation, hydraulic fractures propagate from each other due to the stress interference. The decrease of spacing and reservoir stress difference will promote fracture interaction (3) During sequential fracturing, the subsequent fracture was attracted by the adjacent fracture due to stress reversal, and this phenomenon has been rarely reported in the existing literature. Tight stage spacing caused overlapping of some stages, which had almost no contribution to the production. For the single well fracturing, the stage or cluster spacing should be out of the stress inversion area (4) For the case of the opposite simultaneous fracturing in multiple wells, two fractures with certain spacing were easily attracted to each other due to stress reversal. Overall, it is necessary to appropriately utilize the stress shadow effect, and too small or too large lateral spacing may fail to maximize the shear failure of natural fractures in the overlap region

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declared that they have no conflicts of interest. 14 Geofluids