**1. Introduction**

In recent years, the research into the fluid–solid coupling of hydraulic fracturing has attracted the attention of experts and scholars in many fields around the world. In 1925, Terzaghi [1] regarded the consolidation of deformable saturated porous media as a coupling problem of flow boundary–solid boundary deformation. On the basis that the flow field satisfied Darcy's law, taking the total stress and pore fluid pressure of the saturated rock and soil as the state variables, a relatively complete three-dimensional consolidation theory was established. Barenblatt [2,3] (1959) introduced the definition of the field near the crack boundary and analyzed the linear elasticity by using the relevant provisions and assumptions in elasticity.

Therefore, the rational development of some resources with special mining difficulties in the past two decades has also become the global focus of scientific mining researchers. The research and development of safe and effective mining technology for near-vertical coal seams, deep high-thermal coal seams, short-distance thin coal seams, and deep strongimpact-prone coal seams has entered a new stage. In the 1990s, the hydraulic fracturing

**Citation:** Wang, S.; Luo, J. Study on the Shadow Effect of the Stress Field around a Deep-Hole Hydraulic-Fracturing Top-Cutting Borehole and Process Optimization. *Processes* **2023**, *11*, 367. https:// doi.org/10.3390/pr11020367

Academic Editors: Feng Du, Aitao Zhou and Bo Li

Received: 15 December 2022 Revised: 6 January 2023 Accepted: 10 January 2023 Published: 24 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

process carried out by the American Petroleum System for shale gas development and oil penetration and production was introduced into the mining industry and successfully applied to the seam cutting of the hard roof plate of the mine, the pressure relief of the roadway, and the treatment of the weakening of the top coal. In order to determine that the layout plan of the prescribed hydraulic-fracturing (PHF) boreholes can ensure that the fractures can be successfully expanded and the roof can be cut without premature penetration leading to early pressure reduction, it is important and difficult to reasonably design the spacing of the hydraulic-fracturing boreholes.

On this basis, Hutchinson [4], Sneddon [5], Green [6], Warpinski et al. [7], Wang [8], and Nagel et al. [9] have made some progress in determining the spatial distribution characteristics of stress shadowing. Considering the effects of crack spacing, rock mechanics parameters, and in situ stress in three-dimensional conditions, it was found that the influence range of the increase in the minimum horizontal stress was large and extended in the direction of the crack height, but it had little effect on the crack tip region. Zhao et al. [10] studied the three-dimensional propagation mechanism of the hydraulic fracture in a borehole with multiple angles and different spacings in a three-dimensional rock mass and concluded that the increase in the intermediate principal stress suppressed the propagation of the fracture perpendicular to the direction of the intermediate principal stress; the propagation track of the crack in the direction of the intermediate principal stress became straighter. The change in the intermediate principal stress did not change the main direction of the crack propagation (the direction of the maximum principal stress).

In terms of the development of multi-physical field numerical simulation experimental software, in the past two decades, scholars at home and abroad have mainly carried out engineering scale research on the basis of binary numerical models to explore the law of interaction between porous hydraulic fracturing and fracture propagation. Wang [11], Liu et al. [12], and Michael [13] found in experiments that when designing the perforation spacing and aperture, the stress shadow effect could be weakened by reasonably optimizing the setting of large or nonuniform perforation spacing, so as to make the spatial morphology of the hydraulic fracture propagation more uniform.

To summarize, it is of new significance to study the influence of the borehole size and borehole spacing on the shadow effect between the hydraulic fracturing fractures within the three-dimensional rock mass scale model. Building a three-dimensional mechanical model and grid division of three-dimensional multi-cluster perforated rock materials can more realistically simulate the propagation of fracturing fractures in order to then design the reasonable arrangement of hole spacing and aperture, so that the spatial form of fracture propagation is more conducive to engineering applications. It has been planned to obtain the optimal parameters applicable to roof hydraulic fracturing drilling in the Yushen mining area under the influence of special conditions, such as shadow effect, by means of mechanical model establishment, theoretical calculation and analysis, multi physical field coupling numerical simulation experiments and mutual verification.

#### **2. Simulation Engineering Background and Boundary Condition Setting**

*Geological Condition Background and Working Face Overview*

Hydraulic fracturing and roof cutting are widely used to prevent rock burst or strong ground pressure in deep buried or deep mines from 400 m to 1000 m. With the increase in the buried depth and the change in the importance of the tectonic stress in the in situ stress of the unit coal body, the type of in situ stress changes, as shown in Figure 1a. Therefore, the research background of the shadow effect of deep-hole hydraulic-fracturing jet-fracturing boreholes was set as the type III geostress type with vertical stress as the maximum principal stress, as shown in Figure 1e; in the absence of cutting-induced fractures and the influence of the stress shadow effect between adjacent boreholes, the dominant trend of fracture space expansion gradually moves to the minimum horizontal principal stress σ<sup>3</sup> direction deflection. In most cases, deep-hole hydraulic fracturing is conducted from the leading position of the roadway to the hard roof above the stope, see Figure 1b. The overall development trend of fracturing fractures is perpendicular to the advancing direction of the working face and along the vertical stress direction. This fracture development mode has natural advantages for cutting the hard roof along the mining fracture direction, shortening the fracture step, and weakening the leading support pressure of the stope. Therefore, it is particularly important to reasonably optimize the spacing and aperture of the boreholes to reduce the disturbance of the stress shadow effect on the crack propagation and to avoid the influence of water leakage from the adjacent boreholes on the effect of roof cutting.


(**a**)

(**b**) (**c**) (**d**) (**e**) *ı ı ı* **Working face advancing direction** Transportation roadway Railway roadway Goaf x y

> **Figure 1.** Relationship between the spatial development pattern of the deep-hole hydraulic-fracturing fractures and the evolution of the in situ stress types. (**a**) Comprehensive geological histogram. (**b**) The spatial relationship between the layout of the long-hole hydraulic-fracturing drilling site and the stope roof in the advance roadway (*σ*<sup>V</sup> is the vertical maximum principal stress; *σ*<sup>H</sup> is the maximum horizontal principal stress; and *σ*<sup>h</sup> is the minimum horizontal principal stress); the red dotted line is the possible jet direction when the high-pressure water-jet drill pipe is used to induce the kerf. (**c**) Type I geostress. (**d**) Type II geostress. (**e**) Type III geostress.

We selected the geological conditions in the Xiangning mining area of Hedong Coalfield, Shanxi Province, where the buried depth of the coal seam was 555 m~631 m, which was the Shennan'ao coal mine. In order to facilitate a unified calculation, the buried depth parameter was determined to be 600 m, the working face length was 180 m, and the advancing length was 1286 m. The hard sandstone on the roof was treated to weaken the rock pressure and high-energy dynamic load test pieces, and the multilayer mudstone and fine sandstone on the roof were subject to hydraulic fracturing.

As shown in Figure 1c–e, there were three distinct development stages in the spatial form of hydraulic-fracturing fracture propagation under different in situ stress types: 1. the crack initiation stage along the prefabricated fracture or the jet nozzle; 2. the deflection stage dominated by the ground stress; and 3. the continuous expansion along the direction of the maximum principal stress and perpendicular to the direction of the minimum principal stress.

An excessive size of a three-dimensional fluid–structure coupling experimental model leads to non-convergence; however, too small a size makes it difficult to meet the actual rock-block size at the field. Considering this, a similar model of a normal-cube rock mass with a side length of 10 m was established. The boundary conditions were selected to assign parameters to the model through the measured in situ stress state of the basic roof rock in the Xiangning mining area; by comparing the convergence effects of three adaptive mesh generation methods (coarsening, thinning, and refining), it was determined that the refined mesh better considered the above experimental purposes.

The experimental scheme of borehole spacings of 3 m, 3.5 m, and 4 m was intended to further study the common situation of the borehole spacing for hydraulic fracturing in coal mines in "Detailed rules for prevention and control of rock burst in coal mines" published by the Ministry of Coal Industry of the People's Republic of China.

The pore diameters of the high-pressure fracturing drills were taken to be 0.08 m, 0.10 m, and 0.12 m, which were designed according to the three most commonly used boreholes for hydraulic fracturing in coal mines.

Based on the above two considerations, an experimental scheme was designed for the research object of this manuscript: first, the scheme most suitable for the studied engineering background geological conditions was selected by comparing the fluid–structure coupling boundary surface stress and flow velocity of three different hole diameters; then, the influence of the stress shadow effect on the fracturing effect under 3.0 m, 3.5 m, and 4.0 m hole spacing was compared on the basis of this aperture.

Due to the self-adaptive mesh generation function developed in the multi-physical field coupling numerical simulation software COMSOL 5.6, the mesh at the boundary of the borehole wall of the hydraulic fracturing jet could automatically adapt to the transition state of laminar flow and turbulence, better converge the calculation of the surface velocity and the equivalent stress at the fluid–solid coupling boundary in the turbulent state, where the internal flow velocity exceeded 2300 Reynolds number, and densify and set smooth grids at fewer positions at the borehole wall boundary, as shown in Figure 2. In the practical application of underground hydraulic-fracturing engineering in coal mines, the most common drill pipe sizes are 0.08, 0.10, and 0.12 m. These three sizes of drill pipes can simultaneously meet the rapid and high-pressure passage of the large-flow fracturing fluid and have no destructive impact on the coal pillar in the cross roadway. Three drilling sizes with different diameters of 0.08 m, 0.1 m, and 0.12 m were set, with a hole depth of 3 m, 5 m, and 8 m (5% of different hole spacings) × five × 5 m sandstone-rock material test block. The density was *ρ* 2460 kg/m3; the Young's modulus E was 4.2 Gpa; Poisson's ratio was *μ* 0.75; the vertical stress σ<sup>V</sup> was 20 Mpa; the horizontal maximum principal stress σ<sup>H</sup> was 18.5 Mpa; and the horizontal minimum principal stress σ<sup>h</sup> was 15.5 Mpa. First, we determined the optimal aperture; then, we conducted the stress shadow effect experiment under a fixed aperture size.

**Figure 2.** Boundary conditions and mesh generation. (**a**) Boundary load setting method of type I in situ stress rock specimen. (**b**) Grid of 0.08 m borehole. (**c**) Grid of 0.10 m borehole. (**d**) Grid of 0.12 m borehole.

In the fluid solid coupling simulation experiment, for the fluid mechanics calculation part, the water inlet (the starting point of the fracturing fluid inflow), the wall (the part where the fluid contacts the solid) and the water outlet (the end point of the fracturing fluid outflow) are usually set. The water outlet is free in the jet stage, closed in the fracturing stage, and the model is pressurized.

Figure 1b shows the development trend of fracture propagation direction induced by different geostress types, while Figure 2 shows the experimental model to be consistent with the standard specimen of rock mechanics experiment in the laboratory and the drilling hole layout.

The flow field inside the borehole was determined by the inlet flow velocity, outlet pressure, and wall equation, where the inlet flow velocity was 3 × <sup>10</sup><sup>2</sup> m/s, the outlet was open and free, and the orifice pressure was 0 Pa. The flow field wall was smooth without sliding. The total pressure of the internal flow field was 35 Mpa under the pressure holding state.

Considering that the COMSOL multi physical field simulation software is based on the calculation principle developed by finite element method, the excellent simulation on the boundary of the two physical fields and the relatively accurate calculation of the corresponding stress field and velocity field are highlighted when the hydraulic fracturing process is simulated. However, there is also a defect that the fracture expansion cannot be visualized. Therefore, this manuscript illustrates the scientific problems of the research object by comparing its equivalent force field and surface velocity field, and looks forward to fitting the reflection of crack propagation mode in the laboratory with the simulation results in this manuscript in the future research work.

In consideration of both convergence and experimental effect, it is considered that the finer adaptive mesh generation method is most suitable for the experimental study of fluid structure coupling mechanical model established in this manuscript.

Based on the above boundary conditions and the multi-physical field fluid–solid coupling parameter settings, the mesh division and model size are shown in Figure 2.

The turbulent governing equation of the transient solver inside the borehole was as follows:

$$\begin{cases} \rho \frac{\partial \mathbf{u}}{\partial t} + \rho (\mathbf{u} \cdot \nabla) \mathbf{u} = \nabla \cdot [-p\mathbf{I} + \kappa] + F\\ \rho \nabla \cdot \mathbf{u} = 0 \end{cases} \tag{1}$$

$$\kappa = (\mu + \mu\_T)(\nabla \mu\_2 + (\nabla \mu\_2)^T) \tag{2}$$

$$
\rho \frac{\partial k}{\partial t} + \rho (\boldsymbol{\mu} \cdot \nabla) \boldsymbol{K} = \nabla \cdot \left[ (\boldsymbol{\mu} + \frac{\mu\_T}{\sigma\_\mathcal{K}}) \nabla \boldsymbol{K} \right] + P\_\mathcal{k} - \rho \varepsilon \tag{3}
$$

$$\begin{cases} \rho \frac{\partial \varepsilon}{\partial t} + \rho (\mu \cdot \nabla) \varepsilon = \nabla \cdot \left[ (\mu + \frac{\mu\_T}{\sigma\_t}) \nabla \varepsilon \right] + \mathbb{C}\_{\varepsilon 1} \frac{\varepsilon}{k\_2} P\_k - \mathbb{C}\_{\varepsilon 2} \rho \frac{\varepsilon^2}{k\_2} \\ \varepsilon = \varepsilon p \end{cases} \tag{4}$$

$$
\mu\_T = \rho \mathbb{C}\_{\mu} \frac{k^2}{\varepsilon}, \\
 P\_k = \mu\_T [\nabla \mu : (\nabla \mu + (\nabla \mu)^T)] \tag{5}
$$

where *u* is the velocity field component, m/s; *P* is the pressure, MPa; *F* is the boundary pressure of the solid mechanics, MPa; *K* (*K* is an capital English letter) is the turbulent flow energy N·m; *<sup>ε</sup>* is the turbulent dissipation rate; *<sup>ρ</sup>* is the fluid density, kg·m3;  is the symbol of the nabrakin operator; *κ* (kappa '*κ*' is a Greek alphabet) is the directional gradient matrix, *<sup>μ</sup>* is the dynamic viscosity coefficient, n·s/m3; *<sup>C</sup>ε*<sup>l</sup> is the compressibility of the fluid; here, the incompressible flow was selected, and the parameter was 1; and *k* (*k* is an italic lowercase English letter) is the turbulence auxiliary variable constant.

The solid boundary coupling governing equation is:

$$
\rho \frac{\partial^2 u}{\partial t^2} = \nabla \cdot s + \mathbf{F}\_\mathbf{V} \tag{6}
$$

where *s* is the displacement of the solid boundary, m; and FV is the normal stress at the solid boundary, MPa.

#### **3. Internal Flow Field and Boundary Surface Pressure of Solid Mechanics in Different Boreholes**

Based on the above flow-field transient control equations and solid-mechanics boundary equations, the transient solver was used to solve them. The borehole surface velocity and fluid pressure obtained after convergence and equilibrium are shown in Figures 3–5.

**Figure 3.** Fluid–solid coupling boundary state of 0.08 m borehole. (**a**) Surface velocity of inner wall of 0.08 m borehole. (**b**) Surface velocity of inner wall of 0.08 m borehole (velocity peak position). (**c**) Hole wall pressure of 0.08 m hole diameter. (**d**) Hole wall pressure of 0.08 m hole diameter (stress peak position).

**Figure 4.** *Cont*.

**Figure 4.** Fluid–solid coupling boundary state of 0.10 m borehole. (**a**) Surface velocity of inner wall of 0.10 m borehole. (**b**) Surface velocity of inner wall of 0.10 m borehole (velocity peak position). (**c**) Hole wall pressure of 0.10 m hole diameter. (**d**) Hole wall pressure of 0.10 m hole diameter (stress peak position).

**Figure 5.** Fluid–solid coupling boundary state of 0.12 m borehole. (**a**) Surface velocity of inner wall of 0.12 m borehole. (**b**) Surface velocity of inner wall of 0.12 m borehole (velocity peak position 1). (**c**) Surface velocity of inner wall of 0.12 m borehole (velocity peak position 2). (**d**) Hole wall pressure of 0.12 m hole diameter. (**e**) Hole wall pressure of 0.12 m hole diameter (stress peak position 1). (**f**) Hole wall pressure of 0.12 m hole diameter (stress peak position 2).

According to the analysis of the flow velocity–stress state of the inner wall of the 0.08 m borehole as shown in Figure 3, under the abovementioned computational fluid dynamics (CFD) flow-field parameter setting conditions, the peak flow velocity in the 0.08 m borehole reached about 1.47 Mach number, which met the turbulence state, and the peak velocity appeared at the slot position of the borehole wall. The peak pressure on the surface of the inner wall of the hole reached about 9 MPa, and the position of the peak pressure was basically the same as the position of the peak flow rate. This phenomenon indicated that the crack initiation position of the 0.08 m hole was more likely to occur at the wall of the hole but less likely to occur at the fracturing target rock layer at the bottom of the hole. Therefore, it is not recommended to select the 0.08 m hole for drilling and fracturing hard-rock layers under the working conditions assumed in this paper.

As shown in Figure 4, the flow velocity and boundary surface pressure of the internal flow field of the 0.10 m borehole were analyzed. The peak position of the flow velocity mainly occurred in the area between the jet nozzle and the bottom of the borehole. The peak flow velocity reached about 2 × <sup>10</sup><sup>2</sup> m/s. The peak position of the surface pressure also mainly appeared at the hole bottom, and the peak stress was about 10.2 MPa. Compared with the 0.08 m hole diameter, the peak flow velocity and peak pressure were reduced to some extent, but the flow velocity and pressure in the general area were increased, and the area of the peak position was also significantly increased, mainly at the bottom of the hole, which may weaken the roof by causing multiple cracks.

Figure 5 shows the internal flow velocity and fluid–solid coupling boundary pressure state of the 0.12 m borehole. When the 0.12 m borehole was selected for hydraulic fracturing, the peak flow velocity reached 1.4 × <sup>10</sup><sup>2</sup> m/s, which was widely distributed at the slit between the hole bottom, the hole orifice, and the hole wall. The surface peak pressure was about 7 MPa, which was distributed throughout the hole wall. Through comprehensive analysis of Figures 3–5, with the gradual increase in the borehole diameter, the peak flow velocity and peak surface pressure of the flow field weakened to a certain extent but remained on the same order of magnitude as a whole, while the distribution of the peak area increased significantly, which was not advantageous for the pressure of hard but thin rock layers. However, during the fracturing operation of hard roof plates in deep mines, the target fracturing rock layers are generally thicker basic roof layers. Therefore, in these three schemes, the selection of the 0.12 m borehole diameter for the next step of the stress-shadow effect experiment comparison is of more practical significance for the fracturing optimization under the hypothetical conditions in this paper. When the borehole spacing is small, a stress shadow effect appears, and the influence range of the plastic zone overlaps, resulting in the rock between the boreholes being relatively broken, and the stress in the plastic zone decreases. When the two boreholes move further away, the cracks around the boreholes are developed and connected, but the stress shadow effect does not appear; hence, the stress field value between the boreholes increases.
