© R. Jiang et al., published by IFP Energies nouvelles, 2019

Nomenclature

D_f, θ: Fractal dimension, anomalous diffusion coefficient

C_mt, C_ft, C_vt : Comprehensive compression coefficient of matrix, fracture and cave system, MPa⁻¹

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⁻¹

p_D : Dimensionless bottom hole pressure without considering wellbore reservoir coefficient and skin factor

${\bar{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

μ: Fluid viscosity, mPa s

C_f, α_K : Porosity compressibility coefficient and permeability modulus of fracture system, MPa⁻¹

r₀, L₀ : Capillary radius and length, m

n₀ : Number of capillaries per unit cross-sectional area

δ: Boundary layer thickness, μm

τ₀ : Yield stress value, MPa

l_c : Feature length, m

h, L: Reservoir thickness and half length of the horizontal well, m

z_w : Horizontal well depth, m

c₁ : Nonlinear parameter under the influence of the fluid yield stress and boundary layer

c₂ : 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

β: Fractal index

K_e : Nonlinear coefficient matrix

F_e : Unit load vector

1 Introduction

The fractured-vuggy 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 triple-medium is applied to characterize the complex pores of such reservoirs (Ge and Wu, 1982; Wang Y. et al., 2017). Studying the seepage model of fractured-vuggy 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 fractured-vuggy 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 triple-medium reservoirs considering the variable well storage, innovating the well test theory of triple-medium reservoirs. Camacho-Velazquez et al. (2002) proposed a triple-medium model characterizing a fracture-vuggy 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 fractured-vuggy reservoirs, where caves are embedded in a porous rock and are connected via fracture networks on multiple scales by using the Stokes-Brinkman equations on the fine scale. Nie et al. (2011) established a new flow model in a triple media carbonate reservoir which considers the unsteady inter-porosity 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 fracture-vuggy 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 fracture-vuggy 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. Single-phase 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 multiple-medium 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 dual-porosity model to study the single-phase 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 fractured-vuggy 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 low-speed and non-Darcy 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 low-speed nonlinear seepage phenomenon, the commonly adopted models are mainly divided into three categories: quasi start-up pressure gradient model, section model and continuous model. Quasi start-up pressure gradient model ignores the bending section of the seepage curve, so the model narrows the scope of fluid flow in the low-permeability 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 fractured-vuggy 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 triple-medium reservoir is established. In the conceptual mathematical model, fractured vuggy rock is considered as a triple-continuum 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 fractured-vuggy 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):

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

2. 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.

3. 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 low-permeability reservoirs, small pore throats result in a greater yield stress, and non-Newtonian 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 micro-scale flow effect, resulting in nonlinear seepage strengthening (Cao et al., 2016; Cossio et al., 2013).

This paper combines the micro-seepage mechanism of low-permeability reservoirs, the boundary layer flow and capillary seepage model to derive a new low-speed nonlinear seepage model. It is worth noting that the values of the boundary layer and the starting pressure gradient in low-permeability 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 Hagen-Poiseuille 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. τ₀ is the fluid yield stress, reflecting the influence of fluid non-Newtonian characteristics on seepage flow.

The reservoir core is equivalent to a set of parallel capillary tubes with radius r₀ embedded in solids. If there are n₀ such capillaries per unit of cross-sectional 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}{3r_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₀ = a₁/∇p. If the yield stress value of the same fluid remains unchanged, then 8τ₀/3r₀ can be regarded as a constant value, namely 8τ₀/3r₀ = a₂ (the values of a₁ and a₂ 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{4a_1+a_2}{\nabla p}+\frac{6{a_1}^2+3a_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}{3r_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{4a_1+a_2}{\nabla p}-\frac{6{a_1}^2+3a_1{a}_2}{\nabla p\left(\nabla p-a_1\right)}\right]\nabla p.$$(4)

c₁, c₂, and c₃ are applied to denote the different polynomials of a₁, a₂ $c_1=4a_1+a_2$, c₂ = a₁, $c_3=6{a_1}^2+3a_1{a}_2$). In other words, the value of c₁, c₂, and c₃ can only be obtained by experimental fitting. In order to simplify the formula, c₃ can be represented by c₁ and c₂ through exploring the relationship between c₁, c₂, and c₃. 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₃ can be represented by c₁ and c₂ (c₃ = c₁ c₂). 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₁ reflects the yield stress of the fluid and the influence of the boundary layer while parameter c₂ mainly reflects the influence of boundary layer on the seepage flow. It means c₁/(∇p – c₂) is equivalent to the starting pressure gradient term. When c₁ = 0, the above model is simplified as the Darcy model. But if c₂ = 0, the above model is simplified as the quasi-starting pressure gradient model. When v ≥ 0, then ∇p ≥ c₁ + c₂, so the real minimum starting pressure gradient is ∇p_min = c₁ + c₂.

When c₁ + c₂ = 0, the ∇p_min = 0, if c₁ + c₂ > 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 quasi-linear flow occurs.

2.2 Fractal characteristics of the fracture system

Fractal theory thinks that if fractures have fractal characteristics in fractured-vuggy 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}\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 stress-sensitive 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 }\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 fractured-vuggy 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 fractured-vuggy tight fractal reservoir

3.1 Physical model

A physical model of the horizontal well in a circular triple-medium fractal reservoir just as shown in Figure 1 was established, which met the following assumptions:

1. The model is a three-hole, single-permeability circular reservoir consisting of the matrix system, fracture system and cave system. The quasi-steady-state cross flow generates from the matrix system and cave system to the fracture system.

2. 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.

3. 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.

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

5. The influence of gravity and the capillary force was ignored.

6. The stress sensitivity, fractal characteristics of the fracture system and the nonlinear seepage mechanism of the matrix were considered.

What we study in this paper is a carbonated salt reservoir which is well-developed 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 single-phase transient flow and fractured vuggy rock, conceptualized also as a triple-or multiple-continuum medium, consisting of highly permeable fractures, low-permeability rock matrix, and vugs.

3.2 Mathematical model

In a Cartesian coordinate system, the mass conservation equation of single-phase 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 fractured-vuggy 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³/m³) 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 10^{-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 10^{-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 10^{-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 source-sink 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 source-sink 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}}}^{n-1}+\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\iiint }{\lambda }_{\mathrm{mf}}\mathbf{DN}{\mathbf{N}}^{\mathbf{T}}\mathrm{d}V{{\mathit{P}}_{\mathrm{m}}}^{n-1}+\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\iiint }{\lambda }_{\mathrm{vf}}\mathbf{EN}{\mathbf{N}}^{\mathbf{T}}\mathrm{d}V{{\mathit{P}}_{\mathrm{v}}}^{n-1}.$$(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}}}^{n-1}+\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}}}^{n-1}+\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\iiint }{\lambda }_{\mathrm{vf}}\mathbf{EN}{\mathbf{N}}^{\mathbf{T}}\mathrm{d}V{{\mathit{P}}_{\mathrm{v}}}^{n-1}.$$(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}}^{n-1}+\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}}}^{n-1}+\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.

$${\bar{p}}_{\mathrm{wD}}=\frac{s{\bar{p}}_{\mathrm{D}}+\mathrm{S}{L}_{\mathrm{D}}}{s\left[1+C_{\mathrm{D}}s\left(s{\bar{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{\ln 2}{t_{\mathrm{D}}}\sum _{i=1}^N{V}_i{\bar{p}}_{\mathrm{wD}}.$$(36)

In equation (33): $V_i=(-1{)}^{\left(\frac{N}{2}+i\right)}\sum _{k=\frac{i+1}{2}}^{\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(i-k+1\right)!\left(2k-i+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 triple-medium low-permeability 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.

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).

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 V-shape 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.

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.

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.

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.

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 triple-medium 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.

7 Conclusion

1. A nonlinear seepage model of a horizontal well in a triple-medium 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.

2. 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.

3. The nine stages were summarized in the seepage process namely pure wellbore storage stage, skin effect transition stage, early radial flow, early linear flow, mid-term 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.

4. 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.

5. The triple-medium 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 triple-medium 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

All Tables

All Figures