# Pressure Analysis for Volume Fracturing Vertical Well considering Low-Velocity Non-Darcy Flow and Stress Sensitivity

Authors: Zhongwei Wu; Chuanzhi Cui; Japan Trivedi; Ning Ai; Wenhao Tang
Journal: Geofluids
 (2019)

---

## Abstract

In general, there is stress sensitivity damage in tight reservoirs and fractures. Furthermore, the flow in tight reservoirs is the low-velocity non-Darcy flow. Currently, few researches of pressure analysis for volume fracturing vertical well are conducted simultaneously considering the low-velocity non-Darcy flow and stress sensitivity. In the paper, a novel flow model of a volume fractured vertical well is proposed and solved numerically. Firstly, the threshold pressure gradient, permeability modulus, and experimental data are, respectively, utilized to characterize the low-velocity non-Darcy flow, matrix stress sensitivity, and fracture stress sensitivity. Then, a two-region composite reservoir is established to simulate the vertical well with volume fracturing. After that, the logarithm meshing method is used to discrete the composite reservoir, and the flow model is solved by the method of finite difference and IMPES. Finally, the model verification is conducted, and the effects of the low-velocity non-Darcy flow and stress sensitivity on the pressure and pressure derivative are analyzed. The six flow regimes are identified by the dimensionless pressure and pressure derivative curve. They are, respectively, the fracture linear flow regime, early transition flow regime, radial flow regime, crossflow regime, advanced transition flow regime, and boundary controlling flow regime. The stress sensitivity and threshold pressure gradient have a great effect on the dimensionless pressure and pressure derivative. With the increase of reservoir stress sensitivity, the pressure and pressure derivative are upward at the advanced transition flow and boundary controlling regimes. However, the pressure and pressure derivative are downward at the advanced transition flow and boundary controlling regimes when the fracture sensitivity increases. An increase in the threshold pressure gradient results in a high dimensionless pressure and pressure derivative. This work reveals the effects of low-velocity non-Darcy flow and stress sensitivity on pressure and provides a more accurate reference for reservoir engineers in pressure analysis when developing a tight reservoir by using the volume fracturing vertical well.


---

## Body

## 1. Introduction

Volume fracturing breaks up effective permeable reservoirs to form a fracture network, maximizes the contact area between fracture sides and reservoir matrix, minimizes the fluid seepage distance of oil and gas from the matrix in any directions to fractures, and greatly improves the entire reservoir permeability [1]. Hydraulic fracture technology is an effective method to exploit unconventional reservoirs and is reported in many works [2–5], especially for volume fracturing [6–8]. Many microfractures are developed in the tight reservoir after volume fracturing [9, 10]. During depletion, reservoir pressure changes could lead to stress variations which further alter the fracture apertures, pore, and throats, finally changing microfracture permeability and reservoir permeability [11]. In other words, the fracture and reservoir all are stress sensitivity medium [12–14]. In addition, due to the tiny porous and ultralow permeability of tight reservoirs, the flow in tight reservoirs obeys the law of low-velocity non-Darcy flow instead of Darcy flow [15], which is the fact that the flow velocity in a low pressure gradient regime is lower than what is estimated from Darcy’s law [16].Currently, there are many researches of the volume fracturing. Wang et al. [17] and Lu et al. [18] proposed a semianalytical model to analyze the productivity and pressure of a volume fractured vertical well in tight reservoirs. In their research, the reservoir consists of two regions, one of which is the main fracture. Another is a circular stimulated reservoir volume (SRV), which is idealized to a Warren-Root model [19]. The unstimulated reservoir volume (un-SRV) is not contributed to the flow. That is different from the actual situation, and the stress sensitivity and low-velocity non-Darcy flow are ignored. Zhu et al. [20] gave an analytical pressure solution for the volume fractured vertical well with a rectangular SRV. In their research, the reservoir is composed of multiple linear flow regions, which include SRV and un-SRV. Unfortunately, they do not consider the effect of stress sensitivity and low-velocity non-Darcy flow. Su et al. [21] proposed a fractal analytical model of volume fractured vertical wells. In their research, the whole reservoir is divided into two regions, one of which is an SRV and another is an un-SRV. The flow obeys the Darcy, and the permeability is not a function of effective pressure, which deviates the actual situation in the tight reservoir. Zhang et al. [22] proposed a semianalytical model for vertically fractured wells with stimulated reservoir volumes (SRV). The fractal porosity and permeability are employed to describe the heterogeneous distribution of porous media in the SRV. But the stress sensitivity and low-velocity non-Darcy flow are ignored. The researches mentioned above are all about volume fracturing, but they do not take the effect of the stress sensitivity and low-velocity non-Darcy flow into consideration.Guo et al. [23] proposed a transient mathematical model for horizontal wells with consideration of the threshold pressure and nonlinear flow. In their research, the nonlinear flow is high-velocity nonlinear flow, which is different from that caused by low-velocity mentioned in this paper. And their work is just applied to the multistage fractured horizontal well. Wang et al. [24] conducted research of the transient pressure of the volume fractured vertical well. In his research, the two-region composite model is built. The inner region represents a circle SRV and is idealized as a dual-porosity medium, while the outer zone represents an un-SRV and is idealized as a single porosity reservoir. However, the low-velocity non-Darcy flow and the stress sensitivity of fracture medium are ignored. Ji et al. [25] present a new analytical model for MFHW in tight oil reservoirs where the reservoir is subdivided into five continuous linear flow regions. The effects of non-Darcy flow and stress sensitivity are both considered. However, his research is just applied to the multistage fractured horizontal well. By considering the effect of the stress sensitivity and low-velocity non-Darcy flow, Wu et al. [26] proposed a pressure analysis model for a multistage fractured horizontal well in the tight reservoir. In his research, the tight reservoir is divided into 4 regions, which are region I to region IV. The medium of each region is stress sensitivity, and the flow in the unstimulated reservoir region obeys the low-velocity non-Darcy law. But his research is just for the multistage fractured horizontal well instead of the volume fracturing vertical well. Wu et al. [27] conduct another research of transient pressure. In this work, the fracture is discrete into many subsections. The superposition principle is utilized to build the semianalytical pressure analysis model of a multistage fracturing horizontal well considering the stress sensitivity of fracture medium. However, the reservoir is an ordinary medium instead of stress sensitivity medium in their research, and the low-velocity non-Darcy flow is ignored.As stated above, few works simultaneously considering the effect of the low-velocity non-Darcy flow and stress sensitivity are conducted for the volume fracturing well. This paper fills this gap. To this end, the researches of the flow model of the volume fracturing well and its verification are conducted in Section2. After that, the result and discussion are presented in Section 3. In Section 4, conclusions are presented.


## 2. The Flow Model and Verification

### 2.1. Physical Model of Volume Fracturing Vertical Well

According to microseismic monitoring mapping [24] and the result of indoor core fracturing experiment [28, 29], we can know that the reservoir after volume fracturing can be divided into two regions, one is the SRV and another is the un-SRV. The SRV is circular, and many fractures are developed in the SRV (Figure 1).Figure 1
The reservoir morphology after vertically fractured.
(a)
Microseismic monitoring mapping [24]
(b)
Fracture morphology of unconventional reservoir core [28, 29]As the statement above, the physical model of the vertical well with volume fracturing can be obtained (Figure2). As seen from Figure 2, the tight reservoir is composed of two regions, which are, respectively, the SRV and un-SRV. The un-SRV is treated as a single porosity medium. And the SRV is idealized by double porosity medium [19]. The un-SRV is stress sensitivity medium, and the flow in un-SRV obeys the low-velocity non-Darcy flow. In SRV, the microfracture is developed, and all microfractures are crossed through the vertical well and penetrated the whole reservoir. The fracture permeability is a function of pressure. The relationship of fracture permeability and pressure is obtained from the experimental data of Lei et al.’s work [30]. The flow in the microfracture obeys the Darcy flow law, while the matrix is stress sensitivity medium, and the fluid steadily flows from the matrix to the microfracture and then flows to the well. The top and bottom of the reservoir are the closure. The SRV matrix and un-SRV matrix have the same reservoir properties, which include permeability, porosity, and medium stress sensitivity. The pressure and flow rate in the inner boundary of un-SRV, respectively, equal to that in the outer boundary of SRV. The flow rate of the volume fractured well is constant. The reservoir outer boundary is the closure. The reservoir and fracture are tiny compressible medium, and the process of fluid flow is isothermal.Figure 2
Schematic for the physical model of the volume fracturing vertical well.
(a)
3D schematic
(b)
Horizontal cross-section view
(c)
Vertical cross-section view


### 2.2. Mathematical Model of Volume Fracturing Vertical Well

For un-SRV, the mass conservation equation is as follows:(1)1r∂∂r−rρvM=∂∂tϕMρ,where ρ is the fluid density (g/cm3), r is the radial distance (cm), vM is the flow velocity in the un-SRV region (cm/s), ϕM is the porosity of the un-SRV region (dimensionless), and t is the simulation time (s).Considering low-velocity non-Darcy flow, the flow equation is as follows:(2)vM=−kMμ∂pM∂r−G,where μ is the fluid viscosity (mPa·s), kM is the permeability of the un-SRV region (mD), pM is the pressure in the un-SRV region (105 Pa), and G is the threshold pressure gradient (105 Pa/cm).Substituting Equation (2) into Equation (1), we can obtain the flow controlling equation in un-SRV as follows:
(3)1r∂∂rrρkMμ∂pM∂r−G=∂∂tϕMρ.According to the microcompressible characteristic of un-SRV and fluid, we can obtain the expressions of the compression factor as follows:(4)cM=1ϕM∂ϕM∂pM,(5)cl=1ρl∂ρl∂pM,where cM and cl are compressive factors of the un-SRV region and fluid (10-5 Pa/cm), respectively.Simplifying Equation (3) by using Equations (4) and (5), the expression of the flow controlling equation is
(6)1r∂∂rrρkMμ∂pM∂r−G=ϕMρcM+cl∂pM∂t.Because the tight reservoir is stress sensitivity medium, the permeability of Equation (6) is a function of pressure, rather than a constant. The permeability of the tight reservoir can be calculated by [24, 31]
(7)γM=1kM∂kM∂pM,where γM is the permeability module of un-SRV (10-5 Pa).For SRV, the mass conservation equation flowing in the fracture is(8)1r∂∂r−rρvf−∂∂tϕmρ=∂∂tϕfρ,where vf is the fluid flow velocity (cm/s), ϕf is the fracture porosity (dimensionless), and ϕm is the porosity of SRV medium (dimensionless).The equation describing the fluid flow in the fracture is as follows:(9)vf=−kfμ∂pf∂r,where kf is the fracture permeability (mD) and pf is the fracture pressure (105 Pa).In SRV, the controlling equation for crossflow from the matrix to fracture is(10)αkmμpf−pm=∂∂tϕmρ,where α is the shape factor, which is utilized to characterize the flow between the matrix and the fracture (dimensionless), km is the permeability of SRV medium (mD), and pm is the medium pressure in SRV (105 Pa).According to the microcompressible characteristic of the matrix and fracture in SRV, we can get two expressions as follows:(11)cf=1ϕf∂ϕf∂pf,(12)cm=1ϕm∂ϕm∂pm,where cf and cm are, respectively, the compressive factors of the fracture and medium (10-5 Pa/cm).By Equations (9)–(12), Equation (8) can be simplified to be
(13)1r∂∂rrρkfμ∂pf∂r−ϕmρcl+cm∂pm∂t=ϕfρcl+cf∂pf∂t.Similarly, by Equation (12), Equation (10) can be rewritten as
(14)αkmμpf−pm=ϕmρcl+cm∂pm∂t.The permeability used in Equations (13) and (14) is not a constant, but a function of pressure. The relationship between matrix permeability in SRV and pressure is [24]
(15)γm=1km∂km∂p,where γm are the permeability modulus of the SRV matrix (10-5 Pa).From Lei et al.’s work [30], a polynomial descripting the law of stress sensitivity of fracture is as follows:
(16)kfkfi=c0100peff2+c110peff+c2,where c0, c1, and c2 are the fitting coefficients and peff is the effective stress (105 Pa), the expression of which is peff=pover−pf. pover is overburdened rock stress (105 Pa). By fitting to the experimental data from Lei et al.’s work [30], we can obtain the values of the fitting coefficients: c0=0.25, c1=0.976, and c2=1.0368.The connection conditions between the SRV and un-SRV are as follows:(17)kf∂pf∂rr=rf=kM∂pM∂rr=rf,(18)pfr=rf=pMr=rf,where rf is the outer boundary of SRV (cm).The outer boundary condition is(19)∂pM∂rr=re=0,where re is the outer boundary of un-SRV (cm).The inner boundary condition is(20)kfr∂pf∂rr=rw=μQ2πh,where rw is the well radius (cm).The initial condition is as follows:(21)pf=pm=pM=pin,where pin is the initial pressure (105 Pa).Equations (6), (7), and (13)–(21) compose the mathematical model of the volume fractured vertical well considering the stress sensitivity and low-velocity non-Darcy flow. The next subsection introduces the model solving method.


### 2.3. Model Solution and Verification

#### 2.3.1. Model Solution

From the type of controlling equations (Equations (6) and (13)), we can know that the controlling equations can only be simplified to be strong nonlinear and variable coefficient partial differential equations, although the Pedrosa transformation [32] is used. It is difficult to obtain the analytical solution of the partial differential equation with a strong nonlinear term and variable coefficient. So, the finite difference method [16] is used to solve the model of the volume fracturing vertical well considering the stress sensitivity and low-velocity non-Darcy flow. The flow coefficient term is treated explicitly and changed with the pressure. In order to obtain the discrete grid distribution with characteristics that the grid grows large with being far away from the well, a logarithm transformation is introduced as follows:
(22)lnrrw=x,where x is a constant and is utilized to describe the meshing relationship.By using Equation (22), the simplified Equation (6) is
(23)∂∂xρkMμ∂pM∂x−rρkMμG=rwe2x2ϕMρcM+cl∂pM∂t.By using the finite difference algorithm and upstream weighting method, Equation (23) can be rewritten as
(24)λM,inΔx2pM,i−1n+1−λM,in+λM,i+1nΔx2pM,in+1+λM,i+1nΔx2pM,i+1n+1+rweΔx⋅iλM,in−λM,i+1nΔxG=rweΔx⋅i2ϕMρcM+clinpM,in+1−pM,inΔt,where λM=ρkM/μ; superscripts n and n+1 are the time step (s); Δt and Δx are, respectively, the simulation time interval (s) and discrete unit size (cm); subscript M represents the un-SRV region; and subscripts i−1, i, and i+1 represent the discrete grid position.Similarly, the Equations (13) and (14) can be, respectively, rewritten as
(25)λf,inΔx2pf,i−1n+1−λf,in+λf,i+1nΔx2pf,in+1+λf,i+1nΔx2pf,i+1n+1−rweΔx⋅i2ϕmρcl+cminpm,in+1−pm,inΔt=rweΔx⋅i2ϕfρcf+clinpf,in+1−pf,inΔt,(26)αkmμpf,in+1−pm,in+1=ϕmρcl+cminpm,in+1−pm,inΔt,where λm=ρkm/μ; λf=ρkf/μ; and subscripts f and m are, respectively, the fracture and matrix of the SRV region.The discrete expression of outer boundary conditions is(27)pNn+1=pN−1n+1,where N is the discrete grid number of the reservoir.The discrete expressions of the connection condition between SRV and un-SRV are(28)kfpf,N1n+1−pf,N1−1n+1=kMpM,N1+1n+1−pf,N1n+1,(29)pfx=rweN1⋅Δx=pMx=rweN1⋅Δx,where N1 is the discrete grid number of the SRV.The discrete expression of inner boundary conditions is(30)kfpf,2n+1−pf,1n+1Δx=μQ2πh,where λm, λf, and λM are not a constant, but a function of pressure. In our simulation, we explicitly obtain the value of λm, λf, and λM. In another word, the λm, λf, and λM at (n+1)th time step equal to that at nth time step when calculated the reservoir pressure of (n+1)th time. After the pressure of (n+1)th time is obtained, the λm, λf, and λM at (n+1)th time step are recalculated by using the pressure of (n+1)th time.Thekm,in and kM,in at nth time step are, respectively, calculated by
(31)km,in=kme−γmpin−pm,in,(32)kM,in=kMe−γMpin−pM,in,where subscript is the initial value.The permeability expression of SRV fracturekfn is as follows:
(33)kfn=kfic0100peffn2+c110peffn+c2.Equations (22)–(33) compose the discrete mathematical model for the scenario of the well with a constant flow rate. The above discrete model can be solved by the Gaussian elimination method.


#### 2.3.2. Model Verification

In order to verify our model, we assume that the circle reservoir is not fractured, namely, the reservoir physical properties of SRV are the same to that of un-SRV. And reservoir permeability is constant, and the pseudothreshold pressure gradient (G) is set to zero. Finally, the proposed model can be simplified to be a conventional reservoir with the vertical well. We take a comparison between the results calculated by the simplified model and the commercial software. The simplified model is performed based on MATLAB. The grid spatial distribution of commercial software is the same as the simplified model. The schematic figure of the grid spatial distribution is shown in Figure 3. From Figure 3, we can know that the grid has nonuniform distribution. The closer to the wellbore is, the smaller the grid is.Figure 3
The schematic figure of the grid spatial distribution.The sizes of the simulation reservoir and SRV are, respectively, 100 and 10 m. The discrete grid number of the simulation reservoir and SRV in the radial direction are, respectively, 20 and 13. The flow rate of the volume fracturing vertical and other parameters used in the simulation are in Table1. The definitions of dimensionless parameters are as follows:
(34)tD=kin,ftϕfcf+clμrw2+ϕMcM+clμrw2,(35)pD=2πhkin,fμQpin−p.Table 1
Parameters used in the comparison.
Parameter
Value
Radius of reservoir under vertical well controlling (cm)
10000
Radius of SRV (cm)
1000
Radius of vertical well (cm)
12.4
Reservoir thickness (cm)
600
Reservoir permeability (mD)
0.23
Reservoir porosity
0.15
Porosity of fracture in the SRV
0.02
Reservoir compression factor ((105 Pa)-1)
0.0038
Compression factor of fracture in the SRV ((105 Pa)-1)
0.0056
Fluid compression factor ((105 Pa)-1)
0.006
Fluid viscosity (mPa·s)
5.6
Fluid density (g/cm3)
0.8
Initial pressure (105 Pa)
300
Total grid number
20
Grid number of SRV
13
Permeability modulus of reservoir ((105 Pa)-1)
0
Pseudothreshold pressure gradient (105 Pa/cm)
0
Well flow rate (cm3/s)
300Using the grid distribution and parameters above, the comparing result of our model and commercial software is conducted (Figure4). From Figure 4, we can know that the dimensionless pressure and dimensionless pressure derivative calculated by the simplified model considerately agree with that of the commercial software. That verifies the correctness of our model and its solution method. The next subsection conducts the effects of key parameters, including stress sensitivity and pseudothreshold pressure gradient, on dimensionless pressure and pressure derivative.Figure 4
Comparison result between the paper model and commercial software.


## 2.1. Physical Model of Volume Fracturing Vertical Well

According to microseismic monitoring mapping [24] and the result of indoor core fracturing experiment [28, 29], we can know that the reservoir after volume fracturing can be divided into two regions, one is the SRV and another is the un-SRV. The SRV is circular, and many fractures are developed in the SRV (Figure 1).Figure 1
The reservoir morphology after vertically fractured.
(a)
Microseismic monitoring mapping [24]
(b)
Fracture morphology of unconventional reservoir core [28, 29]As the statement above, the physical model of the vertical well with volume fracturing can be obtained (Figure2). As seen from Figure 2, the tight reservoir is composed of two regions, which are, respectively, the SRV and un-SRV. The un-SRV is treated as a single porosity medium. And the SRV is idealized by double porosity medium [19]. The un-SRV is stress sensitivity medium, and the flow in un-SRV obeys the low-velocity non-Darcy flow. In SRV, the microfracture is developed, and all microfractures are crossed through the vertical well and penetrated the whole reservoir. The fracture permeability is a function of pressure. The relationship of fracture permeability and pressure is obtained from the experimental data of Lei et al.’s work [30]. The flow in the microfracture obeys the Darcy flow law, while the matrix is stress sensitivity medium, and the fluid steadily flows from the matrix to the microfracture and then flows to the well. The top and bottom of the reservoir are the closure. The SRV matrix and un-SRV matrix have the same reservoir properties, which include permeability, porosity, and medium stress sensitivity. The pressure and flow rate in the inner boundary of un-SRV, respectively, equal to that in the outer boundary of SRV. The flow rate of the volume fractured well is constant. The reservoir outer boundary is the closure. The reservoir and fracture are tiny compressible medium, and the process of fluid flow is isothermal.Figure 2
Schematic for the physical model of the volume fracturing vertical well.
(a)
3D schematic
(b)
Horizontal cross-section view
(c)
Vertical cross-section view


## 2.2. Mathematical Model of Volume Fracturing Vertical Well

For un-SRV, the mass conservation equation is as follows:(1)1r∂∂r−rρvM=∂∂tϕMρ,where ρ is the fluid density (g/cm3), r is the radial distance (cm), vM is the flow velocity in the un-SRV region (cm/s), ϕM is the porosity of the un-SRV region (dimensionless), and t is the simulation time (s).Considering low-velocity non-Darcy flow, the flow equation is as follows:(2)vM=−kMμ∂pM∂r−G,where μ is the fluid viscosity (mPa·s), kM is the permeability of the un-SRV region (mD), pM is the pressure in the un-SRV region (105 Pa), and G is the threshold pressure gradient (105 Pa/cm).Substituting Equation (2) into Equation (1), we can obtain the flow controlling equation in un-SRV as follows:
(3)1r∂∂rrρkMμ∂pM∂r−G=∂∂tϕMρ.According to the microcompressible characteristic of un-SRV and fluid, we can obtain the expressions of the compression factor as follows:(4)cM=1ϕM∂ϕM∂pM,(5)cl=1ρl∂ρl∂pM,where cM and cl are compressive factors of the un-SRV region and fluid (10-5 Pa/cm), respectively.Simplifying Equation (3) by using Equations (4) and (5), the expression of the flow controlling equation is
(6)1r∂∂rrρkMμ∂pM∂r−G=ϕMρcM+cl∂pM∂t.Because the tight reservoir is stress sensitivity medium, the permeability of Equation (6) is a function of pressure, rather than a constant. The permeability of the tight reservoir can be calculated by [24, 31]
(7)γM=1kM∂kM∂pM,where γM is the permeability module of un-SRV (10-5 Pa).For SRV, the mass conservation equation flowing in the fracture is(8)1r∂∂r−rρvf−∂∂tϕmρ=∂∂tϕfρ,where vf is the fluid flow velocity (cm/s), ϕf is the fracture porosity (dimensionless), and ϕm is the porosity of SRV medium (dimensionless).The equation describing the fluid flow in the fracture is as follows:(9)vf=−kfμ∂pf∂r,where kf is the fracture permeability (mD) and pf is the fracture pressure (105 Pa).In SRV, the controlling equation for crossflow from the matrix to fracture is(10)αkmμpf−pm=∂∂tϕmρ,where α is the shape factor, which is utilized to characterize the flow between the matrix and the fracture (dimensionless), km is the permeability of SRV medium (mD), and pm is the medium pressure in SRV (105 Pa).According to the microcompressible characteristic of the matrix and fracture in SRV, we can get two expressions as follows:(11)cf=1ϕf∂ϕf∂pf,(12)cm=1ϕm∂ϕm∂pm,where cf and cm are, respectively, the compressive factors of the fracture and medium (10-5 Pa/cm).By Equations (9)–(12), Equation (8) can be simplified to be
(13)1r∂∂rrρkfμ∂pf∂r−ϕmρcl+cm∂pm∂t=ϕfρcl+cf∂pf∂t.Similarly, by Equation (12), Equation (10) can be rewritten as
(14)αkmμpf−pm=ϕmρcl+cm∂pm∂t.The permeability used in Equations (13) and (14) is not a constant, but a function of pressure. The relationship between matrix permeability in SRV and pressure is [24]
(15)γm=1km∂km∂p,where γm are the permeability modulus of the SRV matrix (10-5 Pa).From Lei et al.’s work [30], a polynomial descripting the law of stress sensitivity of fracture is as follows:
(16)kfkfi=c0100peff2+c110peff+c2,where c0, c1, and c2 are the fitting coefficients and peff is the effective stress (105 Pa), the expression of which is peff=pover−pf. pover is overburdened rock stress (105 Pa). By fitting to the experimental data from Lei et al.’s work [30], we can obtain the values of the fitting coefficients: c0=0.25, c1=0.976, and c2=1.0368.The connection conditions between the SRV and un-SRV are as follows:(17)kf∂pf∂rr=rf=kM∂pM∂rr=rf,(18)pfr=rf=pMr=rf,where rf is the outer boundary of SRV (cm).The outer boundary condition is(19)∂pM∂rr=re=0,where re is the outer boundary of un-SRV (cm).The inner boundary condition is(20)kfr∂pf∂rr=rw=μQ2πh,where rw is the well radius (cm).The initial condition is as follows:(21)pf=pm=pM=pin,where pin is the initial pressure (105 Pa).Equations (6), (7), and (13)–(21) compose the mathematical model of the volume fractured vertical well considering the stress sensitivity and low-velocity non-Darcy flow. The next subsection introduces the model solving method.


## 2.3. Model Solution and Verification

### 2.3.1. Model Solution

From the type of controlling equations (Equations (6) and (13)), we can know that the controlling equations can only be simplified to be strong nonlinear and variable coefficient partial differential equations, although the Pedrosa transformation [32] is used. It is difficult to obtain the analytical solution of the partial differential equation with a strong nonlinear term and variable coefficient. So, the finite difference method [16] is used to solve the model of the volume fracturing vertical well considering the stress sensitivity and low-velocity non-Darcy flow. The flow coefficient term is treated explicitly and changed with the pressure. In order to obtain the discrete grid distribution with characteristics that the grid grows large with being far away from the well, a logarithm transformation is introduced as follows:
(22)lnrrw=x,where x is a constant and is utilized to describe the meshing relationship.By using Equation (22), the simplified Equation (6) is
(23)∂∂xρkMμ∂pM∂x−rρkMμG=rwe2x2ϕMρcM+cl∂pM∂t.By using the finite difference algorithm and upstream weighting method, Equation (23) can be rewritten as
(24)λM,inΔx2pM,i−1n+1−λM,in+λM,i+1nΔx2pM,in+1+λM,i+1nΔx2pM,i+1n+1+rweΔx⋅iλM,in−λM,i+1nΔxG=rweΔx⋅i2ϕMρcM+clinpM,in+1−pM,inΔt,where λM=ρkM/μ; superscripts n and n+1 are the time step (s); Δt and Δx are, respectively, the simulation time interval (s) and discrete unit size (cm); subscript M represents the un-SRV region; and subscripts i−1, i, and i+1 represent the discrete grid position.Similarly, the Equations (13) and (14) can be, respectively, rewritten as
(25)λf,inΔx2pf,i−1n+1−λf,in+λf,i+1nΔx2pf,in+1+λf,i+1nΔx2pf,i+1n+1−rweΔx⋅i2ϕmρcl+cminpm,in+1−pm,inΔt=rweΔx⋅i2ϕfρcf+clinpf,in+1−pf,inΔt,(26)αkmμpf,in+1−pm,in+1=ϕmρcl+cminpm,in+1−pm,inΔt,where λm=ρkm/μ; λf=ρkf/μ; and subscripts f and m are, respectively, the fracture and matrix of the SRV region.The discrete expression of outer boundary conditions is(27)pNn+1=pN−1n+1,where N is the discrete grid number of the reservoir.The discrete expressions of the connection condition between SRV and un-SRV are(28)kfpf,N1n+1−pf,N1−1n+1=kMpM,N1+1n+1−pf,N1n+1,(29)pfx=rweN1⋅Δx=pMx=rweN1⋅Δx,where N1 is the discrete grid number of the SRV.The discrete expression of inner boundary conditions is(30)kfpf,2n+1−pf,1n+1Δx=μQ2πh,where λm, λf, and λM are not a constant, but a function of pressure. In our simulation, we explicitly obtain the value of λm, λf, and λM. In another word, the λm, λf, and λM at (n+1)th time step equal to that at nth time step when calculated the reservoir pressure of (n+1)th time. After the pressure of (n+1)th time is obtained, the λm, λf, and λM at (n+1)th time step are recalculated by using the pressure of (n+1)th time.Thekm,in and kM,in at nth time step are, respectively, calculated by
(31)km,in=kme−γmpin−pm,in,(32)kM,in=kMe−γMpin−pM,in,where subscript is the initial value.The permeability expression of SRV fracturekfn is as follows:
(33)kfn=kfic0100peffn2+c110peffn+c2.Equations (22)–(33) compose the discrete mathematical model for the scenario of the well with a constant flow rate. The above discrete model can be solved by the Gaussian elimination method.


### 2.3.2. Model Verification

In order to verify our model, we assume that the circle reservoir is not fractured, namely, the reservoir physical properties of SRV are the same to that of un-SRV. And reservoir permeability is constant, and the pseudothreshold pressure gradient (G) is set to zero. Finally, the proposed model can be simplified to be a conventional reservoir with the vertical well. We take a comparison between the results calculated by the simplified model and the commercial software. The simplified model is performed based on MATLAB. The grid spatial distribution of commercial software is the same as the simplified model. The schematic figure of the grid spatial distribution is shown in Figure 3. From Figure 3, we can know that the grid has nonuniform distribution. The closer to the wellbore is, the smaller the grid is.Figure 3
The schematic figure of the grid spatial distribution.The sizes of the simulation reservoir and SRV are, respectively, 100 and 10 m. The discrete grid number of the simulation reservoir and SRV in the radial direction are, respectively, 20 and 13. The flow rate of the volume fracturing vertical and other parameters used in the simulation are in Table1. The definitions of dimensionless parameters are as follows:
(34)tD=kin,ftϕfcf+clμrw2+ϕMcM+clμrw2,(35)pD=2πhkin,fμQpin−p.Table 1
Parameters used in the comparison.
Parameter
Value
Radius of reservoir under vertical well controlling (cm)
10000
Radius of SRV (cm)
1000
Radius of vertical well (cm)
12.4
Reservoir thickness (cm)
600
Reservoir permeability (mD)
0.23
Reservoir porosity
0.15
Porosity of fracture in the SRV
0.02
Reservoir compression factor ((105 Pa)-1)
0.0038
Compression factor of fracture in the SRV ((105 Pa)-1)
0.0056
Fluid compression factor ((105 Pa)-1)
0.006
Fluid viscosity (mPa·s)
5.6
Fluid density (g/cm3)
0.8
Initial pressure (105 Pa)
300
Total grid number
20
Grid number of SRV
13
Permeability modulus of reservoir ((105 Pa)-1)
0
Pseudothreshold pressure gradient (105 Pa/cm)
0
Well flow rate (cm3/s)
300Using the grid distribution and parameters above, the comparing result of our model and commercial software is conducted (Figure4). From Figure 4, we can know that the dimensionless pressure and dimensionless pressure derivative calculated by the simplified model considerately agree with that of the commercial software. That verifies the correctness of our model and its solution method. The next subsection conducts the effects of key parameters, including stress sensitivity and pseudothreshold pressure gradient, on dimensionless pressure and pressure derivative.Figure 4
Comparison result between the paper model and commercial software.


## 2.3.1. Model Solution

From the type of controlling equations (Equations (6) and (13)), we can know that the controlling equations can only be simplified to be strong nonlinear and variable coefficient partial differential equations, although the Pedrosa transformation [32] is used. It is difficult to obtain the analytical solution of the partial differential equation with a strong nonlinear term and variable coefficient. So, the finite difference method [16] is used to solve the model of the volume fracturing vertical well considering the stress sensitivity and low-velocity non-Darcy flow. The flow coefficient term is treated explicitly and changed with the pressure. In order to obtain the discrete grid distribution with characteristics that the grid grows large with being far away from the well, a logarithm transformation is introduced as follows:
(22)lnrrw=x,where x is a constant and is utilized to describe the meshing relationship.By using Equation (22), the simplified Equation (6) is
(23)∂∂xρkMμ∂pM∂x−rρkMμG=rwe2x2ϕMρcM+cl∂pM∂t.By using the finite difference algorithm and upstream weighting method, Equation (23) can be rewritten as
(24)λM,inΔx2pM,i−1n+1−λM,in+λM,i+1nΔx2pM,in+1+λM,i+1nΔx2pM,i+1n+1+rweΔx⋅iλM,in−λM,i+1nΔxG=rweΔx⋅i2ϕMρcM+clinpM,in+1−pM,inΔt,where λM=ρkM/μ; superscripts n and n+1 are the time step (s); Δt and Δx are, respectively, the simulation time interval (s) and discrete unit size (cm); subscript M represents the un-SRV region; and subscripts i−1, i, and i+1 represent the discrete grid position.Similarly, the Equations (13) and (14) can be, respectively, rewritten as
(25)λf,inΔx2pf,i−1n+1−λf,in+λf,i+1nΔx2pf,in+1+λf,i+1nΔx2pf,i+1n+1−rweΔx⋅i2ϕmρcl+cminpm,in+1−pm,inΔt=rweΔx⋅i2ϕfρcf+clinpf,in+1−pf,inΔt,(26)αkmμpf,in+1−pm,in+1=ϕmρcl+cminpm,in+1−pm,inΔt,where λm=ρkm/μ; λf=ρkf/μ; and subscripts f and m are, respectively, the fracture and matrix of the SRV region.The discrete expression of outer boundary conditions is(27)pNn+1=pN−1n+1,where N is the discrete grid number of the reservoir.The discrete expressions of the connection condition between SRV and un-SRV are(28)kfpf,N1n+1−pf,N1−1n+1=kMpM,N1+1n+1−pf,N1n+1,(29)pfx=rweN1⋅Δx=pMx=rweN1⋅Δx,where N1 is the discrete grid number of the SRV.The discrete expression of inner boundary conditions is(30)kfpf,2n+1−pf,1n+1Δx=μQ2πh,where λm, λf, and λM are not a constant, but a function of pressure. In our simulation, we explicitly obtain the value of λm, λf, and λM. In another word, the λm, λf, and λM at (n+1)th time step equal to that at nth time step when calculated the reservoir pressure of (n+1)th time. After the pressure of (n+1)th time is obtained, the λm, λf, and λM at (n+1)th time step are recalculated by using the pressure of (n+1)th time.Thekm,in and kM,in at nth time step are, respectively, calculated by
(31)km,in=kme−γmpin−pm,in,(32)kM,in=kMe−γMpin−pM,in,where subscript is the initial value.The permeability expression of SRV fracturekfn is as follows:
(33)kfn=kfic0100peffn2+c110peffn+c2.Equations (22)–(33) compose the discrete mathematical model for the scenario of the well with a constant flow rate. The above discrete model can be solved by the Gaussian elimination method.


## 2.3.2. Model Verification

In order to verify our model, we assume that the circle reservoir is not fractured, namely, the reservoir physical properties of SRV are the same to that of un-SRV. And reservoir permeability is constant, and the pseudothreshold pressure gradient (G) is set to zero. Finally, the proposed model can be simplified to be a conventional reservoir with the vertical well. We take a comparison between the results calculated by the simplified model and the commercial software. The simplified model is performed based on MATLAB. The grid spatial distribution of commercial software is the same as the simplified model. The schematic figure of the grid spatial distribution is shown in Figure 3. From Figure 3, we can know that the grid has nonuniform distribution. The closer to the wellbore is, the smaller the grid is.Figure 3
The schematic figure of the grid spatial distribution.The sizes of the simulation reservoir and SRV are, respectively, 100 and 10 m. The discrete grid number of the simulation reservoir and SRV in the radial direction are, respectively, 20 and 13. The flow rate of the volume fracturing vertical and other parameters used in the simulation are in Table1. The definitions of dimensionless parameters are as follows:
(34)tD=kin,ftϕfcf+clμrw2+ϕMcM+clμrw2,(35)pD=2πhkin,fμQpin−p.Table 1
Parameters used in the comparison.
Parameter
Value
Radius of reservoir under vertical well controlling (cm)
10000
Radius of SRV (cm)
1000
Radius of vertical well (cm)
12.4
Reservoir thickness (cm)
600
Reservoir permeability (mD)
0.23
Reservoir porosity
0.15
Porosity of fracture in the SRV
0.02
Reservoir compression factor ((105 Pa)-1)
0.0038
Compression factor of fracture in the SRV ((105 Pa)-1)
0.0056
Fluid compression factor ((105 Pa)-1)
0.006
Fluid viscosity (mPa·s)
5.6
Fluid density (g/cm3)
0.8
Initial pressure (105 Pa)
300
Total grid number
20
Grid number of SRV
13
Permeability modulus of reservoir ((105 Pa)-1)
0
Pseudothreshold pressure gradient (105 Pa/cm)
0
Well flow rate (cm3/s)
300Using the grid distribution and parameters above, the comparing result of our model and commercial software is conducted (Figure4). From Figure 4, we can know that the dimensionless pressure and dimensionless pressure derivative calculated by the simplified model considerately agree with that of the commercial software. That verifies the correctness of our model and its solution method. The next subsection conducts the effects of key parameters, including stress sensitivity and pseudothreshold pressure gradient, on dimensionless pressure and pressure derivative.Figure 4
Comparison result between the paper model and commercial software.


## 3. The Results and Discussions

### 3.1. Flow Regimes

In order to divide flow regimes of the volume fractured vertical well considering stress sensitivity and low-velocity non-Darcy, we set the permeability modulus of the reservoir matrix to 0.001 (105 Pa)-1. The pseudothreshold pressure gradient (G) equals to 0.001 (105 Pa/cm). The fitting factors used to descript fracture sensitivity are, respectively, c0=0.25, c1=0.976, and c2=1.0368. Other parameters used are in Table 1. We obtain the curves of dimensionless pressure and pressure derivative, which can be seen in Figure 5.Figure 5
Dimensionless pressure and pressure derivative curve for dividing the flow regimes.As seen from Figure5, there are six flow regimes presented by the type curves of the volume fractured vertical well considering stress sensitivity and low-velocity non-Darcy in the tight reservoir. Stage 1 is the fracture linear flow regime. In this regime, the fluid flows in the fracture, and there is no interference between fractures. This flow regime is easy to be vanished. Stage 2 is the early transition flow regime [21]. The dimensionless pressure derivative goes downward in this stage. Stage 3 is the radial flow regime, in which the interferences between fractures occur, and the fluid radially flows to the wellbore. Stage 4 is the crossflow regime. With the decrease of pressure in the fracture, the pressure difference between the fracture and the matrix increases. So, the crossflow between the matrix and the fracture occurs. Stage 5 is the advanced transition flow regime. In these flow regimes, the fluid flows both in SRV and un-SRV. Because of poor reservoir physical properties of un-SRV when compared with that of SRV, the pressure derivative goes upward, which is similar to the closed boundary [21]. The stage 6 is the boundary controlling the flow regime. The dimensionless pressure and pressure derivative curves go upward, and its slopes equal to 1 when the boundary is closed.


### 3.2. The Effects of Stress Sensitivity and Low-Velocity Non-Darcy Flow

In this subsection, the effects of stress sensitivity and low-velocity non-Darcy flow on dimensionless pressure and pressure derivative are analyzed.The fitting factors describing fracture stress sensitivity are, respectively,c0=0.25, c1=0.976, and c2=1.0368. The pseudothreshold pressure gradient equals to 0.001 (105 Pa/cm). The effects of reservoir permeability modulus (the value of γm and γM is the same) on dimensionless pressure and pressure derivative are analyzed and shown in Figure 6.Figure 6
The effect of reservoir permeability modulus.As seen from Figure6, the reservoir permeability modulus mainly affects flow regimes five and six, which are, respectively, advanced transition flow regime and boundary controlling flow regime. Due to the high fracture conductivity, the pressure drop in the fracture is tiny, which results that the steady crossflowing rate between fracture and matrix is small. So, the pressure drop in the SRV matrix is small, and the matrix permeability damage is slight. However, in the advanced transition flow regime and boundary controlling flow regime, a greater pressure decreasing in the reservoir is formed. So, the reservoir permeability modulus has great effects on the dimensionless pressure and pressure derivative. The larger reservoir permeability modulus corresponds to a greater damage of the tight reservoir permeability when the same pressure drop occurs. An increase in the tight reservoir permeability damage results in an increase in the pressure drop when the flow rate is a constant. That results the dimensionless pressure and pressure derivative to go upward as shown in Figure 5.In order to investigate the effect of fracture stress sensitivity, we obtain the fracture sensitivity data from Lei et al.’s work [30] (Figure 7). From Figure 7, the fracture stress sensitivity becomes significant from cores 1 to 3. By fitting the experimental data, we can obtain the fitting coefficients of cores 1 to 3. The fitting coefficients of cores 1 to 3 can be seen in Table 2.Figure 7
Fracture stress sensitivity data from Lei et al.’s work [30].Table 2
The fitting coefficients of cores 1 to 3.
Core no.
c0
c1
c2
1
0.15
-0.635
1.0625
2
0.18
-0.705
1.0456
3
0.25
-0.976
1.0368By using the data in Tables1 and 2, the effect of fracture stress sensitivity is conducted. In this sensitive analysis, the reservoir permeability modulus is 0.001 (105 Pa)-1. The pseudothreshold pressure gradient equals to 0.001 (105 Pa/cm). The effects of fracture stress sensitivity on dimensionless pressure and pressure derivative are seen in Figure 8. From Figure 8, we can obtain that the fracture sensitivity mainly influences the advanced transition flow regime and boundary flow regime. With the increase of fracture stress sensitivity (from cases 1 to 3), the dimensionless pressure and pressure derivative decrease at the advanced transition flow regime and boundary flow regime. This is due to the increase of stress sensitivity indicating an increase in the damage of fracture permeability, which results in a large pressure drop when the flow rate is constant.Figure 8
The effect of fracture stress sensitivity.The permeability modulus of the matrix is 0.001 (105 Pa)-1. The fitting factors used to describe fracture stress sensitivity are, respectively, c0=0.25, c1=0.976, and c2=1.0368. And other parameters are in Table 1. The effects of pseudothreshold pressure gradient (G) on dimensionless pressure and pressure derivative are analyzed and shown in Figure 9.Figure 9
The effect of pseudothreshold pressure gradient (G).From Figure9, we get that the pseudothreshold pressure gradient has great influence on the dimensionless pressure and pressure derivative. With the increase of G, the dimensionless pressure and pressure derivative increase. The reason for this case is that an increase in the pseudothreshold pressure gradient is corresponding to large flow resistance, which requires a high driving pressure to maintain a constant flow rate of the volume fracturing well. So, the dimensionless pressure and pressure derivative increase when the pseudothreshold pressure gradient increases.


## 3.1. Flow Regimes

In order to divide flow regimes of the volume fractured vertical well considering stress sensitivity and low-velocity non-Darcy, we set the permeability modulus of the reservoir matrix to 0.001 (105 Pa)-1. The pseudothreshold pressure gradient (G) equals to 0.001 (105 Pa/cm). The fitting factors used to descript fracture sensitivity are, respectively, c0=0.25, c1=0.976, and c2=1.0368. Other parameters used are in Table 1. We obtain the curves of dimensionless pressure and pressure derivative, which can be seen in Figure 5.Figure 5
Dimensionless pressure and pressure derivative curve for dividing the flow regimes.As seen from Figure5, there are six flow regimes presented by the type curves of the volume fractured vertical well considering stress sensitivity and low-velocity non-Darcy in the tight reservoir. Stage 1 is the fracture linear flow regime. In this regime, the fluid flows in the fracture, and there is no interference between fractures. This flow regime is easy to be vanished. Stage 2 is the early transition flow regime [21]. The dimensionless pressure derivative goes downward in this stage. Stage 3 is the radial flow regime, in which the interferences between fractures occur, and the fluid radially flows to the wellbore. Stage 4 is the crossflow regime. With the decrease of pressure in the fracture, the pressure difference between the fracture and the matrix increases. So, the crossflow between the matrix and the fracture occurs. Stage 5 is the advanced transition flow regime. In these flow regimes, the fluid flows both in SRV and un-SRV. Because of poor reservoir physical properties of un-SRV when compared with that of SRV, the pressure derivative goes upward, which is similar to the closed boundary [21]. The stage 6 is the boundary controlling the flow regime. The dimensionless pressure and pressure derivative curves go upward, and its slopes equal to 1 when the boundary is closed.


## 3.2. The Effects of Stress Sensitivity and Low-Velocity Non-Darcy Flow

In this subsection, the effects of stress sensitivity and low-velocity non-Darcy flow on dimensionless pressure and pressure derivative are analyzed.The fitting factors describing fracture stress sensitivity are, respectively,c0=0.25, c1=0.976, and c2=1.0368. The pseudothreshold pressure gradient equals to 0.001 (105 Pa/cm). The effects of reservoir permeability modulus (the value of γm and γM is the same) on dimensionless pressure and pressure derivative are analyzed and shown in Figure 6.Figure 6
The effect of reservoir permeability modulus.As seen from Figure6, the reservoir permeability modulus mainly affects flow regimes five and six, which are, respectively, advanced transition flow regime and boundary controlling flow regime. Due to the high fracture conductivity, the pressure drop in the fracture is tiny, which results that the steady crossflowing rate between fracture and matrix is small. So, the pressure drop in the SRV matrix is small, and the matrix permeability damage is slight. However, in the advanced transition flow regime and boundary controlling flow regime, a greater pressure decreasing in the reservoir is formed. So, the reservoir permeability modulus has great effects on the dimensionless pressure and pressure derivative. The larger reservoir permeability modulus corresponds to a greater damage of the tight reservoir permeability when the same pressure drop occurs. An increase in the tight reservoir permeability damage results in an increase in the pressure drop when the flow rate is a constant. That results the dimensionless pressure and pressure derivative to go upward as shown in Figure 5.In order to investigate the effect of fracture stress sensitivity, we obtain the fracture sensitivity data from Lei et al.’s work [30] (Figure 7). From Figure 7, the fracture stress sensitivity becomes significant from cores 1 to 3. By fitting the experimental data, we can obtain the fitting coefficients of cores 1 to 3. The fitting coefficients of cores 1 to 3 can be seen in Table 2.Figure 7
Fracture stress sensitivity data from Lei et al.’s work [30].Table 2
The fitting coefficients of cores 1 to 3.
Core no.
c0
c1
c2
1
0.15
-0.635
1.0625
2
0.18
-0.705
1.0456
3
0.25
-0.976
1.0368By using the data in Tables1 and 2, the effect of fracture stress sensitivity is conducted. In this sensitive analysis, the reservoir permeability modulus is 0.001 (105 Pa)-1. The pseudothreshold pressure gradient equals to 0.001 (105 Pa/cm). The effects of fracture stress sensitivity on dimensionless pressure and pressure derivative are seen in Figure 8. From Figure 8, we can obtain that the fracture sensitivity mainly influences the advanced transition flow regime and boundary flow regime. With the increase of fracture stress sensitivity (from cases 1 to 3), the dimensionless pressure and pressure derivative decrease at the advanced transition flow regime and boundary flow regime. This is due to the increase of stress sensitivity indicating an increase in the damage of fracture permeability, which results in a large pressure drop when the flow rate is constant.Figure 8
The effect of fracture stress sensitivity.The permeability modulus of the matrix is 0.001 (105 Pa)-1. The fitting factors used to describe fracture stress sensitivity are, respectively, c0=0.25, c1=0.976, and c2=1.0368. And other parameters are in Table 1. The effects of pseudothreshold pressure gradient (G) on dimensionless pressure and pressure derivative are analyzed and shown in Figure 9.Figure 9
The effect of pseudothreshold pressure gradient (G).From Figure9, we get that the pseudothreshold pressure gradient has great influence on the dimensionless pressure and pressure derivative. With the increase of G, the dimensionless pressure and pressure derivative increase. The reason for this case is that an increase in the pseudothreshold pressure gradient is corresponding to large flow resistance, which requires a high driving pressure to maintain a constant flow rate of the volume fracturing well. So, the dimensionless pressure and pressure derivative increase when the pseudothreshold pressure gradient increases.


## 4. Conclusion

(1)
A numerical pressure analysis model of the volume fracturing vertical well considering stress sensitivity and low-velocity non-Darcy flow is proposed, and the numerical solution is obtained by a finite difference method
(2)
The flow regimes for the volume fracturing vertical well considering stress sensitivity and low-velocity non-Darcy flow can be identified by the dimensionless pressure and pressure derivative curves. Six flow stages are identified and include the fracture linear flow regime, early transition flow regime, radial flow regime, crossflow regime, advanced transition flow regime, and boundary controlling flow regime
(3)
The stress sensitivity and pseudopressure gradient have a great effect on the dimensionless pressure and pressure derivative. With the increase of reservoir stress sensitivity, the pressure and pressure derivative are upward at the advanced transition flow and boundary controlling regimes. However, the pressure and pressure derivative are downward at the advanced transition flow and boundary controlling regimes when the fracture sensitivity increases. An increasing pseudopressure gradient is corresponding to high dimensionless pressure and pressure derivative



---

*Source: 2046061-2019-11-20.xml*