Advanced modeling and simulation of flow in subsurface reservoirs with fractures and wells for a sustainable industry
Open Access
Numéro
Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles
Volume 75, 2020
Advanced modeling and simulation of flow in subsurface reservoirs with fractures and wells for a sustainable industry
Numéro d'article 47
Nombre de pages 14
DOI https://doi.org/10.2516/ogst/2020043
Publié en ligne 14 juillet 2020

© L. Mei et al., published by IFP Energies nouvelles, 2020

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Nomenclature

q : Fluid transfer rate between matrix and fracture, [M/(TL3)]

σ : Shape factor, [1/L2]

ρ : Fluid density, [M/L3]

μ : Fluid viscosity, [M/(LT)]

N : Number of sets of fractures (1, 2 or 3)

L : Characteristic length of matrix, [L]

P : Matrix pressure, [M/(LT2)]

P 0 : Initial pressure of matrix, [M/(LT2)]

P ̅ $ \bar{P}$ : Average pressure of matrix, [M/(LT2)]

P f : Fracture pressure, [M/(LT2)]

K m : Matrix permeability, [L2]

ϕ m : Matrix porosity

K a : Permeability constant, [L2]

ϕ a : Porosity constant

D f : Fractal dimension of matrix

D T : Tortuosity fractal dimension of matrix

λ max : Maximum pore diameter of matrix, [L]

A 0 : Unit area, [L2]

c : Total compressibility of reservoir, [(LT2)/M]

t : Transfer flow time, [T]

R : Equivalent radius, [L]

1 Introduction

Naturally fractured reservoirs account for a large proportion of the world’s resources and play an important role in the world’s energy structure. Unconventional resources have significantly transformed the landscape of oil and gas industry in these years. Hydraulic fracturing is a primary technology for economical and effective exploitation of unconventional reservoirs (Li et al., 2015; Wang et al., 2017; Behnoudfar et al., 2019; Vishkai and Gates, 2019), which aims at creating large-scale fracture network to increase production. Therefore, fluids flow mechanism in fractured systems has aroused the strong interest of petroleum engineers, geothermal engineers and geologists, etc.

The dual-porosity model is one of common simulation methods to study the flow properties in fractured reservoirs. The concept of dual-porosity media in fractured reservoirs was firstly proposed by Barenblatt et al. (1960), and then was applied by Warren and Root (1963) to the field of petroleum engineering. As shown in Figure 1, fractured reservoirs are usually supposed to consist of a continuous fracture region that serves as the primary flow channels (high permeability and low porosity), and a discontinuous matrix region which acts as the primary hydrocarbon storage area (low permeability and high porosity). They hypothesized that single-phase of slightly compressible fluid transfer between matrix and fracture occurred under pseudo-steady state conditions. The transfer function was originally formulated as (Warren and Root, 1963):q=σρKmμ(P̅-Pf),$$ q=\sigma \frac{\rho {K}_m}{\mu }\left(\bar{P}-{P}_f\right), $$(1)where q is the fluid transfer rate between matrix and fracture, [M/(TL3)]; σ is the shape factor, [1/L2]; ρ is the fluid density, [M/L3]; μ is the fluid viscosity, [M/(LT)]; Km is the matrix permeability, [L2]; P̅$ \bar{P}$ is the average pressure of matrix, [M/(LT2)]; Pf is the fracture pressure, [M/(LT2)].

thumbnail Fig. 1

Schematic of idealization of fractured reservoirs (Warren and Root, 1963).

The shape factor model for cubic matrix can be written as (Warren and Root, 1963):σ=4N(N+2)L2,$$ \sigma =\frac{4N\left(N+2\right)}{{L}^2}, $$(2)where N is the number of sets of fractures (1, 2 or 3), L is the characteristic length of matrix which is the distance between the fracture surface and the center of matrix, [L]. The shape factors are 12/L2, 32/L2 and 60/L2 for one-dimensional (1D), two-dimensional (2D) and three-dimensional (3D) spaces, respectively.

Furthermore, many scholars have carried out studies on fluids exchange between matrix and fracture based on Warren and Root (1963)’s model. Kazemi et al. (1976) directly extended single-phase transfer of Warren and Root (1963) to two-phase transfer and firstly applied it to the numerical simulation of 3D isotropic dual-porosity porous media. When the length of matrix is the same in each direction, the shape factors are 4/L2, 8/L2 and 12/L2 for 1D, 2D and 3D spaces, respectively. By including the effects of gravity, capillary force and viscous force on transfer function, Thomas et al. (1983) established a 3D three-phase model to simulate fluids flow in fractured reservoirs. By comparing single-porosity fine-grid with the dual-porosity models, they concluded that shape factor was 25/L2 when the two models have good agreement. Coats (1989) solved the diffusion equation of the slightly compressible fluids and derived a transfer function. They found that their shape factor was exactly twice the value of Kazemi et al. (1976). By comparing fine-grid models with experiments, Ueda et al. (1989) believed that the shape factor was twice and three times that of Kazemi et al. (1976) in 1D and 2D spaces, respectively.

However, the application of these shape factors to explain the transient characteristics of matrix pressure in the initial stage of fluids transfer process usually produces large deviation. Therefore, some scholars have conducted many studies on the description of matrix pressure transient behavior. Zimmerman et al. (1993) developed a single-phase dual-porosity model in fractured reservoirs and solved the non-linear ordinary differential equation to obtain matrix-fracture transfer function, which was more accurate than the linear Warren and Root (1963) equation to calculate the flux in the early and late stages. Lim and Aziz (1995) derived shape factor by solving matrix pressure diffusion equation, and used fine-grid numerical simulation to validate the presented shape factor. They considered that both the geometric of matrix and pressure gradient of the whole system have a great impact on shape factor. Noetinger and Estebenet (2000) used continuous time random walks methods to simulate matrix-fracture exchange term and obtained shape factor involved in dual-porosity formulations of fluid flow through fractured reservoirs. Good agreement were found between their results and numerical simulation. Sarma and Aziz (2004) believed that fracture networks were non-orthogonal, and solved the non-orthogonal single-phase pressure diffusion equation to deduce 2D and 3D shape factors. Their results were verified by comparing the discrete fracture model with the dual-porosity model. For different geometries and boundary conditions, Hassanzadeh and Pooladidarvish (2006) introduced the shape factor, and indicated that both the geometric of matrix and pressure change of fracture have an important effect on the shape factor. Ranjbar et al. (2011) solved the single-phase nonlinear diffusivity equation of different pressure depletion regimes to derive the shape factor, and investigated the impact of fracture pressure depletion regime on the shape factor. They further applied the fine-grid numerical simulation technique to verify their shape factor. By solving the saturation diffusion equation in fractured reservoirs, Saboorian-Jooybari et al. (2015) developed the shape factor accounting for both capillary and gravity forces, which was verified using fine-grid numerical simulation. For low permeability fractured reservoirs, He et al. (2017) proposed a new shape factor by solving single-phase pressure diffusion equation. Furthermore, their results shown that the new shape factor could predict production accurately. For tight fractured reservoirs with the existence of the boundary layer and the heterogeneous pressure distribution of matrix, Cao et al. (2019) presented the single-phase matrix-fracture transfer function, which was verified by the numerical simulation and the experimental data.

Fractal geometry theory has been widely used to study the transport properties of porous media. The application of fractal geometry theory in the field of petroleum engineering has received extensive attentions. Both of fracture network (Chang and Yortsos, 1990; Acuna et al., 1995) and matrix (Katz and Thompson, 1985; Krohn, 1988; Adler, 1996) can be represented effectively and conveniently by fractal geometry. Zhou et al. (2000) studied the fractal characteristics of natural fractures, and developed a fractal transfer function model to analyze the effect of pore size distribution of matrix on the fluids transfer between matrix and fracture in fractured reservoirs. Based on the fractal characteristics of fractured reservoirs, Flamenco-Lopez and Camacho-Velazquez (2001) established the fractal dual-porosity equation, and obtained the analytical solution of pressure by Laplace transform method. The importance of transient and pseudo-steady pressure characteristics of fractured reservoirs with fractal geometry in a single-well was discussed. Kong et al. (2009) described the pore characteristics of porous and fractured media through fractal theory and presented the new flow velocity, porosity and permeability models in both porous and fractured media. In their study, the pressure diffusion equations of porous and fractured media were established and their transient pressure characteristics were analyzed. Yao et al. (2012) believed that fractures have fractal characteristics and established two dual-porosity equations based on circular matrix and cylindrical matrix, respectively. They used Laplace method to get the pressure analytical solution, and analyzed the influences of matrix shape and fractal parameters on the transient pressure characteristics of fractured reservoirs. Fan and Ettehadtavakkol (2017) indicated that fractures have fractal characteristics and established induced fractures network distribution model, which was a function of density distribution and permeability/porosity distribution of induced fractures. Lian et al. (2018) used fractal geometry to characterize the different scale distribution of pores and fractures, and presented single-phase radial fluid flow equation by involving nonequilibrium adsorption to study the flow properties in shale gas reservoirs. Wang et al. (2018) combined the dual porosity model with the tri-linear flow model, and used fractal geometry to characterize the heterogeneous distribution of fracture networks and the non-uniform of both matrix/fracture porosity and permeability. By applying Bessel Function and Laplace transform techniques, the flow of fluid through primary hydraulic fractures, stimulated-reservoir-volume and unstimulated reservoir matrices were integrated to derive analytical solutions. The accuracy of the analytical solutions in their study were verified by a synthetic fine-grid numerical simulation example.

In the previous literature, the effect of matrix microstructure on the shape factor is not considered during reservoir development. Fractal geometry theory can relatively accurately characterize the effect of matrix microstructure on fluids transfer between matrix and fracture. Fractal geometry theory has been widely applied in flow characteristics of fractured reservoirs, but its applications in the transfer function and shape factor are still rare. In this paper, we consider the effect of matrix microstructure on fluids transfer between matrix and fracture by means of fractal geometry theory, from which a new shape factor is derived. Single-porosity fine-grid and dual-porosity models are then applied to validate the new shape factor. Finally, the influences of microstructure and characteristic length of matrix on the shape factor are accordingly discussed.

2 Mathematical model

It has been reported that pore surface and pore size distribution of porous media follow fractal features (Katz and Thompson, 1985; Yu and Cheng, 2002). The fractal geometry theory has great advantages in accurate representation of porous media with complex microstructures (Costa, 2006; Erol et al., 2017). Considering the fractal capillary bundle model and Hagen-Poiseuille equation of curved capillary, the fractal permeability and fractal porosity of matrix can be respectively expressed by (Kong et al., 2007):Km=π128A0Dfλmax4DT(3+DT-Df)(xλmax)1-DT=Ka(xλmax)1-DT,$$ {K}_m=\frac{\pi }{128{A}_0}\frac{{D}_f{\lambda }_{\mathrm{max}}^4}{{D}_T\left(3+{D}_T-{D}_f\right)}{\left(\frac{x}{{\lambda }_{\mathrm{max}}}\right)}^{1-{D}_T}={K}_a{\left(\frac{x}{{\lambda }_{\mathrm{max}}}\right)}^{1-{D}_T}, $$(3)ϕm=π4A0DTDfλmax23-DT-Df[1-(λminλmax)3-DT-Df](λmaxx)1-DT=ϕa(λmaxx)1-DT,$$ {\phi }_m=\frac{\pi }{4{A}_0}\frac{{D}_T{D}_f{\lambda }_{\mathrm{max}}^2}{3-{D}_T-{D}_f}\left[1-{\left(\frac{{\lambda }_{\mathrm{min}}}{{\lambda }_{\mathrm{max}}}\right)}^{3-{D}_T-{D}_f}\right]{\left(\frac{{\lambda }_{\mathrm{max}}}{x}\right)}^{1-{D}_T}={\phi }_a{\left(\frac{{\lambda }_{\mathrm{max}}}{x}\right)}^{1-{D}_T}, $$(4)where ϕm is the matrix porosity; A0 is the unit area, [L2]; x is the distance, [L]; Df is the pore fractal dimension, 0 < Df < 2 (in 2D space) and 0 < Df < 3 (in 3D space); λmax is the maximum pore diameter, [L]; DT is the tortuosity fractal dimension, 1 < DT < 2 (in 2D space) and 1 < DT < 3 (in 3D space), DT = 1 represents a straight capillary, DT = 2 represents a plane is completely filled with tortuous line, and DT = 3 represents a 3D space is completely filled with tortuous line; Ka is the permeability constant, [L2]; ϕa is the porosity constant. It is noted that the symbols Df and DT are respectively used for pores and streamlines in matrix in this work.

Substituting equations (3) and (4) into the continuity equation, the matrix fractal pressure diffusion equation can be obtained by (Kong et al., 2007):2Px2+ξ-DTxPx=ϕaoμcKa(xλmax)2(DT-1)Pt,$$ \frac{{\mathrm{\partial }}^2P}{\mathrm{\partial }{x}^2}+\frac{\xi -{D}_T}{x}\frac{\mathrm{\partial }P}{\mathrm{\partial }x}=\frac{{\phi }_{{ao}}\mu c}{{K}_a}{\left(\frac{x}{{\lambda }_{\mathrm{max}}}\right)}^{2\left({D}_T-1\right)}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}, $$(5)where ξ corresponds to 1D flow, 2D flow and 3D flow, which are 1, 2 and 3, respectively. When the capillary is straight, equation (5) can be simplified as the classical pressure diffusion equation.

Before establishing mathematical model, some assumptions need to be considered as follows:

  1. The fractured reservoirs are seen as a dual-porosity model (Fig. 2), where fluids in matrix flow through fractures instead of directly to the bottom of well, and matrix is isotropic.

    thumbnail Fig. 2

    Schematic of dual-porosity media.

  2. The fluid is slightly compressible.

  3. The rock is compressible.

  4. The matrix is represented by the fractal tortuous capillary bundle model.

2.1 Transfer function and shape factor of 1D flow

Supposing matrix has been surrounded by two fractures, which have a spacing L (Fig. 3), fluids flow from matrix to fracture surface, and therefore the fractal pressure diffusion equation in matrix can be given by Kong et al. (2007):2Px2+1-DTxPx=ϕaoμcKa(xλmax)2(DT-1)Pt,$$ \frac{{\mathrm{\partial }}^2P}{\mathrm{\partial }{x}^2}+\frac{1-{D}_T}{x}\frac{\mathrm{\partial }P}{\mathrm{\partial }x}=\frac{{\phi }_{{ao}}\mu c}{{K}_a}{\left(\frac{x}{{\lambda }_{\mathrm{max}}}\right)}^{2\left({D}_T-1\right)}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}, $$(6)where P is the matrix pressure, [M/(LT2)]; ϕao is the value of ϕa under the reference pressure P0; c is the total compressibility of reservoir, [(LT2)/M]; t is the transfer flow time, [T].

thumbnail Fig. 3

Schematic of 1D matrix-fracture system.

Initial and boundary conditions are given by:P=P0,-L/2xL/2, t=0P=Pf,x=±L/2, t>0.$$ \begin{array}{c}P={P}_0,\hspace{0.5em}-L/2\le x\le L/2,\hspace{1em}\enspace t=0\\ P={P}_f,\hspace{0.5em}x=\pm L/2,\enspace \hspace{1em}t>0.\end{array} $$(7)

Using separation variable method and Bessel function and combining equation (6) with equation (7) (see Appendix A for details):P̅-P0Pf-Po=1-n=0(-1)n8×e(2n+1)2Lπ2(2n+1)2δexp(-1M22DT-2π2DT2L2DTt),$$ \frac{\bar{P}-{P}_0}{{P}_f-{P}_o}=1-\sum_{n=0}^{\mathrm{\infty }} \frac{{\left(-1\right)}^n8\times {\mathrm{e}}^{(2n+1{)}^2}}{L{\pi }^2(2n+1{)}^2}\delta \mathrm{exp}\left(-\frac{1}{M}\frac{{2}^{2{D}_T-2}{\pi }^2{D}_T^2}{{L}^{2{D}_T}}t\right), $$(8)where M=ϕaoμcKaλmax2(DT-1)$ M=\frac{{\phi }_{{ao}}\mu c}{{K}_a{\lambda }_{\mathrm{max}}^{2\left({D}_T-1\right)}}$ is a constant.

The partial derivative of equation (8) with respect to time t gives:P̅t1Pf-P0=1M22DT-2π2DT2L2DTn=0(-1)n8×e(2n+1)2Lπ2(2n+1)2δexp(-1M22DT-2π2DT2L2DTt)=1M22DT-2π2DT2L2DT(1-P̅-P0Pf-P0).$$ \frac{\mathrm{\partial }\bar{P}}{\mathrm{\partial }t}\frac{1}{{P}_f-{P}_0}=\frac{1}{M}\frac{{2}^{2{D}_T-2}{\pi }^2{D}_T^2}{{L}^{2{D}_T}}\sum_{n=0}^{\mathrm{\infty }} \frac{{\left(-1\right)}^n8\times {\mathrm{e}}^{(2n+1{)}^2}}{L{\pi }^2(2n+1{)}^2}\delta \mathrm{exp}\left(-\frac{1}{M}\frac{{2}^{2{D}_T-2}{\pi }^2{D}_T^2}{{L}^{2{D}_T}}t\right)=\frac{1}{M}\frac{{2}^{2{D}_T-2}{\pi }^2{D}_T^2}{{L}^{2{D}_T}}\left(1-\frac{\bar{P}-{P}_0}{{P}_f-{P}_0}\right). $$(9)

Equation (9) can be also written as:P̅t=-Kaλmax2(DT-1)ϕaoμc22DT-2π2DT2L2DT(P̅-Pf).$$ \frac{\mathrm{\partial }\bar{P}}{\mathrm{\partial }t}=-\frac{{K}_a{\lambda }_{\mathrm{max}}^{2\left({D}_T-1\right)}}{{\phi }_{{ao}}\mu c}\frac{{2}^{2{D}_T-2}{\pi }^2{D}_T^2}{{L}^{2{D}_T}}\left(\bar{P}-{P}_f\right). $$(10)

The 1D transfer function can be expressed as:q=-(ρϕ)mt=-ρϕaocP̅t(λmaxL)1-DT.$$ q=-\frac{\mathrm{\partial }(\rho \phi {)}_m}{\mathrm{\partial }t}=-\rho {\phi }_{{ao}}c\frac{\mathrm{\partial }\bar{P}}{\mathrm{\partial }t}{\left(\frac{{\lambda }_{\mathrm{max}}}{L}\right)}^{1-{D}_T}. $$(11)

Substituting equation (10) into equation (11) yields:q=ρμKa22DT-2π2DT2L1+DTλmaxDT-1(P̅-Pf).$$ q=\frac{\rho }{\mu }{K}_a\frac{{2}^{2{D}_T-2}{\pi }^2{D}_T^2}{{L}^{1+{D}_T}}{\lambda }_{\mathrm{max}}^{{D}_T-1}\left(\bar{P}-{P}_f\right). $$(12)

Compared to the transfer function (Eq. (1)), shape factor in 1D flow can be expressed as:σ=π222DT-2DT2L1+DTλmaxDT-1.$$ \sigma ={\pi }^2\frac{{2}^{2{D}_T-2}{D}_T^2}{{L}^{1+{D}_T}}{\lambda }_{\mathrm{max}}^{{D}_T-1}. $$(13)

From equation (13), the proposed shape factor considers the effect of flow channel bending and micro pore structure of matrix. Specifically, shape factor is clearly related to tortuosity fractal dimension, maximum pore diameter and characteristic length of matrix.

2.2 Transfer function and shape factor of 2D flow

In 2D model, the matrix is surrounded by four fractures, with two groups of parallel fractures. Fracture space is L. For the 2D matrix-fracture system, the pressure diffusion behavior from matrix to fracture can be similar to a circular region (Fig. 4). The equivalent radius R of matrix is:L2=πR2R=Lπ.$$ {L}^2=\pi {R}^2\Rightarrow R=\frac{L}{\sqrt{\pi }}. $$(14)

thumbnail Fig. 4

Schematic of 2D matrix-fracture system.

The 2D matrix fractal pressure diffusion equation can be given by (Kong et al., 2007):2Pr2+2-DTrPr=ϕaoμcKa(rλmax)2(DT-1)Pt,$$ \frac{{\mathrm{\partial }}^2P}{\mathrm{\partial }{r}^2}+\frac{2-{D}_T}{r}\frac{\mathrm{\partial }P}{\mathrm{\partial }r}=\frac{{\phi }_{{ao}}\mu c}{{K}_a}{\left(\frac{r}{{\lambda }_{\mathrm{max}}}\right)}^{2\left({D}_T-1\right)}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}, $$(15)where r is the distance, [L].

Initial and boundary conditions are given by:P=P0, 0rR, t=0P=Pf, r=R, t>0.$$ \begin{array}{c}P={P}_0,\hspace{1em}\enspace 0\le r\le R,\enspace \hspace{1em}t=0\\ P={P}_f,\hspace{1em}\enspace r=R,\enspace \hspace{1em}t>0.\end{array} $$(16)

Using separation variable method, the analytical solution of equation (15) can be obtained in the form of Bessel function (see Appendix B for details):P=Pf+rDT-12m=12(P0-Pf)ηexp[-1M(DTI1RDT)2t],$$ P={P}_f+{r}^{\frac{{D}_T-1}{2}}\sum_{m=1}^{\mathrm{\infty }} 2\left({P}_0-{P}_f\right)\eta \mathrm{exp}\left[-\frac{1}{M}{\left(\frac{{D}_T{I}_1}{{R}^{{D}_T}}\right)}^2t\right], $$(17)where η=m=11Jv-1(I1)[-DT1-DT2DT(RDTDT)1-3DT2DT(DTI1RDT)-3DT+12DT]×Jv(rDTI1RDT)$ \eta =\sum_{m=1}^{\mathrm{\infty }} \frac{1}{{J}_{v-1}\left({I}_1\right)}\left[-{{D}_T}^{\frac{1-{D}_T}{2{D}_T}}{\left(\frac{{R}^{{D}_T}}{{D}_T}\right)}^{\frac{1-3{D}_T}{2{D}_T}}{\left(\frac{{D}_T{I}_1}{{R}^{{D}_T}}\right)}^{-\frac{3{D}_T+1}{2{D}_T}}\right]\times {J}_v\left(\frac{{r}^{{D}_T}{I}_1}{{R}^{{D}_T}}\right)$, Jv (x) is the Bessel function, v is the order of Bessel function, I1 is the positive zero point of Bessel function (the solution code of I1 is shown in Appendix C), it depends on the value of v, and v=DT-12DT$ v=\frac{{D}_T-1}{2{D}_T}$.

Assuming that the area of matrix Aav is equal to πr2. The average pressure of matrix is as follows:P̅=1πR20RPdAav.$$ \bar{P}=\frac{1}{\pi {R}^2}{\int }_0^R P\mathrm{d}{A}_{{av}}. $$(18)

Substituting equation (17) into equation (18) yields:P̅=Pf+4R2m=1(Pf-P0)0RrDT+12ηdrexp[-1M(DTI1RDT)2t].$$ \bar{P}={P}_f+\frac{4}{{R}^2}\sum_{m=1}^{\mathrm{\infty }} \left({P}_f-{P}_0\right){\int }_0^R {r}^{\frac{{D}_T+1}{2}}\eta \mathrm{d}r\mathrm{exp}\left[-\frac{1}{M}{\left(\frac{{D}_T{I}_1}{{R}^{{D}_T}}\right)}^2t\right]. $$(19)

Both sides of equation (19) are subtracted by P0 and then divided by Pf − P0. Thus, equation (19) can be expressed as:P̅-P0Pf-P0=1+4R2m=10RrDT+12ηdrexp[-1M(DTI1RDT)2t].$$ \frac{\bar{P}-{P}_0}{{P}_f-{P}_0}=1+\frac{4}{{R}^2}\sum_{m=1}^{\mathrm{\infty }} {\int }_0^R {r}^{\frac{{D}_T+1}{2}}\eta \mathrm{d}r\mathrm{exp}\left[-\frac{1}{M}{\left(\frac{{D}_T{I}_1}{{R}^{{D}_T}}\right)}^2t\right]. $$(20)

The partial derivative of equation (20) with respect to t can be expressed as:P̅t1Pf-P0=-1M(DTI1RDT)24R2m=10RrDT+12ηdrexp[-1M(DTI1RDT)2t]=-1M(DTI1RDT)2(P̅-P0Pf-P0-1).$$ \frac{\mathrm{\partial }\bar{P}}{\mathrm{\partial }t}\frac{1}{{P}_f-{P}_0}=-\frac{1}{M}{\left(\frac{{D}_T{I}_1}{{R}^{{D}_T}}\right)}^2\frac{4}{{R}^2}\sum_{m=1}^{\mathrm{\infty }} {\int }_0^R {r}^{\frac{{D}_T+1}{2}}\eta \mathrm{d}r\mathrm{exp}\left[-\frac{1}{M}{\left(\frac{{D}_T{I}_1}{{R}^{{D}_T}}\right)}^2t\right]=-\frac{1}{M}{\left(\frac{{D}_T{I}_1}{{R}^{{D}_T}}\right)}^2\left(\frac{\bar{P}-{P}_0}{{P}_f-{P}_0}-1\right). $$(21)

Equation (21) also can be written as:P̅t=-Kaλmax2(DT-1)ϕaoμc(DTI1RDT)2(P̅-Pf).$$ \frac{\mathrm{\partial }\bar{P}}{\mathrm{\partial }t}=-\frac{{K}_a{\lambda }_{\mathrm{max}}^{2\left({D}_T-1\right)}}{{\phi }_{{ao}}\mu c}{\left(\frac{{D}_T{I}_1}{{R}^{{D}_T}}\right)}^2\left(\bar{P}-{P}_f\right). $$(22)

The 2D transfer function can be expressed as:q=-(ρϕ)mt=-ρϕaocP̅t(λmaxR)1-DT.$$ q=-\frac{\mathrm{\partial }(\rho \phi {)}_m}{\mathrm{\partial }t}=-\rho {\phi }_{{ao}}c\frac{\mathrm{\partial }\bar{P}}{\mathrm{\partial }t}{\left(\frac{{\lambda }_{\mathrm{max}}}{R}\right)}^{1-{D}_T}. $$(23)

Substituting equation (22) into equation (23) yields:q=ρμKaDT2I12R1+DTλmaxDT-1(P̅-Pf).$$ q=\frac{\rho }{\mu }{K}_a\frac{{D}_T^2{I}_1^2}{{R}^{1+{D}_T}}{\lambda }_{\mathrm{max}}^{{D}_T-1}\left(\bar{P}-{P}_f\right). $$(24)

Consequently, the 2D shape factor can be expressed as:σ=DT2I12R1+DTλmaxDT-1.$$ \sigma =\frac{{D}_T^2{I}_1^2}{{R}^{1+{D}_T}}{\lambda }_{\mathrm{max}}^{{D}_T-1}. $$(25)

Substituting equation (14) into equation (25) yields:σ=π1+DT2DT2I12L1+DTλmaxDT-1.$$ \sigma ={\pi }^{\frac{1+{D}_T}{2}}\frac{{D}_T^2{I}_1^2}{{L}^{1+{D}_T}}{\lambda }_{\mathrm{max}}^{{D}_T-1}. $$(26)

The presented 2D shape factor also considers the effect of tortuosity fractal dimension, maximum pore diameter and characteristic length of matrix. Moreover, the value of the positive zero point of Bessel function is related to the value of tortuosity fractal dimension.

2.3 Transfer function and shape factor of 3D flow

In 3D model, the matrix is surrounded by six fractures, with three groups of parallel fractures, which have a spacing L. For the 3D matrix-fracture system, the pressure diffusion behavior from matrix to fracture can be approximately regarded as a sphere (Fig. 5). The equivalent radius R of matrix is:L3=43πR3R=L30.75π.$$ {L}^3=\frac{4}{3}\pi {R}^3\Rightarrow R={L}^3\sqrt{\frac{0.75}{\pi }}. $$(27)

thumbnail Fig. 5

Schematic of 3D matrix-fracture system.

The 3D matrix fractal pressure diffusion equation can be given by (Kong et al., 2007):2Pr2+3-DTrPr=ϕaoμcKa(rλmax)2(DT-1)Pt.$$ \frac{{\mathrm{\partial }}^2P}{\mathrm{\partial }{r}^2}+\frac{3-{D}_T}{r}\frac{\mathrm{\partial }P}{\mathrm{\partial }r}=\frac{{\phi }_{{ao}}\mu c}{{K}_a}{\left(\frac{r}{{\lambda }_{\mathrm{max}}}\right)}^{2\left({D}_T-1\right)}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}. $$(28)

Initial and boundary conditions are given by:P=P0,0rR, t=0P=Pf, r=R, t>0.$$ \begin{array}{c}P={P}_0,\hspace{1em}0\le r\le R,\enspace \hspace{1em}t=0\\ P={P}_f,\enspace \hspace{1em}r=R,\enspace \hspace{1em}t>0.\end{array} $$(29)

Using separation variable method, the analytical solution of equation (28) can be obtained in the form of Bessel function (see Appendix B for details):P=Pf+rDT-22m=12(Pf-P0)ϖexp[-1M(DTI2RDT)2t],$$ P={P}_f+{r}^{\frac{{D}_T-2}{2}}\sum_{m=1}^{\mathrm{\infty }} 2\left({P}_f-{P}_0\right)\mathrm{\varpi exp}\left[-\frac{1}{M}{\left(\frac{{D}_T{I}_2}{{R}^{{D}_T}}\right)}^2t\right], $$(30)where ϖ=1Jv-1(I2)DT-|DT-2|2DT(RDTDT)-|DT-2|+2DT2DT(DTI2RDT)|DT-2|-4DT2DT×Jv(rDTI2RDT)$ \mathrm{\varpi }=\frac{1}{{J}_{v-1}\left({I}_2\right)}{{D}_T}^{-\frac{\left|{D}_T-2\right|}{2{D}_T}}{\left(\frac{{R}^{{D}_T}}{{D}_T}\right)}^{-\frac{\left|{D}_T-2\right|+2{D}_T}{2{D}_T}}{\left(\frac{{D}_T{I}_2}{{R}^{{D}_T}}\right)}^{\frac{\left|{D}_T-2\right|-4{D}_T}{2{D}_T}}\times {J}_v\left(\frac{{r}^{{D}_T}{I}_2}{{R}^{{D}_T}}\right)$, I2 is the positive zero point of Bessel function (the solution code of I2 is shown in Appendix C), it depends on the value of v, and v=|DT-2|2DT$ v=\frac{\left|{D}_T-2\right|}{2{D}_T}$.

Assuming that the volume of matrix Vav is equal to 43πr3$ \frac{4}{3}\pi {r}^3$, the average pressure of matrix is as follows:P̅=143πR30RPdVav.$$ \bar{P}=\frac{1}{\frac{4}{3}\pi {R}^3}{\int }_0^R P\mathrm{d}{V}_{{av}.} $$(31)

Substituting equation (30) into equation (31) yields:P̅=Pf+6R3m=1(Pf-P0)0RrDT+22ϖdrexp[-1M(DTI2RDT)2t].$$ \bar{P}={P}_f+\frac{6}{{R}^3}\sum_{m=1}^{\mathrm{\infty }} \left({P}_f-{P}_0\right){\int }_0^R {r}^{\frac{{D}_T+2}{2}}\mathrm{\varpi d}r\mathrm{exp}\left[-\frac{1}{M}{\left(\frac{{D}_T{I}_2}{{R}^{{D}_T}}\right)}^2t\right]. $$(32)

Both sides of equation (32) are subtracted by P0 and then divided by Pf − P0. Thus, equation (32) can be expressed as:P̅-P0Pf-Po=1+6R3m=10RrDT+22ϖdrexp[-1M(DTI2RDT)2t].$$ \frac{\bar{P}-{P}_0}{{P}_f-{P}_o}=1+\frac{6}{{R}^3}\sum_{m=1}^{\mathrm{\infty }} {\int }_0^R {r}^{\frac{{D}_T+2}{2}}\mathrm{\varpi d}r\mathrm{exp}\left[-\frac{1}{M}{\left(\frac{{D}_T{I}_2}{{R}^{{D}_T}}\right)}^2t\right]. $$(33)

The partial derivative of equation (33) with respect to t can be expressed as:P̅t1Pf-P0=-1M(DTI2RDT)26R3m=10RrDT+22ϖdrexp[-1M(DTI2RDT)2t]=-1M(DTI2RDT)2(P̅-P0Pf-P0-1).$$ \frac{\mathrm{\partial }\bar{P}}{\mathrm{\partial }t}\frac{1}{{P}_f-{P}_0}=-\frac{1}{M}{\left(\frac{{D}_T{I}_2}{{R}^{{D}_T}}\right)}^2\frac{6}{{R}^3}\sum_{m=1}^{\mathrm{\infty }} {\int }_0^R {r}^{\frac{{D}_T+2}{2}}\mathrm{\varpi d}r\mathrm{exp}\left[-\frac{1}{M}{\left(\frac{{D}_T{I}_2}{{R}^{{D}_T}}\right)}^2t\right]=-\frac{1}{M}{\left(\frac{{D}_T{I}_2}{{R}^{{D}_T}}\right)}^2\left(\frac{\bar{P}-{P}_0}{{P}_f-{P}_0}-1\right). $$(34)

Equation (34) can be also written as:P̅t=-Kaλmax2(DT-1)ϕaoμc(DTI2RDT)2(P̅-Pf).$$ \frac{\mathrm{\partial }\bar{P}}{\mathrm{\partial }t}=-\frac{{K}_a{\lambda }_{\mathrm{max}}^{2\left({D}_T-1\right)}}{{\phi }_{{ao}}\mu c}{\left(\frac{{D}_T{I}_2}{{R}^{{D}_T}}\right)}^2\left(\bar{P}-{P}_f\right). $$(35)

The 3D transfer function can be expressed as:q=-(ρϕ)mt=-ρϕaocP̅t(λmaxR)1-DT.$$ q=-\frac{\mathrm{\partial }(\rho \phi {)}_m}{\mathrm{\partial }t}=-\rho {\phi }_{{ao}}c\frac{\mathrm{\partial }\bar{P}}{\mathrm{\partial }t}{\left(\frac{{\lambda }_{\mathrm{max}}}{R}\right)}^{1-{D}_T}. $$(36)

Substituting equation (35) into equation (36) yields:q=ρμKaDT2I22R1+DTλmaxDT-1(P̅-Pf).$$ q=\frac{\rho }{\mu }{K}_a\frac{{D}_T^2{I}_2^2}{{R}^{1+{D}_T}}{\lambda }_{\mathrm{max}}^{{D}_T-1}\left(\bar{P}-{P}_f\right). $$(37)

Consequently, the 3D shape factor can be expressed as:σ=DT2I22R1+DTλmaxDT-1.$$ \sigma =\frac{{D}_T^2{I}_2^2}{{R}^{1+{D}_T}}{\lambda }_{\mathrm{max}}^{{D}_T-1}. $$(38)

Substituting equation (27) into equation (38) gives:σ=(43π)1+DT3DT2I22L1+DTλmaxDT-1.$$ \sigma ={\left(\frac{4}{3}\pi \right)}^{\frac{1+{D}_T}{3}}\frac{{D}_T^2{I}_2^2}{{L}^{1+{D}_T}}{\lambda }_{\mathrm{max}}^{{D}_T-1}. $$(39)

Compared with those shape factors proposed by predecessors, the fractal shape factor accounts for the influences of two additional parameters: tortuosity fractal dimension and maximum pore diameter of matrix. The effect of microstructure on shape factor directly affects fluids transfer in matrix-fracture system. When the flow channels in matrix are approximately straight, tortuosity fractal dimension of matrix DT = 1, our shape factor model has the similar value with Lim and Aziz (1995) model. When the flow channels in matrix are curved, according to the experimental observation, Yin et al. (2017) presented that tortuosity fractal dimension of sandstone, coal and shale media is generally ranging from 1.1 to 2.3. Nelson (2009) reported that pore diameters are usually greater than 2 μm in conventional reservoirs, ranging from 2 to 0.03 μm and from 0.1 to 0.005 μm in tight-gas sandstone and shale, respectively. Within the reasonable value range DT = 1.1 and λmax = 0.15 mm to calculate the positive zero point of Bessel function, we can get I1 = 2.474 and I2 = 3.013, and shape factors are 11.35/L2.1, 20.38/L2.1 and 24.77/L2.1 for 1D, 2D and 3D spaces, respectively. The shape factors predicted by Warren and Root (1963), Kazemi et al. (1976), Lim and Aziz (1995) and our model are presented in Table 1.

Table 1

The shape factors defined by various models.

3 Model verification

Here, the commercial reservoir simulator ECLIPSE is applied instead of the laboratory experiment to verify the proposed fractal shape factor. Three types of single-porosity fine-grid model are developed using ECLIPSE, and offer reference curve for the dual-porosity models with different shape factors. In the single-porosity fine-grid models, matrix is discretized that grids near fractures are finer than others, and fracture spacing is 10 cm (Fig. 6). The dual-porosity models only have one cubic grid with a size of 50 cm, which is mainly controlled by shape factor in transfer function, and fracture spacing is 50 cm. The basic parameters used in these models are shown in Table 2.

thumbnail Fig. 6

Schematic of single-porosity fine-grid models for (a) 1D, (b) 2D, and (c) 3D flow.

Table 2

Basic parameters used in numerical simulation.

Curves of cumulative production over time under 1D, 2D and 3D flows are displayed in Figures 7a7c, respectively. The results of single-porosity fine-grid models are taken as the reference curves. Results show that the fractal shape factors have relatively accurate matching degree with the reference curves. Clearly, the curvature of the fluid channels can’t be ignored in the real rocks. Both tortuosity fractal dimension and maximum pore diameter have an important role in the fluid flow in porous media (Yun et al., 2009; Ye et al., 2019; Huang et al., 2020). Thus, there will be a great error in the fluid exchange term if neglecting the effects of tortuosity fractal dimension and maximum pore diameter of matrix on the shape factor. The proposed fractal shape factor can be used to explain the impacts of tortuosity fractal dimension and maximum pore diameter on the fluid exchange rate, and give a better prediction of production rates and recoveries.

thumbnail Fig. 7

Cumulative production curves of dual-porosity model with different shape factors and single-porosity fine-grid model, (a) 1D flow; (b) 2D flow; (c) 3D flow.

4 Results and discussion

The fractal shape factor model presented in this paper (Eqs. (13), (26) and (39)) is in term of tortuosity fractal dimension, maximum pore diameter and characteristic length of matrix. The impacts of parameters on shape factor are studied in this part to understand the sensitivity to shape factor. It is generally believed that the permeability of low-permeable porous media is less than 10 × 10−3 μm2. With the decrease of permeability, the characteristics of nonlinear flow become obvious. Since the derivation of the presented fractal shape factor model follows Darcy’s law, so according to the experimental data (Yin et al., 2017), tortuosity fractal dimension of matrix is roughly less than 1.5. Six values of tortuosity fractal dimension of matrix are selected as 1.0, 1.1, 1.2, 1.3, 1.4 and 1.5, respectively, which reflect the different curvature of flow channels in porous media. A range of maximum pore diameters of matrix that varies from 0.1 μm to 2 mm, covering conventional rocks, tight sandstone and shale. The fractal shape factor model contains the positive zero point of Bessel function, which is related to tortuosity fractal dimension. So, we first solve the positive zero point of the Bessel function at different values of tortuosity fractal dimension, as shown in Table 3.

Table 3

The positive zero point of Bessel function at different values of tortuosity fractal dimension.

The values of the fractal shape factor decreases as tortuosity fractal dimension increases (Fig. 8). Since the larger tortuosity fractal dimension is, the more curved the flow channel is, which will lead to greater flow resistance, thus reducing the fluid exchange mass between matrix and fracture. However, when tortuosity fractal dimension is larger than 1.25, the shape factor of 2D flow is greater than that of 3D flow. Since the number of fractures around matrix increases, a larger flow cross-sectional area is formed, resulting in a greater fluid exchange rate between matrix and fracture. The shape factor of 2D flow should be less than that of 3D flow. Figure 8 indicates the presented fractal shape factor model is suitable for tortuosity fractal dimension less than 1.25.

thumbnail Fig. 8

The impact of tortuosity fractal dimension of matrix on the fractal shape factor.

The fractal shape factor growths with the increaseof maximum pore diameter of matrix (Fig. 9). The presented fractal shape factor model is applicable in all range of maximum pore diameter of matrix, indicating that the presented fractal shape factor model is suitable for dual-porosity media at different scales and the presented fractal shape factor varies at different scales.

thumbnail Fig. 9

The impact of maximum pore diameter of matrix on the fractal shape factor.

The fractal shape factor decreases with the increase of characteristic length of matrix (Fig. 10). The increase of characteristic length of matrix means that the fluid flow time is longer, reducing the rate of fluids transfer between matrix and fracture.

thumbnail Fig. 10

The impact of characteristic length of matrix on the fractal shape factor.

In order to further understand the fluids transfer mechanism between matrix and fracture, different shape factor models have been presented in the literature. Warren and Root (1963) combined an integral material balance to derive shape factors are 12/L2, 32/L2 and 60/L2 for 1D, 2D and 3D, respectively. Moreover, Lim and Aziz (1995) used approximating analytical solution of matrix pressure diffusion equation to derived shape factors are π2/L2, 18.17/L2 (or 2π2/L2) and 25.67/L2 (or 3π2/L2) for 1D, 2D and 3D, respectively. He et al. (2017) considered the influence of tortuosity (τ) and introduced the modified shape factors are π2/(τL2), 18/(τL2) and 26/(τL2) for 1D, 2D and 3D, respectively. When these shape factors are applied in the transfer function to investigate the fluids transfer between matrix and fracture, the impact of matrix microscopic pore structure is usually neglected, which is deviated from the actual situation. The presented fractal shape factor considers the impacts of tortuosity fractal dimension of matrix, maximum pore diameter and characteristic length. Compared with the previous models, the proposed fractal shape factor model can more accurately characterize the impact of microscopic pore structure on the rate of fluids transfer between matrix and fracture in fractured reservoirs.

5 Conclusion

We established new transfer function and shape factor to investigate the fluids transfer between matrix and fracture in fractured reservoirs. The analytical expression for fractal shape factors in 1D, 2D and 3D spaces are derived via the fractal pressure diffusion equation. The presented fractal shape factor model is authenticated by comparing dual-porosity model with single-porosity fine-grid model. Furtherly, sensitivity analysis of model parameters is carried out to understand their effects on fractal shape factor. Some conclusions can be obtained:

  1. For different porous media, the fractal shape factor can provide corresponding values. Applying the fractal shape factor to transfer function will help us better understand the impact of microscopic pore structure on fluid exchange law between matrix and fracture.

  2. There is a positive correlation between the fractal shape factor and maximum pore diameter of matrix. The fractal shape factor decreases with the growth of tortuosity fractal dimension and characteristic length .

  3. When tortuosity fractal dimension is roughly less than 1.25, the fractal shape factor has practical application value in different scales.

Acknowledgments

This work is supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (No. XDA14010302), the National Natural Science Foundation of China (No. 51904279), and the Fundamental Research Funds for the Central Universities (China University of Geosciences, Wuhan) (No. CUGGC04).

Appendix A

Derivation of 1D fractal pressure diffusion equation

Supposing matrix has been surrounded by two fractures (Fig. 3), equation (A.1) for 1D matrix fractal pressure diffusion equation (Kong et al., 2007):2Px2+1-DTxPx=ϕaoμcKa(xλmax)2(DT-1)Pt.$$ \frac{{\mathrm{\partial }}^2P}{\mathrm{\partial }{x}^2}+\frac{1-{D}_T}{x}\frac{\mathrm{\partial }P}{\mathrm{\partial }x}=\frac{{\phi }_{{ao}}\mu c}{{K}_a}{\left(\frac{x}{{\lambda }_{\mathrm{max}}}\right)}^{2\left({D}_T-1\right)}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}. $$(A.1)

Let W = P − Pf, and initial and boundary conditions can be rewritten as:W=P0-Pf,-L/2xL/2, t=0W=0, x=±L/2, t>0.$$ \begin{array}{c}W={P}_0-{P}_f,\hspace{1em}-L/2\le x\le L/2,\enspace \hspace{1em}t=0\\ W=0,\enspace \hspace{1em}x=\pm L/2,\hspace{1em}\enspace t>0.\end{array} $$(A.2)

Substituting equation (A.2) into equation (A.1), it yields:2Wx2+1-DTxWx=ϕaoμcKa(xλmax)2(DT-1)Wt.$$ \frac{{\mathrm{\partial }}^2W}{\mathrm{\partial }{x}^2}+\frac{1-{D}_T}{x}\frac{\mathrm{\partial }W}{\mathrm{\partial }x}=\frac{{\phi }_{{ao}}\mu c}{{K}_a}{\left(\frac{x}{{\lambda }_{\mathrm{max}}}\right)}^{2\left({D}_T-1\right)}\frac{\mathrm{\partial }W}{\mathrm{\partial }t}. $$(A.3)

Using separation variable method, it must have:W=u(x)exp(-1Mα2t),$$ W=u(x)\mathrm{exp}\left(-\frac{1}{M}{\alpha }^2t\right), $$(A.4)where M=ϕaoμcKaλmax2(DT-1)$ M=\frac{{\phi }_{{ao}}\mu c}{{K}_a{\lambda }_{\mathrm{max}}^{2\left({D}_T-1\right)}}$.

Substituting equation (A.4) into equation (A.3), we can obtain:2ux2+1-DTxux+α2x2(DT-1)u=0.$$ \frac{{\mathrm{\partial }}^2u}{\mathrm{\partial }{x}^2}+\frac{1-{D}_T}{x}\frac{\mathrm{\partial }u}{\mathrm{\partial }x}+{\alpha }^2{x}^{2\left({D}_T-1\right)}u=0. $$(A.5)

According to the solution method of Bessel function, the solution of equation (A.5) can be obtained:μ=m=1[Amsin(αmxDTDT)+Bmcos(αmxDTDT)].$$ \mu =\sum_{m=1}^{\mathrm{\infty }} \left[{A}_m\mathrm{sin}\left(\frac{{\alpha }_m{x}^{{D}_T}}{{D}_T}\right)+{B}_m\mathrm{cos}\left(\frac{{\alpha }_m{x}^{{D}_T}}{{D}_T}\right)\right]. $$(A.6)

Substituting equation (A.6) into equation (A.4), it yields:W=m=1[Amsin(αmxDTDT)+Bmcos(αmxDTDT)]exp(-1Mαm2t).$$ W=\sum_{m=1}^{\mathrm{\infty }} \left[{A}_m\mathrm{sin}\left(\frac{{\alpha }_m{x}^{{D}_T}}{{D}_T}\right)+{B}_m\mathrm{cos}\left(\frac{{\alpha }_m{x}^{{D}_T}}{{D}_T}\right)\right]\mathrm{exp}\left(-\frac{1}{M}{\alpha }_m^2t\right). $$(A.7)

Combining equations (A.2) with (A.7), we can obtain:m=1[Amsin(αmLDT2DTDT)+Bmcos(αmLDT2DTDT)]exp(-1Mαm2t)=0,$$ \sum_{m=1}^{\mathrm{\infty }} \left[{A}_m\mathrm{sin}\left(\frac{{\alpha }_m{L}^{{D}_T}}{{2}^{{D}_T}{D}_T}\right)+{B}_m\mathrm{cos}\left(\frac{{\alpha }_m{L}^{{D}_T}}{{2}^{{D}_T}{D}_T}\right)\right]\mathrm{exp}\left(-\frac{1}{M}{\alpha }_m^2t\right)=0, $$(A.8a)m=1[-Amsin(αmLDT2DTDT)+Bmcos(αmLDT2DTDT)]exp(-1Mαm2t)=0,$$ \sum_{m=1}^{\mathrm{\infty }} \left[-{A}_m\mathrm{sin}\left(\frac{{\alpha }_m{L}^{{D}_T}}{{2}^{{D}_T}{D}_T}\right)+{B}_m\mathrm{cos}\left(\frac{{\alpha }_m{L}^{{D}_T}}{{2}^{{D}_T}{D}_T}\right)\right]\mathrm{exp}\left(-\frac{1}{M}{\alpha }_m^2t\right)=0, $$(A.8b)m=1[Amsin(αmxDTDT)+Bmcos(αmxDTDT)]=P0-Pf.$$ \sum_{m=1}^{\mathrm{\infty }} \left[{A}_m\mathrm{sin}\left(\frac{{\alpha }_m{x}^{{D}_T}}{{D}_T}\right)+{B}_m\mathrm{cos}\left(\frac{{\alpha }_m{x}^{{D}_T}}{{D}_T}\right)\right]={P}_0-{P}_f. $$(A.8c)

Consequently, the solutions of equations (A.8a) and (A.8b) can be expressed as:Am=0, αm=2DT-1(2n+1)πDTLDT.$$ {A}_m=0,\enspace {\alpha }_m=\frac{{2}^{{D}_T-1}(2n+1)\pi {D}_T}{{L}^{{D}_T}}. $$(A.9)

Substituting equation (A.9) into equation (A.8c), equation (A.8c) can be rewritten as:P0-Pf=m=1Bmcos(αmxDTDT).$$ {P}_0-{P}_f=\sum_{m=1}^{\mathrm{\infty }} {B}_m\mathrm{cos}\left(\frac{{\alpha }_m{x}^{{D}_T}}{{D}_T}\right). $$(A.10)

Both sides of equation (A.10) are multiplied by xDT-1cos(αmxDTDT)$ {x}^{{D}_T-1}\mathrm{cos}\left(\frac{{\alpha }_m{x}^{{D}_T}}{{D}_T}\right)$ and integrated with respect to x. Then, the value of Bm can be obtained as:Bm=1-L/2L/2xDT-1cos2(αmxDTDT)dx-L/2L/2(P0-Pf)xDT-1cos(αmxDTDT)dx=(-1)n4(P0-Pf)(2n+1)π.$$ {B}_m=\frac{1}{{\int }_{-L/2}^{L/2} {x}^{{D}_T-1}\mathrm{co}{\mathrm{s}}^2\left(\frac{{\alpha }_m{x}^{{D}_T}}{{D}_T}\right)\mathrm{d}x}{\int }_{-L/2}^{L/2} \left({P}_0-{P}_f\right){x}^{{D}_T-1}\mathrm{cos}\left(\frac{{\alpha }_m{x}^{{D}_T}}{{D}_T}\right)\mathrm{d}x={\left(-1\right)}^n\frac{4\left({P}_0-{P}_f\right)}{\left(2n+1\right)\pi }. $$(A.11)

Substituting equations (A.9) and (A.11) into equation (A.7), thus, equation (A.7) becomes:W=P-Pf=n=0(-1)n4(P0-Pf)(2n+1)πcos[(2n+1)π2xDT(L/2)DT]exp[-1M22DT-2(2n+1)2π2DT2L2DTt].$$ W=P-{P}_f=\sum_{n=0}^{\mathrm{\infty }} {\left(-1\right)}^n\frac{4\left({P}_0-{P}_f\right)}{\left(2n+1\right)\pi }\mathrm{cos}\left[\frac{(2n+1)\pi }{2}\frac{{x}^{{D}_T}}{{\left(L/2\right)}^{{D}_T}}\right]\mathrm{exp}\left[-\frac{1}{M}\frac{{2}^{2{D}_T-2}(2n+1{)}^2{\pi }^2{D}_T^2}{{L}^{2{D}_T}}t\right]. $$(A.12)

Equation (A.12) can be converted to:P=Pf+n=0(-1)n4(P0-Pf)(2n+1)πcos[(2n+1)π2xDT(L/2)DT]exp[-1M22DT-2(2n+1)2π2DT2L2DTt].$$ P={P}_f+\sum_{n=0}^{\mathrm{\infty }} {\left(-1\right)}^n\frac{4\left({P}_0-{P}_f\right)}{\left(2n+1\right)\pi }\mathrm{cos}\left[\frac{(2n+1)\pi }{2}\frac{{x}^{{D}_T}}{{\left(L/2\right)}^{{D}_T}}\right]\mathrm{exp}\left[-\frac{1}{M}\frac{{2}^{2{D}_T-2}(2n+1{)}^2{\pi }^2{D}_T^2}{{L}^{2{D}_T}}t\right]. $$(A.13)

Assuming that the area of matrix Aav is equal to hx, h is the height which is constant, and x is distance. The average pressure of the matrix is as follows:P̅=1Aav-L/2L/2PdAav.$$ \bar{P}=\frac{1}{{A}_{{av}}}{\int }_{-L/2}^{L/2} P\mathrm{d}{A}_{{av}}. $$(A.14)

Substituting equation (A.13) into equation (A.14), it yields:P̅=Pf+n=0(-1)n8(P0-Pf)Lπ2(2n+1)2×δexp[-1M22DT-2(2n+1)2π2DT2L2DTt],$$ \bar{P}={P}_f+\sum_{n=0}^{\mathrm{\infty }} \frac{{\left(-1\right)}^n8\left({P}_0-{P}_f\right)}{L{\pi }^2{\left(2n+1\right)}^2}\times \delta \mathrm{exp}\left[-\frac{1}{M}\frac{{2}^{2{D}_T-2}(2n+1{)}^2{\pi }^2{D}_T^2}{{L}^{2{D}_T}}t\right], $$(A.15)where δ=-L/2L/2cos[(2n+1)π2xDT(L/2)DT]1DT(xL/2)1-DTd[(2n+1)π2xDT(L/2)DT]$ \delta ={\int }_{-L/2}^{L/2} \mathrm{cos}\left[\frac{(2n+1)\pi }{2}\frac{{x}^{{D}_T}}{{\left(L/2\right)}^{{D}_T}}\right]\frac{1}{{D}_T}{\left(\frac{x}{L/2}\right)}^{1-{D}_T}d\left[\frac{(2n+1)\pi }{2}\frac{{x}^{{D}_T}}{{\left(L/2\right)}^{{D}_T}}\right]$.

Both sides of equation (A.15) are subtracted by P0 and then divided by Pf − P0. Thus, equation (A.15) can be expressed as:P̅-P0Pf-Po=1-n=0(-1)n8×e(2n+1)2Lπ2(2n+1)2×δexp(-1M22DT-2π2DT2L2DTt).$$ \frac{\bar{P}-{P}_0}{{P}_f-{P}_o}=1-\sum_{n=0}^{\mathrm{\infty }} \frac{{\left(-1\right)}^n8\times {\mathrm{e}}^{(2n+1{)}^2}}{L{\pi }^2(2n+1{)}^2}\times \delta \mathrm{exp}\left(-\frac{1}{M}\frac{{2}^{2{D}_T-2}{\pi }^2{D}_T^2}{{L}^{2{D}_T}}t\right). $$(A.16)

Appendix B

Derivation of 2D and 3D fractal pressure diffusion equation

The matrix fractal pressure diffusion equation is given by (Kong et al., 2007):2Pr2+ξ-DTrPr=ϕaoμcKa(rλmax)2(DT-1)Pt,$$ \frac{{\mathrm{\partial }}^2P}{\mathrm{\partial }{r}^2}+\frac{\xi -{D}_T}{r}\frac{\mathrm{\partial }P}{\mathrm{\partial }r}=\frac{{\phi }_{{ao}}\mu c}{{K}_a}{\left(\frac{r}{{\lambda }_{\mathrm{max}}}\right)}^{2\left({D}_T-1\right)}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}, $$(B.1)where ξ corresponds to 2D flow and 3D flow, which are 2 and 3, respectively.

Let G = P − Pf, and initial and boundary conditions can be rewritten as:G=P0-Pf,0rR, t=0G=0, r=R, t>0.$$ \begin{array}{c}G={P}_0-{P}_f,\hspace{1em}0\le r\le R,\hspace{1em}\enspace t=0\\ G=0,\hspace{1em}\enspace r=R,\hspace{1em}\enspace t>0.\end{array} $$(B.2)

Substituting equation (B.2) into equation (B.1), it yields:2Gr2+ξ-DTrGr=ϕaoμcKa(rλmax)2(DT-1)Gt.$$ \frac{{\mathrm{\partial }}^2G}{\mathrm{\partial }{r}^2}+\frac{\xi -{D}_T}{r}\frac{\mathrm{\partial }G}{\mathrm{\partial }r}=\frac{{\phi }_{{ao}}\mu c}{{K}_a}{\left(\frac{r}{{\lambda }_{\mathrm{max}}}\right)}^{2\left({D}_T-1\right)}\frac{\mathrm{\partial }G}{\mathrm{\partial }t}. $$(B.3)

Using separation variable method, we can see that:G=w(r)exp(-1Mβ2t).$$ G=w(r)\mathrm{exp}\left(-\frac{1}{M}{\beta }^2t\right). $$(B.4)

Substituting equation (B.4) into equation (B.3), we can obtain:2wr2+ξ-DTrwr=β2r2(DT-1)wt.$$ \frac{{\mathrm{\partial }}^2w}{\mathrm{\partial }{r}^2}+\frac{\xi -{D}_T}{r}\frac{\mathrm{\partial }w}{\mathrm{\partial }r}={\beta }^2{r}^{2\left({D}_T-1\right)}\frac{\mathrm{\partial }w}{\mathrm{\partial }t}. $$(B.5)

According to the solution method of Bessel function, the solution of equation (B.5) can be obtained:w=rDT+1-ξ2m=1[cmJv(βmrDTDT)+bmcos(vπ)Jv(βmrDTDT)-J-v(βmrDTDT)sin(vπ)],$$ w={r}^{\frac{{D}_T+1-\xi }{2}}\sum_{m=1}^{\mathrm{\infty }} \left[{c}_m{J}_v\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)+{b}_m\frac{\mathrm{cos}\left(v\pi \right){J}_v\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)-{J}_{-v}\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)}{\mathrm{sin}\left(v\pi \right)}\right], $$(B.6)where v=|DT+1-ξ|2DT$ v=\frac{\left|{D}_T+1-\xi \right|}{2{D}_T}$, v corresponds to 2D flow and 3D flow, which are v2=DT-12DT$ {v}_2=\frac{{D}_T-1}{2{D}_T}$ and v3=|DT-2|2DT$ {v}_3=\frac{\left|{D}_T-2\right|}{2{D}_T}$, respectively.

Substituting equation (B.6) into equation (B.4), it yields:G=rDT+1-ξ2m=1[cmJv(βmrDTDT)+bmcos(vπ)Jv(βmrDTDT)-J-v(βmrDTDT)sin(vπ)]exp(-1Mβm2t).$$ G={r}^{\frac{{D}_T+1-\xi }{2}}\sum_{m=1}^{\mathrm{\infty }} \left[{c}_m{J}_v\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)+{b}_m\frac{\mathrm{cos}\left(v\pi \right){J}_v\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)-{J}_{-v}\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)}{\mathrm{sin}\left(v\pi \right)}\right]\mathrm{exp}\left(-\frac{1}{M}{\beta }_m^2t\right). $$(B.7)

Combining equations (B.2) and (B.7), we can obtain:RDT+1-ξ2m=1[cmJv(βmRDTDT)+bmcos(vπ)Jv(βmRDTDT)-J-v(βmRDTDT)sin(vπ)]exp(-1Mβm2t)=0,$$ {R}^{\frac{{D}_T+1-\xi }{2}}\sum_{m=1}^{\mathrm{\infty }} \left[{c}_m{J}_v\left(\frac{{\beta }_m{R}^{{D}_T}}{{D}_T}\right)+{b}_m\frac{\mathrm{cos}\left(v\pi \right){J}_v\left(\frac{{\beta }_m{R}^{{D}_T}}{{D}_T}\right)-{J}_{-v}\left(\frac{{\beta }_m{R}^{{D}_T}}{{D}_T}\right)}{\mathrm{sin}\left(v\pi \right)}\right]\mathrm{exp}\left(-\frac{1}{M}{\beta }_m^2t\right)=0, $$(B.8a)rDT+1-ξ2m=1[cmJv(βmrDTDT)+bmcos(vπ)Jv(βmrDTDT)-J-v(βmrDTDT)sin(vπ)]=P0-Pf.$$ {r}^{\frac{{D}_T+1-\xi }{2}}\sum_{m=1}^{\mathrm{\infty }} \left[{c}_m{J}_v\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)+{b}_m\frac{\mathrm{cos}\left(v\pi \right){J}_v\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)-{J}_{-v}\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)}{\mathrm{sin}\left(v\pi \right)}\right]={P}_0-{P}_f. $$(B.8b)

Equation (B.8b) can be converted to:P0-Pf=rDT+1-ξ2m=1Jv(βmrDTDT)[cm+bmcos(vπ)sin(vπ)-bmJ-v(βmrDTDT)Jv(βmrDTDT)]=rDT+1-ξ2m=1Jv(βmrDTDT)[cm+bmcos(vπ)sin(vπ)-bm(βmrDT2DT)-2vΓ(m+1+v)Γ(m+1-v)].$$ {P}_0-{P}_f={r}^{\frac{{D}_T+1-\xi }{2}}\sum_{m=1}^{\mathrm{\infty }} {J}_v\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)\left[{c}_m+{b}_m\frac{\mathrm{cos}\left(v\pi \right)}{\mathrm{sin}\left(v\pi \right)}-{b}_m\frac{{J}_{-v}\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)}{{J}_v\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)}\right]={r}^{\frac{{D}_T+1-\xi }{2}}\sum_{m=1}^{\mathrm{\infty }} {J}_v\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)\left[{c}_m+{b}_m\frac{\mathrm{cos}\left(v\pi \right)}{\mathrm{sin}\left(v\pi \right)}-{b}_m{\left(\frac{{\beta }_m{r}^{{D}_T}}{2{D}_T}\right)}^{-2v}\frac{\mathrm{\Gamma }\left(m+1+v\right)}{\mathrm{\Gamma }\left(m+1-v\right)}\right]. $$(B.9)

In equation (B.9), when r → 0, we have (βmrDT2DT)-2v$ {\left(\frac{{\beta }_m{r}^{{D}_T}}{2{D}_T}\right)}^{-2v}\to \mathrm{\infty }$, but P0 − Pf << ∞, so we can obtain:bm=0.$$ {b}_m=0. $$(B.10)

Then, equation (B.9) can be reduced to:P0-Pf=rDT+1-ξ2m=1cmJv(βmrDTDT).$$ {P}_0-{P}_f={r}^{\frac{{D}_T+1-\xi }{2}}\sum_{m=1}^{\mathrm{\infty }} {c}_m{J}_v\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right). $$(B.11)

Substituting equation (B.10) into equation (B.8a) yields:RDT+1-ξ2m=1[cmJv(βmRDTDT)]exp(-1Mβm2t)=0.$$ {R}^{\frac{{D}_T+1-\xi }{2}}\sum_{m=1}^{\mathrm{\infty }} \left[{c}_m{J}_v\left(\frac{{\beta }_m{R}^{{D}_T}}{{D}_T}\right)\right]\mathrm{exp}\left(-\frac{1}{M}{\beta }_m^2t\right)=0. $$(B.12)

In equation (B.12), due to RDT+1-ξ20$ {R}^{\frac{{D}_T+1-\xi }{2}}\ne 0$ and exp(-1Mβm2t)0$ \mathrm{exp}\left(-\frac{1}{M}{\beta }_m^2t\right)\ne 0$, so we can obtain:Jv(βmRDTDT)=Jv(Im)=0soβm=DTImRDT,$$ \begin{array}{c}{J}_v\left(\frac{{\beta }_m{R}^{{D}_T}}{{D}_T}\right)={J}_v\left({I}_m\right)=0\\ \mathrm{so}\hspace{1em}{\beta }_m=\frac{{D}_T{I}_m}{{R}^{{D}_T}},\end{array} $$(B.13)where Jv (x) is the Bessel function, v is the order of Bessel function, Im is the positive zero point of Bessel function (the solution code of Im is shown in Appendix C), it depends on the value of v, Im corresponds to 2D flow and 3D flow, which are I1 and I2, respectively.

Both sides of equation (B.12) are multiplied by rDTDTJv(βmrDTDT)$ \frac{{r}^{{D}_T}}{{D}_T}{J}_v\left(\frac{{\beta }_m{r}^{{D}_T}}{{D}_T}\right)$ and integrated with respect to x. Then, the value of cm can be obtained as:cm=2(P0-Pf)Jv-1(βmRDTDT)[-DT-v(RDTDT)-v-1βmv-2].$$ {c}_m=\frac{2\left({P}_0-{P}_f\right)}{{J}_{v-1}\left(\frac{{\beta }_m{R}^{{D}_T}}{{D}_T}\right)}\left[-{{D}_T}^{-v}{\left(\frac{{R}^{{D}_T}}{{D}_T}\right)}^{-v-1}{\beta }_m^{v-2}\right]. $$(B.14)

Substituting equations (B.10), (B.13) and (B.14) into equation (B.7), thus, equation (B.7) becomes:G=P-Pf=rDT+1-ξ2m=12(Pf-P0)Jv-1(Im)DT-v(RDTDT)-v-1×(DTImRDT)v-2Jv(rDTImRDT)exp[-1M(DTImRDT)2t].$$ G=P-{P}_f={r}^{\frac{{D}_T+1-\xi }{2}}\sum_{m=1}^{\mathrm{\infty }} \frac{2\left({P}_f-{P}_0\right)}{{J}_{v-1}\left({I}_m\right)}{{D}_T}^{-v}{\left(\frac{{R}^{{D}_T}}{{D}_T}\right)}^{-v-1}\times {\left(\frac{{D}_T{I}_m}{{R}^{{D}_T}}\right)}^{v-2}{J}_v\left(\frac{{r}^{{D}_T}{I}_m}{{R}^{{D}_T}}\right)\mathrm{exp}\left[-\frac{1}{M}{\left(\frac{{D}_T{I}_m}{{R}^{{D}_T}}\right)}^2t\right]. $$(B.15)

Equation (B.15) can be converted to:P=Pf+rDT+1-ξ2m=12(Pf-P0)Jv-1(Im)DT-v(RDTDT)-v-1(DTImRDT)v-2×Jv(rDTImRDT)exp[-1M(DTImRDT)2t]$$ P={P}_f+{r}^{\frac{{D}_T+1-\xi }{2}}\sum_{m=1}^{\mathrm{\infty }} \frac{2\left({P}_f-{P}_0\right)}{{J}_{v-1}\left({I}_m\right)}{{D}_T}^{-v}{\left(\frac{{R}^{{D}_T}}{{D}_T}\right)}^{-v-1}{\left(\frac{{D}_T{I}_m}{{R}^{{D}_T}}\right)}^{v-2}\times {J}_v\left(\frac{{r}^{{D}_T}{I}_m}{{R}^{{D}_T}}\right)\mathrm{exp}\left[-\frac{1}{M}{\left(\frac{{D}_T{I}_m}{{R}^{{D}_T}}\right)}^2t\right] $$(B.16)

Appendix C

Solution code of the positive zero of Bessel function

In this paper, we need to get the positive zero of Bessel function which is solved by MATLAB software. The specific code is as follows:

clear, clc;

format long

x=(0:0.2:1)′;

y_0=fzero(@(x)besselj(0.045,x),3);

References

  • Acuna J.A., Ershaghi I., Yortsos Y.C. (1995) Practical application of fractal pressure transient analysis of naturally fractured reservoirs, SPE Form. Eval. 10, 3, 173–179. [CrossRef] [Google Scholar]
  • Adler P.M. (1996) Transports in fractal porous media, J. Hydrol. 187, 195–213. [CrossRef] [Google Scholar]
  • Barenblatt G.I., Zheltov I.P., Kochina I.N. (1960) Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks, J. Appl. Math. Mech. 24, 5, 1286–1303. [CrossRef] [Google Scholar]
  • Behnoudfar P., Asadi M.B., Gholilou A., Zendehboudi S. (2019) A new model to conduct hydraulic fracture design in coalbed methane reservoirs by incorporating stress variations, J. Pet. Sci. Eng. 174, 1208–1222. [Google Scholar]
  • Cao R., Xu Z., Cheng L., Peng Y., Wang Y., Guo Z. (2019) Study of single phase mass transfer between matrix and fracture in tight oil reservoirs, Geofluids 2019, 1–11. [Google Scholar]
  • Chang J., Yortsos Y.C. (1990) Pressure-transient analysis of fractal reservoirs, SPE Form. Eval. 5, 1, 31–38. [CrossRef] [Google Scholar]
  • Coats K.H. (1989) Implicit compositional simulation of single-porosity and dual-porosity reservoirs, in: SPE Symposium on Reservoir Simulation, Houston, Texas, 6–8 February, Society of Petroleum Engineers. [Google Scholar]
  • Costa A. (2006) Permeability-porosity relationship: A reexamination of the Kozeny-Carman equation based on a fractal pore-space geometry assumption, Geophys. Res. Lett. 33, 2, L02318. [Google Scholar]
  • Erol S., Fowler S.J., Harcouët-Menou V., Laenen B. (2017) An analytical model of porosity–permeability for porous and fractured media, Transp. Porous Media 120, 2, 327–358. [Google Scholar]
  • Fan D., Ettehadtavakkol A. (2017) Semi-analytical modeling of shale gas flow through fractal induced fracture networks with microseismic data, Fuel 193, 444–459. [CrossRef] [Google Scholar]
  • Flamenco-Lopez F., Camacho-Velazquez R. (2001) Fractal transient pressure behavior of naturally fractured reservoirs, in: SPE Annual Technical Conference and Exhibition, 30 September–3 October, New Orleans, Louisiana, Society of Petroleum Engineers. [Google Scholar]
  • Hassanzadeh H., Pooladidarvish M. (2006) Effects of fracture boundary conditions on matrix-fracture transfer shape factor, Transp. Porous Media 64, 1, 51–71. [Google Scholar]
  • He Y., Chen X., Zhang Y., Yu W. (2017) Modeling interporosity flow functions and shape factors in low-permeability naturally fractured reservoir, J. Pet. Sci. Eng. 156, 110–117. [Google Scholar]
  • Huang T., Du P., Peng X., Wang P., Zou G. (2020) Pressure drop and fractal non-Darcy coefficient model for fluid flow through porous media, J. Pet. Sci. Eng. 184, 106579. [Google Scholar]
  • Katz A.J., Thompson A.H. (1985) Fractal sandstone pores: Implications for conductivity and pore formation, Phys. Rev. Lett. 5412, 1325–1328. [CrossRef] [PubMed] [Google Scholar]
  • Kazemi H., Merrill L.S., Porterfield K.L., Zeman P.R. (1976) Numerical simulation of water-oil flow in naturally fractured reservoirs, SPE J. 16, 6, 317–326. [Google Scholar]
  • Kong X., Li D., Lu D. (2007) Basic formulas of fractal seepage and type-curves of fractal reservoirs, Journal of Xi’an Shiyou University (Natural Science Edition) 22, 2, 1–10. [Google Scholar]
  • Kong X., Li D., Lu D. (2009) Transient pressure analysis in porous and fractured fractal reservoirs, Sci. China, Ser. E: Technol. Sci. 52, 9, 2700–2708. [CrossRef] [Google Scholar]
  • Krohn C.E. (1988) Fractal measurements of sandstones, shales, and carbonates, J. Geophys. Res. 93, 3297–3305. [Google Scholar]
  • Li Q., Xing H., Liu J., Liu X. (2015) A review on hydraulic fracturing of unconventional reservoir, Petroleum 1, 1, 8–15. [CrossRef] [Google Scholar]
  • Lian P.Q., Duan T.Z., Xu R., Li L.L., Li M. (2018) Pressure behavior of shale-gas flow in dual porous medium based on fractal theory, Interpretation 6, 4, SN1–SN10. [CrossRef] [Google Scholar]
  • Lim K.T., Aziz K. (1995) Matrix-fracture transfer shape factors for dual-porosity simulators, J. Pet. Sci. Eng. 13, 3, 169–178. [Google Scholar]
  • Nelson P.H. (2009) Pore-throat sizes in sandstones, tight sandstones, and shales, AAPG Bull. 93, 3, 329–340. [CrossRef] [Google Scholar]
  • Noetinger B., Estebenet T. (2000) Up-scaling of double porosity fractured media using continuous-time random walks methods, Transp. Porous Media 39, 3, 315–337. [Google Scholar]
  • Ranjbar E., Hassanzadeh H., Chen Z. (2011) Effect of fracture pressure depletion regimes on the dual-porosity shape factor for flow of compressible fluids in fractured porous media, Adv. Water Resour. 34, 12, 1681–1693. [Google Scholar]
  • Saboorian-Jooybari H., Ashoori S., Mowazi G. (2015) A new transient matrix/fracture shape factor for capillary and gravity imbibition in fractured reservoirs, Energy Sources Part A 37, 23, 2497–2506. [CrossRef] [Google Scholar]
  • Sarma P., Aziz K. (2004) New transfer functions for simulation of naturally fractured reservoirs with dual porosity models, in: SPE Annual Technical Conference and Exhibition, Houston, Texas, 26–29 September, Society of Petroleum Engineers. [Google Scholar]
  • Thomas L.K., Dixon T.N., Pierson R.G. (1983) Fractured reservoir simulation, SPE J. 23, 1, 42–54. [Google Scholar]
  • Ueda Y., Murata S., Watanabe Y., Funatsu K. (1989) Investigation of the shape factor used in the dual-porosity reservoir simulator, in: SPE Asia-Pacific Conference, 13–15 September, Sydney, Society of Petroleum Engineers. [Google Scholar]
  • Vishkai M., Gates I. (2019) On multistage hydraulic fracturing in tight gas reservoirs: Montney Formation, Alberta, Canada, J. Pet. Sci. Eng. 174, 1127–1141. [Google Scholar]
  • Wang W., Yuan B., Su Y., Sheng G., Yao W., Gao H., Wang K. (2018) A composite dual-porosity fractal model for channel-fractured horizontal wells, Eng. Appl. Comput. Fluid Mech. 12, 1, 104–116. [Google Scholar]
  • Wang W., Zheng D., Sheng G., Zhang Q., Su Y. (2017) A review of stimulated reservoir volume characterization for multiple fractured horizontal well in unconventional reservoirs, Adv. Geo-Energy Res. 1, 1, 54–63. [CrossRef] [Google Scholar]
  • Warren J.E., Root P.J. (1963) The behavior of naturally fractured reservoirs, SPE J. 3, 3, 245–255. [Google Scholar]
  • Yao Y., Wu Y., Zhang R. (2012) The transient flow analysis of fluid in a fractal, double-porosity reservoir, Transp. Porous Media 94, 1, 175–187. [Google Scholar]
  • Ye W., Wang X., Cao C., Yu W. (2019) A fractal model for threshold pressure gradient of tight oil reservoirs, J. Pet. Sci. Eng. 179, 427–431. [Google Scholar]
  • Yin S., Xie R., Ding W., Shan Y., Zhou W. (2017) Influences of fractal characteristics of reservoir rocks on permeability, Lithologic Reserv. 29, 4, 81–90. [Google Scholar]
  • Yu B., Cheng P. (2002) A fractal permeability model for bi-dispersed porous media, Int. J. Heat Mass Transf. 45, 14, 2983–2993. [Google Scholar]
  • Yun M., Yu B., Cai J. (2009) Analysis of seepage characters in fractal porous media, Int. J. Heat Mass Transf. 52, 13, 3272–3278. [Google Scholar]
  • Zhou D., Ge J., Li Y., Cai Y. (2000) Establishment of interporosity flow function of complex fractured reservoirs, OGRT 7, 2, 30–32. [Google Scholar]
  • Zimmerman R.W., Chen G., Hadgu T., Bodvarsson G.S. (1993) A numerical dual-porosity model with semianalytical treatment of fracture/matrix flow, Water Resour. Res. 29, 7, 2127–2137. [Google Scholar]

All Tables

Table 1

The shape factors defined by various models.

Table 2

Basic parameters used in numerical simulation.

Table 3

The positive zero point of Bessel function at different values of tortuosity fractal dimension.

All Figures

thumbnail Fig. 1

Schematic of idealization of fractured reservoirs (Warren and Root, 1963).

In the text
thumbnail Fig. 2

Schematic of dual-porosity media.

In the text
thumbnail Fig. 3

Schematic of 1D matrix-fracture system.

In the text
thumbnail Fig. 4

Schematic of 2D matrix-fracture system.

In the text
thumbnail Fig. 5

Schematic of 3D matrix-fracture system.

In the text
thumbnail Fig. 6

Schematic of single-porosity fine-grid models for (a) 1D, (b) 2D, and (c) 3D flow.

In the text
thumbnail Fig. 7

Cumulative production curves of dual-porosity model with different shape factors and single-porosity fine-grid model, (a) 1D flow; (b) 2D flow; (c) 3D flow.

In the text
thumbnail Fig. 8

The impact of tortuosity fractal dimension of matrix on the fractal shape factor.

In the text
thumbnail Fig. 9

The impact of maximum pore diameter of matrix on the fractal shape factor.

In the text
thumbnail Fig. 10

The impact of characteristic length of matrix on the fractal shape factor.

In the text

Les statistiques affichées correspondent au cumul d'une part des vues des résumés de l'article et d'autre part des vues et téléchargements de l'article plein-texte (PDF, Full-HTML, ePub... selon les formats disponibles) sur la platefome Vision4Press.

Les statistiques sont disponibles avec un délai de 48 à 96 heures et sont mises à jour quotidiennement en semaine.

Le chargement des statistiques peut être long.