# Transient Flow Theory of Multiple-Fractured Horizontal Wells with Complex Mechanisms in Shale Gas Reservoirs

Authors: Dianfa Du; Genkai Zhang; Yanwu Zhao; Xiaofei Sun; Bin Zhang
Journal: Geofluids
 (2020)

---

## Abstract

Shale reservoirs have the characterizations of low porosity, low permeability, and hydrocarbon organic matter self-generation and self-storage, resulting in its complex flow mechanisms. Compared with fractured vertical wells, multiple-fractured horizontal wells are widely used due to their advantages of effectively increasing the well control range and further expanding the drainage area. To further study the multiscale flow mechanisms of shale gas, a flow model was established that considered viscous flow in microfractures and inorganic pores, the diffusion of Knudsen in nanoscale porosity, the coexistence of slippage, adsorption-desorption effects under infinity, and closed outer boundary conditions; based on the continuous point source solution, a multiple-fractured horizontal well flow model was established and solved by MATLAB programming. Then, the effects of various factors were investigated. The results show that the Knudsen diffusion and slippage coefficients mainly affect the apparent permeability of the matrix pores. The more the Knudsen diffusion and slippage coefficients are, the earlier the turbulent flow occurs and the higher the gas production is.


---

## Body

## 1. Introduction

The revolution of shale gas in the United States of America (USA) is changing the energy structure of the world [1, 2]. Based on its various forms in the reservoir, the shale gas can be divided into free shale gas, adsorbed shale gas, and dissolved shale gas. Most of the shale gas is deposited in pores and microfractures in a free state, and the proportion of the gas adsorbed on the surface of organic matter and clay mineral particles [3] could generally reach more than 50% [4]. The dissolved gas exists in liquid hydrocarbons or adsorbed on the surface of other substances in the kerogen [5–7] as shown in Figure 1. Through tank degassing experiments, Cheng et al. [8] observed the diffusion of gas from the kerogen or clay particles to the surface of the pores.Figure 1
Schematic diagram of gas molecule in a small part of a kerogen grain pore system of a mudrock.Scholars have committed to the research on shale gas flow models since a decade ago. In 2009, a mathematical model was developed for horizontal wells with triple media branches [9], and the bottom well pressure solution of the horizontal well with the triple media branches in Laplace space was solved by using the point source function theory. Soon, a pressure calculation formula for any point in a stratum with fractures at different inclinations in a fractured horizontal well was developed based on the same theory [10]. Then, Ozkan et al. [11] applied the trilinear flow model to fractured horizontal wells, which is suitable for shale reservoirs with a matrix permeability of only micron-scale Darcy and nanoscale Darcy. This model assumes that the artificial vertical fractures are formed after fracturing, and the surrounding of the fractures is considered to be a dual medium.In the same year, Stalgorova and Mattar [12] proposed a trilinear flow model to analyze the production dynamics of conventional reservoirs and fractured horizontal wells in unconventional gas reservoirs. In this model, the entire reservoir was divided into three regions (the external reservoir, the internal reservoir, and the hydraulic fracture). The flow of gas between the three regions was assumed to be linear in this model. The results indicated that the most effective way to increase the productivity of unconventional gas reservoirs is to increase the fracture density. The effects of fracture permeability on productivity are not obvious. Then, a new trilinear flow model [13] was proposed based on Ozkan’s model. Similarly, the model divided the entire reservoir into three parts: the fractured zone, the unreformed zone, and the artificial fracture. The results studied by this model points out that due to the multiple factors such as rock properties and ground stress, the hydraulic fractures formed during the fracturing process are not necessarily perpendicular to the horizontal wellbore but form a certain fracture network structure. In 2014, Dejam et al. [14] established a fractured horizontal well flow model for shale gas reservoirs based on the triple-porosity four-linear flow model, which considers the effects of gas adsorption and desorption effect of gas in the matrix. The model mainly studies the quasisteady state turbulence of the matrix to the microfracture medium and the large-fracture medium. In 2017, the unsteady flow regimes of a slightly compressible fluid under the linear and radial pre-Darcy flow conditions are modeled and the corresponding highly nonlinear diffusivity equations are solved analytically by aid of a generalized Boltzmann transformation technique [15]. The influence of pre-Darcy flow on the pressure diffusion for homogeneous porous media is studied in terms of the nonlinear exponent and the threshold pressure gradient. In addition, the pressure gradient, flux, and cumulative production per unit area are compared with the classical solution of the diffusivity equation based on the Darcy flow. In 2018, a two-phase flowback model is developed with multiscale diffusion mechanisms, and it is found that the two-phase flowback and the flow consistency between the matrix and fracture network have significant influences on cumulative gas production. The multiscale diffusion mechanisms in different zones should be carefully considered in the flowback model [16]. In 2019, a model considering convective flow, gas diffusion, and surface diffusion was established to investigate gas transport mechanism, and a common practice in modeling shale gas permeability is to use Knudsen diffusion coefficient when calculating diffusive flux, but the use of Knudsen diffusion coefficient would be incorrect if the shale gas flow regime is lying in either the transition diffusion or Fick’s diffusion, in which case the diffusion coefficient must correspond to that regime [17]. Then, a series of shale gas adsorption and desorption experiments are conducted [18]. Desorption and adsorption curves are not coincident, with the former located above the latter, which suggests that adsorption hysteresis also occurs in shale gas.A multiple-fractured horizontal well (MFHW) is extensively applied to develop the shale gas reservoirs, which makes the production of unconventional reservoirs economically and practically feasible [19]. After massive fracturing, the complex fracture network will exist around the horizontal wellbore. To model the flow behavior more rigorously, the fractures within the network should be represented explicitly rather than idealized as dual-porosity media around the wellbore. A mathematical model considering the stress sensitivity of the reservoir permeability, Darcy flow, diffusion, and adsorption and desorption in shale gas reservoirs is developed, and the numerical nonlinear production decline equations are derived and obtained. Subsequently, the model is verified by a simplified model, and production decline curves for a MFHW with discrete fracture networks in shale gas reservoirs are plotted [20]. An analytical solution is developed for the shale gas productivity of a multiple-fractured horizontal well based on a diffusion model and a trilinear flow pattern. The shale gas reservoir is divided into three flow regions: hydraulic-fracture region, microfracture network or dual-porosity region, and pure-matrix region. For the pure-matrix region, a transient diffusion equation is solved based on our previous diffusivity model developed for the shale matrix. For the microfracture network region, a modified dual-porosity model is proposed wherein both the free and adsorbed gases in the shale matrix flow into the microfracture network through a pseudosteady diffusion process. A dimensionless solution is obtained for the bottom-hole pressure in the Laplace domain considering the skin effect. An analytical solution is obtained for the gas production rate in a real-time domain through a partial Taylor series simplification and Laplace inverse transform. This analytical solution is compared with the field data of the shale gas produced from a fractured horizontal well located in southwestern China, and a good agreement is observed [21]. A semianalytical model of finite-conductivity multiple-fractured horizontal wells in LC gas reservoirs is established based on the Laplace-space superposition principle and fracture discrete method. The proposed model is validated against a commercial numerical simulator. Type curves are obtained to study pressure characteristics and identify flow regimes. The effects of some parameters on type curves are discussed [22].Different models are featured in different aspects, so there are many problems at present. (1) The migration mechanisms of shale gas in the pores are not clear. In addition, the equation describing the flow of shale gas in the pores is mostly derived under the conditions of single capillary and ideal gas. (2) The effects of flow law are not considered comprehensively in the existing models, and the models do not fit well with the actual conditions; it is difficult for the model to fully describe the Darcy flow, Knudsen diffusion, adsorption-desorption, Fick diffusion, and other mechanisms. (3) The research results and methods of shale gas in microcosmic are difficult to be used in macroscopic research. The research results between microcosmic flow and macroscopic flow are not closely related.First, the shale gas multiscale flow mechanisms are presented. Then, a flow model was established that considered viscous flow in microfractures and inorganic pores, the diffusion of Knudsen in nanoscale porosity, the coexistence of slippage, and adsorption-desorption effects under infinity and closed outer boundary conditions; based on the continuous point source solution, a multiple-fractured horizontal well flow model was established and solved by MATLAB programming, and the effects of various factors were investigated.


## 2. Shale Gas Multiscale Flow Mechanisms

Unlike conventional reservoirs, shale gas reservoirs have the characteristics of low porosity, low permeability, self-generation, and self-storage of hydrocarbon organics. The spatial structure of the reservoir is a typical multiscale feature. The migration of shale gas in the reservoir is a typical multiscale flow usually including the following three stages: (1) free gas flows from pores to the wellbore; (2) after the pore pressure drops, the adsorbed gas on the pore wall surface begins to desorb into the pores and then become free gas; and (3) dissolved gas in organic matter diffuses to the pore surface [23–25].


### 2.1. Flow in Fractures and Macropores

The flow of shale gas in natural fractures and large matrix pores satisfies Darcy’s law, which is(1)v=−klμ∇p.v is the gas flow velocity in m/s, kl is the medium permeability in μm2, μ is the gas viscosity in mPa·s, and ∇ is the gradient operator, ∇=∂/∂xi→+∂/∂yj→+∂/∂zk→.During the development of shale gas reservoirs, natural fractures and matrix systems exhibit stress-sensitive effects as formation pressures decrease. The permeabilitykl is a function of formation pressure at this time, which is
(2)kl=klpl.pl is the medium system pressure in Pa.


### 2.2. Flow in Nanoscale Pores

In shale reservoirs, nanoscale pores are mainly organic pores, of which the velocity of gas molecules is not zero at the wall surface, and both the slip flow and the Knudsen effect exist in pores, as shown in Figure2.Figure 2
Migration of single component gas in nanoscale pores.There would be slippage when the fluid flowed in porous media [26]. Therefore, the following formula was proposed to calculate the apparent gas permeability of reservoir rocks:
(3)ka=km1+bp¯,where ka is the apparent permeability in μm2, km is the equivalent liquid permeability in μm2, b is the slip factor which depends on the gas type and pore throat characteristics in Pa, and p¯ is the average pressure of the flow path in Pa.There are many models of visual permeability in shale nanoscale pores, which are summarized in Table1, and the slip coefficient expressions of various scholars are obtained. It is noteworthy that these models only consider the slippage mechanism. The Knudsen diffusion terms must be added if these models are applied to nanoscale pores.Table 1
Various calculation formulas for slip factor.
Serial numberExpression ofbAuthor1b=4cλp¯rKlinkenberg [25]2b=8πRT/Mg0.5μgr2/a−1Javadpour [23]RJ/mol/K,TK,Mkg/mol,μgPa⋅s,bPa3b=μgπRTϕτMgK∞Civan [26]RJ/mol/K,TK,Mkg/mol,μgPa⋅s,kmD,bPaIn this paper, the expression of the coefficient of permeability correctionF [27] considering slippage effect and Knudsen is as follows:
(4)F=1+αKnKn1+4Kn1−b0Kn.αKn is the sparse function related to Kn; it can be expressed as follows:
(5)αKn=α02πtan−1c1Knc2.When the flow pattern of shale gas in the nanoscale pores is slippage,b0=−1, α0=64/15π, c1=4, and c2=0.4.


### 2.3. Gas Desorption on Organic Surface

A large amount of shale gas is absorbed on the surface of organic matter under actual formation conditions, and desorption occurs mainly during mining [28–32]. There are many factors affecting the desorption of shale gas, such as organic matter content, organic matter maturity, temperature, and pressure [33–35]. The relationship between various factors is complicated.In the establishment of the flow model, the instantaneous desorption model is used to characterize the shale gas. The expression of the Langmuir isotherm equation is as follows:(6)V=VmppL+p.The mass flow of desorption per unit volume of shale gas reservoir per unit of time during extraction is(7)qdes=ρgsc∂V∂t,where qdes is the shale gas mass flow desorbed from a unit volume of reservoir rock in kg/s and ρgsc is the shale gas density under standard conditions in kg/m3.Substituting Equation (6) into Equation (7),
(8)qdes=ρgscVmpLpL+p2∂p∂t.If the Langmuir equation is expressed in pseudopressure form, the above equation can be changed to(9)qdes=ρgscVmψpLψpL+ψp2∂ψp∂t.At this point, Equations (8) and (9) can be used to express the instantaneous amount of desorption gas, and then the instability flow model of shale gas wells can be studied.


### 2.4. Diffusion in Porous Kerogen (Organic Matter)

There are two forms of diffusion: steady state diffusion and nonsteady state diffusion. When the gas consistence in the porous kerogen is independent of the spatial position and only related to time, it is called pseudosteady state diffusion and can be described by Fick’s first diffusion law; when the gas concentration in organic matter is related to both time and space coordinates, it is called nonsteady state diffusion and can be described by Fick’s second diffusion law.


#### 2.4.1. Pseudosteady State Diffusion

The derivative of shale gas consistence in the matrix with respect to time is proportional to the consistence difference of gas inside and outside the organic matter, which can be expressed as(10)dCodt=DFsCEp−Co.Therefore, the diffusion gas flow from the unit volume of shale organic matter to the nanoscale pores is(11)qF=MgFgdCodt,where Co is the molar concentration of shale gas in organic matter under pseudosteady state conditions in mol/m3; CE is the molar concentration of gas at organic wall in mol/m3; p is the nanoscale pore pressure; Fs is the shape factor in 1/m2, determined by the shape of the organic matter; D is the Fick diffusion coefficient in m2/s; t is the time in s; qF is the Fick mass flow in kg/(m3·s); and Fg is the geometric factor; it is determined by the organic shapes that are listed in Table 2.Table 2
Shale organic form factors and geometric factors of different shapes.
Organic shapeSchematicFeature lengthValues ofFsValues ofFgBlockThickness,2h2π2/4h2CylinderCylindrical radius,R45.7832/R2SphereSphere radius,R6π2/R2


#### 2.4.2. Nonsteady State Diffusion

The diffusion of dissolved gas in organic matter is usually an unsteady state diffusion. If there is a stable and continuous diffusion source, the pseudosteady diffusion rule will be met. Based on theory of uniformity, this paper assumes that the organic matter is spheres of equal size with radius ofRo. The physical model is shown in Figure 3.Figure 3
Matrix part of the quad medium model for shale gas reservoirs.According to the unsteady diffusion theory, the dissolved gas consistenceCo is a function of time and spatial position. For the spherical matrix model, it is assumed that the gas consistence on the outer surface of the organic matter is in dynamic equilibrium with the free gas in the nanoscale pores. The gas consistence change in the organic matter can be described by the following mathematical model:
(12)∂Co∂t=1ro2∂∂roDro2∂Co∂ro,where Co is the volumetric consistence of shale gas in organic matter under unsteady conditions in mol/m3 and ro is the radial coordinates of spherical organic rock mass in m.Before a gas reservoir is put into production, the pressure at each point in the reservoir is the original formation pressurepi, and the gas concentration at each point in the organic matter is related to its internal pressure. Furthermore, it is assumed that the rate of change in gas concentration at the center of the organic matter is zero, and the gas concentration at the outer boundary of the organic matter is the same as the gas concentration in the nanoscale pores. Therefore, the initial boundary conditions for the nonsteady state diffusion in the organic matter satisfy the following process:
(13)Cot=0,ro=Copi,∂Cot,ro=0∂ro=0,Cot,ro=Ro=Cop,where pi is the reservoir original formation pressure in MPa, p is the nanoscale pore pressure connected to the outer surface of organic matter in MPa, and Ro is the spherical organic radius in m.According to the above model and the corresponding definite conditions, the dissolved gas concentration distribution in the organic rock mass can be obtained, and then the mass diffusion of gas per unit volume of organic matter can be determined according to(14)qF=Mg3DRo∂Co∂roro=Ro.qF is the Fick mass flow in kg/(m3·s).


## 2.1. Flow in Fractures and Macropores

The flow of shale gas in natural fractures and large matrix pores satisfies Darcy’s law, which is(1)v=−klμ∇p.v is the gas flow velocity in m/s, kl is the medium permeability in μm2, μ is the gas viscosity in mPa·s, and ∇ is the gradient operator, ∇=∂/∂xi→+∂/∂yj→+∂/∂zk→.During the development of shale gas reservoirs, natural fractures and matrix systems exhibit stress-sensitive effects as formation pressures decrease. The permeabilitykl is a function of formation pressure at this time, which is
(2)kl=klpl.pl is the medium system pressure in Pa.


## 2.2. Flow in Nanoscale Pores

In shale reservoirs, nanoscale pores are mainly organic pores, of which the velocity of gas molecules is not zero at the wall surface, and both the slip flow and the Knudsen effect exist in pores, as shown in Figure2.Figure 2
Migration of single component gas in nanoscale pores.There would be slippage when the fluid flowed in porous media [26]. Therefore, the following formula was proposed to calculate the apparent gas permeability of reservoir rocks:
(3)ka=km1+bp¯,where ka is the apparent permeability in μm2, km is the equivalent liquid permeability in μm2, b is the slip factor which depends on the gas type and pore throat characteristics in Pa, and p¯ is the average pressure of the flow path in Pa.There are many models of visual permeability in shale nanoscale pores, which are summarized in Table1, and the slip coefficient expressions of various scholars are obtained. It is noteworthy that these models only consider the slippage mechanism. The Knudsen diffusion terms must be added if these models are applied to nanoscale pores.Table 1
Various calculation formulas for slip factor.
Serial numberExpression ofbAuthor1b=4cλp¯rKlinkenberg [25]2b=8πRT/Mg0.5μgr2/a−1Javadpour [23]RJ/mol/K,TK,Mkg/mol,μgPa⋅s,bPa3b=μgπRTϕτMgK∞Civan [26]RJ/mol/K,TK,Mkg/mol,μgPa⋅s,kmD,bPaIn this paper, the expression of the coefficient of permeability correctionF [27] considering slippage effect and Knudsen is as follows:
(4)F=1+αKnKn1+4Kn1−b0Kn.αKn is the sparse function related to Kn; it can be expressed as follows:
(5)αKn=α02πtan−1c1Knc2.When the flow pattern of shale gas in the nanoscale pores is slippage,b0=−1, α0=64/15π, c1=4, and c2=0.4.


## 2.3. Gas Desorption on Organic Surface

A large amount of shale gas is absorbed on the surface of organic matter under actual formation conditions, and desorption occurs mainly during mining [28–32]. There are many factors affecting the desorption of shale gas, such as organic matter content, organic matter maturity, temperature, and pressure [33–35]. The relationship between various factors is complicated.In the establishment of the flow model, the instantaneous desorption model is used to characterize the shale gas. The expression of the Langmuir isotherm equation is as follows:(6)V=VmppL+p.The mass flow of desorption per unit volume of shale gas reservoir per unit of time during extraction is(7)qdes=ρgsc∂V∂t,where qdes is the shale gas mass flow desorbed from a unit volume of reservoir rock in kg/s and ρgsc is the shale gas density under standard conditions in kg/m3.Substituting Equation (6) into Equation (7),
(8)qdes=ρgscVmpLpL+p2∂p∂t.If the Langmuir equation is expressed in pseudopressure form, the above equation can be changed to(9)qdes=ρgscVmψpLψpL+ψp2∂ψp∂t.At this point, Equations (8) and (9) can be used to express the instantaneous amount of desorption gas, and then the instability flow model of shale gas wells can be studied.


## 2.4. Diffusion in Porous Kerogen (Organic Matter)

There are two forms of diffusion: steady state diffusion and nonsteady state diffusion. When the gas consistence in the porous kerogen is independent of the spatial position and only related to time, it is called pseudosteady state diffusion and can be described by Fick’s first diffusion law; when the gas concentration in organic matter is related to both time and space coordinates, it is called nonsteady state diffusion and can be described by Fick’s second diffusion law.


### 2.4.1. Pseudosteady State Diffusion

The derivative of shale gas consistence in the matrix with respect to time is proportional to the consistence difference of gas inside and outside the organic matter, which can be expressed as(10)dCodt=DFsCEp−Co.Therefore, the diffusion gas flow from the unit volume of shale organic matter to the nanoscale pores is(11)qF=MgFgdCodt,where Co is the molar concentration of shale gas in organic matter under pseudosteady state conditions in mol/m3; CE is the molar concentration of gas at organic wall in mol/m3; p is the nanoscale pore pressure; Fs is the shape factor in 1/m2, determined by the shape of the organic matter; D is the Fick diffusion coefficient in m2/s; t is the time in s; qF is the Fick mass flow in kg/(m3·s); and Fg is the geometric factor; it is determined by the organic shapes that are listed in Table 2.Table 2
Shale organic form factors and geometric factors of different shapes.
Organic shapeSchematicFeature lengthValues ofFsValues ofFgBlockThickness,2h2π2/4h2CylinderCylindrical radius,R45.7832/R2SphereSphere radius,R6π2/R2


### 2.4.2. Nonsteady State Diffusion

The diffusion of dissolved gas in organic matter is usually an unsteady state diffusion. If there is a stable and continuous diffusion source, the pseudosteady diffusion rule will be met. Based on theory of uniformity, this paper assumes that the organic matter is spheres of equal size with radius ofRo. The physical model is shown in Figure 3.Figure 3
Matrix part of the quad medium model for shale gas reservoirs.According to the unsteady diffusion theory, the dissolved gas consistenceCo is a function of time and spatial position. For the spherical matrix model, it is assumed that the gas consistence on the outer surface of the organic matter is in dynamic equilibrium with the free gas in the nanoscale pores. The gas consistence change in the organic matter can be described by the following mathematical model:
(12)∂Co∂t=1ro2∂∂roDro2∂Co∂ro,where Co is the volumetric consistence of shale gas in organic matter under unsteady conditions in mol/m3 and ro is the radial coordinates of spherical organic rock mass in m.Before a gas reservoir is put into production, the pressure at each point in the reservoir is the original formation pressurepi, and the gas concentration at each point in the organic matter is related to its internal pressure. Furthermore, it is assumed that the rate of change in gas concentration at the center of the organic matter is zero, and the gas concentration at the outer boundary of the organic matter is the same as the gas concentration in the nanoscale pores. Therefore, the initial boundary conditions for the nonsteady state diffusion in the organic matter satisfy the following process:
(13)Cot=0,ro=Copi,∂Cot,ro=0∂ro=0,Cot,ro=Ro=Cop,where pi is the reservoir original formation pressure in MPa, p is the nanoscale pore pressure connected to the outer surface of organic matter in MPa, and Ro is the spherical organic radius in m.According to the above model and the corresponding definite conditions, the dissolved gas concentration distribution in the organic rock mass can be obtained, and then the mass diffusion of gas per unit volume of organic matter can be determined according to(14)qF=Mg3DRo∂Co∂roro=Ro.qF is the Fick mass flow in kg/(m3·s).


## 2.4.1. Pseudosteady State Diffusion

The derivative of shale gas consistence in the matrix with respect to time is proportional to the consistence difference of gas inside and outside the organic matter, which can be expressed as(10)dCodt=DFsCEp−Co.Therefore, the diffusion gas flow from the unit volume of shale organic matter to the nanoscale pores is(11)qF=MgFgdCodt,where Co is the molar concentration of shale gas in organic matter under pseudosteady state conditions in mol/m3; CE is the molar concentration of gas at organic wall in mol/m3; p is the nanoscale pore pressure; Fs is the shape factor in 1/m2, determined by the shape of the organic matter; D is the Fick diffusion coefficient in m2/s; t is the time in s; qF is the Fick mass flow in kg/(m3·s); and Fg is the geometric factor; it is determined by the organic shapes that are listed in Table 2.Table 2
Shale organic form factors and geometric factors of different shapes.
Organic shapeSchematicFeature lengthValues ofFsValues ofFgBlockThickness,2h2π2/4h2CylinderCylindrical radius,R45.7832/R2SphereSphere radius,R6π2/R2


## 2.4.2. Nonsteady State Diffusion

The diffusion of dissolved gas in organic matter is usually an unsteady state diffusion. If there is a stable and continuous diffusion source, the pseudosteady diffusion rule will be met. Based on theory of uniformity, this paper assumes that the organic matter is spheres of equal size with radius ofRo. The physical model is shown in Figure 3.Figure 3
Matrix part of the quad medium model for shale gas reservoirs.According to the unsteady diffusion theory, the dissolved gas consistenceCo is a function of time and spatial position. For the spherical matrix model, it is assumed that the gas consistence on the outer surface of the organic matter is in dynamic equilibrium with the free gas in the nanoscale pores. The gas consistence change in the organic matter can be described by the following mathematical model:
(12)∂Co∂t=1ro2∂∂roDro2∂Co∂ro,where Co is the volumetric consistence of shale gas in organic matter under unsteady conditions in mol/m3 and ro is the radial coordinates of spherical organic rock mass in m.Before a gas reservoir is put into production, the pressure at each point in the reservoir is the original formation pressurepi, and the gas concentration at each point in the organic matter is related to its internal pressure. Furthermore, it is assumed that the rate of change in gas concentration at the center of the organic matter is zero, and the gas concentration at the outer boundary of the organic matter is the same as the gas concentration in the nanoscale pores. Therefore, the initial boundary conditions for the nonsteady state diffusion in the organic matter satisfy the following process:
(13)Cot=0,ro=Copi,∂Cot,ro=0∂ro=0,Cot,ro=Ro=Cop,where pi is the reservoir original formation pressure in MPa, p is the nanoscale pore pressure connected to the outer surface of organic matter in MPa, and Ro is the spherical organic radius in m.According to the above model and the corresponding definite conditions, the dissolved gas concentration distribution in the organic rock mass can be obtained, and then the mass diffusion of gas per unit volume of organic matter can be determined according to(14)qF=Mg3DRo∂Co∂roro=Ro.qF is the Fick mass flow in kg/(m3·s).


## 3. Establishment and Point Source Solutions for Flow Model

### 3.1. Characterization of the Microscopic Flow Model

Strictly speaking, the shale reservoirs should be divided into four parts of organic matter that existed which are microfracture systems, inorganic pores, nanoscale pores, and gas dissolved organic matter [36, 37]. A comprehensive flow model is proposed in this section that considers multiple porous media, including the viscous flow in microfractures and inorganic pores, the diffusion of Knudsen in nanoscale porosity, the coexistence of slippage, adsorption-desorption effects, and the diffusion of dissolved gases in organic matter [38, 39].


#### 3.1.1. Dissolved Gas Diffusion in Porous Kerogen

The diffusion of dissolved gas in porous kerogen can be described by the nonsteady state Fick diffusion law:(15)1ro2∂∂roro2Do∂C∂ro=∂C∂t.


#### 3.1.2. Flow in Nanoscale Pores in Kerogen

There are three sources of gas in organic pores, namely, free gas in pores, desorption gas on the surface of organic matter, and dissolved gas diffused from kerogen. The flow equation in the pores of spherical kerogen is as follows:(16)−1rk2∂∂rkrk2ρkvk+qck+qbk=ϕkfk∂ρk∂t,vk=−kkappμ∂pk∂rk,kkapp=1+α01+AKn−BKn1+4Kn1−bKnk∞,where rk is the radial radius of kerogen in nanoscale pores in m, ρk is the gas density in kerogen nanoscale pores in kg/m3, vk is the gas permeation velocity in kerogen nanoscale pores in m/s, fk is the ratio of kerogen volume to rock volume, k∞ is the inherent permeability of kerogen nanoscale pores in mD, qck is the gas mass rate desorbed from organic pore in kg/s, qbk is the mass rate of gas diffused from organic matter in kg/s, VLk is the Langmuir volume in kerogen in kg/s, and pLk is the Langmuir pressure in kerogen in MPa.
(17)qck=−ρscVLkpLkpk+pLk2∂pk∂t,qbk=−1−ϕkfk3RoρscDo∂C∂roro=Ro,rk.The correction coefficients of permeability can be taken as follows:b=1.000, α0=1.358, A=0.170, and B=0.435.


#### 3.1.3. Flow in Matrix Macropore

The flow equation of matrix macropore in spherical coordinates can be written as(18)−1rm2∂∂rmrm2ρmvm+qkm=ϕm∂ρm∂t,vm=−kmμ∂pm∂rm,qkm=−3Rmfmρmkkappμ∂pk∂rkrk=Rk,rm.


#### 3.1.4. Flow in Microfracture System

The mass conservation equation of the microfracture system is(19)−1rf2∂rf2ρfvf∂rf+qmf=ϕf∂ρf∂t,vf=−kfμ∂pf∂rf,qmf=−3Rmρmkmμ∂pm∂rmrm=Rm.


### 3.2. The Solution of Microscopic Flow Model

The general equation of flow for the corresponding microfracture system under various flow models is as follows:(20)1rD2∂∂rDrD2∂Δm¯f∂rD=fsΔm¯f,where fs is the corresponding parameter and s is the Laplace variable.The developed microscopic flow model divides shale reservoirs into four parts: the microfracture system, inorganic pores, nanoscale porosity pores, and organic matter with dissolved gas. A flow model that is a comprehensive consideration of porous media was established, including viscous flow in microfractures and inorganic pores, Knudsen diffusion in nanoscale porosity, coexistence of slippage, adsorption-desorption effects, and dissolved gas diffusion in organic matter. The solution processes of the various microscopic flow mechanism model are in the appendix.The dimensionless variables are defined in Table3.Table 3
Definitions of dimensionless parameters.
Dimensionless timetD=kftϕμiCtif+m+kLref2Dimensionless radial radius of kerogen in nanoscale poresrkD=rkRkDimensionless matrix radial radiusrmD=rmRmStorage coefficient of the kerogenωk=ϕμCtkiϕμiCtif+m+kStorage coefficient of the matrixωm=ϕμCtmiϕμiCtif+m+kStorage coefficient of the fracture systemωf=ϕμCtfiϕμiCtif+m+kDimensionless Fick diffusion coefficientDoD=ϕμiCtif+m+kRk2kk∞Ro2Transfer coefficient from the matrix to the fractureλmf=kmLref2kfRm2Transfer coefficient from kerogen to matrixλkm=kk∞Rm2kmRk2The resulting expressionfs is as follows:
(21)fs=ωfs+λm−λm2λkωks/λk+ωks+λm+ωms.So far, we have gotten the expressions of dimensionless pseudopressure of microfractures under various flow mechanisms of shale gas reservoirs and given their unified form in spherical coordinates.


### 3.3. Continuous Point Source Solution of Circular Boundary Shale Gas Reservoir

#### 3.3.1. Flow Model

In this section, the continuous point source solution at any point in an anisotropic circular gas reservoir with a closed upper and lower boundary is derived. The physical model is shown in Figure4. The main assumptions in the model are as follows: (1) the gas reservoir is homogeneous, the thickness is h, the initial pressure is pi; (2) the permeability of the gas reservoir in the horizontal and vertical directions are kfh and kfz, respectively; (3) the middle distance of the gas reservoir is the lower boundary zw, where exists a cylindrical source sink. The radius of the microelement cylinder is r and the height is ε. (4) The source and sink strength of the cylinder is qscins, and (5) the gas flow in the gas reservoir meets the Darcy flow, and the effects of capillary pressure and gravity on the flow are ignored.Figure 4
Schematic diagram of a continuous point source physical model in the gas reservoir with closed boundary on top and bottom.According to the flow theory of oil and gas reservoirs, the continuity equation of gas flow is(22)1r∂rρgvr∂r+∂ρgvz∂z=−∂ρgϕf∂t,where ρg is the density of gas in crack systems in kg/m3, vr is the radial flow velocity of gas in m/s, vz is the flow velocity of gas in the vertical direction in m/s, ϕf is the porosity of the fracture system, r/z is the radial and vertical coordinates, and t is the production time in s.The flow equations of gas flow in the transverse and longitudinal directions are as follows:(23)vr=‐kfhμg∂pf∂r,vz=‐kfzμg∂pf∂z,where kfr is the permeability of the reservoir fracture system in the horizontal direction in m2, kfz is the permeability of the reservoir fracture system in the vertical direction in m2, pf is the crack system pressure in Pa, and μg is the gas viscosity in Pa·s.According to the state equation of real gas, the density of gas can be expressed as follows:(24)ρg=pfMgZRT,where Mg is the gas molar mass in kg/mol, T is the absolute temperature of the gas in K, and Z is the gas compression factor.According to the definition of rock compressibility, the following relationship exists between the porosity of the rock fracture system and the pressure of the fracture system:(25)cf=1ϕf∂ϕf∂pf,where cf is the fracture pore compression factor.For the inner boundary conditions, it is assumed that gas enters mainly from the side of the cylinder, and as the radius of the microelement tends to 0, it can be considered as a continuous point sink, and its inner boundary expression is [38]
(26)limε→0limr→0∫zw−ε/2zw+ε/22πrBgkfhμg∂pf∂rdzwr→0=0,z−zw>ε2,qscins,z−zw≤ε2.Among them, the volume coefficient of gasBg is as follows:
(27)Bg=pscZTpfTsc.The upper and lower boundaries of the gas reservoir are all closed. Therefore, the upper and lower boundary conditions can be described as(28)∂pf∂zz=0=0,∂pf∂zz=h=0.The outer boundary conditions can be described as(29)r∂pf∂rr=re=ϑpfr=re.Introducing the following proposed pressure,(30)ψpf=2∫pipfpμgZdp.Substituting the gas density, volumetric coefficient, and pseudopressure definition into the flow equation of the fracture system,(31)1r∂∂rr∂ψf∂r+kfzkfh∂2ψf∂z2=ϕfμgcfgkfh∂ψf∂t,(32)limε→0limr→0∫zw−ε/2zw+ε/2r∂ψf∂rdzwr→0=0,z−zw>ε2,pscTπkfhTschqscins,z−zw≤ε2.For the above continuity equation, the viscosity and compressibility of the gas are all functions of pressure. In order to obtain the analytical solution of the equation, the viscosity and compressibility of the part are usually regarded as the values in the original state.(33)μg=μgi,Cg=Cgi.Several dimensionless variables are defined in Table4.Table 4
Definitions of dimensionless parameters.
Dimensionless microfracture radial distancerD=r/LrefDimensionless outer boundary radiusreD=re/LrefDimensionless timetD=kfht/ϕfμgiCgiLref2Dimensionless distance from the lower boundaryzD=z/hDimensionless height of the gas reservoirhD=h/Lrefkfh/kfzDimensionless pressureΔψf=ψfi−ψfSubstituting the above dimensionless quantities into Equations (30) and (31), the flow equation of the fracture system becomes as follows:
(34)1rD∂∂rDrD∂Δψf∂rD+1hD2∂2Δψf∂zD2=∂Δψf∂tD.The inner boundary condition becomes as follows:(35)limε→0limrD→0∫zwD−ε/2zwD+ε/2rD∂Δψf∂rDdzwDr→0=0,zD−zwD>ε2,−pscTπkfhTschqscins,zD−zwD≤ε2.The upper and lower boundary conditions become as follows:(36)∂Δψf∂zDzD=0=0,∂Δψf∂zDzD=1=0.The outer boundary condition becomes as follows:(37)rD∂Δψf∂rDrD=reD=ϑΔψfrD=reD.


#### 3.3.2. The Solution of the Flow Model

Since there are many derivation formulas for solving similar flow models (appendix), this section does not carry out long derivation but directly gives the general solution of the model; as shown in Equation (38), the solution process requires Laplace transform and Fourier cosine transform.
(38)Δψ¯f=pscTq¯scinsπkfhhTscK0ξ0rD+CnI0ξ0rD+2∑n=1∞K0ξnrD+CnI0ξnrDcosnπzwDcosnπzD,where rD=xD−xwD2+yD−ywD2, ξn=fs+nπ2/hD2n=0,1,2,3⋯, and Cn=reDξnK1ξnreD+ϑK0ξnreD/reDξnI1ξnreD−ϑI0ξnreD.


#### 3.3.3. The Continuous Point Source Solution under Different Outer Boundary Conditions

Above, we have solved the hybrid outer boundary model. However, the outer boundary of a gas reservoir has the following two situations: infinite and closed boundaries. Therefore, the continuous point source solutions were analyzed in these two cases.(1) Infinity Outer Boundary. When the outer boundary of a gas reservoir is infinite, the outer boundary equation is as follows:
(39)limrD→0∂Δψf∂rDrD=reD=0.Therefore,(40)Cn=0.Substituting Equation (40) into Equation (39), we obtain a continuous point source solution with closed top and bottom boundaries and infinite outer boundary:
(41)Δψ¯f=pscTq¯scinsπkfhhTscK0ξ0rD+2∑n=1∞K0ξnrDcosnπzwDcosnπzD.The source intensityqscins for a continuous point is a constant. Thus, the above equation can be written as follows:
(42)Δψ¯f=pscTqscinsπkfhhTscsK0ξ0rD+2∑n=1∞K0ξnrDcosnπzwDcosnπzD.(2) Closed Outer Boundary. When the outer boundary is closed, it is known from the closed boundary conditions:
(43)Cn=K1ξnreDI1ξnreD.When the outer boundary is closed, the continuous point source solution is as follows:(44)Δψ¯f=pscTq¯scinsπkfhhTscK0ξ0rD+K1ξnreDI1ξnreDI0ξ0rD+2∑n=1∞K0ξnrD+K1ξnreDI1ξnreDI0ξnrDcosnπzwDcosnπzD.


## 3.1. Characterization of the Microscopic Flow Model

Strictly speaking, the shale reservoirs should be divided into four parts of organic matter that existed which are microfracture systems, inorganic pores, nanoscale pores, and gas dissolved organic matter [36, 37]. A comprehensive flow model is proposed in this section that considers multiple porous media, including the viscous flow in microfractures and inorganic pores, the diffusion of Knudsen in nanoscale porosity, the coexistence of slippage, adsorption-desorption effects, and the diffusion of dissolved gases in organic matter [38, 39].


### 3.1.1. Dissolved Gas Diffusion in Porous Kerogen

The diffusion of dissolved gas in porous kerogen can be described by the nonsteady state Fick diffusion law:(15)1ro2∂∂roro2Do∂C∂ro=∂C∂t.


### 3.1.2. Flow in Nanoscale Pores in Kerogen

There are three sources of gas in organic pores, namely, free gas in pores, desorption gas on the surface of organic matter, and dissolved gas diffused from kerogen. The flow equation in the pores of spherical kerogen is as follows:(16)−1rk2∂∂rkrk2ρkvk+qck+qbk=ϕkfk∂ρk∂t,vk=−kkappμ∂pk∂rk,kkapp=1+α01+AKn−BKn1+4Kn1−bKnk∞,where rk is the radial radius of kerogen in nanoscale pores in m, ρk is the gas density in kerogen nanoscale pores in kg/m3, vk is the gas permeation velocity in kerogen nanoscale pores in m/s, fk is the ratio of kerogen volume to rock volume, k∞ is the inherent permeability of kerogen nanoscale pores in mD, qck is the gas mass rate desorbed from organic pore in kg/s, qbk is the mass rate of gas diffused from organic matter in kg/s, VLk is the Langmuir volume in kerogen in kg/s, and pLk is the Langmuir pressure in kerogen in MPa.
(17)qck=−ρscVLkpLkpk+pLk2∂pk∂t,qbk=−1−ϕkfk3RoρscDo∂C∂roro=Ro,rk.The correction coefficients of permeability can be taken as follows:b=1.000, α0=1.358, A=0.170, and B=0.435.


### 3.1.3. Flow in Matrix Macropore

The flow equation of matrix macropore in spherical coordinates can be written as(18)−1rm2∂∂rmrm2ρmvm+qkm=ϕm∂ρm∂t,vm=−kmμ∂pm∂rm,qkm=−3Rmfmρmkkappμ∂pk∂rkrk=Rk,rm.


### 3.1.4. Flow in Microfracture System

The mass conservation equation of the microfracture system is(19)−1rf2∂rf2ρfvf∂rf+qmf=ϕf∂ρf∂t,vf=−kfμ∂pf∂rf,qmf=−3Rmρmkmμ∂pm∂rmrm=Rm.


## 3.1.1. Dissolved Gas Diffusion in Porous Kerogen

The diffusion of dissolved gas in porous kerogen can be described by the nonsteady state Fick diffusion law:(15)1ro2∂∂roro2Do∂C∂ro=∂C∂t.


## 3.1.2. Flow in Nanoscale Pores in Kerogen

There are three sources of gas in organic pores, namely, free gas in pores, desorption gas on the surface of organic matter, and dissolved gas diffused from kerogen. The flow equation in the pores of spherical kerogen is as follows:(16)−1rk2∂∂rkrk2ρkvk+qck+qbk=ϕkfk∂ρk∂t,vk=−kkappμ∂pk∂rk,kkapp=1+α01+AKn−BKn1+4Kn1−bKnk∞,where rk is the radial radius of kerogen in nanoscale pores in m, ρk is the gas density in kerogen nanoscale pores in kg/m3, vk is the gas permeation velocity in kerogen nanoscale pores in m/s, fk is the ratio of kerogen volume to rock volume, k∞ is the inherent permeability of kerogen nanoscale pores in mD, qck is the gas mass rate desorbed from organic pore in kg/s, qbk is the mass rate of gas diffused from organic matter in kg/s, VLk is the Langmuir volume in kerogen in kg/s, and pLk is the Langmuir pressure in kerogen in MPa.
(17)qck=−ρscVLkpLkpk+pLk2∂pk∂t,qbk=−1−ϕkfk3RoρscDo∂C∂roro=Ro,rk.The correction coefficients of permeability can be taken as follows:b=1.000, α0=1.358, A=0.170, and B=0.435.


## 3.1.3. Flow in Matrix Macropore

The flow equation of matrix macropore in spherical coordinates can be written as(18)−1rm2∂∂rmrm2ρmvm+qkm=ϕm∂ρm∂t,vm=−kmμ∂pm∂rm,qkm=−3Rmfmρmkkappμ∂pk∂rkrk=Rk,rm.


## 3.1.4. Flow in Microfracture System

The mass conservation equation of the microfracture system is(19)−1rf2∂rf2ρfvf∂rf+qmf=ϕf∂ρf∂t,vf=−kfμ∂pf∂rf,qmf=−3Rmρmkmμ∂pm∂rmrm=Rm.


## 3.2. The Solution of Microscopic Flow Model

The general equation of flow for the corresponding microfracture system under various flow models is as follows:(20)1rD2∂∂rDrD2∂Δm¯f∂rD=fsΔm¯f,where fs is the corresponding parameter and s is the Laplace variable.The developed microscopic flow model divides shale reservoirs into four parts: the microfracture system, inorganic pores, nanoscale porosity pores, and organic matter with dissolved gas. A flow model that is a comprehensive consideration of porous media was established, including viscous flow in microfractures and inorganic pores, Knudsen diffusion in nanoscale porosity, coexistence of slippage, adsorption-desorption effects, and dissolved gas diffusion in organic matter. The solution processes of the various microscopic flow mechanism model are in the appendix.The dimensionless variables are defined in Table3.Table 3
Definitions of dimensionless parameters.
Dimensionless timetD=kftϕμiCtif+m+kLref2Dimensionless radial radius of kerogen in nanoscale poresrkD=rkRkDimensionless matrix radial radiusrmD=rmRmStorage coefficient of the kerogenωk=ϕμCtkiϕμiCtif+m+kStorage coefficient of the matrixωm=ϕμCtmiϕμiCtif+m+kStorage coefficient of the fracture systemωf=ϕμCtfiϕμiCtif+m+kDimensionless Fick diffusion coefficientDoD=ϕμiCtif+m+kRk2kk∞Ro2Transfer coefficient from the matrix to the fractureλmf=kmLref2kfRm2Transfer coefficient from kerogen to matrixλkm=kk∞Rm2kmRk2The resulting expressionfs is as follows:
(21)fs=ωfs+λm−λm2λkωks/λk+ωks+λm+ωms.So far, we have gotten the expressions of dimensionless pseudopressure of microfractures under various flow mechanisms of shale gas reservoirs and given their unified form in spherical coordinates.


## 3.3. Continuous Point Source Solution of Circular Boundary Shale Gas Reservoir

### 3.3.1. Flow Model

In this section, the continuous point source solution at any point in an anisotropic circular gas reservoir with a closed upper and lower boundary is derived. The physical model is shown in Figure4. The main assumptions in the model are as follows: (1) the gas reservoir is homogeneous, the thickness is h, the initial pressure is pi; (2) the permeability of the gas reservoir in the horizontal and vertical directions are kfh and kfz, respectively; (3) the middle distance of the gas reservoir is the lower boundary zw, where exists a cylindrical source sink. The radius of the microelement cylinder is r and the height is ε. (4) The source and sink strength of the cylinder is qscins, and (5) the gas flow in the gas reservoir meets the Darcy flow, and the effects of capillary pressure and gravity on the flow are ignored.Figure 4
Schematic diagram of a continuous point source physical model in the gas reservoir with closed boundary on top and bottom.According to the flow theory of oil and gas reservoirs, the continuity equation of gas flow is(22)1r∂rρgvr∂r+∂ρgvz∂z=−∂ρgϕf∂t,where ρg is the density of gas in crack systems in kg/m3, vr is the radial flow velocity of gas in m/s, vz is the flow velocity of gas in the vertical direction in m/s, ϕf is the porosity of the fracture system, r/z is the radial and vertical coordinates, and t is the production time in s.The flow equations of gas flow in the transverse and longitudinal directions are as follows:(23)vr=‐kfhμg∂pf∂r,vz=‐kfzμg∂pf∂z,where kfr is the permeability of the reservoir fracture system in the horizontal direction in m2, kfz is the permeability of the reservoir fracture system in the vertical direction in m2, pf is the crack system pressure in Pa, and μg is the gas viscosity in Pa·s.According to the state equation of real gas, the density of gas can be expressed as follows:(24)ρg=pfMgZRT,where Mg is the gas molar mass in kg/mol, T is the absolute temperature of the gas in K, and Z is the gas compression factor.According to the definition of rock compressibility, the following relationship exists between the porosity of the rock fracture system and the pressure of the fracture system:(25)cf=1ϕf∂ϕf∂pf,where cf is the fracture pore compression factor.For the inner boundary conditions, it is assumed that gas enters mainly from the side of the cylinder, and as the radius of the microelement tends to 0, it can be considered as a continuous point sink, and its inner boundary expression is [38]
(26)limε→0limr→0∫zw−ε/2zw+ε/22πrBgkfhμg∂pf∂rdzwr→0=0,z−zw>ε2,qscins,z−zw≤ε2.Among them, the volume coefficient of gasBg is as follows:
(27)Bg=pscZTpfTsc.The upper and lower boundaries of the gas reservoir are all closed. Therefore, the upper and lower boundary conditions can be described as(28)∂pf∂zz=0=0,∂pf∂zz=h=0.The outer boundary conditions can be described as(29)r∂pf∂rr=re=ϑpfr=re.Introducing the following proposed pressure,(30)ψpf=2∫pipfpμgZdp.Substituting the gas density, volumetric coefficient, and pseudopressure definition into the flow equation of the fracture system,(31)1r∂∂rr∂ψf∂r+kfzkfh∂2ψf∂z2=ϕfμgcfgkfh∂ψf∂t,(32)limε→0limr→0∫zw−ε/2zw+ε/2r∂ψf∂rdzwr→0=0,z−zw>ε2,pscTπkfhTschqscins,z−zw≤ε2.For the above continuity equation, the viscosity and compressibility of the gas are all functions of pressure. In order to obtain the analytical solution of the equation, the viscosity and compressibility of the part are usually regarded as the values in the original state.(33)μg=μgi,Cg=Cgi.Several dimensionless variables are defined in Table4.Table 4
Definitions of dimensionless parameters.
Dimensionless microfracture radial distancerD=r/LrefDimensionless outer boundary radiusreD=re/LrefDimensionless timetD=kfht/ϕfμgiCgiLref2Dimensionless distance from the lower boundaryzD=z/hDimensionless height of the gas reservoirhD=h/Lrefkfh/kfzDimensionless pressureΔψf=ψfi−ψfSubstituting the above dimensionless quantities into Equations (30) and (31), the flow equation of the fracture system becomes as follows:
(34)1rD∂∂rDrD∂Δψf∂rD+1hD2∂2Δψf∂zD2=∂Δψf∂tD.The inner boundary condition becomes as follows:(35)limε→0limrD→0∫zwD−ε/2zwD+ε/2rD∂Δψf∂rDdzwDr→0=0,zD−zwD>ε2,−pscTπkfhTschqscins,zD−zwD≤ε2.The upper and lower boundary conditions become as follows:(36)∂Δψf∂zDzD=0=0,∂Δψf∂zDzD=1=0.The outer boundary condition becomes as follows:(37)rD∂Δψf∂rDrD=reD=ϑΔψfrD=reD.


### 3.3.2. The Solution of the Flow Model

Since there are many derivation formulas for solving similar flow models (appendix), this section does not carry out long derivation but directly gives the general solution of the model; as shown in Equation (38), the solution process requires Laplace transform and Fourier cosine transform.
(38)Δψ¯f=pscTq¯scinsπkfhhTscK0ξ0rD+CnI0ξ0rD+2∑n=1∞K0ξnrD+CnI0ξnrDcosnπzwDcosnπzD,where rD=xD−xwD2+yD−ywD2, ξn=fs+nπ2/hD2n=0,1,2,3⋯, and Cn=reDξnK1ξnreD+ϑK0ξnreD/reDξnI1ξnreD−ϑI0ξnreD.


### 3.3.3. The Continuous Point Source Solution under Different Outer Boundary Conditions

Above, we have solved the hybrid outer boundary model. However, the outer boundary of a gas reservoir has the following two situations: infinite and closed boundaries. Therefore, the continuous point source solutions were analyzed in these two cases.(1) Infinity Outer Boundary. When the outer boundary of a gas reservoir is infinite, the outer boundary equation is as follows:
(39)limrD→0∂Δψf∂rDrD=reD=0.Therefore,(40)Cn=0.Substituting Equation (40) into Equation (39), we obtain a continuous point source solution with closed top and bottom boundaries and infinite outer boundary:
(41)Δψ¯f=pscTq¯scinsπkfhhTscK0ξ0rD+2∑n=1∞K0ξnrDcosnπzwDcosnπzD.The source intensityqscins for a continuous point is a constant. Thus, the above equation can be written as follows:
(42)Δψ¯f=pscTqscinsπkfhhTscsK0ξ0rD+2∑n=1∞K0ξnrDcosnπzwDcosnπzD.(2) Closed Outer Boundary. When the outer boundary is closed, it is known from the closed boundary conditions:
(43)Cn=K1ξnreDI1ξnreD.When the outer boundary is closed, the continuous point source solution is as follows:(44)Δψ¯f=pscTq¯scinsπkfhhTscK0ξ0rD+K1ξnreDI1ξnreDI0ξ0rD+2∑n=1∞K0ξnrD+K1ξnreDI1ξnreDI0ξnrDcosnπzwDcosnπzD.


## 3.3.1. Flow Model

In this section, the continuous point source solution at any point in an anisotropic circular gas reservoir with a closed upper and lower boundary is derived. The physical model is shown in Figure4. The main assumptions in the model are as follows: (1) the gas reservoir is homogeneous, the thickness is h, the initial pressure is pi; (2) the permeability of the gas reservoir in the horizontal and vertical directions are kfh and kfz, respectively; (3) the middle distance of the gas reservoir is the lower boundary zw, where exists a cylindrical source sink. The radius of the microelement cylinder is r and the height is ε. (4) The source and sink strength of the cylinder is qscins, and (5) the gas flow in the gas reservoir meets the Darcy flow, and the effects of capillary pressure and gravity on the flow are ignored.Figure 4
Schematic diagram of a continuous point source physical model in the gas reservoir with closed boundary on top and bottom.According to the flow theory of oil and gas reservoirs, the continuity equation of gas flow is(22)1r∂rρgvr∂r+∂ρgvz∂z=−∂ρgϕf∂t,where ρg is the density of gas in crack systems in kg/m3, vr is the radial flow velocity of gas in m/s, vz is the flow velocity of gas in the vertical direction in m/s, ϕf is the porosity of the fracture system, r/z is the radial and vertical coordinates, and t is the production time in s.The flow equations of gas flow in the transverse and longitudinal directions are as follows:(23)vr=‐kfhμg∂pf∂r,vz=‐kfzμg∂pf∂z,where kfr is the permeability of the reservoir fracture system in the horizontal direction in m2, kfz is the permeability of the reservoir fracture system in the vertical direction in m2, pf is the crack system pressure in Pa, and μg is the gas viscosity in Pa·s.According to the state equation of real gas, the density of gas can be expressed as follows:(24)ρg=pfMgZRT,where Mg is the gas molar mass in kg/mol, T is the absolute temperature of the gas in K, and Z is the gas compression factor.According to the definition of rock compressibility, the following relationship exists between the porosity of the rock fracture system and the pressure of the fracture system:(25)cf=1ϕf∂ϕf∂pf,where cf is the fracture pore compression factor.For the inner boundary conditions, it is assumed that gas enters mainly from the side of the cylinder, and as the radius of the microelement tends to 0, it can be considered as a continuous point sink, and its inner boundary expression is [38]
(26)limε→0limr→0∫zw−ε/2zw+ε/22πrBgkfhμg∂pf∂rdzwr→0=0,z−zw>ε2,qscins,z−zw≤ε2.Among them, the volume coefficient of gasBg is as follows:
(27)Bg=pscZTpfTsc.The upper and lower boundaries of the gas reservoir are all closed. Therefore, the upper and lower boundary conditions can be described as(28)∂pf∂zz=0=0,∂pf∂zz=h=0.The outer boundary conditions can be described as(29)r∂pf∂rr=re=ϑpfr=re.Introducing the following proposed pressure,(30)ψpf=2∫pipfpμgZdp.Substituting the gas density, volumetric coefficient, and pseudopressure definition into the flow equation of the fracture system,(31)1r∂∂rr∂ψf∂r+kfzkfh∂2ψf∂z2=ϕfμgcfgkfh∂ψf∂t,(32)limε→0limr→0∫zw−ε/2zw+ε/2r∂ψf∂rdzwr→0=0,z−zw>ε2,pscTπkfhTschqscins,z−zw≤ε2.For the above continuity equation, the viscosity and compressibility of the gas are all functions of pressure. In order to obtain the analytical solution of the equation, the viscosity and compressibility of the part are usually regarded as the values in the original state.(33)μg=μgi,Cg=Cgi.Several dimensionless variables are defined in Table4.Table 4
Definitions of dimensionless parameters.
Dimensionless microfracture radial distancerD=r/LrefDimensionless outer boundary radiusreD=re/LrefDimensionless timetD=kfht/ϕfμgiCgiLref2Dimensionless distance from the lower boundaryzD=z/hDimensionless height of the gas reservoirhD=h/Lrefkfh/kfzDimensionless pressureΔψf=ψfi−ψfSubstituting the above dimensionless quantities into Equations (30) and (31), the flow equation of the fracture system becomes as follows:
(34)1rD∂∂rDrD∂Δψf∂rD+1hD2∂2Δψf∂zD2=∂Δψf∂tD.The inner boundary condition becomes as follows:(35)limε→0limrD→0∫zwD−ε/2zwD+ε/2rD∂Δψf∂rDdzwDr→0=0,zD−zwD>ε2,−pscTπkfhTschqscins,zD−zwD≤ε2.The upper and lower boundary conditions become as follows:(36)∂Δψf∂zDzD=0=0,∂Δψf∂zDzD=1=0.The outer boundary condition becomes as follows:(37)rD∂Δψf∂rDrD=reD=ϑΔψfrD=reD.


## 3.3.2. The Solution of the Flow Model

Since there are many derivation formulas for solving similar flow models (appendix), this section does not carry out long derivation but directly gives the general solution of the model; as shown in Equation (38), the solution process requires Laplace transform and Fourier cosine transform.
(38)Δψ¯f=pscTq¯scinsπkfhhTscK0ξ0rD+CnI0ξ0rD+2∑n=1∞K0ξnrD+CnI0ξnrDcosnπzwDcosnπzD,where rD=xD−xwD2+yD−ywD2, ξn=fs+nπ2/hD2n=0,1,2,3⋯, and Cn=reDξnK1ξnreD+ϑK0ξnreD/reDξnI1ξnreD−ϑI0ξnreD.


## 3.3.3. The Continuous Point Source Solution under Different Outer Boundary Conditions

Above, we have solved the hybrid outer boundary model. However, the outer boundary of a gas reservoir has the following two situations: infinite and closed boundaries. Therefore, the continuous point source solutions were analyzed in these two cases.(1) Infinity Outer Boundary. When the outer boundary of a gas reservoir is infinite, the outer boundary equation is as follows:
(39)limrD→0∂Δψf∂rDrD=reD=0.Therefore,(40)Cn=0.Substituting Equation (40) into Equation (39), we obtain a continuous point source solution with closed top and bottom boundaries and infinite outer boundary:
(41)Δψ¯f=pscTq¯scinsπkfhhTscK0ξ0rD+2∑n=1∞K0ξnrDcosnπzwDcosnπzD.The source intensityqscins for a continuous point is a constant. Thus, the above equation can be written as follows:
(42)Δψ¯f=pscTqscinsπkfhhTscsK0ξ0rD+2∑n=1∞K0ξnrDcosnπzwDcosnπzD.(2) Closed Outer Boundary. When the outer boundary is closed, it is known from the closed boundary conditions:
(43)Cn=K1ξnreDI1ξnreD.When the outer boundary is closed, the continuous point source solution is as follows:(44)Δψ¯f=pscTq¯scinsπkfhhTscK0ξ0rD+K1ξnreDI1ξnreDI0ξ0rD+2∑n=1∞K0ξnrD+K1ξnreDI1ξnreDI0ξnrDcosnπzwDcosnπzD.


## 4. Study on Unsteady Flow Rules of Circular Boundary Multiple-Fractured Horizontal Wells

### 4.1. The Establishment of a Physical Model

A model diagram of a multiple-fractured horizontal well in a circular shale gas reservoir is shown in Figure5. The following assumptions are given:
(1)
The gas reservoir is homogeneous, horizontal, and equal in thickness. The upper and lower boundaries are closed. The lateral direction is infinite or closed and the reservoir thickness ish(2)
The horizontal well is located in the center of the gas reservoir and its length isL(3)
The pressure drop caused by fluid flowing in the wellbore is ignored, and the hydraulic fractures supply fluid to the wellbore, regardless of the flow of the reservoir fluid directly to the wellbore(4)
The number of hydraulic fractures isN, and it is a fully cracked fracture with infinite flow conductivity. The fracture is regarded as a 2D surface source rather than a 3D volume source(5)
Fracture cracks are symmetrical bifurcations distributed at equal intervals along the horizontal wellbore(6)
The flow process is an isothermal flow, regardless of the effects of fluid gravity and capillary forcesFigure 5
Schematic diagram of multiple-fractured horizontal well model for shale gas reservoir.


### 4.2. Mathematical Model Establishment and Solution

As shown in Figure5, the direction of the wellbore along the horizontal well is the y-axis, the vertical upward direction is the z-axis direction, and the direction along the fracture surface is the x-axis direction. For multistage fracturing horizontal wells, the flow process of fluids in the reservoir is quite complicated. Not only is the mutual interference between the seams and the fractures, but the flow on the fracture surface is also uneven. Therefore, the solution method for continuing to use the fracturing vertical wells does not work. We can first divide the crack into several units, and it is assumed that the flow on each section is evenly distributed, but the flow in the cracks on different units is different. By solving the pressure response expressions at each fracture cell and combining the bottom hole pressure with the normalized production conditions, the expressions of the bottom hole pressures in multistage fracturing horizontal wells in shale gas reservoirs are solved.In this study, the fractures are infinitely inducing flow cracks that are completely fracturing. According to the continuous point source solution of the circular reservoir obtained, it is integrated along the longitudinal direction that we can obtain a continuous linear source solution. For a continuous line source with an intensity ofqscL, there is an infinite atmospheric boundary:
(45)Δψ¯fL=pscTqscinsπkfhhTscs∫0hK0ξ0rD+2∑n=1∞K0ξnrDcosnπzwDcosnπzDdzw,where rD=xD−xwD2+yD−ywD2 and ξn=fs+nπ2/hD2n=0,1,2,3⋯.The continuous line source intensity and continuous point source intensity have the following relationship:(46)qscL=qscinsh.The definition of the dimensionless distancezwD is as follows:
(47)zwD=zwh.Substituting Equations (46) and (47) into Equation (45), we can obtain a continuous linear source solution for an infinite boundary shale gas reservoir with a circular outer boundary:
(48)Δψ¯fL=pscTqscLπkfhhTscsK0ξ0rD.In the same way, the solution of a continuous line source can be obtained when the gas reservoir boundary is closed:(49)Δψ¯fL=pscTqscLπkfhhTscsK0ξ0rD+K1ξ0reDI1ξ0reDI0ξ0rD.The predecessors have basically formed a set of theoretical systems for the crack dispersion and superposition principle to solve the pressure response. The following is a brief introduction. In the hypothetical part of the model in this section, the number of fractures isN, the length of horizontal wells is L, and each fracture is equally divided into 2M units as shown in Figure 6.Figure 6
Crack discrete grid.From the discrete fracture schematic, the node coordinates of each segment in the fracture can be obtained. The crack node coordinates in the negative direction of thex-axis can be obtained by the following equation:
(50)xi=−M−i−k∗2M+12ΔLfLi,1≤i−k∗2M≤M,yi=yk+1.The coordinates of the crack node in the positive direction of thex-axis can be obtained by the following equation:
(51)xi=i−k∗2M−M−12ΔLfLi,M+1≤i−k∗2M≤2M,yi=yk+1,where xi is the abscissa of the i-th crack cell in m, yi is the ordinate of the i-th crack cell, and ΔLfLi is the section i-th crack length.For the crack unit length, since each crack is divided into2M parts, the length of the i-th crack cell is
(52)ΔLfLi=LfLiM.As mentioned earlier, it has been assumed that the flow distribution of each fracture cell is uniform, so that the pseudopressure drop formed at any point (x,y) in the formation at the i-th crack cell is
(53)Δψ¯fi=∫ΓΔψ¯fLdl.Since the cracks are all parallel to thex-axis, the above curve integrals can be directly converted into coordinate integrals. When the outer boundary of the reservoir is infinite, Equation (49) is substituted into Equation (53), where the coordinates of the i-th node are substituted into, we can get the following equation:
(54)Δψ¯fi=∫xi−ΔLfi/2xi+ΔLfi/2pscTqscLiπkfhhTscsK0ξ0xD−xwDi2+yD−ywDi2dxw.Since the flow is uniformly distributed on the discrete units of the fracture, assuming that thei-th discrete unit is qsci, the following relationship can be satisfied between the line flow and the unit flow:
(55)qscL=qscLiΔLfi.Introducing the following dimensionless pressure and dimensionless production:(56)ψfD=πkfhhTscpscTqscΔψf,(57)qDi=qsciqsc.Substituting Equations (56) and (57) into Equation (54),
(58)Δψ¯fDi=qDisΔLfDi∫−ΔLfDi/2ΔLfDi/2K0ξ0xD−xwDi−α2+yD−ywDi2dα,where ΔLfDi is the dimensionless length of the i-th crack unit.According to the superposition principle of the potential, the total pressure response generated at each position of the crack in theN cracks at any position in the reservoir is
(59)ψ¯DxD,yD=∑i=1Nψ¯fDixD,yD.According to the previous assumption,(60)ψ¯wD=ψ¯DxDj,yDj=∑i=1Nψ¯fDixDj,yDj,where xDj,yDj is the intersection point of the j-th section of a fracture and the wellbore.When a horizontal well is produced at a given production volume, the production of the well is numerically equal to the sum of the flow of all cracks. Therefore, we have the following equation:(61)qsc=∑i=1N∗2MqscLiΔLfi.Substituting Equation (61) into Equation. (57), there are the following normalization conditions:
(62)∑i=1N∗2MqDi=1.In Figure6, we separate the cracks into N×2M units and solve the pseudopressure at each unit, so we can have N×2M equations, which, plus the normalized equation for production, total N×2M+1 equations, while the unknowns include the flow rate and bottom hole pressure of each fracture cell, which is also N×2M+1. The equations are closed and there is a unique solution. For ease of understanding, this set of equations is written in the following matrix form:
(63)A1,1⋯A1,n⋯A1,2M∗N−1⋯⋯⋯⋯⋯−1An,1⋯An,n⋯An,2M∗N⋯⋯⋯⋯⋯⋯−1A2M∗N,1⋯A2M∗N,n⋯A2M∗N,2M∗N−11⋯1⋯10qD1qD2⋯⋯qD,2M∗Nψ¯wD=00⋯⋯01.The expression ofAi,j is as follows:
(64)Ai,j=1sΔLfDj∫−ΔLfDj/2ΔLfDj/2K0ξ0xDi−xDj−α2+yDi−yDj2dα.The above formula is the fracturing horizontal well production and pressure calculation formula under the condition of the infinite outer boundary. If the boundary of the gas reservoir is closed, the formula can be simplified to Equation (64):
(65)Ai,j=1sΔLfDj∫−ΔLfDj/2ΔLfDj/2K0ξ0rD+K1ξ0reDI1ξ0reDI0ξ0rDdα.


### 4.3. Gas Well Bottom Pressure Dynamic Curve Analysis

So far, we have obtained the unstable pressure of multiple-fractured horizontal wells in circular reservoirs considering a variety of complex flow mechanisms. According to the expression, the computer programming and Stehfest numerical inversion techniques [38] were used to obtain pressure and production dynamic curves of shale gas reservoir fractured.Figure7 shows the typical pressure response curves of multistage fracturing horizontal wells in shale reservoirs with different outer boundary conditions. Multiple concaves appear in the dynamic curve, reflecting the effects of various flow mechanisms on the pressure response. As the pressure wave propagates in the formation, different pore media rich in shale gas begin to flow, causing the pressure of the porous media to fall.Figure 7
Pressure response curves of multiple-fractured horizontal wells in shale reservoirs with different outer boundary conditions under double logarithmic coordinates.According to Figure7, we can divide the flow of horizontal wells into the following stages:Stage I: the wellbore storage and subsequent transitional flow stages. In the wellbore storage phase, the pressure and pressure derivative curves coincide with each other with a slope of 1Stage II: the early linear flow stage perpendicular to the fracturing fracture. After the completion of the well storage phase, the free gas in the reservoir flowed first in the direction perpendicular to the fracture surface, and there was no mutual interference between the fractures. In this flow stage, the slope of the pressure derivative curve is 1/2Stage III: early radial flow stage. The pressure derivative curve shows a horizontal line with a slope of 1/2NStage IV: the turbulence stage. This stage corresponds to the linear flow stage from the formation to the entire fracturing zone. The concave character appears on the pressure derivative curveStage V: the boundary control stage. The characteristics of the boundary flow are exhibited after the pressure wave spreads to the boundary. For the closed boundary, when the pressure wave spreads to the boundary, the pressure will increase faster; combined with the continuous line source solution Equation (49) under the closed boundary, the pressure and pressure derivative curves coincide and the slope is 1The parameters used in the mechanistic model are listed in Table5.Table 5
Parameters in the following models.
ParameterValueOriginal pressure,pi (MPa)25Formation thickness,h (m)40Crack number5Initial gas compression factor,Cgi (MPa-1)0.02Microfracture permeability,kf (mD)0.01Horizontal well length,L (m)1000Well storage coefficient,CD10-5Gas reservoir boundary radius,re (m)5000Temperature,T (K)335.15Half length of crack,xf (m)40Relative gas density0.65Microfracture porosity0.02Gas well production,qsc (m3/d)1.2×104Skin factor,Skin0.1The pressure curves for different slippage coefficients under double logarithmic coordinates are shown in Figure8. Making the slippage coefficient b equals 1, 8, and 16, it can be seen that as the slippage coefficient b increases and the turbulence trough appears earlier and earlier. According to the knowledge of flow mechanics, when the shale gas desorbed from the organic pore, slippage appears, so the slippage coefficient b affects the turbulence stage the most. The greater the slippage coefficient b is, the more obvious is the slippage effect, and the larger the apparent permeability of the matrix is, so that the turbulence of matrix pores to the fracture system is more serious, and eventually, the time of turbulence appears earlier. Naturally, a larger slippage coefficient will increase the production of shale gas wells.Figure 8
Pressure dynamic curves of multistage fracturing horizontal wells with different slip coefficients (infinity border).Figure9 shows the influence of different Knudsen diffusion coefficients on the pressure response of multiple-fractured horizontal wells. It is not difficult to see from the figure that the influence of Knudsen diffusion coefficient DK (m2/s) on the pressure dynamics and the gas slippage effect on the well test curve of gas wells are similar. With the increase of DK (m2/s), the turbulence troughs appear earlier and the pressure and pressure derivative curves are getting lower and lower, so the production of gas wells will also increase.Figure 9
Pressure dynamic curves of multiple-fractured horizontal wells with different Knudsen diffusion coefficients under double logarithmic coordinates (infinity border).


## 4.1. The Establishment of a Physical Model

A model diagram of a multiple-fractured horizontal well in a circular shale gas reservoir is shown in Figure5. The following assumptions are given:
(1)
The gas reservoir is homogeneous, horizontal, and equal in thickness. The upper and lower boundaries are closed. The lateral direction is infinite or closed and the reservoir thickness ish(2)
The horizontal well is located in the center of the gas reservoir and its length isL(3)
The pressure drop caused by fluid flowing in the wellbore is ignored, and the hydraulic fractures supply fluid to the wellbore, regardless of the flow of the reservoir fluid directly to the wellbore(4)
The number of hydraulic fractures isN, and it is a fully cracked fracture with infinite flow conductivity. The fracture is regarded as a 2D surface source rather than a 3D volume source(5)
Fracture cracks are symmetrical bifurcations distributed at equal intervals along the horizontal wellbore(6)
The flow process is an isothermal flow, regardless of the effects of fluid gravity and capillary forcesFigure 5
Schematic diagram of multiple-fractured horizontal well model for shale gas reservoir.


## 4.2. Mathematical Model Establishment and Solution

As shown in Figure5, the direction of the wellbore along the horizontal well is the y-axis, the vertical upward direction is the z-axis direction, and the direction along the fracture surface is the x-axis direction. For multistage fracturing horizontal wells, the flow process of fluids in the reservoir is quite complicated. Not only is the mutual interference between the seams and the fractures, but the flow on the fracture surface is also uneven. Therefore, the solution method for continuing to use the fracturing vertical wells does not work. We can first divide the crack into several units, and it is assumed that the flow on each section is evenly distributed, but the flow in the cracks on different units is different. By solving the pressure response expressions at each fracture cell and combining the bottom hole pressure with the normalized production conditions, the expressions of the bottom hole pressures in multistage fracturing horizontal wells in shale gas reservoirs are solved.In this study, the fractures are infinitely inducing flow cracks that are completely fracturing. According to the continuous point source solution of the circular reservoir obtained, it is integrated along the longitudinal direction that we can obtain a continuous linear source solution. For a continuous line source with an intensity ofqscL, there is an infinite atmospheric boundary:
(45)Δψ¯fL=pscTqscinsπkfhhTscs∫0hK0ξ0rD+2∑n=1∞K0ξnrDcosnπzwDcosnπzDdzw,where rD=xD−xwD2+yD−ywD2 and ξn=fs+nπ2/hD2n=0,1,2,3⋯.The continuous line source intensity and continuous point source intensity have the following relationship:(46)qscL=qscinsh.The definition of the dimensionless distancezwD is as follows:
(47)zwD=zwh.Substituting Equations (46) and (47) into Equation (45), we can obtain a continuous linear source solution for an infinite boundary shale gas reservoir with a circular outer boundary:
(48)Δψ¯fL=pscTqscLπkfhhTscsK0ξ0rD.In the same way, the solution of a continuous line source can be obtained when the gas reservoir boundary is closed:(49)Δψ¯fL=pscTqscLπkfhhTscsK0ξ0rD+K1ξ0reDI1ξ0reDI0ξ0rD.The predecessors have basically formed a set of theoretical systems for the crack dispersion and superposition principle to solve the pressure response. The following is a brief introduction. In the hypothetical part of the model in this section, the number of fractures isN, the length of horizontal wells is L, and each fracture is equally divided into 2M units as shown in Figure 6.Figure 6
Crack discrete grid.From the discrete fracture schematic, the node coordinates of each segment in the fracture can be obtained. The crack node coordinates in the negative direction of thex-axis can be obtained by the following equation:
(50)xi=−M−i−k∗2M+12ΔLfLi,1≤i−k∗2M≤M,yi=yk+1.The coordinates of the crack node in the positive direction of thex-axis can be obtained by the following equation:
(51)xi=i−k∗2M−M−12ΔLfLi,M+1≤i−k∗2M≤2M,yi=yk+1,where xi is the abscissa of the i-th crack cell in m, yi is the ordinate of the i-th crack cell, and ΔLfLi is the section i-th crack length.For the crack unit length, since each crack is divided into2M parts, the length of the i-th crack cell is
(52)ΔLfLi=LfLiM.As mentioned earlier, it has been assumed that the flow distribution of each fracture cell is uniform, so that the pseudopressure drop formed at any point (x,y) in the formation at the i-th crack cell is
(53)Δψ¯fi=∫ΓΔψ¯fLdl.Since the cracks are all parallel to thex-axis, the above curve integrals can be directly converted into coordinate integrals. When the outer boundary of the reservoir is infinite, Equation (49) is substituted into Equation (53), where the coordinates of the i-th node are substituted into, we can get the following equation:
(54)Δψ¯fi=∫xi−ΔLfi/2xi+ΔLfi/2pscTqscLiπkfhhTscsK0ξ0xD−xwDi2+yD−ywDi2dxw.Since the flow is uniformly distributed on the discrete units of the fracture, assuming that thei-th discrete unit is qsci, the following relationship can be satisfied between the line flow and the unit flow:
(55)qscL=qscLiΔLfi.Introducing the following dimensionless pressure and dimensionless production:(56)ψfD=πkfhhTscpscTqscΔψf,(57)qDi=qsciqsc.Substituting Equations (56) and (57) into Equation (54),
(58)Δψ¯fDi=qDisΔLfDi∫−ΔLfDi/2ΔLfDi/2K0ξ0xD−xwDi−α2+yD−ywDi2dα,where ΔLfDi is the dimensionless length of the i-th crack unit.According to the superposition principle of the potential, the total pressure response generated at each position of the crack in theN cracks at any position in the reservoir is
(59)ψ¯DxD,yD=∑i=1Nψ¯fDixD,yD.According to the previous assumption,(60)ψ¯wD=ψ¯DxDj,yDj=∑i=1Nψ¯fDixDj,yDj,where xDj,yDj is the intersection point of the j-th section of a fracture and the wellbore.When a horizontal well is produced at a given production volume, the production of the well is numerically equal to the sum of the flow of all cracks. Therefore, we have the following equation:(61)qsc=∑i=1N∗2MqscLiΔLfi.Substituting Equation (61) into Equation. (57), there are the following normalization conditions:
(62)∑i=1N∗2MqDi=1.In Figure6, we separate the cracks into N×2M units and solve the pseudopressure at each unit, so we can have N×2M equations, which, plus the normalized equation for production, total N×2M+1 equations, while the unknowns include the flow rate and bottom hole pressure of each fracture cell, which is also N×2M+1. The equations are closed and there is a unique solution. For ease of understanding, this set of equations is written in the following matrix form:
(63)A1,1⋯A1,n⋯A1,2M∗N−1⋯⋯⋯⋯⋯−1An,1⋯An,n⋯An,2M∗N⋯⋯⋯⋯⋯⋯−1A2M∗N,1⋯A2M∗N,n⋯A2M∗N,2M∗N−11⋯1⋯10qD1qD2⋯⋯qD,2M∗Nψ¯wD=00⋯⋯01.The expression ofAi,j is as follows:
(64)Ai,j=1sΔLfDj∫−ΔLfDj/2ΔLfDj/2K0ξ0xDi−xDj−α2+yDi−yDj2dα.The above formula is the fracturing horizontal well production and pressure calculation formula under the condition of the infinite outer boundary. If the boundary of the gas reservoir is closed, the formula can be simplified to Equation (64):
(65)Ai,j=1sΔLfDj∫−ΔLfDj/2ΔLfDj/2K0ξ0rD+K1ξ0reDI1ξ0reDI0ξ0rDdα.


## 4.3. Gas Well Bottom Pressure Dynamic Curve Analysis

So far, we have obtained the unstable pressure of multiple-fractured horizontal wells in circular reservoirs considering a variety of complex flow mechanisms. According to the expression, the computer programming and Stehfest numerical inversion techniques [38] were used to obtain pressure and production dynamic curves of shale gas reservoir fractured.Figure7 shows the typical pressure response curves of multistage fracturing horizontal wells in shale reservoirs with different outer boundary conditions. Multiple concaves appear in the dynamic curve, reflecting the effects of various flow mechanisms on the pressure response. As the pressure wave propagates in the formation, different pore media rich in shale gas begin to flow, causing the pressure of the porous media to fall.Figure 7
Pressure response curves of multiple-fractured horizontal wells in shale reservoirs with different outer boundary conditions under double logarithmic coordinates.According to Figure7, we can divide the flow of horizontal wells into the following stages:Stage I: the wellbore storage and subsequent transitional flow stages. In the wellbore storage phase, the pressure and pressure derivative curves coincide with each other with a slope of 1Stage II: the early linear flow stage perpendicular to the fracturing fracture. After the completion of the well storage phase, the free gas in the reservoir flowed first in the direction perpendicular to the fracture surface, and there was no mutual interference between the fractures. In this flow stage, the slope of the pressure derivative curve is 1/2Stage III: early radial flow stage. The pressure derivative curve shows a horizontal line with a slope of 1/2NStage IV: the turbulence stage. This stage corresponds to the linear flow stage from the formation to the entire fracturing zone. The concave character appears on the pressure derivative curveStage V: the boundary control stage. The characteristics of the boundary flow are exhibited after the pressure wave spreads to the boundary. For the closed boundary, when the pressure wave spreads to the boundary, the pressure will increase faster; combined with the continuous line source solution Equation (49) under the closed boundary, the pressure and pressure derivative curves coincide and the slope is 1The parameters used in the mechanistic model are listed in Table5.Table 5
Parameters in the following models.
ParameterValueOriginal pressure,pi (MPa)25Formation thickness,h (m)40Crack number5Initial gas compression factor,Cgi (MPa-1)0.02Microfracture permeability,kf (mD)0.01Horizontal well length,L (m)1000Well storage coefficient,CD10-5Gas reservoir boundary radius,re (m)5000Temperature,T (K)335.15Half length of crack,xf (m)40Relative gas density0.65Microfracture porosity0.02Gas well production,qsc (m3/d)1.2×104Skin factor,Skin0.1The pressure curves for different slippage coefficients under double logarithmic coordinates are shown in Figure8. Making the slippage coefficient b equals 1, 8, and 16, it can be seen that as the slippage coefficient b increases and the turbulence trough appears earlier and earlier. According to the knowledge of flow mechanics, when the shale gas desorbed from the organic pore, slippage appears, so the slippage coefficient b affects the turbulence stage the most. The greater the slippage coefficient b is, the more obvious is the slippage effect, and the larger the apparent permeability of the matrix is, so that the turbulence of matrix pores to the fracture system is more serious, and eventually, the time of turbulence appears earlier. Naturally, a larger slippage coefficient will increase the production of shale gas wells.Figure 8
Pressure dynamic curves of multistage fracturing horizontal wells with different slip coefficients (infinity border).Figure9 shows the influence of different Knudsen diffusion coefficients on the pressure response of multiple-fractured horizontal wells. It is not difficult to see from the figure that the influence of Knudsen diffusion coefficient DK (m2/s) on the pressure dynamics and the gas slippage effect on the well test curve of gas wells are similar. With the increase of DK (m2/s), the turbulence troughs appear earlier and the pressure and pressure derivative curves are getting lower and lower, so the production of gas wells will also increase.Figure 9
Pressure dynamic curves of multiple-fractured horizontal wells with different Knudsen diffusion coefficients under double logarithmic coordinates (infinity border).


## 5. Summary and Conclusions

(1)
A flow model for shale gas reservoirs considering flow in microfractures and inorganic pores, the diffusion of Knudsen in nanoscale porosity, the coexistence of slippage, established adsorption-desorption effects, and the continuous point source solution of the flow mechanisms model were obtained under closed and infinity outer boundaries(2)
Based on the continuous point source solution, a circular boundary multiple-fractured horizontal well flow model was established, and the gas well bottom pressure dynamic curve and the relevant sensitivity factors are analyzed with the help of MATLAB programming(3)
The results of the sensitivity analysis show that the Knudsen diffusion and slippage coefficients mainly affect the apparent permeability of the matrix pores. The more the Knudsen diffusion and slippage coefficients are, the earlier the turbulent flow occurs and the higher the gas production is



---

*Source: 7364734-2020-08-03.xml*