one page summary paper

Dr.NoorH

Huang_mines_0052N_10933.pdfsequence1.pdf

Home >English homework help >Article writing homework help >one page summary paper

PORE-SCALE SIMULATION OF MULTIPHASE FLOW USING LATTICE BOLTZMANN

METHOD: DEVELOPMENTS AND APPLICATIONS

Jingwei Huang

A thesis submitted to the Faculty and the Board of Trustees of the Colorado School of Mines

in partial fulfillment of the requirements for the degree of Master of Science (Petroleum

Engineering).

Golden, Colorado Date ________________ Signed: _______________________ Jingwei Huang Signed: _______________________ Dr. Xiaolong Yin Thesis Advisor Golden, Colorado Date ________________ Signed: _____________________ Dr. Erdal Ozkan Professor and Head Department of Petroleum

iii

ABSTRACT

Multiphase flow is an important transport process that is found in many natural phenomena

and engineering applications. In petroleum engineering, sub-surface multiphase flow takes place

throughout the primary, secondary, and tertiary stages of recovery. Advances in high-

performance computing and experimental techniques now allow us to study multiphase flow at

the pore level to better understand pore-scale dynamic and its continuum-level descriptions. In

this work, a parallelized multiphase flow simulator is developed using the lattice Boltzmann

method with the color-gradient multiphase model. Several boundary conditions for the

multiphase flow model are developed, and they allow the simulator to handle a range of

practically relevant multiphase flows that are pressure-driven, rate-driven, or gravity-driven. The

simulator is parallelized, which is critical for simulating flows in computational domains of

meaningful sizes. The parallel performance is good in large computation domains. Validation

tests show that the simulator is able to generate controlled interfacial tension, wettability and can

accurately simulate dynamic problems. As an application, we simulated pressure-driven

displacements in a porous medium made up by packed spheres using fluids with different

viscosity ratios and wettability at different capillary numbers. Our results are consistent with

published numerical and experimental results.

TABLE OF CONTENTS

ABSTRACT ................................................................................................................................... iii

LIST OF FIGURES ....................................................................................................................... vi

LIST OF TABLES ....................................................................................................................... viii

ACKNOWLEDGEMENTS ........................................................................................................... ix

CHAPTER 1 INTRODUCTION .................................................................................................... 1

CHAPTER 2 LITERATURE REVIEW ......................................................................................... 5

2.1 Single-phase LB method and boundary condition ................................................................ 5

2.2 Two-phase LB method .......................................................................................................... 7

2.3 Boundary conditions for two-phase LB method ................................................................. 10

CHAPTER 3 SINGLE- AND MULTIPHASE LATTICE BOLTZMANN METHODS AND PARALLELIZATION .............................................................................................................. 12

3.1 Single-phase LB method and boundary conditions ............................................................. 12

3.1.1 Single-phase LB method .............................................................................................. 12

3.1.2 Pressure and Velocity boundary conditions for single-phase LB method .................... 16

3.2 Multiphase LB model .......................................................................................................... 21

3.3 Parallelization of color-gradient LB code ........................................................................... 24

3.4 Test cases ............................................................................................................................. 31

CHAPTER 4 BOUNDARY CONDITIONS FOR MULTIPHASE LATTICE BOLTZMANN . 39

4.1 Evaluation of existing outflow boundary conditions for two-phase LB method ................ 39

4.2 Modified periodic boundary condition for two-phase LB method ..................................... 43

4.3 Pressure boundary condition for two-phase LB method ..................................................... 45

CHAPTER 5 SIMULATION OF IMMISCIBLE DISPLACEMENTS IN A POROUS MEDIUM .................................................................................................................................. 59

CHAPTER 6 CONCLUSIONS AND FUTURE WORK ............................................................. 74

6.1 Summary and Conclusions .................................................................................................. 74

6.2 Future work ......................................................................................................................... 76

REFERENCES CITED ................................................................................................................. 78

LIST OF FIGURES

Figure 3.1 D3Q19 propagation scheme. ....................................................................................... 13

Figure 3.2 Velocity directions at boundary. .................................................................................. 17

Figure 3.3 Scheme of decomposition: (a) Computation domain before decomposition; (b) Computation domains after decomposition. ................................................................ 26

Figure 3.4 Porous medium used for parallel performance test. .................................................... 29

Figure 3.5 Speedup of multiphase lattice Boltzmann code........................................................... 30

Figure 3.6 Parallel efficiency of multiphase lattice Boltzmann code. .......................................... 31

Figure 3.7 Oscillation and relaxation of a square drop to equilibrium spherical shape. ............... 33

Figure 3.8 Numerical measurement of Laplace law. .................................................................... 35

Figure 3.9 Thickness of interface using color-gradient lattice Boltzmann model. ....................... 36

Figure 3.10 Force balance at the contact line. .............................................................................. 37

Figure 3.11 Contact angles achieved by setting � and � on wall nodes. .............................. 38 Figure 3.12 Droplet deformation in Stokes flow. ......................................................................... 40

Figure 3.13 Deformation parameter as a function of ; at =0.1. ...................................... 41

Figure 4.1 Schematic of outlet boundary. ..................................................................................... 40

Figure 4.2 Phase distribution of a moving droplet in a channel using the convective outflow boundary condition. .................................................................................................... 42

Figure 4.3 Phase distribution of two-phase displacement in a channel using the convective outflow boundary condition. ....................................................................................... 43

Figure 4.4 Phase distribution of a moving droplet in a channel using the modified periodic boundary condition. .................................................................................................... 44

Figure 4.5 Phase distribution of two-phase displacement in a channel using modified periodic boundary condition. .................................................................................................... 45

Figure 4.6 Schematic of inlet and outlet boundaries..................................................................... 46

Figure 4.7 Phase distribution (a) and corresponding dimensionless pressure field (b) of a moving droplet in a channel using the pressure boundary condition. The dimensionless pressure is defined as +⁄ . .............................................. 51 Figure 4.8 Phase distribution (a) and corresponding pressure field (b) of two-phase displacement in a channel using the pressure boundary condition. dimensionless pressure is defined as +⁄ . ..................................................................... 53 Figure 4.9 Simulation setup for capillary filling. The capillary tube has a length L and diameter d. The red fluid is the invading fluid and the blue fluid is the defending fluid. ......... 55

vii

Figure 4.10 The length of tube occupied by the filling fluid as a function of time ∗ = / for (a) drainage, M = 1 (b) imbibition M = 1 (c) drainage M = 0.1 and (d) drainage M = 10. The lines represent the analytical solution based on Eqn. (4.22) and the crosses mark the position (1, 1). ............................................................................... 56

Figure 4.11 Pressure distributions along the flow direction for (a) drainage, M = 1 (b) imbibition M = 1 (c) drainage M = 0.1 and (d) drainage M = 10 at ∗ = / = 0.55, 0.55, 0.50, 0.50 respectively. P* is the dimensionless pressure defined as +⁄ . ...................................................................................................... 57 Figure 5.1 The cylindrical porous medium made up by packed spheres. ..................................... 60

Figure 5.2 Two-phase distribution (a) and wetting-phase pressure distribution (b) in the sphere packing at different times (∗= 0.12, 0.47, 0.71, 1.18), Case A. ................................. 63 Figure 5.3 Two-phase distribution (a) and wetting-phase pressure distribution (b) in the sphere packing at different times (∗= 0.12, 0.47, 0.71, 1.18), Case B. ................................. 64 Figure 5.4 Two-phase distribution (a) and wetting-phase pressure distribution (b) in the sphere packing at different times (∗= 0.12, 0.47, 0.71, 1.18), Case C. ................................. 65 Figure 5.5 Two-phase distribution (a) and wetting-phase pressure distribution (b) in the sphere packing at different times (∗= 0.12, 0.47, 0.71, 1.18), Case D. ................................. 66 Figure 5.6 Saturation distributions at different times in Case A, B, C, and D. Time interval ∆ ∗ between consecutive lines is 0.059. The red line represents the saturation distribution when breakthrough occurs. ...................................................................... 67

Figure 5.7 Phase pressure distribution of Case A, B, C, and D. The open symbols represent pressure profiles at times before breakthrough: ∗ = 0.89, 0.47, 0.30, 0.12 respectively. The lines represent the pressure profiles after breakthrough: ∗ = 1.18 for Case C and D. ........................................................................................................ 68

Figure 5.8 Relative permeability vs. wetting phase saturation of the sphere packing. Red and green points denote breakthrough points. Time interval ∆ ∗ between consecutive points is 0.059. ............................................................................................................. 69

viii

LIST OF TABLES

Table 3.1 Stationary bubble test. �� = � = , �� = � = , = .......................................... 34 Table 3.2 Stationary bubble test. �� = � = , �� = � = , = . ....................................... 34 Table 4.1 Comparison of breakthrough times in the case of capillary filling ............................... 58

Table 5.1 Parameters in different cases of immiscible displacement ........................................... 61

ACKNOWLEDGEMENTS

First of all, I would like to express my gratitude to my advisor, Dr. Xiaolong Yin, for his

guidance, support and help for this work. I also want thank my committee members, Dr. Hossein

Kazemi and Dr. Azra Tutuncu for their valuable comments. I appreciate the help from Dr. Philip

Winterfeld for the help to use the EMG cluster. My appreciation also goes to Dr. Timothy Kaiser,

Director of Research and High-Performance Computing at Colorado School of Mines, for his

help on code testing on BlueM.

I thank my mother for her love, encouragement and support all the time. I want to extend my

appreciations to all my friends at Colorado School of Mines.

Financial supports from the Unconventional Natural Gas and Oil Institute, Coupled Integrated

Multiscale Measurement and Modeling consortium, and the Department of Energy through the

Nuclear Energy University Program DE-NE0000719 are gratefully acknowledged.

CHAPTER 1

INTRODUCTION

In petroleum engineering, multiphase flow through porous media is a key process that

controls production throughout the primary, secondary, and tertiary stages of a reservoir.

Understanding the dynamics of multiphase flow is essential for one to predict and control

reservoir production. A large amount of theoretical, experimental, and numerical studies have

already been dedicated to this area.

On the macroscopic level, Darcy’s law is used to describe fluid flow through porous media. It

is known that Darcy’s law originated from averaging of the Stokes flow (Whitaker 1998).

Description of multiphase flow on the continuum level, on the other hand, usually uses a

multiphase extension of Darcy’s law, which is empirical and does not have a strong theoretical

foundation (Hassanizadeh and Gray 1979). In this extension, relative permeability is introduced

to account for the interaction of the two-phases in the pore space. However, it is still a challenge

to model relative permeability as it is affected by various factors including pore structure,

saturation, viscosity ratio, capillary number, wettability, as well as the way that flow is generated

and maintained in porous medium. Traditionally, relative permeability is determined by

experiment, which is time consuming and often difficult to interpret. It is more difficult to

conduct multiphase flow experiments in tight porous media because of their low permeability

(Karimi et al. 2015). Therefore, it is of great interest to study multiphase flow from the level of

pores with numerical simulations so that one can attain better understanding of pore-scale

dynamics and improve continuum-level descriptions.

Direct numerical simulation of multiphase flow through porous media is an important

approach to study complex multiphase flow processes, such as those in multiphase reactors and

heat/mass exchangers (Kreutzer et al. 2005; Dudukovic et al. 1999), microfluidics (Xu et al.

2014; Wu et al. 2012; Wu et al. 2013), oil and gas production (Orr and Taber 1984; Cueto-

Felgueroso and Juanes 2013), contaminant transport and remediation (Esposito and Thomson

1999; Mackay and Cherry 1989; Hunt et al. 1988), and geological storage of greenhouse gases

(Szulczewski et al. 2012). Compared to continuum-scale simulations and pore network models,

pore-scale direct numerical simulation directly solves the dynamics of fluid motion in pore space,

and therefore is a more fundamental approach to understand the driving and controlling

mechanisms of multiphase flow such as viscous coupling, capillary pressure and the gravity

force (Joekar-Niasar and Hassanizadeh 2012).

Many approaches have been developed to directly simulate multiphase flow in porous media.

These approaches include the phase field (PF) method (Cahn and Hilliard 1958; Cahn and

Hilliard 1959; Badalassi et al. 2003; Lowengrub and Truskinovsky 1998), the volume-of-fl uid

(VOF) method (Hirt and Nichols 1981; Rider and Kothe 1998; Pilliod and Puckett 2004), the

level-set (LS) method (Osher and Sethian 1988; Sussman et al. 1998), the density functional

method (DFM) (Demianov et al. 2011) and the lattice Boltzmann (LB) method (Frisch et al.

1987; Chen et al. 1991; d'Humières 1992; Gunstensen et al. 1991; Gunstensen and Rothman

1993; Grunau et al. 1993; Shan and Chen 1993; Shan and Doolen 1995; Swift et al. 1995). In the

PF model, the interface is described as a thin transition layer where unstable mixtures are

stabilized by non-local energy terms (Cahn and Hilliard 1959; Lowengrub and Truskinovsky

1998). A very thin interface layer is required to improve the accuracy of simulation which may

result in expensive computation cost (Badalassi et al. 2003). The VOF algorithms consist of three

parts: interface reconstruction, advection and surface tension (Gerlach et al. 2006). The VOF

method, as its name implies, tracks the volume of each fluid in the cells that belong to the

interface (Gerlach et al. 2006). The LS method defines the interface by a level-set function

(Osher and Sethian 1988) and has been widely used due to its simplicity (Osher and Fedkiw

2001). Re-initialization techniques are employed to overcome interface smearing and enforce

mass conservation (Sussman et al. 1994). However, VOF and LS method may face numerical

instability when simulating flows in porous media where capillary forces generally dominate

over viscous forces (Scardovelli and Zaleski 1999). In LB methods, phase segregation is

modeled by inter-particle interactions and interface tracking is not needed (He et al. 1999; Reis

and Phillips 2007). Similar to PF, interfaces are diffuse and much thicker than real interfaces in

LB methods. However, for flows where interfacial thickness is not crucial, owning to LB

methods’ simplicity and high computation efficiency (Liu et al. 2012), they have become

attractive numerical options for direct simulation of multiphase flow (Kang et al. 2004; Kang et

al. 2005; Pan et al. 2004; Li et al. 2005; Aidun and Clausen 2010; Chen et al. 2014).

The objective of this study is to develop a multiphase lattice Boltzmann model and a number

of flow driving conditions to facilitate simulations of multiphase flow in porous media. A

parallel, color-gradient based multiphase lattice Boltzmann model is built and several flow

driving conditions have been implemented and compared. With the developed simulator, we

conducted simulations of immiscible displacement in a complex geometry made up by packed

spheres.

The thesis is organized as follows: Chapter 1 introduces motivations and objectives of this

study. Chapter 2 reviews previous work on single-phase lattice Boltzmann, multiphase lattice

Boltzmann and relevant boundary conditions. Chapter 3 first describes the lattice Boltzmann

model applied in this study and the validation tests. Then, the parallel algorithm and parallel

performance of the simulator are presented. Chapter 4 studies the performance of several

boundary conditions for multiphase lattice Boltzmann. In Chapter 5, immiscible displacements

with different viscosity ratios, wettability, and capillary numbers are simulated in a porous

medium made up by packed spheres. Conclusions and recommended future work are presented

in Chapter 6.

CHAPTER 2

LITERATURE REVIEW

In this section, we summarize previous studies on single-phase and multiphase lattice

Boltzmann models and applications. Different boundary conditions for lattice Boltzmann are also

reviewed. Through this review, we show that developing proper boundary conditions for

multiphase lattice Boltzmann method is critical for simulation of many practically relevant flows.

2.1 Single-phase LB method and boundary condition

Traditional computation fluid dynamics (CFD) methods directly solve the Navier-Stokes

equations that are derived from conservation of mass and momentum on the continuum level.

Molecular dynamics, on the other hand, focuses on the motion of individual molecules; they can

determine velocities and locations of molecules at every instant of time. Lattice Boltzmann

method, as a mesoscopic method, deals with the velocity distribution function of molecules. LB

method, compared with traditional CFD, has the advantage that it does not require complex

discretization of partial differential operators; compared with molecular dynamics, it does not

need to focus on the details of every molecule, which significantly reduces computational cost.

With these said, it is also important to keep LB method’s limitations in mind: due to its lattice

nature and connection to the kinetic theory of fluids, it is not as flexible as traditional CFD when

it comes to complex fluids (e.g. non-Newtonian) and complex boundary conditions. Because LB

method still relies on continuous descriptions of molecular velocity distributions, it should not be

applied to flows with significant molecular granularity and fluctuations.

Lattice Boltzmann method was originated from lattice gas (Frisch et al. 1987) and gained a lot

of attention after it was invented. LB method models time evolution of a discrete molecular

velocity distribution function through propagation and collision steps. Propagation moves the

distribution function on a space-filling lattice to simulate the streaming motion of molecules.

Collision, on the other hand, simulates the viscous relaxation of the non-equilibrium part of the

velocity distribution. Many different collision schemes have been developed. Single-Relaxation-

Time scheme (Chen et al. 1991; Qian et al. 1992), which is based on the BGK collision operator

in Boltzmann equation (Bhatnagar et al. 1954), is very simple but is more liable to numerical

instability in complex flows. Multiple-Relaxation-Time (MRT) scheme, first proposed by

d'Humières (1992), applies different relaxation times to different kinetic modes and is able to

achieve a better stability and accuracy with a moderate increase in computational cost (Pan et al.

2006). Other collision models, for example, Two-Relaxation-Time scheme (Ladd 1994) and

Generalized MRT collision scheme (Shan et al. 2007), have their own characteristics but are

similar to MRT in terms of ideas and performance.

Single-phase lattice Boltzmann methods have been successfully applied to simulate flow and

transport in porous media. As in all computational methods, boundary conditions and flow

driving conditions are very important. In most LB applications, a uniform external body force is

applied to drive fluid flow in a computational domain bounded with the periodic boundary

condition. In such systems, the effects of a uniform body force and a pressure gradient are

essentially equal. However, pressure-driven flows are not always equal to body-force-driven

flows. It has been argued that they are equal only in flows with a constant cross-section.

Moreover, the accuracy of body force-approach in a complex geometry like porous media has

not been validated (Zhang and Kwok 2006). Hence, developing pressure boundary condition or

velocity boundary condition is essential for extending capabilities of LB method.

Boundary conditions for solid walls (no-slip boundary condition) of LB method are related to

those from lattice gas automata, e.g., particle bounce-back scheme. Chen et al. (1996) pointed

out that a simple extrapolation of the distributions on the boundary can achieve second-order

accuracy. Maier et al. (1996) also proposed a similar extrapolation method. However, these

methods need distribution’s orientation to the wall and require special treatments for corner

nodes, which cannot be easily extended to general complex geometries (Zou and He 1997). As to

pressure and velocity boundaries that often serve as inlet / outlet conditions to drive fluid flow,

Zou and He (1997) built a pressure and velocity boundary condition by applying a bounce-back

rule for the non-equilibrium distributions on the boundary nodes. At inlet and outlet, zero

transverse velocities were assumed. Zhang and Kwok (2006) proposed a modified periodic

boundary condition to implement a pressure differential. Their model has superior performance

near inlet and outlet but requires periodicity. Kim and Pitsch (2007) proposed yet another

boundary closure scheme for both compressible and incompressible flows and their model shows

a higher accuracy than Zhang and Kwok (2006).

All of the above models have their respective advantages and weaknesses. It is still an active

topic to design simple, steady and accurate boundary conditions for lattice Boltzmann

simulations.

2.2 Two-phase LB method

A number of multiphase LB models have been developed for simulation of multiphase flow,

including the color gradient model (Gunstensen et al. 1991; Gunstensen and Rothman 1993;

Grunau et al. 1993; Liu et al. 2012), the pseudo-potential model (Shan and Chen 1993; Shan and

Doolen 1995), the free energy model (Swift et al. 1995) and the mean-field model (He et al.

1999). The first model, the color gradient model (Gunstensen et al. 1991), was proposed based

on a two-component lattice gas model (Rothman and Keller 1988). Fluid is colored red or blue to

distinguish the two phases. A two-phase perturbation operator is applied to generate interfacial

tension and a recoloring step forces the particle to move to the fluids with the same color.

Variations of density and viscosity, though of limited ranges, were implemented in a subsequent

study (Grunau et al. 1993). To reduce the computation cost of the recoloring step, d’Ortona et al.

(1995) derived a simplified recoloring method. Latva-Kakko and Rothman (2005) proposed a

new algorithm for the recoloring step to reduce lattice pinning by allowing moderate phase

mixing at tangent of interface. In 2007, a new perturbation operator that correctly recovers

Navier-Stokes equations for two-phase flows and can tackle different viscosities and density

ratios was developed (Reis and Phillips 2007). In their work, an explicit expression for the

interfacial tension was derived. Later, Latva-Kakko and Rothman’s recoloring operator was

adopted to Reis and Phillips’ model, which led to reduced spurious currents in two-dimensional

computations (Leclaire et al. 2011). TÖlke et al. (2002) built a three-dimensional 19-velocity

model for immiscible binary fluids with different densities and viscosity ratios. However, this 3D

model does not have an explicit expression for the interfacial tension. Recently, Reis and Phillips’

two-dimensional model was extended to three-dimensional by Liu et al. (2012). Latva-Kakko’s

recoloring operator was also implemented and an explicit expression for the 3D interfacial

tension is available in Liu et al.’s work (2012). In the color gradient model, the procedure to

separate different phases is time-consuming and high spurious velocities can occur near the

interface (Chen and Doolen 1998). One advantage of the model is that the interfacial tension,

contact angle, viscosity ratio and density ratio can all be adjusted independently (Ahrenholz et al.

2008).

The second multiphase LBM model is the pseudopotential model that directly describes

interactions between fluid particles using potentials (Shan and Chen 1993). Shan and Doolen

(1995) then modified the basic model by redefining velocity and significantly reduced errors in

recovering the macroscopic equation. However, the basic pseudo-potential model suffers from

some limitations such as large spurious currents, limited range of density and viscosity ratios,

lack of ability to adjust properties independently, and thermodynamic inconsistency (Chen et al.

2014). A lot of work has been done to improve the pseudopotential model. To name a few, it is

shown that interaction forces that have a higher order of isotropy can suppress spurious currents

(Shan 2006; Sbragaglia et al. 2007). Incorporating a realistic equation of state (EOS) can reduce

spurious currents and increase the range of density ratios (Yuan and Schaefer 2006; Zhang and

Tian 2008). Proper forcing scheme also helps to suppress spurious currents, improve

thermodynamic consistency and increase stability of pseudo-potential models (Kupersthtokh et al.

2009; Guo et al. 2002). Recently, a critical review (Chen et al. 2014) was conducted, which

systematically listed the limitations and the improvements of pseudopotential models.

The third model, the free energy model, was developed to ensure thermodynamic consistency

by introducing a non-ideal pressure tensor directly into the collision operator (Swift et al. 1995).

Although the original free energy model achieves mass and momentum conservation, it lacks

Galilean invariance, despite that this problem can be partially suppressed by adjusting terms in

the pressure tensor (Swift et al. 1996). In 2000, a Galilean invariant free-energy model was

developed by Inamuro et al. (2000) using asymptotic analysis. Kalarakis et al. (2002) also

proposed a two-dimensional free-energy model that recovers Navier-Stokes equations and

preserves Galilean invariance. In the literature, comparisons of the above three multiphase LB

models have been carried out and it was shown that the color gradient model and the free energy

model have similar capabilities (Huang et al. 2011; Yang and Boek 2013).

2.3 Boundary conditions for two-phase LB method

For multiphase LB methods, schemes to recover pressure or velocity boundary conditions are

more difficult to develop than for single-phase LB. In most publications of multiphase LB, a

periodic boundary condition is applied and flow is driven by a body force. Such a set up cannot

simulate multiphase flow in experimental systems that are driven by either pressure or a fixed

injection rate, with non-periodic boundaries. It is therefore important to derive appropriate

boundary conditions and apply them to properly simulate multiphase flow. Some studies have

been undertaken, for example, the stress-free boundary condition (Liu et al. 2012). This method

can ensure mass conservation and improve numerical stability for creeping flows (Liu et al.

2012). Neumann boundary condition has also been formulated and tried in simulations of

immiscible displacement (Dong et al. 2010). A systemic study on outflow boundary conditions

performed by Lou et al. (2013) showed that the convective boundary condition behaved better

than Neumann and extrapolation boundary conditions. Recently, Huang et al. (2014) conducted a

study of immiscible displacement using a velocity boundary condition at inlet and a constant-

pressure boundary condition at outlet based on the non-equilibrium bounce-back method. Later,

in a simulation of liquid CO2 displacing water, constant-flow-rate and constant-pressure

boundary conditions were set at the inlet and outlet, respectively (Liu et al. 2014). With these

developments, multiphase flows that are driven by the rate of injection can now be properly

simulated.

The above brief review shows that currently there is no established method to simulate

multiphase flow driven by controlled inlet and outlet pressure differentials using LB model. It is

therefore important to develop a pressure boundary condition to extend the capability of

multiphase LB methods. In most multiphase LB studies, a uniform body force is used instead of

a pressure difference. However, a pressure differential cannot always be replaced by a body force.

For single-phase flow, Zhang and Kwok (2006) showed that flows driven by pressure differential

and body force are only equal in single-phase flow with constant cross-sections. Kandhai et al.

(1999) also compared body force and pressure boundary approaches in single-phase flow. Their

results show that for flows with small Reynolds numbers and through simple geometries, body

force is equivalent to pressure boundary. However, for flows with high Reynold numbers or

through complex geometries, a pressure boundary cannot always be replaced by a body force. In

view of the importance of the pressure boundary and that few studies have been undertaken to

implement pressure differentials for multiphase LB, in this work, we will present a method to

implement pressure boundary conditions for multiphase lattice Boltzmann simulations.

CHAPTER 3

SINGLE- AND MULTIPHASE LATTICE BOLTZMANN METHODS AND

PARALLELIZATION

The section introduces lattice Boltzmann method for single-phase flow and that for

multiphase flow based on the color gradient model. The parallel performance of the multiphase

lattice Boltzmann code is tested on a super-computer. The simulator is then verified by several

validation tests.

3.1 Single-phase LB method and boundary conditions

Single-phase lattice Boltzmann method has been widely applied to simulate fluid flow and is

the fundamental for multiphase lattice Boltzmann. Here we briefly introduce the single-phase LB

model as well as several boundary conditions.

3.1.1 Single-phase LB method

The Lattice Boltzmann equation can be regarded as a discretization of the continuous

Boltzmann equation:

+ �.∇� + ∇� = Ω (3.1) where is the velocity distribution function, � is particle position. � is particle velocity, is acceleration, Ω is the collision operator between particles. The lattice Boltzmann method has two processes, particle propagation and particle collision. Here, particle is defined as a collection of

molecules on a certain lattice node. The governing transport equation for LB method is

� + ∆t, + ∆t = �, + Ω �, (3.2)

where �, is the particle velocity distribution function, � is the position vector of the fluid node, defined on a cubic lattice, ∆t is the time step, is the discretized velocity possessed by molecules in , Ω represents the collision operator in LB method, �, + Ω �, represents post-collision populations that propagate to location � + ∆t in next time step. In this study, the D3Q19 scheme is used as shown in Figure 3.1. D3Q19 stands for three

dimensions and nineteen discretized velocities, as listed below

= { , , � =± , , , ,± , , , ,± � = ~± ,± , , ± , ,± , ,± ,± � = ~ (3.3)

Figure 3.1 D3Q19 propagation scheme.

Particle collision as denoted by Ω simulates the effect of viscosity, and is very important for LB method. Ω must satisfy mass and momentum conservation. Several collision schemes have been

developed and here we apply the Multiple-Relaxation-Time (MRT) (d'Humieres 2002) collision

method. MRT introduces multiple relaxation parameters to control the rates of relaxation of

various moments of the particle velocity distribution. MRT has better stability and accuracy

compared to the simpler BGK collision.

When applying the MRT scheme, a matrix is introduced into the governing equation

� + ∆t, + ∆t = �, − �− ��[ �, − ] (3.4) = �[ + .� + .� − � ] (3.5)

where � is a transformation matrix (19 × 19) for the D3Q19 scheme (d'Humières 2002). �− is the inverse of �. � is the collision matrix in the moment space. The diagonal collision matrix � in the moment space is defined as (d'Humières 2002)

� ≡ dia� , , , , , , , , , , , , , , , , , , (3.6) where are the relaxation parameters chosen as = = �⁄ , = − / − , == = �⁄ , = . represents the equilibrium distribution. The difference between

and is the non-equilibrium distribution . is the lattice sound speed (√⁄ ) and is a weighting coefficient associated with discretized velocities:

= { / � =/ � = ~/ � = ~ (3.7) Fluid density � and fluid velocity � can be obtained from the first and second moments of the distribution function

� = ∑ = ∑ (3.8) �� = ∑ = ∑ (3.9)

The viscosity in lattice unit is correlated to the relaxation time τ (one of the parameters in the collision matrix �) by

= � τ − . (3.10) where � is the static fluid density. The pressure is given as

= � (3.11) After the collision step, the post-collision populations will either stay at same node or

propagate to neighboring nodes based on Eqn. (3.2). If neighboring node is solid, a boundary

condition must be applied. Here, we use the simple link-bounce-back rule that recovers the no-

slip boundary condition (Ladd 1994). In this model, fluid-solid interface is defined as a

collection of nodes whose positions are halfway along the link that connects fluid and solid

nodes. The populations of fluids that propagate towards solid nodes will be reflected back to the

direction where they came. For a moving interface, the link-bounce-back rule is (Ladd 1994):

− �, + ∆t = �, − � .� � , (3.12) where �, is the incoming population, − �, + ∆t is the bounce-back population in the opposite direction of �, , � = � + ∆t is the location of the interface node, and � is the velocity of the interface.

The mass, time and length units in LB simulation are related to the static fluid density � and the temporal and spatial lattice units that are designated as ∆t and ∆x, respectively. A unit conversion is required to relate simulation units to the actual dimensions. For example, the

resolution of an image of porous medium is usually the real-world equivalence of lattice

Boltzmann ∆x (always 1). By comparing the equilibrium density of the lattice fluid to the actual

density, we can establish conversion in the mass unit. Then, we can take the relation between the

lattice viscosity and the actual viscosity given by the equation below to determine ∆ ( ∆� ∆ ) = ( ∆� ∆ )�� (3.13)

After ∆t and ∆x are determined, other quantities in LB simulation can be converted by using dimensionless numbers. For example, lattice velocity, pressure, and permeability can be

converted to real units by

( ∆∆ ) = ( ∆∆ )�� (3.14) (� ∆∆ ) = (� ∆∆ ∆ )�� (3.15) ( �∆ ) = ( �∆ )�� (3.16)

where is the velocity, is the body force in simulation, ∆ ⁄ is the actual pressure gradient, � is the permeability.

3.1.2 Pressure and Velocity boundary conditions for single-phase LB method

Several boundary conditions for single-phase LB have been proposed. Here we choose to

present the pressure and velocity boundary condition derived by Zou and He (1996), based on

which we developed a pressure boundary condition for multiphase LB. In the original paper, the

method was presented for D2Q9 and D3Q15 LB models. Here, we extend their derivations to the

D3Q19 model.

To implement the pressure boundary, non-equilibrium distributions that are bounced back on

boundary nodes need to be modeled. Take the case of a channel flow as shown in Figure 3.2 for

an example: flow is in the x direction and solid boundaries are applied in the y and z directions.

At the inlet, after propagation, (i = 1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 14, 15, 18, 19) are known; (i = 2, 12, 13, 16, 17) are unknown distributions that need to be determined. Assuming that == on the boundary nodes and the density at boundary is specified as � , can be determined from the two equations below:

Figure 3.2 Velocity directions at boundary.

∑ = �= (3.17) � = + + + + − + + + + (3.18)

Then an explicit expression of is obtained

= � − [ + + + + + + + + + + + + + ]� (3.19) By applying the bounce-back rule for the non-equilibriums of the unknown (i = 2, 12, 13, 16,

17), we have

= + ( − ) (3.20) where is the distributions in the opposite direction of . Using Eqns. (3.5) and (3.20), all

unknown at the boundary nodes can be obtained. In order to ensure that y- and z-momentum

are zero, (i = 2, 12, 13, 16, 17) are modified (Maier et al. 1996)

= + ( + ) i = 2, 12, 13, 16, 17 (3.21) Here we choose:

α = / = − + − + − = − + − + −

(3.22)

This modification makes the x-momentum unchanged while y- and z-momentum zero. Finally,

the unknown distributions at inlet are:

= + � = + � − − + − + − = + � + − + − + − = + � − − + − + − = + � + − + − + −

(3.23)

The same procedure can be applied to the outlet nodes. Density at outlet boundary is specified as � . The unknown distributions at the outlet are: = − �

= − � − − + − + −

= − � + − + − + − = − � − − + − + − = − � + − + − + −

(3.24)

Similarly, we can get the unknown distributions when flow is in y or z direction.

If flow is in the y direction, the unknown distributions at the inlet are:

= + � = + � − − + − + − = + � + − + − + − = + � − − + − + − = + � + − + − + −

(3.25)

The unknown distributions at the outlet are:

= − � = − � − − + − + − = − � + − + − + − = − � − − + − + − = − � + − + − + −

(3.26)

If flow is in the z direction, the unknown distributions at the inlet are:

= + � = + � − − + − + − = + � + − + − + − = + � − − + − + − = + � + − + − + −

(3.27)

The unknown distributions at the outlet are:

= − � = − � − − + − + − = − � + − + − + −

= − � − − + − + − = − � + − + − + −

(3.28)

There are no special treatments for nodes at the corners.

The velocity boundary can also be derived similarly. Take the flow in the x direction as an

example. As is known, from Eqns. (3.17) and (3.18), the expression for density at the inlet is

� = + + + + + + + + + + + + +− (3.29) Then, by following the same procedure, the unknown distributions at the inlet can be obtained.

3.2 Multiphase LB model

In this work, we use a three-dimensional color gradient LB model. As in all multiphase LB

models, hydrodynamics is recovered from the evolution of discrete lattice velocity distribution

functions , where i denotes the index of the discrete velocity distributions, and superscript k

represents fluid components. We assign the values of R and B to k to represent the red and blue

components (� and ). In the color gradient model, for every time step, each component undergoes a collision operation and a propagation operation, just like in single-phase LB model

� + ∆ , + ∆ = �, + Ω �, (3.30) Here, �, is the velocity distribution function of component k at position � and time t. is the lattice velocity carried by �, , ∆ is time step. The lattice velocity � is the D3Q19 model

= { , , � =± , , , ,± , , , ,± � = ~± ,± , , ± , ,± , ,± ,± � = ~ (3.31) Ω is the collision operator that includes three parts (TÖlke et al. 2002) Ω = (Ω ) [(Ω ) + (Ω ) ] (3.32)

(Ω ) is the single-component collision operator. When BGK collision operator is used, it is: (Ω ) = −� ( − ) (3.33)

where is the equilibrium distribution function of , � is the dimensionless relaxation time of component k. Mass and momentum are the first two moments of the distribution

functions:

� = ∑ = ∑ (3.34)

�� = ∑∑ � =∑∑ � (3.35) where � is the local fluid velocity, � is the density of component k, � = �� + � is the total density. In this study, both components are assumed to have equal density. Then, the equilibrium

distribution functions are given by:

= � ( + [ �.� + �.� − � ]) (3.36) where the weight coefficient are

= { / � =/ � = ~/ � = ~ (3.37) The inter-component collision operator (Ω ) is a perturbation operator defined as (Liu et al. 2012)

(Ω ) = |∇� |[ .∇�|∇� | − ] (3.38) where the phase field � is expressed by (Liu et al. 2012)

� �, = �� , − � �,�� , + � �, (3.39) We choose = − / , − = / , − = / as in the work of Liu et al. (2012). Using these parameters, the theoretical expression of the interfacial tension is (Liu et al. 2012)

� = � + � (3.40) The interfacial tension is therefore controlled by three parameters �, and relaxation time �. In this study, we let � = = . The pressure can be obtained via the equation of state and given as

= � = � (3.41)

The relaxation time � is related to the kinematic viscosity of the fluid = (� − ) (3.42) To accommodate large viscosity ratios, an interpolation scheme is used to determine the

relaxation time � within an interface (Grunau 1993). � = {

��, � > � , � > , � − � , � < − (3.43) where is a free positive parameter that controls the zone of application of the relaxation

time model. � = + + , = + + , and = =�� / �� + � , = �� − / , = − / , = − � / , = / . Because the multi-relaxation-time (MRT) scheme can enhance numerical instability

(d'Humières 2002), suppress spurious velocities (Liu et al. 2014; Zhao and Fan 2010) and

minimize the viscosity-dependence of permeability (Li et al. 2005), we use MRT collision

operator instead of BGK in this study. As in the work of Liu et al. (2014), the total distribution

function ( = � + ) is used in the single-component collision and inter-component collision steps. With MRT, (Ω ) and (Ω ) become

(Ω ) = − �− �� [ − ] (3.44) (Ω ) = |∇� |[ .∇�|∇� | − ] (3.45)

where � is an integer transformation tensor (d'Humières 2002). The diagonal collision matrix � in the moment space is defined as (d'Humières 2002)

� ≡ dia� , , , , , , , , , , , , , , , , , , (3.46)

where are the relaxation parameters chosen as = = �⁄ , = − / − , == = �⁄ , = . The recoloring operator Ω� forces the separation of the two components. A widely used recoloring operator is the one proposed by Latva-Kokko and Rothman (2005):

Ω� � = �� + �� , (3.47) Ω� = �� − �� , (3.48) is a free parameter between 0 and 1 that controls interface thickness. � is the angle between the phase field gradient ∇� and . To implement a contact angle, it is assumed that the solid wall is a mixture of two fluids with constant proportions, which means that the solid wall has a

certain value of the phase field � . For instance, if the contact angle of blue phase is , the density of red and blue components on a solid wall are set as � − cos and � + cos , respectively (Zhang et al. 2015).

3.3 Parallelization of color-gradient LB code

In practical applications, simulating three-dimensional flows in a complex geometry often

needs millions to billions of lattices and very high computational cost. A parallelized two-phase

code is therefore of significant value. Parallelization of the color gradient model is based on the

existing single-phase code (Xiao and Yin 2015). Part of the parallel schemes were inherited from

the single-phase code. Extensions to the parallel schemes were built according to the procedures

of the color gradient LB model. In what follows, we present the main steps of the parallelization

algorithm.

1. Parallelism of the single-phase LB code

(1) Build geometry

Read geometry from either experimentally obtained images of porous media or images

generated by numerical reconstruction. Process the geometry into a three-dimensional array that

contains (0, 1) (fluid = 1, solid = 0).

(2) Pre-processing

To parallelize the computation, the computation domain (the three-dimensional array) is

decomposed into a number of subdomains based on the number of processors to be used: 1, 2,

4, …, 2n, with n being the number of cores. Before flow simulation, the following pre-processing

steps are needed.

a) Decompose the computation domain as shown in Figure 3.3 and construct a ghost layer

(green) for each subdomain. Ghost layer contains nodes that are shared between cores and

they are used for communication between cores.

(a)

(b)

Figure 3.3 Scheme of decomposition: (a) Computation domain before decomposition; (b) Computation domains after decomposition.

b) Each core identifies the fluid nodes in its own ghost layer.

c) Each core packs its ghost nodes into an array that contains the coordinates of the nodes and

19 indices to ghost nodes’ neighbors.

d) Each core sends the list of ghost nodes to its neighbors.

e) With ghost layers built, every core constructs its own local lattice Boltzmann node map.

(3) Single-phase LB flow simulation

a) Each core reads input parameters and the local node map.

b) Each core performs collision on every node in the local node map.

c) Each core performs propagation. During a propagation step, if the neighbor is local,

propagation is performed within the local subdomain. If the neighbor is in the ghost layer,

communication is needed. The distribution to be propagation therefore enters a specific

package based on its destination. MPI (Message Passing Interface) is used to communicate

these packages when all of them are ready.

d) A master process collects the information from every core to calculate the average

properties, such as saturation and superficial velocity.

e) Back to b) till all time steps are completed.

2. Parallelism of the color-gradient LB code

The parallel implementation of the two-phase code is an extension of the parallel

implementation of the single-phase code. In the two-phase code, there are two components. Each

component has its own velocity distribution for collision and propagation.

(1) Build geometry

Same as the single-phase code.

(2) Pre-processing

Decompose the computational domain and identify fluid or solid in the same way as the

single-phase code.

(3) Two-phase LB flow simulation

a) Each core reads input parameters and the local node map.

b) Each core calculates the local density of each component and transfer the density of the

fluids nodes that on the sub-domain surface to a density package for communication.

c) All cores exchange density packages by MPI and store received packages to the local ghost

layer such that every fluid node in the local computational domain can access all of its

neighbors’ densities.

d) Perform two-phase collision steps shown in Eqn. (3.38).

e) Propagate the distributions of each component like the single-phase code.

f) A master process calculates the average properties of the entire computation domain.

The performance of the parallelized two-phase LB simulator was tested by simulating

immiscible displacements in a porous medium (Figure (3.4)). Three different size of geometries

were used that contain 6403, 3203 and 1603 nodes respectively.

Figure 3.4 Porous medium used for parallel performance test.

Tests were performed on BlueM, which is a unique high-performance computing system at

Colorado School of Mines (built by IBM) that contains two partitions. One is named Mc2

(Energy) which runs on an IBM BlueGene Q and another is named Aun (Golden) that uses the

iDataplex architecture. We only used Mc2 in the performance test for consistency. Mc2 contains

512 nodes and has 8192 Gbytes of total memory. Each node has 16 cores (PowerPC A2 17 Core).

The performance of the two-phase code is evaluated according to the following benchmarks:

(1) Speedup

Speedup SP is defined as: SP = TTN (3.49) where T is the simulation time using one core, TN is the simulation time using N cores. The speedup of different geometries using different number of cores are shown in Figure 3.5.

Figure 3.5 Speedup of multiphase lattice Boltzmann code.

When comparing speedup in the different geometries, it can be observed that speedup is higher

in the larger system. For larger systems, more computational work is assigned to each core. The

fraction of communication overhead therefore is reduced, leading to a higher speedup and a

higher parallel efficiency as we shall observe below. Though not tested, we expect that the

speedup will continue to increase with increasing system size.

(2) Parallel Efficiency

Parallel efficiency is defined as the ratio of the speed up to the number of cores:

EN = SpN (3.50) From Figure 3.6, we can see that parallel efficiency of the two-phase code decreases with

increasing number of cores used, owning to the increasing ratio of computational load assigned

to each core versus the communication overhead.

Figure 3.6 Parallel efficiency of multiphase lattice Boltzmann code.

As the size of geometry increases, the decreasing rate of parallel efficiency slows down.

To summarize, a parallel LB simulator for two-phase flow in porous media was built based on

the existing single-phase LB simulator. Parallel efficiency is satisfactory. As the size of the

geometry increases, the speedup increases and the parallel efficiency is higher. In these tests,

decomposition was done by consecutive bisections in three directions. However, computational

cost in each subdomain is controlled by the actual number of fluid nodes rather than the

subdomain size. Thus, load balance is not equal, which is an issue that negatively affects the

parallel performance especially for flows through highly heterogeneous porous media geometries.

3.4 Test cases

Several test cases are performed in this part to verify the developed two-phase color gradient

LB simulator. When a drop of one fluid is immersed into another, surface tension will make the

drop relax toward its equilibrium shape – a sphere or a circle. As shown in Figure 3.7, an initial

square red drop ( × ) is centered in a × computation domain, surrounded by the blue fluid. The parameters in this test are �� = � = , �� = � = , � = . × − . The initial configuration is shown in the top left and different snapshots are taken at different

simulation times. As the snapshots show, the square drop deforms gradually under the influence

of the interfacial tension and becomes circular in the end.

Another test case characterizes the capillary pressure within a stationary bubble. A spherical

red bubble is initialized at the center of a × × computation domain, surrounded by the blue fluid. According to the Young-Laplace equation, the pressure difference ∆ across the interface of a spherical bubble is related to the interfacial tension by:

∆ = � (3.51)

where is the radius of bubble and � is the interfacial tension. In our simulations, different pressure differences can be obtained by changing the bubble radius.

Figure 3.7 Oscillation and relaxation of a square drop to equilibrium spherical shape.

As exemplified in Figure 3.8, the slope of ∆ vs. ⁄ is the simulated interfacial tension. Theoretically, the interfacial tension of the color model is given as � = � + �. Without

losing generality, in this study, it is always assumed that � = = . Table 3.1 and Table 3.2 give a comparison between simulated and theoretical results. It can be seen that simulated

interfacial tensions agree well with theoretical results for both = and = . cases. The errors for the interfacial tension, calculated by = |� ℎ − � | � ℎ⁄ × %), increase with A. = generates more accurate interfacial tensions than = . . In many multiphase LB models, spurious current, which refers to non-zero, non-physical fluid velocities

in the vicinity of an interface, is an inevitable problem (Yuan and Schaefer 2006). In this model,

a lower spurious current is obtained (|�| is the maximum spurious current) when = . . At the same time, however, affects the thickness of interface. Figure 3.9 gives the density

distribution along a line that passes through the center of the bubble. When = . , the interfacial thickness is six lattice units. When = . , the thickness decreases to five lattice units. A lower therefore increases the thickness of the interface and suppresses the spurious

current, however at the expense of the accuracy of the interfacial tension.

Figure 3.8 Numerical measurement of Laplace law.

Table 3.1 Stationary bubble test. �� = � = , �� = � = , = � � ℎ % |�|

1.00E-04 4.510E-05 4.44E-05 1.47% 2.79E-05

1.00E-03 4.515E-04 4.44E-04 1.59% 5.08E-05

1.00E-02 4.520E-03 4.44E-03 1.71% 2.02E-04

1.00E-01 4.537E-02 4.44E-02 2.08% 1.75E-03

Table 3.2 Stationary bubble test. �� = � = , �� = � = , = . � � ℎ % |�|

1.00E-04 4.579E-05 4.44E-05 3.03% 2.27E-06

1.00E-03 4.579E-04 4.44E-04 3.03% 1.76E-05

1.00E-02 4.581E-03 4.44E-03 3.06% 1.00E-04

1.00E-01 4.597E-02 4.44E-02 3.43% 9.76E-04

Figure 3.9 Thickness of interface using color-gradient lattice Boltzmann model.

One important ability of the color gradient LB method is that it can produce controlled static

contact angle. As shown in Figure 3.10, the contact angle is determined by the force balance at

the three-phase contact line.

Figure 3.10 Force balance at the contact line.

In the color gradient model, contact angle is achieved by setting � and � on the wall nodes. The contact angle can be analytically determined by (Latva-Kakko and Rothman 2005):

= (� − �� ) (3.52) where � is the initial density. Figure 3.11 shows the different contact angles achieved by this approach. A full range of contact angles can be achieved from wetting to non-wetting, although

those at the extremes are difficult to simulate accurately, owning to the finite thickness of the

fluid-fluid interface.

Figure 3.11 Contact angles achieved by setting � and � on wall nodes.

To examine whether this multiphase LB model can properly simulate the dynamics of two-

phase flow in the absence of solid walls, a Taylor deformation test was performed. As shown in

Figure 3.12, when a spherical drop is placed in the middle of a linear shear flow, viscous effect

tends to deform the drop, and interfacial tension acts to maintain drop’s spherical shape. In this

test, two parallel plates move in the opposite direction to generate a linear shear flow, in the

middle of which a spherical drop is placed. The deformation of the drop is described as a

function of the Reynolds number and the capillary number. The definitions of the Reynold

number and the capillary number are (Liu et al. 2012)

= �, = � (3.53) where = ⁄ is shear flow rate. is the half of the distance between two plates. is the velocity of moving plate. R is the initial radius of droplet which is 10 in lattice unit. = ρν is the dynamic viscosity. For this case, the two fluid has same density and viscosity. At the steady state,

the deformation of droplet is described using a deformation parameter defined as

= −+ (3.54) where a and b are the lengths of major and minor axes, respectively. When flow is Stokes,

should obey the Taylor relation when is low (Taylor 1932)

= ⁄ (3.55) Simulations were performed at different capillary numbers when = . in a × × system. The effect of was also tested.

Figure 3.12 Droplet deformation in Stokes flow.

Figure 3.13 plots the deformation parameter as a function of capillary number. It is can be

seen that the simulation results are in reasonable agreement with the analytical solution. The

agreement is better when = . .

Figure 3.13 Deformation parameter as a function of ; at = 0.1.

CHAPTER 4

BOUNDARY CONDITIONS FOR MULTIPHASE LATTICE BOLTZMANN

For lattice Boltzmann methods, it is not straightforward to implement open boundary

conditions because, during propagation, distributions from the outside of the computational

domain must be speculated. For multiphase LB method, boundary conditions are not well

developed, because the distributions of both components, and sometimes their interactions, must

be modeled at the inlet and outlet boundaries. In this section, we study several boundary

conditions for multiphase LB method. First, we evaluate the performance of an outflow

boundary and a modified periodic boundary for the color gradient lattice Boltzmann model. Then,

a pressure boundary condition that can be applied to simulate pressure-driven flows using

multiphase lattice Boltzmann is derived. As validations, pressure-driven flows through a

capillary tube with different density ratios, viscosity ratios and contact angles are tested.

4.1 Evaluation of existing outflow boundary conditions for two-phase LB method

Several outflow boundary conditions have been tested systematically by Lou et al. (2013)

using a mean-field multiphase LB model (He et al. 1999). Here, we implement the convective

boundary condition, which behaves best among the three outflow boundary conditions tested

(Lou et al. 2013), into the color-gradient multiphase LB model. As shown in Figure 4.1, the

outlet boundary is located at = . Nodes at + and + are the ghost nodes needed for two-phase interactions at the boundary. Nodes at − and − are the two nearest layers of real fluid nodes before the outlet boundary.

Figure 4.1 Schematic of outlet boundary.

For the convective boundary condition, distributions at the outlet boundary obey the

following equation

∂∂t + ∂ = (4.1) where is the phase velocity normal to the outlet boundary (Lou et al. 2013). Here, we choose

as the average velocity at layer − . Using first-order finite difference, the convective boundary condition can be written as

� , + ∆ − ,∆ + � , + ∆ − � − , + ∆∆ = , + ∆ − ,∆ + , + ∆ − − , + ∆∆ = (4.2)

In lattice Boltzmann units, ∆ = , ∆ = , we therefore have � , + = � , + � − , ++

, + = , + − , ++

(4.3)

Then the distribution of outlet boundary nodes , + is known. For the two-phase interaction at boundary, distributions at + are also set according to the convective boundary condition

� + , + = � + , + � , ++ + , + = + , + , ++

(4.4)

To evaluate the performance of the convective boundary condition, two tests were conducted

in a rectangular computational domain with length = , width � = and height = . Two pairs of solid walls are placed in the direction of and �, outflow boundaries are applied in the direction of . At the inlet, velocity boundary condition is applied for both phases. In this

simulation, the parameters are �� = � = , �� = � = , � = . × − , contact angle =°. The first test is a moving drop in a channel flow. A red drop initialized with radius = was placed at ( , , = , , and surrounded by another blue fluid. Density and viscosity are the same for the two phases. The blue fluid is injected continuously with a constant

rate. Performance of the outflow boundary condition is evaluated by checking whether the drop

can leave the outlet boundary smoothly without significant disturbance to the upstream flow nor

to the shape of the drop. Figure 4.2 shows the simulated exit of the drop from the outlet using the

outflow boundary condition. When the red drop reaches the outlet, it keeps its shape and move

out of the system freely without any unphysical phenomenon. This result is consistent with the

report of Lou et al. (2013).

In the second test, we simulated the immiscible displacement in the same computational

domain which is an important topic in many practical problems. Initially, section is

occupied by the blue phase and the rest of the channel is occupied by the red phase. Density and

viscosity are set as unity for both fluids. Also, the blue phase is injected by a constant rate from

the inlet. Figure 4.3 shows the time evolution of the interface during the immiscible displacement.

Similar to the moving drop case, the outflow boundary allows the fluids move out of the domain

freely.

Based on the moving drop case and immiscible displacement case, we can conclude that the

convective outflow boundary condition is implemented to the color-gradient lattice Boltzmann

model successfully.

Figure 4.2 Phase distribution of a moving droplet in a channel using the convective outflow boundary condition.

Figure 4.3 Phase distribution of two-phase displacement in a channel using the convective outflow boundary condition.

4.2 Modified periodic boundary condition for two-phase LB method

In addition to the convective outflow boundary condition, we also tested a modified periodic

boundary condition, designed to combine periodic and outflow boundaries.

The idea of this modified periodic boundary condition is, at the outlet, the convective

boundary condition is applied which has performed well in the previous tests. At the inlet, to

keep the periodicity, fluid is allowed to re-enter from the outlet, yet the color is changed to that

of the displacing fluid. In this way, the velocity field is forced to be continuous between the inlet

and the outlet. Then, both fluids are driven by a uniform body force applied to the entire

computation domain. In the simulations, the parameters are set as �� = � = , �� = � = ,

� = . × − , contact angle = °. With this setup, we conducted the same test of a drop moving through a channel. The snapshots in Figure 4.4 show that the drop can also move out of

the outlet with minimum disturbance to the upstream flow field and to the shape of the drop. The

same test of immiscible displacement in a channel is also conducted. As shown in Figure 4.5, the

two phases can also move out of the domain freely.

As a conclusion, the modified periodic boundary condition works well in the color-gradient

lattice Boltzmann model. This boundary condition should be very useful for simulations of

gravity-driven multiphase flow.

Figure 4.4 Phase distribution of a moving droplet in a channel using the modified periodic boundary condition.

Figure 4.5 Phase distribution of two-phase displacement in a channel using modified periodic boundary condition.

4.3 Pressure boundary condition for two-phase LB method

At the end of this chapter, we present the derivation of pressure inlet and outlet boundary

conditions for two-phase LB simulations, developed based on the single-phase pressure

boundary condition proposed by Zou and He (1997). In their work, the pressure boundary is

applied by letting the non-equilibrium part of distributions bounce back from the boundary nodes.

In order to set up the pressure boundary condition for two-phase flow, the non-equilibrium part

of the total distribution function is bounced back. Then, the unknown distributions of each phase

are determined using its own density defined at the boundary.

Let us take the case of channel flow shown in Figure 4.6 as an example. The boundary is

aligned with the x direction. Inlet and outlet boundary is located at = (I for inlet) and = (O for outlet) respectively. Nodes at = − and = + are the ghost nodes that are needed for two-phase interactions at the boundary. Nodes at x = + and x = − are the nearest layers of real fluid nodes to the inlet and outlet. At the inlet boundary = , a constant density � is set.

Figure 4.6 Schematic of inlet and outlet boundaries.

At the inlet, � (i = 2, 12, 13, 16, 17), (i = 2, 12, 13, 16, 17) are the unknown distributions that need to be determined. Since we have

� = ∑∑ (4.5) � = ∑( + + + + ) − ( + + + + ) (4.6)

the inlet velocity can be obtained as

� = � − ∑[ + + + + + + + + + ( + + + + )] (4.7)

Then, we use the bounce-back rule for the none-equilibrium part of the total distribution function

(i = 2, 12, 13, 16, 17):

= + ( − ) (4.8) where is the distribution function at the opposite direction of , = ∑ , is the total equilibrium function. By using the definitions of the equilibrium distribution functions in Eqn.

(3.36), for each component, we have

, − , = � ,� = (4.9) , − , = � ,� = , , , (4.10)

Then by Eqns. (4.8) and (4.10), we can get

= + � ,� = (4.11) = + � ,� = , , , (4.12) In order to make sure that y- and z-momentum are zero, (i = 2, 12, 13, 16, 17) need to be

modified

= + ( + ),� = , , , (4.13) This modification adds , to the y-, z-momenta and keep the x-momentum unchanged. Here

we choose:

α = / = ∑ − + − + − (4.14)

= ∑ − + − + − Finally, the unknown total distributions at the inlet are:

= + � = + � − ∑ − + − + − = + � + ∑ − + − + − = + � − ∑ − + − + − = + � + ∑ − + − + −

(4.15)

Next, we need to divide to obtain the unknown distribution function for each component.

Although the total density at the inlet � is defined, the densities of component R and component B (�� and � ) are still not known. To solve this problem, we divide the total distribution according to �� and � of the previous time step:

� = �� + � (4.16) = �� + � (4.17)

When calculating phase field gradient in the inter-component collision and recoloring step at

boundary nodes, the information of ghost nodes is needed. We assume the distribution of ghost

nodes at = − is the same as that of inlet boundary nodes: − = . Similarly, at the outlet, by setting a constant density � at outlet boundary nodes and following the same procedure, the unknown distributions at the outlet are:

= − � = − � − ∑ − + − + − = − � + ∑ − + − + − = − � − ∑ − + − + − = − � + ∑ − + − + −

(4.18)

These distributions are then divided to obtain � and according to the densities of component R and component B at the outlet (�� and � ) of the previous time step.

� = �� + � (4.19) = �� + � (4.20)

As at the inlet boundary, we assume the distribution of ghost nodes at = + is the same as that of outlet boundary nodes: + = . To check the performance and accuracy of the pressure boundary condition, three test cases

are carried out. The first two cases show the effect of the pressure boundary condition on two-

phase flow. The third case is the simulation of capillary filling, which is a problem with a

moving contact line that admits an analytical solution.

A. Effect of pressure boundary condition on two-phase flows

In this section, we test the pressure boundary condition to see whether it generates two-phase

flow that are physically correct. The first test is performed in a rectangular computational

domain with length = , width � = and height = . Two pairs of solid walls are placed in the direction of and �, pressure boundaries are applied in the direction of . Initially, the domain is filled with the blue phase, except that a red drop with radius = is placed at ( , , = , , . By making � slightly above the equilibrium value of unity and � slightly below, a pressure difference is established between the inlet and outlet, and fluid flow is

generated. The steady-state flow is characterized by the Reynolds number Re and the capillary

number Ca that defined as

R� = ,Ca = � (4.21) where is the velocity of blue phase along the flow direction. = � is the dynamic viscosity of the blue phase. � is the interfacial tension. In this simulation, the parameters are �� = � = , �� = � = , � = . × − , � = . , � = . , R� = . , Ca = . , contact angle = °. Figure 4.7(a) shows the transient motion of the two phases at different simulation times. The

corresponding pressure fields are presented in Figure 4.7(b). As the figures show, the fluids are

moved by the imposed pressure difference. When the drop approaches the outlet, the interface

becomes perpendicular to the outlet boundary. Compared to previously tested convective

boundary condition, the constant-pressure outlet does affect the shape of the drop. However, this

deformation is reasonable because, based on the Laplace law, pressure difference across a curved

interface is ∆ = � ⁄ + ⁄ where and are the local radii of curvature. At the outlet boundary, the pressure is constant, which means that the pressure difference must be zero across

the interface. The radii of and , therefore, must be infinite. The behavior of the drop is thus

consistent with the constraint imposed by the boundary. Plots of the pressure fields in Figure

4.7(b) also demonstrate that a constant pressure difference is maintained at all times. The

pressure inside the drop is higher than that outside because of the interfacial tension. After the

drop moves out of the domain, the pressure field reverts to that of the single-phase case.

The second test is an immiscible displacement, conducted in the same computational domain

as that used in the moving drop test. Initially, the section is occupied by the blue

phase; the rest of the channel is occupied by the red phase. Fluid flow is driven by a constant

pressure difference. As in pressure-driven flows, varies with time if viscosities of two phases

are not equal, we use the and when the computational domain is fully saturated with the

defending fluid to compute R� and Ca . The parameters for the simulation are: �� = � = , �� = � = , � = . × − , � = . , � = . , R� = . , Ca = . ,

contact angle = . Figure 4.8(a) shows the evolution of the interface in the two-phase displacement test. Similar

to the moving drop case shown in Figure 4.7(a), as the red phase leaves the computational

domain, the interface becomes perpendicular to the outlet boundary. Pressure fields presented in

Figure 4.8(b) show that a constant pressure difference is also maintained. Due to the curved

interface, the pressure in the domain also varies along the cross-flow direction under the effect of

capillary pressure.

To further test the performance of the pressure boundary condition, we varied the viscosity

ratio (ratio of kinematic viscosity between the defending and invading fluids) and capillary

number Ca that are the two most important factors determining the flow patterns (Lenormand et al. 1988). Figure 4.8 also shows the simulation results with = and Ca = . , while keeping other parameters unchanged. In both cases, two-phase flows generated are physically correct,

which demonstrates that the pressure boundary can handle a range of viscosity ratios and

capillary numbers.

Figure 4.8 Phase distribution (a) and corresponding pressure field (b) of two-phase displacement in a channel using the pressure boundary condition. The dimensionless pressure is defined as +⁄ .

B. Capillary filling

Displacement of a wetting fluid with a non-wetting fluid in a capillary tube with a constant

pressure differential, as shown in Figure 4.9, is a problem with a moving contact line that admits

an analytical solution. For a capillary tube that has a diameter , when one fluid (red fluid) is

injected to displace another fluid (blue fluid) under a constant pressure differential, assuming that

the process is quasi-steady and there is no film drainage of the wetting fluid, it is straightforward

to derive that the immiscible displacement can be described by the analytical expression below:

� = − − + √ − − ∙ − + �− ∙

(4.22)

where � is the length occupied by the red phase, � and are, respectively, the viscosity of the red and blue phases, and are the inlet and outlet pressures, � is the interfacial tension,

is the dynamic contact angle (Pooley et al. 2009). Based on this equation, the breakthrough

time ( �⁄ = ) is = ∙ � +− + � (4.23)

To compare the numerical results from simulation using pressure boundary condition to the

analytical solution from Eqns. (4.22) and (4.23), a capillary tube that has length = and diameter = (lattice spacing) is set up. Initially, section . is occupied by the red phase, and the rest of the channel is occupied by the blue phase. Displacement simulations

were run with the following parameters: �� = � = , = × − . To study the effects of viscosity ratio and wettability, four sets of simulations were conducted: (a) = , � = . , � = . , = − . ; (b) = , � = . , � = . , = . ; (c) = . ,

� = . , � = . , = − . ; (d) = , � = . , � = . , = − . .

Figure 4.9 Simulation setup for capillary filling. The capillary tube has a length L and diameter d. The red fluid is the invading fluid and the blue fluid is the defending fluid.

Figure 4.10 compares the simulation results and the analytical solution from Eqn. (4.22). It

can be seen that, in the three drainage cases (case a, c, d), our simulation results are consistent

with the analytical solutions for all the different viscosity ratios / pressure differences. However,

in the imbibition case, the interface in the simulation moves slightly slower than the theoretical

result. The pressure distribution along the flow direction is shown in Figure 4.11. In all four

cases, a constant pressure difference between the inlet and the outlet is maintained and the

pressure gradient within each phase reflects different viscosities correctly.

The breakthrough times from the simulations are compared to the analytical solutions from

Eqn. (4.23) in Table 4.1. Again, in all the three drainage cases (case a, c, d), the simulation

results are in good agreements with analytical results, while there is about 8% error in the

imbibition case.

From the capillary filling tests, the accuracy of the pressure boundary condition is verified

and the simulation results admit analytical solutions.

Table 4.1 Comparison of breakthrough times in the case of capillary filling Type , ℎ , %

drainage 1 7.24E+5 7.2E+5 −0.55%

imbibition 1 7.20E+5 7.8E+5 8.3%

drainage 0.1 7.98E+5 7.7E+5 −3.5%

drainage 10 7.98E+5 7.9E+5 −1.0% % = ( , − , ℎ )/ , ℎ × %

CHAPTER 5

SIMULATION OF IMMISCIBLE DISPLACEMENTS IN A POROUS MEDIUM

The objective of this chapter is to apply the developed parallel multiphase LB simulator with

the color gradient model and the developed pressure boundary condition for multiphase flow to

simulate a few meaningful cases of immiscible displacements in a porous medium made up by

packed spheres.

The sphere packing is constructed by using a DEM (Discrete Element Modeling) method to

fill spherical particles into a cm × cm × cm cubic domain (Petunin et al. 2013). In the sphere packing, the particle diameter obeys a Gaussian distribution and varies from 1.21 mm to

5.97 mm; the number-averaged diameter is 3.34 mm. In simulation, the sphere packing is

discretized into a × × voxel domain using a resolution of 0.2 mm. A column is then cut out from the cubic sphere packing. The length of the column L is 320, and the diameter

of the column d is 100. The slice of = is set as the inlet and the slice of = is set as the outlet. The column (Figure 5.1) is then sealed on the side using solid no-slip walls to only allow

axial flow. To uniformly distribute the pressure and eliminate the effect of solid phase at the

boundaries, we added ten layers of fluid nodes to the inlet and the outlet. The porosity, �, of the column is 22.8%.

Immiscible two-phase flows on the Darcy level are generally described by the multiphase

extension of Darcy’s law

= �� , ∇ , � = , (5.1)

where � is the intrinsic permeability, � , is the relative permeability of fluid i evaluated at a particular wetting phase saturation , is the superficial velocity of fluid i through the porous

medium, ∇ is the pressure gradient of fluid i and subscripts w and nw denote wetting and non- wetting phases, respectively. In this study, the relative permeability is calculated by

� , = = ∇ =∇ , � = , (5.2) With direct numerical simulation, the pressure gradients of both phases can be obtained, which is

a distinctive advantage over experiments where phase pressures are difficult, if not impossible, to

obtain.

Figure 5.1 The cylindrical porous medium made up by packed spheres.

We first apply the pressure boundary to calculate the intrinsic permeability � by setting the sphere packs saturated with only one phase. We note that with MRT, the absolute permeability

obtained, . × − m , shows little viscosity dependence. To simulate immiscible displacement, initially, the entire column (> ) is saturated by one fluid (defending fluid) and the space before the porous column ( ) is filled with the

other fluid (invading fluid). The invading fluid is then pushed into the porous column by the

higher pressure to the inlet. At the outlet, a lower pressure is applied. The densities of both

phases were set at 1. The kinematic viscosity of the defending fluid was fixed at 0.167, and the

kinematic viscosity of the invading fluid was varied. By varying the pressure difference relative

to the viscosity and the interfacial tension of the fluids, flows with different capillary numbers

can be obtained. Results of four cases with different capillary numbers, viscosity ratios, and

wettability, shown in Table 5.1, are reported in this thesis. Note that in all cases, the density ratio

was kept at unity. As is evidenced in Eqns. (4.22) and (4.23), when the Reynolds number is low,

the density ratio does not play a role in the dynamics of displacement. We run all the cases until

a steady saturation is reached in the end. The simulation time for each case is around 70 hours by

using 256 cores. The parallel performance of the code can be found in Chapter 3.

Table 5.1 Parameters in different cases of immiscible displacement Case ∆P A cosθ Ca

A 0.0533 1E−4 −0.5 1 1.5E−1

B 0.0533 1E−4 −0.5 10 2.0E−1

C 0.0533 1E−1 −0.5 1 1.5E−4

D 0.0533 1E−1 0.5 1 1.5E−4

To clearly visualize the immiscible displacement inside the porous medium, Figure 5.2

through Figure 5.5 show the distributions of the two phases and corresponding pressure at the

slice cutting through the center of the column at different times ∗ for Case A, B, C, D respectively. Here the dimensionless time ∗ is defined as the ratio between the simulation time and the residence time of the defending fluid in single-phase flow:

∗ = ��∆ (5.3) In Figure 5.2(a) through Figure 5.5(a), the red color represents the invading fluid, and the yellow

color represents the defending fluid. The blue color represents the solid phase. Figure 5.2(b)

through Figure 5.5(b) show the wetting phase pressure distributions along the column, also at the

center slice at different times. It can be observed that the defending fluid is displaced by the

invading fluid under a constant pressure difference between inlet and outlet. The results

demonstrate that the pressure boundary condition can be applied to simulate multiphase flow in

complex porous media.

In an attempt to quantify the displacement efficiency along the axis of the column, we divide

the column equally into ten sections along the flow direction. In every section, the saturations

and pressures of both phases are calculated. Figure 5.6 shows the saturation distribution along

the column at different times for the four cases. Pressure distributions along the flow direction

are presented in Figure 5.7. Figure 5.8 shows the relative permeability curves using Eqn. (5.2)

for Case A through D.

By comparing the simulation results from different cases, we can illustrate 1) the effect of

capillary number, 2) the effect of viscosity ratio, and 3) the effect of wettability on two-phase

displacement.

Figure 5.2 Two-phase distribution (a) and wetting-phase pressure distribution (b) in the sphere packing at different times (∗= 0.12, 0.47, 0.71, 1.18), Case A.

Figure 5.3 Two-phase distribution (a) and wetting-phase pressure distribution (b) in the sphere packing at different times (∗= 0.12, 0.47, 0.71, 1.18), Case B.

Figure 5.4 Two-phase distribution (a) and wetting-phase pressure distribution (b) in the sphere packing at different times (∗= 0.12, 0.47, 0.71, 1.18), Case C.

Figure 5.5 Two-phase distribution (a) and wetting-phase pressure distribution (b) in the sphere packing at different times (∗= 0.12, 0.47, 0.71, 1.18), Case D.

Figure 5.6 Saturation distributions at different times in Case A, B, C, and D. Time interval ∆ ∗ between consecutive lines is 0.059. The red line represents the saturation distribution when breakthrough occurs.

Figure 5.8 Relative permeability vs. wetting phase saturation of the sphere packing. Red and green points denote breakthrough points. Time interval ∆ ∗ between consecutive points is 0.059.

1. Effect of capillary number

To evaluate the effect of the capillary number, Case A and Case C were conducted with Ca =. and Ca = . × − , respectively. Case A represents situations with ultra-low interfacial tensions that are generally achieved with surfactants or with nearly miscible fluids such as crude

oil and supercritical CO2 at high pressures. Case C represents situations where the two fluids are

immiscible and the interfacial tension is high, e.g., a water-oil displacement. The static contact

angle is °. Case A and Case C therefore represent typical drainage processes. Figure 5.6 shows that as time progresses, the saturation of invading non-wetting phase

increases and a steady saturation profile is reached in the end for both Case A and Case C. As the

porous medium in our simulation is homogeneous, the displacement front takes a Buckley-

Leverett form (Buckley and Leverett 1942) and displacement efficiency is high in both cases. An

obvious earlier breakthrough is observed in the high Ca case (Case A, breakthrough when ∗ =. , = . ) compared to the low Ca case (Case C, breakthrough when ∗ = . , =. ). The increased breakthrough time for the low Ca case is the result of significant adverse capillary pressure. From Figure 5.7(a) and Figure 5.7(c), the pressures of the two phase are

nearly identical when Ca = . , but are very different when Ca = . × − . Another notable observation is the difference between the pressure gradients of the two phases. When Ca is low, at the front of displacement, the non-wetting phase first displaces the wetting phase in larger

pores where the capillary pressure is insignificant. The pressures of the non-wetting phase and

the wetting phase are therefore nearly identical at the front of displacement. In the region behind

the front, the non-wetting phase can enter smaller pores than at the front due to flow pressure;

this leads to a higher difference between the pressures of the non-wetting and wetting phases.

After breakthrough, the pressure gradients of wetting and non-wetting phase become nearly

identical with a nearly constant capillary pressure difference. The relative permeability curves of

the two phases are shown in Figure 5.8(a). We noticed that the relative permeability of both

phases decrease with decreasing capillary number. When the capillary number is small, flow is

blocked in some small pores, which decreases the relative permeability of both invading and

defending fluids. Similar reductions have also been observed in other numerical and

experimental studies (Li et al. 2005; Blunt and King 1991; Avraam and Payatakes 1999).

2. Effect of viscosity ratio

By comparing Case A and Case B, we can observe the effect of viscosity ratio on the

dynamics of multiphase flow. To eliminate the effect of the capillary number, both Case A and

Case B were assigned with a high capillary number so that the flow is dominated by viscous

forces. When ≠ , the interfacial tension can be obtained using Laplace’s law by a bubble test. In Case B, when = and = − , the interfacial tension expressed in lattice unit is . ×

− . According to the saturation profiles in Figure 5.6, in Case B, when the viscosity of the

invading fluid is much lower than the viscosity of the defending fluid, there is a much quicker

breakthrough (breakthrough when ∗ = . , = . ) than Case A (breakthrough when ∗ =. , = . ). At the point of breakthrough, significant amount of wetting fluid is left behind. This phenomenon is known as viscous fingering (Saffman and Taylor 1958; Homsy 1987). The

pressure profile in Figure 5.7(b) concaves up, which is another indication that the less viscous

fluid has developed viscous fingers into the more viscous fluid. In Figure 5.8(b), the shapes of

the relative permeability curves at high M also shows signs of unstable and fingered

displacement. Compared with Case A, the less viscous non-wetting fluid’s relative permeability

concaves down because it breaks into a more viscous fluid, and therefore, has a slower speed

than expected. The more viscous wetting fluid’s relative permeability concaves up, for the same

reason, as it is being penetrated by the front of the less viscous fluid. This trend is in agreement

with results from pore network modeling by Blunt and King (1991), though the flows there were

driven by a rate boundary applied to the inlet and a pressure boundary applied to the outlet. In

some previous LB simulations (Li et al. 2005; Ghassemi and Pak 2011; Huang and Lu 2009),

relative permeability is calculated using steady-state flows through periodic domains, driven by a

body force. In those studies, is defined as the ratio of viscosity between the non-wetting and

wetting fluids. An increased leads to significantly increased non-wetting phase relative

permeability. The relative permeability of the wetting-phase, however, is insensitive to the

change in (Li et al. 2005), which is very different from our observations. This difference is a

clear evidence that relative permeability depends on the way in which multiphase flows are

generated and maintained, as well as the steadiness of the flows.

3. Comparison between drainage and imbibition

Case C and Case D compare a drainage process and an imbibition process in the sphere

packing. Saturation profiles are presented in Figure 5.6(c) and Figure 5.6(d) respectively. A

much quicker breakthrough and higher residual saturations are observed for the imbibition case

(breakthrough when ∗ = . , = . ) than the drainage case (breakthrough when ∗ =. , = . ). The pressure distribution during imbibition is shown in Figure 5.7(d). Compared to drainage, at the front of imbibition, the wetting phase first occupies the small pores.

This preferential occupation leads to a significant capillary pressure at the front of displacement.

In the region behind the front, the wetting phase can gradually enter larger pores; this leads to a

lower difference between the pressures of the non-wetting and wetting phases. After

breakthrough, the pressure gradients of wetting and non-wetting phase become similar. During

imbibition, the wetting phase tends to progress along the small pores as well as the corners of the

pores, leaving the non-wetting phase trapped in large pores. In later stage of imbibition, the non-

wetting phase’s distribution becomes more disconnected in the imbibition process than in the

drainage process. The disconnected non-wetting phase saturation reduces the relative

permeability. Note that saturations that reflect later stages of imbibition correspond to earlier

stages of drainage. As shown in Figure 5.8(c), the relative permeability of the non-wetting phase

in imbibition is initially higher than that in the drainage at high non-wetting phase saturations;

when the non-wetting saturation is low, the relative permeability of the non-wetting phase in

imbibition becomes similar to that in drainage. These relative permeability curves are therefore

reasonable. Similar findings were reported in an experimental study (Naar et al. 1962) and pore

network modeling work (Blunt 1997) for unconsolidated glass beads. For the wetting phase, it is

always stated that there is little hysteresis; compared to drainage, both higher and lower relative

permeability of the wetting phase during imbibition have been reported (Jerauld and Salter 1990).

In this work, we observed that the wetting phase relative permeability in imbibition is lower than

that in drainage.

CHAPTER 6

CONCLUSIONS AND FUTURE WORK

This thesis covers the development of a multiphase flow simulator using the lattice

Boltzmann method, evaluation of several boundary conditions, formulation of pressure boundary

condition, and simulations of pressure driven multiphase flow through a complex porous

medium made up by packed spheres. In this chapter, the main contributions of this work are

summarized, followed by several recommendations for future work.

6.1 Summary and Conclusions

In this work, first we develop a multiphase lattice Boltzmann simulator based on the color-

gradient model. This color gradient model can simulate fluids with large density ratios, viscosity

ratios and produce controlled interfacial tension and contact angle. Simulation of drop

deformation in a simple shear flow shows that this multiphase model is able to accurately

simulate dynamic problems.

To simulate multiphase flows in a meaningful size that often contains up to 109 lattice nodes,

the simulator is parallelized for distributed-memory clusters using MPI (message passing

interface). The parallel performance of the code is tested on BlueM, a high-performance

computing system at Colorado School of Mines. For the test conducted in porous media that

contain 6403, 3203 and 1603 lattice nodes, a maximum speedup of 208 is achieved using 256

cores, yielding a 0.81 parallel efficiency. These numbers are reasonable upon considering that

the domain size and therefore the load of computation assigned to each core is small relative to

the communication cost. As the computation loaded to each core increases, higher speedup and

parallel efficiency are expected.

We note that boundary conditions are important for simulation of multiphase flow through

porous media, and only a few work have been done to implement appropriate boundary

conditions. In this study, we developed / tested several boundary conditions for the color gradient

multiphase LB model, including the convective outflow boundary condition, the modified

periodic boundary condition, and the pressure boundary condition. The convective outflow

boundary condition allows the fluids to exit the computation domain with minimum unphysical

disturbance. The modified periodic boundary also allows fluids to exit freely and can be used to

simulate gravity-driven multiphase flows such as flows in a centrifuge. The pressure boundary

condition, validated in a capillary filling test, enables us to simulate pressure-driven multiphase

flow using the lattice Boltzmann method. The pressure boundary condition, specifically, has not

been developed before this study and therefore is a new contribution.

As an application, the developed simulator is used to study pressure-driven immiscible

displacement in a porous medium made up by packed spheres. By using direct numerical

simulation, pressures of each phase can be obtained, which is a distinctive advantage over

experiments where phase pressures are difficult to measure. Compared to previous LB

simulations, we observe similar effect from the capillary number on multiphase flow. However,

significant difference in the effect of viscosity ratio is observed, owning to the differences in the

way multiphase flow is generated. Using the pressure boundary condition, drainage and

imbibition can be separately evaluated by simulating their respective displacement processes;

these two processes cannot be distinguished in simulations with periodic boundaries with body

force. The relative permeability curves obtained are consistent with published pore network

modeling and experimental results.

In conclusion, we developed a capable and flexible lattice Boltzmann simulator for

multiphase flow. Our simulator is validated and the results are consistent with known

observations. Our parallel simulator allows us to simulate multiphase flow in a cm × cm × cm cubic domain using a 0.2 mm resolution which requires 70 hours’ simulation time running on 256 cores. The increasing parallel efficiency with increasing computational load of every core

indicates that the simulator should be capable of carrying out higher scale simulations on capable

clusters. This lattice Boltzmann simulator is a very useful tool to solve multiphase flow through

porous media; however, it is highly desirable to compare lattice Boltzmann results with other

methods.

6.2 Future work

Development of this multiphase flow simulator together with flexible boundary conditions

opens up many research opportunities. Of the many applications and theoretical studies that can

be conducted, we suggest a few directions that can be pursued in the near future:

 With the different boundary conditions developed in this thesis, we can study the effect of multiphase flow driving mechanisms on relative permeability including the pressure-driven,

rate-driven and gravity-driven flows. This study is a direct application of the methods

developed in this thesis, and the results should be beneficial for improving the understanding

of relative permeability.

 Geomechanical stress plays an important role in controlling single- and multiphase flow in porous media. Stress sensitivity of relative permeability can be modeled by simulating

multiphase flow in a porous medium that has stress-induced deformations. No similar studies

have been reported and therefore it is likely to generate a new research direction.

 Microfluidics is a very useful means to visualize multiphase flow and to study the effects of pore geometry, network topology, and fluid properties. As experimental work on

microfluidics is being carried out in our group, simulations of microfluidic experiments can

be conducted. By comparison with experiments, we will have opportunities for validation of

the simulation as well as the experimental results.

REFERENCES CITED

Aidun, C.K. and Clausen, J.R., 2010. “Lattice-Boltzmann method for complex flows”. Annual

Review of Fluid Mechanics, 42, 439-472.

Ahrenholz, B., Tolke, T., Lehmann, P., Peters, A., Kaestner, A., Krafczyk, M. and Durner, W.,

2008. “Prediction of capillary hysteresis in a porous material using lattice Boltzmann methods

and comparison to experimental data and a morphological pore network model.” Advance in

Water Resources, 31, 1151-1173.

Avraam, D.G., Payatakes, A.C., 1999. “Flow mechanisms, relative permeabilities, and coupling

effects in steady-state two-phase flow through porous media. The case of strong wettability.”

Industrial & Engineering Chemistry Research, 38, 778-786.

Badalassi, V.E., Ceniceros, H.D. and Banerjee, S., 2003. “Computation of multiphase systems

with phase field models.” Journal of Computational Physics, 190, 371-397.

Bhatnagar, P.L., Gross, E.P. and Krook, M. 1954. “A model for collision processes in gases. I.

Small amplitude processes in charged and neutral one-component systems.” Physical Review

94(3): 511.

Blunt, M, 1997. “Effects of heterogeneity and wetting on relative permeability using pore level

modeling.” SPE Journal, 2, 70-87, SPE-36762-PA.

Blunt, M., King, P., 1991. “Relative permeabilities from two- and three-dimensional pore-scale

network modeling.” Transport in Porous Media, 6, 407-433.

Buckley, S.E. and Leverett, M.C., 1942. “Mechanism of Fluid Displacement in Sands.” SPE

Journal, 146, 107-116, SPE-942107-G.

Cahn, J.W. and Hilliard, J.E., 1958. “Free energy of a nonuniform system. I. Interfacial free

energy.” The Journal of Chemical Physics, 28, 258.

Cahn, J.W. and Hilliard, J.E., 1959. “Free energy of a nonuniform system III. Nucleation in a

two-component incompressible fluid.” The Journal of Chemical Physics, 31, 688.

Chen, L., Kang, Q., Mu, Y., He, Y. and Tao, W., 2014. “A critical review of the pseudopotential

multiphase lattice Boltzmann model: Methods and applications.” International Journal of

Heat and Mass Transfer, 76, 210-236.

Chen, S., Chen, H., Martnez, D. and Matthaeus, W., 1991. “Lattice Boltzmann model for

simulation of magnetohydrodynamics.” Physical Review Letters, 67, 3776-3779.

Chen, S., Martínez, D. and Mei., R., 1996. “On boundary conditions in lattice Boltzmann

methods." Physics of Fluids, 8, 2527.

Chen, S. and Doolen, G.D., 1998. “Lattice Boltzmann method for fluid flows.” Annual Review

of Fluid Mechanics, 30, 329-364.

Cueto-Felgueroso, L. and Juanes, R., 2013. “Forecasting long-term gas production from shale.”

Proceedings of the National Academy of Science of the United States of America, 110,

19660-19661.

Demianov, A., Dinariev, O., Evseev, N., 2011. “Density functional modelling in multiphase

compositional hydrodynamics.” The Canadian Journal of Chemical Engineering, 89, 206-

226.

d'Humières, D. 1992. “Generalized lattice-Boltzmann equations.” Rarefied Gas Dynamics:

Theory and Simulations, 450-458.

d'Humières, D., 2002. “Multiple-relaxation-time lattice Boltzmann models in three dimensions.”

Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical

and Engineering Sciences, 360, 437-451.

d'Ortona, U., Salin, D., Cieplak, M., Rybka, R. B. and Banavar, J. R., 1995. “Two-color

nonlinear Boltzmann cellular automata: Surface tension and wetting.” Physical Review E,

51, 3718-3728.

Dong, B., Yan, Y.Y., Li, W. and Song, Y., 2010. “Lattice Boltzmann simulation of viscous

fingering phenomenon of immiscible fluids displacement in a channel.” Computers & Fluids,

39, 768-779.

Dudukovic, M.P., Larachi, F. and Mills, P.L., 1999. “Multiphase reactors-revisited.” Chemical

Engineering Science, 54, 975-995.

Esposito, S.J. and Thomson, N.R., 1999. “Two-phase flow and transport in a single fracture-

porous medium system.” Journal of Contaminant Hydrology, 37, 319-341.

Frisch, U., d'Humières, D., Hasslacher, B., Lallemand, P., Pomeau, Y. and Rivet, J.P., 1987.

“Lattice gas hydrodynamics in two and three dimensions.” Complex System 1, 649-707.

Gerlach, D., Tomar, G., Biswas, G. and Durst, F., 2006. “Comparison of volume-of-fluid

methods for surface tension-dominant two-phase flows.” International Journal of Heat and

Mass Transfer, 49, 740-754.

Ghassemi, A., and Pak, A., 2011. “Numerical study of factors influencing relative permeabilities

of two immiscible fluids through porous media using lattice Boltzmann method.” Journal of

Petroleum Science and Engineering, 77, 135-145.

Grunau, D., Chen, S. and Egger, K., 1993. “A lattice Boltzmann model for multi-phase fluid

flows.” Physics of Fluids A, 5, 2557.

Gunstensen, A.K., Rothman, D.H., Zaleski, S. and Zanetti, G., 1991. “Lattice Boltzmann model

of immiscible fluids.” Physical Review A, 43, 4320.

Gunstensen, A.K. and Rothman, D.H., 1993. “Lattice-Boltzmann studies of immiscible two-

phase flow through porous media.” Journal of Geophysical Research: Solid Earth, 98,

6431-6441.

Guo, Z., Zheng, C. and Shi, B. 2002. “Discrete lattice effects on the forcing term in the lattice

Boltzmann method.” Physical. Review E, 65, 046308.

Hassanizadeh, M. and Gray, W.G., 1979. “General conservation equations for multi-phase

systems: 1. Averaging procedure.” Advance in Water Resources, 2, 131-144.

He, X., Chen, S. and Zhang, R., 1999. “A lattice Boltzmann scheme for incompressible

multiphase flow and its application in simulation of Rayleigh-Taylor instability.” Journal of

Computation Physics, 152, 642-663.

Hirt, C.W. and Nichols, B.D., 1981. “Volume of fluid (VOF) method for the dynamic of free

boundaries.” Journal of Computation Physics, 39, 201-225.

Homsy, G.M., 1987. “Viscous fingering in porous media.” Annual Review of Fluid Mechanics,

19, 271-311.

Huang, H., Huang, J. and Lu, X., 2014. “Study of immiscible displacements in porous media

using a color-gradient-based multiphase lattice Boltzmann method.” Computers & Fluids, 93,

164-172.

Huang, H. and Lu, X., 2009. “Relative permeabilities and coupling effects in steady-state gas-

liquid flow in porous media: A lattice Boltzmann study.” Physics of Fluids, 21, 092104.

Huang, H., Wang, L. and Lu, X., 2011. “Evaluation of three lattice Boltzmann models for

multiphase flows in porous media.” Computers and Mathematics with Applications, 61, 3606-

3617.

Hunt, J.R., Sitar, N. and Udell, K.S., 1988. “Nonaqueous phase liquid transport and cleanup 1.

analysis of mechanism.” Water Resource Research, 24, 1247-1258.

Inamuro, T., Konishi, N. and Ogino, F., 2000. “A Galilean invariant model of the lattice

Boltzmann method for multiphase fluid flows using free energy approach.” Computer Physics

Communications. 129, 32-45.

Jerauld, G.R. and Salter, S.J., 1990. “The effect of pore-structure on hysteresis in relative

permeability and capillary pressure: pore-level modeling.” Transport in Porous Media, 5,

103-151.

Joekar-Niasar, V. and Hassanizadeh, S.M., 2012. “Analysis of fundamentals of two-phase flow

in porous media using dynamic pore-network models: A review.” Critical Reviews in

Environmental Science and Technology, 42, 1895-1976.

Kalarkis, A.N., Burganos, V.N. and Payatakes, A.C., 2002. “Galilean-invariant lattice Boltzmann

simulation of liquid-vapor interface dynamics.” Physical Review E, 65, 056702.

Kandhai, D., Koponen, A., Hoekstra, A., Kataja, M., Timonen, J. and Sloot, P.M.A., 1999.

“Implementation aspects of 3D lattice BGK: Boundaries, accuracy, and a new fast relaxation

method.” Journal of Computation Physics, 150, 482-501.

Kang, Q., Zhang, D. and Chen, S., 2004. “Immiscible displacement in a channel: simulation of

fingering in two dimensions.” Advance in Water Resources, 27, 13-22.

Kang, Q., Zhang, D. and Chen, S., 2005. “Displacement of a three-dimensional immiscible

droplet in a duct.” Journal of Fluid Mechanics, 545, 41-66.

Karimi, S., Saidian, M., Prasad, M. and Kazemi, H. 2015. “Reservoir rock characterization using

centrifuge and nuclear magnetic resonance: A laboratory study of middle Bakken formation.”

SPE-175069-MS, presented at the SPE Annual Technical Conference and Exhibition held in

Houston, Texas, USA.

Kim, S. H. and Pitsch, H., 2007. “A generalized periodic boundary condition for lattice

Boltzmann method simulation of a pressure driven flow in a periodic geometry.” Physics of

Fluids, 19, 108101.

Kreutzer, M.T., Kapteijn, F., Moulijn, J.A. and Heiszwolf, J.J., 2005. “Multiphase monolith

reactors: chemical reaction engineering of segmented flow in microchannels.” Chemical

Engineering Science, 60, 5895-5916.

Kupershtokh, A.L., Medvedev, D.A. and Karpov, D.I., 2009. “On equations of state in a lattice

Boltzmann method.” Computers & Mathematics with. Applications, 58, 965-974.

Ladd, A.J.C., 1994. “Numerical simulations of particulate suspensions via a discretized

Boltzmann equation. Part 1. Theoretical foundation.” Journal of Fluid Mechanics 271, 285.

Latva-Kokko, M. and Rothman, D.H., 2005. “Diffusion properties of gradient-based lattice

Boltzmann models of immiscible fluids.” Physical Review E, 71, 056702.

Leclaire, S., Reggio, M. and Trepaníer, J., 2012. “Numerical evaluation of the two recoloring

operators for an immiscible two-phase flow lattice Boltzmann model.” Applied Mathematical

modeling, 39, 768.

Lenormand, R., Touboul, E. and Zarcone, C., 1988. “Numerical models and experiments on

immiscible displacements in porous media.” Journal of Fluid Mechanics, 189, 165-187.

Li, H., Pan, C. and Miller, C.T., 2005. “Pore scale investigation of viscous coupling effects for

two-phase flow in porous media.” Physical Review E, 72, 026705.

Liu, H., Valocchi, A.J., and Kang, Q., 2012. “Three-dimensional lattice Boltzmann model for

immiscible two-phase flow simulations.” Physical Review E, 85, 046309.

Liu, H., Valocchi, A.J., Zhang, Y. and Kang, Q., 2014. “Lattice Boltzmann phase-field modeling

of thermocapillary flows in a confined microchannel.” Journal of Computation Physics, 256,

334-356.

Liu, H., Valocchi, A.J., Werth, C., Kang, Q. and Oostrom M., 2014. “Pore-scale simulation of

liquid CO2 displacement of water using a two-phase lattice Boltzmann model.” Advance in

Water Resources, 73, 144-158.

Liu, H., Zhang, Y. and Valocchi, A.J., 2012. “Modeling and simulation of thermocapillary flows

using lattice Boltzmann method.” Journal of Computation Physics, 231, 4433-4453.

Lowengrub, J. and Truskinovsky, L., 1998. “Quasi-incompressible Cahn-Hilliard fluids and

topological transitions.” Proceedings of the Royal Society A: Mathematical, Physical and

Engineering Sciences, 454, 2617-2654.

Lou, Q., Guo, Z. and Shi, B., 2013. “Evaluation of outflow boundary conditions for two-phase

lattice Boltzmann equation.” Physical Review E, 87, 063301.

Mackay, D.M. and Cherry, J.A., 1989. “Groundwater contamination: pump-and-treat

remediation.” Environmental Science & Technology, 23, 630-636.

Maier, R. S., Bernard, R.S. and Grunau, D.W., 1996. “Boundary conditions for the lattice

Boltzmann method.” Physics of Fluids, 8, 1788.

Naar, J., Wygal, R.J. and Henderson, J.H., 1962. “Imbibition Relative permeability in

unconsolidated porous media.” SPE Journal, 2, 13-17, SPE-213-PA.

Orr Jr, F.M. and Taber, J.J., 1984. “Use of carbon dioxide in enhanced oil recovery.” Science,

224, 563-569.

Osher, S. and Fedkiw, R.P., 2001. “Level set methods: an overview and some recent results.”

Journal of Computation Physics, 169, 463-502.

Osher, S. and Sethian, J.A., 1988. “Fronts propagating with curvature-dependent speed:

Algorithms based on Hamilton-Jacobi formulations.” Journal of Computation Physics, 79,

12-49.

Pan, C., Hilpert, M. and Miller, C.T., 2004. “Lattice-Boltzmann simulation of two-phase flow in

porous media.” Water Resource Research, 40, W01501.

Pan, C., Luo, L.-S. and C. T. Miller, 2006. “An evaluation of lattice Boltzmann schemes for

porous medium flow simulation.” Computer & Fluids 35, 898-909.

Petunin, V., Labra, C., Xiao, F., Tutuncu, A., Sun, J. and Yin, X., 2013. “Porosity and

permeability change under stress and correlation to rock texture.” Poromechanics, V, 1037-

1046.

Pilliod Jr, J.E. and Puckett, E.G., 2004. “Second-order accurate volume-of-fluid algorithms for

tracking material interfaces.” Journal of Computation Physics, 199, 465-502.

Pooley, C.M., Kusumaatmaja, H. and Yeomans, J.M., 2009. “Modelling capillary filling

dynamics using lattice Boltzmann simulations.” European Physical Journal-Special Topics,

171, 63-71.

Qian, Y. H., d'Humières, D. and Lallemand, P. 1992. “Lattice BGK models for Navier-Stokes

equation.” Europhysics Letters, 17, 479-484.

Reis, T. and Phillips, T.N., 2007. “Lattice Boltzmann model for simulating immiscible two-

phase flows.” Journal of Physics A, 40, 4033.

Rothman, D.H. and Keller, J.M., 1988. “Immiscible cellular-automaton fluids.” Journal of

Statistical Physics ,52, 1119-1127.

Rider, W.J. and Kothe, D.B., 1998. “Reconstructing volume tracking.” Journal of Computation

Physics, 141, 112-152.

Saffman, P.G. and Taylor, G., 1958. “The penetration of a fluid into a porous medium or Hele-

Shaw cell containing a more viscous liquid.” Proceedings of the Royal Society of London.

Series A, Mathematical and Physical Sciences, 245, 312-329.

Sbragaglia, M., Benzi, R., Biferale, L., Succi, S., Sugiyama, K. and Toschi, F., 2007.

“Generalized lattice Boltzmann method with multirange pseudopotential.” Physical Review E,

75, 026702.

Scardovelli, R. and Zaleski, S., 1999. “Direct numerical simulation of free-surface and interfacial

flow.” Annual Review of Fluid Mechanics, 31, 567-603.

Shan, X., 2006. “Analysis and reduction of the spurious current in a class of multiphase lattice

Boltzmann models.” Physical Review E, 73, 047701.

Shan, X. and Chen, H., 2007. “A generalized multiple-relaxation-time Boltzmann collision

model.” International Journal of Modern Physics C, 18, 635-643.

Shan, X. and Chen, H., 1993. “Lattice Boltzmann model for simulating flows with multiple

phases and components.” Physical Review E, 47, 1815.

Shan, X. and Doolen, G., 1995. “Multicomponent lattice-Boltzmann model with interparticle

interaction” Journal of Statistical Physics, 81, 379-393.

http://www.worldscientific.com/worldscinet/ijmpc

Sussman, M., Fatemi, E., Smereka, P. and Osher, S., 1998. “An improved level set method for

incompressible two-phase flows.” Computers & Fluids, 27, 663-680.

Sussman, M., Smereka, P. and Osher, S., 1994. “A level set approach for computing solutions to

incompressible two-phase flow.” Journal of Computational Physics, 114, 146-159.

Swift, M.R., Osborn, W.R. and Yeomans, J.M., 1995. “Lattice Boltzmann simulation of nonideal

fluids.” Physical Review Letters, 75, 830.

Swift, M.R., Orlandini, E., Osborn, W.R. and Yeomans, J.M. 1996. “Lattice Boltzmann

simulation of liquid-gas and binary fluid system” Physical Review E, 54, 5041.

Szulczewski, M.L., MacMinn, C.W., Herzog, H.J., and Juanes, R., 2012. “Lifetime of carbon

capture and storage as a climate-change mitigation technology.” Proceedings of the National

Academy of Science of the United States of America, 109, 5185-5189.

Taylor, G.I., 1932. “The viscosity of a fluid containing small drops of another fluid.”

Proceedings of the Royal Society of London. Series A, 138, 41-48.

TÖlke, J., Krafczyk, M., Schulz, M. and Rank, E., 2002. “Lattice Boltzmann simulations of

binary fluid flow through porous media.” Philosophical Transactions of the Royal Society A,

360, 535.

Whitaker, S., 1998. “The method of volume averaging.” Springer Science & Business Media.

Wu, Q., Ok, J.T., Sun, Y., Retterer, S.T., Neeves, K.B., Yin, X., Bai, B. and Ma, Y., 2013.

“Optic imaging of single and two-phase pressure-driven flows in nano-scale channels.” Lab

on a Chip, 13, 1165-1171.

http://dx.doi.org/10.1039/1473-0189/2001

Wu, M., Xiao, F., Johnson-Paben, R.M., Retterer, S.T., Yin, X. and Neeves, K.B., 2012. “Single-

and two-phase flow in microfluidic porous media analogs based on Voronoi tessellation, Lab

on a Chip, 12, 253-261.

Xiao, F. and Yin, X., 2015. “Geometry models of porous media based on Voronoi tessellations

and their porosity–permeability relations.” Computers & Mathematics with Applications, In

press, available online 19 Oct 2015.

Xu, W., Ok, J.T., Xiao, F., Neeves, K.B. and Yin, X., 2014. “Effect of pore geometry and

interfacial tension on water-oil displacement efficiency in oil-wet microfluidic porous media

analogs.” Physics of Fluids, 26, 093102.

Yang, J. and Boek, E.S., 2013. “A comparison study of multi-component Lattice Boltzmann

models for flow in porous media applications.” Computers and Mathematics with

Applications, 65, 882-890.

Yuan, P. and Schaefer, L., 2006. “Equations of state in a lattice Boltzmann model.” Physics of

Fluids, 18, 042101.

Zhang, J. and Kwok, D.Y., 2006. “Pressure boundary condition of the lattice Boltzmann method

for fully developed periodic flows.” Physical Review E, 73, 047702.

Zhang, L., Kang, Q., Yao, J., Gao, Y., Sun, Z., Liu, H. and Valocchi, A.J., 2015. “Pore scale

simulation of liquid and gas two phase flow based on digital core technology.” Science China

Technological Science, 58, 1375-1384.

Zhang, J. and Tian, F., 2008. “A bottom-up approach to non-ideal fluids in the lattice Boltzmann

Method.” Europhysics Letters, 81, 66005.

Zhao, Y. and Fan, L., 2010. “Multirelaxation-time interaction-potential-based lattice Boltzmann

model for two-phase flow.” Physical Review E, 82, 046708.

http://www.sciencedirect.com/science/journal/08981221

https://www.researchgate.net/researcher/2001399492_Wei_Xu

https://www.researchgate.net/researcher/2006138318_Jeong_Tae_Ok

https://www.researchgate.net/researcher/59121949_Feng_Xiao

Zou, Q. and He, X., 1997. “On pressure and velocity boundary conditions for the lattice

Boltzmann BGK model.” Physics of Fluids, 9, 1591.