Regular Article
Characteristics of transient pressure performance of horizontal wells in fracturedvuggy tight fractal reservoirs considering nonlinear seepage
^{1}
School of Petroleum Engineering, China University of Petroleum (East China), Qingdao 266580, PR China
^{2}
The Pennsylvania State University, University Park, PA 16802, USA
^{*} Corresponding authors: 1215797855@qq.com; 986012825@qq.com
Received:
22
January
2019
Accepted:
10
April
2019
Since the classical seepage theory has limitations in characterizing the heterogeneity of fracturedvuggy tight reservoirs, well test interpretation results are not consistent with actual production by far. Based on the nonlinear percolation theory, a new nonlinear seepage equation considering the boundary layer and yield stress was derived to describe the seepage characteristics of dense matrix blocks and the stress sensitivity and fractal features of fracture systems were characterized by applying the fractal theory. Thus, the nonlinear model of a horizontal well in a fracturedvuggy tight fractal reservoir was established naturally. Then the finite element method was applied to solve the bottom hole pressure based on the processing of internal boundary conditions. After solving the model, the seepage characteristics of different models were summarized by analyzing the bottom hole pressure dynamic curves and the sensitivity analysis of multiple parameters such the nonlinear parameter and fractal index were conducted. Finally, the practicality of the model was proved through a field application. The results show that the pressure dynamic curves can be divided into nine flow stages and the increase of the nonlinear parameter will cause the intensity of the cross flow from matrix blocks to the fracture system to decrease. The fractal index is irrelevant to the intensity of the cross flow while it decides the upwarping degree of the curve at the middle and late flow stages. On the basis of the results of the field application, it can be concluded that the model fits well with actual production and the application of this model can improve the accuracy of well test interpretation.
© R. Jiang et al., published by IFP Energies nouvelles, 2019
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
D_{f}, θ: Fractal dimension, anomalous diffusion coefficient
C_{mt}, C_{ft}, C_{vt} : Comprehensive compression coefficient of matrix, fracture and cave system, MPa^{−1}
r_{w}, r_{e} : Radius of the wellbore and reservoir, m
K_{fh}, K_{fv} : Horizontal and vertical permeability of fracture system, mD
K_{m}, K_{v} : Permeability of matrix and cave system, mD
p_{m}, p_{f}, p_{v} : Pressure of matrix, fracture and cave system, MPa^{−1}
p_{D} : Dimensionless bottom hole pressure without considering wellbore reservoir coefficient and skin factor
${\stackrel{\u0305}{p}}_{\mathrm{wD}}$, p_{wD} : Dimensionless bottom hole pressure considering wellbore storage coefficient and skin factor in Laplace space and real space
ϕ_{fw}, ϕ_{w} : Porosity of fracture system and entire reservoir at the point of the production well
K_{fw}, K_{w} : Permeability of fracture system and entire reservoir at the point of production well, mD
ϕ_{m}, ϕ_{f}, ϕ_{v} : Porosity of matrix, fracture and cave system
C_{f}, α_{K} : Porosity compressibility coefficient and permeability modulus of fracture system, MPa^{−1}
r_{0}, L_{0} : Capillary radius and length, m
n_{0} : Number of capillaries per unit crosssectional area
δ: Boundary layer thickness, μm
τ_{0} : Yield stress value, MPa
h, L: Reservoir thickness and half length of the horizontal well, m
z_{w} : Horizontal well depth, m
c_{1} : Nonlinear parameter under the influence of the fluid yield stress and boundary layer
c_{2} : Nonlinear parameter under the influence of boundary layer
σ_{m}, σ_{v} : Form factor of matrix and cave system
λ_{vf}, λ_{mv}, λ_{mf} : Cross flow factor of cave system to fracture system, matrix system to cave system and matrix blocks to fracture system
ω_{v}, ω_{f} : Elastic storativity ratio of cave system and fracture system
δ_{LVDf}, δ_{LVDv} : Dimensionless intermediate variables of fracture system and cave system
c_{D} : Dimensionless nonlinear parameter
K_{e} : Nonlinear coefficient matrix
1 Introduction
The fracturedvuggy carbonate reservoir has a large number of caves that occur in the interior of the bedrock and its fracture distribution is quite complex. Therefore, the concept of triplemedium is applied to characterize the complex pores of such reservoirs (Ge and Wu, 1982; Wang Y. et al., 2017). Studying the seepage model of fracturedvuggy tight carbonate reservoirs requires full consideration of the nonlinear flow of matrix blocks and the stress sensitivity and fractal characteristics of a fracture system (Guo et al., 2015; Li et al., 2017; Lin and Myers, 2018; Wang W.D. et al., 2017; Zeng et al., 2011). The classical seepage theory has many limitations in characterizing the seepage mode of such reservoirs.
A group of international and domestic academics have done lots of research on the seepage model of fracturedvuggy reservoirs. In order to use simple form and definite physical meaning of mathematical language to solve practical production problems, Yao et al. (2004) conducted the well test interpretation of triplemedium reservoirs considering the variable well storage, innovating the well test theory of triplemedium reservoirs. CamachoVelazquez et al. (2002) proposed a triplemedium model characterizing a fracturevuggy carbonate reservoir and analyzed the model curves in detail by considering the solutions of the situation that the cave system connects with the wellbore. Popov et al. (2009) proposed a unified approach to the problem, which is relevant in flow simulations of fracturedvuggy reservoirs, where caves are embedded in a porous rock and are connected via fracture networks on multiple scales by using the StokesBrinkman equations on the fine scale. Nie et al. (2011) established a new flow model in a triple media carbonate reservoir which considers the unsteady interporosity flow based on the hypothesis that the shape of the cave is spherical. Gomez et al. (2006) applied the global optimization solution to research the well test method of fracturevuggy reservoirs. Wang and Yi (2018) presented a coupling mathematical model for the fractured horizontal well in a triple media carbonate reservoir by conceptualizing caves as spherical shapes in which the infinite conductivity of the acid fractures was taken into account.
The classical seepage theory is based on Euclidean space and has a certain characteristic scale, so Darcy’s law cannot describe the complexity and heterogeneity of the fracturevuggy carbonate reservoir well. It is worth mentioning that applying the fractal theory to petroleum industry is a great innovation, which can effectively characterize the permeability and porosity of heterogeneous reservoirs. Chang and Yortsos (1990) presented a formulation for a fractal fracture network embedded into Euclidean matrix blocks. Singlephase flow in the fractal object was described by an appropriate modification of the diffusivity equation. In addition, the viability of the model in explaining parameters such the transmissibility, storativity and fractal dimension of a complex naturally fractured geothermal reservoir was analyzed through the pressure transient test data (Aprilian et al., 1993). What is more, the fractal theory was applied to real well tests in various fractured reservoirs and the physical meaning of the fractal parameters was presented in the context of well testing in 1995 (Acuna et al., 1995). From then on, seepage models of multiplemedium reservoirs developed rapidly since the application of the fractal theory (Kong et al., 2008; Lu et al., 2017; Razminia et al., 2015; Wang et al., 2015; Xu et al., 2008). In addition to that, Xu et al. (2017) developed a fractal dualporosity model to study the singlephase fluid flow through fractured porous media, which can capture the statistical characteristics of fractures and shed light on the transport mechanism of fractured porous media compared with empirical formulas for effective permeability.
Except applying the fractal theory to deal with the fracture system, studying the seepage model of fracturedvuggy tight fractal reservoirs should also take the nonlinear flow of matrix blocks into consideration (Chen and Yao, 2017; Xu et al., 2013, 2019). In recent years, models of lowspeed and nonDarcy seepage flow developed rapidly (Ezulike and Dehghanpour, 2013; Vijayalakshmi et al., 2009; Wu et al., 2019). When it comes to the mathematical descriptions of a lowspeed nonlinear seepage phenomenon, the commonly adopted models are mainly divided into three categories: quasi startup pressure gradient model, section model and continuous model. Quasi startup pressure gradient model ignores the bending section of the seepage curve, so the model narrows the scope of fluid flow in the lowpermeability reservoirs. Section model does not reflect the nonlinear section of the actual seepage condition, and it is difficult to judge the critical point of linear and nonlinear seepage. Continuous models are all based on the seepage experiment and use different mathematical functions to fit the experimental data, but their physical meanings are not rich enough.
Although the predecessors have carried out a lot of research and summarized the seepage characteristics of fracturedvuggy tight fractal reservoirs, the established seepage models ignore the nonlinear and fractal characteristics, which are too simplified and cannot reflect the heterogeneity characteristics of such reservoirs. This paper derives a new nonlinear seepage equation considering the boundary layer and yield stress to describe the seepage characteristics of dense matrix blocks and considers fractal characteristics of the fracture system, so that a seepage model of the horizontal well in a triplemedium reservoir is established. In the conceptual mathematical model, fractured vuggy rock is considered as a triplecontinuum medium, consisting of fractures, rock matrix, and vugs. There is a strong demand for processing the inner boundary in order to obtain the bottom hole pressure solution based on the finite element principle. After comparing the seepage laws of different models and conducting sensitivity analysis of multiple parameters, the model is applied to conduct a well test interpretation of an actual reservoir.
2 Seepage mechanism of fracturedvuggy tight fractal reservoirs
2.1 Nonlinear seepage mechanism of a matrix system
2.1.1 Microscopic mechanism of the nonlinear seepage
The research shows that the main reasons of nonlinear seepage flow in low permeability reservoirs are described as below (Jiang et al., 2012; Xu et al., 2019):

Low permeability reservoirs have low porosity and permeability, small pores and complex microscopic pore structure. These characteristics cause the fluid to have a microscale flow effect on the lowpermeability reservoir.

The effect of the fluid yield stress and the adsorption of the boundary layer cannot be ignored when the size of pore throat is tiny, belonging to the micron level. The thickness of boundary layer becomes thinner with the displacement pressure rising, the percentage of flowable fluid increasing, as well as the viscosity of the fluid increasing, so the fluid in the pore needs to overcome greater fluid yield stress.

Any kind of the fluid has a certain yield stress, and only when the displacement pressure can overcome the yield stress will the fluid begin to flow. In lowpermeability reservoirs, small pore throats result in a greater yield stress, and nonNewtonian enhancements make the starting pressure gradient more significant.
Yield stress, adsorption of the boundary layer and microscale effects are not alone but mutually affected and coupled. The presence of microscale effects makes the influence of surface force is much larger than the volume force and needs to be regarded as a dominant force. The larger the surface force is, the stronger the reservoir fluid is adsorbed by the boundary constraints. While the presence of the adsorbed boundary layer makes the seepage channel smaller and exacerbates the microscale flow effect, resulting in nonlinear seepage strengthening (Cao et al., 2016; Cossio et al., 2013).
This paper combines the microseepage mechanism of lowpermeability reservoirs, the boundary layer flow and capillary seepage model to derive a new lowspeed nonlinear seepage model. It is worth noting that the values of the boundary layer and the starting pressure gradient in lowpermeability reservoirs are based on experimental data. The new model of this paper can deal with this problem flexibly, which the paper will introduce in detail later.
2.1.2 Derivation of a nonlinear seepage new model
The radius of pore throat in dense matrix blocks is tiny, which can be measured by micrometer, so the effect of the fluid yield stress and the adsorption of the boundary layer cannot be ignored (Wu et al., 2017). Considering the adsorbed boundary layer and the fluid yield stress, the modified form of the HagenPoiseuille law can be expressed as:
$$q=\frac{\pi ({r}_{0}\delta {)}^{4}}{8\mu}\left(\nabla p\frac{8{\tau}_{0}}{3\left({r}_{0}\delta \right)}\right),$$(1)
where δ is the thickness of the boundary layer, reflecting the influence of the boundary layer on seepage flow. τ_{0} is the fluid yield stress, reflecting the influence of fluid nonNewtonian characteristics on seepage flow.
The reservoir core is equivalent to a set of parallel capillary tubes with radius r_{0} embedded in solids. If there are n_{0} such capillaries per unit of crosssectional area, the specific discharge through the porous media block is:
$$v=\frac{{K}_{0}}{\mu}{\left(1\frac{\delta}{{r}_{0}}\right)}^{4}\left(1\frac{8{\tau}_{0}}{3{r}_{0}\left(1\frac{\delta}{{r}_{0}}\right)\nabla p}\right)\nabla p,$$(2)
In equation (2): ${\varphi}_{0}={n}_{0}\mathrm{A\pi}{{r}_{0}}^{2}\frac{{L}_{0}}{A{L}_{0}}={n}_{0}\pi {{r}_{0}}^{2}$; $K={n}_{0}\frac{\pi {{r}_{0}}^{4}}{8}=\frac{{\varphi}_{0}{{r}_{0}}^{2}}{8}$.
In the same capillary tube, the larger the pressure gradient is, the thinner the boundary layer is. Therefore, it is assumed that δ/r_{0} = a_{1}/∇p. If the yield stress value of the same fluid remains unchanged, then 8τ_{0}/3r_{0} can be regarded as a constant value, namely 8τ_{0}/3r_{0} = a_{2} (the values of a_{1} and a_{2} are obtained by experimental fitting) (Celata et al., 2006; Xu and Yue, 2007). Equation (2) can be converted into:
$$v=\frac{{K}_{0}}{\mu}\left[1\frac{4{a}_{1}+{a}_{2}}{\nabla p}+\frac{6{{a}_{1}}^{2}+3{a}_{1}{a}_{2}}{\nabla p(\nabla p{a}_{1})}\left(\frac{4{{a}_{1}}^{3}}{(\nabla p{)}^{3}}+\frac{6{{a}_{1}}^{2}{a}_{2}+6{{a}_{1}}^{3}}{\nabla {p}^{2}(\nabla p{\mathrm{a}}_{1})}\right)+\left(\frac{{{a}_{1}}^{4}}{{\left(\nabla p\right)}^{4}}+\frac{4{{a}_{1}}^{3}{a}_{2}}{\nabla {p}^{3}(\nabla p{a}_{1})}\right)\frac{{{a}_{1}}^{4}{a}_{2}}{\nabla {p}^{4}(\nabla p{a}_{1})}\right]\nabla p$$(3)
Since $\frac{{a}_{1}}{\nabla p}=\frac{\delta}{{r}_{0}}<1$, $\frac{{a}_{2}}{\nabla p{c}_{1}}=\frac{8{\tau}_{0}}{3{r}_{0}(1\frac{\delta}{{r}_{0}})\bullet \nabla p}<1$, the higher order term can be ignored and the formula can be simplified as:
$$v=\frac{{K}_{0}}{\mu}\left[1\frac{4{a}_{1}+{a}_{2}}{\nabla p}\frac{6{{a}_{1}}^{2}+3{a}_{1}{a}_{2}}{\nabla p\left(\nabla p{a}_{1}\right)}\right]\nabla p.$$(4)
c_{1}, c_{2}, and c_{3} are applied to denote the different polynomials of a_{1}, a_{2} ${c}_{1}=4{a}_{1}+{a}_{2}$, c_{2} = a_{1}, ${c}_{3}=6{{a}_{1}}^{2}+3{a}_{1}{a}_{2}$). In other words, the value of c_{1}, c_{2}, and c_{3} can only be obtained by experimental fitting. In order to simplify the formula, c_{3} can be represented by c_{1} and c_{2} through exploring the relationship between c_{1}, c_{2}, and c_{3}. Therefore, the new model can be simplified as:
$$\nu =\frac{K}{\mu}\left(1\frac{{c}_{1}}{\nabla p}\frac{{c}_{3}}{\nabla p(\nabla p{c}_{2})}\right)\nabla p=\frac{K}{\mu}\left(\nabla p{c}_{1}\frac{{c}_{3}}{\nabla p{c}_{2}}\right).$$(5)
Substituting the condition “∇p = 0，v = 0” into equation (5), then c_{3} can be represented by c_{1} and c_{2} (c_{3} = c_{1} c_{2}). Therefore, the new model can be simplified as:
$$\nu =\frac{K}{\mu}\left(1\frac{{c}_{1}}{\nabla p{c}_{2}}\right)\nabla p$$(6)
In equation (5), the parameter c_{1} reflects the yield stress of the fluid and the influence of the boundary layer while parameter c_{2} mainly reflects the influence of boundary layer on the seepage flow. It means c_{1}/(∇p – c_{2}) is equivalent to the starting pressure gradient term. When c_{1} = 0, the above model is simplified as the Darcy model. But if c_{2} = 0, the above model is simplified as the quasistarting pressure gradient model. When v ≥ 0, then ∇p ≥ c_{1} + c_{2}, so the real minimum starting pressure gradient is ∇p_{min} = c_{1} + c_{2}.
When c_{1} + c_{2} = 0, the ∇p_{min} = 0, if c_{1} + c_{2} > 0, then ∇p_{min} > 0. When there is a true starting pressure gradient in the seepage, the seepage curve does not pass through the origin and the model is expressed as:
$$\nu =\frac{{K}_{0}}{\mu}\left(\nabla p\frac{{c}_{1}}{1{c}_{2}\nabla p}\right).$$(7)
Under this circumstance, the seepage curve does not pass through the origin point. As the pressure gradient increases, the curve tends to be straight because when the pressure gradient increases to a certain extent, all pore throats participate in the flow. The thickness of the boundary layer no longer changes, and the quasilinear flow occurs.
2.2 Fractal characteristics of the fracture system
Fractal theory thinks that if fractures have fractal characteristics in fracturedvuggy reservoirs, the permeability or porosity is not only a function of pressure, but also relates to the fractal dimension and anomalous diffusion coefficient. D_{f} denotes fractal dimension which is a symbol of the complexity of a fractal system and it can reflect the geometric characteristics of the reservoirs. θ is anomalous diffusion coefficient which depicts the connectivity and geometric characteristics of a fractal fracture network. It is reasonable to characterize such reservoirs by applying the fractal theory, which can accurately reflect the heterogeneity of the fracture morphology so that the fractal structure can be more in line with the actual reservoir.
The porosity of the reservoir is exponentially related to the effective overburden pressure, so the expression of the porosity with stress sensitivity can be obtained:
$${\varphi}_{\mathrm{f}}\left(r\right)={\varphi}_{\mathrm{fw}}\left(\frac{r}{{r}_{\mathrm{w}}}\right){r}^{{D}_{\mathrm{f}}d}\mathrm{exp}\left[{C}_{\mathrm{f}}\left({p}_{\mathrm{i}}{p}_{\mathrm{f}}\right)\right].$$(8)
The fracture porosity is mainly affected by the aggregation mode of the pore space (characterized by the fractal dimension) and exponentially related to the effective overburden pressure. As for fracture permeability, excepting the aggregation mode of the pore space, the connectivity between fractures (characterized by anomalous diffusion coefficients) also plays an important role (Karimpouli and Tahmasebi, 2016). Pedrosa (1986) found that the permeability and the overburden pressure of the bedrock are exponentially related in the study of unstable test wells in stresssensitive formations. Fractal permeability considering stress sensitivity can be obtained by introducing the anomalous diffusion coefficient based on the definition of fractal porosity:
$${K}_{\mathrm{f}}\left(r\right)={K}_{\mathrm{fw}}\left(\frac{r}{{r}_{\mathrm{w}}}\right){r}^{{D}_{\mathrm{f}}d\theta}\mathrm{exp}\left[{\alpha}_{\mathrm{k}}\left({p}_{\mathrm{i}}{p}_{\mathrm{f}}\right)\right].$$(9)
If we consider the fractal of the fracture in time, it will pose a greater challenge to the calculation of the nonlinear model of horizontal wells in fracturedvuggy tight fractal reservoirs (Albinali et al., 2016; de Swaan, 2016; Noetinger et al., 2001, 2016). The reason why we apply fractal theory to deal with fracture system is to study the initial distribution of fractures in space, so we considered the fractal characteristics of the fracture in space, but not in time.
3 Nonlinear seepage model of the horizontal well in a fracturedvuggy tight fractal reservoir
3.1 Physical model
A physical model of the horizontal well in a circular triplemedium fractal reservoir just as shown in Figure 1 was established, which met the following assumptions:

The model is a threehole, singlepermeability circular reservoir consisting of the matrix system, fracture system and cave system. The quasisteadystate cross flow generates from the matrix system and cave system to the fracture system.

The outer boundary is closed or kept at a constant pressure. The thickness of the reservoir is h, the original pressure is p_{i}, the half length of the horizontal well is L, and the radius of the reservoir is r_{e}.

The model considers the skin factor, wellbore storage coefficient and anisotropy of the permeability. The horizontal and vertical permeability of the fracture system are K_{fh} and K_{fv}, respectively. The permeability of the matrix and cave system are K_{m} and K_{v} respectively.

The horizontal well production is Q, located at the center of the reservoir plane.

The influence of gravity and the capillary force was ignored.

The stress sensitivity, fractal characteristics of the fracture system and the nonlinear seepage mechanism of the matrix were considered.
Fig. 1 Seepage physical model of the horizontal well in a triplemedium fractal reservoir. 
What we study in this paper is a carbonated salt reservoir which is welldeveloped fractures and vugs, the vug in full contact with the fracture is more in line with the actual situation. Lots of current researches on triple medium are based on this model which vugs are in full contact with fractures, ignoring the case that vugs are with limited contact with fracture or within blocks (Kang et al., 2006; Wu et al., 2007). In this study, we focus on singlephase transient flow and fractured vuggy rock, conceptualized also as a tripleor multiplecontinuum medium, consisting of highly permeable fractures, lowpermeability rock matrix, and vugs.
3.2 Mathematical model
In a Cartesian coordinate system, the mass conservation equation of singlephase fluid flow in the model can be combined with the nonlinear motion equation and state equation to obtain a nonlinear seepage model of the horizontal well in a fracturedvuggy carbonate fractal reservoir:
$$\begin{array}{l}\frac{\mathrm{\partial}}{\mathrm{\partial}{x}_{\mathrm{D}}}\left({{r}_{\mathrm{D}}}^{{D}_{\mathrm{f}}d\theta}{e}^{{\alpha}_{\mathrm{KD}}{p}_{\mathrm{fD}}}\frac{\mathrm{\partial}{p}_{\mathrm{fD}}}{\mathrm{\partial}{x}_{\mathrm{D}}}\right)+\frac{\mathrm{\partial}}{\mathrm{\partial}{y}_{\mathrm{D}}}\left({{r}_{\mathrm{D}}}^{{D}_{\mathrm{f}}d\theta}{e}^{{\alpha}_{\mathrm{KD}}{p}_{\mathrm{fD}}}\frac{\mathrm{\partial}{p}_{\mathrm{fD}}}{\mathrm{\partial}{y}_{\mathrm{D}}}\right)+\frac{\mathrm{\partial}}{\mathrm{\partial}{z}_{\mathrm{D}}}\left({{L}_{D}}^{2}{e}^{{\alpha}_{\mathrm{KD}}{p}_{\mathrm{fD}}}\frac{\mathrm{\partial}{p}_{\mathrm{fD}}}{\mathrm{\partial}{z}_{\mathrm{D}}}\right)\\ ={\omega}_{\mathrm{f}}({h}_{\mathrm{D}}{L}_{\mathrm{D}}{)}^{2}\frac{\mathrm{\partial}{p}_{\mathrm{fD}}}{\mathrm{\partial}{t}_{\mathrm{D}}}+{\lambda}_{\mathrm{mf}}{\delta}_{\mathrm{LVDf}}\left({p}_{\mathrm{fD}}{p}_{\mathrm{mD}}\right)+{\lambda}_{\mathrm{vf}}\left({p}_{\mathrm{fD}}{p}_{\mathrm{vD}}\right),\end{array}$$(10)
$${\lambda}_{\mathrm{mf}}{\delta}_{\mathrm{LVDf}}({p}_{\mathrm{fD}}{p}_{\mathrm{mD}})+{\lambda}_{\mathrm{mv}}{\delta}_{\mathrm{LVDv}}({p}_{\mathrm{vD}}{p}_{\mathrm{mD}})=({h}_{\mathrm{D}}{L}_{\mathrm{D}}{)}^{2}\left(1{\omega}_{\mathrm{f}}{\omega}_{\mathrm{v}}\right)\frac{\mathrm{\partial}{p}_{\mathrm{mD}}}{\mathrm{\partial}{t}_{\mathrm{D}}},$$(11)
$${\lambda}_{\mathrm{vf}}({p}_{\mathrm{fD}}{p}_{\mathrm{vD}}){\lambda}_{\mathrm{mv}}{\delta}_{\mathrm{LVDv}}({p}_{\mathrm{vD}}{p}_{\mathrm{mD}})=({h}_{\mathrm{D}}{L}_{\mathrm{D}}{)}^{2}{\omega}_{\mathrm{v}}\frac{\mathrm{\partial}{p}_{\mathrm{vD}}}{\mathrm{\partial}{t}_{\mathrm{D}}},$$(12)
$${\lambda}_{\mathrm{mf}}={\sigma}_{\mathrm{mf}}\frac{{K}_{\mathrm{m}}}{{K}_{\mathrm{fhi}}}{L}^{2};{\lambda}_{\mathrm{mv}}={\sigma}_{\mathrm{mv}}\frac{{K}_{\mathrm{m}}}{{K}_{\mathrm{v}}}{L}^{2};{\lambda}_{\mathrm{vf}}={\sigma}_{\mathrm{vf}}\frac{{K}_{\mathrm{v}}}{{K}_{\mathrm{fhi}}}{L}^{2},$$(13)
where, λ_{vf} is the cross flow factor of cave system to fracture system; λ_{mv} is the cross flow factor of matrix system to cave system; λ_{mf} is the cross flow factor of matrix system to fracture system; σ_{vf}, σ_{mv} and σ_{mf} are the interporosity flow shape factors. The shape factor for m–v or m–f is defined by Warren and Root (1963):
$${\sigma}_{\mathrm{mf}}={\sigma}_{\mathrm{mv}}={\sigma}_{\mathrm{m}}.$$(14)
For v–f interaction, the shape factor for vugs is defined as:
$${\sigma}_{\mathrm{vf}}=\frac{{A}_{\mathrm{vf}}}{{l}_{\mathrm{vf}}},$$(15)
where A_{vf} is the total fracture and vug connection area per unit volume of rock (m^{3}/m^{3}) and l_{vf} is characteristic length.
If (p_{AD} − p_{mD}) < (c_{1D} + c_{2D})/2, δ_{LVDA} = 0, else δ_{LVDA} = 1 − c_{1D}/[2(p_{AD} − p_{mD}) − c_{2D}] (A = f, v)
The initial condition:
$${p}_{\mathrm{lD}}\left({x}_{\mathrm{D}},{y}_{\mathrm{D}},{z}_{\mathrm{D}},{t}_{\mathrm{D}}=0\right)=0\left(l=m,f,v\right).$$(16)
The closed outer boundary condition:
$${\frac{\mathrm{\partial}{p}_{\mathrm{lD}}}{\mathrm{\partial}{r}_{\mathrm{D}}}}_{{r}_{\mathrm{D}}={r}_{\mathrm{eD}}}={\frac{{\mathrm{\partial}p}_{\mathrm{lD}}}{{\mathrm{\partial}z}_{\mathrm{D}}}}_{{z}_{\mathrm{D}}=\mathrm{0,1}}=0.$$(17)
The constant pressure boundary condition:
$${p}_{\mathrm{lD}}\left({r}_{\mathrm{D}}={r}_{\mathrm{eD}}\right)={p}_{\mathrm{lD}}\left({z}_{\mathrm{D}}=\mathrm{0,1}\right)=0.$$(18)
The inner boundary condition:
$$\underset{{\epsilon}_{\mathrm{D}}\to 0}{\mathbf{lim}}{\int}_{{z}_{\mathrm{wD}}\frac{{\epsilon}_{\mathrm{D}}}{2}}^{{z}_{\mathrm{wD}}+\frac{{\epsilon}_{\mathrm{D}}}{2}}\left({{r}_{\mathrm{D}}}^{{D}_{\mathrm{f}}\mathrm{d}\theta}{e}^{{\alpha}_{\mathrm{KD}}{p}_{\mathrm{fD}}}\frac{\partial {p}_{\mathrm{fD}}}{\partial {r}_{\mathrm{D}}}\right)\mathrm{d}{z}_{\mathrm{wD}}{}_{{r}_{D}=\epsilon D}=\frac{1}{2},\left{z}_{\mathrm{D}}{z}_{\mathrm{wD}}\right\le \frac{{\epsilon}_{\mathrm{D}}}{2}.$$(19)
Each dimensionless quantity is defined as follows:
$${L}_{\mathrm{D}}=\frac{L}{h}\sqrt{\frac{{K}_{\mathrm{fvi}}}{{K}_{\mathrm{fhi}}}};{h}_{\mathrm{D}}=\frac{h}{{r}_{\mathrm{w}}}\sqrt{\frac{{K}_{\mathrm{fhi}}}{{K}_{\mathrm{fvi}}}};{x}_{\mathrm{D}}=\frac{x}{L};{y}_{\mathrm{D}}=\frac{y}{L};{z}_{\mathrm{D}}=\frac{z}{h};{\alpha}_{\mathrm{KD}}=\frac{1.842\times 1{0}^{3}\mathrm{Q\mu}B}{{K}_{\mathrm{fhi}}h}{\alpha}_{\mathrm{K}};{\omega}_{\mathrm{f}}=\frac{{\varphi}_{\mathrm{f}}{C}_{\mathrm{ft}}}{{\varphi}_{\mathrm{f}}{C}_{\mathrm{ft}}+{\varphi}_{\mathrm{m}}{C}_{\mathrm{mt}}+{\varphi}_{\mathrm{v}}{C}_{\mathrm{vt}}};{z}_{\mathrm{wD}}=\frac{{z}_{\mathrm{w}}}{h};{\omega}_{\mathrm{v}}=\frac{{\varphi}_{\mathrm{v}}{C}_{\mathrm{vt}}}{{\varphi}_{\mathrm{f}}{C}_{\mathrm{ft}}+{\varphi}_{\mathrm{m}}{C}_{\mathrm{mt}}+{\varphi}_{\mathrm{v}}{C}_{\mathrm{vt}}};{c}_{\mathrm{D}}=\frac{{K}_{\mathrm{fhi}}h{l}_{\mathrm{c}}}{1.842\times 1{0}^{3}\mathrm{Q\mu}B}c;{\delta}_{\mathrm{LVDA}}=1\frac{{c}_{\mathrm{D}}}{2\left({p}_{\mathrm{AD}}{p}_{\mathrm{mD}}\right)+{c}_{D}}\left(A=f,v\right);{C}_{\mathrm{D}}=\frac{0.1592C}{h\varphi {c}_{\mathrm{t}}{{r}_{\mathrm{w}}}^{2}};{r}_{\mathrm{D}}=\frac{r}{L}{p}_{\mathrm{lD}}=\frac{{K}_{\mathrm{fhi}}h\left({p}_{\mathrm{i}}{p}_{\mathrm{l}}\right)}{1.842\times 1{0}^{3}\mathrm{Q\mu}B}\left(l=m,f,v\right);{t}_{\mathrm{D}}=\frac{{K}_{\mathrm{fhi}}t}{({\varphi}_{\mathrm{f}}{C}_{\mathrm{ft}}+{\varphi}_{\mathrm{m}}{C}_{\mathrm{mt}}+{\varphi}_{\mathrm{v}}{C}_{\mathrm{vt}})\mu {{r}_{\mathrm{w}}}^{2}}.$$
3.3 Model solving
Analytic solutions are highly demanding on mathematics and rely heavily on idealized assumptions which are actually very different from the actual engineering problems, and once the engineering problem is slightly more complicated, we cannot directly get the analytical solution, or the analytical solution error is too large.
The finite element method is used to solve the problem of Partial Differential Equations (PDE) with specific boundary conditions in mechanics and mathematics. It requires a large scale of equations, resulting in a large amount of input data and preparation for calculation. The float of the accuracy based on the method is relatively large. The finer the unit division, the more accurate the approximation result, accompanied by a large increase in the amount of calculation. However, the finite element method can simulate infinite complexes with finite, interrelated elements, no matter how complex the geometry can be simplified with the corresponding elements. In that way, results could be calculated and complex engineering problems could be simplified through model analysis.
The process of finite element solution can be described as: total structure discretization – unit mechanics analysis – unit assembly – total structure analysis – application of boundary conditions – total structure response – reaction analysis of a unit inside the structure. The fractal nonlinear seepage model established in the paper is complex and suitable for numerical solution rather than analytical solution.
Considering that the infinite conductivity model is more difficult to solve than the uniform traffic model, this paper selects 0.7 L as the equivalent pressure point under the two models, so that the uniform traffic model can be used to evaluate the bottom hole pressure.
The horizontal well is considered as the sourcesink term of the unit and being integrated:
$${f}_{\mathrm{e}}=\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{N}_{\mathrm{i}}q\mathrm{d}V.$$(20)
The horizontal well is divided into n nodes and the finite element formula of the element sourcesink term can be obtained by applying the Delta function:
$${f}_{\mathrm{e}}=\frac{Q}{n}\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{N}_{\mathrm{i}}\delta \left(x{x}_{0}\right)\delta \left(y{y}_{0}\right)\delta \left(z{z}_{0}\right)\mathrm{d}V=\frac{Q}{n}\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{N}_{\mathrm{i}}\left({x}_{0},{y}_{0},{z}_{0}\right)\mathrm{d}V.$$(21)
The dimensionless form of the formula (18) is:
$${f}_{\mathrm{e}}=\frac{2\pi}{n}\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\overset{}{\iiint}}{N}_{\mathrm{i}}\left({x}_{D0},{y}_{D0},{z}_{D0}\right)\mathrm{d}V.$$(22)
For the fracture system, the finite element equation of the fracture system can be obtained by applying the Galerkin method:
$$\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{N}_{\mathrm{i}}\left[\frac{\partial}{\partial {x}_{\mathrm{D}}}\left({{r}_{\mathrm{D}}}^{{D}_{\mathrm{f}}d\theta}{e}^{{\alpha}_{\mathrm{KD}}{p}_{\mathrm{fD}}}\frac{\mathrm{\partial}{p}_{\mathrm{fD}}}{\mathrm{\partial}{x}_{\mathrm{D}}}\right)+\frac{\mathrm{\partial}}{\mathrm{\partial}{y}_{\mathrm{D}}}\left({{r}_{\mathrm{D}}}^{{D}_{\mathrm{f}}d\theta}{e}^{{\alpha}_{\mathrm{KD}}{p}_{\mathrm{fD}}}\frac{\mathrm{\partial}{p}_{\mathrm{fD}}}{\mathrm{\partial}{y}_{\mathrm{D}}}\right)\right]\mathrm{d}V+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{N}_{\mathrm{i}}\left[\frac{\mathrm{\partial}}{\mathrm{\partial}{z}_{\mathrm{D}}}\left({e}^{{\alpha}_{\mathrm{KD}}{p}_{\mathrm{fD}}}{{L}_{D}}^{2}\frac{\mathrm{\partial}{p}_{\mathrm{fD}}}{\mathrm{\partial}{z}_{\mathrm{D}}}\right)\right]\mathrm{d}V=\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{N}_{\mathrm{i}}{\left({L}_{\mathrm{D}}{h}_{\mathrm{D}}\right)}^{2}{\omega}_{\mathrm{f}}\frac{\mathrm{\partial}{p}_{\mathrm{fD}}}{\mathrm{\partial}{t}_{\mathrm{D}}}\mathrm{d}V+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{N}_{\mathrm{i}}{\delta}_{\mathrm{LVDf}}{\lambda}_{\mathrm{mf}}({p}_{\mathrm{fD}}{p}_{\mathrm{mD}})\mathrm{d}V+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{N}_{\mathrm{i}}{\lambda}_{\mathrm{vf}}({p}_{\mathrm{fD}}{p}_{\mathrm{vD}})\mathrm{d}V.$$(23)
On the basis of Green’s theorem, the finite element equation of the inner element and the outer boundary element under the closed condition can be obtained through integration by parts:
$$\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{{r}_{\mathrm{D}}}^{{D}_{\mathrm{f}}d\theta}{e}^{{\alpha}_{\mathrm{KD}}{p}_{\mathrm{fD}}}\left(\frac{\mathrm{\partial}{N}_{\mathrm{i}}}{\mathrm{\partial}{x}_{\mathrm{D}}}\frac{\mathrm{\partial}{p}_{\mathrm{fD}}}{\mathrm{\partial}{x}_{\mathrm{D}}}+\frac{\mathrm{\partial}{N}_{\mathrm{i}}}{\mathrm{\partial}{y}_{\mathrm{D}}}\frac{\mathrm{\partial}{p}_{\mathrm{fD}}}{\mathrm{\partial}{y}_{\mathrm{D}}}\right)\mathrm{d}V+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{e}^{{\alpha}_{\mathrm{KD}}{p}_{\mathrm{fD}}}\left({{L}_{D}}^{2}\frac{\partial {N}_{i}}{\partial {z}_{D}}\frac{\partial {p}_{\mathrm{fD}}}{\partial {z}_{D}}\right)\mathrm{d}V+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{N}_{\mathrm{i}}{\delta}_{\mathrm{LVDf}}{\lambda}_{\mathrm{mf}}{p}_{\mathrm{fD}}\mathrm{d}V+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{N}_{\mathrm{i}}{\lambda}_{\mathrm{vf}}{p}_{\mathrm{fD}}\mathrm{d}V+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{N}_{\mathrm{i}}{\left({L}_{\mathrm{D}}{h}_{\mathrm{D}}\right)}^{2}{\omega}_{\mathrm{f}}\frac{\mathrm{\partial}{p}_{\mathrm{fD}}}{\mathrm{\partial}{t}_{\mathrm{D}}}\mathrm{d}V=\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{N}_{\mathrm{i}}{\delta}_{\mathrm{LVDf}}{\lambda}_{\mathrm{mf}}{p}_{\mathrm{mD}}\mathrm{d}V+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{N}_{\mathrm{i}}{\lambda}_{\mathrm{vf}}{p}_{\mathrm{vD}}\mathrm{d}V.$$(24)
Thus, the matrix form of the finite element equation is:
$$\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}\left[\mathbf{C}{\mathbf{B}}^{\mathbf{T}}\mathbf{GB}+{\lambda}_{\mathrm{mf}}\mathbf{DN}{\mathbf{N}}^{\mathbf{T}}+{\lambda}_{\mathrm{vf}}\mathbf{EN}{\mathbf{N}}^{\mathbf{T}}+{\omega}_{\mathrm{f}}({L}_{\mathrm{D}}{h}_{\mathrm{D}}{)}^{2}\frac{\mathbf{N}{\mathbf{N}}^{\mathbf{T}}}{\Delta {t}_{\mathrm{D}}}\right]\mathrm{d}V{{\mathit{P}}_{f}}^{n}=\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{\omega}_{\mathrm{f}}({L}_{\mathrm{D}}{h}_{\mathrm{D}}{)}^{2}\frac{\mathbf{N}{\mathbf{N}}^{\mathbf{T}}}{\Delta {t}_{\mathrm{D}}}\mathrm{d}V{{\mathit{P}}_{\mathrm{f}}}^{n1}+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{\lambda}_{\mathrm{mf}}\mathbf{DN}{\mathbf{N}}^{\mathbf{T}}\mathrm{d}V{{\mathit{P}}_{\mathrm{m}}}^{n1}+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{\lambda}_{\mathrm{vf}}\mathbf{EN}{\mathbf{N}}^{\mathbf{T}}\mathrm{d}V{{\mathit{P}}_{\mathrm{v}}}^{n1}.$$(25)
Assuming E is an identity matrix, the expressions of other matrices are as follows:
$$\mathit{B}=\left[\begin{array}{l}\frac{\mathrm{\partial}{N}_{1}}{\mathrm{\partial}{x}_{\mathrm{D}}}\frac{\mathrm{\partial}{N}_{2}}{\mathrm{\partial}{x}_{\mathrm{D}}}...\frac{\mathrm{\partial}{N}_{n}}{\mathrm{\partial}{x}_{\mathrm{D}}}\\ \frac{\mathrm{\partial}{N}_{1}}{\mathrm{\partial}{y}_{\mathrm{D}}}\frac{\mathrm{\partial}{N}_{2}}{\mathrm{\partial}{y}_{\mathrm{D}}}...\frac{\mathrm{\partial}{N}_{n}}{\mathrm{\partial}{y}_{\mathrm{D}}}\\ \frac{\mathrm{\partial}{N}_{1}}{\mathrm{\partial}{z}_{\mathrm{D}}}\frac{\mathrm{\partial}{N}_{2}}{\mathrm{\partial}{z}_{\mathrm{D}}}...\frac{\mathrm{\partial}{N}_{n}}{\mathrm{\partial}{z}_{\mathrm{D}}}\end{array}\right],\mathit{N}=\left[\begin{array}{c}{N}_{1}\\ {N}_{2}\\ .\\ .\\ .\\ {N}_{n}\end{array}\right],{\mathit{P}}_{\mathrm{f}}=\left[\begin{array}{c}{p}_{\mathrm{fD}1}\\ {p}_{\mathrm{fD}2}\\ .\\ .\\ .\\ {p}_{\mathrm{fD}n}\end{array}\right],{\mathit{P}}_{\mathrm{m}}=\left[\begin{array}{c}{p}_{\mathrm{mD}1}\\ {p}_{\mathrm{mD}2}\\ .\\ .\\ .\\ {p}_{\mathrm{mD}n}\end{array}\right],{\mathit{P}}_{\mathrm{v}}=\left[\begin{array}{c}{p}_{\mathrm{vD}1}\\ {p}_{\mathrm{vD}2}\\ .\\ .\\ .\\ {p}_{\mathrm{vD}n}\end{array}\right]\mathit{C}=\left[\begin{array}{ccc}{e}^{{\alpha}_{\mathrm{KD}}{p}_{\mathrm{fD}1}}& & \begin{array}{cc}& \end{array}\\ & {e}^{{\alpha}_{\mathrm{KD}}{p}_{\mathrm{fD}2}}& \begin{array}{cc}& \end{array}\\ \begin{array}{c}\\ \end{array}& \begin{array}{c}\\ \end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \end{array}& \begin{array}{c}\\ {e}^{{\alpha}_{\mathrm{KD}}{p}_{\mathrm{fD}n}}\end{array}\end{array}\end{array}\right],{\mathit{H}}_{\mathrm{A}}=\left[\begin{array}{ccc}{\delta}_{\mathrm{LVDA}1}& & \begin{array}{cc}& \end{array}\\ & {\delta}_{\mathrm{LVDA}2}& \begin{array}{cc}& \end{array}\\ \begin{array}{c}\\ \end{array}& \begin{array}{c}\\ \end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \end{array}& \begin{array}{c}\\ {\delta}_{\mathrm{LVDAn}}\end{array}\end{array}\end{array}\right](A=f,v),\mathbf{G}=\left[\begin{array}{ccc}{{r}_{\mathrm{D}}}^{{D}_{\mathrm{f}}d\theta}& & \\ & {{r}_{\mathrm{D}}}^{{D}_{\mathrm{f}}d\theta}& \\ & & {{L}_{\mathrm{D}}}^{2}\end{array}\right]$$
Further simplification of the matrix equation can be obtained:
$${\mathit{K}}_{\mathrm{e}}{{\mathit{P}}_{\mathrm{f}}}^{n}={\mathit{F}}_{\mathrm{e}}$$(26)
where, the expression forms of K_{e} and F_{e} are followed as:
$${\mathit{K}}_{\mathrm{e}}=\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}\left[\mathbf{C}{\mathbf{B}}^{\mathbf{T}}\mathbf{GB}+{\lambda}_{\mathrm{mf}}{\mathit{H}}_{\mathrm{f}}\mathbf{N}{\mathbf{N}}^{\mathbf{T}}+{\lambda}_{\mathrm{vf}}\mathbf{EN}{\mathbf{N}}^{\mathbf{T}}+{\omega}_{\mathrm{f}}({L}_{\mathrm{D}}{h}_{\mathrm{D}}{)}^{2}\frac{\mathbf{N}{\mathbf{N}}^{\mathbf{T}}}{\Delta {t}_{\mathrm{D}}}\right]\mathrm{d}V,$$(27)
$${\mathit{F}}_{\mathrm{e}}=\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{\omega}_{\mathrm{f}}({L}_{\mathrm{D}}{h}_{\mathrm{D}}{)}^{2}\frac{\mathbf{N}{\mathbf{N}}^{\mathbf{T}}}{\Delta {t}_{\mathrm{D}}}\mathrm{d}V{{\mathit{P}}_{\mathrm{f}}}^{n1}+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{\lambda}_{\mathrm{mf}}{\mathit{H}}_{\mathrm{f}}\mathbf{N}{\mathbf{N}}^{\mathbf{T}}\mathrm{d}V{{\mathit{P}}_{\mathrm{m}}}^{n1}+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{\lambda}_{\mathrm{vf}}\mathbf{EN}{\mathbf{N}}^{\mathbf{T}}\mathrm{d}V{{\mathit{P}}_{\mathrm{v}}}^{n1}.$$(28)
Similarly, the finite element equation of the matrix system can be expressed as:
$${\mathit{K}}_{\mathrm{e}}{{\mathit{P}}_{m}}^{n}={\mathit{F}}_{\mathrm{e}},$$(29)
where, the expression forms of K _{ e } and F _{ e } are followed as:
$${\mathit{K}}_{\mathrm{e}}=\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}\left[{\lambda}_{\mathrm{mf}}{\mathit{H}}_{\mathrm{f}}\mathbf{N}{\mathbf{N}}^{\mathbf{T}}+{\lambda}_{\mathrm{mv}}{\mathit{H}}_{\mathrm{v}}\mathbf{N}{\mathbf{N}}^{\mathbf{T}}+(1{\omega}_{\mathrm{f}}{\omega}_{\mathrm{v}})({L}_{\mathrm{D}}{h}_{\mathrm{D}}{)}^{2}\frac{\mathbf{N}{\mathbf{N}}^{\mathbf{T}}}{\Delta {t}_{\mathrm{D}}}\right]\mathrm{d}V,$$(30)
$${\mathit{F}}_{\mathrm{e}}=\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}(1{\omega}_{\mathrm{f}}{\omega}_{\mathrm{v}})({L}_{\mathrm{D}}{h}_{\mathrm{D}}{)}^{2}\frac{\mathbf{N}{\mathbf{N}}^{\mathbf{T}}}{\Delta {t}_{\mathrm{D}}}\mathrm{d}V{\mathit{P}}_{\mathrm{m}}^{n1}+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{\lambda}_{\mathrm{mf}}{\mathit{H}}_{\mathrm{f}}\mathbf{N}{\mathbf{N}}^{\mathbf{T}}\mathrm{d}V{{\mathit{P}}_{f}}^{n}+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{\lambda}_{\mathrm{mv}}{\mathit{H}}_{\mathrm{v}}\mathbf{N}{\mathbf{N}}^{\mathbf{T}}\mathrm{d}V{{\mathit{P}}_{f}}^{n}.$$(31)
The finite element equation of the cave system can be expressed as:
$${\mathit{K}}_{\mathrm{e}}{{\mathit{P}}_{v}}^{n}={\mathit{F}}_{\mathrm{e}},$$(32)
where, the expression forms of K _{e} and F _{e} are as follows:
$${\mathit{K}}_{\mathrm{e}}=\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}\left[{\lambda}_{\mathrm{vf}}\mathbf{EN}{\mathbf{N}}^{\mathbf{T}}+{\lambda}_{\mathrm{mv}}{\mathit{H}}_{\mathrm{v}}\mathbf{N}{\mathbf{N}}^{\mathbf{T}}+{\omega}_{\mathrm{v}}({L}_{\mathrm{D}}{h}_{\mathrm{D}}{)}^{2}\frac{\mathbf{N}{\mathbf{N}}^{\mathbf{T}}}{\Delta {t}_{\mathrm{D}}}\right]\mathrm{d}V,$$(33)
$${\mathit{F}}_{\mathrm{e}}=\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{\omega}_{\mathrm{v}}({L}_{\mathrm{D}}{h}_{\mathrm{D}}{)}^{2}\frac{\mathbf{N}{\mathbf{N}}^{\mathbf{T}}}{\Delta {t}_{\mathrm{D}}}\mathrm{d}V{{\mathit{P}}_{\mathrm{v}}}^{n1}+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{\lambda}_{\mathrm{mv}}{\mathit{H}}_{\mathrm{v}}\mathbf{N}{\mathbf{N}}^{\mathbf{T}}\mathrm{d}V{\mathit{P}}_{\mathrm{m}}^{n}+\underset{{\mathrm{\Omega}}_{\mathrm{e}}}{\iiint}{\lambda}_{\mathrm{vf}}\mathbf{EN}{\mathbf{N}}^{\mathbf{T}}\mathrm{d}V{\mathit{P}}_{\mathrm{f}}^{n}.$$(34)
Equations (24), (25), and (29), respectively, denote the finite element balance equation of the fracture system, matrix system, and cave system. Where K _{ e } is the nonlinear coefficient matrix, F _{ e } is the unit load vector, the superscript n represents the time variable. The solution process acquires giving the initial value firstly, then conducting linear treatment to obtain the pressure solution.
The wellbore storage coefficient and skin factor are taken into consideration through the Duhamel’s principle.
$${\stackrel{\u0305}{p}}_{\mathrm{wD}}=\frac{s{\stackrel{\u0305}{p}}_{\mathrm{D}}+\mathrm{S}{L}_{\mathrm{D}}}{s\left[1+{C}_{\mathrm{D}}s\left(s{\stackrel{\u0305}{p}}_{\mathrm{D}}+\mathrm{S}{L}_{\mathrm{D}}\right)\right]}.$$(35)
Then, the final form of the bottom hole pressure solution is obtained after conducting the Stehfest numerical inversion:
$${p}_{\mathrm{wD}}=\frac{\mathrm{ln}2}{{t}_{\mathrm{D}}}\sum _{i=1}^{N}{V}_{i}{\stackrel{\u0305}{p}}_{\mathrm{wD}}.$$(36)
In equation (33): ${V}_{i}=(1{)}^{\left(\frac{N}{2}+i\right)}\sum _{k=\frac{i+1}{2}}^{\mathrm{min}\left(i,\frac{N}{2}\right)}\frac{{k}^{\frac{N}{2}}\left(2k+1\right)!}{\left(k+1\right)!k!\left(\frac{N}{2}k+1\right)!\left(ik+1\right)!\left(2ki+1\right)!}$.
4 Model comparison and law summary
The transient pressure and its derivative curves of Darcy model, nonlinear model, fractal model and fractal nonlinear model of triplemedium lowpermeability reservoir are shown in Figure 2. This paper sets nonlinear parameter c_{D} = c_{1D} = c_{2D}, fractal index β = 2 + θ − D_{f}. The values of the known parameters are as follows: C_{D} = 100, S = 0, ω_{f} = 0.1, ω_{v} = 0.5, λ_{m} = 0.001, λ_{v} = 0.1, L_{D} = 5, h_{D} = 200, α_{KD} = 0.02. Darcy model (c_{D} = 0, θ = 0, D_{f} = 2); Nonlinear flow model (c_{D} = 0.3, θ = 0, D_{f} = 2); Fractal nonlinear model (c_{D} = 0.3, θ = 0.3, D_{f} = 1.9). Under the condition of considering the stress sensitivity, each model can adjust the value of the stress sensitivity coefficient α_{KD} according to the actual situation.
Fig. 2 Comparison of well test curves of different models. 
The value of the cross flow factor determines the time of occurrence of the concave section while the value of the elastic storativity ratio determines the depth and width of the concave. If vugs are not considered in the model, then there is only one concave on the pressure derivative curve. In general, the cross flow factor of cave system is larger, so the cross flow stage of cave system to fracture system appears earlier.
The research shows that there are still nine type flow stages in the pressure and its derivative curves even under different models, the specific flow stage division and the comparison of different models are shown in Table 1. The fractal nonlinear model contains all the characteristics of the fractal model and nonlinear model, so the model is more complicated and more factors need to be discussed (Kazemi, 1969).
Flow state division and model comparison.
5 Sensitivity analysis
Firstly, the sensitivity analysis of the nonlinear parameter and fractal index are conducted. The values of the known parameters are as follows: C_{D} = 100, S = 0, ω_{f} = 0.1, ω_{v} = 0.5, λ_{m} = 0.001, λ_{v} = 0.1, L_{D} = 5, h_{D} = 200, α_{KD} = 0.02.
Figure 3 shows that the nonlinear parameter mainly affects the cross flow stage of the matrix system to the fracture system. The parameter decides the ability of cross flow from the matrix system to the fracture system and it can cause a later appearance of the concave of the derivative curve. In other words, the bigger the nonlinear parameter is, the shorter time the cross flow lasts, the narrower and shallower the Vshape concave of the pressure derivative curve. Even so, the nonlinear parameter does not affect the end time of the cross flow. It is obviously shown that the variation of the depth of the concave gets smaller gradually with the increasing nonlinear parameter. Figure 4 indicates that the influence of the fractal index on the pressure and derivative curves are mainly reflected in the stages after the cross flow stage from cave system to fracture system. Starting from this stage, both the pressure and derivative curves turn upward gradually with the increasing time and show a steeper upwarping if the fractal index becomes greater. However, the fractal index does not affect the time of appearance and the duration of the two cross flow stages, so the depth of the concaves of the derivative curves remain unchanged. The fractal parameter also exerts little influence on the early radial flow and early linear flow, and the larger the fractal parameter is, the lower the pressure and derivative curves of the two flow stages are.
Fig. 3 Influence of the nonlinear parameter on well test curves. 
Fig. 4 Influence of the fractal index on well test curves. 
Then, the sensitivity analysis of other factors are conducted and the values of the known parameters are as follows: C_{D} = 100, S = 0, ω_{f} = 0.01, ω_{v} = 0.1, λ_{m} = 0.001, λ_{v} = 0.1, L_{D} = 5, h_{D} = 200, α_{KD} = 0.02.
Figure 5 reflects that if the permeability modulus increases, the upwarping amplitude of the pressure and derivative curves of the stages after the linear flow increase because the stronger the stress sensitivity of the reservoir is, the more difficult the fluid flows, and the larger drawdown pressure is needed. Both the pressure and derivative curves turn upward gradually with the increasing of time. When the value of the permeability modulus is greater than 0.1, the upwarping amplitude increases sharply. Figure 6 shows the influence of the horizontal well position on the well test curves. When other parameters remain unchanged, if the value of z_{wD} is close to 0.5, that is to say, the horizontal well is in the middle of the reservoir, the radial flow stage occurs earlier and the time of the early linear flow stage becomes longer. However, the location of the horizontal well does not influence other flow stages.
Fig. 5 Influence of the permeability modulus on well test curves. 
Fig. 6 Influence of the horizontal well location on well test curves. 
Figure 7 shows that the flow stages after the early linear flow would be delayed because the longer the horizontal well is, the lower the pressure and derivative curves are, the longer time the early radial flow and lasts. Consequently the early linear flow regime will last longer meanwhile the appearance of the next flow periods will be postponed. It is obviously shown that the early radial flow will be covered by the wellbore storage effect if the horizontal well length is short enough. Considering the fractal characteristics and stress sensitivity of the fracture system, the pressure or derivative curves of different horizontal well lengths manifest as several parallel lines instead of converging into a single line in the late radial flow stage.
Fig. 7 Influence of the horizontal well length on well test curves. 
6 Example verification
The standard Genetic Algorithm (GA) is used to fit the measured pressure curve and output the fitting value of the parameter. First, the evolutionary algebra counter, the maximum evolution algebra and the range of parameter values are set to prepare an initial population of several initial solutions for genetic manipulation. Then, the evaluation function value (adaptation) is used to evaluate the merits of the solution or the individual and it is served as the basis for subsequent genetic operations. In addition, selection, replication and mutation operations are used to simulate the natural selection and genetic mechanism of the biological evolution (Yu et al., 2015). Finally, the optimal solution is searched by simulating the natural evolution process, the theoretical curve and the measured curve are drawn and the fitted values are output. The specific processes are as shown in Figure 8. S63 is a production well of X oilfield which was tested in September 2005. The fractal nonlinear model was selected to conduct the well test interpretation and the known parameters of the model are as shown in Table 2.
Fig. 8 Well test fitting flow chart. 
Basic formation parameters.
The analysis of the pressure derivative curve of the well S63 (Fig. 9) shows that the wellbore storage coefficient changes from small to large during well testing, resulting in a poor fit between the actual and theoretical pressure derivative curves in the initial flow stage. The pressure derivative curve is separated from the pressure curve and its slope is less than 1 in wellbore storage stage. In the middle and late flow stages, the pressure derivative curve of S63 well generates two “concave” which clearly reflects the characteristics of the triplemedium reservoir. The actual value fits well with the calculated value of the fractal nonlinear model accept the wellbore storage stage although the well test data is incomplete. The seepage parameters explained by the model agree with the actual situation of the well which are shown in Table 3, indicating that the fractal nonlinear model can be used for well test interpretation and seepage law analysis without considering the variation of wellbore storage coefficient.
Fig. 9 Well test model fitting. 
Well test parameter interpretation.
7 Conclusion

A nonlinear seepage model of a horizontal well in a triplemedium fractal reservoir was derived by orderly adopting a series of treatments, including proposing a new nonlinear seepage model to characterize the seepage situation of matrix blocks and applying the fractal theory to describe the stress sensitivity and fractal characteristics of the fracture system.

To divide the flow stage and summarize the characteristics of different seepage flow models based on the bottom hole pressure curves, the processing of the internal boundary condition was conducted and the finite element principle was applied to solve the models.

The nine stages were summarized in the seepage process namely pure wellbore storage stage, skin effect transition stage, early radial flow, early linear flow, midterm radial flow, cross flow stage of cave system to fracture system, pseudo radial flow of matrix system and cave system, cross flow stage of matrix system to fracture system and late radial flow.

Typical factors influencing the transient pressure performance were analyzed to better understand the flow characteristics and the systematic summary was made based on the sensitivity analysis of the factors.

The triplemedium fractal nonlinear model was applied to conduct the well test interpretation of an actual reservoir. The field application indicates that the measured value and the theoretical calculation value are better fitted if the variation of wellbore storage coefficient ignored. Since the explained seepage parameters are in line with the actual situation, it’s proved that the triplemedium fractal nonlinear model can be used to conduct well test interpretation and seepage law analysis.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant No. 51574265). The authors would like to acknowledge Reviewers and Editors whose critical comments will be very helpful in preparing this article.
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]
 Albinali A., Holy R., Sarak H., Ozkan E. (2016) Modeling of 1D anomalous diffusion in fractured nanoporous media, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 71, 56. [CrossRef] [Google Scholar]
 Aprilian S., Abdassah D., Mucharam L., Sumantri R. (1993) Application of fractal reservoir model for interference test analysis in Kamojang geothermal field (Indonesia), SPE Annual Technical Conference and Exhibition, 3–6, October, Houston, Texas, Society of Petroleum Engineers. [Google Scholar]
 CamachoVelazquez R., VasquezCruz M., CastrejonAivar R., AranaOrtiz V. (2002) Pressure transient and decline curve behaviors in naturally fractured vuggy carbonate reservoirs, SPE Annual Technical Conference and Exhibition, 29 September–2 October, San Antonio, Texas, Society of Petroleum Engineers. [Google Scholar]
 Cao R.Y., Wang Y., Cheng L.S., Ma Y.Z., Tian X.F., An N. (2016) A new model for determining the effective permeability of tight formation, Transp. Porous Media 112, 1, 21–37. [CrossRef] [Google Scholar]
 Celata G.P., Cumo M., McPhail S., Zummo G. (2006) Characterization of fluid dynamic behaviour and channel wall effects in microtube, Int. J. Heat Fluid Flow 27, 1, 135–143. [CrossRef] [Google Scholar]
 Chang J., Yortsos Y.C. (1990) Pressure transient analysis of fractal reservoirs, SPE Form. Eval. 5, 1, 31–38. [CrossRef] [Google Scholar]
 Chen X.J., Yao G.Q. (2017) An improved model for permeability estimation in low permeable porous media based on fractal geometry and modified HagenPoiseuille flow, Fuel 210, 748–757. [CrossRef] [Google Scholar]
 Cossio M., Moridis G.J., Blasingame T.A. (2013) A semianalytic solution for flow in finiteconductivity vertical fractures by use of fractal theory, SPE J. 18, 1, 83–96. [CrossRef] [Google Scholar]
 de Swaan A. (2016) Pressure transients in a fractalcluster model of porous media, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 71, 9. [CrossRef] [Google Scholar]
 Ezulike D.O., Dehghanpour H. (2013) Characterizing tight oil reservoirs using dual and tripleporosity models, SPE Unconventional Resources Conference Canada, 5–7, November, Calgary, Alberta, Canada, Society of Petroleum Engineers. [Google Scholar]
 Ge J.L., Wu Y.S. (1982) The behavior of naturally fractured reservoirs and the technique for well test analysis at constant pressure conduction, Petrol. Explor. Dev. 9, 3, 53–64. [Google Scholar]
 Gomez S., FuentesCruz G., CamachoVelazquez R., VasquezCruz M.A., Otero J., Mesejo A., Castilo N.D. (2006) Application of an evolutionary algorithm in well test characterization of naturally fractured vuggy reservoirs, International Oil Conference and Exhibition in Mexico, 31 August–2 September, Cancun, Mexico, Society of Petroleum Engineers. [Google Scholar]
 Guo J.J., Wang H.T., Zhang L.H., Li C.Y. (2015) Pressure transient analysis and flux distribution for multistage fractured horizontal wells in tripleporosity reservoir media with consideration of stresssensitivity effect, J. Chem. 2015, 1–16. [CrossRef] [Google Scholar]
 Jiang R.Z., Li L.K., Xu J.C., Yang R.F., Zhuang Y. (2012) A nonlinear mathematical model for lowpermeability reservoirs and welltesting analysis, Acta Pet. Sin. 33, 2, 264–268. [Google Scholar]
 Kang Z.J., Wu Y.S., Li J.L., Wu Y.C., Zhang J., Wang G.F. (2006) Modeling multiphase flow in naturally fractured vuggy petroleum reservoirs, SPE Annual Technical Conference and Exhibition, 24–27 September, San Antonio, Texas, Society of Petroleum Engineers. [Google Scholar]
 Karimpouli S., Tahmasebi P. (2016) A hierarchical sampling for capturing permeability trend in rock physics, Transp. Porous Media 116, 3, 1057–1072. [CrossRef] [Google Scholar]
 Kazemi H. (1969) Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution, SPE J. 9, 4, 451–462. [Google Scholar]
 Kong X.Y., Li D.L., Lu D.T. (2008) Transient pressure analysis in porous and fractured fractal reservoirs, Sci. China Ser. E: Technol. Sci. 52, 9, 2700–2708. [Google Scholar]
 Li Y., Wang Q., Li B.Z., Liu Z.L. (2017) Dynamic characterization of different reservoir types for a fracturedcaved carbonate reservoir, SPE Kingdom of Saudi Arabia Annual Technical Symposium and Exhibition, 24–27 April, Dammam, Saudi Arabia, Society of Petroleum Engineers. [Google Scholar]
 Lin Y.Y., Myers M.T. (2018) Impact of nonlinear transport properties on low permeability measurements, J. Nat. Gas Sci. Eng. 54, 328–341. [CrossRef] [Google Scholar]
 Lu H.L., Zhang J.J., Tan X.H., Li X.P., Zhao J.H. (2017) A fractal analysis for net present value of multistage hydraulic fractured horizontal well, Fractals 25, 4. [Google Scholar]
 Nie R.S., Meng Y.F., Yang Z.Z., Guo J.C., Jia Y.L. (2011) New flow model for the triple media carbonate reservoir, Int. J. Comput. Fluid Dyn. 25, 2, 95–104. [CrossRef] [Google Scholar]
 Noetinger B., Estebenet T., Landereau P. (2001) A direct determination of the transient exchange term of fractured media using a continuous time random walk method, Trans. Porous Media 44, 3, 539–557. [Google Scholar]
 Noetinger B., Roubinet D., Russian A., Borgne T.L., Delay F., Dentz M., de Dreuzy J.R., Gouze P. (2016) Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale, Tran. Porous Media 115, 2, 345–385. [Google Scholar]
 Pedrosa Jr O.A. (1986) Pressure transient response in stresssensitive formations, SPE California Regional Meeting, 2–4 April, Oakland, California, Society of Petroleum Engineers. [Google Scholar]
 Popov P., Qin G., Bi L., Efendiev Y., Ewing R.E., Li J. (2009) Multiphysics and multiscale methods for modeling fluid flow through naturally fractured carbonate karst reservoirs, SPE Reserv. Eval. Eng. 122, 218–231. [CrossRef] [Google Scholar]
 Razminia K., Razminia A., Trujilo J.J. (2015) Analysis of radial composite systems based on fractal theory and fractional calculus, Sig. Processing 107, 378–388. [CrossRef] [Google Scholar]
 Vijayalakshmi K., Anoop K.B., Patel H.E., Harikrishna P.V., Sundararajan T., Das S.K. (2009) Effects of compressibility and transition to turbulence on flow through microchannels, Int. J. Heat Mass Trans. 52, 9–10, 2196–2204. [CrossRef] [Google Scholar]
 Wang W.D., Su Y.L., Sheng G.L., Cossio M., Shang Y.Y. (2015) A mathematical model considering complex fractures and fractal flow for pressure transient analysis of fractured horizontal wells in unconventional reservoirs, J. Nat. Gas Sci. Eng. 23, 139–147. [CrossRef] [Google Scholar]
 Wang W.D., Yuan B., Su Y.L., Sheng G.L., Yao W., Gao H., Wang K. (2017) A composite dualporosity fractal model for channelfractured horizontal wells, Eng. Appl. Comp. Fluid Mech. 12, 1, 104–116. [Google Scholar]
 Wang Y., Tao Z.W., Chen L., Ma X. (2017) The nonlinear oil–water twophase flow behavior for a horizontal well in triple media carbonate reservoir, Acta Geophys. 65, 5, 977–989. [CrossRef] [Google Scholar]
 Wang Y., Yi X. (2018) Flow modeling of well test analysis for a multiplefractured horizontal well in triple media carbonate reservoir, Int. J. Nonlin. Sci. Num. Simul. 19, 5, 439–457. [CrossRef] [Google Scholar]
 Warren J.E., Root P.J. (1963) The behavior of naturally fractured reservoirs, SPE J. 3, 3, 245–255. [Google Scholar]
 Wu J.Z., Cheng L.S., Li C.L., Cao R.T., Chen C.C., Cao M., Xu Z.Y. (2017) Experimental study of nonlinear flow in micropores under low pressure gradient, Transp. Porous Media 119, 1, 247–265. [CrossRef] [Google Scholar]
 Wu Y.S., Christine E.E., Guan Q., Kang Z.J., Zhang W.M., Tao Q.F. (2007) A triplecontinuum pressuretransient model for a naturally fractured vuggy reservoir, SPE Annual Technical Conference and Exhibition, 11–14 November, Anaheim, California, Society of Petroleum Engineers. [Google Scholar]
 Wu Z.W., Cui C.Z., Lv G.Z., Bing S.X., Gang Cao (2019) A multilinear transient pressure model for multistage fractured horizontal well in tight oil reservoirs with considering threshold pressure gradient and stress sensitivity, J. Petrol. Sci. Eng. 172, 839–854. [CrossRef] [Google Scholar]
 Xu P., Cheng Y.F., Liu X.Y., Zhang X.C., Shi L.B. (2013) Explosive fracturing simulation experiment for low permeability reservoirs and fractal characteristics of cracks produced, Petrol. Explor. Dev. 40, 5, 682–686. [CrossRef] [Google Scholar]
 Xu P., Liu H.C., Sasmito A.P., Qiu S.X., Li C.H. (2017) Effective permeability of fractured porous media with fractal dualporosity model, Fractals 25, 4. [Google Scholar]
 Xu J.C., Sun B.J., Chen B.L. (2019) A hybrid embedded discrete fracture model for simulating tight porous media with complex fracture systems, J. Petrol. Sci. Eng. 174, 131–143. [CrossRef] [Google Scholar]
 Xu P., Yu B.M., Qiu S.X., Cai J.C. (2008) An analysis of the radial flow in the heterogeneous porous media based on fractal and constructal tree networks, Physica A 387, 26, 6471–6483. [CrossRef] [Google Scholar]
 Xu S.L., Yue X.A. (2007) Experimental research on nonlinear flow characteristics at low velocity, J. China Univ. Pet. 31, 5, 60–63. [Google Scholar]
 Yao J., Dai W.H., Wang Z.S. (2004) Well test interpretation method for triple medium reservoir with variable wellbore storage, J. China Univ. Petrol. 28, 1, 46–51. [Google Scholar]
 Yu H., Pan K., Li S., Guo H., He Y., Xu Y., Zhang T., Du S., Cheng S. (2015) Well testing interpretation method and application in triplelayer reservoirs by polymer flooding, Materwiss. Werksttech. 46, 11, 1133–1141. [CrossRef] [Google Scholar]
 Zeng B.Q., Cheng L.S., Li C.L. (2011) Low velocity nonlinear flow in ultralow permeability reservoir, J. Petrol. Sci. Eng. 80, 1, 1–6. [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Seepage physical model of the horizontal well in a triplemedium fractal reservoir. 

In the text 
Fig. 2 Comparison of well test curves of different models. 

In the text 
Fig. 3 Influence of the nonlinear parameter on well test curves. 

In the text 
Fig. 4 Influence of the fractal index on well test curves. 

In the text 
Fig. 5 Influence of the permeability modulus on well test curves. 

In the text 
Fig. 6 Influence of the horizontal well location on well test curves. 

In the text 
Fig. 7 Influence of the horizontal well length on well test curves. 

In the text 
Fig. 8 Well test fitting flow chart. 

In the text 
Fig. 9 Well test model fitting. 

In the text 