The escalating demand for oil, coupled with the rapid pace of development, has led to the gradual depletion of conventional oil and gas resources. Consequently, unconventional reservoirs like shale oil and gas are gaining prominence. These reservoirs exhibit small pore radii, low permeability, and significant challenges in oil and gas flow. To effectively exploit these resources, hydraulic fracturing techniques become imperative (Al-Rbeawi and Owayed 2020). Hydraulic fracturing technology involves using high-pressure pump systems on the surface to pump fracturing fluid into the well at a volume exceeding the formation's absorption capacity, generating high pressure inside the wellbore. When this pressure surpasses the near-wellbore geo stresses and reaches the rock's tensile strength, fractures are created in the nearby formation layers. By continuing to inject fracturing fluid containing proppants, the fractures are extended and filled with proppants. Eventually, after the pumping is stopped, proppant-supported fractures with specific geometric dimensions and conductivity are formed within the formation, achieving the goal of reservoir stimulation (Baioco et al. 2019; Qian et al. 2023). Since the first field application of hydraulic fracturing technology in 1947, its technical principle has always followed the concept of reservoir stimulation using polymer fracturing fluids and solid proppant particles.

However, traditional methods often bring forth issues like environmental pollution from fracturing fluids and heightened energy consumption due to the extensive use of injection equipment. In response, a novel fracturing technology known as Self-Propping Phase-Transition Fracturing Technology (SPFT) has recently emerged. This technique employs two immiscible fluids (Including phase-transition fluid (PF) and non-phase-transition fluid (NPF)), termed Phase-Transition Fracturing Fluid Systems (PFFSs), to generate hydraulic fractures. Through temperature stimulation, one of these fluids (PF) undergoes an in-situ phase transition, transforming into solid particles (In-situ Generated Proppant or IGP) to support fractures (Fig. 1). This process results in the formation of hydraulic fractures with specific geometric shapes and conductivities (Zhang et al. 2023a). The complete technical process is: (1) PFFS composed of PF and NPF is used to open the reservoir and form artificial fractures. The injected PFFS cools the reservoir near the fracture below the phase transition temperature (Fig. 1a); (2) After the injection of PFFS is stopped, the reservoir heats the PFFS, and the fracture temperature rises above the phase transition temperature. As a result, the PF undergoes a phase transition, turning from a liquid phase into a solid IGP to prop the hydraulic fracture (Fig. 1b); (3) The fracture space formed after NPF flow back becomes a high-conductivity flow channel for oil and gas (Fig. 1c). In this technology, the role of the PF is to undergo a phase transition to generate IGP under the stimulation of reservoir temperature, the NPF occupies part of the fracture space, and after the flow back of the NPF, these spaces serve as high-conductivity flow channels for oil and gas. The function of the IGP is to support the fractures and prevent them from closing.

The SPFT principle (Zhang et al. 2022)

This technology has two primary advantages. Firstly, the flow-back fluid becomes polymer-free as polymers transform into solid proppant particles. This transformation ensures that pollutants do not enter groundwater, thereby preventing environmental pollution. Secondly, the technology eliminates the need to inject solid proppant or crosslinker, resulting in a substantial reduction in injection equipment. Essentially, this new technology provides an effective solution to issues related to environmental pollution and high energy consumption. The SPFT introduces a new reservoir stimulation concept. The formation mechanism of hydraulic fractures and the construction method of conductivity differ significantly from existing fracturing technology. This involves a complex thermal field, hydraulic field, mechanical field, chemical field, and thermal–hydraulic-mechanical-chemical (THMC) coupling process. The thermal field primarily studies the PFFS phase transition heat and temperature distribution(Zhang et al. 2023b); the hydraulic field focuses on fluid flow in the matrix and fractures (Pei et al. 2021); the mechanical field examines rock mass deformation and hydraulic fracture propagation paths; and the chemical field studies PFFS phase transition kinetics, heat release rate, and total heat release (Zhao et al. 2020).

The field variables in SPFT involve multifield and multimedium coupling, with different physical fields interacting. Scholars have extensively researched these complex multifield coupling problems using mathematical models, numerical methods, and other approaches. Concerning the multifield coupling mathematical model, the initial model is the mechanical and hydraulic fields coupling model based on the effective stress principle (Vinci et al. 2014). Scholars have realized the challenge of addressing rock mechanics in a high-temperature environment by considering only the coupling effect between the mechanical and hydraulic fields, as indicated in the studies pertinent to geological carbon dioxide storage, geothermal development, oil and gas development, etc. (Wang et al. 2022). Ignoring the thermal expansion and thermal stress of rock with temperature change and the two-way interaction between thermal, hydraulic, mechanical, and chemical fields makes it difficult to fully understand rock deformation and damage (Ogata et al. 2018; Tao et al. 2020; Zhang et al. 2023d).

Owing to the complexity of the problems, numerical methods are necessary to solve the THMC-coupled mathematical models. Common numerical methods include continuous, discrete, and quasicontinuous mediums (Zhang et al. 2023c). The continuous medium method includes the finite element method, extended finite element method, or boundary element method (Gupta and Duarte 2018; Mukhtar and Duarte 2022; Pariseau et al. 2019), and the discrete medium method includes the discrete element method, numerical manifold method(Kong et al. 2021; Lin et al. 2020; Zhao et al. 2014). The quasicontinuum method regards the domain as composed of numerous particles bonded by virtual bonds (Gao and Klein 1998; Klein and Gao 1998), with representative methods being virtual internal bonds and discretized virtual internal bonds (Peng et al. 2018; Zhang and Chen 2009; Zhu and Zhang 2023).

Progresses in the mathematical models and numerical methods of multifield coupling are significant for revealing the mechanism of multifield coupling(Peng 2023). The representative theory mainly includes the effective stress principle, convective heat transfer in porous media, solute mass transfer, chemical reaction, mineral dissolution, and its influence on physical parameters and damage of solid medium. The main content of the existing THMC model includes conservation equations of mass, energy, and momentum, incomplete coupling model, fully coupled model, multi-phase flow coupling model, and dual porosity model. Its application scenarios include carbon dioxide storage, nuclear waste storage, geothermal development, and oil/gas resource development (Luo et al. 2018; Zhang et al. 2021). Different from existing multifield coupling problems, the thermal, hydraulic, mechanical, chemical, and THMC coupling processes in SPFT are unique:

(1) Thermal field: Existing field models mainly focus on the formation temperature variation under the action of injected cold fluid and acid rock reaction heat (Mou et al. 2023). In addition to studying the formation temperature in SPFT, attention is needed on the effect of PFFS phase transition heat release on temperature distribution.

(2) Hydraulic field: Available models mainly determine the velocity and pressure distribution of single- and multiphase fluids in the flow process (Liu et al. 2016). For SPFT, attention is also necessary on the leak-off of PFFS from the hydraulic fracture to the matrix (Pei et al. 2019).

(3) Mechanical field: Relevant models of the mechanical field mainly describe the stress distribution of rock mass (Tang et al. 2019), deformation, and failure patterns (Taghichian et al. 2018). However, in SPFT, it is also necessary to consider the in-situ generated proppant (IGP) deformation law under closure stress and the variation law of fracture size.

(4) Chemical field: Existing relevant studies mainly focus on acid-rock reactions, mineral dissolution, and precipitation (Gan et al. 2021). The research object of the chemical field in SPFT is the phase transition rate and conversion rate of PFFS, as well as the chemical reaction heat brought by the phase transition process.

In summary, existing models provide a reference for studying multifield coupling problems in SPFT. However, due to the unique characteristics of the thermal, hydraulic, mechanical, and chemical fields, they fail to accurately reflect the multifield coupling process and mechanism in SPFT, thus failing to determine the geometric shape and size variation of hydraulic fractures under different influencing factors. In this regard, based on the characteristics of the SPFT process and the physical and chemical characteristics of PFFS, this paper elaborated a set of THMC multifield coupling models, developed the FEM-DVIB-EPM coupling algorithm, and verified it on a case study. The ultimate goal was to reveal the mechanism of hydraulic fracture propagation caused by SPFT and provide guidance for hydraulic fracturing technology parameters’ optimization.

2.1 Thermal field mathematical model

Temperature is the key condition to control the phase transition time and phase transition position of PFFS(Zhang et al. 2022). On the one hand, the low-temperature PFFS injected into the reservoir reduces the formation temperature. On the other hand, PFFS will release heat during the phase transition process to increase the formation temperature. Based on the PFFS phase transition reaction kinetic model, a mathematical temperature field model was established, including the temperature in the hydraulic fracture and formation.

The temperature model in the fracture considering heat conduction, heat convection, heat exchange between the fluid in the fracture and the rock mass, and phase transition heat release of PFFS can be expressed as follows:

where \(\lambda_{\text{w}}\) is the thermal conductivity of the fluid in the fracture, \({\text{W/(m}} \cdot {\text{K)}}\); \(T_{\text{w}}\) is the temperature of the fluid in the fracture, K; \(c_{\text{w}}\) is the specific heat capacity of the fracturing fluid in the fracture, J/(kg·K); \(\rho_{\text{w}}\) is the density of the fracturing fluid in the fracture, kg/m³; \(K_{\text{f}}\) is the fracture permeability, m²; \(P\) is the pressure of the fluid in the fracture, Pa; \(\lambda_{\text{r}}\) is the heat transfer coefficient between the fluid in the fracture and the rock mass, W/(m²·°C); \(T_{\text{r}}\) is the temperature of rock mass near the fracture, K; \(m_{\text{pf}}\) is the mass of PFFS contained in the unit, kg; \(\Delta H_{\text{R}}\) is the heat release of PFFS in the chemical reaction process, J/kg; \(\alpha (t)\) is the phase transition conversion ratio at the moment \(t\), dimensionless; and \(t\) is the time, s.

The temperature model of the matrix considering the heat conduction, heat convection, and phase transition heat of PFFS leak-off into the matrix is as follows:

where \(\lambda_{x}\), \(\lambda_{y}\), \(\lambda_{z}\) are the rock thermal conductivities along the respective directions,\({\text{W/(m}} \cdot {\text{K)}}\);.. is the specific heat capacity of the fluid in the matrix, J/(kg \(\bullet\) K); \(\rho_{\text{mw}}\) is the density of the fluid in the matrix, kg/m³; \(\phi\) is the porosity of the matrix, dimensionless; \(v_{x}\), \(v_{y}\), and \(v_{z}\) are the flow velocities in the matrix along different directions, m/s; \(m_{\text{pm}}\) is the mass of the PFFS contained in the unit, kg; \(c_{\text{r}}\) is the specific heat capacity of the rock mass, J/(kg \(\bullet\) K); and \(\rho_{\text{r}}\) is the density of the rock mass, kg/m³.

2.2 Chemical field mathematical model

In Eqs. (1) and (2), \(\alpha (t)\) is the phase transition conversion ratio of PFFS at moment \(t\), and \(\Delta H_{\text{R}}\) is the heat release of PFFS during the chemical reaction process. The phase transition process is a chemical reaction process. It is necessary to establish a chemical field mathematical model of PFFS. The heat release \(\Delta H_{\text{R}}\) of PFFS can be measured by DSC (differential scanning calorimetry) experiments. The TA company's DSC Q20 (Fig. 2) was used to test the phase-transition process of PF. Each group of experiments was carried out under a high-purity nitrogen atmosphere with a flow rate of 50 mL/min. The test method is as follows: weigh about 10 mg of a PF sample into the aluminum crucible, quickly seal it, and put it into the DSC sample room; put a sealed empty aluminum crucible in the DSC sample room as a reference; heat the test sample from 313 K (40 °C) to 523 K (250 °C) at a specific heating rate (2.5, 5, 10, 15, or 20 K/min), and measure and record the heat fluxes at different time points (Table 1).

DSC Q20 instrument

The phase transition conversion rate \(\alpha (t)\) at the moment \(t\) can be fitted by the DSC experimental heat fluxes results(Zhang et al. 2023a). The phase transition reaction kinetic model fitted by the DSC experimental test is as follows:

where \(T\) is the ambient temperature, K.

2.3 Hydraulic field mathematical model

The hydraulic field mainly studies the fluid pressure and velocity in the fracture and the matrix, and the pressure field in the fracture will be used as the pressure boundary condition of the matrix seepage field. For the seepage field in the matrix, based on the DVIB theory, the rock matrix is permeable due to the randomly distributed microfractures and micropores inside it. These microfractures are partially connected to form a microfracture network. Each bond cell in the DVIB model represents a small-volume microfracture network (Fig. 3). Each bond bears the seepage characteristics of each microfracture, and the seepage control equation of each bond is (Peng et al. 2017):

where \(k_{\text{b}}\) is the equivalent permeability coefficient; \(P_{\text{m}}\) is the matrix fluid pressure; \(Q_{{\text{int}}}\) is the source or sink; \(s_{\text{b}}\) is the equivalent water storage rate; \(\rho\) is the fluid density; and \(g\) is the gravitational acceleration.

Fracture medium model via the DVIB theory (Peng et al. 2017)

By assembling the seepage equation of the bond, the matrix form of the seepage equation of each bond cell can be obtained, and then the matrix form of the seepage equation of the whole computational domain can be further assembled (Peng 2018).

During the hydraulic fracturing, high-pressure fracturing fluid is injected into the reservoir from the wellbore. Part of the high-pressure fracturing fluid entering the reservoir or fracture leaks off into the matrix, causing changes in the seepage field, and the other part is retained in the fracture, forcing the fracture to propagate. When studying the fluid velocity and pressure in the fracture, due to the narrow fracture width, the difference between the fluid velocity and pressure in the fracture along the direction perpendicular to the fracture wall is not large, and the flow in the fracture can be simplified to two dimensions. In this paper, the Navier‒Stokes (N‒S) equation is used to establish the fluid flow equation in the fracture considering the leak-off of fracturing fluid:

where \(u\) and \(v\) are the flow velocities in the x and z directions in the fracture, respectively, m/s; \(w\) is the fracture width, m; and \(v_{\text{lm}}\) is the fluid leak-off velocity, m/s.

The fluid leak-off velocity \(v_{\text{lm}}\) can be expressed as (Xue et al. 2018):

where \(k\) is the matrix permeability, m²; \(\mu_{\text{a}}\) is the fluid viscosity, Pa·s; \(p\) is the fracture fluid pressure, Pa; \(p_{\text{e}}\) is the matrix pore pressure, Pa; and \(l_{\text{m}}\) is the leak-off depth, m.

The change rate of fluid momentum in the fracture with time depends on the resultant force acting on the fluid in the unit. The momentum equation of incompressible single-phase non-Newtonian fluid is (Bi 2013):

where \(\tau = \left( {\begin{array}{*{20}c} {\tau_{xx} } & {\tau_{\begin{subarray}{l} xz \\ \end{subarray} } } \\ {\tau_{zx} } & {\tau_{zz} } \\ \end{array} } \right) = \mu_{a} \left( {\begin{array}{*{20}c} {2\frac{\partial u}{{\partial x}}} & {\frac{\partial u}{{\partial z}} + \frac{\partial v}{{\partial x}}} \\ {\frac{\partial v}{{\partial z}} + \frac{\partial u}{{\partial x}}} & {2\frac{\partial v}{{\partial z}}} \\ \end{array} } \right).\)

2.4 Mechanical field mathematical model

The mechanical field is described by the DVIB constitutive model (Zhang and Gao 2012). The DVIB constitutive model is directly derived from the energy potential between microscopic particles. Therefore, it contains the microscopic fracture mechanism and does not need the external fracture criterion (Zhang et al. 2016). According to the principle of DVIB, the stress tensor and elastic tensor are as follows (Wang et al. 2012):

where

where \(\psi\) is the strain energy density of the unit; \(\Phi\) is the energy potential of the bond; \(V\) is the volume of the unit; \(\Delta_{\text{n}}\) is the normal deformation; \(\Delta_{\text{t}}\) is the tangential deformation; and \(D\left( {\theta ,\varphi } \right)\) is the distribution density function of the bond.

The energy potential of the bond \(\Phi\) can be expressed as in (Xu and Alan 1999):

The normal deformation \(\Delta_{\text{n}}\) and tangential deformation \(\Delta_{\text{t}}\) are expressed as:

where \(\delta_{\text{n}}\) and \(\delta_{\text{t}}\) are the normal and tangential displacement eigenvalues, respectively; \(\phi_{\text{n}}\) is the work related to the normal deformation \(\Delta_{\text{n}}\); \(\phi_{\text{t}}\) is the work related to the tangential deformation \(\Delta_{\text{t}}\); \(l_{0}\) is the length of the bond before the deformation; and \({\varvec{\xi}}\) is the initial direction vector of the virtual internal bond.

The DVIB-based model is strain-softening; when the element deformation is too large, the stiffness is degraded to a very small value. Once this occurs, it is equivalent to forming a small fracture surface at the element of stiffness degradation. It is assumed that the small fracture surface passes through the geometric center of the tetrahedral element and is perpendicular to the first principal strain direction. The arrangement of all these small fracture surfaces is regarded as a three-dimensional macroscopic fracture surface (Wang et al. 2012). The following test conditions are selected to test these elements whose stiffness degrades to a very low level:

where \(\varepsilon_{1}\) is the maximum principal strain of the element; \(\lambda\) is a coefficient greater than 1, where \(\lambda\) is equal to 10; and \(\varepsilon_{\text{t}}\) is the strain value corresponding to the peak stress in the uniaxial tensile stress–strain curve.

2.5 THMC coupling process

According to the interaction between field variables, the multiphysics coupling analysis method can be divided into two categories: one-way coupling and two-way coupling. According to the treatment of the control equations of each field variable, it can be divided into direct coupling and indirect coupling. Among them, one-way coupling refers to only analyzing the influence of variable A on variable B, without considering the influence of variable B on variable A. Two-way coupling analyzes the interaction between variable A and variable B. Direct coupling refers to putting the thermal field, hydraulic field, mechanical field, and chemical field in the same set of control equations and then solving the control equations with boundary conditions. Indirect coupling calculates the thermal field variables, hydraulic field variables, mechanical field variables, and chemical field variables, which are output separately through iterative calculation until the field variables meet a certain calculation accuracy (Yan et al. 2020).

The direct coupling solution method, where each field variable uses the same control equations and multiple degrees of freedom are solved at each node, leads to high computational difficulty, low efficiency, and large computational cost. Especially in the case of a three-dimensional model and four sets of field variable coupling, the degrees of freedom will increase exponentially, and it is challenging to achieve a direct coupling solution. Therefore, the indirect one-way and two-way coupling methods are used to solve the THMC coupling model. The model coupling principle is illustrated in Fig. 4.

THMC coupling process

Equations (1) to (3) describe the coupling of the thermal and chemical fields. The thermal field and mechanical field are mainly coupled by thermal stress. Because the influence of thermal stress is small, it is not the focus of this paper. Still, the influence of the temperature shown in Eq. (21) on the PFFS rheology is considered, indirectly affecting the fracture size. The hydraulic field affects the chemical field through the flow of PFFS. The mechanical field changes the permeability through the deformation of the rock mass and the formation of artificial fractures, changes the distribution range of PFFS, and finally acts on the chemical field. The influence of the hydraulic field on the mechanical field and that of the thermal field on the hydraulic field are emphasized below.

2.5.1 The hydraulic field effect on the mechanical field

The fluid pressure acting on the fracture surface will cause the expansion of the hydraulic fracture. In turn, the expansion of the hydraulic fracture will increase the fracture opening, change the permeability of the fracture, and affect the hydraulic field. The fluid pressure is applied to the partition element node as an equivalent node force. In the four-node tetrahedral element, according to the shape of the contact surface between the bond cell and the fracture, the partition element can be divided into two types of (Huang and Zhang 2010): triangular partition element (Fig. 5a) and quadrilateral partition element (Fig. 5b).

The method of applying equivalent nodal force

The fluid pressure is directly applied to the node in the calculation process. The force acting on the fracture surface \({\varvec{F}}\) is expressed as:

where \({\varvec{P}}\) is the fluid pressure and \(A\) is the partition surface area.

In the triangular partition element, the force components \(\left( {F_{mx^{\prime}} ,F_{my^{\prime}} ,F_{mz^{\prime}} } \right)\) of the nodes in the x-, y-, and z-directions are equal to \({\varvec{F}}\), that is:

The relationships between the force components \(\left( {F_{ix^{\prime}} ,F_{iy^{\prime}} ,F_{iz^{\prime}} } \right)\), \(\left( {F_{jx^{\prime}} ,F_{jy^{\prime}} ,F_{jz^{\prime}} } \right)\), and \(\left( {F_{kx^{\prime}} ,F_{ky^{\prime}} ,F_{kz^{\prime}} } \right)\) of the \(I\), \(J\) and \(K\) nodes in the three x-, y-, and z-directions are derived as in (Huang et al. 2013):

Similarly, in the quadrilateral partition element, the equivalent nodal force can be expressed as:

2.5.2 The thermal field effect on the hydraulic field

The thermal field can change the fluid’s viscosity, affecting the hydraulic field. Therefore, the coupling model of the thermal field and the hydraulic field can be characterized by the empirical viscosity-temperature model:

The mathematical relationship between the density and temperature of PFFS is:

The mathematical relationship between the specific heat capacity and temperature of PFFS has the following form:

where \(\mu_{\text{PF}} \left| {_{T} } \right.\) is the viscosity of PFFS at temperature \(T\); \(\text{Den}\left| {_{T} } \right.\) is the density of PFFS at temperature \(T\); \(\text{SHC}\left| {_{T} } \right.\) is the specific heat capacity of PFFS at temperature \(T\).

3.1 Solution method

As shown in Sect. 2.5, solving the THMC coupling model directly is problematic. Therefore, this paper uses indirect one-way and two-way coupling methods to solve the THMC coupling model. The finite element method (FEM) solves the thermal field to reflect the complex boundary of the domain better. The hydraulic field is solved based on the discretized virtual internal bonds (DVIB) theory and the FEM to reduce the solution workload and improve the calculation accuracy fully. The DVIB theory and the element partition method (EPM) solve the mechanical field. The chemical field model via Eq. (3) is a simple polynomial, and the value of the field variable can be obtained by substituting the temperature and time parameters, which is unnecessary to elaborate. This section mainly establishes the solution method of the thermal, hydraulic, and mechanical fields.

3.1.1 Algorithm for solving the thermal field model

The FEM is used to solve the thermal field. The finite element format of the temperature field in the fracture derived via Eq. (1) is as follows:

where \({\varvec{R}}\) is the heat conduction matrix; \({\varvec{S}}\) is the heat capacity matrix; \({\varvec{T}}_{{\varvec{\text{w}}}}\) is the heat capacity matrix; and \({\varvec{F}}\) is the temperature load vector.

The domain is discretized into a unit cell, and the temperature field of the unit cell \({\varvec{T}}_{{\varvec{\text{w}}}}^{\text{e}}\) is expressed as the interpolation relationship of the node temperature, that is:

where \({\varvec{N}}\left( {x,y,z} \right)\) is the interpolation matrix, \({\varvec{N}}\left( {x,y,z} \right) = \left[ {N_{1} ,N_{2} ,N_{3} ,N_{4} } \right]\), and \({\varvec{q}}_{T}^{\text{e}}\) is the node’s temperature matrix.

The interpolation matrix \({\varvec{N}}\left( {x,y,z} \right)\) can be expressed as:

where \(V\) is the unit’s volume, while \(a_{i}\), \(b_{i}\), \(c_{i}\), and \(d_{i}\) are the coefficients related to the node’s geometric position.

In Eq. (22), it is necessary to integrate time to derive the unknown quantity \({\varvec{T}}_{{\varvec{\text{w}}}} \left( {t_{{n{ + }1}} } \right)\) at the next moment \(t_{{n{ + }1}}\) from the known quantity \({\varvec{T}}_{{\varvec{\text{w}}}} \left( {t_{n} } \right)\) at the current moment \(t_{n}\). In a short time step, \({\varvec{F}}\) is also known. The integration of time continues until the time point is sought. Therefore, time integration is also a recursive relationship. The time interval \(\Delta t\) can be regarded as a time unit. The integration of time also discretizes Eq. (22) in time. Equation (22) is processed by the trapezoidal method, and the following results are obtained:

The time recursive relationship is derived via Eq. (25), thus avoiding derivative calculation and discretizing the time term. The matrix temperature field obtained via Eq. (2) can be calculated discretely by the same method.

3.1.2 Algorithm for solving the hydraulic field model

The hydraulic field includes the seepage equation in the matrix shown in Eq. (4) and the N‒S equation in the fracture shown in Eqs. (5) to (8). The DVIB theory is used to solve the seepage equation in the matrix. Each bond is regarded as a unit; then, the discrete form of Eq. (4) is (Peng et al. 2017):

where \({\mathbf{H}}^{{\text{b}}} = \frac{{k_{\text{b}} }}{l}\left[ {\begin{array}{*{20}c} 1 & { - 1} \\ { - 1} & 1 \\ \end{array} } \right]\), \({\mathbf{S}}^{{\text{b}}} = \frac{{s_{\text{b}} l}}{\rho g}\left[ {\begin{array}{*{20}c} \frac{1}{3} & \frac{1}{6} \\ \frac{1}{6} & \frac{1}{3} \\ \end{array} } \right]\), \({\mathbf{Q}}^{{\text{b}}} = \left[ {\begin{array}{*{20}c} {q_{i} } \\ {q_{j} } \\ \end{array} } \right]\), \({\mathbf{P}}^{\text{b}} = \left[ {\begin{array}{*{20}c} {p_{i} } \\ {p_{j} } \\ \end{array} } \right]\), \(\mathop {\mathbf{P}}\limits^{{.}{\text{b}}} = \left[ {\begin{array}{*{20}c} {\dot{p}_{i} } \\ {\dot{p}_{j} } \\ \end{array} } \right]\).

By assembling the matrix form seepage equation of the bond, the matrix form of the seepage equation of each bond cell can be obtained, and then the seepage eqsuation of the whole discrete domain can be further assembled.

The number of virtual internal bonds in the ideal bond cell is sufficient, and the virtual internal bonds are evenly distributed in the bond cell. Therefore, the calculation formulas of \(k_{\text{b}}\) and \(s_{\text{b}}\) in the ideal state derived from the ideal bond cell have a certain deviation from the actual situation. In the practical application of the bond’s equivalent seepage parameters, the bond cell’s geometric correction parameters should be introduced to improve the model’s accuracy. After correction, \(k_{\text{b}}\) and \(s_{\text{b}}\) are derived as follows:

where \(k_{\text{m}}\) is the matrix permeability; a is the fracture opening; \(\Omega\) is the total number of virtual internal bonds; \(V\) is the volume of the bond cell; \(S_{\text{m}}\) is the water storage rate of the matrix; \(S_{\text{f}}\) is the water storage rate of the fracture; \(L\) is the sum of the bond length in the bond cell; \(\lambda_{1}\), \(\lambda_{2}\), \(\lambda_{3}\), and \(\lambda_{4}\) are 0.83, 0.43, 1, and 0.45, respectively.

The fracture width is narrow, generally within 1 cm, the fracture height can attain the order of several ten meters, and the fracture length even exceeds 100 m. Therefore, the fluid velocity and pressure difference along the fracture width direction can be neglected. The flow in the fracture is simplified into two dimensions. In the two-dimensional case, the FEM is used to solve the fluid flow equation in the fracture shown in Eqs. (5) ~ (8), which can better reflect the complex boundary of the domain and improve the calculation rate. Therefore, the four-node quadrilateral mesh is used to solve the fluid flow model in the fracture. The velocity and pressure at any point in the element can be expressed as:

where \({\varvec{u}}_{I}^{\text{e}}\) and \({\varvec{v}}_{I}^{\text{e}}\) are the vectors composed of the velocity of each node in the x and y directions, respectively; \({\varvec{p}}_{I}^{\text{e}}\) is the vector composed of the pressure of each node; and \({{\varvec{\Psi}}}\) is the interpolation function.

Therefore, Eq. (5) can be written as:

The matrix form of Eq. (29) is:

where

The numerical integration of Eq. (31) can be solved using the Gauss integral. The same method can be used to solve the fluid momentum equations in the fractures shown in Eqs. (7) and (8).

3.1.3 Algorithm for solving the mechanical field model

The EPM (Huang and Zhang 2010; Wang et al. 2012) is used to construct the stiffness matrix of the element through which the fracture passes (i.e., partition element). For the tetrahedral element, the displacement vector and the force vector of the element node can be expressed as:

For the triangular partition element and quadrilateral partition element, the element stiffness matrices are (Huang and Zhang 2010):

where \(A_{{im^{\prime } }}\), \(A_{{jm^{\prime \prime } }}\), \(A_{km\prime \prime \prime }\), \(A_{im\prime }\), \(A_{i\prime k}\), \(A_{j\prime m}\) and \(A_{jk\prime }\) are the areas of the control area of the contact point pair; \(k_{\text{n}}\) and \(k_{\text{s}}\) are the normal and tangential stiffnesses of the contact surface, respectively; and \(\delta_{i,j}\) is the Kronecker symbol.

When the EPM is used to simulate the initiation and propagation of fractures in a jointed rock mass, the whole research object is directly meshed instead of changing the meshing scheme for setting joint elements separately. After meshing, the fracture surface and the element node position coordinates are used to determine whether the element intersects with the fracture. If they intersect, the element is automatically converted into a partition element, and the stiffness matrix of the partition element is directly used during the calculation. If they do not intersect, the stiffness matrix of the conventional tetrahedral element is used.

3.2 Solving procedure

According to the above solution method, the proposed coupling calculation process is shown in Fig. 6. The solution process mainly includes the following steps:

1. (1)

   Meshing nodes, units, and input basic parameters.

2. (2)

   The mechanical field and hydraulic field calculation programs are called in turn to calculate the stress, strain, fracture geometry, size, fluid pressure in the fracture, and pore pressure.

3. (3)

   The chemical field calculation program was used to calculate the conversion ratio and heat release of PFFS based on the kinetic equation of the chemical reaction.

4. (4)

   The thermal field calculation program derives the fracture and matrix temperatures.

5. (5)

   Calling the mechanical field, hydraulic field, chemical field, and thermal field calculation program again, update the stress, strain, fracture geometry/size, fluid pressure, pore pressure, conversion ratio/heat release of PFFS, fracture temperature, and matrix temperature.

6. (6)

   If the difference in the thermal field, hydraulic field, mechanical field, and chemical field calculated by steps (2) ~ (4) and that calculated by step (5) is less than 5%, the calculation of the current time step will be ended. Otherwise, step (5) will be repeated until the difference is less than 5%.

7. (7)

   If the current calculation time equals or exceeds the total time, the calculation will end. Otherwise, the next time step will be calculated, and steps (2) ~ (6) will be repeated.

THMC coupling model solving process

3.3 Model verification

Because the current THMC coupling model cannot fully reflect the thermal, hydraulic, mechanical, and chemical field coupling processes involved in the SPFT, it is difficult to use a single case to verify the model established in this paper. Therefore, the thermal, hydraulic, mechanical, and chemical fields are verified in four cases.

3.3.1 Thermal Field Model Verification

The established thermal field model is used to simulate the temperature in the fracture and matrix. During the fluid injection process, the bottom hole temperature gradually decreases. The bottom hole temperature at different times is calculated by the conventional wellbore temperature field model (Luo et al. 2020b). The bottom hole temperature value at each time is taken as the boundary condition for the solution of the thermal field in the fracture and matrix. The main parameters used in the calculation are shown in Table 2.

The temperature in the fracture and matrix calculated by the thermal field mathematical model established in this paper is compared with the results calculated by Luo et al. (Luo et al. 2020b). Based on the principle of energy conservation, Luo et al., established the heat transfer equation in hydraulic fractures, and coupled with the pseudo-three-dimensional fracture propagation model to construct a pseudo-three-dimensional fracture temperature field model, and introduced the influence of temperature on fluid viscosity. The biggest difference between the model established in this paper and the research results of Luo et al. is that the two-dimensional flow of fracturing fluid in the fracture and the influence of temperature on fluid density and specific heat capacity are considered, which is more in line with the actual situation of the fracturing process. The temperature in the fracture varies greatly at different times. The temperature of the fracture front is equal to the original reservoir temperature, and the temperature at the fracture mouth is close to the bottom hole temperature (Fig. 7). The mathematical model established in this paper is a three-dimensional model that considers the influence of leak-off and two-dimensional flow of fluid in the fracture. The isotherm appears semi-ellipse, while the two-dimensional temperature field model assumes that the fracturing fluid flows in the fracture in one dimension. The temperature is equal on the cross-section of the central axis perpendicular to the length direction of the fracture (Fig. 8). Therefore, the calculation results of the two models are different. When the established model is reduced to a one-dimensional flow state, the calculation results are close to the reference model’s (Fig. 9).

Temperature distribution in the fracture at different moments calculated by the proposed model

Temperature distribution in the fracture at the end of fluid injection calculated by the one-dimensional model

Comparison of the results calculated via different models

For the temperature of the matrix, the comparison model only considers the heat exchange between the fracturing fluid and the matrix but does not consider the thermal convection caused by the flow of the seepage fluid. During the fluid injection process, the fluid pressure in the fracture is much larger than the matrix fluid pressure. After the fracturing fluid leaks into the matrix, the pore pressure of the reservoir increases, and the seepage velocity in the vicinity of the fracture is significantly accelerated. The thermal convection caused by the seepage velocity cannot be ignored (Ma et al. 2022). Under the combined influence of the seepage velocity of the matrix fluid and the heat conduction of the matrix, the temperature influence range in the direction of the leak-off is significantly increased, reaching approximately 1 m, exceeding the influence range of about 0.5 m calculated via the reference model (Fig. 10).

Temperature of the matrix along the direction of leak-off at the fracture mouth position

3.3.2 Chemical field model validation

The chemical field model is verified by the phase transition conversion ratio and temperature data measured by DSC experiments. Figure 11 shows the measured data of the PFFS conversion rate at a heating rate of 2.5 K/min and the conversion rate calculated by Eq. (3). The relationship between the phase transition conversion ratio and the time calculated by the kinetic model meets the experimental phenomenon. The phase transition conversion ratio calculated by the model is faster at the initial stage, which is different from the DSC test data. Still, the difference is small, which can meet the needs of numerical calculation and field engineering applications.

PFFS conversion rate versus time (heating rate 2.5 K/min)

3.3.3 Hydraulic field model validation

In computational fluid dynamics (CFD), the cavity lid-driven fluid model is usually used to test the accuracy of the N‒S equation numerical method. The numerical results of Ghia et al. (Ghia et al. 1982) are currently recognized benchmark solutions. Figure 12 shows the physical model of the lid-driven fluid cavity model. The cavity’s left, lower, and right sides are static walls and the upper wall moves to the right at a driving rate \(u_{\text{sc}}\). Under the action of the top driving rate \(u_{\text{sc}}\), several circulations are generated in the cavity.

Physical model of the square cavity lid-driven fluid model

For the mathematical model’s verification, the side length of the square cavity was set to 1 m, the fluid viscosity was 1 mPa·s, the density was 998 kg/m³, the upper wall driving rate was 1 m/s, and the Reynolds number was 1000. The calculated velocity components u on the vertical centerline and the velocity component v on the horizontal centerline were compared with the results of the benchmark solution (Ghia et al. 1982) and the finite volume method (FVM) (Dai 2014). The obtained calculation results of the velocity component u on the vertical centerline (Fig. 13) and the velocity component v on the horizontal centerline (Fig. 14) of the square cavity only slightly differed from the benchmark solution, being also close to the results calculated via the FVM. On the vertical centerline of the square cavity, the minimum velocity had a certain gap with the reference solution. The minimum velocity calculated by the model in this paper was approximately -0.2 m/s, the minimum velocity calculated by the FVM was about -0.3 m/s, and the minimum velocity in the benchmark solution was approximately -0.4 m/s. The difference was within the acceptable range in the whole square cavity, and the calculation results of the velocity component on the horizontal centerline were smaller than those of the benchmark solution and the FVM, meeting the accuracy requirements.

The calculation results of the velocity component u on the vertical centerline via the proposed model, FVM, and benchmark solutions

The calculation results of the velocity component v on the horizontal centerline via the proposed model, FVM, and the benchmark solutions

3.3.4 Validation of the mechanical field model

The 3D physical model with preset fractures shown in Fig. 15 was used to simulate the three-dimensional mechanical field by using this mechanical model and the XFEM model to verify the accuracy of the established mechanical field model. The size of the three-dimensional physical model with preset fractures was 1 m × 1 m × 1 m. The fracture was located in the xz plane and 0.5 m in the y-direction. The strain calibration results corresponding to the stress peak of the model material during uniaxial tension and compression were 4 × 10^–4 and 3.3 × 10^–3, respectively. The size of the preset fracture was 1.0 m × 0.1 m. The stresses applied in the x-, y-, and z-directions were 45, 40, and 50 MPa, respectively. The viscosity of the injected fluid was 50 mPa·s, and the injection flow rate was 5 × 10^–6 m³/s.

Three-dimensional physical model with preset fractures

The fracture propagation length calculated by the model established in this paper was similar to that of the XFEM model (Luo et al. 2020a) (Fig. 16). However, the XFEM model treated the fluid flow in the fracture as one-dimensional flow, while the proposed model considered the fluid flow in the two directions of fracture length and fracture height. Therefore, the fluid pressure distribution in the hydraulic fracture differed, resulting in different fracture propagation velocities, the fracture propagation rate calculated by the proposed model was higher in the initial stage and lower in the later one.

Fracture length evolution

4.1 Coupling analysis of THMC under different injection modes

According to the existing fracturing technology, combined with the physical properties of PFFS, three injection modes suitable for SPFT were designed as in (Zhang et al. 2022):

(1) PFFS is directly used for fracturing and prop, and finally, the displacement fluid is injected to displace all PFFS into the reservoir.

(2) The conventional fracturing fluid forms a fracture, and then the PFFS is used to prop the fracture. Finally, the displacement fluid is injected to displace all PFFS into the reservoir.

(3) Using conventional fracturing fluid to form a fracture, then using PFFS to prop the fracture, and finally injecting sand-carrying fluid to plug the fracture mouth to prevent the flow back of lighter density IGP, and finally injecting displacement fluid.

According to the above three injection modes, the main parameters shown in Table 3 are used to calculate the temperature, stress, fracture geometry, and other parameters under different injection modes. The bottom of the established physical model is a fixed displacement in the three directions of x, y, and z. The physical model size is 200 m (x-direction) × 100 m (y-direction) × 100 m (z-direction), the preset fracture is in the y–z plane, located at y = 50, and the z-direction is in the range of 30 to 60 m. The vertical stress is applied at the top, the maximum horizontal principal stress is applied in the x-direction, and the minimum horizontal principal stress is applied in the direction.

In each injection mode, the volume of the injected fluid is shown in Table 4. Mode A is calculated in two cases: the total fluid volume of mode A-1 is equal to that of mode B and mode C, and the volume of PFFS of mode A-2 is equal to that of mode B and mode C to eliminate the interference of total fluid volume and PFFS volume on the result analysis.

Under different injection modes, the displacement difference in the x-direction is slight (Fig. 17), i.e., between -3 and 3 cm. Under the action of compressive stress in the x direction, the displacement at the fracture mouth is positive, and that at the far end of the fracture is negative. Because the bottom surface of the physical model is fixed, the displacement of the bottom surface is 0, and the top surface has the largest x-direction displacement value. Near the fracture mouth, the fluid pressure is large, causing rapid changes in displacement. The displacement cloud diagram shown in Fig. 17 is sliced, and the x-direction displacement slices under different injection modes shown in Fig. 18 are obtained. Generally speaking, stress concentration occurs at the tip of the fracture, whether it is the tip of the length direction or the tip of the height direction, so the position where the displacement changes rapidly is the location of the fracture. As shown in Fig. 19, the stress concentration range shows the fracture location.

The x-direction displacements under different injection modes (unit: cm)

The x-direction displacement slices under different injection modes

The y-direction stress distributions under different injection modes

The displacement variation law in the y and z directions is similar to that in the x-direction. Under compressive stress, the displacement value in the y-direction is positive on one side and negative on the other, and the displacement near the fracture changes rapidly (Fig. 20a). The z-direction is slightly different. Because the bottom surface displacement is fixed, the displacement of the bottom of the model in the z-direction is zero, and the maximum displacement value in the z-direction is at the top of the model. The displacement value is negative because its direction points to the negative z-axis (Fig. 20b).

Displacements in the y- and z-directions and displacement slices under injection mode B

Under different injection modes, the fracture length and height are different even when the total injection volume is the same (Figs. 21 and 22). When mode A-1 and mode A-2 are used, only PFFS is injected except for the displacement fluid. The injection process is relatively simple. However, due to the leak-off of PFFS (Zhao et al. 2020), the efficiency of the fracturing fluid is low, and the fracture length and height are small. When mode B is adopted, the volume of injected PFFS is reduced, the amount of conventional fracturing fluids increases, and the injected conventional fracturing fluids help reduce the PFFS leak-off and improve the fracturing fluid’s efficiency. The fracture length is 20 m higher than that of mode A-1 and 32 m higher than that of mode A-2, reaching 120 m. In mode C, sand-carrying fluid is added to prevent the flowback of IGP with low density in addition to conventional fracturing fluids. The workload of fluid preparation is the largest, and the injection process is the most complicated. The volume of conventional fracturing fluids is smaller than that in mode B, the leak-off reduction effect is also slightly worse, and the fracture length is approximately 116 m.

Fracture geometries under different injection modes

Fracture sizes under different injection modes

In the injection process, the fluid with a lower temperature enters the fracture and exchanges heat with the matrix. The temperature of the fluid in the fracture gradually increases, and the temperature of the matrix near the fracture gradually decreases (Zhou et al. 2022). Under different injection modes, the temperature variation in the fracture is similar (Fig. 23). The distance from the fracture length to the fracture mouth (i.e., zero point) is lower, and the closer to the fracture tip, the higher the temperature. Figure 24 shows the temperature distribution in the direction of the dimensionless fracture length; that is, the fracture coordinates of each injection mode in Fig. 23 are divided by the fracture length, and the dimensionless fracture length with an interval of 0, 1 is obtained. The temperature distribution on the dimensionless fracture length is drawn. Because the PFFS chemical reaction process will release a certain amount of heat and the PFFS volume injected under the condition of mode A-1 is the largest, the temperature in the fracture is the highest, which is approximately 4 ~ 5 °C higher than that of other injection modes. Injection mode A-2 injects the least liquid, the cooling effect on the reservoir is weak, and the reservoir temperature is slightly higher than in modes B and C.

Temperature in the fracture under different injection modes

Temperature distribution in the dimensionless length direction of the fracture under different injection modes

In summary, from the perspective of the hydraulic fracture size under different modes, the fracture formed by mode B is the longest, and the injected fluid is only conventional fracturing fluid and PFFS, without any solids, which is the best choice for fracturing. The fracture size formed by mode C is the second, but mode C needs to inject sand-carrying fluid, which increases the difficulty of fluid preparation and the complexity of the fluid injection process. The density of IGP, fracture closure rate, and phase transition time should be controlled as much as possible to avoid injecting sand-carrying fluid. For mode A-1 and mode A-2, except for the displacement fluid, only the PFFS is injected in the whole process, and the fluid preparation and injection procedures are the simplest. However, due to the leak-off of PFFS, the fracture length is small, which can be used as an alternative when matching the fracture size requirements.

4.2 Coupling analysis of THMC under different PFFS volumes

To study the effect of the amount of PFFS on the shape and size of fractures, based on the basic parameters shown in Table 3, injection mode B was adopted, that is, using conventional fracturing fluid to form the fractures, then using PFFS to prop the fractures, and finally injecting displacement fluid. The PFFS volume changed to 50, 90, 130, and 170 m³, and the coupling analysis of THMC was carried out.

According m the fracture geometry (Fig. 25) and fracture sizes (Fig. 26) under different PFFS volumes, larger volumes of the injected fluid corresponded to longer fractures. However, under the influence of leak-off, when the fracture height was constant, longer fractures corresponded and larger leak-off areas and volumes, the fracture length increment with increasing fluid volume gradually dropped from 12 m to 8 m, and then to 4 m. The change in the fracture height was less pronounced, and the difference from the reservoir thickness was not large. At PFFS volumes of 50 m³ and 90 m³, the fracture height was approximately 30 m. At PFFS volumes of 130 m³ and 170 m³, it was about 34 m.

Fracture geometries under different volumes of PFFS

Fracture length and height variations with PFFS volume

The temperature was found to be the key factor controlling the phase transition time and position of the PFFS, and it was crucial to the success or failure of the SPFT implementation (Zhang et al. 2023b). Figure 27 shows the temperature in the fracture central axis for different PFFS volumes. Figure 28 depicts the temperature distribution in the dimensionless length direction of the fracture at different PFFS volumes. On the one hand, the injected low-temperature fluid cooled down the reservoir; conversely, it increased the reservoir temperature by chemical reaction heat. From the calculation results shown in Figs. 27 and 28, the cooling effect of the injected PFFS was stronger than the heating effect of the reaction heat. As the amount of injected PFFS increased, the temperature in the fracture decreased. The temperature in the fracture at the PFFS volume of 50 m³ was more than 6 °C higher than that at 170 m³. The temperature at the fracture mouth was approximately 60 °C lower than that at the fracture tip, greatly influencing the PFFS phase transition. Therefore, the fracture can be divided into three sections according to the temperature. In each section, a PFFS system with different phase transition times and temperatures should be used to realize rapid phase transition in the fracture and prevent fracture closure.

Temperature in the fracture central axis for different PFFS volumes

Temperature vs. the dimensionless distance in the fracture length direction for different PFFS volumes

4.3 Coupling analysis of THMC under different displacements

To study the influence of injection displacement on the size and shape of fractures, injection mode B was adopted based on the basic parameters listed in Table 3. The displacements were set to 3, 5, 7, and 9 m³/min, and the coupling analysis of THMC was carried out. Figure 29 and Fig. 30 illustrate the fracture geometries and sizes under different displacements, respectively. When the displacement was low (3 m³/min), the injection time was long, and the leak-off volume was large. Therefore, a large amount of fracturing fluid leaked into the formation, and the fracturing fluid efficiency was low(Jafari et al. 2021), resulting in smaller fracture lengths and heights, which were about 96 and 30 m, respectively, at a displacement of 3 m³/min. With increased displacement, the injection time and the total leak-off volume were reduced, the fracturing fluid efficiency was high, and the length and height of the formed fracture were both large. When the flow rate increased to 9 m³/min, the length and height of the fracture reached approximately 132 m and 40 m.

Fracture geometries under different displacements

Fracture length and height variations with displacement

Figures 31 and 32 show the temperature distributions in the fracture central axis and dimensionless length directions under different displacements, respectively. For equal total amounts of injected fluid, the greater the displacement, the more fracturing fluid enters the reservoir at the same time, the better the cooling effect on the reservoir, and the lower the heating rate of the fracturing fluid. The temperature in the fracture was lower for large displacements: at 9 m³/min, it was lower than that at 3 m³/min by more than 7 °C.

The temperature in the fracture central axis under different displacements

Temperature vs. dimensionless distance in the fracture length direction under different displacements

Self-Propping Phase-Transition Fracturing Technology (SPFT) involves the coupling of thermal, hydraulic, mechanical, and chemical (THMC) components. In this study, a comprehensive FEM-DVIB-EPM method was proposed, integrating thermal, hydraulic, mechanical, and chemical field models. The method's effectiveness was verified in a case study, shedding light on the multifield coupling process under varying injection modes, fluid volumes, and displacements. The key findings include:

1. (1)

   The utilization of indirect one- and two-way coupling methods in solving the THMC coupling model proved advantageous, reducing solution complexity, computational workload, and enhancing accuracy and efficiency.

2. (2)

   The optimal injection mode entails initiating fractures using conventional fracturing fluid, propping fractures with PFFS, and displacing all PFFS into the reservoir through injection of displacement fluid.

3. (3)

   The study revealed that an increase in the PFFS amount corresponds to a higher fracture length. However, the incremental impact on fracture length diminishes with an increase in phase transition fluid. Specifically, at PFFS amounts of 130 and 170 m³, the difference in fracture lengths was negligible. Under the case study parameters, 130 m³ of PFFS was sufficient to form a longer fracture.

4. (4)

   The larger the displacement, the greater the fracture height, although the incremental effect on height diminishes with increased displacement. Excessive displacements may result in an uncontrollable increase in fracture height. Optimal fracture geometry, allowing control over height and preventing longitudinal reservoir penetration, was achieved at a displacement of 5 m³/min.