In the simulation, the reservoir is assumed to be a porous medium. Only a single phase of fluid is saturated into the solid skeleton and pores. The total stress within the reservoir consists of two parts: the total effective stress σ ¯ and the pore pressure p _w . According to the principle of virtual work, the finite element method equilibrium equations can be written as (DASSAULT, 2017): ∫ σ ¯ − p w I δ ε d V = ∫ S t ⋅ δ v d S + ∫ V f ⋅ δ v d V (1)

Where σ ¯ and p _w are the total effective stress and pore pressure, respectively; δε and δv are the virtual strain rate and virtual velocity, respectively; t and f are the surface displacement per unit area and body force per unit volume, respectively; and I is the unit matrix.

During the hydraulic fracturing process, the coal reservoir is always in a saturated state. In a saturated state, the fluid in the porous medium should comply with the continuity equation. The total mass change rate of the fluid in the control body is equal to the mass of the fluid passing through the surface of the control body per unit time (Sun and Zheng, 2020), and the specific equation is as follows: 1 J ∂ ∂ t J ρ w n w + ∂ ∂ χ ⋅ ρ w n w v w = 0 (2) where J represents the rate of volume change of the porous medium; ρ _w , n _w , and v _w stand for the fluid mass density, medium porosity, and the average velocity of the liquid relative to the solid phase, respectively; χ is a space vector.

Figure 2 shows that the flow of fracturing fluid inside the fracture is composed of tangential flow inside the fracture and normal filtration flow on the fracture surface. The fracturing fluid is assumed to be an incompressible Newtonian fluid. The tangential flow equation within the fracture (Carrier and Granet, 2012) is as follows: q = − w 3 12 μ ∇ p f (3) where q represents the volumetric flow velocity passing through the cross-section of the fracture, w stands for the width of the fracture, μ denotes the viscosity of the fracturing fluid, and ∇ p f signifies the fluid pressure gradient in the direction of the fracture.

The definition of normal flow representing the filtration behavior from fractures to porous coal matrix (Guo and Zhao, 2015) is as follows: q t = c t p f − p t q b = c b p f − p b (4) where q _t and q _b indicate the normal flow rates at the top and bottom of the cohesive element, respectively, c _t and c _b symbolize the filtration coefficients of the top and bottom surfaces of the cohesive element, respectively, p _f signifies the fracture flow pressure in the cohesive element, and p _t and p _b are the pore pressures at the top and bottom surfaces of the cohesive element, respectively.

In the construction of the traction separation constitutive model of the cohesive element, we used the quadratic nominal stress criterion (Sun and Zheng, 2020) as the initial damage criterion of the element. When the element stress state satisfies the following functional equation, the element stress state reaches its extreme strength value and initial damage begins to occur: σ n N max 2 + σ s S max 2 + σ t T max 2 = 1 (5) where ⟨σ _s ⟩ is the normal stress of the cohesive element; σ _s and σ _t are the shear stress in two orthogonal directions along the plane in a three-dimensional state; N _max is the ultimate tensile strength of the cohesive element; S _max and T _max are the ultimate shear strength of the cohesive element in two orthogonal directions, respectively.

By introducing the damage variable D to describe the evolution process of element damage, the actual stress state of the cohesive element at a certain instant is defined as follows: σ n = 1 − D σ ¯ n σ ¯ n > 0 σ ¯ n σ ¯ n < 0 σ s = 1 − D σ ¯ s σ t = 1 − D σ ¯ t (6) where D stands for the damage variable, and D=1 means that the cohesive element completely fails (i.e., fracture propagation occurs); σ ¯ n , σ ¯ s and σ ¯ t denote the yield stress states corresponding to the current strain based on the separation displacement in the undamaged state. To describe the evolution of damage variable D with tensile separation displacement, the concept of effective displacement (Camanho and Davila, 2003) needs to be introduced, and its specific definition formula is as follows: δ m = δ n 2 + δ s 2 + δ t 2 (7) where δ _m represents the effective displacement, ⟨δ _n ⟩ signifies the normal displacement; δ _s and δ _t symbolize the tangential displacement in two orthogonal directions.

The constitutive curve of the effective traction force of the cohesive element as a function of effective displacement is shown in Figure 3, where E represents the element stiffness related to geometric parameters; δ m 0 denotes the effective displacement at the initial damage; δ m f signifies the effective displacement at complete separation state.

The mixed fracture energy G ^c varies with the different mixing ratios of type I and type II fracture forms. We adopted a mixed fracture energy model based on the B-K criterion (Sun and Zheng, 2020), and the formula is as follows: G c = G n c + G s c − G n c G s G t η (8) where G n c and G s c are the critical fracture energies for pure mode I fracture and pure mode II fracture, respectively; G ⁿ is the work carried out by the normal traction force; G ^s is the work performed by the tangential traction force; η is a material parameter, which is considered 1.5 in this study.

Given that the reliability of FDEM in simulating single hydraulic fracture propagation has been proven by many scholars (Carrier and Granet, 2012; Sun and Zheng, 2020), the reliability of single-fracture model will not be discussed here. In this study, we focused on verifying the modeling effectiveness of multiple hydraulic fracture propagation.

The propagation direction of multiple hydraulic fractures was compared with Wu’s numerical solutions (Wu and Olson, 2015). According to Wu’s case, a 100 m × 111 m model was established with two horizontal perforations spaced at 10 m intervals. The Young’s modulus of the reservoir is 30 GPa, the Poisson’s ratio is 0.35, the maximum and minimum horizontal stresses are both 46.78 MPa, and the minimum principal stress direction is in the x-direction. The viscosity of the fracturing fluid is 0.001 Pa·s, and the water injection rate is 0.0018 m³/s. According to Wu’s description, we have established the model shown in Figure 4. In addition to the parameters mentioned above, it is also necessary to provide the parameters of the cohesive elements, as shown in Table 1. Finally, a two-dimensional homogeneous reservoir model was established, consisting of 16320 matrix elements and 33000 cohesive elements.

Figure 5 exhibits that the numerical results of fracture propagation simulated in our study are in good agreement with the numerical calculation results of Wu. The expansion of two fractures tends to move away from each other, which may be attributed to the strong stress interference under a low-stress difference. Meanwhile, we found that the two hydraulic fractures in our model are not completely symmetrical. This is because the cohesive elements in the FDEM are inserted between the matrix elements, resulting in a certain dependence of the propagation direction of hydraulic fractures on element meshing.

In rock mechanics, the Weibull distribution function is widely employed to characterize the heterogeneity of rock (Lei and Gao, 2018; Li and Guo, 2020). The density function expression of the Weibull distribution is as follows: ϕ x = m a x a m − 1 e − x a m (9) where ϕ(x) represents the probability density where the strength of the coherent element is equal to x; a stands for the statistical average of material strength; m denotes the heterogeneity coefficient, with a value range of (0, ∞). However, it is worth mentioning that a smaller m represents a stronger degree of heterogeneity in the coal reservoir and vice versa.

Compared to other studies that considered reservoir heterogeneity (Wu and Gao, 2022), we did not heterogenized the mechanical properties of the coal matrix, but instead, we assigned a Weibull distribution function to the strength of the cohesive element. That is, the cohesive elements inserted between matrix in the FDEM macro-scale coal reservoir model represent various natural defects present in the coal reservoir to some extent. Heterogenization of the mechanical properties of the cohesive elements can endow the entire coal reservoir with heterogeneity.

The implementation of the Weibull distribution was achieved through self-programmed Python scripts. To verify the correctness of the program, three samples were taken at m=2, and the sampling results are displayed in Figure 6. From this figure, it can be observed that the sampling results have a high degree of agreement with the Weibull distribution, verifying the correctness of the script.

Here, we established a two-dimensional hydraulic fracture propagation model. The model is exhibited in Figure 7. The size of the model is 20 m × 20 m, which have three clusters of perforations in the middle. The perforation direction is parallel to the maximum horizontal principal stress direction, and the spacing between each cluster of perforations is 0.5 m, with a perforation length of 0.4 m. After dividing the model into element networks using rectangular elements, zero-thickness cohesive elements were embedded between the divided matrix elements. The final model included 49214 matrix elements and 98028 cohesive elements.

The rock mechanics parameters, injection parameters, and in-situ stress conditions used in the simulation are presented in Table 2. The heterogeneity is assigned to the cohesive element according to the method described in Section 3.1. Considering that when m=20, the coal reservoir is approximately completely homogeneous, three different heterogeneous levels of coal reservoirs are set: strong, medium, and weak. Each heterogeneous level of coal reservoir is set with three different heterogeneity coefficients. Therefore, a total of nine coal reservoir models with different degrees of heterogeneity are established as follows: strong heterogeneous coal reservoir (m=1.5, 1.8, 2), medium heterogeneous coal reservoir (m=3, 5, 8), and weak heterogeneous coal reservoir (m=12, 16, 20).

The fracture network morphology of different heterogeneous coal reservoirs after 20 s is depicted in Figure 8. To display the fracture morphology more clearly, all cloud maps are the results of deformation amplification by 30 times. It can be seen that there are significant differences in the morphology of fracture networks under different degrees of heterogeneity. When the heterogeneity level of the coal reservoir was strong (m=1.5, 1.8, 2.0), the fracture network morphology was tortuous but the deflection angle was small, and more secondary fractures were generated. When the heterogeneity level of the coal reservoir was medium (m=3.0, 5.0, 8.0), the continuity of the main fractures became stronger, and the deflection angle of some main fractures increased. When the heterogeneity level of the coal reservoir was weak (m=12.0, 16.0, 20.0), the morphology of the fracture network underwent substantial changes, which resulted in major fractures with larger deflection angles. Some major fractures even extended perpendicular to the direction of maximum principal stress. As the heterogeneity of the coal reservoir weakens, the reason for the gradual increase in the deflection angle of the main fracture may be that when the heterogeneity of the coal reservoir is strong, there are more fracture elements with weak mechanical properties in the coal reservoir, which are controlled by the maximum horizontal principal stress during the fracturing process leading to the formation of secondary fractures with smaller deflection angle. These secondary fractures can easily connect with the main fracture during the fracturing process, thus they control the deflection direction of the main fracture. Therefore, in coal reservoirs with strong heterogeneity, the deflection angle of the main fracture was small. After the heterogeneity of the coal reservoir weakened, the mechanical properties of the fracture elements in the reservoir were almost uniform. Hence, the local stress disturbance caused by the expansion of the main fractures makes it difficult for the secondary fractures to initiate, and the main fracture gradually deviates from the maximum principal stress direction due to the strong stress interference during multiple hydraulic fractures.

Although the final fracture network morphology varied under different heterogeneous coal reservoir conditions, there were still similarities in the evolution process. In Figure 9, m=1.8, 5, and 16 are analyzed as representatives of strong, medium, and weak heterogeneous coal reservoirs. It can be observed that in the early stage of fracturing, a large number of secondary fractures are generated in the middle area of the fractures. This is due to the heterogeneity of the mechanical properties of coal reservoir. Strong stress interference causes the elements around the perforation with weak mechanical properties to initiate, resulting in many isolated secondary fractures around the main fracture. Note that some secondary fractures may not necessarily be in the area where the fracturing fluid flows through, indicating that the formation of fractures is not only due to changes in pore pressure caused by fracturing fluid injection but also possibly due to stress interference between fractures. Some secondary fractures were unable to connect with the main fracture during the fracturing process, leading to discontinuous fracture phenomena (as marked in Figure 9). In coal reservoirs with stronger heterogeneity, this phenomenon was more obvious. This was because the heterogeneity of the coal reservoir was stronger, and the proportion of elements with weaker mechanical properties was also higher, which made it easier to initiate and extend under stress interference. Although all three initial fractures could initiate and form a network of fractures in the early stage of fracturing in different heterogeneous coal reservoirs, as the main fractures on both sides expanded, the middle main fracture could not further expand and form a main fracture under stress interference.

The pressure-time curves at the injection points of different heterogeneous coal reservoirs were extracted from the numerical simulation results (Figure 10). From Figure 10A, it can be seen that in reservoirs with different degrees of heterogeneity, the evolution trend of pore pressure curves is almost similar, which can be divided into the following three stages: 1) The stage of sudden decrease after the pressure rises to its peak: With the continuous injection of fluid, the pressure in the initial fracture continuously accumulates, and under the action of pressure, the width of the fracture continuously widens until the fracture pressure is reached, and the fracture begins to expand. Because the pressure required for fracture propagation is smaller than the pressure required for fracture initiation, a substantial pressure drop is observed. 2) The stage of small amplitude fluctuation of pressure curve: Owing to the presence of coal reservoir heterogeneity, a small amplitude fluctuation is observed in the pressure curve. At this stage, hydraulic fractures constantly turn and connect with each other, and the shape of the fracture network becomes more complex. 3) The stage of stable pressure rise: With the continuous increase of injection time, the continuous injection of fluid leads to a continuous increase in net pressure within the expanded fracture, which further expands the fracture network.

However, in reservoirs with different degrees of heterogeneity, significant changes occurred in fracture pressure. Figure 10B shows the fracture pressure curve under different heterogeneity coefficient. Overall, as m increases (coal reservoir heterogeneity decreases), the fracture pressure gradually increases and finally stabilizes at a certain value. This may be due to the fact that when m is small (with strong reservoir heterogeneity), there are more fracture elements with weak mechanical properties, and these weak elements are more likely to be destroyed and connected during the fracturing process. Hence, it is difficult to maintain pressure at the injection point in reservoirs with strong heterogeneity, resulting in lower fracture pressure; As m increases, the heterogeneity of the coal reservoir weakens, and the weak fracture elements in the reservoir decrease. Thus, the fracture pressure at the injection point exhibits an upward trend; As m further increases, the coal reservoir becomes closer to a completely homogeneous reservoir, and the mechanical properties of the fracture elements in the reservoir are consistent. Therefore, the fracture pressure ultimately stabilizes at a certain value. Moreover, it should be noted that the fracture pressure shows a stepwise upward trend. According to the classification in Section 3.2, the fracture pressure is similar in reservoirs with the same heterogeneity level, and when the heterogeneity level of the reservoir changes, the fracture pressure will increase in a stepped manner.

To investigate the impact of injection rate on multiple hydraulic fracture propagation in heterogeneous coal reservoirs, numerical simulations were conducted with a heterogeneity coefficient of m=5 for injection rate equal to 0.001 m³/s, 0.002 m³/s, 0.003 m³/s, 0.004 m³/s, and 0.005 m³/s. Figure 11 depicts the morphology of the fracture network under different injection rate when the injection time is 15 s. It can be seen that as the injection rate increases, the morphology of the fracture network changes from multiple fractures to single fractures and then to multiple fractures. Hence, it is reasonable to infer that an excessively high injection rate does not necessarily mean a complex fracture network.

Figure 12 exhibits the curve of the fracture pressure at the injection point under different injection rate, indicating a positive correlation between injection rate and fracture pressure. This is because as the injection rate increases, the energy of the fracturing fluid increases, but the filtration loss and resistance loss along the way of the fracturing fluid flow inside the fracture are fixed value. Therefore, as the flow rate increases, it is easier to maintain pressure inside the fracture and the fracture pressure is also larger.

According to previous research results, stress disturbances greater than 1 MPa may have an impact on the propagation of hydraulic fractures (Zhao et al., 2015; Yu et al., 2017). Thus, we defined the area with an increase in σ _h (the minimum horizontal principal stress) greater than 1 MPa as the fracturing disturbance area, and the fracturing disturbance areas of heterogeneous coal reservoirs with different injection rate were obtained accordingly (see Figure 13). Overall, the injection rate had a substantial effect on the fracturing disturbance area. As the injection rate increased, the fracturing disturbance area of the reservoir gradually increased. The fracturing disturbance area of the reservoir had the following characteristics: 1) The fracturing disturbance area was always concentrated at the tip of the main fracture; 2) The fracturing disturbance areas at the tips of different main fractures were prone to converge and thus to form a large area of fracturing disturbance; 3) The fracturing disturbance area generated by the main fracture tip had a double wing shape, and the fracturing disturbance area formed by the convergence of multiple fracture tips had a quasi square shape.

Overall, when conducting multi-stage fracturing in horizontal well in a coal reservoir with a heterogeneity coefficient m=5, a moderate injection rate should be used. Although increasing the injection rate will increase the fracturing disturbance area, it will lead to a decrease in the complexity of the fracture network and an increase in the fracture pressure.

To investigate the effect of perforation spacing on multi-cluster fracturing of horizontal wells in heterogeneous coal reservoirs, numerical simulations were conducted under conditions of perforation spacing of 0.5 m, 1.0 m, 1.5 m, and 2.0 m. Figure 14 illustrates the fracture network morphology under different perforation spacings with an injection time of t=15 s. As the perforation spacing increased, the shape of the fracture network gradually transformed from a complex curved shape to a double-wing curved fracture. This demonstrates that an increase in perforation spacing will greatly decrease the degree of mutual interference between multiple fractures. The double-wing curved fractures obtained in heterogeneous coal reservoirs were not completely symmetrical, and the expansion of one side of the fracture was hindered when the perforation spacing was 1.5 m and 2.0 m. This is due to the discreteness of the mechanical properties of fracture elements in heterogeneous coal reservoirs, which causes the fracturing fluid to be hindered by fractures with strong mechanical properties during the injection process, ultimately resulting in an asymmetric fracture network structure. At the same time, the stress interference on the middle fracture is too large, and it is still difficult for the middle fracture to expand by increasing the perforation spacing.

The perforation spacing also has a certain impact on the fracture pressure at the injection point. Figure 15 shows that as the perforation spacing increases, the fracture pressure at the injection point initially decreases and then stabilizes. When the perforation spacing increased from 0.5 m to 1.5 m, the fracture pressure decreased linearly, demonstrating that the increase in perforation spacing effectively decreased stress interference between fractures. However, as the perforation spacing continued to increase, the fracture pressure did not further decrease. This is because the stress interference of hydraulic fractures has a certain range of influence. When this range is exceeded, the mutual interference between fractures will be very limited. At this time, it will be more difficult to affect the fracture pressure at the injection point by increasing the perforation spacing.

Following the method explained in Section 5.1, we obtained the fracturing disturbance areas under different perforation spacings (see Figure 16). It can be observed that the fracturing disturbance area still has similar characteristics as described in Section 5.1. However, under different perforation spacing conditions, the area and shape of the fracturing disturbance area do not show significant changes.

In a word, when conducting hydraulic fracturing in coal reservoirs with a heterogeneity coefficient m=5, it is recommended to increase the perforation spacing appropriately, as increasing the perforation spacing will reduce the stress interference between the fractures, and at the same time, the fracture pressure will decrease.