Simulation of Hydraulic Fracturing Using Different Mesh Types Based on Zero Thickness Cohesive Element

Hydraulic fracturing is a significant technique in petroleum engineering to enhance the production of shale gas or shale oil reservoir. The process of hydraulic fracturing is extremely complicated, referring to the deformation of solid formation, fluid flowing in the crack channel, and coupling the solid with fluid. Simulation of hydraulic fracturing and understanding the course of the mechanism is still a challenging task. In this study, two hydraulic fracturing models, including the Khristianovic–Geertsma–de Klerk (KGD) problem and the hydraulic fracture (HF) intersection with the natural fracture (NF), based on the zero thickness pore pressure cohesive zone (PPCZ) element with contact friction is established. The element can be embedded into the edges of other elements to simulate the fracture initiation and propagation. However, the mesh type of the elements near the PPCZ element has influences on the accuracy and propagation profile. Three common types of mesh, triangle mesh, quadrangle mesh, and deformed quadrangle mesh, are all investigated in this paper. In addition, the infinite boundary condition (IBC) is also discussed in these models. Simulation indicates that the results of pore pressure for zero toughness regime simulated by the triangle mesh are much lower than any others at the early injection time. Secondly, the problem of hydraulic fracturing should be better used with the infinite boundary element (IBE). Moreover, suggestions for crack intersection on the proper mesh type are also given. The conclusions included in this article can be beneficial to research further naturally fractured reservoirs.


Introduction
There are a large number of shale gas reserves in the world. They usually exist in the form of adsorption state and free state under the rock formation. America is the pioneer of shale gas engineering and makes a great breakthrough. In 2011, it was estimated that the production of shale gas in America accounted for 21.69% of the output of the total natural gas [1]. Since the huge success is achieved, shale gas and shale oil are widely known around the world, including China, Europe, India, etc. [2,3]. In the meantime, the utilization of shale resources cannot just relieve the pressure on the yield of crude oil but also reduce the pollution induced by fossil fuel. However, shale gas or shale oil is difficult to be exploited for engineers because the permeability of the shale formation is extremely low. It is reported that only a few wells that have perfect natural fracture format can be brought into production and 90% of wells need reservoirs reconstruction such as acidizing and fracturing in order to obtain the ideal yield [4]. The cohesive zone model (CZM) is the predominant tool in the FEM, which dates from the 1960s. Barenblatt [24] and Dugdale [25] originally introduced this model and investigated fracture propagation in the brittle materials and ductile materials with the small scale of plasticity, respectively. CZM postulates that there is a damage evolution zone ahead of the crack tip, which avoids the singularity at the crack tip in the linear elastic fracture mechanics (LEFM). Moreover, the extra definition of the crossing criterion is not specified under the framework of the cohesive method. Additionally, CZM can solve the problem that the lubrication equation contains a degenerate nonlinear partial differential equation that poses a considerable challenge for numerical modeling [26]. Owing to these properties, a variety of simulations of hydraulic fracturing based on CZM occur rapidly. Chen [27] realized the process of single fluid-driven fracture propagation in the ABAQUS software and validated its correctness by virtue of comparing it with the zero toughness solution. Guo [28,29] established a coupling stress-seepage-damage model to simulate the interaction between hydraulic fracture and natural fracture without crossing the criterion, as well as study the effect of approaching angle and geologic parameters on the penetration of hydraulic fracture. Nguyen [30] provided a zero thickness flow cohesive interface element and inserted it into the rock formation, and analyzed the relationship between the order of element and modeling issues such as convergence and numerical integration. Li [31], Wang [32], and Xiang [33] started to develop global zero thickness pore pressure cohesive zone (PPCZ) elements in order to carry on the fluid-driven fracture propagation or intersection between HF and NF without predefining the propagation path.
Although many scholars succeeded in investigating hydraulic fracturing for either single crack or multiple cracks, they seldom illuminated the reason they utilized a triangle or quadrangle mesh and their advantages or drawbacks. Furthermore, the infinite boundary condition (IBC) played an incredibly important role in fracture mechanics, most of the asymptotic solutions were founded on the infinite boundary. However, IBC was hardly mentioned in previous works [27]. In this paper, a hydro-mechanical model, which is based on the zero thickness PPCZ element implemented into the UEL of ABAQUS, is established to simulate the HF propagation and intersection with NF. Three types of mesh shapes are discussed in this study, and the influence of IBC on the results is presented as well. Finally, the systematic suggestion will be given for further researches in the future.

Methodology of Pore Pressure Cohesive Zone Model
Since the coupling problem is very complex leading to a challenging simulation, for simplicity, some following assumptions will be made in this model: (1) Isotropic rock media; (2) small rock deformation; (3) two-dimensional plane strain problem; (4) incompressible fluid and laminar flow; (5) no fluid lag, no leakoff, and no fluid gravity. According to the achievement of Adachi, zero fluid lag is reasonable if the confining stress is huge [34].

Governing Equation of Solid Medium and Fluid in the Fracture
Problems of hydraulic fracturing are the coupling process, which consists of the solid part and the fluid part. Considering the rock matrix is an elastic and homogeneous material, the momentum balance equation of the formation can be described as: where σ is the total stress tensor, ρ s is the density of the solid phase, b and u are the body force and solid displacement vector, respectively. The solid deformation can be expressed by Hooke's law: where D is the material parameter matrix, ε is the strain vector whose component can be written as ε ij = 1 2 u i,j +u j,i . Figure 1 presents the force condition of the fracture surface and fracturing liquid. p f (s,t) represents the fluid pressure, Q is the injection rate, s is the coordinate of the crack; w(s, t) and l(t) are the width and half length of the crack, respectively. σ v and σ h are in-situ stresses representing a vertical and horizontal direction. In terms of Terzaghi's theory [35], the behavior of the solid is related to effective stress. Therefore, the effective stress σ ij around the fracture is given by: where δ ij is the Kronecker delta function, α is the biot's coefficient that can be equal to the unit with regard to incompressible solids.
Processes 2020, 8, x 4 of 20 Figure 1 presents the force condition of the fracture surface and fracturing liquid. pf(s,t) represents the fluid pressure, Q is the injection rate, s is the coordinate of the crack; w(s, t) and l(t) are the width and half length of the crack, respectively. v and h are in-situ stresses representing a vertical and horizontal direction. In terms of Terzaghi's theory [35], the behavior of the solid is related to effective stress. Therefore, the effective stress  ij σ around the fracture is given by: where δ ij is the Kronecker delta function, α is the biot's coefficient that can be equal to the unit with regard to incompressible solids. Fluid in the fracture should satisfy the mass conservation equation that can be given by where f ρ is the fluid density, q is the volume fluid rate.
Since fluid flow is regarded as laminar flow, the discipline of flow inside the crack can be described by the Poiseuille equation between two parallel crack surfaces. The equation is also called cubic law: where μ is the fluid viscosity. This equation can be considered as the tangential flow in the crack. Inserting Equation (  Fluid in the fracture should satisfy the mass conservation equation that can be given by where ρ f is the fluid density, q is the volume fluid rate.
Since fluid flow is regarded as laminar flow, the discipline of flow inside the crack can be described by the Poiseuille equation between two parallel crack surfaces. The equation is also called cubic law: where µ is the fluid viscosity. This equation can be considered as the tangential flow in the crack. Inserting Equation (5) into (4), the lubrication equation can be derived: Equation (6) can be solved by combining an appropriate boundary condition as follows: q(0, t) = Q(t) w(s = l, t) = 0, q(s = l, t) = 0 (7) where l is the fracture length.

Fracture Initiation and Propagation of PPCZ Element
The simulation of hydraulic fracture with CZM is shown in Figure 2. Fracturing fluid is injected into the crack through the wellbore, and the fracture is opened under the fluid pressure. When the fluid pressure at the tip reaches the strength of the formation, the fracture will initiate and start to propagate. where l is the fracture length.

Fracture Initiation and Propagation of PPCZ Element
The simulation of hydraulic fracture with CZM is shown in Figure 2. Fracturing fluid is injected into the crack through the wellbore, and the fracture is opened under the fluid pressure. When the fluid pressure at the tip reaches the strength of the formation, the fracture will initiate and start to propagate. Generally, damage of the PPCZ element mainly conforms to traction separation rules which includes the normal mode and shear mode. In this paper, the mix-mode damage law [36] shown in Figure 3 is used. From Figure 3, the condition of the fracture initiated can be controlled by the quadratic nominal stress law: where Tn is the traction in the normal mode, while Tt is the shear mode traction, N and T are the prescribed maximum traction in the normal mode and shear mode, respectively. The symbol <> is the Macauley operator, avoiding the damage caused by compressive traction, which can be given by: After the initiation, fracture propagation is followed by the Benzeggagh-Kenane criterion that is described by: where GI and GII are the fracture energy in the normal mode and shear mode, respectively. Subscript C indicates the critical value of the fracture energy, η is the material constant.
Equation (10) corresponds to the area of the cyan triangle representing the total fracture energy of the PPCZ element in Figure 3, which can be divided into the yellow part (shear part) and blue part (normal part). If the fracture energy arrives at the critical energy, the material will completely fail and cannot support any traction. In the process of damage evolution, the damage is measured by the total relative displacement of the crack surfaces: Generally, damage of the PPCZ element mainly conforms to traction separation rules which includes the normal mode and shear mode. In this paper, the mix-mode damage law [36] shown in Figure 3 is used. From Figure 3, the condition of the fracture initiated can be controlled by the quadratic nominal stress law: where T n is the traction in the normal mode, while T t is the shear mode traction, N and T are the prescribed maximum traction in the normal mode and shear mode, respectively. The symbol <> is the Macauley operator, avoiding the damage caused by compressive traction, which can be given by: After the initiation, fracture propagation is followed by the Benzeggagh-Kenane criterion that is described by: where G I and G II are the fracture energy in the normal mode and shear mode, respectively. Subscript C indicates the critical value of the fracture energy, η is the material constant. Equation (10) corresponds to the area of the cyan triangle representing the total fracture energy of the PPCZ element in Figure 3, which can be divided into the yellow part (shear part) and blue part (normal part). If the fracture energy arrives at the critical energy, the material will completely fail and cannot support any traction. In the process of damage evolution, the damage is measured by the total relative displacement of the crack surfaces: where ∆ is the total relative displacement, ∆ n and ∆ t are the components in the normal and shear direction. Damage factor d can be determined by: where ∆ max is the maximum of total relative displacement in the fracturing history, ∆ o is the displacement at the initiation, ∆ f is the displacement at the complete damage.
Processes 2020, 8, x 6 of 20 where Δ ̅ is the total relative displacement, Δ ̅ n and Δ ̅ t are the components in the normal and shear direction. Damage factor d can be determined by: where Δ ̅ max is the maximum of total relative displacement in the fracturing history, Δ ̅ o is the displacement at the initiation, f  is the displacement at the complete damage.

Formulation of PPCZ Element
The vital step of simulation is inserting PPCZ elements into edges of other elements. A PPCZ element with friction contact is shown in Figure 4. There are six nodes in one element. The first four black nodes only include the displacement degree of freedom and the rest, describing the fluid flow, merely contain the pore pressure degree of freedom. The xy stands for the global coordinate system and nt for the local coordinate system. Γ is the surface outside the fracture, Γ fc is the internal fracture surface where the friction f is applied.
The displacement and pore pressure of the PPCZ element needs to be converted into a variable at the node. The nodal displacement vector of the 2D PPCZ element in the global coordinate can be written as: T   1  1  2  2  3  3  4  4 , , , , , , , where u and v are the displacement components in the x and y direction. By virtue of the transformation matrix, local nodal displacement u ̅ can be obtained: Figure 3. Mix-mode damage criterion.

Formulation of PPCZ Element
The vital step of simulation is inserting PPCZ elements into edges of other elements. A PPCZ element with friction contact is shown in Figure 4. There are six nodes in one element. The first four black nodes only include the displacement degree of freedom and the rest, describing the fluid flow, merely contain the pore pressure degree of freedom. The xy stands for the global coordinate system and nt for the local coordinate system. Γ is the surface outside the fracture, Γ fc is the internal fracture surface where the friction f is applied.
The displacement and pore pressure of the PPCZ element needs to be converted into a variable at the node. The nodal displacement vector of the 2D PPCZ element in the global coordinate can be written as: where u and v are the displacement components in the x and y direction. By virtue of the transformation matrix, local nodal displacement u can be obtained: where R is the transformation matrix.
In the local coordinate, the displacement separation is given by: Introducing an operator L, the equation above can be expressed as: where the expression of the matrix L is: where R is the transformation matrix.
In the local coordinate, the displacement separation is given by: Introducing an operator L, the equation above can be expressed as:  u Lu (16) where the expression of the matrix L is: When the displacement and pore pressure are assumed to be varying linearly over the cohesive element, linear interpolation functions can be applied to connect the node value with the field value. Therefore, the displacement jump vector at the middle face is formulated as:  Δ Nu (18) where matrix N is the interpolation function that is followed as: where N1 and N2 are the components, N 1 =  (14) and (16) into Equation (18), the relationship between local separation field and global nodal displacement can be derived as: From Equation (20), the width of the crack w is attained: In the same way, pore pressure in the crack can be discretized by interpolation as follows: T ff p  Np (22) where  When the displacement and pore pressure are assumed to be varying linearly over the cohesive element, linear interpolation functions can be applied to connect the node value with the field value. Therefore, the displacement jump vector at the middle face is formulated as: where matrix N is the interpolation function that is followed as: where N 1 and N 2 are the components, . Substituting Equations (14) and (16) into Equation (18), the relationship between local separation field and global nodal displacement can be derived as: From Equation (20), the width of the crack w is attained: Processes 2020, 8,189 8 of 20 In the same way, pore pressure in the crack can be discretized by interpolation as follows: where

Discretization and FEM Implementation
In order to be embedded in the FEM, the governing equation in the strong form should be turned into weak form by the principle of the virtual work. With regard to one element, the total virtual work of the internal force is equal to that of the external force.
where T c is the vector of the cohesive force, T f . is the friction traction, F ext is the external force, δ∆, δw, and δu . are the virtual displacement separation, virtual width, and virtual displacement, respectively. Substituting Equations (20), (21) and (22) into (23), Equation (23) can be rewritten as: Combining Equation (6) with (7), the fluid flow of the equation can be discretized by: where Q is the injection rate. The partial derivative with respect to time is replaced by the first order difference, Equation (25) can be rewritten in the matrix form: where U is the difference of nodal displacement between the current step and previous step. The PPCZ element is programmed with the Fortran code and embedded in the UEL of ABAQUS software [37]. In this subroutine, the element stiffness matrix and right-hand-side vector should be specified, which can be solved by Equations (24) and (26).

Numerical Integral Scheme
According to Camanho's analysis [36], care must be taken in the integral aspect. The Newton-Cotes integral techniques have advantages over Gaussian integral techniques because the Gaussian integral will result in undesired spurious oscillations of the traction field under the large traction gradient condition. Nguyen [30] also pointed out that the Newton-Cotes quadrature must be utilized at least for the KGD problem if the cohesive element was quadratic. Hence, the integral method of PPCZ in this study is the Newton-Cotes integral.
Equations (24) and (26) both involve the curvilinear coordinate s. However, the element local coordinate ξ must be used to implement integral equations in the FEM framework. The relationship between s and ξ is expressed by: where s 1 = 0, s 2 = l e , l e is the length of one PPCZ element. By means of the chain rule, the derivative of N f with respect to s can be described by:

Results and Discussion
In this part, two numerical examples, including the single HF propagation and HF intersection with NF, are presented. These two instances are both simulated with a triangle, quadrangle, and deformed quadrangle mesh. The deformed quadrangle mesh is designed by Cordero [38] for accurate and fracturing needs. The infinite boundary element (CINPE4) is added as well to analyze the effect on the results. All fractures, including HFs and NFs, in these models are constituted by the PPCZ element, and the medium element around the PPCZ element is the 4-node plane strain element (CPE4) or 3-node plane strain element (CPE3). All models are solved by the ABAQUS standard solver that is very efficient for coupling problems in hydraulic fracturing. Finally, a comprehensive suggestion of building up a hydraulic fracturing model will be given.

Simulation with Different Mesh Types
The figure of the KGD problem is shown in Figure 1. The model can be considered in half due to its symmetry. The dimension of the FEM model that is restricted to the simple support is 10 m × 20 m. As shown in Figure 5, the solid red line denotes the hydraulic fracture, the solid red circle is the source of fluid injection. The FEM model is surrounded by the infinite boundary element (IBE) of which the length should be equal to 10 m at least. The mesh near the injection point and the HF is much finer than that in the far-field. The size of the smallest mesh is approximately 0.16 m. All the results of the FEM with different mesh types are compared with the analytical results that are proposed by Detournay [11]. Two categories of analytical results, zero toughness (ZT) solution and large toughness (LT) solution, are both illuminated in Detournay's work and employed here. Parameter K is introduced to distinguish zero toughness and large toughness solutions. The expression of K is given by: Processes 2020, 8,189 10 of 20 where (30) where K Ic is the fracture toughness, E is the Young's modulus, v is the Poisson's ratio. If K > 1, the solution of the KGD problem belongs to the large toughness, otherwise the zero toughness is dominated. The input parameters of the KGD for numerical simulation are shown in Table 1. The value of K is 0.244 and 5.8 for the zero toughness regime and large toughness regime, respectively.  A comparison between numerical results and analytical results with regard to large toughness is shown in Figure 6. DQuad stands for deformed quadrangle mesh, Quad for the normal quadrangle mesh, and Tri for the triangle mesh. The subplot (a)-(b) represents the fracture length, width, and injection point pressure, respectively. Numerical results with different mesh strategy hardly make a difference and agree well with the analytical results. Figure 7 shows the solution of zero toughness, which indicates that numerical consequences with different mesh types are nearly the same in the aspect of crack length and width. For the pore pressure, the performance of the triangle mesh is slightly poorer than the other two mesh at the early injection time. The reason is that the accuracy of the triangle mesh is lower than that of the quadrangle mesh, which has an effect on the PPCZ element. Therefore, the triangle mesh should be avoided dealing with the problem where the effect of pore pressure is significant during the early simulation period.

Simulation with Different Boundary Sizes
The problem of boundary sizes is non-negligible when fractures are involved. The larger model should be established in FEM to approximate the infinite boundary, resulting in the burden of the computer duo to quantities of the mesh. Infinite boundary elements can tackle this question. Parameter RHL is brought in order to measure the ratio of the model height to length. For simplicity, this part only considers quadrangle mesh and large toughness conditions.
Results simulated by different boundary sizes are shown in Figure 8. There is a slight difference in Figure 8a. The tendency of all numerical curves is consistent with the analytical curve except for the curve of RHL = 1.0. This is because that boundary has an influence on the field near the crack tip leading to the poor consequences. This characteristic is more obvious in Figure 8b,c. The discrepancy between the analytical solution and numerical solution with low RHL tends to be large as the time increases. Hence, from Figure 8, the optimal value of the RHL is 3.0 at least if the IBE is not used in the FEM model.

HF Intersection with NF
In terms of complicated fracture networks, the capability of treating fracture intersection for one numerical method is extremely essential. The model of intersection between HF and NF emulated by PPCZ is shown in Figure 9. The intersection angle of the HF and NF is α. HF is injected into injection rate Q to approach the NF. The initial in-situ stress in the horizontal direction is σ h and σ v is in the vertical direction. The infinite boundary is subjected to a simple support. All parameters the example needs are shown in Table 2. In order to better compare with experimental results [39,40], the parameters are selected as reasonable as possible. For convenience, the large toughness category is only taken into account. A series of numerical cases with different mesh types, different stress contrasts, and intersection angle are presented as follows:     Figure 10 shows the fracture profile and the contour of displacement magnitude of intersection between HF and NF under different conditions. The intersection angle of 45 degrees is presented from Figure 10a1 to Figure 10a6 where the first three plots denote the results with the low stress contrast and the rest with the high stress difference. From these figures, we can see that all the HFs propagate along the NFs. The bi-wing of the NF is activated when the stress contrast is low, while only one wing is extended if the stress contrast is high. This is the reason that small stress differences can hardly squeeze the NFs and allow fluid to flow inside. HFs are likely to be opened and prevent the NF from cracking along the down-wing under the large stress differences. According to the figure, it seems that the small stress contrast will lead to tensile mode crack while the large stress contrast results in shear damage. Although there are some differences in the value of displacement, the consequences of the intersection are coincident.
Results of 60 degrees are exhibited from Figure 10b1 to Figure 10b6. If stress differences are relatively low, the phenomenon of HF propagation along the NF will occur. However, HF will penetrate the NF once the stress differences are large. The reason is that the high stress differences restrict the opening of NF and stresses on the other side of NF become much larger leading to crossing the NF. Similar conclusions are arrived at the intersection of 75 and 90 degrees. Concerning the magnitude of displacement, the value of triangle mesh is much larger than the other two mesh types when the intersection angle is small and stress contrast is low. However, quadrangle and deformed quadrangle value turn huge as the intersection angle and stress contrast increases. This is because quadrangle mesh is more sensitive to the variation of circumstances.
In order to guarantee the validity and reasonability of the numerical results, comparison with experimental results is performed in Figure 11. The dashed line is an asymptotic line to distinguish the case of opening and crossing. The data points below the dashed line represent the condition where the HF cannot penetrate the NF while points above the dashed line denote the penetration. From Figure 11, we can find that numerical results are in great agreement with the experiments. Additionally, there are few differences among the consequences simulated by the three mesh types. Hence, they are all effective at least for the HF intersection with NF. restrict the opening of NF and stresses on the other side of NF become much larger leading to crossing the NF. Similar conclusions are arrived at the intersection of 75 and 90 degrees. Concerning the magnitude of displacement, the value of triangle mesh is much larger than the other two mesh types when the intersection angle is small and stress contrast is low. However, quadrangle and deformed quadrangle value turn huge as the intersection angle and stress contrast increases. This is because quadrangle mesh is more sensitive to the variation of circumstances.  In order to guarantee the validity and reasonability of the numerical results, comparison with experimental results is performed in Figure 11. The dashed line is an asymptotic line to distinguish the case of opening and crossing. The data points below the dashed line represent the condition where the HF cannot penetrate the NF while points above the dashed line denote the penetration. From Figure 11, we can find that numerical results are in great agreement with the experiments. Additionally, there are few differences among the consequences simulated by the three mesh types. Hence, they are all effective at least for the HF intersection with NF.

Conclusions
In this paper, the pore pressure cohesive zone (PPCZ) element with contact friction based on the framework of the UEL subroutine is established. Hydraulic fracturing problems can be solved by inserting PPCZ elements into the edges of formation elements. PPCZ elements cannot just simulate the single HF propagation but also the intersection between HF and NF even the complex networks without definition of the crossing criterion. Plenty of numerical examples are emulated for hydraulic fracturing with three types of mesh, including triangle mesh, quadrangle mesh, and deformed quadrangle mesh. In addition, the influence of the infinite boundary is also referred. Finally, the following conclusions can be reached.
(1) Two types of hydraulic fracturing, large toughness and zero toughness, can be correctly captured by the PPCZ element. For the large toughness problem, numerical results solved by three types of mesh make little difference. However, with respect to the zero toughness regime, the pore pressure calculated by the triangle mesh is lower than any other mesh type at the early injection time. Therefore, the triangle mesh should be avoided treating the problem where the significance of pore pressure is terrifically striking during the early simulation period. (2) The infinite boundary element (IBE) should be utilized in solving the hydraulic fracturing problem, which cannot only enhance the computational accuracy but also alleviate the burden of the computer. In case that IBE is not available in some conditions, the ratio of model height to length (RHL) is at least equal to 3.0. (3) HF will propagate along the NF when the stress contrast is small, while HF will cross the NF once ∆σ becomes large. The results of HF intersection with NF simulated by three types of mesh make little difference, which all conform with the experimental results. Hence, the investigation of fracture intersection can be used by any type of mesh mentioned above. Owing to the property of the PPCZ element, only cracking along the edges of the element, triangle mesh, and deformed quadrangle mesh should be applied. In this way, there are more alternatives to cracking direction to increase the validity of the fracture pattern.
Author Contributions: M.C. produced the majority of work on numerical simulation, data reduction, and completion of the manuscript. M.L. mainly designed the methodology of this paper. Y.W. and B.K. verified the numerical results and checked the manuscript. All authors have read and agreed to the published version of the manuscript.