© Y.-M. Yao et al., published by IFP Energies nouvelles, 2019

Nomenclature

B: Capillary pressure constant, MPa

h: Thickness of the reservoir, m

K: Absolute permeability, mD, 10⁻¹⁵m²

K_r : Relative permeability

K_rg,eq1 : Equivalent relative permeability

L_eq : Equivalent interblock distance, m

m: Pseudo-pressure, kg/(m³ · s)

P: Pressure, Pa

P_c : Capillary pressure, Pa

P_init : Initial pressure of gas reservoir, Pa

P^(e) : Pressure at equivalent radius

q: Phase flow rate, m³/d

r_e : Equivalent radius, m

r_w : Wellbore radius, m

S: Water phase normalized saturation

S_o : Outlet saturation in the condition of only water flowing into wellbore

S_g : Gas phase saturation

S_rw : Water saturation at the entrance of wellbore

T: Temperature, K

T_eq : Equivalent transmissibility, m⁴ · s kg⁻¹

Δx: Grid length, m

ϕ: Porosity

λ: Mobility, m · s kg⁻¹

λ_w,eq1 : Equivalent mobility, m · s kg⁻¹

μ: Viscosity, Pa · s

ρ: Density, kg/m³

Subscripts and superscripts

g: Gas phase

w: Water phase

(i, j): Physical quantity in well-contained grid

(adj): Physical quantity in adjacent grid

inlet: Physical quantity at the inlet

outlet: Physical quantity at the outlet

+: Left-hand limit

−: Right-hand limit

1 Introduction

Numerical reservoir simulation for oil or gas production must account for the presence of wells. It is necessary to use grid blocks whose horizontal dimensions are much larger than the diameter of a well (Williamson and Chappelear, 1981). As a result, the pressure calculated for a well block is greatly different from the bottom-hole pressure of the well. In the reservoir simulation, it is necessary to build a numerical well model to calculate the well flow rate associating with the well block pressure and the bottom-hole pressure. Peaceman well model, constructed from the analytical solution of One-Dimensional (1D) radially symmetric single-phase flow, is widely used in the general reservoir simulator. Initially, Peaceman’s finite difference well model was for the well located at the center of square grid for single phase flow (Peaceman, 1978). Then, it was extended into various applications, including horizontal wells, anisotropic reservoirs, and also multiphase flow, etc. (Babu and Odeh, 1988; Babu et al., 1991; Peaceman, 1983; Shiralkar, 1989; Su, 1995). Chen and Zhang (2009) extended the well model for other numerical methods such as standard finite element, control volume finite element, and mixed finite element methods.

By analogy to single-phase flow, numerical well models are directly extended to multi-phase flow without the solid support of mathematical analysis. The single-phase flow in the neighborhood of the well is described by the linear equations, in contrast to that the multi-phase flow should be described by the nonlinear equations. Only when the phase saturations remain constant, the direct extension to multi-phase flow is accurate. It has been shown that, at the interface of two different type rocks which have the different relative permeability and capillary pressure curves, multiphase flows in porous media face the problem of the discontinuities of fluid properties and transport properties (Ding et al., 2016; Wang et al., 2006). This implies that the physical quantities might vary quickly in the neighborhood of the production wells.

It has been observed that in the core experiment, when two-phase fluid flows into the non-capillary force container from the core sample containing capillary force, the wetting phase saturations at the outlet end of the core increase significantly compared with the rest of the area (Richardson et al., 1952). This phenomenon was called as the capillary end effect. Ignoring the capillary end effect, the measurements of relative permeability will produce big errors. Assuming that water-phase saturation is 1 at the outlet end of the water-wet core, a saturation distribution equation near the outlet end was established (Hadley and Handy, 1956). Then the way to correct the errors in the relative permeability measurements, considering the capillary end effect, has been proposed (Huang and Honarpour, 1998; Virnovsky et al., 1995).

In the low-permeability reservoir, the capillary pressure has an important role on multi-phase flows. The capillary pressure can be ignored inside the wellbore while it should be considered for the multi-phase flow inside the reservoir. Analogous to the core experiment described above, the capillary end effect may also have a significant impact on fluid flows in the vicinity of production well. Actually, it has already been shown that the capillary end effect in water-wet reservoir creates a water layer near the wellbore with residual saturation of gas or oil. This layer will cause the hydraulic resistance for gas or oil flow (Barenblatt et al., 1989; Bedrikovetsky, 2013). Holditch (1979) found out that water phase will be completely trapped in reservoir if the pressure drop between the formation pressure and the bottom-hole pressure cannot overcome the capillary end effect. Even when the pressure drop is large enough to overcome the capillary end effect, severe reduction in gas production may also occur, depending on the level of capillary pressure inside the reservoir. Settari and Aziz (1974) analyzed the incompressible two-phase flow in the vicinity of the production well. Ding (1995) presented a scaling-up procedure for the radial flow near the well which consists in scaling up transmissibility and numerical PI in the vicinity of the well, then he extended this technique to the near-well upscaling and formation damage modeling which can give quite satisfactory results for long-term productivity prediction (Ding, 2004, 2009, 2011). Naik et al. (2018) established a semi-analytical modeling of steady-state gas-water flow in the rock with constant and piecewise-constant contact angles accounting for gas compressibility and found out that, water saturation significantly rises at the end-effect zone, creating a hydraulic resistance for gas flow. The capillary end effect may have the significant impact on the production well. As far as we know, this capillary end effect has not been considered in the existing well model for reservoir simulation. In this article, we will construct an alternative numerical well model to describe the gas recovery for the low-permeability water-wet reservoir.

For the convenience of generalization to non-uniform grids and flexible grids, Ding proposed a more accurate and effective numerical well model for the single-phase flow (Ding and Jeannin, 2004; Ding and Renard, 1994; Ding et al., 1998). The technique proposed by Ding et al. is applied. The article is organized as follows. The numerical well models are introduced in Section 2. In Section 3, the 1D radial steady flow for immiscible two-phase fluid toward well is investigated. Section 4 describes the procedure to establish the well model for multi-phase flow. In Section 5, the accuracy of the proposed well model is verified by numerical examples and the capillary end effect on the gas and water production is analyzed. Conclusion and discussion are given in Section 6.

2 Brief review of the previous numerical well models

In the neighborhood of the production well, the single-phase flow can be approximated by 1D radially symmetric steady flow:

$$\frac{\mathrm{d}}{\mathrm{d}r}\left(r\frac{\mathrm{d}P}{\mathrm{d}r}\right)=0.$$(1)Assuming the isolated well is at (x₀, y₀), the pressure distribution is given by:

$$P\left(x,y\right)=P_{\mathrm{wf}}+\frac{\mu q}{2\pi \rho \mathrm{Kh}}\ln (r/r_{\mathrm{w}}),$$(2)where $r=\sqrt{(x-x_0{)}^2+(y-y_0{)}^2}$.

Consider the single phase flow in two dimensional grids with Δx = Δy (see Fig. 1). For the incompressible single-phase flow, the exchange flow rate between the control volume (i, j) and its neighbors can be obtained from the discretization of Darcy’s equation as follows:

$$q=\frac{\rho \mathrm{Kh}}{\mu }\left[P^{\left(i-1,j\right)}+P^{\left(i+1,j\right)}+P^{\left(i,j-1\right)}+P^{\left(i,j+1\right)}-4P^{\left(i,j\right)}\right],$$(3)where the variables ρ, μ and K represent the fluid density, fluid viscosity and permeability, respectively. If supposing P^(i+1,j) = P^(i−1,j) = P^{(i, j+1)} = P^(i,j-1) = P^(adj), the relation equation (3) becomes:

$$q=\frac{4\rho Kh}{\mu }\left[P^{\left(\mathrm{adj}\right)}-P^{\left(i,j\right)}\right].$$(4)

To establish the numerical well model, an equivalent radius r_e is introduced, where the flowing pressure equals P^(i,j). With the boundary condition $P{|}_{r=r_{\mathrm{e}}}=P^{(i,j)}$ and P|_r=Δx = P^(adj), we can integrate equation (1) to get the pressure distribution. The flow rate $q=2\pi r\frac{\rho Kh}{\mu }\frac{\mathrm{d}P}{\mathrm{d}r}$ is then calculated to be

$$q=\frac{2\pi \rho Kh}{\mu }\frac{\left[P^{\left(\mathrm{adj}\right)}-P^{\left(i,j\right)}\right]}{\ln \left(\Delta x/r_{\mathrm{e}}\right)}.$$(5)

Combing equations (4) and (5), we can get

$$r_{\mathrm{e}}=e^{-\pi /2}\Delta x\approx 0.208\Delta x.$$(6)

Integrating equation (1) with the boundary condition $P{|}_{r=r_{\mathrm{w}}}=P_{\mathrm{wf}}$ and ${P|}_{r=r_{\mathrm{e}}}=P^{(i,j)}$ again, we have

$$q=\frac{2\pi \rho Kh}{\mu }\frac{\left[P^{\left(i,j\right)}-P_{\mathrm{wf}}\right]}{\ln \left(r_{\mathrm{e}}/r_{\mathrm{w}}\right)},$$(7)where P_wf and r_w are respectively the bottom-hole pressure and the radius of the production well. The relation equation (7) with the definition of equivalent radius equation (6) is just the well-known Peaceman well model. The above derivations can be found out in the literature (Peaceman, 1978).

Without the solid mathematical analysis, Peaceman well model equation (7) for single-phase flow is directly extended to multi-phase flow and then get

$$q_{\alpha }=\frac{2\pi {\rho }_{\alpha }K{K}_{\mathrm{r\alpha }}h}{{\mu }_{\alpha }}\frac{P_{\alpha }^{\left(i,j\right)}-P_{\mathrm{wf}}}{\ln \left(r_{\mathrm{e}}/r_{\mathrm{w}}\right)},$$(8)where the subscript α represents the different phase, K_rα is the relative permeability and h is the layer thickness. Here, r_e is still calculated by the relation equation (6). The Peaceman well models equations (6) and (7) are widely used in reservoir simulators.

As pressure varies rapidly logarithmically near the well, its gradient cannot be approximated by linear relation. To improve the accuracy and the efficiency for the prediction, Ding and Renard (1994) proposed an alternative numerical well model based upon the analytical solution equation (2). The finite-volume method applied to the flow equation for the wellblock yields:

$$q=\frac{Kh}{\mu }\left({\int }_{\mathrm{\Gamma }}\frac{\partial P}{\partial x}\mathrm{d}y+{\int }_{\mathrm{\Gamma }}\frac{\partial P}{\partial y}\mathrm{d}x+{\int }_{\mathrm{\Gamma }}\frac{\partial P}{\partial x}\mathrm{d}y+{\int }_{\mathrm{\Gamma }}\frac{\partial P}{\partial y}\mathrm{d}x\right).$$(9)

Take the interface Γ1 between wellblock and adjacent block (i + 1, j) as example, the flow-term at the interface is

$$q_{\mathrm{\Gamma }1}={\int }_{\mathrm{\Gamma }}\frac{\partial P}{\partial x}\mathrm{d}y=\frac{\mu q}{\pi Kh}\arctan \frac{\Delta y}{\Delta x}.$$(10)

Ding et al. (1998) defined an equivalent interblock distance L_eq1 to replace the original interblock distance Δx_1/2. The pressure P^(e) is located at the equivalent radius r_e. The numerical flow-term at the interface can be expressed as

$$q_{\mathrm{\Gamma }1}=T_{\mathrm{eq},1}\left(P^{\left(i+1,j\right)}-P^{\left(e\right)}\right)=\frac{Kh}{\mu }\Delta y\frac{P^{\left(i+1,j\right)}-P^{\left(e\right)}}{L_{\mathrm{eq}1}},$$(11)where T_eq,1 is the equivalent transmissibility to relate P^(e) and adjacent grid pressure through equivalent interblock distance. From equations (2), (10) and (11), the equivalent interblock distance and equivalent transmissibilities can be respectively given as:

$$L_{\mathrm{eq}1}=\frac{\Delta y}{2}\frac{\ln \frac{\Delta x_{1/2}}{r_{\mathrm{e}}}}{\arctan \frac{\Delta y}{\Delta x}} \hspace{1em}\mathrm{and}\hspace{1em}{T}_{\mathrm{eq},1}=\frac{Kh}{\mu }\frac{\Delta y}{L_{\mathrm{eq}1}}.$$(12)

For the other three directions, equivalent interblock distances and equivalent transmissibilities can be calculated using the same way.

Substituting the flow-terms at four directions into equation (9), the pressure P^(e) located at equivalent radius r_e can be given by:

$$q=\left[T_{\mathrm{eq},1}{P}^{\left(i+1,j\right)}+T_{\mathrm{eq},2}{P}^{\left(i,j+1\right)}+T_{\mathrm{eq},3}{P}^{\left(i-1,j\right)}+T_{\mathrm{eq},4}{P}^{\left(i,j-1\right)}-\left(T_{\mathrm{eq},1}+T_{\mathrm{eq},2}+T_{\mathrm{eq},3}+T_{\mathrm{eq},4}\right)P^{\left(e\right)}\right].$$(13)

Substituting the P^(e) into the analytical solution equation (2), the equation relating pressure, wellbore pressures and the flow rate of the well can be expressed as:

$$q=\frac{2\pi \rho Kh}{\mu }\frac{\left[P^{\left(e\right)}-P_{\mathrm{wf}}\right]}{\ln \left(r_{\mathrm{e}}/r_{\mathrm{w}}\right)}.$$(14)

Equations (13) and (14) constitute the new numerical well model for the incompressible single-phase flow. It should be noted that the equivalent radius r_e can be any value in this model.

3 Two flow patterns for the one-dimensional radial steady flow for immiscible two-phase fluid toward well

In this section, the 1D radial steady flow for immiscible two phase fluid (gas and water) toward well is investigated. A condition of constant temperature is supposed in this paper. And the governing equations can be written as:

$$\{\begin{array}{l}\frac{\mathrm{d}{q}_{\mathrm{w}}}{\mathrm{d}r}=\frac{\mathrm{d}}{\mathrm{d}r}\left(-2\pi r{\rho }_{\mathrm{w}}\frac{K{K}_{r\mathrm{w}}}{{\mu }_{\mathrm{w}}}\frac{\mathrm{d}{P}_{\mathrm{w}}}{\mathrm{d}r}\right)=0\\ \frac{\mathrm{d}{q}_{\mathrm{g}}}{\mathrm{d}r}=\frac{\mathrm{d}}{\mathrm{d}r}\left(-2\pi r{\rho }_{\mathrm{g}}\frac{K{K}_{r\mathrm{g}}}{{\mu }_{\mathrm{g}}}\frac{\mathrm{d}{P}_{\mathrm{g}}}{\mathrm{d}r}\right)=0\end{array}.$$(15)

Only the gas compressibility is considered here, and it varies with the change of gas phase pressure:

$${\rho }_{\mathrm{g}}={\rho }_{\mathrm{g}}\left(P_{\mathrm{g}}\right).$$(16)

Denote S_wi as the irreducible water saturation and S_gr as the residual saturation of gas phase, then the water saturation S_w can be normalized as:

$$S=\frac{S_{\mathrm{w}}-S_{\mathrm{wi}}}{1-S_{\mathrm{wi}}-S_{\mathrm{gr}}}.$$(17)

The relative permeability K_rα of the α-phase can be expressed as the function of the normalized water saturation. That is:

$$K_{\mathrm{rw}}=K_{\mathrm{rw}}\left(S\right),\hspace{1em}{K}_{\mathrm{rg}}=K_{\mathrm{rg}}\left(S\right),$$(18)with K_rw being a strictly increasing function, but K_rg being a decreasing function. The capillary pressure, defined as the difference between the gas and water pressures, is also a decreasing function of the normalized water saturation S:

$$P_{\mathrm{c}}\left(S\right)=P_{\mathrm{g}}-P_{\mathrm{w}}.$$(19)

Capillary pressure is relatively large inside the reservoir but can be ignored inside the wellbore. There are two different capillary pressure curves on the two sides of the entrance of the production well. Analogous to fluid flows through the interface of the different type rocks, physical quantities may change significantly around the interface. Settari and Aziz (1974) pointed out that both water and gas can flow into the production well when the normalized water saturation reaches the value of 1 on the wall of the wellbore and otherwise water phase will be trapped in the reservoir. That means there are two different flow patterns at the vicinity of the wellbore.

Define P_inlet and S_inlet as the gas-phase pressure and normalized water-phase saturation of the point r = r_inlet in the vicinity of the production well. When P_inlet ∈ (P_wf + P_c(S_inlet), ∞), both two phase can flow into wellbore, i.e., q_g > 0 and q_w > 0 . When P_inlet ∈ (P_wf, P_wf + P_c(S_inlet)], only the gas can flow into wellbore, i.e., q_g > 0 and q_w = 0. Details of these two flow patterns are provided in the following sections.

3.1 Flow pattern 1: both two phases flowing into wellbore

In the case where both water and gas flow into the well, we will have:

$${S|}_{r=r_{\mathrm{w}}^{+}}=1.$$(20)

The reason that the relation equation (20) holds is as follow. If ${S|}_{r=r_{\mathrm{w}}^{+}}<1$, then ${K_{\mathrm{rg}}|}_{r=r_{\mathrm{w}}^{+}}>0$ and ${P_{\mathrm{c}}|}_{r=r_{\mathrm{w}}^{+}}>0$. Thus, ${P_{\mathrm{g}}|}_{r=r_{\mathrm{w}}^{+}}>{P_{\mathrm{w}}|}_{r=r_{\mathrm{w}}^{+}}$. Since there is no capillary pressure in the well, we have ${P_{\mathrm{w}}|}_{r=r_{\mathrm{w}}^{-}}={P_{\mathrm{g}}|}_{r=r_{\mathrm{w}}^{-}}=P_{\mathrm{wf}}$. The water pressure must be continuous at r = r_w, so ${P_{\mathrm{g}}|}_{r=r_{\mathrm{w}}^{+}}>{P_{\mathrm{w}}|}_{r=r_{\mathrm{w}}^{+}}={P_{\mathrm{w}}|}_{r=r_{\mathrm{w}}^{-}}={P_{\mathrm{g}}|}_{r=r_{\mathrm{w}}^{-}}$. Combined with ${K_{\mathrm{rg}}|}_{r=r_{\mathrm{w}}^{+}}>0$, the relation ${P_{\mathrm{g}}|}_{r=r_{\mathrm{w}}^{+}}>{P_{\mathrm{g}}|}_{r=r_{\mathrm{w}}^{-}}$ will result in infinite velocity of gas phase, which is impossible. Thus, the relation equation (20) holds. Under the condition ${S|}_{r=r_{\mathrm{w}}^{+}}=1$, we have ${K_{\mathrm{rg}}|}_{r=r_{\mathrm{w}}^{+}}=0$. Considering the gas flow rate q_g > 0, so ${\frac{\mathrm{d}{P}_{\mathrm{g}}}{\mathrm{d}r}|}_{r=r_{\mathrm{w}}^{+}}\to +\mathrm{\infty }$. This will lead to the relation ${\frac{\mathrm{d}S}{\mathrm{d}r}|}_{r=r_{\mathrm{w}}^{+}}\to -\mathrm{\infty }$ (see Fig. 2).

3.2 Flow pattern 2: only gas phase flowing into wellbore

If P_inlet ∈ (P_wf, P_wf + P_c(S_inlet)], then we have ${P_{\mathrm{w}}|}_{r=r_{\mathrm{inlet}}}\le P_{\mathrm{wf}}$ and ${P_{\mathrm{g}}|}_{r=r_{\mathrm{inlet}}}>P_{\mathrm{wf}}$. Thus, the water phase is trapped in the reservoir and only the gas phase can flow into the wellbore. In this article, we only consider the production well case, so we have:

$$q_{\mathrm{w}}=0,\hspace{1em}\frac{\mathrm{d}{P}_{\mathrm{w}}}{\mathrm{d}r}=0.$$(21)Then, we obtain:

$$P_{\mathrm{inlet}}-P_{\mathrm{c}}\left(S_{\mathrm{inlet}}\right)=P_{\mathrm{wf}}-P_{\mathrm{c}}\left(S_{\mathrm{rw}}\right)$$(22)where S_rw represents saturation at the entrance of wellbore. The saturation and pressure distributions are shown in Figure 3.

Given the pressure P_inlet and saturation S_inlet, the pressure and saturation in the vicinity of the production well for the two flow patterns can be determined, which will be given in the next section.

4 Modified well model

In Section 2, numerical well models for the incompressible single-phase flow are reviewed. In the following, we will use Ding’s finite-volume method to derive an alternative numerical well model for multi-phase flow based upon the analytical equation of the two flow patterns described in the above section.

The finite-volume method applied to the flow equation for the wellblock yields:

$$\{\begin{array}{l}{q}_{\mathrm{w}}=\frac{Kh}{\mu }\left({\int }_{\mathrm{\Gamma }}{\lambda }_{\mathrm{w}}\frac{\partial P}{\partial x}\mathrm{d}y+{\int }_{\mathrm{\Gamma }}{\lambda }_{\mathrm{w}}\frac{\partial P}{\partial y}\mathrm{d}x+{\int }_{\mathrm{\Gamma }}{\lambda }_{\mathrm{w}}\frac{\partial P}{\partial x}\mathrm{d}y+{\int }_{\mathrm{\Gamma }}{\lambda }_{\mathrm{w}}\frac{\partial P}{\partial y}\mathrm{d}x\right)\\ q_{\mathrm{g}}=\frac{Kh}{\mu }\left({\int }_{\mathrm{\Gamma }}{K}_{\mathrm{rg}}\frac{\partial m}{\partial x}\mathrm{d}y+{\int }_{\mathrm{\Gamma }}{K}_{\mathrm{rg}}\frac{\partial m}{\partial y}\mathrm{d}x+{\int }_{\mathrm{\Gamma }}{K}_{\mathrm{rg}}\frac{\partial m}{\partial x}\mathrm{d}y+{\int }_{\mathrm{\Gamma }}{K}_{\mathrm{rg}}\frac{\partial m}{\partial y}\mathrm{d}x\right)\end{array},$$(23)where λ_α = K_rα/μ_α is the mobility of the α-phase fluid and m is the pseudo-pressure defined by:

$$m=\int \frac{{\rho }_{\mathrm{g}}\left(P_{\mathrm{g}}\right)}{{\mu }_{\mathrm{g}}}\mathrm{d}{P}_{\mathrm{g}}.$$(24)

Take the interface Γ1 between wellblock and adjacent grid (i + 1, j) as example. Considering the axisymmetric properties of the phase pressure distributions, from equation (15), the flow-term between the adjacent grid and wellblock grid can be obtained

$$\{\begin{array}{l}{q}_{\mathrm{w}\mathrm{\Gamma }1}={\int }_{\mathrm{\Gamma }}Kh{\lambda }_{\mathrm{w}}\frac{\partial P_{\mathrm{w}}}{\partial x}\mathrm{d}y=q_{\mathrm{w}}\frac{2\arctan \frac{\Delta y}{\Delta x}}{2\pi }\\ q_{\mathrm{g}\mathrm{\Gamma }1}={\int }_{\mathrm{\Gamma }}Kh{K}_{\mathrm{rg}}\frac{\partial m}{\partial x}\mathrm{d}y=q_{\mathrm{g}}\frac{2\arctan \frac{\Delta y}{\Delta x}}{2\pi }\end{array}.$$(25)

To let the equivalent interblock distance be phase-independent, we choose its value L_eq1 = [Δyln(Δx_1/2/r_e)]/[2arctan(Δy/Δx)], just the same as that in the single-phase flow given by equation (12). Considering that the different phases have the different pressure distributions, we put the flux corrections into the mobility and relative permeability terms, which is expressed as following:

$$\{\begin{array}{l}{q}_{\mathrm{w}\mathrm{\Gamma }1}=T_{\mathrm{eqw},1}(P_{\mathrm{w}}^{(i+1,j)}-P_{\mathrm{w}}^{\left(e\right)})=Kh{\rho }_{\mathrm{w}}{\lambda }_{\mathrm{w},\mathrm{eq}1}\Delta y\frac{[P_{\mathrm{w}}^{(i+1,j)}-P_{\mathrm{w}}^{\left(e\right)}]}{L_{\mathrm{eq}1}}\\ q_{\mathrm{g}\mathrm{\Gamma }1}=T_{\mathrm{eqg},1}(m^{(i+1,j)}-m^{\left(e\right)})=Kh{K}_{\mathrm{rg},\mathrm{eq}1}\Delta y\frac{[m^{\left(i+1,j\right)}-m^{\left(e\right)}]}{L_{\mathrm{eq}1}}\end{array},$$(26)where λ_w,eq1 and K_rg,eq1 are the equivalent mobility and equivalent relative permeability, respectively.

With the help of the boundary condition ${P_{\alpha }|}_{r=r_{\mathrm{e}}}=P_{\alpha }^{\left(e\right)}$ and ${P_{\alpha }|}_{r=\Delta x}=P_{\alpha }^{(i+1,j)}$, we integrate equation (15) from r = r_e to r = Δx_1/2 and then obtain

$$\{\begin{array}{l}{q}_{\mathrm{w}}=\frac{\left[P_{\mathrm{w}}^{\left(i+1,j\right)}-P_{\mathrm{w}}^{\left(e\right)}\right]h}{{\int }_{\Delta x_{1/2}}^{r_e}\frac{1}{2\pi \mathrm{rK}{\rho }_{\mathrm{w}}{\lambda }_{\mathrm{w}}}\mathrm{d}r}\\ q_{\mathrm{g}}=\frac{\left[m^{\left(i+1,j\right)}-m^{\left(e\right)}\right]h}{{\int }_{\Delta x_{1/2}}^{r_e}\frac{1}{2\pi \mathrm{rK}{K}_{r\mathrm{g}}}\mathrm{d}r}\end{array}.$$(27)

From equations (25)–(27), the equivalent mobility (λ_w,eq1), equivalent relative permeability (K_rg,eq1) and equivalent transmissibilities for the water and gas phases can be obtained:

$$\{\begin{array}{l}{\lambda }_{\mathrm{w},\mathrm{eq}1}=\frac{2L_{\mathrm{eq}1}\arctan \frac{\Delta y}{\Delta x}}{\Delta y{\int }_{r_{\mathrm{e}}}^{\Delta x_{1/2}}\frac{1}{r{\lambda }_{\mathrm{w}}}\mathrm{d}r}\\ K_{\mathrm{rg},\mathrm{eq}1}=\frac{2L_{\mathrm{eq}1}\arctan \frac{\Delta y}{\Delta x}}{\Delta y{\int }_{r_{\mathrm{e}}}^{\Delta x_{1/2}}\frac{1}{r{K}_{\mathrm{rg}}}\mathrm{d}r}\end{array}\hspace{1em}\mathrm{and}\hspace{1em}\{\begin{array}{l}{T}_{\mathrm{eq}1,\mathrm{w}}=\frac{Kh{\rho }_{\mathrm{w}}{\lambda }_{\mathrm{w},\mathrm{eq}1}\Delta y}{L_{\mathrm{eq}1}}\\ T_{\mathrm{eq}1,\mathrm{g}}=\frac{Kh{K}_{\mathrm{rg},\mathrm{eq},1}\Delta y}{L_{\mathrm{eq}1}}\end{array}.$$(28)

The other three adjacent grid equivalent mobilities (or relative permeability) and equivalent transmissibilities can be calculated as the same method. Substituting flow-terms at four directions into equation (23), the pressure P^(e) located at r_e can be determined by:

$$\{\begin{array}{l}{q}_{\mathrm{w}}=\left[T_{\mathrm{eqw},1}{P}_{\mathrm{w}}^{\left(i+1,j\right)}+T_{\mathrm{eqw},2}{P}_{\mathrm{w}}^{\left(i,j+1\right)}+T_{\mathrm{eqw},3}{P}_{\mathrm{w}}^{\left(i-1,j\right)}+T_{\mathrm{eqw},4}{P}_{\mathrm{w}}^{\left(i,j-1\right)}-\left(T_{\mathrm{eqw},1}+T_{\mathrm{eqw},2}+T_{\mathrm{eqw},3}+T_{\mathrm{eqw},4}\right)P^{\left(e\right)}\right]\\ q_{\mathrm{g}}=\left[T_{\mathrm{eqg},1}{m}^{\left(i+1,j\right)}+T_{\mathrm{eqg},2}{m}^{\left(i,j+1\right)}+T_{\mathrm{eqg},3}{m}^{\left(i-1,j\right)}+T_{\mathrm{eqg},4}{m}^{\left(i,j-1\right)}-\left(T_{\mathrm{eqg},1}+T_{\mathrm{eqg},2}+T_{\mathrm{eqg},3}+T_{\mathrm{eqg},4}\right)m^{\left(e\right)}\right]\end{array}.$$(29)

After obtaining the pressure P^(e) located at equivalent radius, integrating equation (15) on the interval [r_w, r_e] again in the wellblock grid, we obtain the calculation formula for two-phase flow expressed as:

$$\{\begin{array}{l}{q}_{\mathrm{w}}=\frac{\left[P_{\mathrm{w}}^{\left(e\right)}-P_{\mathrm{wf}}\right]h}{{\int }_{r_{\mathrm{w}}}^{r_{\mathrm{e}}}\frac{1}{2\pi \mathrm{rK}{\rho }_{\mathrm{w}}{\lambda }_{\mathrm{w}}}\mathrm{d}r}\\ q_{\mathrm{g}}=\frac{\left[m^{\left(e\right)}-m_{\mathrm{wf}}\right]h}{{\int }_{r_{\mathrm{w}}}^{r_{\mathrm{e}}}\frac{1}{2\pi \mathrm{rK}{K}_{\mathrm{rg}}}\mathrm{d}r}\end{array}.$$(30)

Equations (29) and (30) constitute the two-phase flow well model based on finite-volume method. And the equivalent radius can also be any value.

To get the gas and water production from the formula, the saturation distribution on the interval [r_w, r_e] should be calculated. In the following, we present the procedure to calculate the saturation distribution for two different flow patterns.

4.1 Both of two phase flowing into wellbore (when $P_{\mathrm{g}}^{(i,j)}-P_{\mathrm{c}}\left(S^{(i,j)}\right)>P_{\mathrm{wf}}$)

The phase pressures can be given by the following equation:

$$\{\begin{array}{l}\frac{\mathrm{d}}{\mathrm{d}r}\left(2\pi r{\rho }_{\alpha }h\frac{K{K}_{\mathrm{r}\mathrm{\alpha }}}{{\mu }_{\alpha }}\frac{\mathrm{d}{P}_{\alpha }}{\mathrm{d}r}\right)=0\\ P_{\mathrm{c}}\left(S\right)=P_{\mathrm{g}}-P_{\mathrm{w}}\end{array},$$(31)with the boundary condition:

$${P_{\mathrm{\alpha }}|}_{r=r_{\mathrm{w}}}=P_{\mathrm{wf}}\hspace{1em}\mathrm{and}\hspace{1em}{P_{\mathrm{\alpha }}|}_{r=r_{\mathrm{e}}}=P_{\mathrm{\alpha }}^{\left(i,\mathrm{j}\right)}\hspace{1em}\left(\alpha =\mathrm{w},\mathrm{g}\right).$$(32)

4.2 Only gas phase flowing into wellbore (when $P_{\mathrm{g}}^{(i,j)}-P_{\mathrm{c}}\left(S^{(i,j)}\right)\le P_{\mathrm{wf}}$)

Substituting equation (21) into equation (15), we obtain the following saturation equation:

$$\frac{\mathrm{d}}{\mathrm{d}r}\left(-2\pi r{\rho }_{\mathrm{g}}\left(P_{\mathrm{g}}\right)h\frac{K{K}_{\mathrm{rg}}\left(S_{\mathrm{g}}\right)}{{\mu }_{\mathrm{g}}}\frac{\mathrm{d}{P}_{\mathrm{c}}\left(S\right)}{\mathrm{d}S}\frac{\mathrm{d}S}{\mathrm{d}r}\right)=0,$$(33)with the boundary condition:

$${S|}_{r=r_{\mathrm{w}}}=S_{\mathrm{rw}}\hspace{1em}\mathrm{an}\mathrm{d}\hspace{1em}{S|}_{r=r_{\mathrm{e}}}=S^{\left(i,j\right)},$$(34)where the saturation S_rw is determined by equation (22) and the phase pressures are given by equation (35):

$$P_{\mathrm{g}}^{\left(i,j\right)}-P_{\mathrm{c}}\left(S^{\left(i,j\right)}\right)=P_{\mathrm{wf}}-P_{\mathrm{c}}\left(S_{\mathrm{rw}}\right).$$(35)

Equations (31)–(35) can be solved numerically. In order to describe the dramatic changes of saturation and pressure near the wellbore, the logarithmic coordinate grids are adopted. In the following numerical tests, the solution interval [r_w, r_e] is divided into enough equal parts on the logarithmic coordinate. Newton-Raphson iteration method is applied in the solution of the discretized nonlinear equations for equations (31)–(35) to obtain the saturation and pressure distributions for two-phase flow and only gas flow, respectively. Practically, a short subroutine can be added into the numerical simulator. Due to solving the ordinary differential equations, the computation speed is very fast. The additional modification has little influence on the program speed.

5 Numerical results

In this section, numerical tests are performed to validate the accuracy of the proposed model. In the Two-Dimensional (2D) examples of Test 1 and Test 2, the layer thickness is assumed to be unit thickness, i.e., h = 1. For comparison, the numerical results calculated from the original Peaceman well model equation (7) are also provided. Based upon the calculations, the capillary end effects on the predictions of the gas recovery for the low-permeability reservoir are analyzed. In the proposed numerical well model the equivalent radius (r_e) can be any value. Without loss of generality, in the following calculations, we choose r_e = 0.208Δx, just the same as that in the Peaceman well model.

5.1 Test 1

In this numerical test, the simulation area is a 2D square reservoir with the size 500 × 500 m. There are five production wells in the gas reservoir, whose distribution is shown in Figure 4. The outer boundary is imposed with impermeable boundary conditions. Two-phase immiscible flow for water and natural gas is considered here. Water is incompressible with invariable density 10³ kg/m³ and gas is compressible with its density varying with its pressure. The relevant parameters of the reservoir and fluids are shown in Table 1. The capillary pressure curve takes the form P_c = −BlnS, where B = 1 MPa. The capillary pressure curves and relative permeability are shown in Figure 5. In this test, we aim to explore the accuracy of proposed well model for the multi-well problem. Since the grid size is much bigger than wellbore radius, the impermeable boundary has little effect on the well model performance.

To test the accuracy of the original Peaceman and the proposed well models, different grid systems of 50 × 50, 150 × 150 and 250 × 250 are adopted in the computations. Due to that the locations of four production wells in the corner have the geometric symmetry, only the calculations of Well 1 and Well 2 are analyzed. For the calculations of the proposed well model, it is shown that the water and gas productions weakly depend upon the grid size, i.e., the calculated results obtained from the different grid sizes are very close (see Fig. 6). However, if adopting the original Peaceman model, there are the obvious deviations among the results for water and gas productions of the different grid sizes.

The water saturation increasing phenomenon will cause the hydraulic resistance for gas flow and increase the water production. All the grid sizes in this example are much larger than the size of the water layer near the wellbore. Therefore, without considering the capillary end effects, Peaceman well model overestimates the gas production and underestimates the water production under the coarse grid system. But it is indicated that, under the mesh refinement, the calculations of the Peaceman well model gradually approach to those of the proposed well model.

In Test 2, a further analysis of the capillary end effect on the well productions is presented.

5.2 Test 2

In this numerical test, the simulation area is a 2D square reservoir with the size 500 m × 500 m. Parameters are the same as Table 1. A production well is located at the center of the reservoir. In this example, the two-phase flow can be simplified to be 1D radial flow since the wellbore radius is much smaller than the computational domain size. It means that we can get the reference solution by solving the ordinary difference equation of equation (15) for two different flow patterns, respectively. Though actually the simulation area is not radial symmetry, the radially symmetric solution can be seemed as a good reference in the neighborhood of the well. To test the accuracy of the original Peaceman and the proposed well models, the different size grids of 50 × 50, 100 × 100, 200 × 200 and 300 × 300 are adopted in the computations. In the following, first we consider the capillary pressure curve with B = 1 MPa in the calculations for Figures 7 and 8 and then investigated the influence the capillary pressure curve with different values of the constant B in the calculations for Figures 9 and 10.

Figure 7 shows that all the calculations with different grid size are agreed with the reference solution very well when adopting the proposed well model. Nevertheless, adopting Peaceman well model, the calculations of the coarse grid system 50 × 50 deviate from the reference solution even up to 40%. With the mesh subdivision, the result of Peaceman well model is gradually approaching to the reference solution. The calculations indicate that extremely fine grid is needed if using the original Peaceman well model. The capillary end effect has the important influence for the prediction of the low-permeability reservoir recovery. Due to the capillary end effect, the saturation and the gas-phase pressure change rapidly in the vicinity of the production well and the saturation and the gas-phase pressure gradients may even diverge at the entrance of the wellbore (see Fig. 8), which is the same as Figure 2. In the original Peaceman well model, these dramatic changes are not considered and therefore may cause a poor prediction for the low-permeability reservoir recovery.

To investigate the influence of the strength of the capillary pressure upon the predictions of the original Peaceman well model, we perform the calculations with the grid number of 100 × 100 for the different values of the constant B = 0.1 MPa, 1 MPa and 3 MPa. When the capillary pressure is much smaller than the pressure difference between the formation pressure and the well bottom-hole pressure, for instance, in the case of B = 0.1 MPa, Peaceman well model can provide a rather accurate result (see Fig. 9). However, along with the increase of the strength of the capillary pressure, zone affected by capillary pressure is also increasing (see Fig. 10) and therefore the gap between the calculations of Peaceman well model and the proposed model becomes larger and larger. When B = 3 MPa, the water phase is blocked inside the reservoir and then cannot flow out of the production well. In this case, the original Peaceman well model yield an incorrect calculation if the grid size is not enough small.

6 Conclusion

In the recovery of low permeability reservoirs, the capillary pressure has an important effect, which may reduce gas production. When pressure difference between the formation pressure and the well bottom-hole pressure is not large enough to overcome the capillary end effect, only the gas phase can flow into the wellbore and cause the gas production significantly reduced. When the pressure drop is enough to overcome the capillary end effect, both water and gas phases can flow into the wellbore. In this case, the saturation and gas-phase pressure gradients diverge in the vicinity of wellbore. To construct an effective numerical well model for gas reservoir simulation is the main objective of this article. In this article, other near-wellbore processes, such as perforation damage, fines movement, near-wellbore fracturing and wettability changes introduced by stimulation, are not considered. The conclusion can be summarized as follows:

1. Considering the capillary end effect, there are two flow patterns in the vicinity of the production well: both water and gas phase can flow into the production well and only gas phase can enter the production well. With the help of the solution of the 1D radially symmetric flow for two-phase immiscible flow, an alternative well model for low permeability gas reservoir is constructed, where two different flow patterns are considered.

2. When the capillary pressure is much smaller than the pressure difference between the formation pressure and the well bottom-hole pressure, the original Peaceman well model can provide a rather accurate result. The Peaceman well model of multi-phase flow is extended directly from that of single-phase flow, where the capillary end effect is not considered. As a result, along with the increase of the strength of the capillary pressure, the error of the predictions of Peaceman well model becomes larger and larger. The saturation and the gas-phase pressure changes rapidly in the vicinity of the production well and their gradients may even diverge at the entrance of the production well. As a result, the original Peaceman well model can provide good prediction only when the grid size is enough small. This problem is out of control since it is difficult to estimate the suitable grid size for the original Peaceman well model.

3. In contrast to the original Peaceman well model, the proposed well model can provide rather accurate prediction for gas and water productions. The prediction weakly depends on the grid size adopting in the computation, since the dramatic change for saturation and gas-phase pressure in the vicinity of the wellbore has been already considered in the proposed well model. In the recovery of the low-permeability gas reservoirs, the capillary end effect is the main reason for water-phase trapping. The increasing water phase saturation reduces gas mobility near wellbore and then affects well productivity and gas recovery. The proposed well model can accurately describe the capillary end effect in the recovery procedure for low-permeability recovery.

For the single-phase flow, Ding et al. (1998) established the well model based upon the finite-volume method by using the axisymmetric analytical solution. In this article, flow patterns in low-permeability gas reservoirs are obtained from the near well axisymmetric equation. With the help of the solution equations (31)–(35) of the 1D radially symmetric flow for the two-phase flow, the well model for low permeability gas reservoir in cartesian mesh is established, based upon Ding’s method. It is expected that the proposed method can be extended to flexible grids based upon the near-well axisymmetric solutions. This remains the future research.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (No. 11572315).

References

All Tables

All Figures