Scholarly article on topic 'An Iterative Implicit Scheme for Nanoparticles Transport with Two-Phase Flow in Porous Media'

An Iterative Implicit Scheme for Nanoparticles Transport with Two-Phase Flow in Porous Media Academic research paper on "Mechanical engineering"

CC BY-NC-ND
0
0
Share paper
Academic journal
Procedia Computer Science
OECD Field of science
Keywords
{Nanoparticles / "Two-phase flow" / "Porous media" / "Oil reservoir" / "Iterative implicit method"}

Abstract of research paper on Mechanical engineering, author of scientific article — Mohamed F. El-Amin, Jisheng Kou, Shuyu Sun, Amgad Salama

Abstract In this paper, we introduce a mathematical model to describe the nanoparticles transport carried by a two-phase flow in a porous medium including gravity, capillary forces and Brownian diffusion. Nonlinear iterative IMPES scheme is used to solve the flow equation, and saturation and pressure are calculated at the current iteration step and then the transport equation is solved implicitly. Therefore, once the nanoparticles concentration is computed, the two equations of volume of the nanoparticles available on the pore surfaces and the volume of the nanoparticles entrapped in pore throats are solved implicitly. The porosity and the permeability variations are updated at each time step after each iteration loop. Numerical example for regular heterogenous permeability is considered. We monitor the changing of the fluid and solid properties due to adding the nanoparticles. Variation of water saturation, water pressure, nanoparticles concentration and porosity are presented graphically.

Academic research paper on topic "An Iterative Implicit Scheme for Nanoparticles Transport with Two-Phase Flow in Porous Media"

^^ Procedia Computer Science

Volume 80, 2016, Pages 1344-1353

ICCS 2016. The International Conference on Computational

Science

An Iterative Implicit Scheme for Nanoparticles Transport with Two-Phase Flow in Porous Media

Mohamed F. El-Amin1, Jisheng Kou2, Shuyu Sun3 , and Amgad Salama3

1 Effat University, Jeddah 21478, Kingdom of Saudi Arabia 2 School of Mathematics and Statistics, Hubei Engineering University, Xiaogan 432000, Hubei, China 3 King Abdullah University of Science and Technology (KAUST), Thuwal 23955-6900, Kingdom of

Saudi Arabia momousa@effatuniversity.edu.sa

Abstract

In this paper, we introduce a mathematical model to describe the nanoparticles transport carried by a two-phase flow in a porous medium including gravity, capillary forces and Brownian diffusion. Nonlinear iterative IMPES scheme is used to solve the flow equation, and saturation and pressure are calculated at the current iteration step and then the transport equation is solved implicitly. Therefore, once the nanoparticles concentration is computed, the two equations of volume of the nanoparticles available on the pore surfaces and the volume of the nanoparticles entrapped in pore throats are solved implicitly. The porosity and the permeability variations are updated at each time step after each iteration loop. Numerical example for regular heterogenous permeability is considered. We monitor the changing of the fluid and solid properties due to adding the nanoparticles. Variation of water saturation, water pressure, nanoparticles concentration and porosity are presented graphically.

Keywords: Nanoparticles, Two-phase flow, Porous media, Oil reservoir, Iterative implicit method

1 Introduction

Recently, applications of nanoparticles (1-100 nm) have been reported in petroleum industry such as oil and gas exploration and production that become a promising field of research. In general, the sizes of PN are in the range of 10 — 500 nm, while pore radii of a porous medium (sandstone) are from 6 to 6.3 x 104 nm. However, if a particle larger than a pore throat may block at the pore throat during nanoparticles transport with flow in the porous medium. There are certain types of nanoparticles, such as polysilicon nanoparticles (PN), can be used in oilfields to enhance water injection by changing the wettability through their adsorption on porous walls. The PN may be classified based on wettability of their surfaces into two types,namely:

* Corresponding Author

1344 Selection and peer-review under responsibility of the Scientific Programme Committee of ICCS 2016

© The Authors. Published by Elsevier B.V.

doi:10.1016/j.procs.2016.05.423

(1) Lipophobic and hydrophilic PN and exists in water phase only.

(2) Hydrophobic and lipophilic PN and exists only in the oil phase.

In [18, 17] authors have founded a mathematical model of nanoparticles transport in two-phase flow in porous media based on the formulation of fine particles transport in two-phase flow in porous media provided in Refs. [21]. El-Amin et al. [9, 8, 11] introduced a model to simulate nanoparticles transport in two-phase flow in porous media. In [12, 10], authors extended the model to include a negative capillary pressure and mixed relative permeabilities correlations to fit with the mixed-wet system. The model of flow and transport of nanoparticles in porous media consists mainly of five PDEs, namely, pressure, saturation, nanoparticles concentration, volume of the nanoparticles available on the pore surfaces , and finally the volume of the nanoparticles entrapped in pore throats. Moreover, the model includes variations of both porosity and permeability due to nanoparticles precipitation on the pores walls of the medium.

The model of two-phase fluid flow in porous media is a coupled system of nonlinear time-dependent partial differential equations. Two different types of time discretization schemes are often used to solve this coupled system. The first one is the fully implicit scheme [2, 6, 7] that implicitly treats with all terms including capillary pressure. This scheme results in a system of nonlinear equations and has unconditional stability and maintains the inherent coupling of two-phase flow model. The second scheme is the IMplicit-EXplicit (IMEX) [1, 3, 13, 16] which generally treats the linear terms implicitly and evaluates the others explicitly, and consequently. This scheme is conditionally stable, however, it has advantage that is to eliminate the nonlinearity of original equations. The IMplicit Pressure Explicit Saturation (IMPES) approach is viewed as an IMEX method, solves the pressure equation implicitly and updates the saturation explicitly. The IMPES method is conditionally stable, and hence it must take very small time step size, especially for highly heterogeneous permeable media where the capillary pressure affects substantially on the path of fluid flow. The instability of the IMPES method [5] results from the decoupling between the pressure equation and the saturation equation as well as the explicit treatment of the capillary pressure. The IMPES for two-phase flow has been improved in several versions (e.g. [4]). As an iterative method, the computational cost and memory required by iterative IMPES method is smaller than the fully coupled approach at each iterative step, which is more pronounced for very large size computational problems. The main disadvantage of iterative IMPES method is the decoupling of pressure and saturation equations, which results from the explicit treatment for capillary pressure. A linear approximation of capillary function is introduced to couple the implicit saturation equation into pressure equation [19]. Kou and Sun [20] presented an iterative version of their previous scheme proposed in [19]. Unlike iterative IMPES, capillary pressure is not computed by the saturations at the previous iteration, but the linear approximation of capillary function at the current iteration is used, which is constructed by the saturations at the current and previous iterations.

In this work, we use the iterative IMPES scheme introduced in [20, 19] to solve the flow equation of the model of nanopaticles transport in porous media. Then, we used the saturation and pressure calculated at the current iteration step to calculate the transport equation implicitly. Therefore, once we compute the nanoparticles concentration, the two equations of volume of the nanoparticles available on the pore surfaces and the volume of the nanoparticles entrapped in pore throats are solved implicitly. In our scheme, we update the porosity and the permeability variations at each time step after each iteration loop.

2 Modeling and Mathematical Formulation

2.1 Two-phase flow Model

In this section, a mathematical model is developed to describe the nanoparticles transport carried by two-phase flow in porous media. Let us consider two-phase immiscible incompressible flow in a heterogeneous porous medium domain governed by the Darcys law and the equations of mass conservation for each phase as,

d (S) , ^ m

——--+ V- ua = qa, a = w,n. (1)

Ua =--—K (Vpa + pag) , a = W,n. (2)

where Sa is the saturation, ua is the velocity of the phase a. w stands for the wetting phase, and n stands for the nonwetting phase. $ is the porosity of the medium, qa is the external mass flow rate. K is the absolute permeability tensor is chosen as K = kI, where I is the identity matrix and k is a positive real number. kra is the relative permeability, pa is the density, and pa is the pressure of the phase a, g = (0, -g)T is the gravitational acceleration. pa is the viscosity and ka = kraK is the effective permeability. The fluid saturations for the wetting and non-wetting are interrelated by,

Sw + S„ = 1. (3)

Now, we describe the governing equations used in [15] and [20] as,

V- (un + uc) = -V • AtKV$w - V • A„KV$C = qw + qn. (4)

- qw = -V- (fw ua) = —V • AtKV$w. (5)

where fw = Aw/At is the flow fraction, Aa = kra/pa is the mobility, $w = pw + pwg is the water pressure potential, and = pc + (pn - pw) g is the capillary pressure potential. The total velocity ut = uw + un = ua + uc is defined as the sum of the two velocity variables ua = -AtKV$w and uc = -AnKV$c. The wetting-phase velocity may be expressed by, uw = fwua. The two-phase capillary pressure can be expressed by, pc (Sw) = pn - pw.

2.2 Nanoparticles Transport Model

Ju and Fan [18] reported that there are two types of polysillicon nanoparticles (PN) can be used in oil fields to improve oil recovery and enhance water injection, respectively. The polysillicon nanoparticles are classified based on wettability of the surface of the PN. The first type is called lipophobic and hydrophilic polysillicon nanoparticles (LHPN) and exists in water phase only, while the second type is called hydrophobic and lipophilic polysillicon nanoparticles (HLPN) and exists in the oil phase only. The sizes of PN are in the range of 10 to 500 nm, therefore, Brownian diffusion is considered. Assuming that we have a number m of size interval of the nanoparticles in water phase, the transport equation for each size interval i of the nanoparticles in the water/oil phase can be written as,

d ($SatCi,a) + ua • VC\a = V • - Ri,a + Qi,a. (6)

where i = 1, 2,,m. Cia is the volume concentrations of nanoparticles in size interval i in the phase a.

Nanoparticles that have sizes smaller than microns have strong Brownian motion, and can bring nanoparticles very close to the pore wall. Therefore, nanoparticle may retention to decrease as flow velocity increases. The diffusion coefficient Di of the nanoparticle can be calculated using the Stokes-Einstein equation,

i 3njdp, ( )

where kB is the Boltzmann constant, T absolute temperature, dp is the particle diameter, and j is the fluid viscosity. In Eq. (7), particle diffusion constant is inversely proportional to particle diameter dp, which increases with decreasing particle size. For instance, for water at 20oC the Brownian diffusivity for a 1-jm-diameter particle is 4.3 x 10~9cm2 , which is small compared to solute diffusion but potentially significant over the small distances within pore spaces. Qia is the rate of change of particle volume belonging to a source/sink term. Ria is the net rate of loss of nanoparticles in size interval i in the phase a. The net rate of loss of nanoparticles may be written as [18, 17, 21],

d (54>).

Ria = (8)

where (5^)i a = v1iia + v2,i,a is the porosity change due to release or retention of nanoparticles of interval i in the phase a. v1,i,a is the volume of the nanoparticles of interval size i in contact with the phase a available on the pore surfaces per unit bulk volume of the porous medium. v2,i,a is the volume of the nanoparticles of interval size i entrapped in pore throats from the phase a per unit bulk volume of porous medium due to plugging and bridging. Also, v1ia and v2,i,a may be defined in terms of the mass of particles per unit fluid volume deposited at the pore bodies a1,i,a and pore throats a2,i, a of the porous medium as,

a1 ,ia &2,i,a /„N

v1 ,i, a = -, v2 ,i, a = -. (9)

where pb is the density of particulate suspensions.

At the critical velocity of the surface deposition only particle retention occurs while above it retention and entrainment of the nanoparticles take place simultaneously (Gruesbeck and Collins [14]). A modified Gruesbeck and Collins's model for the surface deposition is expressed by [18],

I \ld,i,a |ua| Ci,a, |ua| < Un

Yd,i,a |ua | Ci,a Ye,i,av1,i,a |ua uac| , |ua| ^ uo

where Yd}i,a is the rate coefficients for surface retention of the nanoparticles in interval i in the phase a. Ye,i,a is the rate coefficients for entrainment of the nanoparticles in interval i in the phase a. uac is the critical velocity for the phase a. Similarly, the rate of entrapment of the nanoparticles in interval i in the phase a is,

dv2, i, a

where jpt^, a is the pore throat blocking constants.

= Ypt,i,a |ua| Ci,a, (11)

2.3 Porosity and Permeability Variations

Porosity may be changed because nanoparticles deposition on the pore surfaces or blocking of pore throats. The porosity variation may be by [18, 21],

$ = $0 - £ (^)i a (12)

where $0 is the initial porosity. Also, the permeability variation due to nanoparticles deposition on the pore surfaces or blocking of pore throats may be expressed as [17],

K = Ko

(1 - f ) kf + f f

where K0 is the initial permeability, kf is constant for fluid seepage allowed by the plugged pores. The flow efficiency factor expressing the fraction of unplugged pores available for flow is given by,

f =1 - E Yfi(E A (14)

i \ a /

Yf,i is the coefficient of flow efficiency for particles i. The value of the exponent l has range from 2.5 to 3.5. For the nanoparticles transport carried by fluid stream in the porous media, deposition on pore surfaces and blockage in pore throats may occur. The retained particles on pore surfaces may desorb for hydrodynamic forces, and then possibly adsorb on other sites of the pore bodies or get entrapped at other pore throats.

2.4 Initial and Boundary Conditions

The initial variables at the beginning of the flow displacing process in the computational domain Q are defined by,

Sw = SW, Ci a a = 0, vi, i a a =0, V2, i, a = 0 in Q at t = 0. (15)

The boundary dQ of the computational domain Q is subjected to both Dirichlet and Neumann conditions such that dQ = rD U rN and rD n =0, where r^ is the Dirichlet boundary and r^ is the Neumann boundary. The boundary conditions considered in this study are summarized as follow,

pw (or pn)= pD on rD, (16)

ut • n = qN, Sw (or Sn) = SN, Ci,w = C0w, on Viia = 0, Vi^a = 0 , (17)

where n is the outward unit normal vector to dQ, pD is the pressure on rD and qN the imposed inflow rate on rN, respectively.

3 Iterative Implicit Scheme

In this study, we consider only one interval size in the wetting phase. So, for example we drop the subscript from Cj, , w, v1,i , a and v2,i., a to become C, v1 and v2, respectively. Similarly, we will drop the subscript i, w, a from all symbols including constants. Define the time step

length Atn = tn+1 _ tn, the total time interval [0,T] may be divided into NT time steps as 0 = t° < t1 < ■ ■ ■ < tNT = T. The current time step is represented by the superscript n + 1. Also, the iteration loop consists of a number Nj iterations of each time step. The current iteration step is denoted by k +1. The iterative backward Euler time discretization is used for the equations of pressure, saturation, concentration and the two volumes to obtain,

-V ■ A, (S^+1,fc) K (vn+1,k, vn+1,k) V^l+W Xn (SW+1,k) K (vn+1,k, vn+1,k) V$c (Sl+1,k+1) = ql+1 + qn+\

n+1,k n+1,k

n+1,k+1

w nn+1

= —V ■ A, (SW+1,k)

(vvn+1,kV$W+1,k+1,

n+1,k n+1,k >1 , ,v2 ,

Sn+1, k + 1 cn+1, k+1 _ sn Cn

V ■ {^[v^,vn+1k) SW+1k+1 DVCn+1,k+1} + R (sw+1,k+1, $W+1,k+1, Cn+1,k, vn+1k) + QW+1,

+ V ■ (uW+1,k+1Cn+1„k+^ =

n+1,k+1

v1 ' _ v

Yd |uW+1,k+1| Cn+1,k+1,

Yd |uW+1,k+11 Cn+1,k+1 _ Ye |uW+1,k+1 _ ucl v1

n+1,k+1

(c1) , (C2)

n+1,k+1 vn , _v

,n+1,k + 11 s~in+1,k + 1 C

From now upon we will refer to the condition |uW+1,k+1| < uc by (c1) and to the condition |uW+1,k+1| > uc by (c2). The superscripts k and k + 1 represent the iterative steps within the current time step n +1. For each iteration, the variables AW,An,At and fW is calculated using the saturation from the previous iteration. The pressure equation is solved firstly to obtain the wetting-phase pressure at the current iteration and then the Darcys velocity can be calculated. Therefore, the saturation at the current iteration is computed explicitly in the current iteration. Then, the concentration and values are computed implicitly at the current time step. Finally, the permeability, porosity, and other parameters such as AW,An,At and fW are updated. This procedure is repeated until the convergence criterion of errors has been satisfied. In this iterative scheme the capillary potential $c is linearized as follows [20],

$c {sW+1,k+1) = {sW+1,k) + (SW+1k) [srW+1,k+1 _ sW+1,k]

where is derivative of $c. The changes of saturation in a time step are often very small, and hence the linear approximation is reasonable. On the other hand, we use the relaxation

approach to control the convergence of nonlinear iterative solvers. Therefore, the iterative scheme of pressure and saturation equations may be rewritten as,

-V • At (Sn+1'k) Kn+1'fcV$w+1'fc+1 -V^An (Sn+1'k) Kn+1'k V$ c (Sn+1'k) = Qn+1, (24) $c (Zw+1'k+1) = (Sw+1'k) + (Sw+1'k) [sw+1'k+1 - Sw+1'k], (25)

on+1,k+1 _ on

$n+1'k w —w + V^At (Sw+1'k) Kn+1-kV$w+1'k+1 = qw+1, (26)

where the relaxation equation is,

Sn+1'k+1 _ Sn+1'k + ^ ^Sn+1'k+1 _ Sn+1'^ (27)

where = (0,1] is a relaxation factor for saturation. In a similar manner we may write implicit iterative scheme of the equation of concentration as follow,

$n+1'k Sw_C_°wC + v^ ^un+1'k+1cn+1'k + 1^ = v^ $n+1'kSn+1'k+1

Sn+1'k + 1cn+1'k+1 _ Sn Cn

DVCn+1'k+1) + R Sw+1'k+\$w+1'k+1, Cn+1'k, vn+1'k^j + Qn+1. (28)

In the above equation, we use v1 from the previous iteration step in the term R. Once Cn+1'k+1 is obtained one can get v!+1'k+1 as follows,

n+1'k+1 n 1 - v1

Yd\uw+1'k+1\ Cn+1'k+1, (a)

Yd\uw+1'k+1\Cn+1'k+1 - Ye\ul+1'k+1 - uc\Zn+1'k+1, (C2) Finally, the equation of v2 may be written as,

n+1'k+1 2

= Ypt \uw+1'k+1\ Cn+1'k+1, (30)

4 Numerical Experiments

In this section, we test some examples to show the performance of the presented scheme. Before presenting the numerical examples, let us introduce the necessary physical parameters used in the computations.

In this study, we consider the following capillary pressure formula,

pc = -Bc log(S),

and the normalized wetting phase saturation are correlated by,

S = SS - SZ , 0 < S < 1,

1 Snr Swr

where Bc is the capillary pressure parameter, Swr is the irreducible (minimal) water (wetting phase) saturation, and Snr is the residual (minimal) oil (nonwetting phase) saturation after water flooding.

Also, the expressions of the relation between the relative permeabilities and the normalized wetting phase saturation S is given as,

k = k0 Sa k = k0 (1 — S)b

krw krws , krn krn (1 s) '

where a and b are positive real numbers, k°w = krw (S = 1) is the endpoint relative permeability to the wetting phase, and k°n = krn (S = 0) is the endpoint relative permeability to the nonwetting phase. The capillary pressure function and relative permeabilities are chosen to be zero for the residual saturations of water and oil; that is, S = Sw. In computation, we take the minimum of saturation as Sw,min = 10-4. Moreover, since the relaxation factor has some effects on the stability and efficiency of iterative methods. The choice strategy of relaxation factor [20], is applied in these calculations. The relaxation factor 0S is computed where we take ISW+1,0 _ SW+1~11 = 1. The iteration loop continues until |SW+1,k+1 _ SW+1,k| < Sw^in. In our tests, we use 2-norm for vectors and matrices.

Moreover, we consider one size interval of nanoparticles suspension in the water phase. The following parameters values are used in the calculations, Yd = 16 m-1, Ypt = 1.28 m-1, Ye = 30 m-1, Yd = 16 m-1, uc = 4.6 x 10-6 m/s, and D = 5.6 x 10-8 m2/s. The nanoparticles diameter is taken as 40 nm and concentration C0 = 0.1. Swr = Snr = 0.001, = 0.3, kf = 0.6, Yf = 0.01. The viscosities of water and oil are 1 cP and 0.45 cP, respectively. The injection rate is 1.5 PV/year. The relative permeabilities are quadratic, krw0 = kro0 = 1, a = b = 2, and the capillary pressure parameter is Bc = 50 bar. The domain dimension is taken as 0.3 m x 0.2 m. The computational domain is divided into 2400 uniform rectangles. The choice of relaxation factor given in Ref. [20] is applied for the iterative method, and the three parameters are taken as 9s,min = 0.1,9s,max = 0.9 and p = 0.2. The time step used in this example is

0.0075 day. In the given example, the tested medium consists of two subdomains with different configurations for the distribution of permeability as shown in the lower left graph of Figure

1. Distributions of water saturation (upper left), water pressure (upper right), nanoparticles concentration (middle left), deposited nanoparticles concentration (middle right) and porosity reduction (lower right) are shown in the Figure 1 too.

References

[1] U. Ascher, S. J. Ruuth, and B. R. Wetton. Implicit-explicit methods for time-dependent partial differential equations. SIAM J. Numer. Anal., 32:797-823, 1995.

[2] K. Aziz and A. Settari. Petroleum reservoir simulation. Science Pub., London, 1979.

[3] S. Boscarino. Error analysis of imex runge-kutta methods derived from differential-algebraic systems. SIAM J. Numer. Anal., 45:1600-1621, 2007.

[4] Z. Chen, G. Huan, and Y. Ma. Computational methods for multiphase flows in porous media. SIAM Comp. Sci. Eng., Philadelphia, 2006.

[5] K. H. Coats. Impes stability: selection of stable time steps. SPE-84924, SPE Reservoir-Simulation Symposium, Houston, TX., 2001.

[6] D. A. Collins, L. X. Nghiem, Y. K. Li, and J. E. Grabenstetter. An efficient approach to adaptive implicit compositional simulation with an equation of state. SPE Reservoir Engineering, 7:259264, 1992.

[7] C. N. Dawson, H. Kle, M. M. Wheeler, and C. S. Woodward. A parallel, implicit, cell-centered method for two-phase flow with a preconditioned newton- krylov solver. Computational Geo-sciences, 1:215-249, 1997.

Figure 1: Variation of saturation (upper left), pressure (upper right), concentration (middle left), deposited concentration (middle right), permeability (lower left), and porosity (lower right).

[8] M. F. El-Amin, A. Salama, and S. Sun. Modeling and simulation of nanoparticle transport in multiphase flows in porous media: Co2 sequestration. SPE-163089, Mathematical Methods in Fluid Dynamics and Simulation of Giant Oil and Gas Reservoirs, 2012.

[9] M. F. El-Amin, A. Salama, and S. Sun. Modeling and simulation of nanoparticles transport in a two-phase flow in porous media. SPE-154972, SPE International Oilfield Nanotechnology Conference and Exhibition, Noordwijk, The Netherlands, 2012.

[10] M. F. El-Amin, A. Salama, and S. Sun. Numerical and dimensional analysis of nanoparticles transport with two-phase flow in porous media. Journal of Petroleum Science and Engineering, 128:53 - 64, 2015.

[11] M. F. El-Amin, A. Salama, and S. Sun. Numerical and dimensional analysis of nanoparticles transport with two-phase flow in porous media. Journal of Petroleum Science and Engineering, 128:53-64, 2015.

[12] M. F. El-Amin, S. Sun, and A. Salama. Enhanced oil recovery by nanoparticles injection: modeling and simulation. SPE-164333, SPE Middle East Oil and Gas Show and Exhibition held in Manama,, Bahrain, 10-13 March, 2013.

[13] J. Frank, W. Hundsdorfer, and J. G. Verwer. On the stability of implicit-explicit linear multi-step methods. Appl. Numer. Math., 25:193-205, 1997.

[14] C. Gruesbeck and R. E. Collins. Entrainment and deposition of fines particles in porous media. SPE J, 24:847-855, 1982.

[15] H. Hoteit and A. Firoozabadi. Numerical modeling of two-phase flow in heterogeneous permeable media with different capillarity pressures. Advances in Water Resources, 31:56 - 73, 2008.

[16] W. Hundsdorfer and S. J. Ruuth. Imex extensions of linear multistep methods with general monotonicity and boundedness properties. J. Comput. Phys., 225:2016-2042, 2007.

[17] B. Ju and T. Fan. A study of wettability and permeability change caused by adsorption of nanometer structured polysilicon on the surface of porous media. SPE-77938, SPE Asia Pacific Oil and Gas Conference and Exhibition, Melbourne, Australia, 2002.

[18] B. Ju and T. Fan. Experimental study and mathematical model of nanoparticle transport in porous media. Powder Technology, 192:195 - 202, 2009.

[19] J. Kou and S. Sun. A new treatment of capillarity to improve the stability of impes two-phase flow formulation. Comp. Fluids, 39:1923-1931, 2009.

[20] J. Kou and S. Sun. On iterative impes formulation for two phase flow with capillarity in heterogeneous porous media. International Journal of Numerical Analysis and Modeling. Series B, 1(1):20-40, 2010.

[21] Liu X. H. and F. Civian. Characterization and prediction of formation damage in two-phase flow systems. SPE-25429, Production Operations Symposium,, Oklahoma City, OK, U.S.A., 1993.