Pressure transient analysis for multi-wing fractured vertical well in coalbed methane reservoir based on fractal geometry and fractional calculus

Abstract

In this paper, a semi-analytical mathematical model of pressure transient analysis (PTA) for multi-wing fractured vertical well (MWFV) in coalbed methane (CBM) reservoir is proposed, which considers the complexity of porous media by fractal geometry, the anomalous diffusion based on fractional calculus and the stress sensitivity represented by the exponential expression. Then through line source theory, dimensionless transformation, Pedrosa transformation, and other methods, the solution of the bottom hole pressure is obtained. At last, the PTA curve is presented by the Stehfest inversion method, and seven flowing regimes are identified. It can be observed that after introducing fractal geometry and fractional calculus, the PTA curve is quite different from the traditional curve in the slope of the derivative curve. The influences of related parameters are discussed, including mass fractal dimension, anomalous diffusion coefficient, number of hydraulic fractures, fracture angle, and stress sensitivity factor. Relevant results can provide better guidance to understand the CBM production performance in MWFV.

Keywords

Pressure transient analysis multi-wing fractured vertical well coalbed methane complexity of porous media anomalous diffusion stress sensitivity

Introduction

It was reported in BP ENERGY OUTLOOK in 2016 that to maintain the stability and increment of the global economic (Wei et al., 2021a), the energy demand will be in dramatic growth and fossil fuel will account for 80% of the whole energy supplement on a global scale (Wei et al., 2021b). For the aim to alleviate the burden of traditional petroleum energy supply, the unconventional natural gas resource is gradually getting people's attention (Nie et al., 2012). More importantly, the unconventional natural gas resource has an abundant distribution worldwide. For example, in China, a large number of unconventional natural gas resources such as tight sandstone gas (Sulige Gas Field) and coalbed methane (Taiyuan Formation) are developed in the Ordos Basin. Also, in the United States, numerous tight gas resources are stored in the western part of the Rocky Mountains, and in the San Juan Basin, there is also a large amount of coalbed methane (Miskimins and Jennifer, 2009). Coalbed methane (CBM), which is a member of the family of unconventional natural gas resources, is beginning to get much exploitation due to its high thermal efficiency and low pollution (Dong et al., 2017).

The permeability of CBM reservoirs is usually low and vertical fracturing technology is widely introduced to the exploitation of CBM to achieve high productivity (Zhao et al., 2016). Meanwhile, during the fracturing period, multiple fractures are formed (Figure 1) in the effect of stress instead of ideal single planner fracture which is known as multi-wing fractured vertical well (MWFV).

Figure 1.

Multiple fractures connected to the vertical wellbore (Chen et al., 2016).

Pressure transient analysis (PTA) is an efficient method for the monitoring of the performance of production well, and many researchers studied the PTA curve of MWFV. Choo and Wu (1987) researched MWFV in the gas reservoir. Through numerical method, they obtained the PTA curve and noticed that there is an obvious difference in the PTA curve between multi-wing fractures and single fracture when the dimensionless conductivity of fractures is small than 10. But the numerical method is complicated in representation and time-consuming in the calculation, a semi-analytical method was proposed later according to the research achievement of Cinco-Ley H (Cinco et al., 1978; Cinco and Meng, 1988). Craig and Blasingame (2006) constructed a vertical fractured well model with two fractures (one is secondary fracture) and proposed the solution of bottom hole pressure at the constant production rate. They analyzed the impacts of the length of two fractures, the conductivity of fractured, and the angle between two fractures. Luo and Tang (2015) studied the PTA curve of MWFV in detail based on the semi-analytical method. They found that there is a slight ‘hump’ in the pressure derivative curve which is called ‘interference regime’. Research on the PTA curve of MWFV in CBM reservoirs is also within the scope of researchers. Chen et al. (2016) established the corresponding model for MWFV in CBM reservoirs which took absorption, diffusion, and stress sensitivity into consideration. They also drew the corresponding PTA curve and analyzed the influence of various factors in detail. Zhang et al. (2018) proposed a model for MWFV in CBM reservoirs which coupled the stimulated region through the composite model. A transition regime in the PTA curve between the stimulated region and unstimulated region was put forward and the radius of the stimulated region was analyzed.

However, all research mentioned above is based on a uniform research background: the Euclidean geometric space, in which the reservoir is a homogeneous porous media and natural fractures (or cleats in CBM reservoirs) are distributed uniformly. While this is opposite to the fact that the distribution of cleats is extremely random in reservoirs (Wang et al., 2016). That is to say, the complexity of reservoirs is ignored in the above research. In order to modify this problem, the fractal theory is applied to revise. Based on fractal theory, Chang and Yortsos (1990) derived the formula of permeability and porosity (fractal permeability and porosity) which are the power-law function of reference distance. This method can better reflect the complexity of porous media and the reservoirs described is closer to reality and is widely used in much research. Through fractal permeability and porosity, Beier (1994) conducted research on the PTA curve of a vertical fractured well and found this method can match the field data better. Flamenco-Lopez and Camacho-Velazquez (2001) investigated the PTA curve of naturally fractured reservoirs which took fractal characteristics into consideration. Cossio et al. (2012) utilized the tri-linear model (Lee and Brockenbrough, 1986) to study the pressure performance of the fractal vertical fractured well. Also, Wang et al. (2015a, 2015b) introduced fractal theory into the tri-linear model to represent the stimulated region in multi-staged horizontal wells and analyzed the properties of the stimulated region on the PTA curve. Wang et al. (2019) coupled the fractal theory with the five-linear model (Stalgorova and Mattar, 2013) to study the imbibition effect on production performance of the multi-staged horizontal well.

Meanwhile, it is approved that due to the extremely random distribution of natural fractures, the characteristic of transport in porous media is transformed into anomalous diffusion (Metzler et al., 1994; Park et al., 2000). Aimed at anomalous diffusion, fractional calculus is proposed to describe this phenomenon, that is to say, modifying the darcy equation by fractional calculus. This method is widely used in the study of PTA considering anomalous diffusion (Raghavan, 2012; Raghavan and Chen, 2013; Razminia et al., 2015). In recent years, Ozcan et al. (2014) combined anomalous diffusion with the tri-linear model to obtain an analytical solution of a fractured horizontal well in tight formation and suggested that the new tri-linear anomalous-diffusion model is useful to make a transient analysis. Aimed at shale gas reservoirs, Ren and Guo (2015) studied the pressure response of multiple fractured horizontal wells considering anomalous diffusion, but the conductivity of hydraulic fractures was not accounted. Gu et al. (2019) derived a FFDM model (fractally fractional diffusion model) to study the PTA curve of multiple fractured horizontal wells (infinite conductivity) with the stimulated region. Jiang et al. (2021) established the PTA curve for multi-fractured horizontal wells in CBM reservoirs based on fractal and fractional calculus, but they also assumed that hydraulic fractures are infinite conductivity.

Through extensive literature research, lots of work had been done on the pressure transient of MWFV and fractal reservoir, but there is no study about the PTA research of MWFV that both consider the fractal geometry, anomalous diffusion, and stress sensitivity in CBM reservoir. In this paper, we derive a corresponding mathematical model based on fractal theory and fractional calculus and conduct some relevant analysis.

The following of this paper is arranged as:

The physical model of MWFV in CBM reservoirs is introduced and some assumptions are made to simplify this problem.

Based on the physical model, the mathematical model is established, and then the mathematical model is solved to obtain the pressure response of MWFV.

The verification of the results of this paper is conducted and a detailed analysis of the PTA curve is performed.

Physical model

The physical model for MWFV in coalbed methane is illuminated in Figure 2. The basic hypotheses are:

During the fracturing period, M hydraulic fractures centered on the wellbore are formed, which penetrate the reservoir and each of which has the same length (L_F), permeability (K_F), and with (W_F). The pressure drop along the fracture is considered. The gas firstly flows into the fractures, then flows into the wellbore through fractures.

The MWFV produces at the constant rate (q_sc), the whole reservoir shares the initial pressure (p_i).

The cleat system which is embedded in the Euclidean matrix system is represented as fractal geometry to consider the complexity of porous media, and the flow in the cleat system obeys anomalous diffusion. The gas flows from the matrix system to the cleat system by pseudo diffusion.

Ignoring the temperature change, capillary force, and gravity. Only a single gas phase is considered.

Figure 2.

The physical model for MWFV.

Establishment and solution of mathematical model

Model establishment

The whole system consists of two elements: reservoir and hydraulic fractures, therefore, the model establishment also includes two parts.

Model of reservoir

According to Appendix B, the flowing equation of the gas in the cleat system considering fractal geometry, anomalous diffusion, and stress sensitivity can be formulated as:

\begin{aligned} [\frac{\partial^{2} ψ_{c}}{\partial r^{2}} + \frac{d - θ - 1}{r} \frac{\partial ψ_{c}}{\partial r} + γ_{p} {(\frac{\partial ψ_{c}}{\partial r})}^{2}] \cdot {(\frac{r}{L_{ref}})}^{d - θ - 2} \\ = e^{γ_{p} (ψ_{i} - ψ_{c})} [\frac{c_{tc} ϕ_{cr}}{k_{cri}} {(\frac{r}{L_{ref}})}^{d - 2} \frac{\partial^{α} ψ_{c}}{\partial t^{α}} + \frac{2 p_{sc}}{T_{sc}} \frac{\partial^{α} C_{m}}{\partial t^{α}}] \end{aligned}

(1)

Based on the dimensionless variable listed in Appendix A, the dimensionless form of equation (1) is:

[\frac{1}{r_{D}^{θ}} \frac{\partial^{2} ψ_{cD}}{\partial r_{D}^{2}} + \frac{d - θ - 1}{r_{D}^{θ + 1}} \frac{\partial ψ_{cD}}{\partial r_{D}} - \frac{γ_{pD}}{r_{D}^{θ}} {(\frac{\partial ψ_{cD}}{\partial r_{D}})}^{2}] = e^{γ_{pD} ψ_{fD}} [ω \frac{\partial^{α} ψ_{cD}}{\partial t_{D}^{α}} + (1 - ω) \frac{\partial^{α} C_{mD}}{\partial t_{D}^{α}}]

(2)

The gas from the matrix system to the cleat system is described by Fick's first law (Xu et al., 2015):

(1 - ω) \frac{\partial^{α} C_{mD}}{\partial t_{D}^{α}} = - λ (C_{mD} - β ψ_{cD})

(3)

The corresponding initial condition, outer boundary condition, and inner boundary are:

ψ_{cD} |_{t = 0} = 0

(4)

ψ_{cD} |_{r_{D} \to \infty} = 0

(5)

r_{D}^{d - θ - 1} {\frac{\partial ψ_{cD}}{\partial r_{D}} |}_{r_{D} \to 0} = - e^{γ_{pD} ψ_{cD}} q_{D}

(6)

Model of hydraulic fractures

There are M hydraulic fractures in total and have the same properties, so it is reasonable to write the flowing equation for the i-th fracture only, and the other fractures’ flowing equations are the same. According to Appendices A and B, the dimensionless flowing equation can be described as:

\frac{\partial^{2} ψ_{FD i}}{\partial x_{D}^{2}} + \frac{2 e^{- γ_{pD} ψ_{cD}}}{C F D_{i}} r_{D}^{d - θ - 2} {\frac{\partial ψ_{cD}}{\partial x_{D}} |}_{y_{D} = W_{FD} / 2} = 0

(7)

The relevant boundary conditions are:

{\frac{\partial ψ_{FD i}}{\partial x_{D}} |}_{x_{D}^{^{'}} = 0} = - \frac{2 π}{C F D_{i}} q_{FD i}

(8)

ε_{FD i} = - \frac{e^{- γ_{pD} ψ_{cD}}}{π} r_{D}^{d - θ - 2} {\frac{\partial ψ_{fD}}{\partial y_{D}} |}_{y_{D} = W_{FD i} / 2}

(9)

Model solution

Solution for reservoir model

Equation (2) is strongly nonlinear due to the consideration of the stress sensitivity, therefore, Pedrosa transformation and Perturbation transformation are adopted to handle this problem.

Pedrosa transformation is shown as:

ψ_{cD} = - \frac{1}{γ_{pD}} \ln (1 - γ_{pD} σ_{cD})

(10)

Through equation (10), equations (2)–(6) is transformed into:

\frac{1}{r_{D}^{θ}} \frac{\partial^{2} σ_{cD}}{\partial r_{D}^{2}} + \frac{d - θ - 1}{r_{D}^{θ + 1}} \frac{\partial σ_{cD}}{\partial r_{D}} = \frac{ω}{{(1 - γ_{pD} σ_{cD})}^{α}} \frac{\partial^{α} ψ_{cD}}{\partial t_{D}^{α}} + (1 - ω) \frac{\partial^{α} C_{mD}}{\partial t_{D}^{α}}

(11)

(1 - ω) \frac{\partial^{α} C_{mD}}{\partial t_{D}^{α}} = - λ [C_{mD} + \frac{β}{γ_{pD}} \ln (1 - γ_{pD} σ_{cD})]

(12)

σ_{cD} |_{t_{D} = 0} = 0

(13)

σ_{cD} |_{r_{D} \to \infty} = 0

(14)

r_{D}^{d - θ - 1} {\frac{\partial σ_{cD}}{\partial r_{D}} |}_{r_{D} \to 0} = - q_{D}

(15)

Then, according to Perturbation transformation:

\begin{aligned} σ_{cD} = & σ_{cD 0} + γ_{pD} σ_{cD 1} + γ_{pD}^{2} σ_{cD 2} + \dots \\ - \frac{1}{γ_{pD}} \ln (1 - γ_{pD} σ_{cD}) = σ_{cD} + \frac{1}{2} γ_{pD} σ_{cD}^{2} + \dots \\ \frac{1}{1 - γ_{pD} σ_{cD}} = & 1 + γ_{pD} σ_{cD} + γ_{pD}^{2} σ_{cD} + \dots \end{aligned}

(16)

Due to the tiny value of γ_pD, equation (16) can be also formulated as:

\begin{aligned} σ_{cD} \approx & σ_{cD 0} \\ - \frac{1}{γ_{pD}} \ln (1 - γ_{pD} σ_{cD}) \approx & σ_{cD} \\ \frac{1}{1 - γ_{pD} σ_{cD}} \approx & 1 \end{aligned}

(17)

So, equations (11)–(15) can be furtherly written as:

\frac{1}{r_{D}^{θ}} \frac{\partial^{2} σ_{cD 0}}{\partial r_{D}^{2}} + \frac{d - θ - 1}{r_{D}^{θ + 1}} \frac{\partial σ_{cD 0}}{\partial r_{D}} = ω \frac{\partial^{α} ψ_{cD 0}}{\partial t_{D}^{α}} + (1 - ω) \frac{\partial^{α} C_{mD}}{\partial t_{D}^{α}}

(18)

(1 - ω) \frac{\partial^{α} C_{mD}}{\partial t_{D}^{α}} = - λ (C_{mD} - β σ_{cD 0})

(19)

σ_{cD 0} |_{t_{D} = 0} = 0

(20)

σ_{cD 0} |_{r_{D} \to \infty} = 0

(21)

r_{D}^{d - θ - 1} {\frac{\partial σ_{cD 0}}{\partial r_{D}} |}_{r_{D} \to 0} = - q_{D}

(22)

At last, taking Laplace transformation (in regard to t_D) to equations (18)–(22) and combing equations (18) and (19):

\frac{1}{r_{D}^{θ}} \frac{\partial^{2} {\bar{σ}}_{cD 0}}{\partial r_{D}^{2}} + \frac{d - θ - 1}{r_{D}^{θ + 1}} \frac{\partial {\bar{σ}}_{cD 0}}{\partial r_{D}} = s^{α} f (s) {\bar{σ}}_{cD 0}

(23)

{\bar{σ}}_{cD 0} |_{r_{D} \to \infty} = 0

(24)

r_{D}^{d - θ - 1} {\frac{\partial {\bar{σ}}_{cD 0}}{\partial r_{D}} |}_{r_{D} \to 0} = - {\bar{q}}_{D} (s)

(25)

Where,

f (s) = [ω + \frac{β (1 - ω)}{1 + λ s}] s^{α}

(26)

Making n = d − θ−1, m = d − 1, equation (23) is reshaped to:

\frac{\partial^{2} {\bar{σ}}_{cD 0}}{\partial r_{D}^{2}} + \frac{n}{r_{D}} \frac{\partial {\bar{σ}}_{cD 0}}{\partial r_{D}} - r_{D}^{θ} s^{α} f (s) {\bar{σ}}_{cD 0} = 0

(27)

The general solution of equation (27) is:

{\bar{σ}}_{cD 0} = r_{D}^{(1 - n) / 2} [a I_{v} (\frac{2 \sqrt{s^{α} f (s)}}{2 + m - n} r_{D}^{\frac{2 + m - n}{2}}) + b K_{v} (\frac{2 \sqrt{s^{α} f (s)}}{2 + m - n} r_{D}^{\frac{2 + m - n}{2}})]

(28)

Where,

v = \frac{1 - n}{2 + m - n}

(29)

Taking the outer boundary equation (24) and inner boundary equation (25) into equation (29), the unknown coefficient a and b are:

{\begin{matrix} a = 0 \\ b = \frac{2 {\bar{q}}_{D} (s)}{Γ (1 - v) {(2 + θ)}^{1 - v} {[\sqrt{s^{α} f (s)}]}^{v}} \end{matrix}

(30)

Where,

Γ (1 - v) = \int_{0}^{\infty} e^{- τ} τ^{- v} d τ

(31)

Finally,

{\bar{σ}}_{cD 0} = r_{D}^{(1 - n) / 2} B {\bar{q}}_{D} (s) [K_{v} (\frac{2 \sqrt{s^{α} f (s)}}{2 + m - n} r_{D}^{\frac{2 + m - n}{2}})], B = b / {\bar{q}}_{D} (s)

(32)

Equation (32) is the line-source solution for the reservoir model.

Solution for hydraulic fractures model

The method to solve hydraulic fractures model is analogous to solution for reservoir model, after Pedrosa transformation, Perturbation transformation, and Laplace transformation, Equations (7)–(9) are transformed into:

\frac{\partial^{2} {\bar{ψ}}_{FD i}}{\partial x_{D}^{2}} + \frac{2}{C F D_{i}} r_{D}^{d - θ - 2} {\frac{\partial {\bar{σ}}_{cD}}{\partial x_{D}} |}_{y_{D} = W_{FD} / 2} = 0

(33)

{\frac{\partial {\bar{ψ}}_{FD i}}{\partial x_{D}} |}_{x_{D} = 0} = - \frac{2 π}{C F D_{i}} {\bar{q}}_{FD i} (s)

(34)

{\bar{ε}}_{FD i} = - \frac{1}{π} r_{D}^{d - θ - 2} {\frac{\partial {\bar{σ}}_{cD 0}}{\partial y_{D}} |}_{y_{D} = W_{FD i} / 2}

(35)

Substituting equation (35) into equation (33) and taking double integral:

\int_{0}^{r_{D}} \int_{0}^{δ} \frac{\partial^{2} {\bar{ψ}}_{FD i}}{\partial x_{D}^{2}} d r_{D} d δ = \frac{2 π}{C F D_{i}} \int_{0}^{r_{D}} \int_{0}^{δ} {\bar{ε}}_{FD i} d r_{D} d δ

(36)

The double integral in the left side of equation (36) is:

\begin{aligned} \int_{0}^{r_{D}} \int_{0}^{δ} \frac{\partial^{2} {\bar{ψ}}_{FD i}}{\partial x_{D}^{2}} d r_{D} d δ = & {\bar{ψ}}_{FD i} (r_{D}) - {\bar{ψ}}_{FD i} (0) + \frac{2 π}{C F D_{i}} {\bar{q}}_{FD i} r_{D} \\ = & {\bar{ψ}}_{FD i} (r_{D}) - {\bar{ψ}}_{wD} + \frac{2 π}{C F D_{i}} {\bar{q}}_{FD i} r_{D} \end{aligned}

(37)

Then, equation (36) can be written as:

{\bar{ψ}}_{wD} - {\bar{ψ}}_{FD i} (r_{D}) = \frac{2 π}{C F D_{i}} ({\bar{q}}_{FD i} r_{D} - \int_{0}^{r_{D}} \int_{0}^{δ} {\bar{ε}}_{FD i} d r_{D} d δ)

(38)

Assuming that the pressure is equal at the interface between fractures and reservoir, equation (38) can also be expressed as:

{\bar{ψ}}_{wD} - {\bar{ψ}}_{cD} (r_{D}) \approx {\bar{ψ}}_{wD} - {\bar{σ}}_{cD 0} (r_{D}) = \frac{2 π}{C F D_{i}} ({\bar{q}}_{FD i} r_{D} - \int_{0}^{r_{D}} \int_{0}^{δ} {\bar{ε}}_{FD i} d r_{D} d δ)

(39)

Pressure response of MWFV

To obtain the pressure response of MWFV, discretization and superposition principles are utilized to cope with this problem. Each hydraulic fracture is split into N segments, as shown in Figure 3. Therefore, the length of each segment can be calculated by the length of the fracture, and the midpoint coordinate of the j-th segment in the i-th fracture can also be calculated. Applying the superposition principle, the total pressure response at any point r_D or (x_D, y_D) can be obtained by superimposing pressure responses of all segments at that point:

{\bar{σ}}_{cD 0 i, j} (r_{D}, u) = \sum_{j = 1}^{N} \sum_{i = 1}^{M} \int_{(j - 1) Δ L_{FD}}^{j Δ L_{FD}} {\bar{ε}}_{FD i, j} B r_{D i, j}^{(1 - n) / 2} K_{v} [\frac{2 \sqrt{s^{α} f (s)}}{2 + m - n} r_{D i, j}^{\frac{2 + m - n}{2}}] d r_{D}

(40)

Figure 3.

The discretization of i-th hydraulic fracture.

By the relation between radial coordinates and cartesian coordinates, equation (40) can be written as:

{\bar{σ}}_{cD 0 i, j} (x_{D}, y_{D}, u) = \sum_{j = 1}^{N} \sum_{i = 1}^{M} {\bar{ε}}_{FD i, j} E_{i, j} (x_{D}, y_{D})

(41)

Where,

E_{i, j} (x_{D}, y_{D}) = \int_{x_{D j}}^{x_{D j + 1}} {\begin{matrix} B {(\sqrt{\begin{matrix} {(x_{D} - l_{D} \cos θ_{i})}^{2} + \\ {(y_{D} - l_{D} \sin θ_{i})}^{2} \end{matrix}})}^{(1 - n) / 2} \cdot \\ K_{v} [\frac{2 \sqrt{s^{α} f (s)}}{2 + m - n} {(\sqrt{\begin{matrix} {(x_{D} - l_{D} \cos θ_{i})}^{2} + \\ {(y_{D} - l_{D} \sin θ_{i})}^{2} \end{matrix}})}^{\frac{2 + m - n}{2}}] \end{matrix}} d l_{\begin{matrix} D \end{matrix}}

(42)

Point (x_D, y_D) can be placed at the midpoint of the segment, and equation (40) is further written as:

{\bar{ψ}}_{wD} - \sum_{j = 1}^{N} \sum_{i = 1}^{M} {\bar{ε}}_{FD i, j} E_{i, j} (X_{D a, b}, Y_{D a, b}) = \frac{2 π}{C F D_{i}} ({\bar{q}}_{FD i} r_{D} - \int_{0}^{r_{D}} \int_{0}^{δ} {\bar{ε}}_{FD i} d r_{D} d δ)

(43)

Combining partial integration and numerical integration, equation (43) is finally transformed into:

\begin{aligned} {\bar{ψ}}_{wD} - \sum_{j = 1}^{N} \sum_{i = 1}^{M} {\bar{ε}}_{FD i, j} E_{i, j} (X_{D a, b}, Y_{D a, b}) \\ = & \frac{2 π}{C F D} ({\bar{q}}_{FDa} r_{D a, b} - \sum_{v = 1}^{b - 1} {\bar{ε}}_{FD a, v} (i - v) Δ L_{FD a, b}^{2} + {\bar{ε}}_{FD a, b} \frac{Δ L_{FD a, b}^{2}}{8}) \end{aligned}

(44)

By taking all midpoints into equation (44), an M × N order matrix formula is obtained, but there are M × N + M + 1 unknowns. Two other restrictions are listed as follows.

The sum of the output of all segments in each fracture is equal to the whole output of the fracture:

\sum_{b = 1}^{N} {\bar{ε}}_{FD a, b} Δ L_{FD a, b} = {\bar{q}}_{FDa}

(45)

The sum of the output of all fractures is equal to the whole output of the well:

\sum_{a = 1}^{M} {\bar{q}}_{FDa} = \frac{1}{s}

(46)

Combing equations (44)–(46), the bottomhole pressure in the Laplace domain can be calculated. And through Duhamel's principle (Van Everdingen and Hurst, 1949), the wellbore storage and skin effect are taken into consideration:

{\bar{ψ}}_{wD} (s, S, C_{D}) = \frac{S + u {\bar{ψ}}_{wD} (s)}{u + u^{2} C_{D} (S + u {\bar{ψ}}_{wD} (s))}

(47)

Then by Stehfest numerical inversion (Stehfest, 1970), we can get the solution in the real-time domain.

Simulation results

Verification

In order to verify the accuracy of the fractal solution of this paper, a comparison between the fractal solution (solved by the semi-analytical method in this paper) and the numerical solution (solved by COMSOL). As shown in Figures 4 and 5, the red solid line has good consistent with the green hollow dotted line and the relative error is acceptable, which approve that the fractal solution is reliable.

Figure 4.

Comparison with numerical solution.

Figure 5.

Relative error between semi-analytical solution and numerical solution.

Also, making d = 2, θ = 0, the model in this paper can degenerate to the model analogous to Ren's model. As shown in Figure 6, there is a good agreement between these two results, which is explaining the correctness of the model in this paper.

Figure 6.

Comparison with Ren's model.

Typical PTA curve

According to different features of the PTA curve, seven flowing regimes are split:

I: Wellbore storage and skin effect regime. II: Bi-linear flowing regime. The slope of the derivative curve is bigger than 1/4 compared to the traditional model. III: Interference regime. A slight hump in the derivative curve is a feature in this regime. IV: Linear flowing regime. The derivative curve tends to have a slope greater than 1/2, which is also different from the traditional model. V: Transition regime. VI: Diffusion regime. For pseudo steady diffusion, there is an obvious dip in the derivative curve, which illustrates that CBM is flowing from matrix to cleat under the action of concentration difference. But for transient diffusion, there is also a dip which is less obvious than that if the diffusion is pseudo steady. VII: Later radial flowing regime. After considering fractal geometry and anomalous diffusion, the derivative curve is no more a horizontal line. (Figure 7)

Figure 7.

Typical PTA curve of the pressure response curve of MWFV.

Discussions

Mass fractal dimension and anomalous diffusion index

As shown in Figures 8 and 9. Mass fractal dimension and anomalous diffusion index have an opposite impact on the PTA curve. When the mass fractal dimension becomes smaller, the position of curves in the Bi-linear flowing regime, Interference regime, and Linear flowing regime become up, while the position of curves in the Diffusion regime and Later radial flowing regime becomes down. Instead, the PTA curve shows an inverse phenomenon as the anomalous diffusion index increase. In addition, the increment of anomalous diffusion index also makes flowing regimes after flowing regime I delay.

Figure 8.

The impact of mass fractal dimension.

Figure 9.

The impact of the anomalous diffusion index.

Number of hydraulic fractures

In Figure 10, the impact of the number of hydraulic fractures is mainly in Bi-linear flowing, Interference, and Linear flowing regime. First, the position of the PTA curve in these regimes goes down as the number of fractures increases. Additionally, the hump of the derivative curve in the Interference regime becomes more obvious, which means more fractures represent worse interference.

Figure 10.

The impact of the number of hydraulic fractures.

Hydraulic fracture length

With the increment of hydraulic fracture length, the stimulated region becomes larger, which means the energy wastage required for gas flowing becomes smaller. So, it can be observed from Figure 11 that a bigger fracture length makes the position of the PTA curve go down in Bi-linear flowing, Interference, and Linear flowing regimes. Meanwhile, the duration of Bi-linear flowing is longer, and the duration of the Linear flowing regime is the opposite.

Figure 11.

The impact of fracture length.

Angle between adjacent hydraulic fracture

When the angle between adjacent hydraulic fractures decreases, the Interference regime appears early and the hump in the derivative curve is more obvious (Figure 12). A smaller angle represents that the adjacent fractures are very close to each other, which means the interference between adjacent fractures is easier to appear and the greater energy loss in the whole production. Therefore, the relevant phenomenon is formed.

Figure 12.

The impact of adjacent hydraulic fracture angle.

Hydraulic fracture conductivity

In Figure 13, it can be seen that with the increment of fracture conductivity, the duration of the Bi-linear flowing regime becomes shorter, but the duration of the Linear flowing regime becomes longer. If the value of fracture conductivity is large enough, which means the fracture has infinite conductivity, the Bi-linear flowing regime will be covered up, leaving the Linear flowing regime only.

Figure 13.

The impact of hydraulic fracture conductivity.

Stress sensitivity factor

With the increment of stress sensitivity factor, the PTA curve in Later radial flowing regime upwarps (Figure 14). This is because the stress sensitivity effect will be reflected when the pressure difference increases to a certain extent, which leads to the increment of flowing resistance and energy wastage and is not conducive to production. Therefore, it shows the above phenomenon.

Figure 14.

The impact of stress sensitivity factor.

Field application

In this part, one field application is conducted based on the model above. Well X3 is an MWFV of CBM reservoir in China. Through the micro-seismogram result, there are five fracture wings formed during the fracturing period. The matching of the PTA curve is shown in Figure 15 and relevant interpretation results are listed in Table 1. It should be noted that due to the short time of shut-in, the whole flowing regimes are not complete. Also, the interpreted adsorption coefficient, transfer coefficient, and storage coefficient are not necessarily accurate because the Diffusion regime is not identified.

Figure 15.

The matching result of X3 well.

Table 1.

The interpretation results of X3 well.

Parameters	Value
Mass Fractal Dimension	1.92
Anomalous Diffusion index	0.05
Hydraulic fracture length	51.3 m
Fracture conductivity	44.5 mD·m
Wellbore storage coefficient	2.14 × 10⁻⁷ m³/Pa
Skin factor	0.01
Adsorption coefficient	0.78
Transfer coefficient	0.02
Storage coefficient	0.21

Conclusions

A new mathematical model for multi-wing fractured vertical well in coalbed methane reservoir is established which takes the complexity of porous media, anomalous diffusion, and stress sensitivity into consideration. Then the solution is obtained by a semi-analytical method.

The corresponding PTA curve is drawn and divided into seven different flowing regimes. The main difference compared to the traditional one is reflected in the slope of the derivative curve. The slope is bigger than 1/4 in the Bi-linear flowing regime or 1/2 in the Linear flowing regime. And the derivative is no more a horizontal line in the Later radial flowing regime.

Sensitivity analysis reveals that mass fractal dimension and anomalous diffusion index have a great influence on the whole PTA curve, also the increasing anomalous diffusion index leads to the delay of flowing regimes after regime I. Parameters related to the geometry of hydraulic fractures mainly affect the middle flowing regime. Stress sensitivity factor cause the upward of later flowing regimes.

Field application shows that based on the proposed model, the PTA curve has a good fitting result which represents the practicality of this paper's model.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Key S&T Special Projects, The National Natural Science Foundation of China, (grant number 2017ZX05008-004-004-001, 41972129).

ORCID iD

Haoyuan Li

Appendix

References

Beier

(1994) Pressure transient model for a vertically fractured well in a fractal reservoir. SPE Formation Evaluation 9(02): 122–128.

Chang

Yortsos

(1990) Pressure transient analysis of fractal reservoirs. SPE Formation Evaluation 5(01): 31–38.

Chen

Liao

Zhao

, et al. (2016) A semi-analytical mathematical model for transient pressure behavior of multiple fractured vertical well in coal reservoirs incorporating with diffusion, adsorption, and stress-sensitivity. Journal of Natural Gas and Science Engineering 2016: 570–582.

Choo

(1987) Transient pressure behavior of multiple-fractured gas wells. SPE/DOE Joint Symposium on Low Permeability Reservoirs, Denver, Colorado, USA, 18–19 May, Society of Petroleum Engineers.

Cinco

Meng

(1988) Pressure transient analysis of wells with finite conductivity vertical fractures in double porosity reservoirs. SPE Annual Technical Conference and Exhibition, Houston, Texas, USA, 2–5 October, Society of Petroleum Engineers.

Cinco

Samaniego

Dominguez

(1978) Transient pressure behavior for a well with a finite-conductivity vertical fracture. SPE Journal 18(04): 253–264.

Cossio

(2012) A Semi-analytic solution for flow in finite conductivity vertical fractures using fractal theory. SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA, 8–10 October, Society of Petroleum Engineers.

Craig

Blasingame

(2006) Constant-rate drawdown solutions derived for multiple arbitrarily oriented uniform-flux, infinite-conductivity, or finite-conductivity fractures in an infinite-slab reservoir. SPE Gas Technology Symposium, Calgary, Alberta, Canada, 15–17 May, Society of Petroleum Engineers.

Dong

Cheng

Jin

, et al. (2017) Effects of diffusion and suction negative pressure on coalbed methane extraction and a new measure to increase the methane utilization rate. Fuel 197: 70–81.

10.

Flamenco-López

Camacho-Velázquez

(2001) Fractal transient pressure behavior of naturally fractured reservoirs. SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, USA, 30 September–31 October, Society of Petroleum Engineers.

11.

Ding

Gao

, et al. (2019) A fractally fractional diffusion model of composite dual-porosity for multiple fractured horizontal wells with stimulated reservoir volume in tight gas reservoirs. Journal of Petroleum Science and Engineering 173: 53–68.

12.

Jiang

Liu

Wang

, et al. (2021) A semi-analytical fractal-fractional mathematical model for multi-fractured horizontal wells in coalbed methane reservoirs. Journal of Energy Resource and Technology 143(1): 1–18.

13.

Lee

Brockenbrough

(1986) A new approximate analytic solution for finite-conductivity vertical fractures. SPE Formation Evaluation 1(01): 75–78.

14.

Luo

Tang

(2015) Pressure-transient analysis of multiwing fractures connected to a vertical wellbore. SPE Journal 20(02): 360–367.

15.

Metzler

Glockle

Nonnenmacher

(1994) Fractional model equation for anomalous diffusion. Physica A Statistical Mechanics & Its Applications 211(1): 13–24.

16.

Miskimins

Jennifer (2009) Design and life-cycle considerations for unconventional reservoir wells. SPE Production & Operations 24(02): 353–359.

17.

Nie

Meng

Guo

, et al. (2012) Modeling transient flow behavior of a horizontal well in a coal seam. International Journal of Coal Geology 92(2): 54–68.

18.

Ozcan

Sarak

Ozkan

, et al. (2014) A trilinear flow model for a fractured horizontal well in a fractal unconventional reservoir. SPE Annual Technical Conference and Exhibition, Amsterdam, The Netherlands, 27–29 October, Society of Petroleum Engineers.

19.

Park

Choe

Kang

(2000) Pressure behavior of transport in fractal porous media using a fractional calculus approach. Energy Sources 22(10): 881–890.

20.

Raghavan

(2012) Fractional diffusion: performance of fractured wells. Journal of Petroleum Science and Engineering 92–93: 167–173.

21.

Raghavan

Chen

(2013) Fractional diffusion in rocks produced by horizontal wells with multiple, transverse hydraulic fractures of finite conductivity. Journal of Petroleum Science and Engineering 109: 167–173.

22.

Razminia

Baleanu

(2015) Investigation of the fractional diffusion equation based on generalized integral quadrature technique. Applied Mathematical Modelling 39(1): 86–98.

23.

Ren

Guo

(2015) Anomalous diffusion performance of multiple fractured horizontal wells in shale gas reservoirs. Journal of Natural Gas and Science Engineering 26: 642–651.

24.

Stalgorova

Mattar

(2013) Analytical model for unconventional multifractured composite systems. SPE Reservoir Evaluation & Engineering 16(03): 246–256.

25.

Stehfest

(1970) Algorithm 368: numerical inversion of Laplace transforms (D5). Communications of the ACM 13(1): 47–49.

26.

Van Everdingen

Hurst

(1949) The application of the Laplace transformation to flow problems in reservoirs. Transactions of AIME 186(305): 97–104.

27.

Wang

Shahvali

(2015a) A semi-analytical fractal model for production from tight oil reservoirs with hydraulically fractured horizontal wells. Fuel 158: 612–618.

28.

Wang

Sheng

, et al. (2015b) A mathematical model considering complex fractures and fractal flow for pressure transient analysis of fractured horizontal wells in unconventional reservoirs. Journal of Natural Gas and Science Engineering 23: 139–147.

29.

Wang

Wei

(2016) Semi-analytical modeling of flow behavior in fractured media with fractal geometry. Transport in Porous Media 112: 707–736.

30.

Wang

Yang

Ding

, et al. (2019) A multi-linear fractal model for pressure transient analysis of multiple fractured horizontal wells in tight oil reservoirs including imbibition. Fractals 27: 1940004.

31.

Wei

Cheng

She

, et al. (2021a) Numerical study on the heat transfer behavior in naturally fractured reservoirs and applications for reservoir characterization and geothermal energy development. Journal of Petroleum Science and Engineering 202: 108560.

32.

Wei

Cheng

Wang

, et al. (2021b) Numerical study on temperature behavior in naturally fractured geothermal reservoirs and analysis methodology for geothermal reservoir characterization and development. SPE Annual Technical Conference and Exhibition, Dubai, UAE, 21–23 September, Society of Petroleum Engineers.

33.

Guo

Wei

, et al. (2015) Production performance analysis for composite shale gas reservoir considering multiple transport mechanisms. Journal of Natural Gas Science and Engineering 26: 382–395.

34.

Zhao

Zhang

Feng

, et al. (2016) Performance analysis of fractured wells with stimulated reservoir volume in coal seam reservoirs. Oil & Gas Science Technology 71: 1.

35.

Zhang

Kou

Wang

, et al. (2018) Performance analysis for a model of a multi-wing hydraulically fractured vertical well in a coalbed methane gas reservoir. Journal of Petroleum Science and Engineering 166: 104–120.