International Journal of Multiphase Flow xxx (2010) xxx–xxx
Contents lists available at ScienceDirect
International Journal of Multiphase Flow
j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / i j m u l fl o w
Effect of fluids properties on non-equilibrium capillarity effects: Dynamic
pore-network modeling
Vahid Joekar-Niasar ⇑, S. Majid Hassanizadeh
Department of Earth Sciences, Utrecht University, PO Box 80021, 3508 TA Utrecht, The Netherlands
a r t i c l e
i n f o
Article history:
Received 22 January 2010
Received in revised form 14 September
2010
Accepted 16 September 2010
Available online xxxx
Keywords:
Porous media
Non-equilibrium capillarity effect
Two-phase flow
Dynamic pore-network model
Drainage
Imbibition
a b s t r a c t
We have developed a Dynamic Pore-network model for Simulating Two-phase flow in porous media
(DYPOSIT). The model is applicable to both drainage and imbibition processes. Employing improved
numerical and geometrical features in the model facilitate a physically-based pore-scale simulator. This
computational tool is employed to perform several numerical experiments (primary and main drainage,
main imbibition) to investigate the current capillarity theory. Traditional two-phase flow formulations
state that the pressure difference between the two phase is equal to the capillary pressure, which is
assumed to be a function of saturation only. Many theoretical and experimental studies have shown that
this assumption is invalid and the pressure difference between the two fluids is not only equal to the capillary pressure but is also related to the variation of saturation with time in the domain; this is referred to
as the non-equilibrium capillarity effect. To date, non-equilibrium capillarity effect has been investigated
mainly under drainage. In this study, we analyze the non-equilibrium capillarity theory under drainage
and imbibition as a function of saturation, viscosity ratio, and effective viscosity. Other aspects of the
dynamics of two-phase flow such as trapping and saturation profile are also studied.
Ó 2010 Elsevier Ltd. All rights reserved.
area and they proposed the following equation for capillary
pressure:
1. Introduction
1.1. Theory of non-equilibrium capillarity effect
Classical equations for multiphase flow in porous media are
based a central equation – capillary pressure–saturation relationship. This equation, which is written based on the thermodynamically-equilibrium assumption, is commonly written as:
Pn —Pw ¼ Pc ðSw Þ
ð1Þ
in which, Pc is the capillary pressure, Sw is the saturation of the wetting fluid, and Pn and Pw are the nonwetting and wetting fluid pressures, respectively.
In fact, there are two major assumptions in this equation: (a)
fluid distribution and consequently capillary pressure is a function
of wetting fluid saturation only and (b) fluids pressure difference is
equal to capillary pressure (at all times and under all conditions).
Regarding the first assumption, it is known that there is hysteresis in Pc–Sw curves obtained for drainage and imbibition processes. Using a thermodynamic approach, Hassanizadeh and Gray
(1990, 1993b) have suggested that the non-uniqueness in the Pc–
Sw relationship is indeed due to the absence of specific interfacial
⇑ Corresponding author.
E-mail address: joekar@geo.uu.nl (V. Joekar-Niasar).
Pc ¼ P c ðSw ; anw Þ
ð2Þ
A number of computational and experimental works have shown
that Pc–Sw–anw surfaces for drainage and imbibition more or less
coincide under a wide range of drainage and imbibition history. This
means that inclusion of anw leads to the removal, or significant
reduction, of hysteresis in Pc–Sw relationship (e.g. Reeves and Celia,
1996; Held and Celia, 2001; Cheng et al., 2004, 2007; Joekar-Niasar
et al., 2008, 2009, 2010b; Porter et al., 2009).
Regarding the second assumption underlying Eq. (1), it is now
an established fact that Pn–Pw is equal to capillary pressure but
only under equilibrium conditions (see Hassanizadeh et al.
(2002) for an extended review of experimental evidences). According to Entov (1980) capillary pressure–saturation relationship is
not unique and, even though it is obtained under equilibrium conditions, it is a function of the history of fluids movements. In fact, it
depends not only on the volume fraction of each phase, but also on
their microscale distribution and change of saturation with time.
Non-equilibrium effects in capillary pressure can be of significant
importance in industrial porous media, such as paper pulp drying
process, where gradients of fluids pressure and fluids velocities
are large (Lewalle et al., 1994). For non-equilibrium conditions,
the following equation for the difference in fluids pressure has
0301-9322/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved.
doi:10.1016/j.ijmultiphaseflow.2010.09.007
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
2
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
been suggested (Stauffer, 1978; Kalaydjian and Marle, 1987;
Hassanizadeh and Gray, 1990):
Pn —Pw ¼ Pc s
@Sw
@t
ð3Þ
where s, a non-equilibrium capillarity coefficient, is a material
property that may still be a function of saturation, and fluid–fluid
properties (e.g. interfacial area, contact line). Potential dependencies of s are discussed in the next paragraph. In this work, we focus
on investigating the validity of Eq. (3) by means of numerical experiments and determining the dependencies of the coefficient s.
1.2. Non-equilibrium capillarity coefficient
Eq. (3) has been the subject of many studies in recent years,
computationally using Darcy-scale models (see e.g. Manthey
et al., 2005; Das et al., 2006), and pore-scale models (see e.g. Dahle
et al., 2005; Gielen et al., 2005; Joekar-Niasar et al., 2010a) as well
as experimentally see e.g. Hassanizadeh et al., 2004; Oung et al.,
2005; O’Carroll et al., 2005; Bottero and Hassanizadeh, 2006;
Bottero, 2009; Berentsen and Hassanizadeh, 2006; Camps-Roach
et al., 2010; O’Carroll et al., 2010). Hassanizadeh et al. (2002)
reviewed extensively the experimental works in which nonequilibrium effects have been observed.
Up to now, most studies have studied the variation of s with
saturation. However, based on pore-scale insights Hassanizadeh
et al. (2002) conjectured that s may be related to phase trapping,
capillary blockage (and consequently interfacial area and contact
lines), and contact angle. Furthermore, in a recent study, O’Carroll
et al. (2010) stated that s may be a function of fluid/fluid/solid contact line friction. Dahle et al. (2005) developed a bundle-of-tube
model to investigate the variation of s with variance of radii distribution and with saturation under drainage process. They found
that s increases with decrease of wetting fluid saturation and with
increase of variance. In another study, Joekar-Niasar et al. (2010a)
investigated variation of s with viscosity ratio (M) and saturation
during drainage using a dynamic pore-network model. Viscosity
ratio is defined as M = ln/lw, where ln and lw are the viscosities
of nonwetting and wetting fluid, respectively. Joekar-Niasar et al.
(2010a) found a similar trend for variation of s with saturation
as Dahle et al. (2005) did. In addition, they found that larger values
of viscosity ratio results in larger values of s. In another study,
Manthey et al. (2005) studied variation of s with permeability, heterogeneity and entry capillary pressure during drainage. They employed a continuum-scale two-phase simulator in which local
capillary pressure was defined based on Brooks-Corey relation.
Based on volumetric phase averaging, they determined the magnitude and behavior of s.
Although some aspects of non-equilibrium capillarity effect has
been studied, there is still much room for further studies. There is
still no study on the combined effect of saturation, effective viscosity, and viscosity ratio on the variation of s or on the presence of
hysteresis on s–Sw relationship.
1.3. Dynamic pore-network modeling
We employ pore-network modeling as an upscaling technique
from pore scale to macroscale. This approach is computationally
not too expensive, compared to other pore-scale simulators such
as Lattice Boltzmann or Smoothed Particle Hydrodynamics
(see Joekar-Niasar et al., 2010a, for a comparison between
pore-network modeling and Lattice Boltzmann).
Dynamic pore-network models can simulate transient behavior
of flow with time for various capillary numbers and viscosity ratios. Capillary number (Ca) is traditionally defined as the ratio of
qinv
viscous forces of the invading fluid to capillary forces linrvnw
.
The first dynamic pore-network model reported in the literature
was developed by Koplik and Lasseter (1985) to simulate imbibition process in a two-dimensional pore network. Later, several dynamic pore-network models were developed for various
applications, such as simulating pressure field evolution with time
in a two-phase drainage process (see e.g. Aker et al., 1998a; Aker
et al., 1998b; Van der Marck et al., 1997; Dahle and Celia, 1999),
effect of capillary number on residual water saturation as well as
fractional flow behavior during drainage (see e.g. Singh and
Mohanty, 2003; Al-Gharbi and Blunt, 2005), effect of flow regime
on residual nonwetting fluid in two-phase imbibition processes
(see e.g. Nguyen et al., 2006; Mogensen and Stenby, 1998; Koplik
and Lasseter, 1985; Hughes and Blunt, 2000; Thompson, 2002),
evaporation (e.g. Prat, 2002), three-phase flow (see e.g. Pereira
et al., 1996) ganglia movement see e.g. Constantinides and Payatakes, 1996; Dias and Payatakes, 1986a,b) and its contribution to
relative permeability (e.g. Avraam and Payatakes, 1995a), and
fluid–fluid interface velocity under drainage (Nordhaug et al., 2003).
Strong nonlinearity at the pore-scale causes severe numerical
stability problems in dynamic pore-network models, such that
simulation for unfavorable viscosity ratios or for capillary-dominated flow have often proven troublesome. For instance,
Al-Gharbi and Blunt (2005) could not obtain a full consistency
between fluid occupancy at equilibrium resulted from quasi-static simulations and dynamic simulations for the same boundary
conditions.
In order to avoid these problems and also to be able to simulate imbibition as well as drainage under a wide range of viscosity ratios and different flow regimes, a new computational
algorithm and a robust pore-network model was developed in
Joekar-Niasar et al. (2010a). This new model has the following
main features:
(a) In contrast to all other dynamic pore-network models
(except for Thompson, 2002), we assign two separate pressure fields to the two fluids within pore bodies or pore
throats, and thus assign local capillary pressure to pore
bodies.
(b) Pore-scale flow mechanisms such as snap-off and cooperative pore filling (Lenormand and Zarcone, 1983, 1984) during both drainage and imbibition are included.
(c) A semi-implicit approach is used for the saturation update in
order to obtain numerical stability in simulations even for
capillary-dominated and viscous fingering regimes. Thus,
the resulting set of equations for fluid pressures contain both
advection-type terms (corresponding to viscous forces) and
diffusion-type terms (corresponding to capillary forces).
(d) To mimic the geometry of pore space in granular media, pore
bodies and pore throats are represented as truncated octahedrons and parallelepipeds, respectively. Local capillary pressure is determined from local interface curvature, which is
related to local pore body saturation.
Details of the computational algorithm employed in our
dynamic pore-network model are given in Joekar-Niasar et al.
(2010a). They developed a pore-network model with cubic pore
bodies and parallelepiped pore throats to simulate drainage
processes under different pressure gradients imposed at boundaries. They analyzed effect of viscosity ratio on non-equilibrium
capillarity coefficient for different saturations. In addition, nonequilibrium effects on specific interfacial area were studied. Their
results were only for primary drainage.
1.4. Objectives
The main objectives of this paper are as follows:
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
Developing a dynamic pore-network model (called DYPOSIT)
for simulating two-phase flow under both drainage and imbibition conditions. The model employs improved numerical algorithms for computing pressure fields and updating the
saturation. Using these algorithms, the transient behavior of
capillary pressure under different flow regimes and viscosity
ratios for drainage and imbibition processes are simulated; such
studies under imbibition have not been reported in the literature. To mimic larger pore spaces among the solid grains in a
granular medium, the shape of pore bodies has been selected
as truncated octahedrons. This is a modification of the geometry
compared with Joekar-Niasar et al. (2010a), in which pore
bodies were represented as cubes.
Studying the evolution of average fluids pressure difference (Pn–
Pw) with time and saturation for M = 0.1, 1, 10 for primary and
main drainage as well as main imbibition processes. As a result,
variation of s with viscosity ratio (and consequently effective
viscosity) and saturation under drainage and imbibition has
been analyzed.
Analysis of the validity of Eq. (3) during imbibition as well as
drainage. Up to now, behavior of s has been studied only during
primary drainage. We are interested in analyzing existence of a
unique s under conditions where the local equilibrium at interfaces is still assumed (so, for example, dynamics of contact
angle is neglected).
3
(a)
(b)
2. Model description
(c)
2.1. Model features
2.1.1. Structure and geometry
Coordination number is a topological property of the network,
which is the number of pore throats connected to a given pore
body. Although there are many irregular quasi-static pore networks (e.g. Bakke and Øren, 1997; Joekar-Niasar et al., 2009,
2010b; Raoof and Hassanizadeh, 2010), the majority of dynamic
pore networks have a regular structure. We have also developed
a regular three-dimensional lattice network with fixed coordination number of six. Pore bodies and pore throats are presented
by ‘‘truncated octahedron” and ‘‘parallelepiped”, respectively. This
allows simultaneous existence of both fluids in a single pore element. The octahedron pore bodies can be unequally truncated
since they are connected to pore throats of different sizes, as
shown in Fig. 1a. Truncated parts of the octahedron have the shape
of square pyramids with base width of 2rij (where rij is equal to the
radius of inscribed circle of the pore throat ij) as shown in Fig. 1a.
The radius of the inscribed sphere of a pore body, Ri, and the radius
0
of the inscribed circle of the diagonal cross sectionpof
ffiffi pore body, Ri ,
are shown in Fig. 1b. It should be noted that Ri ¼ 36 R0i . A cross section through the vertices of a pore body and the pore throat connecting them is shown in Fig. 1c.
The size distribution of pore bodies is specified by a truncated
log-normal distribution, with no spatial correlation, expressed by:
f ðRi ; rnd Þ ¼ qffiffiffiffiffiffiffiffiffiffiffi
"
2 #
R
pffiffiffi
ln i
2 exp 12 r Rm
nd
"
!
Rmax
Rmin
2r
2r
!#
ln R
ln R
m
m
erf pffiffiffiffiffiffiffi
pr2nd Ri erf pffiffiffiffiffiffiffi
2
2
nd
nd
ð4Þ
in which Rmin is the lower range of truncation, Rmax is the upper
range of truncation, Rm is the mean of inscribed sphere radii, and
rnd is the standard deviation. Radius and length of pore throats connecting the pore bodies are determined based on the size of neighboring pore bodies. Spacings between the pore body layers of the
network in x-, y- and z-directions are chosen to be variable. Let
Fig. 1. (a) Schematic presentation of a pore body and its connected pore throats.
Truncated parts of the pore body have the width of 2rij, which is the inscribed radius
of pore throat ij. (b) Cross sections of the pore body along the vertices and through
the edges. Radius of inscribed sphere is denoted by R and radius of inscribed circle
0
in the cross section along vertices is denoted by R . (c) Cross section of two pore
bodies and connected pore throats. Geometrical configuration
for determining
the
pffiffi
pffiffi
pore throat radius (rij) based on pore bodies radii, R0i ¼ 26 Ri and R0j ¼ 26 Rj .
the spacings between layers k and k + 1 in the three directions be
denoted by kx,k, ky,l, and kz,m. Then, designating each pore body by
its lattice indices, namely k, l, and m, lattice spacings are given as
follows:
kx;k ¼
pffiffiffi
2 maxfR0 ðk; l; mÞ þ R0 ðk þ 1; l; mÞ : l ¼ 1; ny ; m ¼ 1; nz g;
k ¼ 1; nx
ð5aÞ
ky;l ¼
pffiffiffi
2 maxfR0 ðk; l; mÞ þ R0 ðk; l þ 1; mÞ : k ¼ 1; nx ; m ¼ 1; nz g;
ð5bÞ
pffiffiffi
2 maxfR0 ðk; l; mÞ þ R0 ðk; l; m þ 1Þ : k ¼ 1; nx ; l ¼ 1; ny g;
ð5cÞ
l ¼ 1; ny
kz;m ¼
m ¼ 1; nz
where nx, ny, and nz denote the total number of pore bodies in x, y,
and z directions, respectively. If there is no spacing between vertices
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
4
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
of square cross-sections shown in Fig.
pffiffiffi 1c, the distance between centers of pore bodies will be equal to 2 times the summation of cross
section inscribed radii. Based on the length of a pore throat and
sizes of its neighboring pore bodies, radius of the inscribed circle
in a pore throat cross section is determined (for more detailed
explanation, see (Joekar-Niasar et al., 2008). Consider two pore
bodies i and j, with a center-to-center distance dij (see Fig. 1c),
and cross-sections inscribed radii R0i and R0j , respectively. We define
the dimensionless radii Rei and Rej as follows:
Rei ¼ R0i =dij ;
Rej ¼ R0j =dij
n
ð6Þ
Then, we can calculate dimensionless inscribed radius of the pore
throat ij, reij , as follows:
reij ¼ .i .j
.i
.i1=n þ .1=n
j
;
n>0
ð7Þ
e i sinðp=4Þ
R
¼
e i cosðp=4ÞÞn
ð1 R
.j ¼
ð8Þ
e j sinðp=4Þ
R
e j cosðp=4ÞÞn
ð1 R
ð9Þ
where n is a positive number, which can control the ratio between
radii of a pore body and pore throat, referred to as ‘‘aspect ratio”. A
larger n results in smaller pore throats (larger aspect ratios).
2.1.2. Fluids and network parameters and specifications
Table 1 shows fluids and network properties used in the simulations. Viscosity of the nonwetting fluid is kept constant and equal
n
to 0.001 and viscosity ratio is defined as the M ¼ llw . Statistical
properties of pore bodies, pore throats, and aspect ratio distributions are shown in Table 2. Aspect ratio is defined as pore body inscribed sphere radius divided by pore throat radius. Corresponding
to Table 2, Fig. 2a and b show the aspect ratio distribution as well
as pore body-pore throat size distributions, respectively.
Aspect ratio is controlled by the parameter n given in Eq. (9). In
this work, we have set n = 1.0.
2.2. Governing equations
ð10Þ
A flux Q aij is assigned to a pore throat ij for each fluid separately. A
volume balance for fluid a in pore body i is written as:
X a
Dsa
Vi i ¼
Q ij ;
Dt
j2N
i
a ¼ w; n
Specifications
Ri (mm)
rij (mm)
Rasp
min
max
mean
st. deviation
0.0077
0.0200
0.0125
0.0028
0.0048
0.0162
0.0084
0.0017
1.55
4.00
2.25
0.38
Q aij ¼ K aij Dpaij ;
a ¼ w; n
ð12Þ
a
where K ij is a function of geometry and fluid occupancy of pore
throats and Dpaij is the pressure difference in fluid a across the pore
throats. Summation of Eq. (11) for the two fluids, which are assumed to be incompressible, yields:
X
j2Ni
Q tot
ij ¼
X n
Q ij þ Q w
¼0
ij
ð13Þ
j2Ni
Note that compared to one single pressure per body and one effective conductivity per pore throat (see e.g. (Al-Gharbi and Blunt,
2005; Mogensen and Stenby, 1998), assigning different pressures
and conductivities to each fluid has major advantages. For example,
it allows us to include mechanisms related to the local capillary
pressure (such as snap-off, counter-current flow) in simulations.
w
n
Eqs. (10)–(12) form a determinate set to be solved for sw
i ; pi , and pi .
2.2.2. Pressure field solver
To reduce the computational demand, the general equations are
i , defined as the saturareformulated in terms of a total pressure p
tion-weighted average of fluid pressures in a pore body:
w
n n
i ¼ sw
p
i pi þ si pi
ð14Þ
Combining Eqs. (10) and (14) and using sni þ sw
i ¼ 1, we get the following equations for pressures of wetting and nonwetting fluids:
n c
pw
i ¼ pi si pi
n
w
i þ si pci
pi ¼ p
ð15Þ
ð16Þ
Substitution of Eqs. (15) and (16) in Eq. (12) and subsequently in Eq.
i :
(13) results in an equation for p
2.2.1. General equations for two-phase flow
The local capillary pressure for pore body i is defined as:
w
pci ¼ pni —pw
i ¼ f si
Table 2
Statistical properties of the radii of inscribed spheres of pore bodies (Ri), and of
inscribed circle of pore throats rij, and aspect ratio distribution Rasp.
ð11Þ
where Ni is the set of all pore throats connected to pore body i, Vi is
the volume of pore body i, and sai is the saturation of fluid a in pore
body i. The volumetric flux of fluid a in pore throat ij is given by an
equation similar to Washburn formula:
X
j2Ni
X h n
c
n
w
w
Kw
K ij sw
pi
i K ij 1 si
ij þ K ij ðpi pj Þ ¼
j2Ni
i
w
þ Kw
K nij sw
pcj
j
ij 1 sj
ð17Þ
In this equation, the right-hand side as well as the coefficients of the
left-hand side depend on local saturation only. This linear system of
i by diagonally-scaled biconjugate gradiequations was solved for p
ent method using SLATEC mathematical library (Fong et al., 1993).
2.2.3. Saturation update
i , pressures of fluids can be back-calculated
After calculating p
from Eqs. (15) and (16), based on saturation values. Then,
Table 1
Fluid and network properties used in the simulations.
Parameter
Symbol
Value
Unit
Contact angle
Interfacial tension
Wetting fluid viscosity
Nonwetting fluid viscosity
Total no. of pore bodies in flow direction
Total no. of pore bodies in lateral directions
Domain size
Permeability
h
0.0
0.0725
0.0001, 0.001, 0.01
0.001
45
35
1.9 1.9 2.37
1.43 1012
Degree
kg s2
kg m1 s1
kg m1 s1
–
–
mm3
m2
rnw
lw
ln
nz
nx,ny
–
K
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
(a)
5
6
Frequency (%)
5
4
3
2
1
3
3.
1
3.
2
3.
3
3.
4
3.
5
3.
6
3.
7
3.
8
3.
9
4
2
2.
1
2.
2
2.
3
2.
4
2.
5
2.
6
2.
7
2.
8
2.
9
1.
5
1.
6
1.
7
1.
8
1.
9
0
Aspect Ratio
(b)
20
Pore throat
18
Pore Body
Frequency (%)
16
14
12
10
8
6
4
2
0.
00
9
0.
01
0.
01
1
0.
01
2
0.
01
3
0.
01
4
0.
01
5
0.
01
6
0.
01
7
0.
01
8
0.
01
9
0.
02
7
0.
00
8
6
00
0.
0.
00
5
4
00
0.
3
00
00
0.
0.
0.
00
2
0
Radius (mm)
Fig. 2. Network geometry properties (a) Aspect ratio distribution, (b) pore body and pore throat distributions.
commonly, Eq. (12) is used to calculate Q aij , and using Eq. (11) new
saturation values can be calculated explicitly. This procedure, however, will result in numerical problems for a capillary-dominated
flow regime, as mentioned in Koplik and Lasseter (1985). He found
that the explicit saturation update was not numerically stable for
very small capillary numbers and he could not successfully simulate a capillary-dominated flow. Therefore, we have developed a
semi-implicit approach, analogous to a fractional flow formulation,
to control the nonlinearities under such flow conditions. Details of
the semi-implicit saturation update are given in Joekar-Niasar et al.
(2010a). The resulting discretized equation reads:
!
!
n w
c
c kþ1
X K nij K w
V i X K ij K ij @pij w kþ1
ij @pij
si
þ
sw
j
tot
Dt j2N K tot
@sw
@sw
ij
ij
ij
j2Ni K ij
i
¼
n
V i w k X K ij tot
si þ
tot Q ij
Dt
j2Ni K ij
where superscript k denotes time step level, and
ð18Þ
@pcij
@sw
ij
is calculated
pci —sw
i
relationship of the upstream pore body (given in Section
from
a
2.3.1). One should note that since Q tot
ij and K ij are calculated from
time step k, this scheme is not fully implicit. Here also, the diagonally-scaled biconjugate gradient method from SLATEC mathematical library (Fong et al., 1993) was used to solve Eq. (18).
to a minimum saturation. On the other hand, if local imbibition occurs, the wetting fluid saturation in a pore body can go back to 1.
So, we calculate Dti for all pore bodies, depending on the local process, from the following formulas:
8
w
< qVni sw
for local drainage; Q ni > 0
i si;min
i
Dt i ¼
: Vni 1 sw
for local imbibition; Q ni < 0
i
q
ð19Þ
Dtglobal ¼ minfDti g
ð20Þ
i
where the accumulation rate of the nonwetting fluid is defined as
P
Q ni ¼ j2Ni Q nij . Then, the global time step is chosen to be the minimum of all Dti s.
It should be noted that we have imposed a truncation criterion of
w
w
106 for sw
i si;min and 1 si in order to reduce the number of times
the equations are solved. Also, note that in Eq. (19), there is a
correspondence between saturation change (numerator) and the
accumulation rate of nonwetting fluid (denominator). That is, when
local saturation is close to the limits, the accumulation rate of
nonwetting fluid is also very small. This means that Dti always
remains finite and nonzero.
2.3. Local rules
2.2.4. Time step
The time step is determined based on the time of filling of pore
bodies by the nonwetting fluid (or wetting fluid), denoted by Dti.
The wetting fluid saturation of a pore body varies between 1 and
sw
i;min , as we assume that a pore body may be drained down only
2.3.1. Capillary pressure curves for pore bodies and pore throats
Local capillary pressure within a pore is a function of the curvature of fluid–fluid interface through Young–Laplace equation,
regardless of whether drainage or imbibition occurs. However,
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
6
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
during imbibition, topology of fluid–fluid interfaces is much more
complicated.
Interfaces may be categorized in two different types: ‘‘arc menisci (AM)”, formed along the edges, and ‘‘main terminal menisci
(MTM)”, spanning the pore throat cross section during drainage
(Mason and Morrow, 1987). Behavior of ‘‘MTM” during drainage
and imbibition is different. During drainage, a MTM is formed at
the entrance to a pore throat within a pore body filled with the
nonwetting fluid. However, during imbibition a MTM can be
formed such that it spans over a pore body and connected pore
throats; this is referred to as ‘‘cooperative filling” (Lenormand
and Zarcone, 1984).
In any case, for a given fluid–fluid interface position within a
pore body, we can determine corresponding capillary pressure
and saturation. Therefore, for a given pore body, a unique relationship between capillary pressure and local saturation can be obtained. Local pci —sw
curves for drainage and imbibition are
i
discussed and derived below.
2.3.1.1. Drainage. During drainage, the invasion of a pore throat by
the nonwetting fluid, and consequently the movement of main terminal menisci, is controlled by the entry capillary pressure. In
drainage, an interface is located within a pore body. If a pore body
is filled with both fluids, the wetting fluid is residing along edges of
the pore bodies (see Fig. A.1b and c in Appendix A). The saturation
of the pore body (i.e. volume of the wetting fluid divided by the
volume of the pore body) depends on the prevailing capillary pressure. For a given capillary pressure, the curvature of interfaces in
the edges of the pore body can be calculated and, consequently,
the corresponding saturation can be estimated. In Appendix A, details of derivation of the (local) pci —sw
i relationship for an octahedron pore body are presented.
The resulting pci —sw
i relationship in terms of the radius Ri of the
inscribed sphere of the pore body i and other geometrical parameters is:
pci ¼ 2rnw ji ; ji
8
sw sdr 3:5
>
>
1
1
dr
i
i
>
þ R1
sw
< rij Ri
i P si
1sdr
i
i
¼
a
>
>
sw
dr
>
: R1 sdri ; a ¼ 2:98sdr13:85 smin
< sw
i
i < si
i
i
i
ð21Þ
in which, sdr
i is the wetting fluid saturation corresponding to the inis the
scribed sphere of the pore body, given by Eq. (A.2), and smin
i
minimum possible saturation in a simulation given by Eq. (22).
The resulting curve is shown in Fig. A.2a.
Obviously, it is impossible to completely displace the wetting
fluid from the corners of a pore body. We assume that each pore
body has a minimum saturation sw
i;min , which depends on the im(Pcglobal
posed global pressure difference
defined in Section 3.2) as
well as the blockage of the invading fluid. The capillary pressure
blockage of the invading fluid Pceblock is also a global variable, de-
2.3.1.2. Imbibition. During imbibition, locations of interfaces are
not only controlled by the geometry of pore bodies and pore
throats, but also by local fluid configuration in pore throats connected to a pore body (Lenormand and Zarcone, 1984). One of
the mechanisms that controls configuration of fluid–fluid interfaces in a pore body during imbibition is ‘‘cooperative pore filling”.
Lenormand and Zarcone (1984) observed in micro-model experiments that, depending on the number of pore throats filled with
the nonwetting fluid, the topology of local main terminal interface
can change. They defined a so-called I-mechanism for pore filling,
where I represents number of pore throats filled with the nonwetting fluid. They showed that with the decrease of I, radius of interface may decrease, resulting in an increase of local capillary
pressure as shown schematically in Fig. 3. Some researchers such
as Øren et al. (1998), Hughes and Blunt (2000) suggested simple
algebraic equations for pore filling mechanism, regardless of the
geometry of pore bodies and pore throats. These relations do not
apply to all pore spaces. For example, they do not apply to highporosity domains with small aspect ratio, such as micro-model
experiments of Pyrak-Nolte (2007), as illustrated in simulations
of the imbibition process by Joekar-Niasar et al. (2009).
We define local capillary pressure–saturation relationship as a
function of saturation as well as pore throats filling. Details of derivation of local pci —sw
i relationship during imbibition are given in
Appendix B. Finally, we have suggested the following relationships
which results in a trend consistent with experimental observations
of Lenormand and Zarcone (1984):
pci ¼ 2rnw ji ; ji
8
a
3:5
a
>
simb
sw
simb
simb
>
1
i
i
i
>
þ
< r1ij R1i si dr
Ri
1simb
sdr
i
i
i
¼
a
>
sw
>
1
1
i
>
: R sdr ; a ¼ 2:98sdr 3:85
i
i
i
imb
sw
i P si
imb
smin
< sw
i
i < si
ð23Þ
where simb
is the pore body saturation at which one pore throat is
i
still filled with the nonwetting fluid. It should be noted that this saturation is not predetermined in the simulations. Its value for each
pore body will be determined during the simulation.
2.3.2. Entry capillary pressure for a pore throat
We assume that a pore throat will be invaded by the nonwetting fluid when the capillary pressure in a neighboring pore body
becomes larger than the entry capillary pressure of the pore throat.
For a pore throat with square cross-section, the entry capillary
pressure can be calculated as follows (due to (Mayer and Stowe,
1965; Princen, 1969a,b, 1970; Ma et al., 1996):
pce;ij
¼
rnw
r ij
h þ cos2 h p=4 sin h cos h
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
cos h p=4 h þ sin h cos h
!
ð24Þ
fined to be the minimum entry capillary pressure of all pore
throats neighboring the vicinity of the nonwetting fluid but not invaded by it yet. Thus, using the pci —sw
i relationship given by Eq.
(21), the local minimum wetting fluid saturation in a pore body
may be determined as follows:
dr
sw
i;min ¼ si
Ri
2rnw
n
o2:98sdri 3:85
min Pcglobal ; Pceblock
ð22Þ
A capillary pressure should be also assigned to a pore throat once it
is invaded and both fluids are present. We assume that the capillary
pressure of a pore throat pcij is equal to the capillary pressure of the
upstream pore body.
Fig. 3. Pore filling mechanisms during imbibition referred to I-mechanism. Wetting
fluid is shown as hashed areas. With decrease of number of pore throats filled with
the nonwetting fluid, radius of the interface may decrease. Consequently capillary
pressure assigned to the interface can increase (due to Lenormand and Zarcone,
1984).
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
where h is the contact angle.
2.3.3. Conductivities of pore throats
Conductivities of pore throats are determined based on their
size and fluid occupancy. One of the following two cases may
occur.
(a) A pore throat is occupied by the wetting fluid only. For this
case, the following equation was obtained by Azzam and
Dullien (1977):
Kw
ij ¼
p eff 4
rij
8l w
ij l
ð25Þ
K nij ¼ 0
where lw is the viscosity of the wetting fluid, lij is the length
of pore throat, and
reff
ij
rffiffiffiffi
4
r ij
¼
p
ð26Þ
(b) A pore throat is occupied by both fluids. Then, following
Ransohoff and Radke (1988), we can write:
4 p c 4
r ij
blw
ij l
p eff 4
K nij ¼
r
8ln lij ij
Kw
ij ¼
ð27Þ
ð28Þ
reff
ij
r
nw
pcij
0sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1
2
c2
1 @ rij ð4 pÞr ij
¼
þ r ij A
2
p
ever, we have added five buffer layers at each boundary to
reduce the boundary effect. Thus the network has a length of 45
pore bodies along the main flow direction. The buffer layers are
not included in the averaging window.
3.2. Boundary conditions
For our simulations, we assume that the network is connected
to a nonwetting fluid reservoir on one side and a wetting fluid reservoir on the other side. Fluid pressures are specified at these
boundaries. Side boundary conditions are assumed to be periodic.
For drainage and imbibition simulations the following procedure is followed:
(a) Drainage: Pressure at the nonwetting fluid reservoir is fixed
to Pntop and pressure at the wetting fluid reservoir is fixed to
zero. The difference between the two boundary pressures
during drainage is referred to as ‘‘global pressure difference”
P cglobal ¼ Pntop . After the nonwetting fluid breaks into the wetting fluid reservoir a boundary condition for its pressure is
needed. One may choose to set it equal to the wetting fluid
boundary pressure. But, this will cause an unrealistic relaxation of the nonwetting fluid pressure field throughout the
network. Instead, as an assumption we have chosen to set
the gradient of capillary pressure within an invaded pore
throat (connected
@pc
to the wetting fluid boundary) to be equal
to zero @lijij ¼ 0 .
(b) Imbibition: Similar to the drainage process, pressure at the
nonwetting fluid reservoir is fixed to P ntop and pressure at
where
rcij ¼
7
ð29Þ
the wetting fluid reservoir is fixed to zero. However, P ntop
should be so small that imbibition process can occur continuously. During imbibition similar to the drainage, it is
ð30Þ
assumed that
In Eq. (27), b is a resistance factor that depends on geometry, surface roughness, crevice roundness and other specifications of the
cross section. Detailed explanation about b can be found in Zhou
et al. (1997). As mentioned earlier, the pore throat capillary pressure pcij is set equal to the capillary pressure of the upstream pore
body.
@pcij
@lij
¼ 0 as long as a pore throat at the nonwett-
ing fluid boundary is filled with both fluids.
3.3. Drainage simulations
3. Simulations and analysis
For primary drainage, the network is initially fully saturated with
the wetting fluid. Simulation of drainage starts with raising the pressure of the nonwetting fluid reservoir, and establishing a global pressure difference, P cglobal , across the network. When the pressure
difference is larger than the entry pressure of the largest pore throat
connected to the nonwetting fluid reservoir, drainage starts. In quasi-static simulations, the nonwetting fluid reservoir pressure is increased in incremental steps so that the network is invaded in
steps. At the end of each step, when there is no flow (static conditions), the overall saturation of the network and average specific
interfacial area are determined. Then, global pressure difference is
increased again. In dynamic simulations, the imposed P cglobal is so
large that the whole network (or a large part of it) could be flooded
in one step. The simulations are continued till the change of average
saturation in a selected averaging window is negligible.
For main drainage simulations, the initial saturation occupancy
is based on the last snapshot of the quasi-static main imbibition
experiment. The simulation procedure is otherwise similar to the
primary drainage simulation.
3.1. Network size
3.4. Imbibition simulations
To analyze Darcy-scale equations using pore-network models,
the size of the pore-network should be at least one REV. REV size
was determined by performing quasi-static drainage simulations
in networks with different sizes but with the same statistical
parameters. Our simulations show that the REV size for these statistical parameters is a cube with length of 35 pore bodies. How-
Consider a pore network filled by the nonwetting fluid at the
end of a (main or primary) drainage experiment. The wetting fluid
is still present along edges of pore bodies and pore throats. Starting
from an equilibrium condition, all pore bodies have the same
capillary pressure. The global capillary pressure is decreased by
reducing the nonwetting fluid reservoir pressure or increasing
2.3.4. Snap-off
During imbibition, if the local capillary pressure in a pore throat
becomes smaller than a critical value (defined below), the corner
interfaces become unstable and snap-off will occur. The criterion
for snap-off in a square cross-sectional pore throat has been defined as follows (Vidales et al., 1998):
pcij 6
rnw
r ij
ðcos h sin hÞ
ð31Þ
We assume that as soon as snap-off occurs, the nonwetting fluid retreats instantaneously into the two neighboring pore bodies, and
the pore throat is filled up with the wetting fluid.
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
8
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
Our simulations result in local-scale variables such as pressure,
saturation, and fluxes. These have to be averaged over an averaging
window to obtain macroscopic variables. Average saturation is
simply defined as the ratio of the wetting fluid to the total pore volume of the network:
Pnpb w
Vw
i¼1 si V i
n ¼ Pnpb
V þV
i¼1 V i
ð32Þ
w
Sn ¼ 1 S w
in which npb is the total number of pore bodies within the averaging
window.
The total flux across any given surface is equal to the sum of
fluxes of all pore throats intersecting that surface. The averaging
of pressure is, however, less straightforward. Commonly, average
pressure is obtained using an intrinsic phase average operator
(see e.g. Whitaker, 1977). However, recently it has been shown
that the intrinsic phase average pressure introduces numerical
artifacts when both pressure and saturation are spatially variable
(see Nordbotten et al., 2007, 2008; Korteland et al., 2010; Bottero,
2009). Instead, a centroid-corrected averaging operator has been
suggested by Nordbotten et al. (2008) to alleviate problems associated with intrinsic phase averaging. However, it has been shown
numerically that this approach also suffers from artifacts at large
wetting saturations (Korteland et al., 2010). Thus, we use the
intrinsic phase average since it is still the most commonly used
operator:
Pa ¼
1
dV a
Z
Pnpb
pai sai V i
Pa dV ¼ Pi¼1
;
npb a
dV a
i¼1 si V i
a ¼ n; w
4.1. Non-equilibrium capillarity effects
As mentioned in the Introduction, in standard models of twophase flow in porous media, the difference in average fluid pressures, Pn–Pw, is assumed to be an algebraic function of saturation.
In this section, we show that the curves of Pn–Pw versus Sw depend
strongly on boundary conditions and dynamics of the system, as
shown by Yang et al. (2009) for bundle-of-tube model. We performed a large number of dynamic primary drainage, main drainage, and main imbibition simulations for various viscosity ratios
and under a range of boundary pressure values. For every simulation, we calculated the variations of average saturation and the
‘‘fluids pressure difference” with time. From those results, plots
of Pn–Pw versus saturation were made.
For primary and main drainage simulations, three different global pressure differences were imposed over the domain: 45, 60 and
90 kPa. The entry capillary pressure is 15 kPa. For main imbibition,
four different global pressure differences were imposed over the
domain, equal to +5, 0, 5, 10 kPa. Also, all simulations were performed for three different viscosity ratios, M = 0.1, 1, and 10 (ratio
of viscosity of the nonwetting fluid to that of the wetting fluid).
Computational time increased with decrease of viscosity ratio
M=0.1
30
25
PD
20
MD
MI
15
10
0
0.2
0.4
0.6
0.8
1
Saturation (Sw)
(b) 35
M=1.0
30
25
PD
20
MD
MI
15
10
5
ð33Þ
4. Results and discussion
35
5
0
0.2
0.4
0.6
0.8
1
Saturation (Sw)
(c)
Average Pressure Difference
(kPa)
Sw ¼
(a)
Average Pressure Difference
(kPa)
3.5. Averaging procedure
and decrease of global pressure difference. A simulation took from
24 to 72 h on Intel(R) CPU 6600, 2.4 GHz with 2GB RAM.
Plots of Pn–Pw as well as Pc versus average saturation, Sw, are
shown in Fig. 4. It is clear that Pn–Pw curves strongly depend on
boundary conditions and significantly deviate from the Pc curves.
Their behavior is in agreement with Eq. (3): at any given saturation, the fluids pressure difference is larger than the capillary pressure during drainage and smaller during imbibition. The
differences are larger for larger nonwetting fluid boundary pressure during drainage and for large wetting fluid boundary pressure
Average Pressure Difference
(kPa)
the wetting fluid reservoir pressure. A decrease in the global capillary pressure causes the interfaces to relax gradually and main
imbibition experiment will start. The imbibition simulation will
stop when all pore throats at the outflow boundary (nonwetting
fluid reservoir) are fully filled with the wetting fluid.
For the quasi-static imbibition simulations, an approach similar
to the drainage has been employed.
40
M=10
35
30
25
PD
20
MD
15
MI
10
5
0
0.2
0.4
0.6
0.8
1
w
Saturation (S )
PD:+45 kPa
PD:+60 kPa
MD:+45 kPa
MD:+60 kPa
MI:-100.2
kPa
MI:+5kPa
MI:-5 0.6
kPa
w
Saturation
Pc-S(S
Equil.)
0.4
PD:+90 kPa
MD:+90 kPa
0.8
1
MI:0 kPa
Fig. 4. Average fluids pressure difference curves, Pn–Pw for three different global
pressure differences during primary and main drainage (PD, MD), and four different
global pressure differences during main imbibition (MI), for viscosity ratios (a)
M = 0.1, (b) M = 1 and (c) M = 10 and for ln = 0.001 Pa s. These curves should be
compared to the equilibrium Pc–Sw curves resulted from quasi-static simulations.
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
9
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
Fluid entrapment at pore scale, and consequently fluid distribution at macroscale is controlled by pore-scale invasion mechanisms, such as piston-like movement and snap-off. Snap-off
refers to the mechanisms in which residual fluid can cause a local
rupture in the bulk fluid. The importance of these mechanisms varies depending on the process (drainage or imbibition). For instance,
snap-off is more important during imbibition than drainage, since
the nonwetting fluid (bulk fluid) is the receding one and it is more
vulnerable to rupture. As a result, at the end of imbibition process,
a significant amount of the nonwetting fluid remains behind,
known as the residual saturation. Due to the importance of imbibition in reservoir engineering, the dependence of residual saturation on capillary number, viscosity ratio, contact angle and pores
0.4
0.5
According to Eq. (3), the deviation of non-equilibrium fluids
pressure difference from the macroscopic capillary pressure is
0.3
4.2. Effect of viscosity ratio on fluid distribution
4.3. Dependence of s on fluid viscosities and saturation
0.2
It must be noted that the effective viscosity (or effective resistance)
of a system is related to fluid viscosities as well as spatial distribution of fluids. The definition given in Eq. (34) is only suitable for a
stable invading front where a piston-like movement is dominant,
as the viscosities are linearly weighted by saturation. This equation
is analogous to Kirchhoff’s law for electrical resistors. If the resistors
are in series, the effective resistance is the summation of all resistors (similar to piston-like movement). However, unstable flow
can be analyzed as a circuit with mixed combination of resistors
in parallel and series. Every single capillary blockage increases the
resistance of the system against flow. Consequently, resistance
due to capillary blockage is more pronounced under unstable flow
compared to stable one. For unstable fronts, e.g. in a viscous fingering regime, some researchers, such as Koval (1963), have suggested
other empirical relationships for effective viscosity (see e.g. Fayers
et al., 1990), which are not linearly weighted by saturation. During
both drainage and imbibition, the effective viscosity changes considerably if fluid viscosities are not equal. During primary drainage,
effective viscosity is initially equal to lw, which was chosen to be
different in simulations with different viscosity ratios (i.e. nonwetting fluid saturation was kept constant and lw is varied to get various M values). As the drainage proceeds, effective viscosity
decreases for M = 0.1 and increases for M = 10. At the end of drainage, the effective viscosity is almost equal to ln.
According to Fig. 4 the effect of viscosity ratio on the transient
behavior of fluids pressure difference is significant. During drainage, initially (i.e. at higher saturations), Pn–Pw is larger for small
M and smaller for larger M (for equal ln). This is due to the fact that
the pressure drop over the whole domain for M = 10 is smaller than
the other cases. The decline in the fluids pressure difference during
drainage as residual saturation is approached, as an artefact of the
intrinsic phase averaging. As the nonwetting fluid invades the network, the wetting fluid pressure behind the invasion front increases along with the nonwetting fluid pressure. As a result,
after breakthrough, the average fluid pressure difference may even
decrease. Consequently, the average pressure difference approaches the Pc–Sw curve. These results illustrate that under dynamic conditions, phase pressures difference is not equal to
capillary pressure–saturation curve obtained under equilibrium
conditions. In the other words, the classical two-phase flow simulators, which do not include dynamic effects in their capillary pressure–saturation curve, may not hold under dynamic conditions.
0.1
ð34Þ
Residual Nonwetting Phase (Sor)
leff ¼ ln Sn þ lw Sw
aspect ratio has been studied significantly (see e.g. Dias and
Payatakes, 1986a,b; Vizika et al., 1994; Mogensen and Stenby,
1998; Hughes and Blunt, 2000). The typical curve relating residual
(trapped) nonwetting fluid saturation to the capillary number is
shown in Fig. 5. This curve was obtained from imbibition
n
simulations at M ¼ llw ¼ 1:0 under constant flux boundary conditions. The residual nonwetting fluid saturation was determined
as a function of capillary number at various flow rates. This figure
shows that with the increase of capillary number, the residual
nonwetting fluid saturation decreases. The largest decrease occurs
for capillary numbers between 106 and 104. At higher capillary
numbers, large flow rate cause a piston-like movement and
suppresses the snap-off mechanism.
As mentioned before, all previous studies were performed
under a constant capillary number. In this section, we show that
under Dirichlet boundary condition, a variable spatial distribution
of the trapped phase can be resulted, since the flow rate and
capillary number change continuously in time and space. Our
imbibition results (not presented here) show that, for M = 0.1, with
the increase of wetting fluid saturation, effective viscosity
increases, and consequently flow rate decreases. But, for M = 1
and M = 10 with the increase of wetting fluid saturation, flow rate
increases. This increase is very pronounced for M = 10, which
changes about two order of magnitude.
As shown in Fig. 5, change of flow rate can influence the trapped
nonwetting fluid saturation. This effect can be observed when
examining the average saturation profile along the network, shown
in Fig. 6a and b. Here, the saturation is averaged over a cross section of the network located at position x and then plotted against
x/l at different times, and for two different viscosity ratios. An
interesting result here is the non-monotonic behavior of saturation
for M = 10. Moreover, the saturation front for the case of M = 10 is
steeper than for M = 1.0. With the significant increase of flow rate
for M = 10, snap-off is suppressed. Suppression of snap-off helps
increasing the flooding efficiency. To show the effect of viscosity
ratio on saturation profile during drainage, saturation profiles for
M = 1.0 and 0.1 for DP = 45 kPa have been shown in Fig. 7. As it
can be seen, for M = 0.1, the invading front is unstable and the
slope of the front is much smaller than for M = 1.0. The comparison
between invasion front for these two cases shows how significant
the effect of local heterogeneities on saturation profile is. Consequently, the flooding efficiency is smaller for M = 0.1 compared
with M = 10.
M=1
0
during imbibition. These curves are used later to calculate nonequilibrium capillarity coefficient, s.
Resulting Pn–Pw curves also depend strongly on viscosity ratio
and effective viscosity. In order to explain this, let us define the
effective viscosity, leff, as follows:
1.E-08
1.E-07
1.E-06
1.E-05
1.E-04
1.E-03
Capillary number
Fig. 5. Effect of capillary number on residual nonwetting saturation for
n
M ¼ llw ¼ 1:0.
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
10
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
(a)
1
(a)
M=10
M=0.1
0.8
0.6
Saturation (Sn)
0.8
Saturation (Sw)
1
t=0.018s
0.4
0.6
0.4
0.2
0.2
t=0.31 s
0
0
0
0
0.2
0.4
0.6
0.8
0.2
1
(b)
0.6
0.8
1
Length (X/L)
Length (X/L)
(b)
1
0.4
1
M=1
M=1.0
0.8
0.6
Saturation (Sn)
Saturation (Sw)
0.8
t=0.048 s
0.4
0.6
0.4
0.2
t=0.11 s
0.2
0
0
0
0
0.2
0.4
0.6
0.8
1
Length (X/L)
Fig. 6. Wetting fluid saturation profiles during imbibition for DP = 10 kPa for (a)
M = 10 and (b) M = 1. Dark-colored curves show the wetting fluid saturation
distribution at the time of breakthrough of the wetting fluid. Light-colored curves
show the saturation profile at earlier times. In all cases ln = 0.001 Pa s.
related to the time rate of change of saturation. Hassanizadeh and
Gray (1993a) have suggested that the non-equilibrium capillarity
coefficient, s, is a function of saturation. It is also believed to depend on other medium and/or fluid properties.
To compute s at a given saturation, a plot of Pn–Pw versus @Sw/@t
at that saturation must be made. This is done by starting from one
of the graphs in Fig. 4 for a given viscosity ratio. At a given saturation, different values of Pn–Pw–Pc can be obtained pertaining to different global pressure differences. Corresponding values of @Sw/@t
are known for that saturation. Note that such plots are made for
drainage and imbibition separately, whereby the corresponding
(drainage or imbibition) Pc–Sw curve is used to calculate Pn–Pw–
Pc. By plotting Pn–Pw–Pc versus @Sw/@t at various values of Sw and
fitting a linear relationship to the data points, s is found as the
slope of the fitted line. Results for primary and main drainage as
well as for main imbibition are plotted in Fig. 8a and b,
respectively.
Under drainage, the fitted line passed through the origin for different M values, as it should. However, under imbibition, there was
an intercept. This means that for the curve to go through the origin
a nonlinear relationship should be employed. However, with the
increase of the pressure gradient, this nonlinearity decreases, as
the capillary forces are suppressed by the viscous forces.
Two different aspects of these data are important; (a) order of
magnitude of the non-equilibrium capillarity coefficient s and (b)
its variation with saturation for different viscosity ratios as well
as for drainage and imbibition.
0.2
0.4
0.6
0.8
1
Length (X/L)
Fig. 7. Effect of viscosity ratio on invasion process during drainage. Nonwetting
fluid saturation profiles during drainage for DP = 45 kPa for (a) M = 0.1 and (b)
M = 1. For M = 0.1, the invasion front is unstable and the role of porous medium
geometry in viscous fingering is obvious. In all cases ln = 0.001 Pa s.
Values of s obtained here (ranging from 100 to 1000 Pa s for
M = 1.0) are in agreement with results obtained in other pore-network modeling studies as shown in Table 3. The dimensions of our
network correspond to a porous medium sample size of 1.9
1.9 1.9 mm3 with an intrinsic permeability of 1.43 1012 m2.
While the permeability corresponds to a fine sand, the sample size
is very small. The general understanding is that the magnitude of s
increases with the size of the observation or averaging window and
is inversely correlated with permeability. Dahle et al. (2005), using
a bundle-of-tube model, concluded that s value can be proportional to L2, where L is the length of averaging window. A similar
result was found by Manthey et al. (2005) based on simulations
at continuum scale. In laboratory experiments by Hassanizadeh
et al. (2004), the value of s for a fine sand sample of 3 cm in height
was found to be 5 105 Pa s. The pressure measurements were
actually done by transducers with a diameter of around one centimeter. Using similar transducers in experiments with the same
fine sand, Bottero and Hassanizadeh (2006) found a s value of
around 105 Pa s. But, when Bottero (2009) upscaled the results to
the column scale (18 cm), the average s value was found to be
around 106 Pa s.
Our results shown in Fig. 8a and b illustrate that for the same
nonwetting fluids, smaller viscosity ratio is, larger non-equilibrium
capillarity coefficient will be. Since effective viscosity is defined as
the sum of linear saturation-weighted viscosities, for the same saturation smaller viscosity ratio results in larger effective viscosity.
Thus, we conclude that larger effective viscosity results in larger
non-equilibrium capillarity coefficient (s(M = 0.1) > s(M = 1) >
s(M = 10), ln = constant). This trend is obvious for drainage
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
11
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
(a) 1.E+04
(b)
MI: M=0.1
MI: M=1.0
MI: M=10
1.E+03
(Pa.s)
1.E+03
(Pa.s)
1.E+04
1.E+02
PD:M=0.1
PD:M=1.0
PD:M=10
MD:M=0.1
MD:M=1.0
MD:M=10
1.E+01
1.E+00
0
0.2
0.4
0.6
0.8
1.E+02
1.E+01
1.E+00
0
1
0.2
w
(c)
(d)
1.E+07
1.E+06
1.E+05
1.E+05
PD:M=0.1
PD:M=1.0
PD:M=10
MD:M=0.1
MD:M=1.0
MD:M=10
1.E+04
1.E+03
0.2
0.4
0.6
0.8
1.E+07
1.E+06
0
0.4
1
Saturation (Sw)
Saturation (S )
0.6
0.8
MI: M=0.1
MI: M=1.0
MI: M=10
1.E+04
1.E+03
1
0
0.2
Saturation (Sw)
0.4
0.6
0.8
1
Saturation (Sw)
Fig. 8. Variation of non-normalized (a and b) and normalized (c and d) non-equilibrium capillarity coefficient as a function of saturation for different viscosity ratios
M = 0.1,1,10 for (a and c) primary and main drainage, (b and d) main imbibition simulations. s is normalized by effective viscosity defined as leff = (1 Sw)ln + Swlw. The
vertical axis is in logarithmic scale and in all cases ln = 0.001 Pa s.
Table 3
Values of s in literature for computational and experimental works.
Reference
Exp. type
Process
Fluids
M
Pd (kPa)
K (m2) 1012
s (Pa s)
@s
@Sw
Domain dimensions (cm)
Dahle et al. (2005)
Gielen et al. (2005)
Gielen (2007)
Joekar-Niasar et al. (2010a)
Joekar-Niasar et al. (2010a)
Das et al. (2007)
Manthey et al. (2004)
BTM
PNM
Drain
Drain
–
Oil–water
1
10
0.8
5
4.7
–
274
1.2 105
0
<0
0.1 (L)
0.3 0.3 1
PNM
PNM
CM
Lab
Drain
Drain
Drain
Drain
–
–
Oil–PCE–water
PCE–water
1,10
0.1
>0.6
0.9
4
4
1.2
5.58
150
150
5
3
900,2750
350
105–107
2–7 104
0
0
0
<0
>0
0
0
0
0.5 0.5 0.5
0.5 0.5 0.5
O’Carroll et al. (2005)
O’Carroll et al. (2005)
Bottero and Hassanizadeh (2006)
Bottero (2009)
Camps-Roach et al. (2010)
Camps-Roach et al. (2010)
7
Lab
Lab
Lab
Drain
Drain
Drain
PCE–water
PCE–water
PCE–water
0.8
0.8
0.9
2.1
2.4
6
15.8
12.6
–
5.64 10
1.99 107
105–107
Lab
Lab
Drain
Drain
Air–water
Air–water
0.02
0.02
2.2
4.6
53
14.7
2 105–106
105–8 105
0
0
–
5.07(D), 9.62(L)
5.07(D), 9.62(L)
9.8 (D), 19(L)
10 (D), 20(L)
10 (D), 20(L)
BTM: Bundle of tubes, PNM: Pore-network model, CM: Continuum model, Lab: Laboratory experiment D: Diameter, L: Length.
(Fig. 8a) and imbibition (Fig. 8b). This agrees with the explanation
given by Barenblatt et al. (2003) and Entov (1980), who stated that
dynamic effects in capillary pressure are related to the finite time
required for the fluids in the pore structure to rearrange themselves. Indeed, for larger effective viscosity values, more time is required for fluids to reach the equilibrium condition, which
corresponds to a larger s value.
As it is known, dimension of s is the same as viscosity. Since
resistance of the system is related to viscosities of both fluids
and spatial distribution of fluids, effective viscosity can be used
to normalize the non-equilibrium capillarity coefficient. The normalized s is expected to be less dependent on saturation and
discrepancies in Fig. 8a and b can be reduced. But, results given
in Fig. 8c and d show that the dimensionless s is still a strong function of saturation specially for main drainage and main imbibition
curves. As it can be seen in Fig. 8c, there is not any significant difference among the curves obtained under main drainage. But the
primary drainage curves are not identical, although the order of
magnitude of them are the same. The difference among the curves
is even larger under main imbibition (Fig. 8d). This difference is
probably due to the fact that saturation alone is not sufficient to
characterize the spatial distribution of fluids in the pores. As mentioned earlier, capillary processes and related coefficients (and thus
s coefficient) may not only depend on saturation. As pointed out by
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
12
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
O’Carroll et al. (2010), the magnitude of the observed dynamic effect may be a function of fluid/fluid/solid contact line friction in
addition to fluid viscosities. It is known that contact lines are the
most detailed description of two-phase flow invasion and its
fluid–fluid spatial distribution that include all pore-scale fluid–
fluid properties such as equilibrium contact angle, interfacial tension, displacement distance, and front topology. These pore-scale
issues are too difficult to be measured and simulated. As a macroscopic property, specific interfacial area may be an effective
parameter in addition to effective viscosity and saturation. This
dependence, however, is not being investigated in this paper, as
it requires major computational effort. That is, in addition to main
drainage and imbibition simulations carried out here, one has to
perform a large number of (primary and secondary) scanning
drainage and imbibition simulations to create capillary pressure–
saturation–interfacial area relationship. Such a computational
work is beyond the scope of the present work.
Furthermore, as explained in Eq. (34), effective viscosity defined
as linearly saturation-weighted viscosity does not account for all
nonlinearities related to the fluid distribution as well as macroscopic interface front; the fluids can be distributed as disconnected
blobs and macroscopic interface represent fingering or stable front.
It is interesting to look at the variation of s versus two system
parameters, namely the invading phase saturation (Sinv) and
dimensionless effective viscosity (l), which is defined as the effective viscosity divided by the viscosity of invading fluid. Fig. 9 shows
the contour plot of variation of log(s/leff) versus saturation and l.
The contour plot is obtained by fitting a polynomial function to the
data points of Fig. 8 resulted from drainage and imbibition simulations. The variation ranges for fitted parameters have been calculated for 95% confidence level based on t-distribution. The
trajectories of the primary drainage and main imbibition simulations for different viscosity ratios have been also shown. Since
the effective viscosity is divided by invading fluid’s viscosity, starting points (S) for M = 0.1 and M = 10 are different. But endpoints (E)
are the same, almost equal to unity. As it is known, changing the
invasion process from drainage to imbibition, may change the stability. This figure can be interpreted in analogy to the bifurcation
diagrams. The upper branch is the unstable one, which is the tra-
jectory of drainage
for M = 0.1 (S0.1 ? E) and imbibition for
M = 10 S010 ! E . The lower branch is the stable one, which is the
trajectory
for M = 10 (S10 ? E) and imbibition for
of drainage
M = 0.1 S00:1 ! E . The stable and unstable branches are separated
from each other by the line associated with the trajectory of M = 1.
Fig. 8a shows that during stable primary drainage (M P 1), for a
wide range of saturation, non-equilibrium capillarity coefficient (s)
increases with the decrease of wetting fluid saturation, which is
similar to most findings as reported in Table 3. For M < 1, the reverse trend has been observed; namely s decreases with the decrease of wetting fluid saturation. An empirical formula relating s
to medium and fluid properties was proposed by Stauffer (1978).
In that formula, however, no dependence on saturation was included. We propose to modify Stauffer formula as follows:
s¼
aleff Pcd
kk
qg
2
ð35Þ
where l is replaced by leff, which is a function of saturation. In this
equation, a is a constant, is porosity, k and Pcd are the coefficients
in Brooks-Corey formula (Brooks and Corey, 1964), k is the saturated
permeability, q is the water mass density, and g is the gravity.
Eq. (35) suggests that the change of s with saturation is proportional to the change of leff with saturation. Recalling the simple linear algebraic equation for leff, we can write:
@ leff
@Sw
¼ lw ln
ð36Þ
@l
which means that @S@ sw / @Seff
is negative for M > 1 and positive for
w
M < 1:
@s
/
@Sw
lw ln 6 0; M P 1
lw ln > 0; M < 1
This trend is in full agreement with the observations reported in the
literature, as shown in Table 3. It was found in Mirzaei and Das
(2007) in column-scale drainage simulations, in Joekar-Niasar
et al. (2010a) in pore-network model simulations and in Bottero
and Hassanizadeh (2006) in PCE–water experiments (M 0.9) that
s increased with the decrease of wetting fluid saturation. For
10
4.5
4.75
5
5.25
5.5
5.75
6
6.25
above
S0.1 , S' 10
unstable branch
8
6
4
2
S1.0 , S' 1.0
E
S10 ,S' 0.1
0
0.0
stable branch
0.2
0.4
0.6
0.8
1.0
Saturation (invading)
Fig. 9. Contour plot of log(s/leff) versus l and saturation of the invading phase obtained from dynamic simulations. l is defined as the ratio of leff to the viscosity of the
invading phase. The arrows show the trajectory of drainage simulations from starting point, S, to end point, E, for different viscosity ratios. The corresponding starting points
for imbibition are denoted by S0 . A simple polynomial surface fitted to the drainage and imbibition data points has the following form:
2
s=leff ¼ 6:945 1:526E5 þ 1:623 4:592E5Sinv 5:500 1:906E5l þ 3:591 1:615E5Sinv l 3:494 1:177E5Sinv þ 5:117 1:804E4l2 .
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
13
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
M = 0.1, the inverse trend of s as a function of saturation has been
reported in pore-network modeling results of Joekar-Niasar et al.
(2010a) and in air-water experiments by Camps-Roach et al. (2010).
Main drainage and imbibition curves shown in Fig. 8a and b
show a strong non-monotonic behavior in variation of s versus saturation. This might be due to contribution of spatial distribution of
disconnected phase over the domain to the average phase pressures. While the disconnected phase does not contribute to the
change of saturation, it does contribute to the averaged fluid pressure. The contribution of the trapped phase is obviously more
important under imbibition than drainage. This can be more visibly
seen in main imbibition curves shown in Fig. 8b. As it can be seen
for all viscosity ratios, there is a sharp increase of s with increase of
wetting phase saturation and afterwards s value decreases with
the increase of wetting phase saturation and there is not a dominant trend over the whole range of saturation.
5. Summary and conclusion
We have developed a DYnamic POre-network SImulator for
Two-phase flow (DYPOSIT) to investigate non-equilibrium effects
in fluid pressure fields under dynamic imbibition and drainage
conditions. The pore-network model consists of truncated octahedron pore bodies and parallelepiped pore throats. The angularity in
cross-sections allows for simultaneous flow of both nonwetting
and wetting fluids. Two pressure fields, one for each fluid, are computed separately. Thus, counter-current flow within pore throats
can be simulated. This means that capillary diffusion in the network is properly taken into account. To improve numerical stability of the model under capillary-dominated flow, a semi-implicit
algorithm is employed. This allows us to simulate flow dynamics
for different flow regimes and viscosity ratios for drainage as well
as imbibition.
We find that in spite of the fact that we have not included localscale dynamics of the interface (such as dynamic contact angle),
the average fluids pressure difference (under large pressure gradients) is significantly distinct from the capillary pressure obtained
under equilibrium conditions. Our analysis shows that the invasion
regime has a significant effect on deviation of fluids pressure difference from capillary pressure. This deviation can be much larger
for stable front regime compared to the viscous fingering regime. It
means that in capillary fingering regime or viscous fingering regime less dynamic effect is resulted compared with stable front
regime.
Several numerical simulations for primary drainage, main imbibition, and main drainage are implemented. We have calculated
non-equilibrium capillarity coefficient (s) as a function of saturation, for three different viscosity ratios under Dirichlet boundary
conditions. We found that there is hysteresis in s–Sw relationship;
i.e. although the order of magnitude of s during drainage and imbibition for a given viscosity ratio does not change significantly, the
curves for primary drainage and main imbibition are different. This
is probably due to the fact that saturation alone is not sufficient to
characterize the distribution of fluids in the pores. As mentioned
earlier, capillary processes and related coefficients (and thus s
coefficient) may not only depend on saturation but also on specific
interfacial area. This dependence, however, is not investigated in
this paper, as it requires major computational effort. That is, in
addition to main drainage and imbibition simulations carried out
here, one has to perform a large number of (primary and secondary) scanning drainage and imbibition simulations. Such a computational work is beyond the scope of the present work.
Our analysis also shows a strong dependence of s on viscosity
ratio as well as effective viscosity. Variation of s with saturation
is strongly dependent on the viscosity ratio. This is in agreement
with laboratory experiments reported in Table 3. Given the definition of viscosity ratio as the ratio of nonwetting fluid viscosity to
the wetting fluid’s, if it is smaller than one under drainage or larger
than one under imbibition can create unstable invading fronts.
Since different viscosity ratios lead to different distribution of
phases, interfacial area between the two phases will be different.
Thus, the dependence of s on fluid viscosities will be probably also
modeled by the inclusion of specific interfacial area.
Moreover, we have shown that snap-off is highly related to the
viscosity ratio and capillary number. A viscous fingering during
imbibition under Dirichlet boundary conditions creates non-monotonic distribution of trapped nonwetting fluid. With the invasion of
wetting fluid during imbibition, flooding efficiency increases and
less nonwetting fluid remains in the domain.
Acknowledgments
Comments by one of the anonymous reviewers were most valuable and helped to improve the manuscript. The authors are members of the International Research Training Group NUPUS, financed
by the German Research Foundation (DFG) and The Netherlands
Organization for Scientific Research (NWO). This research is
financed by Utrecht Centre of Geosciences.
Appendix A. Local pci —sw
i relationship for drainage
Consider the octahedron shown in Fig. 1 to be associated with
pore body i. As mentioned earlier, the size of a pore body is specified by the radius of inscribed sphere of the p
octahedron,
Ri. So
ffiffiffi
the length of a typical edge O1A is equal to ai ¼pffiffi 6Ri . Total
pffiffiffi volume
of an an octahedron is then given by V oct:hedi ¼ 32 a3i ¼ 4 3R3i . Next,
consider the corners to be truncated by pore throats connecting to
the pore body. A pore throat ij results in a pyramid with square
base and edge length 2rij to be cut off the octahedron (rij is the radius of inscribed circle p
offfiffiffithe pore throat). The height of the pyramid will be equal to
2rij . Thus, the volume of the truncated
pffiffi
section is equal to 4 3 2 r3ij . So, the volume of the truncated octahedron forming pore body i is given by:
pffiffiffi
pffiffiffi
4 2X
V i ¼ 4 3R3i
r3
j2Ni ij
3
ðA:1Þ
During primary drainage, the wetting fluid saturation of a pore body
specified below. An important
varies from 1 to residual value Smin
i
intermediate saturation corresponds to the situation that the nonwetting fluid fills up the inscribed sphere of the pore body. Denoted
by sdr
i , the corresponding wetting fluid saturation is given by:
sdr
i ¼ 1
4pR3i
p R3
¼ 1 pffiffiffi 3 piffiffiffiP
3V i
3 3Ri 2 j2N r 3ij
i
ðA:2Þ
For the derivation of local capillary pressure–saturation curves, we
identify two different ranges as described below.
w
A.1. Capillary pressure for the range sdr
i 6 si 6 1 (Zone 1)
We assume that the entrance of nonwetting fluid into a pore
body and the position of fluid–fluid interface may be idealized by
sequences shown in Fig. A.1a. Position 1 shows the moment of
invasion of nonwetting fluid into the pore body. The interface
has the same radius of curvature as when it was in the pore throat
ij and the local capillary pressure is the same as the (entry) capillary pressure of the pore throat, denoted by pce;ij (the starting point
w
of pci —sw
i curve in Fig. A.2a at si ¼ 1). We assume that as the interface moves into the pore body, its radius remains unchanged until
position 2 is reached. At this point, the saturation of nonwetting
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
14
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
r 3ij
2
3
r3
p Viji . Thus, for the range
2
1 P sw
i P 1 3 p V i , the local capillary pressure will be constant
equal to the entry capillary pressure of the pore throat; pci ¼ pce;ij
(see the plateau in Fig. A.2a). As the interface moves into the pore
body, it will expand and its capillary pressure decreases (interface
3 in Fig. A.1a). To simplify the geometry of the interface, we
assume that from this saturation onwards, the nonwetting fluid
fills up a sphere of increasingly larger radius until it fills up the
inscribed sphere of the pore body, at which point the wetting fluid
(a)
ðA:3Þ
r ci
In this range, as
increases, the local capillary pressure, given by
pci ¼ 2rnw =r ci decreases until sdr
i is reached (see Fig. A.2).
The irregular shape of the resulting curves is approximated by
the following formula for the range 1 to sdr
i for sake of computational efficiency.
sw
i
¼
sdr
i
þ 1
sdr
i
1 sdr
i
1 sdr
i
4prci 3
3V i
ðA:4Þ
4pr 3ij
3V i
2
1.5
1
Si dr
0
0
0.2
0.4
0.6
0.8
1
Saturation (Sw)
(b)
Imbibition_Analytical
Imbibition_Fitted
3
2.5
Mean Curvature ( )
4pr ci 3
3V i
Drainage_Analytical
Drainage_Fitted
0.5
saturation is sdr
i given by Eq. (A.2).
If the radius of such an intermediate sphere is r ci , then the
wetting fluid saturation in this range is given by the following
relation:
sw
i ¼ 1
3
2.5
Mean Curvature ( )
fluid in the pore body is equal to
2
1.5
1
0.5
A.2. Capillary pressure for the range 0 <
sw
i
6
sdr
i
(Zone 2)
Once the nonwetting fluid fills the inscribed sphere, and after
that, the wetting fluid will be present only along the edges of the
pore body. So, interfaces will be formed along pore body’s 12 edges
as well as in the opening of pore throats not invaded by the nonwetting fluid yet. The latter interfaces are formed in the vertices of
the truncated octahedron (where non-invaded pore throats are
connected). Edge and vertex interfaces are shown in Fig. A.1b
and d, respectively.
Edge interfaces form part of a cylinder (with only one finite radius of curvature), and vertex interfaces form part of a spherical
surface. When all pore throats have been invaded by the nonwetting fluid, there will be no vertex interface in a pore body.
Si imb
0
0
0.2
0.4
0.6
0.8
1
w
Saturation (S )
Fig. A.2. Local pci —sw
i curves during drainage (a) and imbibition (b). Thick curves
show the analytical calculations and thin curves show the fitted curves.
Let the radius of edge interfaces be denoted by rci . We assume
that the radius of curvature of interfaces along the edges is negligible compared r ci . As shown in Fig. A.1b, we denote the half
corner angle of all edges by b and half angle of interfaces by
a. We can show that
NMO = b = 0.5cos1(1/3) = 0.9553, and
0
0
NOM = a = (0.61548 h). For given b and rci , area of the MNM N
(a)
(b)
(c)
(d)
Fig. A.1. Schematic presentation of interfaces in a pore body (a) expansion of an interface into a pore body before filling the inscribed sphere (b and c) interface geometry
along the edges (d) interface geometry in a vertex of pore body when the corresponding pore throat is filled with wetting fluid.
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
is equal to 2 ðAMMNO A
follows:
AMMNO ¼ 0:5r ci
A
ONM0
2
ONM0 Þ,
which are themselves defined as
B.1. All pore throats are (partially) filled by the nonwetting fluid, I = Z
(Zone 1)
ðA:5Þ
In this range all pore throats are (partially) filled by the nonwetting fluid and pci —sw
relationship is the same as Zone 2
i
dr
for drainage process. Thus, pci —sw
0 < sw
i 6 si
i curve follows the
same curve shown in Fig. A.2(a) till one of the pore throats is fully
filled by the wetting fluid.
sin a cos h
sin b
2
¼ 0:5r ci a
ðA:6Þ
Consequently, the cross-sectional area of residual wetting fluid may
be written as:
Awet
¼ r ci
i
2
15
sin a cos h
a
sin b
ðA:7Þ
Total length of edges filled with the wetting fluid is given by:
X
pffiffiffi
Li ¼ 12 6Ri 8
r ij
j2Ni
ðA:8Þ
With the total pore body volume given by Eq. (A.1), the saturation of
wetting fluid in the pore body is determined as:
pffiffiffi pffiffiffi
P
a cos h a
3 2rci 2 3 6Ri 2 j2Ni r ij sinsin
b
pffiffiffi
sw
P
i ¼
3 6R3i 2 j2Ni r 3ij
ðA:9Þ
This relation is valid only when all pore throats connected to the
pore body i are invaded by the nonwetting fluid, so that there is
no vertex interface in that pore body. This may not be the case when
the nonwetting
fluid has filleddronly the inscribed sphere of the pore
body r ci ¼ Ri for which sw
i ¼ si given by Eq. (A.2). In fact, substituting r ci ¼ Ri in Eq. (A.9), results in a saturation different from sdr
i . We
denote this saturation by si . Forcing saturation to vary between sdr
i
and smin
(the latter one was defined by Eq. (22)), results in the foli
lowing equation:
min
sdr
i si
si smin
i
2 pffiffiffi pffiffiffi
3
P
a cos h a
3 2 3 6Ri 2 j2Ni r ij sinsin
b
5
pffiffiffi
smin
4
P
i
4j2i 3 6R3i 2 j2Ni r 3ij
min
sw
þ
i ¼ si
B.2. More than one of the pore throats and not all of them are
(partially) filled by the nonwetting fluid, 1 < I < Z (Zone 2)
As long as all the pore throats are (partially) filled by the nonwetting fluid, the wetting fluid will be present only along the 12
edges of the pore body. When nonwetting fluid recedes completely
from one of the pore throats, an interface will be created in the
opening of that pore throat. This vertex interfaces is shown in
Fig. A.1d. For simplicity, we assume that all wetting fluid’s volume
existing in that pore body will remain in that vertex. If nonwetting
fluid recedes from two or more pore throats, the wetting fluid will
be distributed among them. Knowing the volume of the wetting
fluid and the number of pore throats nonwetting fluid has receded
from, radius of the curvature of the interface can be estimated. Re0
0
0
0
fer to Fig. A.1d, let’s denote the length of edge of pyramid O1A B C D
0
0
base by a . We can calculate a based on the geometry of pore body,
pore throats, and the number of pore throats partially filled with
nonwetting fluid (I), as follows:
a0 ¼
" pffiffiffi
#1=3
6I
3 2 w
8 X
r3ij
si V i þ
6I
6 I j¼1
0
in which Vi can be calculated from Eq. (A.1). Using the value of a ,
mean radius of curvature of the interface is given by
0
a0
pffi . As it can be seen with decrease of I, a
rci ¼
will
2 cos hþsin1
ðA:10Þ
This relationship is shown as a thick curve in Fig. A.2a.
A.3. Local pci —sw
i relationship for the full range of saturation
Finally, to reduce the computational efforts, we have chosen to
fit one single curve to pci for the full range of saturation. Eq. (21)
shows the fitted pci —sw
i relationship for the full range of saturation
in a pore body. This relationship is shown by a thin curve in
Fig. A.2a.
Appendix B. Local pci —sw
i relationship for imbibition
Local variation of capillary pressure is much more complex under imbibition than drainage. As mentioned in Section 2.3.4, interface topology, and consequently capillary pressure in a pore body
are controlled by the number of pore throats which are still (partially) filled by the nonwetting fluid. Introducing Z as the total
number of pore throats connected a pore body, and I as the number
of pore throats (partially filled by the nonwetting fluid), the following zones are introduced to defining local pci —sw
i relationship:
All pore throats are (partially) filled by the nonwetting fluid,
I = Z (Zone 1).
More than one of the pore throats and not all of them are (partially) filled by the nonwetting fluid, 1 < I < Z (Zone 2).
Only one pore throat is (partially) filled by the nonwetting fluid,
I = 1 (Zone 3).
ðB:1Þ
3
3
decrease and consequently capillary pressure will increase. This
trend is shown schematically in Fig. A.2b as a step-wise curve.
The jumps in the curve occur when the nonwetting fluid recedes
from one of the pore throats (I decreases), and refer to Eq. (B.1), radius of curvature decreases. We follow this trend till only one of the
pore throats remains filled with the nonwetting fluid.
There is a qualitative consistency between the trend of this local
pci —sw
i curve with the experimental observations of Lenormand and
Zarcone (1984), who observed with decrease of number of pore
throats (partially) filled with the nonwetting fluid, local capillary
pressure may increase.
B.3. Only one pore throat is (partially) filled by the nonwetting fluid,
I = 1 (Zone 3)
When only one of the pore throats is (partially) filled with the
nonwetting fluid, the behavior of the interface follows the explanations given in Zone I during drainage.
B.4. Local pci —sw
i relationship for the full range of saturation
As mentioned before, the variation of capillary pressure with
saturation is very irregular and it depends on the fluid occupancy
of pore throats. The thick curve in Fig. A.2b shows an example of
this curve for a pore body. But, handling such a discontinuous
curve in the pore-network model is nontrivial. Thus, to reduce
the computational efforts, we have chosen to fit one single curve
to pci for the full range of saturation. Eq. (23) shows the fitted
pci —sw
i relationship for the full range of saturation in a pore body.
This relationship is shown by a thin curve in Fig. A.2b. The fitted
curve results decrease of capillary pressure with increase of
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
16
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
saturation to a minimum capillary pressure. This minimum capillary pressure occurs at a saturation corresponding to the start of
Zone 3 during imbibition, where only one pore throat is (partially)
filled by the nonwetting fluid. We denote this saturation by simb
as
i
shown in Fig. A.2b.
References
Aker, E., Maloy, J., Hansen, A., Batrouni, G.G., 1998a. A two-dimensional network
simulator for two-phase flow in porous media. Transport Porous Media 32,
163–186.
Aker, E.K., Maloy, K.J., Hansen, A., 1998b. Simulating temporal evolution of pressure
in two-phase flow in porous media. Phys. Rev. E 58, 2217–2226.
Al-Gharbi, M.S., Blunt, M.J., 2005. Dynamic network modeling of two-phase
drainage in porous media. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 71,
016308.
Avraam, D.G., Payatakes, A.C., 1995a. Generalized relative permeability coefficients
during steady-state, two-phase flow in porous media and correlation with the
flow mechanisms. Transport Porous Media 20, 135–168.
Azzam, M.I.S., Dullien, F.A.L., 1977. Flow in tubes with periodic step changes in
diameter: a numerical solution. Chem. Eng. Sci. 32, 1445–1455.
Bakke, S., Øren, P.-E., 1997. 3-d pore-scale modelling of sandstones and flow
simulations in the pore networks. SPE J. 2, 136–149.
Barenblatt, G., Patzek, T., Silin, D., 2003. The mathematical model of nonequilibrium
effects in water–oil displacement. SPE J. 8, 409–416.
Berentsen, C.W.J., Hassanizadeh, S.M., 2006. On the equivalence of measured and
averaged pressures in porous media. In: 16th International Conference on
Computational Methods in Water Resources.
Bottero, S., 2009. Advances in Theories of Capillarity in Porous Media. Ph.D. Thesis,
Utrecht University.
Bottero, S., Hassanizadeh, S.M., 2006. Experimental and numerical study to
investigate dynamic capillary pressure effect in two-phase flow in porous
media. In: 16th International Conference on Computational Methods in Water
Resources.
Brooks, R.H., Corey, A.T., 1964. Hydraulic Properties of Porous Media. Tech. Rep.
Hydrol. Pap. No. 3, Colo. State Univ.
Camps-Roach, G., O’Carroll, D.M., Newson, T.A., Sakaki, T., Illangasekare, T.H., 2010.
Experimental investigation of dynamic effects in capillary pressure: grain size
dependency and upscaling. Water Resour. Res. 46, W08544.
Chen, D.Q., Pyrak-Nolte, L.J., Griffin, J., Giordano, N.J., 2007. Measurement of
interfacial area per volume for drainage and imbibition. Water Resour. Res. 43,
W12504.
Cheng, J.T., Pyrak-Nolte, L.J., Nolte, D.D., 2004. Linking pressure and saturation
through interfacial area in porous media. Geophys. Res. Lett. 31, L08502.
Constantinides, G.N., Payatakes, A.C., 1996. Network simulation of steady-state twophase flow in consolidated porous media. AIChE J. 42, 369–382.
Dahle, H.K., Celia, M.A., 1999. A dynamic network model for two-phase immiscible
flow. Comput. Geosci. 3, 1–22.
Dahle, H.K., Celia, M.A., Hassanizadeh, S.M., 2005. Bundle-of-tubes model for
calculating dynamic effects in the capillary–pressure–saturation relationship.
Transport Porous Media 58, 5–22.
Das, D., Gauldie, R., Mirzaei, M., 2007. Dynamic effects for two-phase flow in porous
media: fluid property effects. AIChE J. 53, 2505–2520.
Das, D.B., Mirzaei, M., Widdows, N., 2006. Non-uniqueness in capillary pressure–
saturation–relative permeability relationships for two-phase flow in porous
media: interplay between intensity and distribution of random microheterogeneities. Chem. Eng. Sci. 61, 6786–6803.
Dias, M.M., Payatakes, A.C., 1986a. Network models for two-phase flow in porous
media part 1. Immiscible microdisplacement of non-wetting fluids. J. Fluid
Mech. 164, 305–336.
Dias, M.M., Payatakes, A.C., 1986b. Network models for two-phase flow in porous
media part 2. Motion of oil ganglia. J. Fluid Mech. 164, 337–358.
Entov, V.M., 1980. Theory of nonequilibrium effects associated with the flow of
nonuniform fluids in porous media. Fluid Dynam. 15, 365–369.
Fayers, F.J., Blunt, M.J., Christie, M.A., 1990. Accurate calibration of empirical viscous
fingering models. In: 2nd European Conference on the Mathematics of Oil
Recovery, 75015 Paris, pp. 45–55.
Fong, K.W., Jefferson, T.H., Suyehiro, T., Walton, L., 1993. Guide to the slatec
common mathematical library. <http://www.netlib.org/slatec>.
Gielen, T., Hassanizadeh, S.M., Leijnse, A., Nordhaug, H.F., 2005. Dynamic effects in
multiphase flow: a pore-scale network approach. In: Das, D.B.M.H.S. (Ed.),
Upscaling Multiphase Flow in Porous Media. Springer, The Netherlands, pp.
217–236.
Gielen, T.W.J., 2007. Upscaling Multiphase Transport Processes in Porous Media.
Ph.D. Thesis, Delft University of Technology.
Hassanizadeh, S., Celia, M., Dahle, H., 2002. Dynamic effects in the capillary
pressure–saturation relationship and their impacts on unsaturated flow. Vadose
Zone J. 1, 38–57.
Hassanizadeh, S.M., Gray, W.G., 1990. Mechanics and thermodynamics of
multiphase flow in porous media including interphase boundaries. Adv.
Water Resour. 13, 169–186.
Hassanizadeh, S.M., Gray, W.G., 1993a. Thermodynamic basis of capillary pressure
in porous media. Water Resour. Res. 29, 3389–3405.
Hassanizadeh, S.M., Gray, W.G., 1993b. Toward an improved description of the
physics of two-phase flow. Adv. Water Resour. 16, 53–67.
Hassanizadeh, S.M., Oung, O., Manthey, S., 2004. Laboratory experiments and
simulations on the significance of non-equilibrium effect in the capillary
pressure–saturation relationship. In: Schanz, T. (Ed.), Unsaturated Soils:
Experimental Studies. Proceedings of the International Conference: From
Experimental Evidence towards Numerical Modeling of Unsaturated Soils, vol.
1. Springer Verlag, Berlin, pp. 3–14.
Held, R.J., Celia, M.A., 2001. Modeling support of functional relationships between
capillary pressure, saturation, interfacial area and common lines. Adv. Water
Resour. 24, 325–343.
Hughes, R.G., Blunt, M.J., 2000. Pore scale modeling of rate effect in imbibition.
Transport Porous Media 40, 295–322.
Joekar-Niasar, V., Hassanizadeh, S., Leijnse, A., 2008. Insights into the
relationships among capillary pressure, saturation, interfacial area and
relative permeability using pore-network modeling. Transport Porous Media
74, 201–219.
Joekar-Niasar, V., Hassanizadeh, S.M., Dahle, H.K., 2010a. Non-equilibrium effects in
capillarity and interfacial area in two-phase flow: dynamic pore-network
modelling. J. Fluid Mech. 655, 38–71.
Joekar-Niasar, V., Hassanizadeh, S.M., Pyrak-Nolte, L.J., Berentsen, C., 2009.
Simulating drainage and imbibition experiments in a high-porosity
micromodel using an unstructured pore network model. Water Resour. Res.
45, W02430.
Joekar-Niasar, V., Prodanović, M., Wildenschild, D., Hassanizadeh, S.M., 2010b.
Network model investigation of interfacial area, capillary pressure and
saturation relationships in granular porous media. Water Resour. Res. 46,
W06526.
Kalaydjian, F., Marle, C.M., 1987. Thermodynamic aspects of multiphase flow in
porous media. Collection colloques et séminaires – Institut français du pétrole
45, 513–531.
Koplik, J., Lasseter, T.J., 1985. Two-phase flow in random network models of porous
media. Soc. Petrol. Eng. J. 25, 89–110.
Korteland, S., Bottero, S., Hassanizadeh, S.M., Berentsen, C.W.J., 2010. What is the
correct definition of macroscale pressure? Transport Porous Media 84, 153–
175.
Koval, E.J., 1963. A method for predicting the performance of unstable miscible
displacements in heterogenous media. Trans. AIME 219, 145–150.
Lenormand, R., Zarcone, C., 1983. Mechanism of the displacement of one fluid by
another in a network of capillary ducts. J. Fluid Mech. 135, 337–353.
Lenormand, R., Zarcone, C., 1984. Role of roughness and edges during imbibition in
square capillaries. In: 59th Annual Technical Conference of the SPE. Houston,
Richardson, TX, No. 13264.
Lewalle, J., Singh, K., Bambacht, J., 1994. Analysis of the wet pressing of paper pulp.
Int. J. Multiphase Flow 20, 415–437.
Ma, S., Mason, G., Morrow, N.R., 1996. Effect of contact angle on drainage
and imbibition in regular polygonal tubes. Colloids Surfaces A 117, 273–
291.
Manthey, S., Hassanizadeh, S.M., Helmig, R., 2005. Macro-scale dynamic effects in
homogeneous and heterogeneous porous media. Transport Porous Media 58,
121–145.
Manthey, S., Hassanizadeh, S.M., Oung, O., Helmig, R., 2004. Dynamic capillary
pressure effects in two-phase flow through heterogeneous porous media. In:
15th International Conference on Computational Methods in Water Resources.
Mason, G., Morrow, N.R., 1987. Meniscus configurations and curvatures in nonaxisymmetric pores of open and closed uniform cross section. Proc. Roy. Soc.
Lond. Series A, Math. Phys. Sci. 414, 111–133. http://www.jstor.org/stable/
2398111.
Mayer, R.P., Stowe, R.A., 1965. Mercury porosimetry-breakthrough pressure for
penetration between packed spheres. J. Colloid Sci. 20, 891–911.
Mirzaei, M., Das, D.B., 2007. Dynamic effects in capillary pressure–saturations
relationships for two-phase flow in 3d porous media: implications of microheterogeneities. Chem. Eng. Sci. 62, 1927–1947.
Mogensen, K., Stenby, E.H., 1998. A dynamic two-phase pore-scale model for
imbibition. Transport Porous Media 32, 299–327.
Nguyen, V.H., Sheppard, A.P., Knackstedt, M.A., Pinczewski, W., 2006. The effect of
displacement rate on imbibition relative permeability and residual saturation. J.
Petrol. Sci. Eng. 52, 54–70.
Nordbotten, J., Celia, M., Dahle, H., Hassanizadeh, S., 2007. Interpretation of
macroscale variables in Darcy’s law. Water Resour. Res. 43, W08430.
Nordbotten, J., Celia, M., Dahle, H., Hassanizadeh, S., 2008. On the definition of
macroscale pressure for multiphase flow in porous media. Water Resour. Res.
44, W06S02.
Nordhaug, H.F., Celia, M., Dahle, H.K., 2003. A pore network model for calculation of
interfacial velocities. Adv. Water Resour. 26, 1061–1074.
O’Carroll, D.M., Phelan, T.J., Abriola, L.M., 2005. Exploring dynamic effects in
capillary pressure in multistep outflow experiments. Water Resour. Res. 41,
W11419.
O’Carroll, D.M., Mumford, K.G., Abriola, L.M., Gerhard, J.I., 2010. Influence of
wettability variations on dynamic effects in capillary pressure. Water Resour.
Res. 46, W08505.
Øren, P.E., Bakke, S., Arntzen, O.J., 1998. Extending predictive capabilities to network
models. SPE J. 3, 324–336.
Oung, O., Hassanizadeh, S.M., Bezuijen, A., 2005. Two-phase flow experiments in a
geocentrifuge and the significance of dynamic capillary pressure effect. J. Porous
Media 8, 247–257.
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007
V. Joekar-Niasar, S. Majid Hassanizadeh / International Journal of Multiphase Flow xxx (2010) xxx–xxx
Pereira, G.G., Pinczewski, W.V., Chan, D.Y.C., Paterson, L., Øren, P.E., 1996. Pore-scale
network model for drainage-dominated three-phase flow in porous media.
Transport Porous Media 24, 167–201.
Porter, M.L., Schaap, M.G., Wildenschild, D., 2009. Lattice-boltzmann simulations of
the capillary pressure–saturation–interfacial area relationship for porous
media. Adv. Water Resour. 32, 1632–1640.
Prat, M., 2002. Recent advances in pore-scale models for drying of porous media.
Chem. Eng. J. 86, 153–164.
Princen, H.M.J., 1969a. Capillary phenomena in assemblies of parallel
cylinders i. capillary rise between two cylinders. J. Colloid Interface Sci.
30, 69–75.
Princen, H.M.J., 1969b. Capillary phenomena in assemblies of parallel cylinders ii.
capillary rise in systems with more than two cylinders. J. Colloid Interface Sci.
30, 359–371.
Princen, H.M.J., 1970. Capillary phenomena in assemblies of parallel cylinders iii.
liquid columns between horizontal parallel cylinders. J. Colloid Interface Sci. 34,
171–184.
Pyrak-Nolte, L.J., 2007. Url of Website with Original Data. <http://www.physics.
purdue.edu/rockphys/dataimages/>. <http://www.physics.purdue.edu/rockphys/
DataImages/>.
Ransohoff, T.C., Radke, C.J., 1988. Laminar flow of a wetting liquid along the corners
of a predominantly gas-occupied noncircular pore. J. Colloid Interface Sci. 121,
392–401.
Raoof, A., Hassanizadeh, S., 2010. A new method for generating pore-network
models of porous media. Transport Porous Media 81, 391–407.
17
Reeves, P.C., Celia, M.A., 1996. A functional relationship between capillary pressure,
saturation, and interfacial area as revealed by a pore-scale network model.
Water Resour. Res. 32, 2345–2358.
Singh, M., Mohanty, K.K., 2003. Dynamic modelling of drainage through threedimensional porous materials. Chem. Eng. Sci. 58, 1–18.
Stauffer, F., 1978. Time dependence of the relations between capillary pressure,
water content and conductivity during drainage of porous media. IAHR.
Thessaloniki, pp. 3.35–3.52.
Thompson, K.E., 2002. Pore-scale modeling of fluid transport in disordered fibrous
materials. AIChE J. 48, 1369–1389.
Van der Marck, S.C., Matsuura, T., Glas, J., 1997. Viscous and capillary
pressures during drainage: network simulations and experiments. Phys.
Rev. E 56, 5675–5687.
Vidales, A.M., Riccardo, J.L., Zgrabli, G., 1998. Pore-level modelling of wetting on
correlated porous media. J. Phys. D – Appl. Phys. 31, 2861–2868.
Vizika, O., Avraam, D.G., Payatakes, A.C., 1994. On the role of viscosity ratio during
low-capillary-number forced imbibition in porous media. J. Colloid Interface Sci.
165, 386–401.
Whitaker, S., 1977. Simultaneous heat, mass, and momentum transfer in porous
media: a theory of drying. Adv. Heat Transfer 13, 119–200.
Yang, D., Currier, R.P., Zhang, D.Z., 2009. Ensemble phase averaged equations for
multiphase flows in porous media. Part 1: the bundle-of-tubes model. Int. J.
Multiphase Flow 35, 628–639.
Zhou, D., Blunt, M.J., Orr, F.M., 1997. Hydrocarbon drainage along corners of
noncircular capillaries. J. Colloid Interface Sci. 187, 11–21.
Please cite this article in press as: Joekar-Niasar, V., Majid Hassanizadeh, S. Effect of fluids properties on non-equilibrium capillarity effects: Dynamic porenetwork modeling. Int. J. Multiphase Flow (2010), doi:10.1016/j.ijmultiphaseflow.2010.09.007