In this study, the physical fields involved in the fracturing process of reservoir rock include the temperature field, fluid field, and solid field (Wang et al., 2019; Wang, 2020). The matrix deformation governing equation of reservoir rock is as follow: L T ( σ ′ − α m p l ) + ρ B g = 0 , (1) where L is the differential operator, σ ′ is the effective stress tensor, α is the Biot coefficient, m is the element tensor, p l is rock mass pore fluid pressure, p B is saturated bulk density of rock, and g is the gravity vector.

The following Eqs. 2, 3 are the governing equations of rock matrix seepage and fluid flow in fractures, respectively: d i v [ k μ l ( ∇ p l − ρ l g ) ] = ( φ K l + α - φ K s ) d p l d t + α d ε v d t , (2) ∂ ∂ x [ k f r μ n ( ∇ p n − ρ f n g ) ] = S f r d p n d t + α ( Δ e ˙ ε ) , (3) where k is the inherent permeability of pore structure, μ l is the fluid velocity in the pore, K l is pore fluid stiffness, K s is solid skeleton stiffness, ε v is the volumetric strain, t is the current moment, k f r is fracture inherent permeability, μ n is fluid velocity in fracture, p n is fluid pressure in fracture, ρ f n is fluid density in fracture, S f r is the parameter describing rock compressibility under fluid action, and Δ e ˙ ε is the fracture strain rate. Some detailed contents, such as numerical implementation, proppant model, elastic constitutive equation, and fracture criterion based on fracture energy, were omitted to avoid redundancy, which could be found in the related references (Wang et al., 2019; Wang, 2020).

The governing equation of thermal effects between the rock matrix and fluid in pores and fractures is as follows: d i v [ k b ∇ T f ] = ρ b c b ∂ T f ∂ t + ρ f c f q f ∇ T f , (4) where k b is the thermal conductivity coefficient, T f is the fluid temperature, ρ b is the volume density, c b is the specific heat coefficient, ρ f is the fluid density, c f is the specific heat coefficient of fluid, and q f is the Darcy fluid flux.

The differential governing Eqs. 1–3 of solid deformation of the rock matrix, fluid flow in pores and fractures, and thermal effects are discretized using the conventional finite element method. The form of heat transfer between element nodes is shown in Figure 1, and the temperature and heat flux at the nodes are as follows: q c 1 = α c ( T N ) ( T N − T f 1 ) , (5a) q c 2 = α c ( T N ) ( T N − T f 2 ) , (5b) where q c 1 and q c 2 are the heat flux transmitted at the fracture plane node, T f 1 and T f 2 are the temperature value at the fracture plane node, T N is the temperature value of the node within the fracture, and α c is the contact thermal conductivity. Temperature changes in rock cause volume expansion and contraction: Δ V V = α T Δ T , (6) where Δ T is the temperature change of the rock element, Δ V is the volume change, V is the initial volume, and α T is the linear thermal expansion coefficient of the rock matrix.

The engineering-scale 3D model of multistage hydrofracturing of horizontal well with multiple perforation clusters in deep tight rock is established, as shown in Figure 2. There are five perforation cluster locations, sequentially numbered from 1 to 5. The basic physical parameters and perforation cluster spacing settings of the model are shown in Table 1. According to the different initiation sequences of perforation clusters, the fracturing scheme can be divided into the sequential (fractures are fractured in the order of 1→2→3→4→5), simultaneous (fractures are fractured simultaneously), and alternate fracturing (fractures are fractured in the order of 1→3→2→5→4 (Wang, 2020), and the model is analyzed by finite element-discrete element methods in program package ELFEN (Rockfield Software Ltd., 2016). In this study, the three-dimensional deflection and induced shear stress evolution of alternate fracturing with different perforation cluster spacing are analyzed.

Figure 3 shows the thermal gradient on the 3D fracture surfaces and surrounding rock. The value of thermal gradient can represent temperature difference and thermal diffusion trend. It can be seen that there is a significant temperature gradient at the fracture surface due to the temperature difference, which makes the temperature domain rapidly diffuse. The temperature gradient of the rock mass from the fracture area to the periphery gradually decreases until the far-field temperature gradient disappears.

Figure 4 shows the final fracture morphology and stress field results under different perforation cluster spacing in alternate fracturing. When the cluster spacing is large, the fractures are nearly parallel and stable. With the decrease of the perforation cluster spacing, the mutual interference between fractures gradually increases, and the deflection intensifies. The spacing of perforation clusters is an important factor affecting the spatial deflection behavior of spatial fractures.

Figure 5 shows the results of fracture propagation at each stage of alternate fracturing. The first fracture propagates in space close to the plane because there is no interference from other fractures. The fracture 3 propagates in space close to the plane due to the spacing between fracture 1 and fracture 3 is 100 m (This is equivalent to increasing the spacing between sequentially activated fractures). The fracture 2 is disturbed by fractures 1 and 3, and the deflection is intensified. Similarly, fracture 4 propagates after fracture 5 has completed fracturing and is deflected by disturbance from fracture 3 and fracture 5.

A comparison of induced shear stress disturbance under thermal-hydraulic (HM) coupling and THM coupling is shown in Figure 6. The upper left side of the fracture is positive shear stress and the right side is negative shear stress. At the same time, the shear stress is negative at lower left of the fracture, andshear stress is positive at the lower right of the fracture. By comparing Figures 6A,B, it is found that the induced shear stress field area around the fracture is larger when the thermal effect is considered, and the shear stress disturbance caused by the reduction of induced shear stress superposition is stronger. The degree of fracture propagation and spatial deflection also increase. Therefore, thermal diffusion has a great influence on the induced shear stress field disturbance around the fracture, which needs to be considered.

Figure 7 shows the shear stress disturbance of alternate fracturing at different spacing of perforation clusters. At 100 m of perforation cluster spacing, the induced shear stress field around the fractures are slightly superimposed and reduced, and the fracture propagation is less affected by stress shadow effect. With the decrease of the spacing between perforating clusters, the superposition and disturbance of stress field between fractures gradually intensify. At 50 m of perforation cluster spacing, the positive and negative shear stress areas near the fractures superimpose obviously, resulting in the shear stress field on the right side of fracture 4 reduces obviously, and the fracture deflected to the left side of the larger shear stress field. At 25 m of perforation cluster spacing, the induced shear stress field around the fractures are superimposed and reduced seriously, resulting in the shear stress field on the right side of fracture 4 reduces obviously, and the deflection degree to the larger shear stress area on the left side becomes larger.

Figure 8 shows the shear stress disturbance of alternate fracturing. In the first stage, the stress around the first fracture is almost symmetrical and the fracture propagates in a plane direction. In the second stage, fracture 3 begins to propagate, and since the two fractures are separated by uninitiated fracture 2, the stress shadow area around fracture 1 does not cover the shear stress area caused by fracture 3, so fracture 3 is not affected by stress shadow and propagates in plane. In the third stage, the shear stress field around fracture 2 is superimposed and reduced with the stress shadow around fracture 1 and fracture 3, causing fracture 2 to propagate and deflect towards the larger stress area. Similarly, in the fifth stage, the induced shear stress field around fracture 4 is superimposed and reduced with the stress shadow around fracture 3 and fracture 5, and the fracture deflection is intensified.