Numerical Simulations of Radial Well Assisted Deflagration Fracturing Based on the Smoothed Particle Hydrodynamics Method

Abstract

The technology of radial-well-assisted hydraulic fracturing is applied in the stimulation of low-permeability hydrocarbon reservoirs where commercial production cannot be achieved by the conventional fracturing method. Here, a study on the reservoir stimulation effect and the fracture propagation pattern of radial-well-assisted deflagration fracturing was carried out. Based on smooth particle hydrodynamics (SPH), rock mechanics theory, and finite element theory, a numerical model of radial-well-assisted deflagration fracturing was established by integrating the JWL state equation. Research on the effects of the deflagration position, radial well azimuth and horizontal principal stress difference on the fracture propagation was carried out. The results show that the deflagration position, radial well azimuth and horizontal principal stress difference have significant effects on the fracture area in deflagration fracturing. The closer distance from the deflagration position is, the larger the radial well azimuth and the smaller the horizontal stress difference are, leading to a larger fracture area, which is conducive to reservoir stimulation. During fracturing, both shear fractures and tensile fractures are formed. The formation and conversion of shear fractures and tensile fractures are related to the deflagration position, radial well azimuth, horizontal principal stress difference, etc.

1. Introduction

Ultrashort radius radial horizontal wells (which drilling professionals also call radial wells for short), a type of horizontal well, have superiority in the development of low-permeability or fractured reservoirs, and they are also applied in the stimulation of old oil wells. Nevertheless, the increment achieved in oil and gas production by drilling the radial well to connect it to the reservoir is always not satisfactory. Thus, the technology of radial-well-assisted fracturing is feasible in the stimulation of these reservoirs. Efforts have been devoted to the study of this technology, and some progress has been made [1,2,3,4]. In radial-well-assisted hydraulic fracturing, stress concentration occurs mostly around the radial well root, where the fracture is first initiated, and the reservoir stimulation degree is similar to that of vertical well fracturing [5]. Based on this method, a new stimulation technology based on the use of explosives to operate deflagration fracturing in radial wells is proposed here, with the aim of reservoir volume stimulation. Compared with the perforated hole, the radial well has a larger length and provides a path for the transportation of high-energy fracturing materials to the far end. Deflagration fracturing leads to significant improvements in oil and gas production and the generation of network fractures around the radial well used to communicate with the reservoir. The drainage area and reservoir connectivity are enhanced, and the oil and gas production is increased.

The deflagration fracturing technology has been applied widely. A.V. Dubovik carried out a study on the mechanical impact properties of explosive materials by the shell fracture method and mathematical characterization of the ignition model through theoretical analysis [6]. In 2016, V. Vershinin established the temperature–stress coupling equation of the explosion of liquid explosives in the near-wellbore region, analyzed the effects of the explosion parameters on post-stimulation production, and discussed the generation mechanism of the micro-fracture network near the wellbore [7]. Progress has been made in deflagration fracturing technology. V. n. Odintsev carried out simulations of the explosive effect on the gas dynamic state in an outburst-hazardous coal bed, and the mechanism of the gas dynamic fracture of the outburst-hazardous coal bed during the explosion in the coal bed was considered [8]. The conditions of fractured reservoirs in the layer furthest from the deflagration hole were studied. A.K. Raina carried out a study on the complex interactions between the explosion and rock fracture [9]. In China, a great deal of effort has been put into the study of deflagration fracturing technology. For example, some progress has been made in the study of deflagration fracturing technology by the China National Petroleum Corporation (CNPC) High Energy Gas Fracturing Technology Center, Xi’an Shiyou University, which has been studying this technology since the end of 20th century. In recent years, the study of deflagration fracturing technology has focused on the simulation of multi-stage pulsed combustion–fracturing coupling and the optimization of high-energy fracturing parameters [10,11], fracture dynamic propagation in high-energy gas fracturing in casing perforated wells [12], and deflagration transition models [13].

However, studies on the mechanisms of the deflagration and fracture propagation in radial-well-assisted deflagration fracturing have not been carried out because of the different borehole lengths, borehole sizes, and perforation modes in radial wells. The application of radial-well-assisted deflagration fracturing techniques in the oil and gas fields is limited. In this study, a new radial-well-assisted reservoir stimulation technique was proposed. The technology feasibility was analyzed through numerical simulation. Moreover, research on the mechanism of radial-well-assisted deflagration fracturing provides an understanding of the rock failure mode, fracture morphology, and fracture parameters in the radial well and a theoretical basis for subsequent experimental studies and field practice.

2. Deflagration Fracturing Theory

2.1. Principle of the SPH Method

In 1977, Lucy and Gingold et al. proposed the SPH method. SPH is a pure Lagrange theory and is solved numerically without meshing, which avoids the problem of the interface between Euler mesh and material in the Euler description. The SPH method enables the simple computation of the interface between media and avoids the sharp drop in the computation precision caused by mesh distortion in an extremely deformed structure. Thus, the SPH method is applicable for solving the large deformation of dynamic structures in cases of high-speed collision. Many achievements have been made in research on the SPH method. The research on the explosion theory based on the SPH method targets the simulation of soil breakage or large deformation caused by stress waves using the SPH method or improved SPH method [14,15,16,17]; the simulation of rock failure and fracture propagation under stress [18,19,20]; and the simulation of underwater deflagration processes [21,22]. The related theories of rock stress failure provided a basis for the study described here.

The SPH method is conducted in two steps. The first is kernel approximation, which is the core step, and here, the function is expressed as the integral form. The second is particle approximation, which is discretization. The integral operation is discretized into summation computation, and the volume domain is discretized into finite particles with the properties of the mass, velocity, pressure, etc. [23].

2.1.1. Kernel Approximation

The first step of the SPH method is kernel approximation, where a function is expressed as an integral form with the property of the Dirac smooth kernel function. The functions f(x) and x are set as the spatial vectors, and an integral form is expressed as:

$$f\left(x\right)={\displaystyle {\int }_{\mathsf{\Omega }}f\left(x^{\prime }\right)\delta \left(x-x^{\prime }\right)d{x}^{\prime }}$$

(1)

As the δ function is continuous but not differentiable, and it is not computable, the kernel function W is used to substitute the δ function. A function of a particle in a domain Ω can be approximated as the Dirac smooth kernel function W as follows:

$$f\left(x\right)={\displaystyle {\int }_{\mathsf{\Omega }}f\left(x^{\prime }\right)W\left(x-x^{\prime },h\right)d{x}^{\prime }}$$

(2)

where f(x) is a function at the 3D coordinate x, $x-x^{\prime }$ is the particle spacing, and $h$ is the smooth length of the particles.

Using the Gaussian divergence theorem, the equation is transformed into the spatial derivative of the kernel function as follows:

$$\nabla f\left(x_i\right)={\displaystyle {\int }_{\mathsf{\Omega }}f\left(x^{\prime }\right){\nabla }_{i}W\left(x-x^{\prime },h\right)d{x}^{\prime }}$$

(3)

2.1.2. Particle Approximation

The second step of the SPH method is particle approximation, where the integral operation in the first step is discretized into the summation operation. The domain is discretized into a finite number of particles with the properties of the mass, velocity, pressure, density, etc. The approximate derivation is expressed as follows:

$$f\left(x\right)={\displaystyle {\int }_{\mathsf{\Omega }}f\left(x^{\prime }\right)W\left(x-x^{\prime },h\right)d{x}^{\prime }}={\displaystyle \sum _{j=1}^{N}f\left(x_j\right)W\left(x-x_j,h\right)\Delta V_j}={\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}f\left(x_j\right)W\left(x-x_j,h\right)}$$

(4)

where j is a SPH particle in the domain i, m_j is the mass of the SPH particle, ρ_j is the density of the SPH particle, and N is the number of particles in the smooth length.

Thus, the function approximation of the final particle is expressed as:

$$\nabla f\left(x_i\right)={\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}f\left(x_j\right){\nabla }_i{W}_{ij}}$$

(5)

2.2. JWL State Equation

Deflagration is a process of violent chemical reaction, where a large amount of energy is released and accompanied by a high temperature, high pressure, and high-speed stress wave. In deflagration fracturing, a deflagration point is required to ignite the explosive or gunpowder. Here, TNT with the standard parameters is used as the explosive, and the explosion reaction process is simulated using the JWL state equation. The JWL state equation [24] has been widely used in engineering applications, and accurate physical property parameters have been reported [25]. These parameters are obtained by a large number of cylinder experiments and fitting. The JWL state equation is expressed as follows [26]:

$$P=A\left(1-\frac{\omega }{R_{1}V}\right)e^{-R_{1}V}+B\left(1-\frac{\omega }{R_{2}V}\right)e^{-R_{2}V}+\frac{\omega e}{V}$$

(6)

where P is the pressure of the deflagration product, MPa, V is the relative specific volume of the deflagration product, and e is the specific internal energy of the deflagration product. The unit of R₁ and R₂ is GPa, A, B, and ω are dimensionless, and they are all constants.

2.3. Criterion for Fracture Propagation in Cohesive Unit

The principle of deflagration fracturing is that the high pressure is generated within the wellbore through deflagration, and rock failure occurs when the pressure exerted on the rock by the high-energy gas is higher than the maximum bearing pressure of the rock on the borehole wall. According to the damage mechanics principle and the stiffness decay method, the cohesive damage unit describes the process of fracture initiation and propagation in rock under stress. Before damage in the cohesive unit, the stress–strain relation is expressed as follows [27]:

$$t=\left\{\begin{array}{l}{t}_n\\ t_s\\ t_t\end{array}\right\}=K\delta =\left[\begin{array}{ccc}{K}_{nn}& 0& 0\\ 0& K_{ss}& 0\\ 0& 0& K_{tt}\end{array}\right]\left\{\begin{array}{l}{\delta }_n\\ {\delta }_s\\ {\delta }_t\end{array}\right\}$$

(7)

where t is the nominal stress, MPa; t_n and t_s are the normal stress, MPa; t_t is the shear stress, MPa; δ is the nominal strain; δ_n and δ_s are the normal strain; δ_t is the shear strain; K is the stiffness matrix; K_nn is the directional stiffness; and K_ss and K_tt are the shear stiffness.

Rock failure occurs when the stress is higher than the limit of the rock strength. Rock failure leads to reductions in the mechanical properties, and the stress is negatively correlated with the displacement. As the rock failure continues, the relative displacement of the fracture surface increases, which leads to reductions in the rock stiffness and the stress on the rock. When the stress drops to zero, the rock is completely damaged, and the fracture is initiated.

Different formation conditions lead to differences in the rock physical properties. The fractures are generally composed of tensile fractures caused by tensile stress and shear fractures caused by shear stress. Here, fracture initiation is expressed by the maximum tensile stress criterion as follows:

$${\sigma }_{\max }\ge P_t$$

(8)

where σ_max is the maximum principal stress on the rock, MPa, and P_t is the rock tensile strength, MPa.

Fracture initiation in the cohesive unit is characterized by the BK criterion [28] as follows:

$$G_n^c+\left(G_s^c-G_n^c\right){\left(\frac{G_s+G_t}{G_n+G_s+G_t}\right)}^{\eta }=G^c$$

(9)

where G_c is the cohesion energy; G_n is the normal fracture energy; G_s is the fracture energy in the first shear direction; G_t is the fracture energy in the second shear direction; and η is the material parameter.

3. Establishment of the Numerical Model

The radial well is completed by an open hole, and no casing is tripped in the radial well. No perforation hole is required, and open hole fracturing bullets are reapplied. The mechanisms of fracture initiation and propagation in deflagration fracturing are related to the mechanism of rock fracture on the borehole wall caused by explosive deflagration.

3.1. Model Characterization

3.1.1. Assumptions

The underground rock is heterogeneous and anisotropic. The stress field around the borehole is unstable. To facilitate the calculation and analysis of the model, the assumptions are made as follows: a homogeneous and isotropic reservoir; only a single radial well in the formation and a negligible effect of the vertical section; a negligible effect of deflagration at a high temperature on the rock properties; the pore pressure and seepage of the reservoir are not considered due to the short deflagration process; and there are no natural fractures around the wellbore.

3.1.2. Model Settings

Taking a well in the Changguangmiao Block in the Yanchang Oilfield as an example, the parameters of the wellbore and rock mechanics are set as shown in

Table 1. Reservoir parameters.

Parameters	Values	Parameters	Values
Rock density	2500 kg/m3	Radial well diameter	0.03 m
Horizontal maximum principal stress	30 MPa	Horizontal minimum principal stress	28 MPa
Vertical stress	35 MPa	Rock tensile strength	3 MPa
Poisson’s Ratio	0.25	Elastic modulus	12.9 GPa

:

The parameters, such as A, B, R₁, R₂, ω, ρ, and D, of the explosive materials are referenced to the parameters of standard TNT, measured by the cylinder experiment, as shown in

Table 2. Parameters in the JWL state equation.

Parameters	Constants					Deflagration Rate	Explosive Density	Initial Specific Internal Energy
Unit	A/GPa	B/GPa	R1	R2	ω	D/m·s−1	ρ/kg·m−3	e/GPa
Values	371.2	3.231	4.15	0.95	0.3	6930	1600	7.0

[29]:

3.2. Model Establishment

In this study, the JWL state equation is used to simulate the explosion’s initiation. The standard TNT explosive model is established as a cylindrical model using the finite element software. The model has the size dimensions of R_B = 0.03 m and L_B = 18 m and is meshed into 4680 C3D8R hexahedrons, as shown in Figure 1. These meshes are converted into particles according to time steps. After 10⁻⁸ s, the computation starts, and one mesh is converted into one particle. A total of 4680 particles are generated.

In the JWL state equation, the physical parameters are assigned to the TNT model. To simulate the underground stress, the reservoir model is established as a cuboid with the size dimensions of 10 m × 10 m × 20 m and 34,803 hexahedral C3D8R elements, as shown in Figure 2. The horizontal maximum principal stress, the horizontal minimum principal stress, and the vertical stress are assigned in x, y, and z directions, respectively. The radial well extends in the y direction. Cohesive insertion unities are programmed in the Python language, and the global cohesive unit is inserted to simulate fracture initiation and propagation. The infinite element theory is introduced in the computation of the bounce-back of the stress at the boundary. The infinite boundary is defined by the CIN3D8 element. It is assumed that the tensile stress is positive and the compressive stress is negative. The displacement boundary conditions are U_x = 0, U_y = 0 and U_z = 0.

The radial well has the diameter of D_rad = 0.03 m and the length of H_rad = 18 m and is in the center of the reservoir. The radial well has closed ends. The standard TNT explosive is placed in the middle of the radial well, and the center of the charge column is ignited.

3.3. Model Verification

According to the deflagration parameters proposed by Hu Zhaoying et al. [30], the shock wave over-pressure was calculated using the model developed in this study. The calculated shock wave over-pressure reached the peak at the horizontal distances of 3.5 m, 5.0 m, 7.0 m, and 9.5 m from the explosion center. The calculated values were compared with the measured values, as shown in

Table 3. Comparison of the numerical simulation results with the experimental results.

Explosion Center/m	Over-Pressure Peak/MPa		Error/%
	Measured	Calculated
3.5	0.257	0.263	2.3%
5.0	0.124	0.132	6.5%
7.0	0.099	0.106	7.1%
9.5	0.051	0.056	9.8%

. The numerical simulation results are close to the experimental results, and the errors are less than 10%, which verifies the rationality of the numerical model.

4. Analysis of Results

4.1. Effect of the Deflagration Position

A numerical model of single-point deflagration was established in this study. Near-wellbore deflagration, middle deflagration, and toe deflagration were simulated, respectively, by setting the deflagration position at the left end (0, −5, 0), the midpoint (0, 0, 0), and the right end (0, 6, 0), as shown in Figure 3. It is interesting to note that, in a short time, the deflagration of the explosive materials started and extended preferentially along the wellbore until all the explosives were ignited, and then the energy was rapidly radiated as a spherical form and transmitted to the reservoir, as shown in Figure 4.

To obtain the real fracture morphology, the XOY and YOZ sections and the cohesive units corresponding to the deflagration points (as shown in Figure 5) were extracted, and the fracture propagation pattern across the radial well was analyzed. In the cohesive units, three basic fracture morphologies are simulated: type I—open fracture, type II—slip fracture, and type III—rip fracture. The MMIXDMI parameter represents the field output in cases of tensile and shear failure. When the value is between −1~0, rock failure does not occur (blue area); when the value is 0~0.5, tensile failure is dominant (green area); and when the value is 0.5~1, shear failure is dominant (red area). Moreover, the area of the cohesive fracture surface is calculated by post-processing in the Python language, which facilitates the quantitative analysis of the data.

During deflagration, under three deflagration points, part of the reservoir in the area around the radial well (XOY and YOZ planes) in the direction of the horizontal minimum principal stress is fractured completely, and tensile failure is dominant, accompanied by some shear failure. Tensile failure occurs mainly in the area parallel to the direction of the horizontal maximum principal stress (XOZ plane), and almost no shear failure occurs. The fracture propagation mechanism is analyzed quantitatively in five typical interfaces.

The simulation shows that when deflagration occurs in the middle of the radial well, a higher deflagration pressure and larger fracture area are generated in the middle area, and a smaller fracture area is generated at both ends (Figure 6). The fracture is analyzed at the cross-section of the radial well. Around the deflagration point, tensile failure leads to the fracture area of 41.58 m², and the fracture extends in a cross-shape from the deflagration point (green area). Shear failure leads to the fracture area of 58.37 m² (red area) around the tensile fracture. Comparatively, the reservoir stimulation near the deflagration point is dominated by shear failure. As the distance from the deflagration point increases, the shear failure decreases, the total fracture area decreases, and the reservoir stimulation efficiency decreases, as shown in

Table 4. Fracture quantitative characterization at different deflagration points.

Positions	Fracture Area at the Left Deflagration Point/m2			Fracture Area at the Middle Deflagration Point/m2			Fracture Area at the Right Deflagration Point/m2
	Shear Failure Area	Tensile Failure Area	Fracture Area	Shear Failure Area	Tensile Failure Area	Fracture Area	Shear Failure Area	Tensile Failure Area	Fracture Area
Left	57.26	40.13	97.39	16.88	47.91	64.79	4.23	10.34	14.57
Mid	17.26	50.13	67.39	58.37	41.58	99.95	13.47	47.75	61.22
Right	3.288	13.51	16.79	13.15	45.09	58.24	56.54	40.22	96.76

. In the sections that are 5 m and 6 m away from the deflagration point, the shear fracture areas of 17.26 m² and 13.47 m² and the tensile fracture areas of 50.13 m² and 47.75 m² are generated, respectively. The reservoir stimulation areas are reduced by 32.58% and 38.75%, respectively, compared with those at the deflagration point.

When deflagration occurs at the radial well ends, the fracture is also initiated in the same way as the abovementioned mechanism. This indicates that the deflagration point directly affects the fracture propagation pattern. A larger deflagration pressure wave and larger fracture area are generated at the position closer to the deflagration point, and this is favorable for local reservoir stimulation. With the increase in the distance from the deflagration point, the reservoir stimulation efficiency is reduced significantly, and this indicates the characteristics of deflagration fracturing, in which a greater reservoir stimulation area is generated near the radial well, and a lower stimulation volume is generated far from the deflagration point.

4.2. Effect of the Radial Well Azimuth

Assuming that deflagration occurs in the middle of the radial well and that the abovementioned rock mechanic parameters are unchanged, the effect of the radial well azimuth on the deflagration fracturing was simulated by setting the radial well azimuth α (the angle between the radial well axis and the horizontal maximum principal stress) as 0°, 15°, 45°, 75°, and 90°.

The results show that with the increase in the radial well azimuth, the fracture area at the deflagration point increases from 60.09 m² to 99.95 m², with an increment of 66.33%. As the radial well azimuth changes, the stress distribution on the borehole wall and deflagration stress field change accordingly, which affects the fracture propagation pattern. The larger radial well azimuth indicates lower horizontal stress on the section across the wellbore, and the fractures are prone to initiation and propagation, resulting in a larger stimulation area. Moreover, with the increase in the radial azimuth, the shear fracture area increases gradually, as shown in Figure 7. The tensile fracture area reaches a maximum of 61.15 m² when α = 75° and decreases to 41.58 m² when α = 90°. With the increase in the radial well azimuth (α > 75°), the stress in the rock changes, and the fracture mechanism changes from tensile failure to shear failure, as shown in

Table 5. Fracture characterization under different radial well azimuths.

Radial Well Azimuths	Fracture Area in the Middle of the Radial Well/m2
	Shear Failure Area	Tensile Failure Area	Fracture Area
0°	12.64	47.45	60.09
15°	17.17	50.94	68.11
45°	27.43	57.59	85.02
75°	35.68	61.15	96.83
90°	58.37	41.58	99.95

. This provides a reference for understanding the fracture mechanism of deflagration fracturing.

4.3. Effect of the Horizontal Stress Difference

Assuming that deflagration occurs in the middle of the radial well, the abovementioned rock mechanic parameters are unchanged, and the horizontal minimum principal stress is 28 MPa, the effect of the horizontal stress difference on the fracture morphology was simulated by setting the horizontal maximum principal stress as 33 MPa and 36 MPa (Figure 8). When the horizontal stress difference increased from 2 MPa to 8 MPa, the fracture area decreased from 99.95 m² to 71.33 m², with a reduction of 28.63%. The large stress difference had a significant negative effect on the fracture initiation and propagation, as shown in

Table 6. Fracture characterization under variable horizontal stress differences.

Stress Difference/MPa	Fracture Area in the Middle of the Radial Well/m2
	Shear Failure Area	Tensile Failure Area	Fracture Area
2	58.37	41.58	99.95
5	22.98	51.60	74.58
8	16.8	54.53	71.33

. However, in deflagration fracturing, there is still a good capacity for fracture initiation in the rock. With the increase in the horizontal stress difference, the shear failure area decreased rapidly, and the tensile failure area increased to a certain extent. A larger stress difference is conductive to tensile failure.

It should be noticed that with the increase in the horizontal maximum principal stress, the difference between the horizontal maximum principal stress and the vertical stress decreased. The horizontal maximum principal stress of 36 MPa is greater than the vertical stress. According to the theory of solid mechanics, variation in the maximum stress directions of the X, Y, and Z axis leads to variation in the fracture morphology in deflagration fracturing. Thus, further study of fracture initiation and propagation in the cohesive plane is needed. An example in a cohesive unit on the X-Y plane was analyzed here, as shown in Figure 9.

The numerical simulation shows that when the horizontal maximum principal stress is lower than the vertical stress, the fractures on the X-Y plane are dominated by tensile failure. When the horizontal maximum principal stress increases to 33 MPa, a small number of shear fractures occur on the X-Y plane, with the shear fracture proportion of 9.61%. When the horizontal maximum principal stress is larger than the vertical stress, shear failure occurs obviously on the X-Y plane, with the shear fracture proportion of up to 56.2%, as shown in

Table 7. Quantitative characterization of the fracture morphology on the X-Y plane under variable stress differences.

Stress Difference/MPa	Fracture Area in the Middle of the Radial Well/m2
	Shear Failure Area	Tensile Failure Area	Fracture Area
2	0	20.00	20.00
5	1.92	18.08	20.00
8	11.25	8.75	20.00

. It is confirmed that in all the stress distributions, the deflagration horizontal fracture occurs on the plane of the radial well axis, indicating that the stress distribution has little effect on the horizontal fractures, and fracture initiation and propagation in deflagration fracturing are mainly affected by the stress wave generated by deflagration. The stress distribution has an obvious effect on the rock failure type. When the horizontal maximum principal stress and the vertical stress are reversed, fracture initiation in the horizontal direction is more likely to be subject to shear failure.

5. Discussion

In this study, the SPH theory was integrated into the finite element model, and the deflagration reaction was simulated using the JWL state equation. The stress field was introduced to simulate the effects of stress on the radial wells in different directions. The process of fracture initiation in deflagration fracturing was simulated with the cohesive units and solid mechanics theory. This study provided a simulation of the process and mechanism of deflagration fracturing in reservoirs. The previous studies focused on the calculation of the deflagration parameters, including the deflagration load, ignition model, deflagration spacing, correction of the SPH model, and the mechanism of deflagration’s damage to the soil, water, brittle materials, and steel. This study focused on practical applications in well stimulation. Non-convergence in the calculation of 3D fracture propagation is avoided by the cohesive unit. The model provides a higher calculation efficiency and better understanding of the stress around the radial well and the fracture morphology parameters. The simulation shows that the deflagration position, the radial well azimuth, and the horizontal principal stress difference have a significant effects on the fracture formation area and generation of shear and tensile fractures. In deflagration fracturing, the high temperature generated by deflagration has a certain effect on the reservoir fluid and then affects the pore pressure and seepage process. Nevertheless, considering that the deflagration is completed in a very short time (less than 0.1 s) and the stress between the deflagration gas and rock mainly stems from the high-speed pressure wave, the pore pressure and seepage and the effect of the temperature on the rock are not calculated in the mathematical model. Furthermore, the explosive volume, the rock mineral composition, and the radial well diameter also affect the reservoir stimulation effect. Subsequently, in the future, we will study these factors to reveal more information about the mechanism of radial-well-assisted deflagration fracturing.

6. Conclusions

The numerical model of deflagration fracturing was established by the SPH method, and the effects of the deflagration position, radial well azimuth, and horizontal principal stress difference on fracture propagation were simulated. The conclusions were obtained as follows:

• In deflagration fracturing, the deflagration position, the radial well azimuth, and the horizontal principal stress difference have significant effects on the fracture formation area. The closer the distance from the deflagration position is, the larger the radial well azimuth and the smaller the horizontal stress difference are, leading to a larger fracture area, which is conducive to reservoir stimulation. As the distance from the deflagration point increases from 5 m to 11 m, the fracture area decreases by 33.5% and 85.1%, respectively. With the increase in the distance from the deflagration point, the reservoir stimulation efficiency decreases exponentially. When the radial well azimuth increases from 0° to 90°, the fracture area at the deflagration point increases by 66.33%. The larger azimuth is favorable for reservoir stimulation. As the horizontal principal stress difference increases from 2 MPa to 8 MPa, the fracture area at the deflagration point decreases by 28.63%. A higher stress difference has a significant negative effect on fracture initiation and propagation. However, in deflagration fracturing, there is still a good capacity for fracture initiation in the rock.
• In deflagration fracturing, both shear slip fractures and tensile fractures occur. The formation and conversion of shear slip fractures and tensile fractures are related to the deflagration position, radial well azimuth, horizontal principal stress difference, etc. Shear failure is sensitive to the deflagration distance, the radial well azimuth, and the horizontal principal stress. With the increase in the distance from the deflagration point, decrease in the radial well azimuth, and increase in the horizontal principal stress, the shear failure area decreases by 93.39%, 78.31%, and 71.23%, respectively. The tensile failure area is closely related to the deflagration distance, radial well azimuth, and the horizontal principal stress.
• The process of deflagration fracturing can be characterized based on the JWL constitutive equation, the SPH method, and infinite element numerical simulation. The error of the numerical simulation is less than 10%, and the problems of the stress wave rebound and poor calculation accuracy are mitigated.