Excerpt
Implementation and validation of
gravity effects
in a porenetwork solver
Bilal Özkan Lafci
ETH BSc Mechanical and Process Engineering
Bachelor Thesis  Spring 2012
Institute of Fluid Dynamics &
Separation Process Laboratory
ETH Zurich
With the aid of a gravity implemented pressure solver and a pore network model, we have been
able to simulate and analyze the mechanism of fluid displacement in porous media. Exemplarily,
the pressure solver is used to simulate the injection of ethyl acetate and heptane into the water
and glass beads filled CCSshowcase, which is an experimental setup to understand the process
behind carbon dioxide storage. In a first step, gravity effects are implemented into the flow rate
and thus change the local capillary pressure by altering the wetting and nonwetting saturation.
The validity of this approach is shown by a series of Test Cases which are stepwise expanded to
more complex pore networks. A parameter study is done to understand the spreading behavior
of less dense, less viscous and with high interfacial tension afflicted nonwetting fluid injections
of different flow rates. Results of this parameter study allow comprehending the influence of
capillary, gravity and viscous forces during fluid displacement in the CCSshowcase. For this
reason the showcase glass beads are modeled by the conceptual approach of the general site bond
correlated pore network model hypon, which relates pore throat and pore body size distribution
to generate a regular isotropic pore network as simulation basis. Finally, the similarity in the
plume of the experiments and simulations is investigated and shows consistency in the case of
fluids with low interfacial tension like ethyl acetate, but requires a more precise capture of the
network properties for fluids with high interfacial tension and low viscosity like heptane.
Contents
1
Introduction
5
2
Pore network solver with gravity
6
2.1
Physical displacement mechanism . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.2
Local capillary pressure p
c
i
(s
w
i
) (Lee et al. [2009]) . . . . . . . . . . . . . . . . . .
8
2.3
Gravity affected flow rate between pore bodies . . . . . . . . . . . . . . . . . . . .
10
2.4
Mathematical Equations and algorithm . . . . . . . . . . . . . . . . . . . . . . . .
11
2.4.1
Assumptions
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.4.2
Boundary Conditions
. . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.4.3
Governing equations
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.4.4
Pressure Field Solver
. . . . . . . . . . . . . . . . . . . . . . . . . . .
13
2.4.5
Timestep Calculation
. . . . . . . . . . . . . . . . . . . . . . . . . .
14
2.4.6
Algorithm
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
3
Validation and code development
17
3.1
Test Case 1  Macroscale pore column (single phase) . . . . . . . . . . . . . . . .
17
3.2
Test Case 2  Microscale pore column (single phase & low interfacial tensions) . .
18
3.3
Test Case 3  Microscale pore column (single phase & high interfacial tensions) .
18
3.4
Test Case 4  Macroscale pore column (multiphase) . . . . . . . . . . . . . . . . .
19
3.5
Test Case 5  Microscale pore column (multiphase & low interfacial tension) . . .
20
3.6
Test Case 6  Microscale pore column (multiphase & high interfacial tension) . .
21
3.7
Test Case 7 & 8  Microscale pore column with injection of less dense fluid . . . .
21
3.8
Test Case 9  2DPore field with equal sized pore bodies (low interfacial tension)
22
3.9
Test Case 10  2DPore field with equal sized pore bodies (high interfacial tension) 22
3.10 Test Case 11  2DPore field with poly sized pore bodies (low interfacial tensions)
24
3.11 Test Case 12  2DPore field with poly sized pore bodies (high interfacial tensions) 24
4
Parameter study on fluid displacement
27
4.1
Drainage Phasediagram for immiscible displacement after Lenormand et al. [1988] 28
4.2
Influence of Interfacial Tension . . . . . . . . . . . . . . . . . . . . . . . . . . . .
30
4.3
Influence of the viscosity ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
4.4
Influence of inlet pressure and injection rate . . . . . . . . . . . . . . . . . . . . .
34
4.4.1
Displacement evolution
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
4.4.2
Influence of flow rate variations . . . . . . . . . . . . . . . . . . . . . . . .
35
5
Simulation of the CCSShowcase
39
5.1
Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
3
5.2
Pore network modeling of CCSShowcase . . . . . . . . . . . . . . . . . . . . . . .
40
5.2.1
HYPONmodel after Acharya et al. . . . . . . . . . . . . . . . . . . . . . .
41
5.2.2
Assumptions for the pore network
. . . . . . . . . . . . . . . . . . . . . .
43
5.2.3
Generation of the pore network . . . . . . . . . . . . . . . . . . . . . . . .
43
5.3
Simulation of the experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
44
5.3.1
Physical properties of injected fluids . . . . . . . . . . . . . . . . . . . . .
44
5.3.2
Ethyl acetate with a flow rate of 150
ml
min
and 2mm glass beads . . . . . .
44
5.3.3
Heptane with a flow rate of 150
ml
min
and 2mm glass beads . . . . . . . . .
46
5.3.4
Comparison of ethyl acetate and heptane
. . . . . . . . . . . . . . . . . .
47
5.4
Simulation of a High Pressure Showcase . . . . . . . . . . . . . . . . . . . . . . .
48
6
Conclusion
50
7
Outlook
51
4
1 Introduction
The industrial revolution and its greatest invention the steam machine marked a turning point
of human existence in societies of developed countries. An agricultural work dominated life
after the circle of the year changed within two centenaries to an energy consumption orientated
urban living environment. The comprehension and mode of life in today's industrial nations
lead to the depletion of fossil fuels, which are the dominating transportable and combusted
form of energy. However, the depletion of oil fields and coal beds signifies a notable enrichment
in atmospheric carbon dioxide which remains one of the prime reasons for global warming.
Figure 1.1:
Saline aquifer scaled down to the dimen
sions of the CCSShowcase with a fluo
rescently dyed plume spreading
Carbon dioxide Capture and Storage abbrevi
ated as CCS is a possibility to encounter the
global energetic mismatch and to enable de
veloping countries the access to sufficient en
ergy for the need and purpose of their citi
zens. In the process of CCS carbon dioxide
is captured as waste product of fossil fuel ox
idation through adsorption in energy plants
and stored in deep and shallow saline aquifers.
These saline aquifers consist of porous and
permeable lime or sandstone filled with saline
ambient water and are isolated from envi
ronmental soil layers through an imperme
able caprock. Carbon dioxide is fluidized and
pumped as supercritical CO
2
into the aquifer
where it is embedded by different mechanism
like structural, residual, solubility and finally
mineralization trapping. Due to public`s lack
of understanding and thus the general disacceptance of CCS a commercialization of this tech
nology cannot be realized. Hence, a showcase for demonstration and teaching purpose is built
in the Focusproject "Carbon Dioxide Capture and Storage Showcase". The figure 1.1 shows
a plume of fluorescently dyed ethyl acetate spreading towards the top of the water and glass
beads filled showcase. Experiments with ethyl acetate and heptane are undertaken to find a
fluid with comparable spreading properties like supercritical CO
2
. The purpose of this study is
to complement the experimental analysis by simulations of ethyl acetate and heptane through
implementing gravity effects into a pressure solver in §2, which is afterwards validated in a series
of Test Cases in §3. To understand the fluid migration a parameter study is done and analyzed
after influencing factors in §4. In chapter §5 we are finally able to simulate and predict the fluid
migration of heptane and ethyl acetate by using the pore network model Hypon.
5
2 Pore network solver with gravity
Porous media like sandstone or limestone, which can be found in a saline aquifer for carbon
dioxide storage or packed glass beads in the CCSShowcase are simulated on a microscale using
different computational techniques which are selected through a compromise between informa
tive value and computational demand.
In this sense a more computational expensive possibility in multiphase simulation is the Lattice
Boltzmann method, which solves the NavierStokes equation in an arbitrary pore geometry and
topology, but is consequently limited in computational time and network size. On the other hand,
stochastic rulebased models are developed and used for certain spreading pattern like Invasion
Percolation for capillary fingering, Diffusionlimited aggregation (DLA) for viscous fingering and
AntiDiffusionlimited Aggregation (AntiDLA) for a stable displacement of the injected fluid.
Lenormand et al. [1988] showed in a comparison with experiments that these methods work in
their limiting cases very well and even quickly, but are insufficient in a general twophase flow
situation where both viscous and capillary forces are present.
A possible compromise in this question is a porescale network solver, which is able to solve
the pressure field based on the simplified momentum equation in form of the HagenPoiseuille
flow on a geometrically simplified pore network out of pore bodies and pore throats. Quasistatic
porenetwork solver calculate the pressure field by incrementally increasing the capillary pres
sure so that the invading fluid is driven from one to another equlibrium state in a capillary
dominated pattern without viscous effects. Due to this lack a dynamic pore network solver is
firstly advanced by Koplik and Lasseter [1985] which calculates the pressure distribution in both
fluids using a saturation dependent local capillary pressure. Both quasistatic and dynamic pore
network algorithms provide sufficient information on macroscopic transport properties like rela
tive permeability, average saturation and average capillary pressure. Provided that the network
morphology is appropriately captured the simulations are able to predict the spreading behavior
of the injected fluid which is realized by implementing gravity as external force into the pressure
solver.
The following dynamic pore network is based on the work of Thompson [2002] and Koplik and
Lasseter [1985], which is improved due to stability criteria in the local capillary pressure by the
work of Lee et al. [2009] and expanded through including gravity effects into the flow behavior.
2.1 Physical displacement mechanism
Lenormand et al. [1983] analyzed in their paper the basic mechanism of immiscible displacement
of incompressible fluids in an etched grid. Due to simplifications of the pore morphology the
fundamental results of their work are used to model the physical principles behind the displace
ment pattern in a pore network.
In a first ascertainment, the displacement is divided into drainage, where the nonwetting fluid is
6
pushing the wetting fluid, and imbibition, where the reverse action happens.
Figure 2.1: (a) Cubical pore bodies i with nonwetting phase and j with wetting phase connected
over a pore throat ij with square crosssection M . (b) Square crosssection M with
central nonwetting phase flow and wetting phase corner flow
Drainage injection of a nonwetting fluid into a pore structure means that the fluid has a higher
nonwetting pressure and is thus able to displace the wetting bulk, which is called Meniscus
movement or Piston displacement. Regarding a single pore i (figure 2.1a) the pressure
difference, which is resulting from this observation should be high enough to overcome the entry
pressure
p
c
e,ij
=
nw
r
ij
+ cos
2

4
 sin cos
cos 
4
 + sin cos
(2.1)
of the pore throat ij, which is a function of the inscribed pore throat radius r
ij
, the contact
angle and the interfacial tension
nw
. The equation (2.1) is a geometrical fit to the wellknown
YoungLaplace equation p
c
=
2
nw
r
cos for noncircular pore throats and is used by JoekarNiasar
et al. [2010] for their pore network solver as well. The displaced wetting phase can be pressed
in another pore body j or into the corners of pore throat ij, resulting in a wetting corner flow
(crosssection view M in figure 2.1b). Pore throats with circular crosssection are only developing
a thin wetting film along the walls. Both mechanism result from the wettability property of
the wetting ambient fluid towards the rigid network, which is assumed to be perfect in our case
(contact angle
0).
Figure 2.2: Snapoff effect in a pore throat: (a) separation process between two pore bodies
(Lenormand et al.); (b) crosssectional view M on the instable nonwetting flow with
swelling wetting phase corners
7
Accumulation of wetting phase in a pore neck or in certain pore bodies causes a capillary
pressure difference due to a gradient in interfacial curvature after the YoungLaplace equation,
which is followed by an instable interface (figure 2.2b) of the nonwetting bulk and consequently
drives spontaneously accumulated corner wetting phase back into the center of the pore neck.
The outcome of this Snapoff effect is that clusters of nonwetting phase are separated from
the invading nonwetting bulk (figure 2.2a) and thus trapped in the porous media. The snapoff
effect is strongly linked to the aspect ratio of the pore body to the pore throat size and occurs
if the capillary pressure falls under a the certain value of
p
c
ij
p
c
snap,ij
=
nw
r
ij
(cos  sin )
(2.2)
(JoekarNiasar et al. [2010], Rossen [2003]'s work on steadystate foam generation, Kovscek and
Radke [1996] on snapoff in constricted noncircular capillaries). Moreover, each phase flows
simply in the domain already occupied by itself, called Bulk flow.
Imbibition as reinvasion of wetting fluid is controlled by the pore body radii instead of the pore
throat radii. The flow mechanism apart from secondary effects like corner and thin film flow
are the same as in the drainage process. However, the meniscus movement is now an instability
conducted rapid Retraction of the nonwetting phase.
In all considerations, a more physical representation of sinusoidal pore throats, allowing realistic
Haines jumps by slowing the nonwetting down and reaccelerating it, is like in Lenormand et al.
[1983] basic work neglected. The reason is simply based on the computational intensity (see
Reeves and Celia [1996] pore network). Furthermore, the saturation dependent local capillary
pressure p
c
i
inside a pore body i is derived by a filling procedure out of one pore throat instead of
multiple ones. Although Lenormand et al. [1983] presents an empirical expression for a filling out
of multiple pore throats, called "Cooperative filling" , the local capillary pressure p
c
i
in pore
body i is derived by a filling procedure out of one pore. While the dynamic capillary pressure
is changing during the discharge of nonwetting phase, entry and snapoff pressure are both fixed
flow conditions for imbibition and drainage, depending on pore network geometry and interfacial
tension
nw
between the nonwetting and wetting fluid.
2.2 Local capillary pressure
p
c
i
(s
w
i
) (Lee et al. [2009])
The original YoungLaplace equation (2.3) describes the pressure difference between the wetting
and nonwetting fluid interface in the case of a static mechanical equilibrium, where it relates the
capillary pressure p
c
to the interfacial tension
nw
, the contact angle and the spherical shape
curvature r. Figure 2.3a and 2.3b explain the idea behind equation (2.3) on the basis of a stress
balance.
p
c
= p
nw
 p
w
=
2
nw
r
cos
(2.3)
The strong dependance on interfacial curvature is recognized during the exemplary nonwetting
phase invasion of pore body i, which is resulting in an altering local capillary pressure curve.
The invasion of pore body i is coupled with a change of the nonwetting (s
nw
i
) and wetting phase
saturation (s
w
i
), which is respectively defined as the occupied fraction of pore body i.
8
Figure 2.3: (a) Spherical interface between two immiscible fluids in a cylindrical tube [Tchelepi
and Horne, 2007]; (b) Visualization of a stress balance on the fluid interface [Tchelepi
and Horne, 2007]; (c) Nonwetting invasion of a pore body: (1) entry capillary pres
sure; (2) transition capillary pressure (poresidetouch); (3) capillary pressure by
pinned wetting corners; [Lee et al., 2009] ; (d) resulting local capillary pressure curve
The following computation method for a capillary pressure p
c
i
(s
w
i
) depending on the local wetting
phase saturation is based on the work of Lee et al. [2009] and visualized in figure 2.3c and 2.3d.
1. When entering the pore body i the local capillary pressure will be equal to the entry
capillary pressure of pore throat ij (figure 2.3c&d  1st period).
p
c
i
(s
w
i
= 1) = p
c
e,ij
(2.4)
2. The fluid interface is expanding with an increase in curvature while invading pore body i,
which is resulting in a decreasing local capillary pressure until the nonwetting front touches
the pore body walls. This configuration with maximum inscribed sphere radius r
max
is the
infimum of the local capillary pressure curve and thus characterized by the transition
saturation s
w,tr
i
= 1 
4
3V
i
r
3
max
with V
i
the total volume of pore body i . The capillary
pressure for a growing nonwetting bubble is given by equation (2.5) and can be derived
using a volume based instead of curvature based YoungLaplace Equation (figure 2.3c&d 
2nd period).
p
c
i
(s
w
i
) = 2
nw
cos (
3
4
V
i
(1  s
w
i
))

1
3
f or
s
w
i
s
w,tr
i
(2.5)
3. After surpassing the transition saturation s
w,tr
i
the fluid interfaces are pinned into the
corners and the interfacial curvature is rising again. An approximate functional form
for the capillary pressure is given by Lee et al. [2009] assuming that the wetting phase
volume becomes proportional to r
3
. Approaching the limit of zero wetting saturation and
thus infinite local capillary pressure evokes numerical instabilities which are damped by a
relaxation parameter in the code (figure 2.3c&d  3rd period).
p
c
i
(s
w
i
) = 2
nw
cos (a V
i
s
w
i
)

1
3
with
a =
3
4
1  s
w,tr
i
s
w,tr
i
f or
0 < s
w
i
s
w,tr
i
(2.6)
Furthermore, when the pore throat is invaded by both fluids, the capillary pressure in a
pore throat p
c
ij
equals the capillary pressure of the upstream pore body.
9
The stress balance for a static equilibrium is used to relate wetting and nonwetting phase pressure
over the dynamic capillary pressure p
c
i
. This should be understood in a discrete manner of
successive equilibrium states.
p
c
i
(s
w
i
) = p
nw
 p
w
(2.7)
2.3 Gravity affected flow rate between pore bodies
The driving local capillary pressure p
c
i
(s
w
i
) as well as the snapoff effect p
c
snap,ij
and the threshold
entry pressure p
c
e,ij
are principal ideas behind fluid displacement and hence for the development
of the pore network solver algorithm. JoekarNiasar et al. [2010] mention in their work that
"no gravity effect has been considered in [their] simulations . . . " and that "adding gravity does
not constitute any major complication in the code". In this subsection we discuss the way of
including gravity effects in our equations by employing a different approach than in common
codes developed for simulating drying (Huinink et al. [2002], Segura and Toledo [2005]), gas
drive (Bondino et al. [2007], Ezeuko et al. [2010]) and gravity included Invasion percolation
algorithms (Wilkinson [1984]).
· Gravity weighted snapoff and threshold pressure
Wilkinson's algorithm is the first porescale code (1984) for invasion percolation which
considers buoyancy effects by using a height dependent entry potential instead of the entry
threshold pressure. Ezeuko et al. [2010] adapted this idea to their dynamic porenetwork
solver and defined the entry condition p
min
c
= min
i
(
2
nw
r
ij
 gH
i
), where is the density
difference of nonwetting and wetting phase, g the gravitational constant and H
i
the height
from the top of the network to the pore body i. A similar term can be matched to our
equations to decrease the entry pressure from upper pore throats so that a vertical invasion
is preferred. The local capillary pressure from previous subsection §2.2 can be used without
modification due to a sufficient change of porous media properties.
· Gravity weighted dynamic capillary pressure
Another possibility is to alter the capillary pressure by adding a height and direction
dependent gravity term. Here, the change directly influences the driving force of the
nonwetting displacement. Consequently, there should be no need to modify entry and
threshold pressure.
The hydrostatic pressure profile should be flexible and selfgenerated without defining a reference
height in this underlying porenetwork solver. Therefore the second option is chosen and the
gravity term is included into the direction dependant flow rate
Q
ij
= K
ij
(p
ij

g l
ij
sin(
ij
)) with p
ij
= p
i
 p
j
and
= nw, w
(2.8)
between two pore bodies i and j, where l
ij
is the pore throat length,
ij
the angle to the horizontal
plain and K
ij
the conductivity between pore bodies i and j . The flux Q
ij
is correlated with
the nonwetting and wetting phase saturation and therefore changes indirectly the local capillary
pressure p
c
i
.
Q
ij
s
i
s
w
i
p
c
i
(s
w
i
)
10
Figure 2.4: Two pore bodies i and j filled with a single wetting phase and connected over the
pore throat l
ij
The figure 2.4 consists of two pore bodies i and j connected with a pore throat of magnitude
l
ij
and filled with a single wetting phase. The upper pore body pressure p
w
j
is set zero while the
lower pore body pressure p
w
i
is unknown. We are assuming that there is no flux Q
w
ij
and thus
the conductivity K
w
ij
is constant. Solving equation (2.8) after p
w
i
is resulting in the expected
hydrostatic pressure of a magnitude of
w
gl
ij
and is therefore confirming our implementation
idea in a first step. Indeed, an algorithmic validation will show rigorously the confirmation of
our assumption.
2.4 Mathematical Equations and algorithm
Before capturing the pore network of the experimental CCSShowcase, the gravity need to be
implemented and validated in the already existing pressure solver after JoekarNiasar et al. [2010]
and Thompson [2002].
The network utilized in their work is a structured isotropic regular lattice, which consists of
cubical volume carrying unresistant pore bodies and cuboid volumeless resistant pore throats.
This assumption allows less computational and memory demanding simulations and is legitimate
because of an exact fluid tracking in pore throats is not necessary. The size of each pore body is
specified after a truncated lognormal distribution, which is related by geometrical considerations
to the pore throat distribution [JoekarNiasar and Majid Hassanizadeh, 2011]. A presentation
of an exemplary two dimensional pore field with a regular coordination number of 4 (number of
adjacent pore bodies) is shown in figure 2.5. A detailed analysis and capture of the pore network
properties of the showcase will be presented in chapter §5 and should be focused on the goal of
predictive simulations of the fluid displacement.
Based on the physical displacement mechanism and the angular network the assumptions for
the pressure solver are summarized as follow.
2.4.1 Assumptions
1. The fluids are immiscible, incompressible and Newtonian, separated by a sharp interface
guaranteeing no mixing or diffusion.
2. Flux in the pore throats are laminar and thus modeled as HagenPoiseuille flow.
3. While the volumeless pore throats are modeled as resistors between two pore bodies, all
11
Excerpt out of 54 pages
 Quote paper
 Bilal Özkan Lafci (Author), 2012, Implementation and Validation of Gravity Effects in a PoreNetwork Solver, Munich, GRIN Verlag, https://www.grin.com/document/276760
Publish now  it's free
Comments