An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method

Computers and Mathematics with Applications xxx (xxxx) xxx

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method ∗

Shuoliang Wang a , , Chunlei Yu b , Guoqiang Sang c , Rongze Yu d , Feng Cheng d a

School of Energy, Faculty of Engineering, China University of Geosciences, Beijing 100083, PR China Shengli Oilfield Exploration and Development Research Institute, Dongying, Shandong 257000, PR China State Key Laboratory of Enhanced Oil Recovery, Beijing 100083, PR China d PetroChina Research Institute of Petroleum Exploration and Development, Beijing 100083, PR China b c

article

info

Article history: Received 25 October 2018 Received in revised form 3 November 2019 Accepted 21 November 2019 Available online xxxx Keywords: Dynamic capillary force Reservoir numerical simulation Finite difference Sweep coefficient Water cut

a b s t r a c t Capillary pressure is a key parameter in reservoir numerical simulation. It is generally considered that the capillary force is an only function of water saturation in the current traditional reservoir numerical simulation methods. However, large amount of indoor experiments have shown that capillary force is not only dependent on the water saturation, but also the velocity of fluid flow. In this work, we introduce the dynamic capillary pressure to the oil–water two-phase flow equations. Then, we discrete the modified oil and water phase basic governing equation with difference methods. The full implicit method is selected to solve the pressure equation. Finally, a numerical simulation method coupled with the dynamic capillary force based on the full-implicit method is established. We discussed the effect of dynamic capillary pressure on the oil–water two-phase flow in two-dimensional homogeneous and heterogeneous porous media. Comparing to static capillary pressure model, the dynamic capillary pressure model predicts higher water-saturation near the injection and production wells. In addition, a ‘‘funnel’’ like pressure profile appears in the water saturation profile. At low water cut stage, the rate of water cut increase is relatively high, while opposite is true at high water cut stage. We also calibrated our newly proposed model by comparing to the analytical method and the actual oil field production data. Our work should provide important insights into the effect of dynamic capillary force on the oil–water two-phase flow and better estimation of flow behavior in oil and gas reservoir. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Capillary pressure (Pc) is an important parameter in petro-physics of oil and gas development [1]. The capillary pressure curve of reservoir reveals the pore size and the pore throat volume during the immiscible flow. Thus, it can be used to analyze the microscopic heterogeneity of reservoir [2]. In addition, the depth of the fluid interface (oil–water interface and oil–gas interface) and the shape of water flooding front which play an important role in the development plan in oil field can be determined by the capillary force curve [3]. Capillary pressure is also an important parameter in reservoir numerical simulation. It can be used to analyze remaining oil distribution in the reservoir and predict the recovery factor ∗ Correspondence to: Petroleum Engineering China University of Geosciences, PR China. E-mail address: [email protected] (S. Wang). https://doi.org/10.1016/j.camwa.2019.11.013 0898-1221/© 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

2

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

of oil field. Therefore, accurate representation of capillary pressure curve plays a key role in reservoir description and oil field development. In the conventional reservoir simulation method, the common belief is that such capillary pressure is the pressure difference between the wetting phase and non-wetting phase, considered as static capillary pressure (SCP). The pressure difference between the wetting phase and non-wetting phase is equal to capillary pressure under equilibrium conditions. SCP is constant when the oil–water interface reaches equilibrium. Within reservoir numerical simulation and classic petrophysics, it is usually considered that the capillary pressure is only dependent on water saturation [4,5]. The capillary pressure measured by mercury intrusion method, centrifugal method and semipermeable membrane method is actually SCP. Flow velocity is not considered in the aforementioned capillary pressure measurement methods. There have been a number of experimental studies on the SCP. Bear and Verruijt [6] indicated that SCP is only related to the relative saturation of the wetting phase. The most popular methods to calculate SCP are Brooks and Corey method [7] and van Genuchten method [8]. The capillary force used in reservoir numerical simulation method has been commonly considered to be a function of the wetting-phase saturation only [9–13], without considering fluid flow velocity [14,15]. However, dynamic motion behavior was detected during the fluid flow [16]. In actual reservoir, the oil and water seepage is always in motion. In such unsteady state, the oil–water interface equilibrium is disturbed. Under the nonequilibrium conditions, the pressure difference between the wetting phase and non-wetting phase is the dynamic capillary pressure (DCP) that is completely different from the conventional static capillary pressure [17]. There have been a number of studies on DCP between gas and water two-phase models in hydrodynamics [18–20]. Experimentally, Kalaydjian [21] conducted a two-phase (oil and water) flow experiment and found that the DCP in water–oil phase flow in porous medium is very different from the SCP (when water saturation = 28%, SDP = 2 KPa, CDP = 10 KPa). Theoretically, Hassanizadeh and Gray [22] proposed a new dynamic empirical coefficient to describe the difference between SCP and DCP. Later on, Hassanizadeh [23] summarized the range of dynamic empirical coefficient (5×104 ∼6×104 Pa · s) through numerous experiments, indicating large deviations of DCP from SCP. Camps-Roach [24] showed that the DCP is larger than the SCP at a given saturation. To account for such DCP, many researchers [25,26] developed numerical approaches based on the material balance. They found that DCP is not only dependent on static factors, such as porosity, permeability and fluid properties, but also dynamic properties (flow fluid velocity. etc.). There have been a number of works incorporating DCP into numerical simulation [27–29]. Sandoval-Torres [30] used finite element method to study wood vacuum drying process considering DCP, showing excellent agreement with experimental results. Although there has been a number of studies of dynamic capillary pressure model in hydrodynamics, its effect on the reservoir numerical simulation in petroleum is rarely studied. Shu [31] used the same method in the oil– water two-phase flow in one-dimensional numerical simulation and studied the effect of the DCP on the oil and water migration law. Tian [32] incorporated dynamic capillary force into the one-dimensional numerical simulation by using IMPES (The implicit pressure–explicit saturation method) method considering DCP with one injection and one production well. They found that DCP could greatly affect the water breakthrough time, i.e., 279 day under SCP condition comparing to 417 day under DCP condition. In summary, the capillary force involved in the process of fluid mechanics in porous medium and reservoir numerical simulation should be a DCP associated with the fluid flow velocity. However, there are many production wells and injection wells in actual oil and gas reservoirs. The velocity near the production wells and injection wells is much larger than that in inter-well region. The numerical simulation model in hydrodynamics cannot be used to simulate actual oil and gas reservoir. The reservoir numerical simulation model introduced earlier does not take into account the effect of sudden increase of flow velocity near production and injection wells on reservoir saturation and pressure, and the solution method is not full implicit. For this reason, the calculation results are only compared with the results of laboratory experiments, have not been calibrated by comparing to actual productions. In this study, we couple the dynamic capillary pressure to three-dimensional reservoir scale numerical simulation to study oil–water two-phase flow. DCP is characterized by the method developed by Hassanizadeh and Gray [22], which demonstrates that the capillary pressure constantly changes during the unsteady state motion, which is not only a function of wet-phase fluid saturation, but also influenced by flow rate. Initially, the dynamic capillary pressure is firstly introduced into the fundamental governing equation in a reservoir scale and fully implicit method is used to calculate saturation and pressure. Then, we also compare the numerical simulations with static and dynamic capillary pressure in two-dimensional homogeneous and heterogeneous model. In the last section, the calculation results are compared not only with the analytical method but also with the production data from an offshore field in China. Our newly proposed model shows excellent reliability and accuracy. 2. Mathematical model In this section, we present the mathematical model describing the two-phase flow in a three-dimensional porous medium. Based on the oil–water two-phase flow model, SCP is replaced by DCP in the governing equation. The preconditions of the oil–water two-phase flow model are listed as: (1) There are exactly two phases (water, oil) in the reservoir. (2) There is no mass exchange between oil and water phases. Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

3

(3) Fluid flow is regarded as Darcy flow. (4) The process is isothermal. (5) Instantaneous phase equilibrium is observed during the two-phase flow. 2.1. Darcy’s equation Oil phase:

vo = −

kro

µo

k (∇ Po − γo ∇ D)

(2.1)

Water phase:

vw = −

kr w

µw

k (∇ Pw − γw ∇ D)

(2.2)

where vo is volumetric velocity of oil phase, kro is the oil phase relative permeability, k is the absolute permeability tensor,

µo , Po , and γo is the viscosity, pressure, and specific gravity of oil phase, respectively, vw is volumetric velocity of water phase, kr w is the water phase relative permeability, µw , Pw , and γw is the viscosity, pressure, and specific gravity of water phase, D is the depth. 2.2. Mass conservation equation Oil phase:

( →) ∂ (φρo So ) −∇ ρo − vo + qo = ∂t

(2.3)

Water phase:

( →) ∂ (φρw Sw ) −∇ ρw − vw + qw = ∂t where ρo , So , qo is the density, saturation, and the flow rate of oil phase, respectively, φ is the porosity.

(2.4)

2.3. Governing equation Oil phase:

[ ∇·

kkro ρo

µo

]

(∇ po − γo ∇ D) − qo =

∂ (ρo φ So ) ∂t

(2.5)

Water phase:

[ ∇·

kkr w ρw

µw

]

(∇ pw − γw ∇ D) − qw =

∂ (ρw φ Sw ) ∂t

(2.6)

where, γo = ρo g, γw = ρw g 2.4. Auxiliary equation

{ { { {

So + Sw = 1 Po − Pw = Pcowd

(2.7)

kro = kro (Sw ) kr w = kr w (Sw )

(2.8)

ρo = ρo (Po ) ρw = ρw (Pw )

(2.9)

µo = µo (Po ) µw = µw (Pw )

(2.10)

where Pcowd is the dynamic capillary pressure. In the traditional numerical simulation method, the capillary force Pc is given by a function of saturation of wetting phase, and we call it SCP (Pcow ). In this section, we replace the SCP with DCP (Pcowd ). Pcowd = Po − Pw = Pcow − τ (

∂ Sw ) ∂t

(2.11)

Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

4

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

In 1978, Stauffer [33] proposed the empirical expression of dynamic parameter:

τ=

α s µ w φ Pe 2 ( ) Kλ ρw g

(2.12)

Where αs is empirical constant, which is equal to 0.1; µw is flow velocity of the wetting phase, (m/s); φ is the porosity, (%); K is the absolute permeability, (10−3 µm2 ); λ is factors of capillary pressure–saturation relationships in the Brook–Corey model; Pe is the capillary pressure of pressure–saturation relationship in the Brook–Corey model, (MPa); ρw is density of the wetting phase, (kg/m3 ); g is gravity acceleration, (m/s2 ). 3. Numerical model In this section, the governing equations are initially discretized, while full-implicit iteration is used to calculate pressure and saturation. 3.1. Discretization of the governing equation In this work, DCP is introduced into the discretization of governing equation. The specific discrete procedures are listed in Appendix A. The relationship between DCP and SCP can be written as, n Pon+1 − Pwn+1 = Pco w − τ(

n+1 Sw − Swn ∂ Sw n ) = Pco w −τ ∂t ∆t

(3.1)

The numerical simulation model developed in this paper is an oil–water two-phase model. The absolute permeability used in this paper is anisotropy (the absolute permeability of XYZ directions is different). The governing equation of this simulator is anisotropy. Therefore, the velocity, water saturation, pressure calculated by this simulator of XYZ directions are also different. Although the dynamic capillary pressure equation is one-dimensional formula, the dynamic capillary pressure of XYZ directions is different. The dynamic capillary pressure is a function of absolute permeability, water saturation, pressure and velocity. The basic parameters in dynamic capillary pressure equation are three-dimensional. Therefore, the dynamic capillary pressure calculation results are three-dimensional. The oil phase governing equation is discretized as a sample in this section, which is shown as follows,

∆Ton+1 ∆P n+1 − ∆Ton+1 γog ∆D + qno+1 Vijk =

Vijk

∆t n+1

[(φρo So )n+1 − (φρo So )n ]

(3.2)

The governing equation of the water phase is expressed with the same format. S n +1 − S n

w n+1 w ∆Twn+1 ∆P n+1 − ∆Twn+1 ∆(Pco ) − ∆Twn+1 γw ∆Dn+1 w −τ ∆t V ijk + qnw+1 Vijk = ∆t n+1 [(φρw Sw )n+1 − (φρw Sw )n ]

Φ is used to present the capillary pressure during the calculation. { n Φo = pno − γo D Φwn = pnw − γw D = pno − (pncow − τ

n −S n−1 Sw w ) ∆t

− γw D

(3.3)

(3.4)

To summarize, the governing equation discretization of the oil phase is expressed as follows,

∆Ton+1 ∆Φon+1 + qno+1 Vijk =

Vijk

∆t n+1

[(φρo So )n+1 − (φρo So )n ]

(3.5)

Similarly, the governing equation discretization of the water phase can be written as:

∆Twn+1 ∆Φwn+1 + qnw+1 Vijk =

Vijk

∆t n+1

[(φρw Sw )n+1 − (φρw Sw )n ]

(3.6)

3.2. Approximation of pressure and saturation In this study, X is set as an arbitrary variable during multiple iterations. In random iteration from Step l to Step l + 1, the X value can be assigned as X l and X l+1 . When liml→n X l+1 − X l < δ , X l+1 → X n+1 is the solution of Step n. The increment between two iteration steps is δ X = X l+1 − X l , while the step length is δ X = X l+1 − X n = X l + δ X − X n . Substituting these equations into Eqs. (3.5) and (3.6), the full-implicit iteration equations of oil phase (3.7) and water phase (3.8) are given as follows,

∆Tol+1 ∆Φol+1 + qno+1 VM = ∆Twl+1 ∆Φwl+1 + qnw+1 VM =

VM

∆t n+1 VM

∆t n+1

[ ] δ (ϕρo So ) + (ϕρo So )l − (ϕρo So )n

(3.7)

[ ] δ (ϕρw Sw ) + (ϕρw Sw )l − (ϕρw Sw )n

(3.8)

Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

5

The relationship between pressure and saturation is expressed by conductivity To as follows,

δ To =

∂ To ∂ To δ Po + δ Sw ∂ Po ∂ Sw

(3.9)

γo is calculated by explicit difference schemes, while γol+1 = γol .δ Φo can be written as follow, δ Φo = Φol+1 − Φol = (Po − γo D)l+1 − (Po − γo D)l = Pol+1 − Pol = δ Po

(3.10)

The unknowns in the right side of Eq. (3.7) can be expressed as follows,

δϕ = ϕ0 CR δ Po δρo =

∂ρo So δ Po ∂ Po

δ So = −δ Sw

(3.11) (3.12) (3.13)

The full-implicit iteration Eq. (3.7) of oil phase can be written as,

∂ To ∂ To δ Po ∆Φol + ∆ δ Sw ∆Φol + qno+1 VM ∂ Po ∂ Sw VM ∂ρol VM = (ρol Sol ϕ0 CR + ϕ l Sol )δ Po − ϕ l ρol δ Sw n + 1 ∆t ∂ Pol ∆t n+1 ] VM [ + n+1 (ϕρo So )l − (ϕρo So )n ∆t ∆Tol ∆Φol + ∆Tol ∆δ Po + ∆

(3.14)

The discretization Eq. (3.14) of oil phase using full-implicit iteration method is listed in Appendix B. Three-dimension grid (i, j, k) is replaced with m., while upstream scheme is used to induce all the auxiliary points (i + 21 , j + 12 , k + j + 21 ) to the corresponding node. Eq. (3.14) can be written as,

∂ Tom−1 l ∂ Tom−1 l ∂ Tom−1 l ∂ Tom−1 l δ Pom−1 (Tom−1 − Φom + Φom−1 ) + δ Swm−1 (− Φom + Φ ) ∂ Pom−1 ∂ Pom−1 ∂ Swm−1 ∂ Swm−1 om−1 [ ] ∂ Tom l ∂ Tom l VM ∂ρol l l + δ Pom −Tom − Tom−1 + Φom+1 − Φom − (ρol Som ϕ0 CR + ϕ l Som ) n + 1 l ∂ Pom ∂ Pom ∆t ∂ Pom ∂ Tom l VM ∂ Tom l l Φ − Φ + ϕ l ρom ) + δ Pom+1 (Tom ) + δ Swm+1 (0) + δ Sw m ( ∂ Swm om+1 ∂ Swm om ∆t n+1 ] VM [ (ϕρom Som )l − (ϕρom Som )n − ∆Tol ∆Φol − qno +1 VM = ∆t n+1

(3.15)

In this study, we assume that,

⎧ ∂ Tom−1 l ∂ Tom−1 l ⎪ ⎪ b1om = Tom−1 − Φom + Φ ⎪ ⎪ ∂ Pom−1 ∂ Pom−1 om−1 ⎪ ⎪ ⎪ ⎪ ∂ Tom−1 l ∂ Tom−1 l ⎪ ⎪ b2om = − Φom + Φ ⎪ ⎪ ⎪ ∂ Swm−1 ∂ Swm−1 om−1 ⎪ ⎪ ⎪ ∂ Tom l ∂ Tom l ⎪ 1 ⎪ ⎪ coim = −Tom − Tom−1 + Φom+1 − Φ ⎪ ⎪ ∂ P ∂ Pom om ⎪ om ⎪ ⎨ VM ∂ρol l l − n+1 (ρol Som ϕ0 CR + ϕ l Som ) l ⎪ ∆t ∂ Pom ⎪ ⎪ ⎪ ⎪ ∂ Tom l ∂ Tom l VM ⎪ 2 l l ⎪ ⎪ ⎪ coim = ∂ S Φom+1 − ∂ S Φom + ∆t n+1 ϕ ρom ⎪ wm wm ⎪ ⎪ ⎪ d1 = T ⎪ ⎪ om om ⎪ ⎪ ⎪ 2 ⎪ d = 0 ⎪ om ⎪ ⎪ [ ] ⎪ ⎩ gom = VM (ϕρom Som )l − (ϕρom Som )n − ∆T l ∆Φ l − qn+1 VM o o o n + 1 ∆t

(3.16)

So, the solution matrix of the oil phase equation can be simplified as: 1 2 b1om δ Pom−1 + b2om δ Swm−1 + com δ Pom + coim δ Swm + d1om δ Pom+1 + d2om δ Swm+1 = gom

(3.17)

Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

6

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

According to the discretization process of the oil phase equation, the water phase equation considering dynamic capillary pressure can be written as follows:

∆Twl ∆Φwl + ∆Twl ∆δ Po + ∆

∂ Pcowe ∂ Tw δ Po ∆Φwl − ∆Twl ∆ δ Sw ∂ Po ∂ Sw

τ ∂ Tw ∆T l ∆δ Sw + ∆ δ Sw ∆Φwl + qnw+1 VM ∆t w ∂ Sw l VM l l l ∂ρw = (ρw Sw ϕ0 CR + ϕ l Sw )δ Pw n + 1 l ∆t ∂ Pw ] VM l l VM [ − ϕ ρw δ S w + (ϕρw Sw )l − (ϕρw Sw )n ∆t ∆t −

(3.18)

In this study, we assume that,

⎧ ∂ Twm−1 l ∂ Twm−1 l ⎪ Φwm + Φ b1wm = Twm−1 − ⎪ ⎪ ⎪ ∂ P ∂ Pom−1 wm−1 om−1 ⎪ ⎪ ⎪ ∂ T ∂ Pcowe α l ∂ T wm−1 l w m−1 l ⎪ ⎪ Φwm + Φwm−1 − Twl m−1 − T b2wm = − ⎪ ⎪ ∂ Swm−1 ∂ Swm−1 ∂ Swm−1 ∆t wm−1 ⎪ ⎪ ⎪ ⎪ ∂ Twm l ∂ Tw m l ⎪ 1 ⎪ cw Φwm+1 − Φ ⎪ m = −Tw m − Tw m−1 + ⎪ ⎪ ∂ Pom ∂ Pom wm ⎪ ⎪ l ⎪ VM ∂ρ ⎪ ⎪ ⎨ − n+1 (ρwl Swl m ϕ0 CR + ϕ l Swl m l w ) ∆t ∂ Pw m ∂ Tom l ∂ Tom l VM ⎪ 2 ⎪ Φwm+1 − Φwm + ϕ l ρwl m ⎪ cwm = n+1 ⎪ ∂ S ∂ S ∆ t ⎪ w m w m ⎪ ⎪ ⎪ ∂ Pcowe ∂ Pcowe α l α l ⎪ ⎪ + Twl m + Twl m−1 + T + T ⎪ ⎪ ⎪ ∂ Swm ∂ Sw m ∆t wm ∆t wm−1 ⎪ ⎪ 1 ⎪ dw m = Tw m ⎪ ⎪ ⎪ ⎪ ⎪ d2 = −T l ∂ Pcowe − α T l ⎪ wm wm ⎪ ⎪ ∆t w m ⎪ [∂ Swm+1 ] ⎩ VM l gwm = ∆t n+1 (ϕρwm Swm )l − (ϕρwm Swm )n − ∆Twl ∆Φw − qnw+1 VM

(3.19)

So, the solution matrix of the water phase equation can be simplified as: 1 2 1 2 b1wm δ Pom−1 + b2wm δ Swm−1 + cw m δ Pom + cw m δ Sw m + dw m δ Pom+1 + dw m δ Sw m+1 = gw m

(3.20)

The order of nodes in the equation is arranged according to the order of nodes in the numerical simulation model. The order of unknown quantities at nodes is arranged according to the order of (δ Pom , δ Sw ). A group of differential equations on any node can be obtained after rearranging the differential equation. All nodes and unknown quantities are sorted naturally. u = δ Pom1 , δ Sw1 , δ Pom2 , δ Sw2 , δ Pom3 , δ Sw3 , . . . . . . . . . . . . , δ PomN , δ SwN

[

]T

(3.21)

The difference equations on all grids in the entire solution area are as follows:

− →

− →

Au = G.

(3.22)

− →

The coefficient matrix A is a diagonal form, G is a constant column vector.

⎡

1 co1

2 co1

d1o1

d2o1

1

2

d1o1 1 co2 1 cw2 b1o3 b1w3

d2o1 2 co2 2 cw2 b2o3 b2w3

⎢ cw1 ⎢ ⎢ b1 ⎢ o2 ⎢ b1 ⎢ w2 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

cw1 b2o2 2

bw 2

d1o2

d2o2

1

2

dw 2

dw2 2 co3

d1o3

d2o3

1

2

1

dw3

d2w3

1 coN

1 coN

d1oN

d1oN

1 cw N

1 cw N

d1wN

d1wN

1 co3

cw3

.. .

cw3

.. .

.. .

.. .

⎤ ⎡ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥·⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

δ Po1 δ Sw1 δ Po2 δ Sw2 δ Po3 δ Sw3 .. .

δ PoN δ Sw N

⎤

⎡ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥=⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

go1 gw 1 go2 gw2 go3 gw3

.. .

goN gw N

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

(3.23)

Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

7

Fig. 1. Dependence of relative permeability on the water saturation.

Fig. 2. Dependence of dynamic and static capillary force on the oil saturation.

Table 1 Simulation model parameters. Grid node

40*40*1

Dx (m)

2

Dy (m) Top deep (m) Porosity (f) Water viscosity (mPa · s) Exploit scheme

2 1000 0.2 0.5 Water flooding

Dz (m) Initial water saturation (f) Permeability (10−3 µm2 ) Oil viscosity (mPa · s) Well pattern

1 0.365 500 10 One injection well and one production well

4. Numerical experiments 4.1. Numerical simulation parameters This section mainly discusses the difference between SCP simulation method and DCP numerical simulation method. The basic parameters of the numerical simulation model are as shown in Table 1, and the relative permeability curve used in the numerical simulation model is shown in Fig. 1. The measured experimental capillary pressure conducted by Fuík R [34] is used to analyze the results by considering DCP. Fig. 2 shows DCP and SCP obtained in the experiments. Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

8

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 3a. Contour plot of water saturation before water break through from SCP.

In this subsection, we compare the results from the numerical simulator which uses SCP and that uses DCP in the homogeneous model. The saturation field are obtained in three different time steps: (1) before water break through; (2) water break through to the bottom of production well; (3) 95% water cut. 4.2. Planar homogeneous model numerical simulation results 4.2.1. Water saturation comparison before the water break through In Figs. 3, we present the water saturation field before the water break through from SCP and DCP. Before water breakthrough (Figs. 3, 4), the difference after taking DCP into consideration is mainly reflected near water well. The characteristic of radial fluid flow is that the pressure dropping speed is fast around the injection well and the production well, and the large pressure gradient will lead to a faster fluid velocity, which will change the capillary pressure. The area near the water well, the water saturation value calculated by DCP method is higher; The interwell area, due to the slower fluid flow rate between the oil well and water well, there are not too many differences between two methods; The area near the oil well, as the injection water does not breakthrough, there is no change in the water saturation. For clarity, we also present the water saturation field along the diagonal direction between the injection and production well as shown in Fig. 4. 4.2.2. Water saturation comparison when water break through to the bottom of production well After water breakthrough (Figs. 5, 6) near the water well, the water saturation and saturation gradient of DCP method are both high nearby the water well; In the interwell area, the water saturation is almost as same as the value before water breakthrough; In the area near the oil well, the water saturation value calculated by DCP method is higher. It is indicated from the water saturation field that the trade of water saturation value near the oil well is slightly upturn. 4.2.3. Water saturation comparison when water cut of production well is 95% When the water cut arrives 95% in oil well (Figs. 7, 8), there is a relatively large difference in water saturation between DCP and SCP because the water saturation profile curve upwards near the oil well. The water saturation and saturation gradient are both high around the water well and oil well. The water saturation value along the connection of water well and oil well is lower than that near water well and oil well. A ‘‘funnel’’ similar to the pressure profile appears in the water saturation profile. Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

9

Fig. 3b. The same as Fig. 3a, but from DCP.

Fig. 4. Water saturation profile along the diagonal direction between the injection and production wells before water break through from SCP and DCP.

4.2.4. Water cut curve comparison We will compare the water cut curve of SCP method and DCP numerical simulation method in this part. The water breakthrough time calculated by the SCP method is the same as that calculated by the DCP method. The main difference of the two water cut curves lies in the fast-rising water cut period. The water cut calculated by the DCP numerical simulation method is slightly higher than that obtained by SCP method. In high water cut stage, the water cut calculated by DCP numerical simulation method is rising slowly, which is closer to the water cut calculated by the SCP method, but it is always higher than that of the SCP method. (Fig. 9). Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

10

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 5a. Contour plot of water saturation field at water break through from SCP.

Fig. 5b. The same as Fig. 5a, but from DCP.

Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

11

Fig. 6. Water saturation profile when water break through from SCP and DCP.

Fig. 7a. Contour plot of water saturation at 95% water cut from SCP.

4.3. Planar heterogeneous model numerical simulation results In order to analyze the effect of DCP on the simulation results under plane heterogeneous conditions, a high permeability channel is set up between the water well and the production well. The permeability of the channel is 500×10−3 µm2 . The permeability horizontal distribution is as shown in Fig. 10. The simulator is set up in two patterns to analyze the influence of DCP on numerical simulation results. One uses the DCP simulator and the other uses SCP simulator. The saturation field is obtained in three different time steps. The three different time steps are: before water break through, the water break through point and 95% water cut point. Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

12

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 7b. The same as Fig. 7a, but from DCP.

Fig. 8. Water saturation profile at 95% water cut from SCP and DCP.

Water saturation profile value is extracted along the connection direction of the injection well and the production well in different recovery time as follows.

4.3.1. Water saturation comparison before the water break through Before water breakthrough (Figs. 11, 12), in the area near the water well, the water saturation value calculated by the DCP simulation method is higher than that calculated by the SCP simulation method near the water well. In the interwell area and the area near oil well, there are not too many differences. Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

13

Fig. 9. Water cut calculated from SCP and DCP.

Fig. 10. Permeability (the unit is 10−3 µm2 ) horizontal distribution in the heterogeneous model.

4.3.2. Water saturation comparison when water break through to the bottom of production well After water breakthrough (Figs. 13, 14), in the area near water well, the water saturation and water saturation gradient of DCP simulation method are both higher than that calculated by SCP simulation method; In the area near oil well, the water saturation profile curve is slightly upturn. 4.3.3. Water saturation comparison when water cut of production well is 95% When the water cut arrives 95%, (Figs. 15, 16), there is a relatively large difference in water saturation because the water saturation profile curve upwards in the vicinity of the oil well. Both of the water saturation and saturation gradient of DCP numerical simulation method are higher near water well and oil well than SCP numerical simulation method. The saturation value in the interwell area is lower than that near injection well and production well. A ‘‘funnel’’ similar to the pressure profile appears in the water saturation profile. 4.3.4. Water cut curve comparison As Fig. 17 shown, the water breakthrough time calculated by the SCP numerical simulation method is the same as that calculated by the DCP method. The main difference performs at the stage of rapid rising in water cut rate. Compared with Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

14

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 11a. Contour plot of water saturation in heterogeneous model before water breakthrough from SCP.

Fig. 11b. The same as Fig. 11a, but from DCP.

homogeneous reservoir, the water cut difference is quite small in heterogeneous reservoir. The reason is that the injection water flows mainly along the high permeability channel in heterogeneous reservoir. The saturation changes dramatically Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

15

Fig. 12. Water saturation profile before water break through from SCP and DCP.

Fig. 13a. Contour plot of water saturation in heterogeneous model at water breakthrough from SCP.

in the high permeability region. The influence of saturation on the capillary force in high permeability region is greater than that in other region. 5. Validation of simulation results In order to calibrate our simulation results, analytical method (fractional flow model) is used to study the influence of dynamic capillary on the water cut in this research. The calculation formula of water cut is shown as follows: fw =

{ [ ]/ } λw ∂ pc ∂ Sw 1 + λo A − (ρo − ρw )g sin α qt λw + λo ∂ Sw ∂ x

(5.1)

Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

16

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 13b. The same as Fig. 13a, but from DCP.

Fig. 14. Water saturation profile when water break through from SCP and DCP.

where fw is the water cut rate, λw is the water mobility, k · kr w /µw, md/mPa · s, λo is the oil mobility, k · kro /µo , md/mPa · s, A is core cross sectional area, m2 , Pc is capillary pressure, MPa, Sw is water saturation, %, qt is the total volume flow of water and oil, m3 /s, ρo is oil density, g/cm3 , ρw is water density, g/cm3 . From Eq. (5.1), it can be seen that the water cut is related to the slope of the capillary force curve. As Fig. 18 shown, there flooding velocities are discussed in this study. The slope of DCP curves shows different characteristics at various saturation conditions. When So > 50%, the higher the displacement velocity is, the greater the slope of the DCP curve and the quicker water cut rise. When So < 50%, it is opposite. The slope of dynamic capillary is always larger than the static capillary when the velocity is less than 0.08ml/min. When So > 57%, the slope of the dynamic capillary (velocity = 0.25 ml/min) is lower than the static one. Elsewhere, the slope of the DCP curve is much greater than the slope of the SDP curve. From the analytical method (fractional flow Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

17

Fig. 15a. Contour plot of water saturation in heterogeneous model at 95% water cut from SCP.

Fig. 15b. The same as Fig. 15a, but from DCP.

Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

18

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 16. Water saturation profile at 95% water cut from SCP and DCP.

Fig. 17. Heterogeneous model water cut curve calculated from SCP and DCP.

Fig. 18. Capillary force slope calculated by analytical method (fractional flow model).

Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

19

Fig. 19. Permeability distribution in A6H well.

Fig. 20. Water cut rate in the single well of China offshore oil field. Table 2 Specific reservoir parameters of production well A6H. Parameter

Value

Parameter

Value

Thick of reservoir Average porosity Viscosity of water Viscosity of oil

15 m 0.28 0.65 mPa · s 1.2mPa · s

Average permeability Permeability variation coefficient Average water saturation Lithology

850md 0.68 0.54 Unconsolidated sand

model), it is indicated that water cut of the DCP method is higher than that of the SDP method, which is consistent with the simulation result. We also compare our model with an actual production well in China offshore oil field. The specific reservoir parameters of well A6H are shown in Table 2. The reservoir has high water aquifer energy. In the stage of medium and high water cut, large extraction liquid pump is used to enhance the liquid of production well. After the single well enhanced liquid treatment, the flow velocity in the bottom of the well increased rapidly, and the capillary force has changed accordingly. Therefore, production well A6H in China offshore oil field is selected to validate the accuracy of DCP numerical simulation method. Based on the geological background of A6H well, three-dimensional geological model is established as Fig. 19 shown. Two different methods are used to simulate in A6H well. The water cut rate in the single well is shown in Fig. 20. It can be seen from Fig. 20 that the water cut is equal to 0 at multiple time points. The water cut is equal to 0 indicating shutting down well. Three enhanced liquid treatments were carried out in November 2011, December 2012 and December 2014. The liquid rate was enhanced from 942 m3 /day to 1254 m3 /day for the first time, from 1200 m3 /day to 1568m3 /day for the second time, from 1511 m3 /day to 1785 m3 /day. The first two measures were very effective, and the water cut curve showed a marked decline. The results of the third enhanced liquid treatment was not good, and the water cut curve increased significantly. The SCP considers many enhanced liquid mechanisms, such as starting pressure. From the Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

20

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

comparison of the actual and simulation water cut curve in A6H well, the water cut curve calculated by SCP numerical simulation method reaches plateau and not affected by the enhanced liquid treatment. The water cut curve calculated by the DCP numerical simulation method is closer to the actual production curve than SCP numerical simulation method. It shows that the DCP numerical simulation method has better accuracy and reliability comparing to field data. 6. Conclusion In this study, the dynamic capillary force is coupled with the numerical simulation. New governing equation is established, and finite-difference algorithm is used to solve water phase and oil phase governing equation. Full implicit method is utilized in solving saturation and pressure. We studied the difference between reservoir numerical simulation with static and dynamic capillary pressure in twodimensional homogeneous and heterogeneous model. Comparing to static capillary pressure model, the dynamic capillary pressure model predicts higher water-saturation near the injection and production wells. A ‘‘funnel’’ like pressure profile appears in the water saturation profile due to dynamic capillary pressure. In order to validate the accuracy of simulation results, the effect of capillary on water cut is analyzed in fractional flow model. It shows that numerical simulation results considering the dynamic capillary force are consistent with the fractional flow theory. Moreover, our model shows better agreement with the water cut from an actual production well than traditional numerical simulation method. Our works have shown that it is necessary to consider dynamic capillary pressure effect in numerical simulation and provided important insights into accurate estimation of well performance. Acknowledgment This work was supported by the National Science and Technology Major Project (Grant Nos. 2017ZX05009001 and 2017ZX05035-004-005). And this work was Funded Fundamental Research Funds for the Central Universities (2-9-2018-217). Appendix A. Discretization of the governing equation The oil phase governing equation is expanded initially into a rectangular coordinate component as follows,

∂ ∂D ∂D ∂ po ∂ ∂ po [λo ( − γo )] + [λo ( − γo )] ∂x ∂x ∂x ∂y ∂y ∂y ∂D ∂ po ∂ (φρo So ) ∂ − γo )] + qo = + [λo ( ∂z ∂z ∂z ∂t where, λl =

(A.1)

ρl kkrl µl

λl is the flow coefficient. The point (i, j, k, n + 1) could be abbreviated as the form of i + 1

∆xi

+ +

[λoi+ 1 (

pni++11 − pni +1

∆xi+ 1

2

− γoi+ 1

Di+1 − Di

∆xi+ 1

2

2

1

∆yj 1

∆zk

[λoj+ 1 (

pnj++11 − pnj +1

2

∆yj+ 1

[λok+ 1 (

+ qno+1 =

2

1

∆t

∆yk+ 1

pni−+11 − pni +1

∆xi− 1

2

2

− γoj+ 1

∆yj+ 1

2

) + λoj− 1 (

− γok+ 1

2

∆xi− 1

∆xj− 1

∆yk+ 1

) + λok− 1 ( 2

2

− γoj− 1

)]

Dj−1 − Dj

∆ xj− 1

2

+1 n+1 pnk− 1 − pk

∆xk− 1

= i + 21 , j, k.

2

2

Dk+1 − Dk

2

Di−1 − Di

2

pnj−+11 − pnj +1

2

2

− γoi− 1

2

Dj+1 − Dj

2

+1 n+1 pnk+ 1 − pk

) + λoi− 1 (

1 2

)] (A.2)

2

− γok− 1 2

2

Dk−1 − Dk

∆xk− 1

)]

2

[(φρo So )n+1 − (φρo So )n ]

Multiply both sides by Vijk = ∆xi ∆yj ∆zk and the following coefficient of conduction is defined as follows,

⎧ Vijk λoi+ 12 ∆yj ∆zk Vijk λoi− 21 ∆yj ∆zk ⎪ ⎪ TXoi+ 1 = = λoi+ 1 , TXoi− 1 = = λ 1 ⎪ ⎪ 2 2 2 ⎪ ∆xi ∆xi+ 1 ∆xi+ 1 ∆xi ∆xi− 1 ∆xi− 1 oi− 2 ⎪ ⎪ 2 2 2 2 ⎪ ⎨ Vijk λoj+ 12 ∆xi ∆zk Vijk λoj− 21 ∆xi ∆zk TXoj+ 1 = = λoj+ 1 , TXoj− 1 = = λ 1 2 2 2 ∆yj ∆yj+ 1 ∆yj+ 1 ∆yj ∆yj− 1 ∆yj− 1 oj− 2 ⎪ ⎪ 2 2 2 2 ⎪ ⎪ ⎪ ⎪ V λok+ 12 ∆xi ∆yj Vijk λok− 21 ∆xi ∆yj ⎪ ⎪ TXok+ 1 = ijk = λok+ 1 , TXok− 1 = = λ 1 ⎩ 2 2 2 ∆zk ∆zk+ 1 ∆zk+ 1 ∆zk ∆zk− 1 ∆zk− 1 ok− 2 2

2

2

(A.3)

2

The two-order difference operator is defined as follows,

⎧ ⎪ ⎨ ∆x TXo ∆x P = TXoi+ 21 (pi+1 − pi ) + TXoi− 12 (pi−1 − pi ) ∆y TXo ∆y P = TXoj+ 1 (pj+1 − pj ) + TXoj− 1 (pj−1 − pj ) 2 2 ⎪ ⎩ ∆z TXo ∆z P = TX 1 (pk+1 − pk ) + TX 1 (pk−1 − pk ) ok− ok+ 2

(A.4)

2

Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

21

Eq. (A.4) can be simplified to Eq. (A.5),

∆x TXo ∆x P n+1 + ∆y TYo ∆y P n+1 + ∆z TZo ∆z P n+1 − ∆x TXo γog ∆x D

− ∆y TYo γog ∆y D − ∆z TZo γog ∆z D + qno+1 Vijk =

Vijk

∆t

[(φρo So )n+1 − (φρo So )n ]

(A.5)

Eq. (A.5) can be written as Eq. (A.6),

∆To ∆P n+1 − ∆To γog ∆D + qno+1 Vijk =

Vijk

∆t

[(φρo So )n+1 − (φρo So )n ]

(A.6)

Appendix B. Full-implicit discretization Eq. (3.14) of oil phase

The point (i, j, k, n+1) is selected as the discretization item, Eq. (3.18) is expanded as follows,

[ ] ∆Tol ∆Φol = TXoi+ 1 (Φoil +1 − Φoil ) − TXoi− 1 (Φoil − Φoil −1 ) 2 2 ] [ l l l + TYoj+ 1 (Φoj+1 − Φoj ) − TYoj− 1 (Φoj − Φojl −1 ) 2 2 [ ] l l l l + TZok+ 1 (Φok+1 − Φok ) − TZok− 1 (Φok − Φok ) −1

(B.1)

[ ] ∆Tol ∆δ Po = TXoi+ 1 (δ Poi+1 − δ Poi ) − TXoi− 1 (δ Poi − δ Poi−1 ) 2 2 [ ] + TYoj+ 1 (δ Poj+1 − δ Poj ) − TYoj− 1 (δ Poj − δ Poj−1 ) 2 2 [ ] + TZok+ 1 (δ Pok+1 − δ Pok ) − TZok− 1 (δ Pok − δ Pok−1 ) 2 2 { [ l+1 [ ] ]} l+1 l = TXoi+ 1 (Poi+1 − Poi+1 ) − (Poi − Poil ) − TXoi− 1 (Poil+1 − Poil ) − (Poil+−11 − Poil −1 ) 2 2 { [ l+1 ] [ ]} l+1 l l + TYoj+ 1 (Poj+1 − Poj+1 ) − (Poj − Poj ) − TYoj− 1 (Pojl+1 − Pojl ) − (Pojl+−11 − Pojl −1 ) 2 2 { ] ]} [ l+1 [ l+1 l+1 l+1 l l l l + TZok+ 1 (Pok+1 − Pok+1 ) − (Pok − Pok ) − TZok− 1 (Pok − Pok ) − (Pok −1 − Pok−1 )

(B.2)

[ ] ∂ Toi+ 1 ∂ Toi− 1 ∂ To 2 2 l l l l l δ Po ∆Φo = δ P 1 (Φ − Φoi ) − δ P 1 (Φ − Φoi−1 ) ∆ ∂ Po ∂ Poi+ 1 oi+ 2 oi+1 ∂ Poi− 1 oi− 2 oi 2 2 [ ] ∂ Toj+ 1 ∂ Toj− 1 2 2 l l l l + δ P 1 (Φ − Φoj ) − δ P 1 (Φ − Φoj−1 ) ∂ Poj+ 1 oj+ 2 oj+1 ∂ Poj− 1 oj− 2 oj 2 2 [ ] ∂ Tok+ 1 ∂ Tok− 1 2 2 l l l l − Φok ) − + δ P 1 (Φ δ P 1 (Φ − Φok−1 ) ∂ Pok+ 1 ok+ 2 ok+1 ∂ Pok− 1 ok− 2 ok

(B.3)

[ ] ∂ Toi+ 1 ∂ Toi− 1 ∂ To 2 2 l l l l l ∆ δ Sw ∆Φo = δ S 1 (Φ − Φoi ) − δ S 1 (Φ − Φoi−1 ) ∂ Sw ∂ Swi+ 1 wi+ 2 oi+1 ∂ Swi− 1 wi− 2 oi 2 2 [ ] ∂ T w j+ 1 ∂ T w j− 1 2 2 l l l l + δ S 1 (Φ − Φwj ) − δ S 1 (Φ − Φoj−1 ) ∂ Swj+ 1 wj+ 2 wj+1 ∂ Swj− 1 wj− 2 oj 2 2 [ ] ∂ Tok+ 1 ∂ Tok− 1 2 2 l l l l + δS − Φok ) − δS − Φok−1 ) 1 (Φ 1 (Φ ∂ Swk+ 1 wk+ 2 ok+1 ∂ Swk− 1 wk− 2 ok

(B.4)

2

2

2

2

2

2

2

2

Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

22

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Combine the aforementioned equations, the following equation is given,

[ ] ∆Tol ∆Φol + TXoi+ 1 (δ Poi+1 − δ Poi ) − TXoi− 1 (δ Poi − δ Poi−1 ) 2 2 [ ] + TYoj+ 1 (δ Poj+1 − δ Poj ) − TYoj− 1 (δ Poj − δ Poj−1 ) 2 2 [ ] + TZok+ 1 (δ Pok+1 − δ Pok ) − TZok− 1 (δ Pok − δ Pok−1 ) 2 2 [ ] ∂ Toi+ 1 ∂ Toi− 1 2 2 l l l l δ P 1 (Φ δ P 1 (Φ − Φoi−1 ) + − Φoi ) − ∂ Poi+ 1 oi+ 2 oi+1 ∂ Poi− 1 oi− 2 oi 2 2 ] [ ∂ Toj− 1 ∂ Toj+ 1 2 2 l l l l − Φoj ) − + δ P 1 (Φ δ P 1 (Φ − Φoj−1 ) ∂ Poj+ 1 oj+ 2 oj+1 ∂ Poj− 1 oj− 2 oj 2 2 ] [ ∂ Tok− 1 ∂ Tok+ 1 2 2 l l l l + − Φok ) − δ P 1 (Φ δ P 1 (Φ − Φok−1 ) ∂ Pok+ 1 ok+ 2 ok+1 ∂ Pok− 1 ok− 2 ok 2 2 [ ] ∂ Toi+ 1 ∂ Toi− 1 2 2 l l l l + − Φoi ) − δ S 1 (Φ δ S 1 (Φ − Φoi−1 ) ∂ Swi+ 1 wi+ 2 oi+1 ∂ Swi− 1 wi− 2 oi 2 2 [ ] ∂ Twj− 1 ∂ Twj+ 1 2 2 l l l l δ S 1 (Φ δ S 1 (Φ − Φoj−1 ) + − Φwj ) − ∂ Swj+ 1 wj+ 2 wj+1 ∂ Swj− 1 wj− 2 oj 2 2 ] [ ∂ Tok− 1 ∂ Tok+ 1 2 2 l l l l δS − Φok ) − δS − Φok−1 ) + qno+1 VM + 1 (Φ 1 (Φ ∂ Swk+ 1 wk+ 2 ok+1 ∂ Swk− 1 wk− 2 ok 2

VM

=

∆t n+1

+

∆t n+1

VM

(ρ

(B.5)

2

ϕ

l l o So 0 CR

+ϕ

l l So

∂ρol VM ) δ Po − ϕ l ρol δ Sw ∂ Pol ∆t n+1 ] n

(ϕρo So )l − (ϕρo So )

[

References [1] L.K. Abidoye, D.B. Das, Scale dependent dynamic capillary pressure effect for two-phase flow in porous media, Adv. Water Resour. 74 (2014) 212–230. [2] W. Rose, W.A. Bruce, Evaluation of capillary character in petroleum reservoir rock, J. Petrol. Sci. Eng. 1 (05) (1949) 127–142. [3] E.O. Udegbunam, D.T. Numbere, Model for locating fluids’ contracts in petroleum reservoirs, in: SPE Eastern Regional Meeting, Charleston, West Virginia, USA, 8–10 November, 1994. [4] R. Ehrlich, E.L. Etris, D. Brumfield, L.P. Yuan, S.J. CRABTREE, Petrography and reservoir physics III: physical models for permeability and formation factor (1), AAPG Bull. 75 (10) (1991) 1579–1592. [5] R.A. Dawe, Reservoir Physics at the Pore Scale, Commemoration volume of 75, 1990. [6] J. Bear, A. Verruijt, Theory and Applications of Transport in Porous Media, 1987, France. [7] R.H. Brooks, A.T. Corey, Hydraulic properties of porous media, Hydro Paper 3 (27) (1964). [8] M.T. Van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils 1, Soil Sci. Am. J. 44 (5) (1980) 892–898. [9] R. Helmig, Multiphase Flow and Transport Processes in the Subsurface: A Contribution to the Modeling of Hydro Systems, Springer Verlag, 1997. [10] J. Bear, A. Verruijt, Modeling Groundwater Flow and Pollution, Springer Science & Business Media, 2012, p. 2. [11] J. Mikyska, M. Benes, T.H. Illangasekare, Numerical investigation of nonaqueous phase liquid behavior at heterogeneous sand layers using VODA multiphase flow code, J. Porous Media 12 (7) (2009) 685–694. [12] R. Fučík, M. Beneš, J. Mikyška, T.H. Illangasekare, Generalization of the benchmark solution for the two-phase porous-media flow, in: Int. Conf. on Finite-Element Models, MODFLOW, and More: Solving Groundwater Problems, Carlsbad, Czech Republic, 2004. [13] R. Fučík, M. Beneš, J. Mikyška, T.H. Illangasekare, An improved semi-analytical solution for verification of numerical models of two-phase flow in porous media, Vadose Zone J. 6 (2007) 93–104. [14] Y. Muneta, M.I. Mubarak, H.H. Alhasani, K. Arisaka, Formulation of ‘‘capillary force barriers’’ in moderately-oil wet systems and its application to reservoir simulation, SPE Reserv. Eval. Eng. 8 (5) (2005) 388–396. [15] M. Sheffield, Three-phase fluid flow including gravitational, viscous and capillary forces, Soc. Pet. Eng. J. 9 (2) (1969) 255–269. [16] D.A. Weitz, J.P. Stokes, R.C. Ball, et al., Dynamic capillary pressure in porous media: Origin of the viscous-fingering length scale, Phys. Rev. Lett. 59 (26) (1987) 2967–2970. [17] S. Bottero, S.M. Hassanizadeh, P.J. Kleingeld, et al., Nonequilibrium capillarity effects in two — phase flow through porous media at different scales, Water Resour. Res. 47 (10) (2011) 599–609. [18] R.G. Bentsen, Effect of hydrodynamic forces on capillary pressure and relative permeability, Transp. Porous Media 17 (2) (1994) 121–132. [19] B. Agostini, B. Watel, A. Bontemps, B. Thonon, Liquid flow friction factor and heat transfer coefficient in small channels: an experimental investigation, Exp. Therm Fluid Sci. 28 (2) (2004) 97–103. [20] S.M. Hassanizadeh, M.A. Celia, H.K. Dahle, Dynamic effect in the capillary pressure-saturation relationship and its impacts on unsaturated flow, Vadose Zone J. 1 (1) (2002) 38–57. [21] F.J.M. Kalaydjian, Dynamic capillary pressure curve for water/oil displacement in porous media: theory vs. experiment, in: SPE Annual Technical Conference and Exhibition, Washington, D.C. USA, 4–7 October, 1992.

Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.

S. Wang, C. Yu, G. Sang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

23

[22] S.M. Hassanizadeh, W.G. Gray, Thermodynamic basis of capillary pressure in porous media, Water Resour. Res. 29 (10) (1993) 3389–3405. [23] S.M. Hassanizadeh, O. Oung, S. Manthey, Laboratory experiments and simulations on the significance of non-equilibrium effect in the capillary pressure-saturation relationship, in: Unsaturated Soils: Experimental Studies, Springer, Berlin, Heidelberg, 2005, pp. 3–14. [24] G. Camps-Roach, D.M. O’Carroll, T.A. Newson, T. Sakaki, T.H. Illangasekare, Experimental investigation of dynamic effects in capillary pressure: Grain size dependency and upscaling, Water Resour. Res. 46 (8) (2010) 1–13. [25] B. Johannesson, Modeling of a Viscous Fluid Percolating a Porous Material Due to Capillary Forces, Report TVBM. 3095, 2000. [26] R. Fučík, J. Mikyška, Numerical investigation of dynamic capillary pressure in two-phase flow in porous medium, Math. Bohem. 136 (4) (2011) 395–403. [27] V. Joekar-Niasar, S.M. Hassanizadeh, Effect of fluids properties on non-equilibrium capillarity effects: Dynamic pore-network modeling, Int. J. Multiph. Flow. 37 (2) (2011) 198–214. [28] C.J.V. Duijn, M.J.D. Neef, Similarity solution for capillary redistribution of two phases in a porous medium with a single discontinuity, Adv. Water Resour. 21 (6) (1998) 451–461. [29] L. Hou, B.E. Sleep, T.C.G. Kibbey, The influence of unavoidable saturation averaging on the experimental measurement of dynamic capillary effects: A numerical simulation study, Adv. Water Resour. 66 (2) (2014) 43–51. [30] S. Sandoval-Torres, J. Rodríguez-Ramírez, L.L. Méndez-Lagunas, Modeling plain vacuum drying by considering a dynamic capillary pressure, Chem. Biochem. Eng. Q. 25 (3) (2011) 327–334. [31] W.B. Shu, X.U. He-Hua, J.Y. Wan, L.I. Yan-Zhen, Numerical research on dynamic effect of capillary pressure in low permeability reservoirs, Chin. J. Hydrodyn. 29 (2) (2014) 189–196. [32] S. Tian, G. Lei, S. He, L. Yang, Dynamic effect of capillary pressure in low permeability reservoirs, Petrol. Explor. Dev. 39 (3) (2012) 405–411. [33] F. Stauffer, Time dependence of the relations between capillary pressure, water content and conductivity during drainage of porous media, in: IAHR Symposium on Scale Effects in Porous Media, vol. 29, 1978, pp. 335–353. [34] R. Fučík, J. Mikyška, T. Sakaki, Significance of dynamic effect in capillarity during drainage experiments in layered porous media, Vadose Zone J. 9 (3) (2010) 697–708.

Please cite this article as: S. Wang, C. Yu, G. Sang et al., An oil–water two-phase reservoir numerical simulation coupled with dynamic capillary force based on the full-implicit method, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.11.013.