# Direct Numerical Simulation of Transport Phenomena on Pore-space Images

Urban and Civil

Nov 16, 2013 (4 years and 5 months ago)

107 views

Peyman Mostaghimi, Martin Blunt, Branko Bijeljic

11
th

January 2010, Pore
-
scale project meeting

Direct Numerical Simulation of Transport
Phenomena on Pore
-
space Images

1

Flow at pore scale

In petroleum science and engineering, scales of interest may vary from
molecular level to a mega level.

The pore scale is of the order of a typical pore which is in the range few
microns. Modelling fluid flow at the pore scale can provide a predictive tool for
estimating rock and flow properties at larger scales.

One of the most used ways to capture the

morphology of a porous medium as the

main input for pore scale modelling is

micro
-
CT imaging.

Motivation

Network

modelling

the

representation

of

the

pore

space

by

an

equivalent

representation

of

pores

and

throats

has

been

successful
:

we

now

understand

trends

in

recovery

with

wettability

and

can

predict

single

and

multi
-
phase

properties
.

However
……
.
the

extraction

of

networks

involves

ambiguities

and

there

are

some

cases

where

the

method

does

not

work

so

well
.

Now

have

direct

three
-
dimensional

imaging

of

pore

spaces
.

Why

not

simulate

multiphase

flow

directly

on

these

images?

Micro
-
CT imaging AND DIRECT SIMULATION

Post

processing

Micro
-
CT

images,

a

matrix

can

be

generated

for

a

core

which

shows

whether

there

is

a

solid

inside

the

voxel

or

a

pore
.

Zero

means

that

voxel

is

a

pore

and

one

means

it

is

a

solid

phase
.

Two

methods

to

simulate

fluid

flow

in

porous

media

directly

without

the

need

for

simplified

geometries
:

-

the

lattice

Boltzmann

method

(Edo)

-

conventional

computational

fluid

dynamics

algorithms

based

on

the

relevant

flow

and

conservation

equations
.

Governing equations

Conservation of mass:

Navier
-
Stokes equation:

-
state and incompressible flow:

0

v

t
v
v
v
v
2
.



P
t
0

v
v
v
v
2
.



P
Dimensionless analysis

0
~
v
v
v

0
0
0
/
~
l
v
P
P
P

0
~
l
2
2
0
2
~

l

0
0
Re
v
l

v
P
v
v
~
~
~
~
)
~
~
.
~
Re(
2

1
10
/
.
10
/
10
/
10
10
Re
5
2
3
3
3
6
5


m
s
N
m
kg
s
m
m
-
state Navier
-
Stokes
equation:

Reynolds number for flow in porous media:

v
2
0



P
Stokes Equation :

FORMULATION

0
1
1
1
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2

z
w
y
v
x
u
z
P
z
w
y
w
x
w
y
P
z
v
y
v
x
v
x
P
z
u
y
u
x
u

u
Equation

(u: velocity in x
-
direction)

v
Equation

(v: velocity in y
-
direction)

w
Equation

(w: velocity in z
-
direction)

p Equation

Gridding

Marker
-
and
-
cell grid:

Existence of solid phase in each grid causes six velocity components be
zero in the 3 dimensional models.

Discretized Form

x
P
P
z
u
u
u
y
u
u
u
x
u
u
u
K
J
I
K
J
I
K
J
i
K
J
i
K
J
i
K
J
i
K
J
i
K
J
i
K
J
i
K
J
i
K
J
i

,
,
1
,
,
2
1
,
,
,
,
1
,
,
2
,
1
,
,
,
,
1
,
2
,
,
1
,
,
,
,
1
1
)
(
2
)
(
2
)
(
2

y
P
P
z
v
v
v
y
v
v
v
x
v
v
v
K
J
I
K
J
I
K
j
I
K
j
I
K
j
I
K
j
I
K
j
I
K
j
I
K
j
I
K
j
I
K
j
I

,
1
,
,
,
2
1
,
,
,
,
1
,
,
2
,
1
,
,
,
,
1
,
2
,
,
1
,
,
,
,
1
1
)
(
2
)
(
2
)
(
2

z
P
P
z
w
w
w
y
w
w
w
x
w
w
w
K
J
I
K
J
I
k
J
I
k
J
I
k
J
I
k
J
I
k
J
I
k
J
I
k
J
I
k
J
I
k
J
I

1
,
,
,
,
2
1
,
,
,
,
1
,
,
2
,
1
,
,
,
,
1
,
2
,
,
1
,
,
,
,
1
1
)
(
2
)
(
2
)
(
2

0
,
,
1
,
,
,
,
,
1
,
,
,
,
,
1

z
w
w
y
v
v
x
u
u
k
J
I
k
J
I
K
j
I
K
j
I
K
J
i
K
J
i

i
K
J
I
K
J
I
nb
nb
K
J
i
K
J
i
A
p
p
u
a
u
a
,
,
,
,
1
,
,
,
,

j
K
J
I
K
J
I
nb
nb
K
j
I
K
j
I
A
p
p
v
a
v
a
,
,
,
1
,
,
,
,
,

k
K
J
I
K
J
I
nb
nb
k
J
I
k
J
I
A
p
p
w
a
w
a
,
,
1
,
,
,
,
,
,

The momentum equations can be rewritten as:

SIMPLE algorithm

The SIMPLE (Semi
-
Implicit Method for pressure
-
Algorithm:

1. Guess the pressure field p*

2. Solve the momentum equations to obtain u*,v*,w* by algebraic
mutigrid

solver

3. Solve the p’ equation (The pressure
-
correction equation) by algebraic
mutigrid

solver

4. p=p*+p’

5. Calculate u, v, w from their starred values using the

velocity
-
correction equations

6. Solve the
discretization

equation for other variables, such as

temperature, concentration, and turbulence quantities.

7. Treat the corrected pressure p as a new guessed pressure p*,

converged solution is obtained.

Storing matrices in CRS and AMG for solving all linear systems of equations.

Boundary Condition

0

p
u
0

N
u

i
K
J
I
K
J
I
nb
nb
K
J
i
K
J
i
A
p
p
u
a
u
a
,
,
,
,
1
,
,
,
,

j
K
J
I
K
J
I
nb
nb
K
j
I
K
j
I
A
p
p
v
a
v
a
,
,
,
1
,
,
,
,
,

k
K
J
I
K
J
I
nb
nb
k
J
I
k
J
I
A
p
p
w
a
w
a
,
,
1
,
,
,
,
,
,

. . . Boundary Condition

)
(
8
)
(
2
3
2
2
2
y
O
y
y
u
y
y
u
u
u
p
n

)
(
2
)
(
3
2
2
2
y
O
y
y
u
y
y
u
u
u
p
s

2
2
2
)
(
3
4
4
3
8
y
u
u
u
y
u
s
p
n

COMPARISON OF THE Three METHODS FOR BC FOR
FLOW BETWEEN TWO INFINITE PARALEL PLATES

first method

second method

third method

VELOCITY PROFILE FOR FLOW BETWEEN TWO INFINITE PARALEL PLATES

We see non
-
zero velocity even for one block within the channel and for more
than one we see agreement to within machine accuracy with the analytical
solution.

Lid
-
driven Cavity

DISPERSION MODELLING

When a miscible fluid is injected in a flowing fluid in a saturated porous
diffusion.

In brief, dispersion is the spread or mixing of flowing fluids due to all these
mechanisms.

To model advection term we use stream tracing algorithm and for diffusion
we apply random walk method.

Diffusion
X
X
t
x
t
t
x

)
(
)
(

STREAMLINE TRACING

Interpolation to estimate the velocity vectors
at a point within the grid block

K
J
i
i
K
J
i
K
J
i
u
x
x
x
u
u
u
,
,
,
,
,
,
1
)
(

K
j
I
j
K
j
I
K
j
I
v
y
y
y
v
v
v
,
,
,
,
,
1
,
)
(

k
J
I
k
k
J
I
k
J
I
w
z
z
z
w
w
w
,
,
,
,
1
,
,
)
(

)
ln(
,
,
,
,
1
i
e
K
J
i
K
J
i
x
u
u
u
u
x

)
ln(
,
,
,
1
,
i
e
K
j
I
K
j
I
y
v
v
v
v
y

)
ln(
1
,
,
1
,
,
i
e
k
J
I
k
J
I
z
w
w
w
w
z

)
,
,
min(
z
y
x

)
1
.(
,
,
,
,
1
,
,
,
,
1
0

K
J
i
K
J
i
u
u
x
i
K
J
i
K
J
i
e
u
u
u
x
x
x
)
1
.(
,
1
,
,
1
,
,
1
,
,
1
,
0

K
j
I
K
j
I
v
v
y
i
K
j
I
K
j
I
e
v
v
v
y
y
y
)
1
.(
1
,
,
1
,
,
1
,
,
1
,
,
0

k
J
I
k
J
I
w
w
z
i
k
J
I
k
J
I
e
w
w
w
z
z
z
The time

of flight:

The coordinates of exit location
:

Diffusion

t
D
m

2

Cos
Sin
X
X
.
.
0

Sin
Sin
Y
Y
.
.
0

Cos
Z
Z
.
0

Random walking method for both advection and diffusion:

t
D
t
x
V
t
x
t
t
x
m

2
)
(
)
(
)
(
Random walking method just for diffusion part of

flow

:

diffusion

t
D
X
t
x
t
t
x
m

2
)
(
)
(
. . . Diffusion

Gridding:

Resolution:

0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
1
2
3
4
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
1
2
3
4
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
m

851
.
8
50
50

m

4255
.
4

m

702
.
17

0

3
.
6

Pe
4
.
0

Pe
m
avg
D
L
u
Pe

t
D
m

2

. . . Particle tracking

Gridding:

Resolution:

Sandpack

LV60B

338
338

m

851
.
8
. . . Particle tracking

Gridding:

Resolution:

Sandpack LV60B

338
338

m

851
.
8
. . . Particle tracking

Gridding:

Resolution:

Sandpack

LV60B

338
338

m

851
.
8
DISPERSION COEFFIEICENT

The average of positions of particles:

Variance of X can be calculated:

And the longitudinal dispersion coefficient:

For showing the importance of diffusion, dispersion is modelled for a range
Peclet number:

N
t
X
t
X
X
N
i
i
i



1
)
(
)
(
N
t
Y
t
Y
Y
N
i
i
i



1
)
(
)
(
N
t
Z
t
Z
Z
N
i
i
i



1
)
(
)
(
N
X
t
X
X
t
X
N
i
i
i




1
2
2
2
)
)
(
(
)
)
(
(

dt
d
D
L
2

m
avg
D
L
u
Pe

(Bijeljic et al. 2004)

Multiphase flow at the pore scale

Having the interface at different saturations, the flow of each phase can be
modelled by the Stokes solver and the relative permeability can be
predicted.

Also for reactive transport (Branko), the code can be used to simulate the
flow at each time step.

Courtesy of
Masa

Prodanovic