The XFEM is applied to model arbitrary discontinuities and the discontinuity in a displacement filed along the crack path without the requirement of re-meshing. The displacement field u is approximated to discrete the equilibrium equation by the special enriched functions with additional node degrees of freedom (Daux et al., 2004; Salimzadeh and Khalili, 2015; Shi et al., 2017; Somnath and Vaibhav, 2024; Marzok and Waisman, 2024). In the XFEM, the displacement u vector can be approximated as Equation 1: u = ∑ I ∈ S a l l N I u x u I + ∑ I ∈ S f r a c N I u x H x a I + ∑ I ∈ S t i p N I u x ∑ l = 4 4 F l x b I l , (1) where Sall is the set of all ordinary nodes; Sfrac and Stip are the set of Heaviside enrichment nodes and the set of fracture-tip enrichment nodes, respectively; u _I is node degrees of freedom in a regular finite element; a _I and b ^l _I represent the enrichment nodal degrees of freedom for the fracture and crack tip, respectively; N _I ^u is the standard finite shape function of node I; H( x ) is the enrichment shape jump functions; and F _l is the singular displacement field around the fracture tip.

The jump function H can be written as Equation 2: H x = 1 , x ≥ 0 − 1 , x < 0 , (2) and x is the signed distance function vector, which determines the crack side. The values H(x)=1 or −1 represent the opposite sides on the crack wall.

The asymptotic crack-tip function F _l (x) is used to calculate the field displacement around the crack tip. This function based on the asymptotic features can be given by Equation 3: F l = r sin θ 2 , r sin θ 2 sin ⁡ θ , r cos θ 2 , r cos θ 2 cos ⁡ θ , (3) where r and θ are the polar coordinates at the crack tip with its origin coordinate system.

The cohesive zone method (CZM) is a general method for modeling the fracture initiation and propagation by using the traction–separation relation (Carrier and Granet, 2012; Chen, 2012; Guo et al., 2015; Nguyen et al., 2017; Dandi et al., 2023). The cohesive zone (shown in Figure 1) can describe the fracture initial loading, the initial damage, and the evolution damage at the failure fracture surface.

Damage initiation refers to the stiffness degradation of the cohesive zone in an enriched element. The process begins when the initial stress or strain conforms to special crack initiation criteria we chose. The maximum principle stress criterion is adopted in this simulation (Equation 4): f = σ max σ max a , (4) where σ ⁰ _max is the maximum allowable principal stress. The symbol⟨⟩represents the Macaulay bracket used to indicate that pure compressive stress cannot initiate the damage. It assumes that damage will be initiated when the stress ratio f reaches the value of 1, which means that the stress magnitude is bigger than the maximum allowable principal stress of the objective.

Once crack initiation occurs, the fracture energy law can be used to evaluate the fracture evolution in the whole process. The Benzeggagh–Kenane (BK law) fracture criterion is introduced in our model as the most acceptable criterion in modeling fracture propagation. The BK law can be described as Equation 5: G n c + G s c − G n c G s + G t G n + G s + G t η = G c , (5) where G represents the energy release rate and the superscript c expresses the critical energy release rate. The subscripts n, s, and t denote normal and two shear directions, and η is a constant for the corresponding material property. In this criterion, the energy release rate of the two shear directions should be equal.

The relationship between damage and stress (normal and shear) in the traction–separation law can be represented as Equations 6–8: t n = 1 − D t ¯ n , t ¯ n ≥ 0 t ¯ n , t ¯ n < 0 , (6) t s = 1 − D t ¯ s , (7) t t = 1 − D t ¯ t , (8) where t _n , t _s , and t _t are the stress components in traction separation behavior without damage. Here, the indexes n, s, and t represent the direction of normal stress and the two shear stresses, respectively. D is the damage variable that represents the average overall damage. In particular, no damage occurs (D=0) at the beginning of the simulation, and the cohesive element also satisfies the condition without damage under pure compressive stress. In addition, the complete degradation of the cohesive element at all integration points occurs when D=1, which means zero load-carrying capacity for these elements.

The evolution of the damage variable, D, is expressed in Equation 9: D = δ n f δ n max − δ n 0 δ n max δ n f − δ n 0 , (9) where δ n f and δ n 0 are the displacement at the complete failure and the initial opening before the damage occurs, respectively. δ n max is the maximum displacement during the loading history.

The flow patterns in the pore cohesive model are shown in Figure 2, in which the tangential flow and normal flow are both taken into consideration. The continuity flow is assumed to be incompressible Newtown fluid, with the normal flow permeating into the porous medium, as well as the tangential flow across the gap opening. The pore pressure on the fracture surface can be simulated by introducing the pore pressure node in the enrichment elements. This allows the fully coupled seepage–stress hydraulic fracturing to be modeled.

The tangential flow formulated from the Poiseuille law can be described by Equation 10: q d = w 3 12 μ ∇ p f , (10) where q is the fluid volume through the fracture, w is the fracture width, μ is the viscosity of the fracturing fluid, ∇ p f is the gradient of pressure along the fracture, and d is the gap opening.

The normal flow representing the fluid leak-off can be described as Equation 11: q t = c t p i − p t q b = c b p i − p b , (11) where q t and q b are the fluid rate permeating into the up and side surfaces, respectively; c t and c b are the leak-off coefficients for the top and bottom cohesive layers, respectively; p t and p b are the pore pressure on the top and bottom surfaces, respectively; and p i is the pore pressure in the middle of the fracture.

As a two-dimensional fluid structure coupling model, all vertical displacements of the grid are constrained. Meanwhile, node displacement at the edge of the model is constrained to prevent the overall deviation of the model. Other related boundary conditions are combined with practice and determined in the following section.

To investigate the effect of the perforation angle on the near-wellbore fracture propagation, the simulation models of different angles θ=15°, 30°, 45°, 60°, 75°, and 90° are presented in this section. The sensitive parameters remain constant, i.e., elasticity modulus E=10 GPa, Poisson’s ratio υ=0.2, horizontal stress difference ▽σ=2 MPa, tensile strength σ _t =1.2 Mpa, and fluid injection rate Q=0.0005 m³/s. Meanwhile, the other parameters for reservoirs and fracturing treatment are identical to those given in Table 2. Figure 5 shows the complicated trajectory of hydraulic fracture with different perforation angles. It can be seen that sharper fracturing reorientation happens when the perforation angle is low, especially when the value is less than 60°. In this situation, the fracture initiates from the perforation and sharply reorients in the direction of maximum horizontal stress with relatively low propagation pressure in the fracture, while a long-curving distance fracture can be captured when the perforation angle exceeds 60°. The largest fracture propagation pressure occurs at θ=90°, which attains a high value of 33.53 MPa, as shown in Figure 5F. This is because the severe reorientation may lead to an excessive curve of the fracture morphology, which needs more fracturing fluid and higher injection pressure to maintain the propagation of fracture. Figure 6A shows the relation between the perforation angle and reorientation angle of horizontal distance dx and vertical distance dy. A larger perforation angle can result in more reorientation of the hydraulic fracture, which represents a longer reorientation distance in the vertical and horizontal directions. The largest value of reorientation dx=12.21 m and dy=7.24 m can be achieved at the perforation angle θ=90°. However, the sharpest change in the reorientation distance occurs when the angle increases from θ=60° to θ=75°, corresponding to the distance change from dx=2.83 m to dx=11.83 m and dy=1.21 m to dy=6.64 m. In order to obtain easier fracture propagation and proppant placement, a lower curving distance of fracture is needed in the filed application. Therefore, a lower perforation angle adjusted to be less than 60° should be adopted to obtain a plane fracture according to our simulation model.

The elasticity modulus is an important property affecting the fracture propagation in hydraulic fracturing treatment. In all four cases, the perforation angle has the same value of θ = 60° and the fluid injection rate remains constant at 5 × 10^-4 m³/s. Then, tensile strength, Poisson’s ratio, and horizontal stress difference are fixed at 1.2 MPa, 0.2, MPa and 2 MPa, respectively. The other parameters are given in Table 2. As expected, the fracture initiates from the perforation and rotates to the direction of maximum horizontal stress with different values of elasticity modulus. However, a larger elasticity modulus results in a larger curving distance and more reorientation before the fracture rotates to the direction of maximum horizontal stress. Meanwhile, a larger elasticity modulus can also lead to a higher fracture propagation pressure, which can make the fracture grow longer and narrower in the same injection condition. The maximum propagation pressure and fracture length among all the four cases reach up to 36.88 MPa and 13.4 m, respectively. Nevertheless, the fracture width decreases to the lowest value of 6 mm in all four scenarios. Figure 6B shows the transformation law of reorientation distance with different elasticity moduli. The larger elasticity modulus exhibits longer reorientation distance in both the vertical and horizontal directions, but the distance scope changes less when the elasticity modulus increases from 10 GPa to 40 GPa. For example, the vertical reorientation distance only changes from 1.21 m to 1.52 m, and the horizontal reorientation distance changes from 2.83 m to 4.36 m. Therefore, the elasticity modulus can influence the morphology of fracture with a higher elasticity modulus, resulting in a long curving distance and more reorientation for hydraulic fracture. In addition, a larger elasticity modulus contributes to generating a longer and narrower fracture in the same injection condition. Partly in the same injection condition, the hydraulic fracture can penetrate deep into the reservoir with a large elasticity modulus when the fracture finally returns to the direction perpendicular to the minimum horizontal stress.

The effect of different Poisson’s ratios υ = 0.20, 0.22, 0.24, 0.26, 0.28, and 0.30 on the propagation trajectories of hydraulic fracture is studied in this section. The fluid injection has the same value of 0.0005 × 10^-4 m³/s, and the elasticity modulus remains constant at 10 GPa. Furthermore, tensile strength, horizontal stress difference, and perforation angle are 1.2 MPa, 2 MPa, and 60°, respectively. The other simulation parameters are listed in Table 2. The fracture tends to initiate in the direction of the original perforation and finally turns to the direction of maximum horizontal stress. A larger value of Poisson’s ratio can lead to a larger propagation pressure in the fracture, increasing from 22.91 MPa to 24.49 MPa in the all six models. Second, the change in fracture morphology can be hardly observed in the different schemes of different Poisson’s ratios. It is obvious that the change in Poisson’s ratio can influence the propagation condition and fracture morphology with limited scope. Figure 6C shows the relation between the reorientation distance and Poisson’s ratio. Both vertical and horizontal reorientation distances increase as Poisson’s ratio increases; yet the values of the distance change less. In particular, the vertical reorientation distance only increases from 1.42 m to 1.45 m, increasing by only 0.03 m. It can be seen that the reorientation and trajectory of hydraulic fracture are generally slightly influenced by the change in Poisson’s ratio.

Tensile strength should be a key crack parameter of rock mass in the hydraulic fracturing treatment. To analyze the reorientation mechanism of fracture, the effect of different tensile strengths σ _t = 1.2 MPa, 2.2 MPa, 3.2 MPa, 4.2 MPa, 5.2 MPa, and 6.2 MPa is discussed in this section. The other five sensitive parameters should be kept constant, i.e., elasticity modulus E = 10 GPa, perforation angle θ=60°, Poisson’s ration υ=0.2, horizontal stress difference ▽σ = 2 MPa, and injection rate Q = 5 × 10^-4 m³/s. Other parameters, including fracturing treatment and rock mechanics, are consistent with those given in Table 2. The fracture trajectory and behavior of hydraulic reorientation can be significantly influenced by transforming the tensile strength. A long curving distance and complicate trajectory can be captured in the larger-tensile strength scheme. It is worth noting that more reorientation occurs in the propagation of hydraulic fracture with the increase in tensile strength. Simultaneously, a great increase in propagation pressure within the fracture can be observed when the tensile strength increases, for instance, the propagation pressure can reach up to approximately 33.83 MPa when σ _t = 6.2 MPa. Furthermore, the geometry of fracture, including the fracture length and width, also changes with the increase in tensile strength. The maximum fracture opening attains an extremely large value of 1.77 cm when tensile strength is 6.2 MPa, while a low fracture width of 1.07 cm is obtained when tensile strength is 1.2 MPa. The vertical and horizontal reorientation distances are influenced by tensile strength, as shown in Figure 6D. Both vertical and horizontal reorientation distances increase with the increase in tensile strength; in particular, the horizontal reorientation distance increases more than the vertical reorientation distance. Therefore, a long curving distance is associated with a larger tensile strength, which can cause more serious reorientation and complicated trajectory of fracture propagation.

Horizontal stress difference controls the main mechanism of fracture reorientation, as proven by several researchers (Zou et al., 2018). The perforation angle and injection rate are fixed at θ = 60° and Q = 5 × 10^-4 m³/s, respectively. Reservoir parameters including elasticity modulus, Poisson’s ratio, and tensile strength have the same values of E = 10 GPa, υ = 0.2, and σ _t = 1.2 MPa, respectively. Other parameters are consistent with those given in Table 2. The fracture trajectory with increasing horizontal stress is significantly different in all six simulation cases. A lower horizontal stress difference can result in a more complicated trajectory. Lower propagation pressure within the fracture causes the increase in horizontal stress difference, which will result in quick reorientation of the initiated fracture. In particular, the trajectory of the hydraulic fracture shows a sharp reorientation after initiating from the perforation when the horizontal stress difference exceeds 8 MPa. Figure 6E shows the revolution of the reorientation distance with the increase in horizontal stress difference. The horizontal reorientation distance shows a more obvious decrease than the vertical reorientation distance. The vertical reorientation stress changes from 0.05 m to 1.21 m, which can hardly be ignored. Therefore, without sufficient knowledge about the direction of in situ stress, a larger stress difference may cause a timely reorientation of the near-wellbore fracture, for which less reorientation and simple fracture trajectory will be obtained to benefit for the fracture propagation and proppant placement.

The fluid injection rate is an important factor that influences the propagation behavior of the hydraulic fracture. In this case, different trajectories of hydraulic fracture for different injection rates (Q = 3 × 10^-4 m³/s, 5 × 10^-4 m³/s, 7 × 10^-4 m³/s, 9 × 10^-4 m³/s, 1.1 × 10^-3 m³/s, and 1.3 × 10^-3 m³/s) with a perforation angle θ=60° are established. The reservoir parameters have the same values in all six cases, i.e., elasticity modulus E = 10 GPa, Poisson’s ratio υ = 0.2, tensile strength σ_t = 1.2 MPa, and horizontal stress difference ▽σ = 2 MPa. Other parameters are consistent with those presented in Table 2. It can be concluded that a larger propagation pressure within the fracture should be provided for extension with the increase in the injection rate. Furthermore, a higher injection rate creates a wider fracture along the direction of the perforation angle after the fracture initiates. It should also be noted that a more complicated fracture trajectory that takes up more fracturing time in reorienting to the direction of maximum horizontal stress is formed because of the large injection rate. The reorientation distance increases with the increase in the injection rate, as shown in Figure 6F. The larger horizontal reorientation distance dx can be observed when the injection rate increases to Q = 1.3 × 10^-3 m^-3/s compared with the vertical reorientation distance dy. Second, the vertical reorientation distance dy changes in a relatively less range, just increasing from 0.94 m to 6.34 m with the increase in the injection rate in all six cases. Therefore, a larger curving distance can be achieved at a higher injection rate, which needs more fluid volume to extend a longer reorientation fracture. This is because the larger reorientation fracture along the perforation angle is more likely to be generated because of the large propagation pressure within the fracture, which is provided by the high injection rate. For a longer and preferred fracture plane, a lower injection rate is recommended for easier reorientation at the beginning of the fracturing; however, a larger injection rate should be adopted to generate a longer and wider fracture after the initiated fracture finally rotates to the direction of maximum horizontal stress.

The above research conclusions have a reference value for perforation schemes in practical engineering applications. The degree of distortion of fractures under different geological conditions can be foreseen to a certain extent, which is beneficial for the utilization of remaining oil.