A Lattice Boltzmann Model for Diffusion of Binary
Gas Mixtures
University of Cambridge
Department of Engineering
This dissertation is submitted for the degree of
Doctor of Philosophy
by:
Sam Bennett
King’s College
2010
A Lattice Boltzmann Model for Diffusion of Binary Gas Mixtures
Sam Bennett
This thesis describes the development of a Lattice Boltzmann (LB) model for a binary gas
mixture. Specifically, channel flow driven by a density gradient with diffusion slip occurring
at the wall is studied in depth.
The first part of this thesis sets the foundation for the multi-component model used in
the subsequent chapters. Commonly used single component LB methods use a non-physical
equation of state, in which the relationship between pressure and density varies according
to the scaling used. This is fundamentally unsuitable for extension to multi-component
systems containing gases of differing molecular masses that are modelled with the ideal gas
equation of state. Also, existing methods for implementing boundary conditions are unsuit-
able for extending to novel boundary conditions, such as diffusion slip. Therefore, a new
single component LB derivation and a new method for implementing boundary conditions
are developed, and validated against Poiseuille flow. However, including a physical equa-
tion of state reduces stability and time accuracy, leading to longer computational times,
compared with ‘incompressible’ LB methods. The new method of analysing LB boundary
conditions successfully explains observations from other commonly used schemes, such as
the slip velocity associated with ‘bounce-back’.
The new model developed for multi-component gases avoids the pitfalls of some other
LB models, a single computational grid is shared by all the species and the diffusivity is
independent of the viscosity. The Navier-Stokes equation for the mixture and the Stefan-
Maxwell diffusion equation are both recovered by the model. However, the species mo-
mentum equations are not recovered correctly and this can lead to instability. Diffusion
slip, the non-zero velocity of a gas mixture at a wall parallel to a concentration gradient,
is successfully modelled and validated against a simple one-dimensional model for channel
flow. To increase the accuracy of the scheme a second order numerical implementation
is needed. This can be achieved using a variable transformation method which does not
result in an increase in computational time.
i
Simulations were carried out on hydrogen and water diffusion through a narrow channel,
with varying total pressure and concentration gradients. For a given value of the species
mass flux ratio, the total pressure gradient was dependent on the species concentration
gradients. These results may be applicable to fuel cells where the species mass flux ratio
is determined by a chemical reaction and the species have opposing velocities. In this case
the total pressure gradient is low and the cross-channel average mass flux of hydrogen is
independent of the channel width.
Finally, solutions for a binary Stefan tube problem were investigated, in which the
boundary at one end of a channel is permeable to hydrogen but not water. The water has
no total mass flux along the channel but circulates due to the slip velocity at the wall. The
cross-channel average mass flux of the hydrogen along the channel increases with larger
channel widths. A fuel cell using a mixture of gases, one being inert, will experience similar
circulation phenomena and, importantly, the width of the pores will affect performance.
This thesis essentially proves the viability of LB models to simulate multi-component
gases with diffusion slip boundaries, and identifies the many areas in which improvements
could be made.
ii
Declaration
This dissertation is my own work and contains nothing which is the outcome of work done
in collaboration with others, except as specified in the text and Acknowledgements. This
dissertation contains approximately 31,000 words and 43 figures. This dissertation has not
been submitted for a degree or diploma or other qualification at any other University.
Sam Bennett
Hopkinson Laboratory, Cambridge
21st of July, 2010
iii
Acknowledgements
I’m not sure how, or if, I would have completed the work presented in this thesis without two
people, Pietro Asinari of Politecnico di Torino and Paul Dellar of the University of Oxford.
They both recognised a student out of his depth in the quagmire of Lattice Boltzmann
literature and tossed him a lifeline, for which I am profoundly grateful. They invested their
time and energy in me and my work, their reward is only this small paragraph. I would
also like to thank their respective deputies, Salvador Izquierdo and Tim Reis, for ensuring
that every spare moment of my time in Turin and Oxford (and China) was spent eating
delicious food and having fun.
Whenever I have needed reminding that there is more to engineering than Chapman-
Enskog, Joy Ward and everyone involved in outreach (too numerous to mention) were
there; thank you all for cheering me up.
Without Peter Benie coaxing my computer, ‘Ortley’, through its old age I would never
have finished this thesis, thank you for that and a million other small things.
I would like to thank my supervisor and advisor, John Young and Graham Pullan, and
also Nick Collings.
Antonios Triantafyllidis, Andy Garmory and Robert Gordon all tried their best to an-
swer my impossible questions, thank you, and stick with the CFD. Cathy Elks kindly
pointed out the parts of this thesis which were unnecessarily obfuscated.
Finally I thank Rolls Royce Fuel Cells Ltd. and the EPSRC for their generous financial
support in exchange for the very occasional one page report.
iv
Contents
Contents v
List of Symbols ix
1 Introduction 1
1.1 Notes on style . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 An Introduction to LB Methods 5
2.1 A one-dimensional scheme (D1Q2) . . . . . . . . . . . . . . . . . . . . . . 5
2.2 D1Q2 example code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.3 Analysing the D1Q2 scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.4 One dimension, three velocity (D1Q3) scheme . . . . . . . . . . . . . . . . 12
2.5 Concluding remarks on the 1D examples . . . . . . . . . . . . . . . . . . . 14
3 Literature Review 17
3.1 Lattice Boltzmann Methods (LBM) . . . . . . . . . . . . . . . . . . . . . . 17
3.1.1 Lattice Boltzmann from the Lattice Gas Automaton . . . . . . . . 17
3.1.2 The scope of Lattice Boltzmann . . . . . . . . . . . . . . . . . . . . 19
3.1.3 The rise of Lattice BGK . . . . . . . . . . . . . . . . . . . . . . . . 20
3.1.4 Boundaries for LBM . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.1.5 Improvements on LBGK . . . . . . . . . . . . . . . . . . . . . . . . 23
3.1.6 LBM for multi-component fluids . . . . . . . . . . . . . . . . . . . . 24
3.1.7 LBM in porous media and solid oxide fuel cells . . . . . . . . . . . 26
v
3.2 Diffusion slip . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.2.1 Origin of diffusion slip . . . . . . . . . . . . . . . . . . . . . . . . . 29
4 Single-Component Lattice Boltzmann 32
4.1 The D2Q9 grid and the Boltzmann equation . . . . . . . . . . . . . . . . . 33
4.1.1 The D2Q9 grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.1.2 Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.1.3 The truncated nine moment system . . . . . . . . . . . . . . . . . . 35
4.1.4 The scaled evolution equation . . . . . . . . . . . . . . . . . . . . . 37
4.1.5 Moments of the evolution equation . . . . . . . . . . . . . . . . . . 38
4.1.6 The equation of state . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.2 The equilibrium distribution . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.2.1 Mass conservation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.2.2 Momentum conservation . . . . . . . . . . . . . . . . . . . . . . . . 40
4.2.3 Higher moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.2.4 Recovering the equilibrium distribution function . . . . . . . . . . . 42
4.3 The MRT collision operator . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.3.1 Evolution of moments . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.3.2 The momentum flux density tensor (ΠXX ,ΠXY and ΠY Y ) . . . . . . 46
4.3.3 The momentum equation . . . . . . . . . . . . . . . . . . . . . . . . 48
4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5 Boundary Conditions and Numerical Implementation 52
5.1 Boundary conditions for LB . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.1.1 Kinetic style boundaries . . . . . . . . . . . . . . . . . . . . . . . . 52
5.1.2 Moment grouping . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.1.3 Analysis of Zou and He boundary conditions . . . . . . . . . . . . . 55
5.1.4 Analysis of ‘bounce-back’, specular and diffusive boundary conditions 58
5.1.5 Treatment of corner nodes . . . . . . . . . . . . . . . . . . . . . . . 60
vi
5.2 Numerical implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
6 Solution of Poiseuille Flow 66
6.1 Theoretical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
6.1.1 Moment based boundary conditions for Poiseuille Flow . . . . . . . 68
6.2 Non-dimensional results for Poiseuille flow . . . . . . . . . . . . . . . . . . 71
6.2.1 Higher moment relaxation frequencies . . . . . . . . . . . . . . . . . 76
6.3 A dimensional problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
7 LB for Multi-Component Gases 80
7.1 Multi-component nomenclature . . . . . . . . . . . . . . . . . . . . . . . . 81
7.2 Binary channel flow theory . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
7.3 Multi-component Lattice Boltzmann scheme . . . . . . . . . . . . . . . . . 84
7.4 Mass conservation equation . . . . . . . . . . . . . . . . . . . . . . . . . . 85
7.5 Momentum equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
7.5.1 Mixture Centred Equilibrium (MCE) . . . . . . . . . . . . . . . . . 87
7.5.2 Species Centred Equilibrium (SCE) . . . . . . . . . . . . . . . . . . 90
7.5.3 Corrected Species Centred Equilibrium (CSCE) . . . . . . . . . . . 92
7.6 Expressions for diffusion slip . . . . . . . . . . . . . . . . . . . . . . . . . . 93
7.6.1 Implementing diffusion slip in the LB method . . . . . . . . . . . . 94
7.6.2 Intersections between walls and open boundaries . . . . . . . . . . . 96
7.7 Modified relaxation times . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
7.8 Self-diffusion test case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
7.9 A dimensional test case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
7.9.1 Alternative boundary conditions . . . . . . . . . . . . . . . . . . . . 107
7.10 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
8 Improvements to Multi-component Methods 110
8.1 Second order scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
vii
8.1.1 Crank-Nicolson integration . . . . . . . . . . . . . . . . . . . . . . . 110
8.1.2 Diffusion slip boundary . . . . . . . . . . . . . . . . . . . . . . . . . 112
8.1.3 Other boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
8.1.4 Intersection between a wall and another boundary . . . . . . . . . . 115
8.1.5 Binary channel flow results . . . . . . . . . . . . . . . . . . . . . . . 117
8.2 Pre-conditioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
8.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
9 Calculations of Binary Flow 121
9.1 Fully developed binary channel flow . . . . . . . . . . . . . . . . . . . . . . 122
9.2 Binary flow with one inert component . . . . . . . . . . . . . . . . . . . . . 130
9.2.1 The Stefan tube problem . . . . . . . . . . . . . . . . . . . . . . . . 130
9.2.2 A LB model of a Stefan tube type problem . . . . . . . . . . . . . . 131
9.2.3 The effect of channel width . . . . . . . . . . . . . . . . . . . . . . 136
9.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
10 Conclusions 141
10.1 Summary of results from this work . . . . . . . . . . . . . . . . . . . . . . 141
10.2 Suggestions for future research . . . . . . . . . . . . . . . . . . . . . . . . . 143
References 146
viii
List of Symbols
The page on which the following symbols are first used is listed.
A . . . . . . . . . . . 4
A . . . . . . . . 111
α . . . . . . . . . . . 4
B . . . . . . . . . . . 4
β . . . . . . . . . . . 4
Bσ,ς . . . . . . . . . 82
c . . . . . . . . . 33
Dσ,ς . . . . . . . . . 28
e . . . . . . . . . 11
Kn . . . . . . . . . 35
Ma . . . . . . . . . 35
f . . . . . . . . . . . 6
gσ . . . . . . . . 111
γ . . . . . . . . 118
H . . . . . . . . . 68
i . . . . . . . . . . . 4
k . . . . . . . . . 89
κσ . . . . . . . . . 95
κ¨σ . . . . . . . . 113
L . . . . . . . . . 34
λ . . . . . . . . . . . 7
Λ . . . . . . . . . 44
λδ . . . . . . . . . 86
λν . . . . . . . . . 44
λξ . . . . . . . . . 44
M . . . . . . . . . 81
m . . . . . . . . . . . 6
M . . . . . . . . . . . 6
Mσ . . . . . . . . . 28
mσ . . . . . . . . . 82
nσ . . . . . . . . . 28
ν . . . . . . . . . 50
Nx . . . . . . . . 105
Ω . . . . . . . . . . . 7
p . . . . . . . . . 13
Φσ . . . . . . . . 112
Π . . . . . . . . . 36
ρ . . . . . . . . . . . 6
ρ0 . . . . . . . . . 35
RT . . . . . . . . . 39
R0 . . . . . . . . . 81
σ . . . . . . . . . . . 4
ς . . . . . . . . . . . 4
S(Mσ) . . . . . . . . 93
t . . . . . . . . . . . 6
tf . . . . . . . . . 34
Θσ . . . . . . . . . 97
Θ¨σ . . . . . . . . 116
U . . . . . . . . . 34
u . . . . . . . . . . . 6
u¨ . . . . . . . . 111
x . . . . . . . . . . . 4
ξ . . . . . . . . . 50
Xσ . . . . . . . . . 81
y . . . . . . . . . . . 4
Yσ . . . . . . . . . 81
ix
Chapter 1
Introduction
Solid Oxide Fuel Cells (SOFCs) are being developed by Rolls Royce Fuel Cells, with the
hope that they will replace gas turbines in stationary power applications. The study of
fuel cells involves many engineering problems which are not yet well understood, one of
which is the diffusion of the fuel through the electrodes of the cell. The porous electrodes
have irregular geometries, and currently no satisfactory method exists for modelling the
gas flow through them. This work is the first step towards developing a Lattice Boltzmann
(LB) model suitable for that task.
The origin of the Lattice Boltzmann method is the work in the 19th century by Maxwell
and Boltzmann on the kinetic theory of gases. The first success of the kinetic theory was
Maxwell’s prediction that the viscosity of a gas was independent of its density, which was
later experimentally confirmed. Chapman and Enskog, working independently in 1916,
obtained a class of slowly varying or ‘normal’ solutions (see Grad (1963)) to the Boltz-
mann equation and found expressions for the transport coefficients which again showed
good agreement with experiments. The work of Bhatnagar, Gross and Krook significantly
simplified the Boltzmann Equation, yet key results could still be obtained with sufficient
accuracy. However most of kinetic theory remains unconfirmed, as measurements of the
velocity of individual molecules of a gas are currently infeasible. It is also a subject which
is particularly labour intensive to study, Chapman reportedly described reading his own
1
book as ‘like chewing glass’ (According to Brush (1976) citing the Observer newspaper
from 1957).
The concept that lies behind Lattice Boltzmann (and before that the Lattice Gas Au-
tomaton) is to store the velocity of every particle that makes up a gas, and track the
particles as they collide with each other. Of course, to store every particle would be far
too computationally expensive, and so simplifications are made. This is where the two
ideologies of Lattice Boltzmann split. On the one hand, efforts can be made to ensure
the computational particles behave as closely as possible to the behaviour predicted by
the kinetic theory. On the other hand, the particle’s behaviour is not of concern, as long
as the behaviour of certain averages (called ‘moments’) such as density and momentum
agree with the results from kinetic theory. Both viewpoints have merit, but they often get
confused and mixed together, as the discussions on boundary conditions in Section 5.1.1
and mixture collision integrals in Section 3.1.6 will show. This work considers only the
behaviour of the moments, and the particles in the computational scheme do not behave
like real particles. An introduction to the way Lattice Boltzmann models are constructed
is given in Chapter 2, and then a review of the Lattice Boltzmann literature relevant to
this work is in Chapter 3.
Aside from Lattice Boltzmann this work also covers the implementation of a diffusive
slip boundary condition for a binary gas mixture. Diffusion slip has an interesting history,
but one which has mostly been ignored. There is a very simple reason for this, in most
applications diffusion slip itself can be safely ignored. However in fuel cells and some
other applications this is not the case as the flow through the porous material is driven
by concentration gradients not total pressure gradients, and diffusion slip becomes a key
factor. The story of diffusion slip, and some more details about its kinetic origin, are in
Chapter 3.
The main difficulty in extending the standard Lattice Boltzmann model to multi-
component gas mixtures is the non-physical equation of state commonly used. Therefore,
in Chapter 4, a model for a compressible single component fluid with a barotropic equation
of state is developed. This method can easily be extended to gas mixtures, but is not the
2
most efficient tool for single component gases. The boundary conditions commonly used
are based on specifying the behaviour of the particles, and not the moments. However,
for diffusion slip boundary conditions to be implemented a set of conditions on the mo-
ments must be imposed. Therefore, a general new approach to boundaries is conceived
in Chapter 5. This approach is not specific to any problem and represents a significant
shift in the way boundary conditions for Lattice Boltzmann are considered, and sheds new
light on established boundary condition methods. Both the new equation of state and new
boundary treatment are validated in Chapter 6 by simulating Poiseuille flow through a
channel.
The extension to multi-component gases, building on the foundations of the preced-
ing chapters, is described in Chapter 7. A key capability of the model developed is the
separation of diffusive and viscous effects, allowing flow with any Schmitt number (the
ratio of viscous diffusion to mass diffusion, ν/Dσ,ς) to be modelled. Diffusive slip bound-
ary conditions are developed, and binary channel flow is simulated and validated against
a simple one-dimensional scheme for solving the macroscopic flow equations. To obtain
accurate results in reasonable computational time improvements to the numerical scheme
are required, as described in Chapter 8.
Finally, in Chapter 9 two binary problems are investigated, channel flow and the Stefan
tube, and some interesting behaviour of binary gases is revealed. These problems highlight
the capability of the method and its potential future applications.
1.1 Notes on style
The principle complaint amongst engineers reading the LB literature is the use of non-
dimensional variables. In this work there are two scales by which variables may be non-
dimensionalised, indicated with a˜or aˆover the variable. However, for clarity in Chapters
2 and 5 the˜is dropped, so that a variable without a˜in those chapters is assumed to be
non-dimensional.
Another feature common to LBM is the use of the momentum flux density tensor,
3
following the approach of Landau and Lifshitz (1959). The symbol Π is used for all moments
of the velocity distribution function f , see Section 4.1.3. Παβ (i.e. Π with two subscripts)
is the momentum flux density tensor, and consists of the second order moments of f . The
momentum flux density tensor is composed of its equilibrium values minus the viscous
stress tensor. In this work the viscous stress tensor is −Π(neq)αβ , and so Παβ = Π(eq)αβ + Π(neq)αβ .
The subscript i indicates an element of a vector. For example fi is element number i of
f where the use of bold font indicates a vector. Matrices are represented by the use of a
bold capital letter, for example M .
The subscript σ is used to identify any species in a mixture, then ς identifies the other
species. These are replaced by the subscripts A and B when the expression would not be
true for any species, but applies only to specific species.
Similarly, the subscript α represents any Cartesian direction, and β represents the other
direction. All examples are two dimensional. The subscripts x and y are used when
referring to a direction in particular.
The superscript (eq) is used to indicate the variable is evaluated at equilibrium, the su-
perscript (neq) indicates the non-equilibrium component of that value, for example f (neq) =
f − f (eq).
The superscript † is used to indicate a variable after the collision step, while the super-
script + indicates a variable after the collision and streaming step.
The ′ superscript indicates that the value of the variable has been corrected to include
certain numerical effects, as described in Section 5.2. The γ superscript indicates the vari-
ables used in conjunction with pre-conditioning, described in Section 8.2. The ∗ superscript
is used to describe multi-component collision operators in Section 3.1.6, and is not used
again.
Double dots, like ¨ , indicate a variable derived from the variable g, introduced and
explained in Chapter 8.
4
Chapter 2
An Introduction to LB Methods
This section is intended as an introduction to general LB methods. In order to keep the
chapter concise the example chosen is particularly basic. The solution to a one-dimensional
problem will be examined, a fluid subject to a ‘square wave’ initial density profile in an
infinite domain. As time progresses the fluid in the region of initially high density disperses,
until the density of the fluid everywhere is the same and the velocity is zero. A periodic
problem was chosen so that boundary conditions need not be considered.
For ease of explanation all variables in this chapter are non-dimensional.
2.1 A one-dimensional scheme (D1Q2)
Figure 2.1: The D1Q2 grid. The arrows show the velocities associated with each element of f .
If the variables of interest are the density and velocity of the fluid, at least two values
must be stored at each grid point. To use a LB method the situation can be considered in
this way: At any grid point let the mass of particles per unit volume having unit velocity
to the left be f1 and the mass of particles per unit volume having unit velocity to the
5
right be f2. Then, the mass density is (f1 + f2) and the net mass flux to the right is
(f2− f1). If an equal number of particles are travelling in each direction then the net mass
flux (equivalently, momentum) is zero. The mass density of particles is represented as the
variable f(i, x, t), where i is a directional index, in this example i = 1 points left and i = 2
points right as shown in Figure 2.1. x is the spatial position of the node and t is time. At
each grid point the particle mass density can be expressed in vector form, and a vector of
the macroscopic variables, m, is defined as,
fD1Q2 =
 f1
f2
 mD1Q2 =
 ρ
ρu
 (2.1)
where the subscript of the scalar f is the directional index i. The mass density is ρ and
the macroscopic net velocity is u. The shorthand ‘DaQb’ is used for describing a Lattice
Boltzmann scheme, where ‘a’ is the number of dimensions and ‘b’ is the number of discrete
velocities. The vector m contains the moments (weighted summation) of f , and this
relationship can be expressed mathematically using a transform matrix, M , as follows,
mD1Q2 =
 1 1
−1 1
 f1
f2
 = MD1Q2fD1Q2 (2.2)
There are two different bases in which information about the fluid can be expressed. The
‘moment’ basis contains quantities which are weighted functions of the particle f values,
whereas the ‘particle’ basis (or the ‘f ’ values) contains values for the mass of particles
travelling in each direction. Clearly the particle basis is well suited to computation, whereas
the problem is defined in terms of the moment basis. Switching between the two is the key
concept in analysing LBM.
Having defined the representation of the fluid in terms of the two populations at each
grid point, the unsteady behaviour of the fluid will now be examined. Each time step is split
into two distinct processes, namely streaming and collisions. The streaming process will
involve all the population’s (the ‘f ’ values) moving to the next grid space in the direction of
6
that population’s discrete velocity. Each node will effectively swap populations with each
of its nearest neighbours. This is the principle method of communicating disturbances
through the grid, and the speed of these exchanges must be much higher than the bulk
velocity of the fluid. After the streaming is completed, there is a local collision step. This
may be thought of as those particles arriving at each grid point physically colliding, and
exchanging momentum, with each other. In fluid dynamics problems the total mass at
each grid point must be conserved during the collision step. The modelling of this collision
step is vital in ensuring the correct macroscopic behaviour.
Before modelling the collision step an equilibrium solution must be considered. Equi-
librium occurs when the spatial gradients of velocity are zero; in this specific case that will
be when there is zero velocity everywhere. Equilibrium values for the populations for each
discrete velocity in terms of an equilibrium density may be defined, using Equation 2.2, as,
f (D1Q2,eq) = (MD1Q2)−1
 ρ(eq)
0
 =
 ρ(eq)2
ρ(eq)
2
 (2.3)
To model the collision process, a linear relaxation to the equilibrium state is used. That is
to say that when the particles from neighbouring sites collide, they do so in such a way as to
move closer to an equilibrium state. The rate at which this relaxation occurs is controlled by
a parameter λ. The model for this relaxation process is taken from the Bhatnagar-Gross-
Krook (BGK) [Bhatnagar et al. (1954)] equation commonly used in continuum kinetic
methods. The BGK equation states that, in a collision, particles within a population
all having the same velocity leave that population at a rate proportional to the current
population, whereas those particles entering that population do so at a rate proportional
to the equilibrium value of that population. The collision term is written as Ω(f) and the
post collision state as f †.
fD1Q2,† = fD1Q2 + Ω(f) = fD1Q2 + λ(fD1Q2,eq − fD1Q2) (2.4)
Remembering that the density is defined as the sum of the elements of fD1Q2, the vectors
7
in the previous equation may be summed over their elements to give,
ρ† = ρ+ λ(ρeq − ρ) (2.5)
It has already been specified that the collision process must conserve mass and so ρ† = ρ.
Therefore ρeq = ρ. Equivalently
∑
i Ω(f)i = 0.
2.2 D1Q2 example code
The MATLAB code for the D1Q2 example outlined above is included here. Initial, and
arbitrary, values are defined as follows:
nx=31; %number of grid points
nt=30; %number of time steps
M = [1,1;-1,1];
lambda =1/3;%This controls the rate of relaxation
%initialise values of f
for x=1:round(nx/3)
f(1,x,1)=0.5;
f(2,x,1)=0.5;
end
for x=(round(nx/3)+1):(round(2*nx/3)-1)
f(1,x,1)=1;
f(2,x,1)=1;
end
for x=round(2*nx/3):nx
f(1,x,1)=0.5;
8
f(2,x,1)=0.5;
end
The collision and stream rules are applied at every time step. At each collision step the
density is calculated from the existing fluid populations, then the equilibrium function is
calculated using this density. Finally the f values are relaxed towards equilibrium. The
stream step moves the f values to the nearest neighbour, but also moves them to the next
time step, to save a separate copy operation. The node points at each end must be dealt
with separately, and a periodic condition has been enforced in this particular example for
simplicity.
for t=1:nt
%Collide
for x=1:nx
rho(x) = f(1,x,t)+f(2,x,t);
feq(:,x) = [0.5*rho(x), 0.5*rho(x)];
f(:,x,t)=f(:,x,t)+lambda.*(feq(:,x)-f(:,x,t));
end
%stream
for x=2:nx
f(1,x-1,t+1)=f(1,x,t);
end
for x=1:(nx-1)
f(2,x+1,t+1)=f(2,x,t);
end
%periodic boundary
f(2,1,t+1)=f(2,nx,t);
f(1,nx,t+1)=f(1,1,t);
end
9
The results of this calculation, the density and velocity, are shown in Figure 2.2. The
behaviour is not as expected from a viscous fluid. The total fluid in the domain, the sum
over the nodes of the density remains constant over time at 41 non-dimensional units.
The scheme therefore conserves mass. However, there is no reduction in the magnitude of
the peak velocity of the fluid, and even over much larger time intervals the peak velocity
remains the same so that no steady-state is reached. Therefore, it appears that the viscous
behaviour of the fluid is not correctly modelled, a conclusion which is supported by the
analysis in the subsequent section.
t = 20
u
t = 15
u
t = 10
u
t = 5
u
t = 1
u
ρ
t = 20
ρ
t = 15
ρ
t = 10
ρ
t = 5
ρ
t = 1
10 20 30
10 20 30
10 20 30
10 20 30
10 20 30
10 20 30
10 20 30
10 20 30
10 20 30
10 20 30
−0.5
0
0.5
−0.5
0
0.5
−0.5
0
0.5
−0.5
0
0.5
−0.5
0
0.5
1
1.5
2
1
1.5
2
1
1.5
2
1
1.5
2
1
1.5
2
Figure 2.2: Density and velocity profiles at different time steps for D1Q2 example. Periodic
boundary conditions were used. Non-dimensional units.
10
2.3 Analysing the D1Q2 scheme
An analytical expression for the system outlined above is now developed. Combining the
streaming and collision step yields an update equation:
fD1Q2i (x+ e
D1Q2
i , t+ 1) = f
D1Q2
i (x, t) + λ(f
D1Q2,eq
i (x, t)− fD1Q2i (x, t)) (2.6)
Where eD1Q2i is introduced to represent the velocity with which f
D1Q2
i streams. Each
element moves one grid space per time step, and in this example eD1Q2 = [−1 1]T . A
truncated Taylor expansion,
fD1Q2i (x+ e
D1Q2
i , t+ 1) = f
D1Q2
i (x, t) +
∂fD1Q2i
∂t
+ eD1Q2i
∂fD1Q2i
∂x
(2.7)
is used to transform Equation 2.6 into a partial differential equation,
∂fD1Q2i
∂t
+ eD1Q2i
∂fD1Q2i
∂x
= λ(fD1Q2,eqi − fD1Q2i ) (2.8)
Only the first order terms in the Taylor expansion have been included, and there is an
error associated with this, but for this simple example this is ignored.
Summing the elements of the vectors in this equation, noting that
∑
i eifi = ρu, gives,
∂ρ
∂t
+
∂(ρu)
∂x
= λ(ρ− ρ) = 0 (2.9)
This is the mass conservation equation. Multiplying Equation (2.8) by eD1Q2i and summing
the vector elements, noting
∑
i e
D1Q2
i e
D1Q2
i f
D1Q2
i = ρ, yields,
∂(ρu)
∂t
+
∂ρ
∂x
= λ(ρequeq − ρu) = −λρu (2.10)
It is now clear that this simple scheme does not recover the correct momentum equation. It
is worth exploring the reasons for this. The initial restriction of only two discrete velocities
means that there can be only two meaningful macroscopic variables per grid space. The
11
Figure 2.3: The D1Q3 grid. The velocity of the particles associated with f2 is zero.
momentum equation requires a third macroscopic quantity, the momentum flux density.
Therefore, the error in the left hand side of the momentum equation can be attributed
to having too few lattice velocities. The error on the right hand side of the momentum
equation is due to the choice of using a BGK style collision model and the choice of f (eq).
There are an infinite number of ways of modelling the collision process, and each one will
give different macroscopic equations. If the collision process were no longer local to each
grid node it would be possible to compute the derivatives of the macroscopic variables
and perhaps recover the correct momentum equation. However this would sacrifice the
simplicity of the current code, and move towards a scheme which is more representative of
traditional CFD methods.
2.4 One dimension, three velocity (D1Q3) scheme
To recover the correct momentum equation from the scheme, whilst keeping a local collision
operator, another particle velocity is introduced, in this case a population with zero velocity
as shown in Figure 2.3. The velocity vector becomes eD1Q3 = [−1 0 1]T , and the transform
matrix is,
MD1Q3 =

1 1 1
−1 0 1
1 0 1
 (2.11)
The new third moment was selected following the general structure that each row of M
is e transposed with each element raised to the power of the row number minus one.
For example row 3 is [(−1)2 (0)2 (1)2]. This approach is in keeping with traditional
kinetic theory, however there are different strategies for generating moment bases based
on orthogonalisation. Any third moment that is linearly independent of the first two
12
is acceptable. However, the choice made here reduces the algebraic complexity of what
follows. The new moment vector is,
mD1Q3 =

ρ
ρu
ΠXX
 (2.12)
Where ΠXX is the momentum flux in the x direction. The equilibrium state is
mD1Q3,eq =

ρ
ρu
p+ ρu2
 (2.13)
where the Euler equation was used to provide a value for ΠXX at equilibrium, and p is the
pressure of the gas. An equation of state is now required. For simplicity this is taken to
be p = 1/3ρ. The issue of equations of state will be discussed in considerable length later
on, and it is unhelpful to burden this introduction with that discussion.
Using the same technique described in Section 2.3 the same mass conservation equation
is recovered, however there is a different momentum equation,
∂(ρu)
∂t
+
∂ΠXX
∂x
= λ(ρequeq − ρu) = 0 (2.14)
The value for ΠXX is unknown, and so the final independent equation in the moment basis
is used, which may be obtained by multiplying the D1Q3 version of Equation (2.8) by
eD1Q3i e
D1Q3
i and summing over i.
∂ΠXX
∂t
+
∂(ρu)
∂x
= λ(Π(eq)xx − ΠXX) (2.15)
Where the fact that the magnitudes of the eD1Q3i are either unity or zero has been used to
give the result (eD1Q3i )
3 = eD1Q3i . The approximation that must be employed to arrive at
an expression for ΠXX is that it is not too far from its equilibrium value. Therefore the
13
terms on the left hand side of Equation (2.15) may be given their equilibrium values,
∂(p+ ρu2)
∂t
+
∂(ρu)
∂x
= λ(p+ ρu2 − ΠXX) (2.16)
Using the mass continuity equation, 2.9, the equation of state (p = ρ/3) and the fact that
at equilibrium ∂(ρu)/∂t = −∂(p+ρu2)/∂x, after some manipulation the expression for the
momentum flux in the x direction becomes,
ΠXX = p+ ρu
2 − 2
3λ
ρ
∂u
∂x
(2.17)
For details of exactly how this is reached, consult Dellar (2002b). It is only when the
equation of state p = ρ/3 is used that such a simple expression is obtained. Inserting this
expression into Equation (2.14) yields,
∂(ρu)
∂t
+
∂(ρu2)
∂x
= −∂p
∂x
+
∂
∂x
(
2ρ
3λ
∂u
∂x
)
(2.18)
which is the compressible one-dimensional Navier-Stokes equation with the viscosity given
by ν = λ/3. Graphs of the same problem solved with D1Q2 are shown in Figure 2.4.
The D1Q3 results show the peak velocity decreasing after 20 time step, and this should be
compared with the results from the D1Q2 model. The scheme reaches a steady state with
zero velocity and a constant density, as would be expected.
2.5 Concluding remarks on the 1D examples
A relatively simple code, D1Q2, was developed and although there are problems with the
momentum equation it is remarkable that the mass conservation equation is recovered.
Such a model could be used to model the transport of a passive scalar in an established
flow field. The introduction of one extra velocity was sufficient to solve the problem with
the correct momentum equation.
The one-dimensional code will now be abandoned in favour of a two-dimensional scheme
14
t = 20
u
t = 15
u
t = 10
u
t = 5
u
t = 1
u
ρ
t = 20
ρ
t = 15
ρ
t = 10
ρ
t = 5
ρ
t = 1
10 20 30
10 20 30
10 20 30
10 20 30
10 20 30
10 20 30
10 20 30
10 20 30
10 20 30
10 20 30
−0.5
0
0.5
−0.5
0
0.5
−0.5
0
0.5
−0.5
0
0.5
−0.5
0
0.5
1
1.5
2
1
1.5
2
1
1.5
2
1
1.5
2
1
1.5
2
Figure 2.4: Density and velocity profiles at different time steps with the D1Q3 model. Periodic
boundary conditions were used. Non-dimensional units.
15
with nine velocities. However, the method that will be used is based on the same principles
as the one-dimensional method developed here. By choosing the moment system wisely
the macroscopic equations will appear in conservation form, as they did in the example
above. With nine velocities per grid node it will be shown that the correct momentum
equation in two dimensions can be recovered, but not the energy equation.
16
Chapter 3
Literature Review
This chapter contains a review of the literature on the Lattice Boltzmann method and the
phenomenon of diffusive slip. With a topic as broad as Lattice Boltzmann even the most
extensive literature review would not be comprehensive, and hence only selected topics have
been examined. The literature review on diffusive slip is more comprehensive, however,
because the subject has not been so extensively researched.
3.1 Lattice Boltzmann Methods (LBM)
3.1.1 Lattice Boltzmann from the Lattice Gas Automaton
A search for ‘Lattice Boltzmann’ on the ‘Web of Knowledge’ journal search engine returns
4,941 results, an indication of the size of the following that LBM has enjoyed in recent years.
The earliest entries date back to 1989, when LBM emerged from the Lattice Gas Cellular
Automaton (LGCA or LGA) computational scheme. The intention of both methods was
to model the behaviour of individual particles, and thus model the macroscopic behaviour
of a collection of particles. In the LGA the domain was split into nodes, each of which
had a population of particles moving with its own discrete velocities. The key difference
between this and LBM is that the number of particles at each node with a certain velocity
is either one or zero; in LBM this could be any real number. The LGA was computationally
17
intensive, and the LB method is better suited to the architecture of most computers as
the local collision step may be easily parallelised. For problems where the Navier-Stokes
equations are not applicable, for example in the study of rarefied gas flow, there was no
established computational method. Since then, both LBM and Monte Carlo methods have
been developed.
The earliest review of LBM was conducted by Succi et al. (1991) and it is worthwhile
examining his work in some depth as it lays the foundation for many misconceptions that
follow. The description of the Lattice Boltzmann Equation (LBE) begins, ‘The LBE is a
nonlinear finite difference equation (even though it does not result from the discretisation of
any partial differential equation!)’. This is interesting as almost exactly the same method is
presented as the discretisation of the continuous Bhatnagar Gross Krook (BGK)equation,
Bhatnagar et al. (1954), in later work.
The Succi paper described a Lattice Boltzmann method which will not be examined
here as it has since been superseded. Results were then presented for laminar flow through
porous material, bifurcations of a two-dimensional Poiseuille flow and fully developed two-
dimensional forced turbulence. This general approach, presenting a model with very com-
plicated and specific applications, is common in the LBM literature. The scope of papers
generally far exceeds their depth. In the results presented by Succi for the flow through
a porous material there are several immediate problems. The pressure force is applied to
each grid site, meaning that a region of fluid almost entirely encircled by impermeable
material will have a high velocity. The boundaries are not examined in depth. A reader
more experienced in traditional computational fluid dynamics techniques would be left
wondering what time steps and grid spacing were used, how the method converged and
whether it was stable. Instead, what Succi presents is a series of particular cases where
the results look qualitatively impressive.
Most of the early work on LBM was dedicated to finding model expressions for the
collision term. As seen in Chapter 2, the collision term is essential in determining which
macroscopic equations the scheme describes. A particularly interesting paper on this sub-
ject is that of Chen et al. (1992). They correctly identified that if one has an LB model
18
with only one speed (not one velocity), an increase in the density of the fluid is equivalent
to an increase in the momentum flux (in Chapter 2 in the D1Q2 model
∑
i eieifi =
∑
i fi).
Therefore there must be at least one population having a velocity with a different speed,
the simplest being a zero velocity population. Unfortunately, the authors chose to cloud
the solution to this simple mathematical problem with talk of a ‘particle reservoir’. The
simplest form of the collision integral was used, a linear relaxation to equilibrium,
Ω(fi) = −λ(fi − f (eq)i ) (3.1)
It should be noted that Ωi is a general expression for the collision integral, such that
f+ = f + Ω(f). A general expression for the equilibrium distribution was also given,
and then constants within that expression were derived by consideration of mass and
momentum conservation, and by the correct form of the three viscous stress components.
This approach is misleading however, because the form of the general expression had
already imposed some symmetry constraints which were not explicitly mentioned. Finally,
they noted that there was still one degree of freedom in the definition of the equilibrium
function and so an arbitrary constraint was introduced. Nevertheless, the algorithm they
described is still widely used today in an essentially unchanged form.
3.1.2 The scope of Lattice Boltzmann
LBM, trading on its computational simplicity and stability, has been applied to a bewil-
dering array of problems: magnetohydrodynamics, immiscible fluids, compressible fluids,
incompressible fluids, nonlinear reaction-diffusion equations, thermo-hydrodynamics, non-
Newtonian flow, multi-phase, multi-components and short-time motion of colloidal particles
had all received the Lattice Boltzmann treatment as early as 1993. Although obviously
very different from each other, these papers shared a common theme, the use of LBM as a
substitution for the Navier-Stokes or similar equations, with some extra physics bolted on.
It is possible that this work does quickly, and in some cases successfully, simulate these
complex phenomena, which traditional CFD struggles to do, but it is beyond the scope of
19
the present work to study all the attempts in detail. The results of the multi-component
models are examined, however, particularly because the early models do not recover the
correct physics.
3.1.3 The rise of Lattice BGK
The most popular LBM ‘starter’ paper is that of Chen and Doolen (1998), with an as-
tounding 1,501 (as of July 2010) citations. This paper defines the ‘second generation’ of
LBM, abandoning the derivation from the lattice cellular automaton and instead rewriting
the history of LBM as a discretisation of the continuous BGK model. The two-dimensional
nine-velocity model (called D2Q9) is the standard model. The Lattice Boltzmann equation
became the now common Lattice BGK (LBGK) equation,
fi(α + eα,i, β + eβ,i, t+ 1) = fi(α, β, t)− λ[fi(α, β, t)− f (eq)i (α, β, t)] (3.2)
The equilibrium function which was previously so carefully constructed is now the contin-
uous Maxwellian distribution confined to the discrete velocity set using a Gauss-Hermite
quadrature. This change in derivation surprisingly does not lead to a significant change
in the value of the equilibrium function itself. Essentially, although the method has not
changed the interpretation is quite different and this is a dangerous step. It is no longer
clear that the method has evolved from considerations of the conservation of mass and
momentum in terms of macroscopic variables and, instead, it is viewed as being derived
purely from microscopic considerations.
The early papers of He and Luo (1997b,c,a) strongly supported the LBGK philosophy,
although Luo’s later papers highlighted shortcomings of the method, Luo (2004). More
recently he has turned to the Multiple Relaxation Time (MRT) model, see Lallemand and
Luo (2000). However, the damage had already been done in establishing the microscopic
interpretation of LBGK. Not that this was the fault of Luo, rather the lax referencing of
his early papers by other authors. Succi et al. (2002) reviewed LBGK in the context of
other microscopic lattice based methods.
20
Because of this implied microscopic nature of LBGK, efforts have been made to uncover
behaviour beyond Navier-Stokes, especially in micro-channel flow. Verhaeghe et al. (2009)
proved that all the attempts beforehand, and they listed 15 such papers, have been exercises
more in curve fitting than modelling the physics. They modelled Navier-Stokes flow with
Knudsen velocity slip at the wall. In order to enforce this slip they used an approach
mathematically similar to the approach in this work, but the explanation is confused by
the unnecessary use of microscopic boundaries, as will be discussed in Chapter 5.
This thesis, building on the work of Asinari and Ohwada (2009), hopefully puts beyond
doubt the fact that with a two-dimensional nine-velocity LBGK method it is not possible
to recover the rarefied gas equations beyond the Navier-Stokes level.
3.1.4 Boundaries for LBM
The confusion about the microscopic or macroscopic nature of LBM has probably had the
biggest impact on boundary conditions, a subject in which there is a lack of good quality
work. By far the most common approach to modelling solid-fluid boundaries is the so-
called ‘bounce-back’ rule, which is often said to be derived from physical considerations of
the particle collision with the wall; when the particle impacts with the wall its velocity is
reversed. In fact, there has never been a suggestion that this is how real particles behave
and this method originated from the LGA. When real particles impact with a wall they
reflect either in a specular way, with the component of their velocity normal to the wall
reversing, or in a diffusive way, returning from the wall in a random direction. Attempts
have been made to replicate both these processes in LBM, notably by Sofonea and Sekerka
(2005a).
Efforts have also been made to improve on the poor performance of the bounce-back
method. Ziegler (1993) showed that the results of the bounce-back boundary condition
more accurately represented a no-slip wall if the position of the wall was interpreted as
half way between the last fluid node and the first solid node, a view shared by He et al.
(1997). Another correction to the bounce-back method was proposed by Junk and Yang
21
(2005) who were particularly concerned with geometry that varied on a scale less than
one grid spacing. Inamuro et al. (1995) devised a method where the unknown distribution
functions at a wall were set to their equilibrium values with a counter slip-velocity such
that the overall velocity was zero.
For open (inlet and outlet) boundary conditions there is no simple solution to be found
in continuum theory, equivalent to bounce-back, and in some respects this has been an
advantage in forcing the invention of more creative schemes. The situation can be simply
stated by considering that at a straight boundary there are 3 unknown elements of f ,
and usually only two hydrodynamic constraints. Noble et al. (1995b,a) considered closure
methods for the extra unknown introduced at a boundary. Zou and He (1997) proposed
what is probably the most popular scheme, the so-called ‘non-equilibrium bounce-back’
method (this scheme will be examined in Chapter 5). Yu et al. (2005) tried to modify the
Zou and He method to reduce the reflection of pressure waves at the boundary.
For setting a non-zero velocity along the walls of a channel Latt et al. (2008) conducted a
review of several methods, and found no significant differences between their performances.
Hollis et al. (2008) used the viscous shear stress as a closure condition for boundaries, and
went on to try and deal with nodes at the intersection between intersecting straight bound-
aries, something conspicuously absent from most other papers. Ginzburg and d’Humieres
(2003) conducted a review of the different methods of reflecting the particles from a surface.
An interesting approach was taken by Skordos (1993), who used finite difference approxi-
mations of the derivatives of the hydrodynamic variables to recreate the values of the higher
moments at a boundary. Maier et al. (1996) and Chen et al. (1999) used extrapolation of
the known elements of f for the same purpose.
In general, schemes which extrapolate the boundary condition from the interior should
be treated, in the opinion of the author, with suspicion. The boundary conditions are an
essential part of the solution of partial differential equations, and they must be imposed
on the PDE, not the other way around.
22
3.1.5 Improvements on LBGK
The Multiple Relaxation Time (MRT) collision term, proposed initially by Higuera et al.
(1989) and later by Lallemand and Luo (2000), is an extension of the LBGK model. Each
moment is relaxed at a different rate, and the collision terms are given by,
Ω(f) = M−1ΛM(f eq − f) (3.3)
where Λ is a diagonal matrix of relaxation rates, explained in Section 4.3. This offers sev-
eral advantages over LBGK, see for example Dellar (2003) and Reis and Phillips (2008).
There are two different moment systems commonly used. Continuum kinetic methods use
a moment basis of increasing powers of the molecular velocity, as explained in Section 2.4.
The advantage of this basis is that it is easily extendable, and discrepancies between con-
tinuum kinetic theory and Lattice Boltzmann may be easily observed, as in Section 4.1.3.
However, Lallemand and Luo (2000) used a moment basis derived from a purely mathe-
matical background, such that the moments are orthogonal. This allows the identification
and isolation of non-hydrodynamic ‘ghost’ moments, at the expense of a more conceptually
complicated scheme. In terms of computational expense, it is not clear that either basis
offers an advantage.
The single relaxation time model may be seen as a special case of the MRT model with
all relaxation times set to the same value. Another special case is the two relaxation time
scheme proposed by Ginzburg et al. (2008).
The motivation for the development of the MRT approach was to reduce the effect that
higher order (unphysical) moments have on the lower order moments of interest. Latt and
Chopard (2006) developed a novel solution, instead of MRT, to this problem. After each
time step they ‘regularise’ f to remove any influence from the higher order moments. This
is computationally expensive however, and MRT offers a much neater solution. In fact, the
Latt and Chopard (2006) method is equivalent to MRT with a particular choice of collision
times that sets the higher moments to equilibrium after each collision, as in Higuera et al.
23
(1989) and McNamara et al. (1995).
3.1.6 LBM for multi-component fluids
In the modelling of mixtures of gases, there are two important differences from the single
component method. Firstly, the effect of the different molecular masses of the gases must
be taken into consideration. Secondly, the collision term must be modified to take into
account the momentum exchanged between the species.
The first problem proves difficult to overcome using LBGK, as conventionally a non-
physical equation of state is used and the grid spacing is proportional to the molecular
masses, as will be explained in Section 4.1.6.
McCracken and Abraham (2005) proposed a method of interpolating between two dif-
ferently spaced grids. However, this is too computationally intensive for real applications.
If one considers that the ratio between the molecular masses of a binary gas could be, say,
20, the grid spacing ratio would have to be the same. Therefore a 10× 10× 10 grid for one
species would require a 200× 200× 200 grid for the other species.
Yu and Zhao (2000) proposed a model that allowed a change in the equation of state
using a finite difference based modification to the pressure term in the momentum equation.
This was done with the use of an ‘attractive force term’ added to the collision term which
then appeared in the momentum equation. The modified equation of state did not appear
in the viscous terms however. Essentially, this was the first attempt at a model which
would be developed more successfully by Dellar (2002b), who removed the need for the
finite difference scheme. He developed a more robust model for generalised equations of
state, based on the shallow water equations but applicable to any barotropic equation of
state. Guo and Zhao (2003) applied essentially the same method to multiple gas mixtures.
To address the second problem it seems the natural choice to turn to continuum kinetic
theory as described by Hamel (1966) and Sirovich (1962). The collision term in both these
24
models may be grossly simplified, and written as,
Ωσ = Ω
σ,σ + Ωσ,ς (3.4)
Where Ωσ,σ represents self-collisions between particles of the same species, and Ωσ,ς rep-
resents cross collisions between particles of different species. This approach was adopted
by Sofonea and Sekerka (2001, 2005b) who realised that if BGK-style models are used to
describe the collision process, such as,
Ωσ = λ
σ,σ(fσ − f (eq),σ,σ) + λσ,ς(fσ − f (eq),σ,ς), (3.5)
then this expression may be rearranged into the form,
Ωσ = λ
∗(fσ − f (eq),∗) (3.6)
where the values with a superscript ∗ are functions of the self-collision and cross-collision
relaxation frequencies and equilibrium functions. Sofonea and Sekerka set f (eq),∗ to be a
function of the mixture velocity, the mass mean of the species velocities. The value of λ∗
may then be set to model diffusion correctly, but the viscosity recovered is incorrect.
To circumvent this shortcoming, Facin et al. (2004), Luo and Girimaji (2002) and Asinari
(2005) separately proposed ‘three relaxation time’ models, of the form of equation (3.5).
Arcidiacono et al. (2006b,a, 2007) also proposed a similar model, but derived from an
entropic point of view. Asinari (2006) demonstrated that when these models are written in
the form of Equation (3.6) the equilibrium function does not have a simple form, implying
that when a mixture composed of identical components is considered, it does not behave
in the same way as a single gas.
Asinari (2008, 2009) developed a lattice Boltzmann model based on the Andries et al.
(2002) (AAP) model for gas mixtures. Two vital points are established, the first is that
the diffusive behaviour of the gases can be modelled by an exchange of momentum between
the species. The second is that the MRT method can be used to apply different relaxation
25
times to the diffusive behaviour and the viscous behaviour. In this manner the diffusivity
and viscosity can be made entirely independent of each other.
To summarise, the differentiation between self and cross collisions is appropriate for
continuous kinetic theory. However, it is difficult to separate the diffusivity and viscosity
of the gases in a LB model from this starting point. It is better to start from a moment
basis, where diffusivity and viscosity are already separated.
3.1.7 LBM in porous media and solid oxide fuel cells
A review of the early work using LB for flow in porous materials can be found in Chen
and Doolen (1998). Ahrenholz et al. (2006) provide an insight into how more complicated
porous media geometries may be reconstructed in a computational model.
Several people have attempted to model flow through the porous material of a Solid
Oxide Fuel Cell (SOFC). Asinari et al. (2007) used electron microscopy to image the SOFC
electrolyte and then took average properties of grain size to reconstruct the geometry.
Joshi et al. (2007b,a) thought that they modelled the non-continuum diffusion equations
in porous media and SOFCs, despite only using a D2Q9 scheme. Both authors used bounce-
back boundary conditions and, most importantly, did not refer to diffusion slip at the solid
walls or even the slip between the species. The scope of these papers is very large, and it
is difficult to validate their results quantitatively.
Wang and Afsharpoya (2006) considered much larger channels in SOFCs, and claimed
to predict turbulence at high Reynolds numbers, again using a bounce back wall condition.
Although not related to LBM, interesting work was done by Lee et al. (2002, 2003, 2005)
on solid oxide fuel cell geometry and composition, although no attempt was made to model
flow. Instead, permeability was deduced purely from geometric considerations.
Most recently, Rama et al. (2010) conducted a study of polymer electrolyte fuel cells,
and integrated the flow computations with the electrochemistry. Unfortunately, they used
the Luo and Girimaji (2002) model, in which each species’ grid spacing is different. The
extra interpolation needed because of this adds an unnecessary computational overhead.
26
Furthermore, there was no detailed consideration of fluid-solid boundaries, and the bounce-
back boundary condition was used without justification or explanation.
3.2 Diffusion slip
The unusual history of diffusion slip is now described. The story is characterised by
painstaking work over many years and a reluctance by the wider engineering community
to accept clear experimental evidence. This should be contrasted with the often rushed
work of the LBM community and their eagerness to accept new theories despite an absence
of evidence.
The story begins with the work of Sir Thomas Graham, a Scottish chemist who lived
between 1805 and 1869. His original papers, Graham (1876), although fascinating, are
difficult to read now as concepts such as viscosity had not been formalised. However,
Mason and Evans (1969) and Mason et al. (1969) have produced a modern summary,
which should be considered essential reading for any student of this subject. Graham was
responsible for discovering the relation which has become known as ‘Graham’s law’; when
two gases diffuse into each other, under a condition of constant total pressure (defined as
the sum of the partial pressures), the ratio of their molecular fluxes is inversely proportional
to the square root of the ratio of their molecular masses. In the experiments conducted by
Graham, a stucco plug was used to separate the two gases. Over time, the results of these
experiments became confused, and it was wrongly assumed that the pores in the stucco
must have been very small so that the result was due to non-continuum effects. In fact,
stucco has rather large pores, and the gas flows were well within the continuum regime.
The scale of the misunderstanding regarding the zero total pressure drop problem is
eloquently outlined by Mills (2003). If a diffusing gas mixture in a channel has constant
total pressure, the mass mean velocity is not zero. Almost all standard text books claim
the opposite, however, including Bennett and Myers (1982), Welty et al. (1984), White
(1988), Holman (1988), Geankoplis (1993), Baehr and Stephan (1998), Incropera et al.
(2001) and Cengel (1998). This is an astounding error in the history of the analysis of
27
diffusing gas mixtures.
The discrepancy has been highlighted and attempts have been made to address it from
both a mathematical and an experimental point of view. Dullien and Scott (1962) and
Evans et al. (1961b) both confirmed mathematically that Graham’s law applies for both
continuum and non-continuum conditions in a channel. Excellent experiments were carried
out by Knaff and Schlunder (1985), Hoogschagen (1955), Remick and Geankoplis (1973)
and Evans et al. (1961a, 1962a,b, 1963). All these experiments confirmed Graham’s law
over a wide range of pressures.
The reason for at least some of the confusion is due to the treatment of the solid wall
boundary condition. Kramers and Kistemaker (1943) proposed that the mass mean velocity
of a diffusing gas mixture is not zero at the wall. Instead they identified a ‘diffusion slip’
of the mixture directed, in the case of a binary gas, away from the higher concentration
of the heavier gas. The magnitude of the diffusion slip velocity derived by Kramers and
Kistemaker (1943) is,
uslip =
(
MA −MB
nAMA + nBMB
−
√
MB −
√
MA
nA
√
MA + nB
√
MB
)
DAB
dnA
dx
(3.7)
where uslip is the mixture slip velocity in the x direction. In this equation the light gas is
A and the heavy gas is B, Mσ is the molar mass of gas σ, nσ the number of molecules per
unit volume of gas σ, DAB the diffusion coefficient and x is the distance along the channel.
The expression for the diffusion slip velocity was extended to multiple gas mixtures by
Jackson (1977) and studied in some depth by Noever (1990a,b,c), who referred to it as as
the ‘baroeffect’.
More recently, Young and Todd (2005) and then Mills (2007) revisited the topic, but
otherwise it remains relatively unexplored. In any application which has a small or zero
total pressure change, but large changes in gas composition, correct modelling of diffusion
slip at solid wall boundaries is essential.
28
Slip Velocity
Ve
loc
ity
Distance from wall
Figure 3.1: Species and mixture velocities tangential to the wall, close to a solid wall. A charac-
teristic mean free path length is l, while a characteristic grid spacing may be ∆x. The flat velocity
profiles away from the wall indicate that viscous effects are small compared to the diffusive effect.
The slip velocity of the mixture is predominantly due to diffusion slip.
3.2.1 Origin of diffusion slip
The question of exactly how diffusion slip occurs is now examined in more depth. Figure
3.1 shows the velocity profiles of a binary gas near a wall. In the bulk of the flow the total
pressure gradient is small and a concentration gradient parallel to the wall and which is
constant across the channel causes diffusion of the two species in a direction parallel to
the wall. Some suitable average mean free path l is defined such that the particles less
than a distance l from the wall collide with the wall without making any other collisions.
In the region between the wall and a distance l from the wall the continuum equations,
Navier-Stokes and Stefan-Maxwell, are not applicable. Unless this region is to be modelled
correctly with the appropriate rarefied gas equations, boundary conditions appropriate to
the continuum modelling level must be applied. Looking again at Figure 3.1, it is suggested
that the appropriate continuum boundary conditions at the wall are indicated by the dotted
lines. Finally, although the grid spacing ∆x drawn on Figure 3.1 is larger than l, even if
∆x < l the no-slip condition would not be appropriate for the model developed in this
work, as the governing equations cannot resolve physical effects on that scale.
29
Slip Velocity
Ve
loc
ity
Distance from wall
Figure 3.2: Species and mixture velocities tangential to the wall, close to a solid wall. A char-
acteristic mean free path length is l, while a characteristic grid spacing may be ∆x. The total
pressure gradient is large, therefore the slip velocity at the wall is predominantly due to Knudsen
slip (which can be ignored for continuum flow).
Figure 3.2 shows a binary gas subject to Knudsen slip, and is included for comparison
with diffusion slip. The two are similar, however there is one very important distinction,
Knudsen slip has a large effect when the pressure of the gas is so low, or the channel is so
narrow, that the mean free path of the gas is no longer negligibly small. In such a case the
Knudsen number (l/H where H is the channel width) is greater than, say, 0.01. Diffusion
slip has a large effect when the diffusion fluxes of the species are much greater then the
convective fluxes. Both types of slip occur at all times, however for low Knudsen numbers
Knudsen slip can be ignored and for flows with high pressure gradients diffusion slip can
be ignored.
The argument of Kramers and Kistemaker (1943), for a Maxwell gas with perfect diffu-
sive reflection of molecules from the wall, is essentially as follows. The number of molecules
of each gas passing through a surface parallel to the concentration gradient is nσ c¯σ/4, where
c¯σ is the mean molecular velocity in the gas and nσ is the number of molecules. The mo-
mentum carried through this surface is therefore nσ c¯σmσuσ/4, where mσ is the mass of a
30
molecule of species σ and uσ is the velocity of the species σ. If the surface is a stationary
physical wall, the molecules returning from the wall have zero average velocity, and zero
average momentum. Therefore, the total momentum crossing the surface is,
nAc¯AmAuA + nB c¯BmBuB = 0 (3.8)
From the equipartition of energy
√
mAc¯A =
√
mB c¯B, and Equation (3.8) becomes,
ρAuA√
MA
= − ρBuB√
MB
(3.9)
Loyalka (1971) extended this analysis to include the effect of different intermolecular force
laws and gas-surface interactions. When the species have identical molecular masses, the
analysis of Kramers and Kistemaker (1943) predicts zero slip, although a small slip velocity
actually occurs due to these additional considerations. For species with very different
molecular masses however, the Kramers and Kistemaker (1943) expression gives reasonably
good agreement with Loyalka’s full expression. This analysis was extended to a higher level
of accuracy by involving a second order Chapman-Enskog expansion by Ivchenko et al.
(2002). They found the diffusion slip coefficient was then highly sensitive to the model
used for intermolecular forces.
A study of the layer adjacent to the wall was conducted using the numerical technique of
molecular-dynamics by Mo and Rosenberger (1991). The results were, however, inconclu-
sive. More recently, Kim et al. (2009) also studied this subject with both Direct Simulation
Monte Carlo and D2Q9 and D2Q16 Lattice Boltzmann models. The results reported were
very limited and did not recover Graham’s law. Clearly, further work in this area is needed.
To the best of the author’s knowledge, no attempt has been made to model macroscopic
diffusive slip using the Lattice Boltzmann method.
31
Chapter 4
Single-Component Lattice Boltzmann
In this chapter, a kinetic-style scheme to solve the single component isothermal compress-
ible Navier-Stokes equations is derived. The nine velocity two dimensional (D2Q9) grid is
used. To make this model easily extendable to multi-component gas mixtures the correct
Ideal Gas equation of state is used. This is not the norm in LBM and it does mean that
a time accurate solution cannot be obtained, as will be demonstrated. However, unless a
physical equation of state is used when multi-component gases are considered, the effect
of the different molecular weights of the species cannot be incorporated into the model.
The system of moments is chosen, and then the equilibrium distribution function is
defined using the principles of mass and momentum conservation. The higher equilibrium
moments do not enter the mass and momentum equations and are selected on the grounds
of stability. The results of Dellar (2002b) are used to ensure this. The relaxation process
using the ‘Multiple Relaxation Time’ (MRT) model is investigated. An expression for the
momentum flux density tensor including viscous terms leads to the momentum equation,
which turns out to be the steady-state compressible Navier-Stokes equation.
The results of this chapter are not new, but the presentation and interpretation are
unique. Chapman and Cowling (1970) introduced a very successful method of analysis for
continuum kinetic theory, and this has been used with some success in Lattice Boltzmann
methods. The principle idea is that the velocity distribution function is expanded in
32
terms of a small parameter, such that f = f (0) + f (1) + f (2)..., and the time and space
derivatives are similarly expanded. By substituting the expanded velocity distribution
function into the Boltzmann equation, taking moments, and equating terms of similar
powers in , equations for the macroscopic variables are recovered. The so-called zeroth
approximation represents a flow at equilibrium. The first approximation leads to the
Navier-Stokes equations and subsequent approximations lead to the Burnett equations
and increasingly rarefied flow.
When the Boltzmann equation is represented in terms of just nine particle velocities the
use of the Chapman-Enskog expansion is unnecessary. With such a limited set of particle
velocities the third order moments are not independent, and effects past the Navier-Stokes
level cannot be recovered. The parameter  used in the Chapman-Enskog analysis is
meaningless and instead will be replaced with small parameters (Ma and Kn) defined
by the computational grid spacing and time step used to simulate the physical problem.
These parameters should not be confused with the physical Mach and Knudsen numbers,
although, in general, the literature does not make this distinction and can therefore be
confusing.
4.1 The D2Q9 grid and the Boltzmann equation
4.1.1 The D2Q9 grid
This work is limited to the D2Q9 grid shown in Figure 4.1, with velocities given by;
ex = c
(
0 1 0 −1 0 1 −1 −1 1
)T
(4.1)
ey = c
(
0 0 1 0 −1 1 1 −1 −1
)T
(4.2)
where c is the speed at which the values f are moved across the lattice, the so-called
‘grid speed’. Models can be constructed with any number of velocities, but D2Q9 has the
advantages of being compatible with a Cartesian grid and having independent zeroth and
33
Figure 4.1: The D2Q9 grid. The velocity associated with f0 is zero.
second order moments (like D1Q3 but unlike D1Q2). The velocity distribution function is
expressed as a vector of the nine components,
f =
(
f0 f1 f2 f3 f4 f5 f6 f7 f8
)T
(4.3)
As seen in Chapter 2 the Lattice Boltzmann method, when analysed with the use of a
Taylor expansion, leads to a partial differential equations of the form,
∂f
∂t
+ ex
∂f
∂x
+ ey
∂f
∂y
= Ω(f) (4.4)
In this case the resulting equation can be compared with the continuous Boltzmann equa-
tion, indeed the form is the same. However, such severe restrictions on f exist that
the comparison is almost meaningless. In the continuous Boltzmann equation, Ω(f) is a
function that describes the effect that collisions between molecules have on the velocity
distribution function. In the Lattice Boltzmann method, Ω(f) describes the effect of the
‘collision’ step in the calculation.
4.1.2 Scaling
Two scales are used throughout this work. A characteristic length, time and velocity for
the physical problem L, tf and U (where Utf = L) are defined. This establishes the non-
34
dimensional variables in the fluid dynamic scale (macroscopic), denoted by the use of a
caret, like .ˆ
xˆα =
xα
L
tˆ =
t
tf
uˆα =
uα
U
(4.5)
Another characteristic length, time and velocity ∆x, ∆t and c (where ∆tc = ∆x) are also
defined. This is the computational scale. ∆x and ∆t are the grid length and time step
and c will be called the ‘grid speed’. Variables in this scale are denoted with a˜.
x˜α =
xα
∆x
t˜ =
t
∆t
u˜α =
uα
c
(4.6)
The key advantage of the computational scale is that the scaled grid spacing and time step
are unity. Two small parameters, Ma and Kn, are defined as,
Kn =
∆x
L
Ma =
U
c
(4.7)
The same scale is used for the density in the computational scale and in the fluid dynamic
scale, such that ρ˜ = ρˆ = ρ/ρ0. ρ0 is a reference density, the value of which should be chosen
so that ρ˜ is close to unity, to reduce numerical errors. The velocity distribution function f
is scaled in the same way, f˜ = fˆ = f/ρ0 .
For the pressure, sufficient scales are already defined thus,
p˜ =
p
ρ0c2
pˆ =
p
ρ0U2
(4.8)
4.1.3 The truncated nine moment system
The continuum velocity distribution function can be described by an infinite set of integrals,
or moments. However once the velocity set is reduced to nine velocities, just nine moments
are needed to exactly reconstruct the discrete velocity distribution function. The definition
of these moments is only subject to the condition that they must be independent, and
different sets are commonly used. This discussion is restricted to the kinetic style set, and
35
Π0
ΠX ΠY
ΠXX ΠXY ΠY Y
ΠXXX ΠXXY ΠXY Y ΠY Y Y
ΠXXXX ΠXXXY ΠXXY Y ΠXY Y Y ΠY Y Y Y
Table 4.1: Moments of the D2Q9 velocity set
the results of other authors are converted into this set when appropriate.
The general notation for each moment, Π, is,
ΠXXX︸ ︷︷ ︸,Y Y Y︸ ︷︷ ︸ = ∑i emx,ieny,ifi
m n
(4.9)
A special case is Π0 =
∑
i fi. The first nine moments of this infinite set can be arranged
as shown in Table 4.1. Due to the fact that the magnitude of the velocity components are
unity in the computational scale, e˜3α = e˜α. This means that the moments that are ‘grayed
out’ in Table 4.1 are not independent, and,
ΠXXX =
∑
i
e3x,ifi =
∑
i
c3e˜3x,ifi =
∑
i
c3e˜x,ifi =
∑
i
c2ex,ifi = c
2ΠX (4.10)
Similarly,
ΠY Y Y = c
2ΠY , ΠXXXY = c
2ΠXY , ΠXY Y Y = c
2ΠXY (4.11)
The flux of each moment is related to the moments on the line beneath it in Table 4.1,
as will be shown later. The nine independent moments may be arranged in a vector of
moments m as,
m = (Π0 ΠX ΠY ΠXX ΠY Y ΠXY ΠXY Y ΠXXY ΠXXY Y )
T (4.12)
36
In the computational scale the relation between m and f may be expressed as,
m˜ = M˜f˜ (4.13)
where M˜ is the transform matrix,
M˜ =

1 1 1 1 1 1 1 1 1
0 1 0 −1 0 1 −1 −1 1
0 0 1 0 −1 1 1 −1 −1
0 1 0 1 0 1 1 1 1
0 0 1 0 1 1 1 1 1
0 0 0 0 0 1 −1 1 −1
0 0 0 0 0 1 −1 −1 1
0 0 0 0 0 1 1 −1 −1
0 0 0 0 0 1 1 1 1

(4.14)
4.1.4 The scaled evolution equation
The dimensions of the collision integral, Ω(f), is density time−1, and it is non-dimensionalised
by Ω˜(f) = Ω(f)∆t/ρ0. Picking a mix of scalings which will later prove useful one can
rewrite the evolution equation in non-dimensional form,
1
tf
∂f˜i
∂tˆ
+
c
L
e˜x,i
∂f˜i
∂xˆ
+
c
L
e˜y,i
∂f˜i
∂yˆ
=
Ω˜(f˜i)
tc
(4.15)
Then, substituting the definitions of the small parameters Ma and Kn yields,
MaKn
∂f˜i
∂tˆ
+ Kne˜x,i
∂f˜i
∂xˆ
+ Kne˜y,i
∂f˜i
∂yˆ
= Ω˜(f˜i) (4.16)
This equation will provide the vital link between the computational scheme and the equa-
tions for the macroscopic variables, as it links the quantities computed in the Boltzmann
scale to the fluid dynamic length and time scales. In the D2Q9 model there are 9 values for
37
i, and nine equations, which translate into nine independent equations for the macroscopic
variables.
4.1.5 Moments of the evolution equation
Nine unique moments were defined in Section 4.1.3 and now an equation for the evolution of
these moments will be derived. This may easily be achieved using the evolution equation for
the distribution function. Taking the vector form of Equation (4.16) and pre-multiplying
by M˜ , results in the evolution equation for the moments,
MaKn
∂
∂tˆ

Π˜0
Π˜X
Π˜Y
Π˜XX
Π˜Y Y
Π˜XY
Π˜XY Y
Π˜XXY
Π˜XXY Y

+ Kn
∂
∂xˆ

Π˜X
Π˜XX
Π˜XY
Π˜X
Π˜XY Y
Π˜XXY
Π˜XXY Y
Π˜XY
Π˜XY Y

+ Kn
∂
∂yˆ

Π˜Y
Π˜XY
Π˜Y Y
Π˜XXY
Π˜Y
Π˜XY Y
Π˜XY
Π˜XXY Y
Π˜XXY

=

∑
i Ω˜i(f)∑
i e˜x,iΩ˜i(f)∑
i e˜y,iΩ˜i(f)∑
i e˜x,ie˜x,iΩ˜i(f)∑
i e˜y,ie˜y,iΩ˜i(f)∑
i e˜x,ie˜y,iΩ˜i(f)∑
i e˜x,ie˜y,ie˜y,iΩ˜i(f)∑
i e˜x,ie˜x,ie˜y,iΩ˜i(f)∑
i e˜x,ie˜x,ie˜y,ie˜y,iΩ˜i(f)

(4.17)
38
This provides further proof that the moment system is truncated, as the evolution of each
moment can be described without reference to any moment outside the chosen set of nine.
Thus there is no direct equivalent to Newton’s viscosity law; the momentum equation
cannot be expressed purely in terms of density and velocity. Instead, the viscosity is
controlled by the relaxation frequency of the second order moments.
4.1.6 The equation of state
The Ideal Gas equation links the pressure of the gas to the density,
p˜ = R˜T ρ˜ (4.18)
Where R˜T is the gas constant divided by c2. Several authors, For example Reis and Phillips
(2008) and Lallemand and Luo (2000), propose ˜(RT ) = 1/3 for maximum stability and
accuracy. However in doing so they are abandoning the physical equation of state in
favour of a new, unphysical, link between pressure and density that is a function of the
scaling used. In the incompressible limit the equation of state may be ignored, but for
multicomponent applications it is vital that the equation of state is correct. The correct
Ideal gas equation is used in this chapter to ensure the results may be extended to multi-
component gas mixtures.
4.2 The equilibrium distribution
Conventionally, equilibrium has been defined when all the moments of Ω(f) are zero. For
the continuous Boltzmann equation this is a sensible proposition; it is not the same as
saying that there are no collisions, just that the collisions on average have no effect on the
macroscopic properties. In the Lattice Boltzmann method equilibrium is defined in the
same way, noting the higher moments will still be affected by the finite number of velocities,
and so the particle distribution function cannot be compared to its continuum equivalent.
Therefore, in this section the right hand side of equation (4.17) will be set to zero and the
39
moments given the superscript (eq) to indicate that they refer to the equilibrium state.
4.2.1 Mass conservation
The first moment of the evolution equation, (4.17), at equilibrium is
KnMa
∂Π˜
(eq)
0
∂tˆ
+ Kn
∂Π˜
(eq)
X
∂xˆ
+ Kn
∂Π˜
(eq)
Y
∂yˆ
= 0 (4.19)
The mass conservation equation in the fluid dynamic scaling is,
∂ρˆ
∂tˆ
+
∂(ρˆuˆx)
∂xˆ
+
∂(ρˆuˆy)
∂yˆ
= 0 (4.20)
By comparison between the form of these equations, the lower order moments of the equi-
librium function may be defined as,
Π˜
(eq)
0 = ρˆ (4.21)
Π˜
(eq)
X = Maρˆuˆx (4.22)
Π˜
(eq)
Y = Maρˆuˆy (4.23)
Using these relations the D2Q9 model exactly recovers the mass conservation law with a
discrete velocity set. Numerical effects due to discretizing space and time have yet to be
considered.
4.2.2 Momentum conservation
The second and third moments of the evolution equation, (4.17), at equilibrium are:
Kn
2
Ma
∂(ρˆuˆx)
∂tˆ
+ Kn
∂Π˜eqXX
∂xˆ
+ Kn
∂Π˜eqXY
∂yˆ
= 0 (4.24)
Kn
2
Ma
∂(ρˆuˆy)
∂tˆ
+ Kn
∂Π˜eqXY
∂xˆ
+ Kn
∂Π˜eqY Y
∂yˆ
= 0 (4.25)
40
where Equation (4.22) and (4.23) have been used. The equilibrium form of the momentum
conservation equation is
∂(ρˆuˆx)
∂tˆ
+
∂(pˆ+ ρˆuˆ2x)
∂xˆ
+
∂(ρˆuˆyuˆx)
∂yˆ
= 0 (4.26)
∂(ρˆuˆy)
∂tˆ
+
∂(ρˆuˆyuˆx)
∂xˆ
+
∂(pˆ+ ρˆuˆ2y)
∂yˆ
= 0 (4.27)
From comparison between these two sets of equations the second order moments must be,
Π˜eqXX = 
2
Ma(pˆ+ ρˆuˆ
2
x) (4.28)
Π˜eqY Y = 
2
Ma(pˆ+ ρˆuˆ
2
y) (4.29)
Π˜eqXY = 
2
Ma(ρˆuˆyuˆx) (4.30)
4.2.3 Higher moments
The higher moments present a problem for the nine velocity model. In continuous kinetic
theory the first order moments are the mass flux, the second order moments are the mo-
mentum flux and the third order moments are the energy flux. Therefore, the appropriate
equation to use for defining the equilibrium values of the third order moments would be the
conservation of energy. However, it is in the third order moments that the D2Q9 lattice
starts to show its deficiencies (Π˜XXX = Π˜X and Π˜Y Y Y = Π˜Y ), and so it is impossible to re-
cover the energy equation correctly. The scheme is often referred to as isothermal because
the BGK collision operator is intentional changed to conserve temperature rather than
energy under collisions, see Dellar (2001). D’Humie`res (1992) formulated a thermal D2Q9
scheme, but it only correctly reproduces the Euler equations. The non-equilibrium part of
the third order moments should not be regarded as meaningful, although the equilibrium
parts have to take the correct values that yield the Navier-Stokes equations.
There is another problem with the higher moments: they may cause instability. Insta-
bility of this type is normally associated with an oscillation on a wavelength equal to the
grid spacing, and cannot be identified using the long wavelength analysis used here. The
41
priority when assigning equilibrium values to the higher moments is therefore to suppress
this instability. Dellar (2002b) has analysed these unstable modes and recommends values
for the higher moments. Because he has used a different moment base, they are translated
here,
Π˜eqXXY = Ma
1
3
ρˆuˆy (4.31)
Π˜eqXY Y = Ma
1
3
ρˆuˆx (4.32)
Π˜eqXXY Y = 
2
Ma
1
3
(p+ ρˆuˆ2x + ρˆuˆ
2
y) (4.33)
The third order equilibrium moments also appear in the expression for the momentum
flux density tensor, however the fourth order moment has no impact on the hydrodynamic
variables. As such, any interest in it is small; if it is stable that is sufficient.
4.2.4 Recovering the equilibrium distribution function
Once the moments of the equilibrium distribution function are known, it is a trivial task
to recover the actual function, using the linear mapping matrix M˜ .
f˜ eq = M˜−1m˜eq (4.34)
42
Where the equilibrium moments are,
m˜eq =

Π˜
(eq)
0
Π˜
(eq)
X
Π˜
(eq)
Y
Π˜
(eq)
XX
Π˜
(eq)
Y Y
Π˜
(eq)
XY
Π˜
(eq)
XY Y
Π˜
(eq)
XXY
Π˜
(eq)
XXY Y

=

ρˆ
Maρˆuˆx
Maρˆuˆy
2Ma(pˆ+ ρˆuˆ
2
x)
2Ma(pˆ+ ρˆuˆ
2
y)
2Ma(ρˆuˆxuˆy)
Ma
1
3
ρˆuˆx
Ma
1
3
ρˆuˆy
2Ma
1
3
(pˆ+ ρˆuˆ2x + ρˆuˆ
2
y)

=

ρ˜
ρ˜u˜x
ρ˜u˜y
(p˜+ ρ˜u˜2x)
(p˜+ ρ˜u˜2y)
(ρ˜u˜xu˜y)
1
3
ρ˜u˜x
1
3
ρ˜u˜y
1
3
(p˜+ ρ˜u˜2x + ρ˜u˜
2
y)

(4.35)
Using the equation of state to remove references to pressure, f˜ (eq) can be written as,
f˜ eq =

f˜
(eq)
0
f˜
(eq)
1
f˜
(eq)
2
f˜
(eq)
3
f˜
(eq)
4
f˜
(eq)
5
f˜
(eq)
6
f˜
(eq)
7
f˜
(eq)
8

=

−1
3
ρ˜(5 ˜(RT )− 3 + 2(u˜2x + u˜2y))
1
6
ρ˜(2u˜2x − u˜2y + 2u˜x + 2 ˜(RT ))
1
6
ρ˜(−u˜2x + 2u˜2y + 2u˜y + 2 ˜(RT ))
1
6
ρ˜(2u˜2x − u˜2y − 2u˜x + 2 ˜(RT ))
1
6
ρ˜(−u˜2x + 2u˜2y − 2u˜y + 2 ˜(RT ))
1
12
ρ˜(u˜2x + u˜
2
y + 3u˜xu˜y + u˜x + u˜y +
˜(RT ))
1
12
ρ˜(u˜2x + u˜
2
y − 3u˜xu˜y − u˜x + u˜y + ˜(RT ))
1
12
ρ˜(u˜2x + u˜
2
y + 3u˜xu˜y − u˜x − u˜y + ˜(RT ))
1
12
ρ˜(u˜2x + u˜
2
y − 3u˜xu˜y + u˜x − u˜y + ˜(RT ))

(4.36)
4.3 The MRT collision operator
The Multiple Relaxation Time (MRT) collision operator involves a linear relaxation to
equilibrium. Each moment of the distribution function will relax towards its equilibrium
value at a different rate, controlled by a relaxation frequency. These relaxation frequencies
43
are assembled in a matrix Λ,
Λ =

λ0 0 0 0 0 0 0 0 0
0 λ1 0 0 0 0 0 0 0
0 0 λ2 0 0 0 0 0 0
0 0 0 λ
ξ+λν
2
λξ−λν
2
0 0 0 0
0 0 0 λ
ξ−λν
2
λξ+λν
2
0 0 0 0
0 0 0 0 0 λν 0 0 0
0 0 0 0 0 0 λ6 0 0
0 0 0 0 0 0 0 λ7 0
0 0 0 0 0 0 0 0 λ8

(4.37)
λν controls shear viscosity while λξ is related to bulk viscosity. Splitting the terms like this
will enable these two effects to be separated. Isotropy requires that λ1 = λ2 and λ6 = λ7.
The collision operator for the velocity distribution function is constructed by pre-
multiplying the distribution function evolution equation, (4.16), by M to get the moment
evolution equation, then setting the right hand side to be a linear relaxation to equilibrium.
MΩ(f) = Λ(meq −m) (4.38)
Premultiplying by the inverse of M , the collision operator becomes,
Ω(f) = M−1Λ(meq −m) = M−1ΛM(f eq − f) (4.39)
In the applications considered in this work mass is always conserved during a collision, and
hence the value for λ0 will always be multiplied by zero and is irrelevant. If mass is not
conserved, for example due to chemical reactions, this would not be the case. Likewise,
for single component gases, momentum is conserved and λ1 and λ2 are irrelevant. In
multi-component gases momentum is exchanged between the species.
44
4.3.1 Evolution of moments
Taking moments of the evolution Equation (4.16) results in the evolution equation for the
nine moments, and if the MRT collision operator is used the result is,
MaKn
∂
∂tˆ

Π˜0
Π˜X
Π˜Y
Π˜XX
Π˜Y Y
Π˜XY
Π˜XY Y
Π˜XXY
Π˜XXY Y

+ Kn
∂
∂xˆ

Π˜X
Π˜XX
Π˜XY
Π˜X
Π˜XY Y
Π˜XXY
Π˜XXY Y
Π˜XY
Π˜XY Y

+ Kn
∂
∂yˆ

Π˜Y
Π˜XY
Π˜Y Y
Π˜XXY
Π˜Y
Π˜XY Y
Π˜XY
Π˜XXY Y
Π˜XXY

=

λ˜0(Π˜eq0 − Π˜0)
λ˜1(Π˜eqX − Π˜X)
λ˜2(Π˜eqY − Π˜Y )
λ˜ξ+λ˜ν
2
(Π˜eqXX − Π˜XX) + λ˜
ξ−λ˜ν
2
(Π˜eqY Y − Π˜Y Y )
λ˜ξ+λ˜ν
2
(Π˜eqY Y − Π˜Y Y ) + λ˜
ξ−λ˜ν
2
(Π˜eqXX − Π˜XX)
λ˜ν(Π˜eqXY − Π˜XY )
λ˜6(Π˜eqXY Y − Π˜XY Y )
λ˜7(Π˜eqXXY − Π˜XXY )
λ˜8(Π˜eqXXY Y − Π˜XXY Y )

(4.40)
To obtain the correct mass conservation equation Π˜0 = Π˜
eq
0 = ρˆ and Π˜α = Π˜
(eq)
α =
Maρˆuˆα. The momentum equation involves the terms, Π˜XX , Π˜Y Y and Π˜XX , which are not
at equilibrium. It will require some assumptions to simplify these terms.
45
4.3.2 The momentum flux density tensor (ΠXX,ΠXY and ΠY Y )
The components of the momentum flux density tensor appear in the 4th, 5th and 6th rows
of Equation (4.40), where they are expressed in terms of the derivatives of the other mo-
ments. An expression for the higher moments is currently unknown, but could be calculated
from the final three rows of Equation (4.40). Such a calculation would be very complicated
and to simplify this problem the Maxwellian iteration of Ikenberry and Truesdell (1956)
is used. VanKampen (1987) and Dellar (2002a) proved that this approach is equivalent
to the Chapman-Enskog method for recovering the Navier-Stokes equations. Equilibrium
values (from Equation (4.35)) are substituted into the left hand side of Equation (4.40),
while retaining the real distribution function (that is, not at equilibrium) on the right hand
side. The 4th, 5th and 6th rows become:
MaKn
∂
∂tˆ

2Ma(pˆ+ ρˆuˆ
2
x)
2Ma(pˆ+ ρˆuˆ
2
y)
2Ma(ρˆuˆxuˆy)
+ Kn ∂∂xˆ

Maρˆuˆx
Ma
1
3
ρˆuˆx
Ma
1
3
ρˆuˆy
+ Kn ∂∂yˆ

Ma
1
3
ρˆuˆy
Maρˆuˆy
Ma
1
3
ρˆuˆx

≈

λ˜ξ+λ˜ν
2
(Π˜eqXX − Π˜XX) + λ˜
ξ−λ˜ν
2
(Π˜eqY Y − Π˜Y Y )
λ˜ξ+λ˜ν
2
(Π˜eqY Y − Π˜Y Y ) + λ˜
ξ−λ˜ν
2
(Π˜eqXX − Π˜XX)
λ˜ν(Π˜eqXY − Π˜XY )
 (4.41)
The Ideal Gas Equation, (4.18), is used to replace pressure terms with density terms,
and when the expressions for the equilibrium moments is included on the left hand side
46
Equation (4.41) becomes,
MaKn
∂
∂tˆ

R˜T ρˆ+ 2Ma(ρˆuˆ
2
x)
R˜T ρˆ+ 2Ma(ρˆuˆ
2
y)
2Ma(ρˆuˆxuˆy)
+ Kn ∂∂xˆ

Maρˆuˆx
Ma
1
3
ρˆuˆx
Ma
1
3
ρˆuˆy
+ Kn ∂∂yˆ

Ma
1
3
ρˆuˆy
Maρˆuˆy
Ma
1
3
ρˆuˆx

≈

λ˜ξ+λ˜ν
2
(2Ma(pˆ+ ρˆuˆ
2
x)− Π˜XX) + λ˜
ξ−λ˜ν
2
(2Ma(pˆ+ ρˆuˆ
2
y)− Π˜Y Y )
λ˜ξ+λ˜ν
2
(2Ma(pˆ+ ρˆuˆ
2
y)− Π˜Y Y ) + λ˜
ξ−λ˜ν
2
(2Ma(pˆ+ ρˆuˆ
2
x)− Π˜XX)
λ˜ν(2Maρˆuˆxuˆy − Π˜XY )
 (4.42)
With any physical (that is, independent of scaling) equation of state the correct unsteady
equations cannot be recovered, only the steady-state Navier-Stokes equations can be re-
covered. To illustrate this point the terms that lead to this error will be included, but
removed at a later stage. However, to simplify the algebra the terms of the type ∂(ρuu)/∂t
will be set to zero at this stage. The mass conservation equation, (4.20), is used to change
the remaining time derivative terms to spatial derivatives.
KnMa
∂
∂xˆ

(1− R˜T )ρˆuˆx
(1
3
− R˜T )ρˆuˆx
1
3
ρˆuˆy
+ KnMa ∂∂yˆ

(1
3
− R˜T )ρˆuˆy
(1− R˜T )ρˆuˆy
1
3
ρˆuˆx

≈

λ˜ξ+λ˜ν
2
(2Ma(pˆ+ ρˆuˆ
2
x)− Π˜XX) + λ˜
ξ−λ˜ν
2
(2Ma(pˆ+ ρˆuˆ
2
y)− Π˜Y Y )
λ˜ξ+λ˜ν
2
(2Ma(pˆ+ ρˆuˆ
2
y)− Π˜Y Y ) + λ˜
ξ−λ˜ν
2
(2Ma(pˆ+ ρˆuˆ
2
x)− Π˜XX)
λ˜ν(2Maρˆuˆxuˆy − Π˜XY )
 (4.43)
47
After rearranging, the normal components of the momentum flux density tensor are given
by,
 Π˜XX
Π˜Y Y
 ≈
 2Ma(pˆ+ ρˆuˆ2x)
2Ma(pˆ+ ρˆuˆ
2
y)
−
MaKn
[
∂
∂xˆ
 (13( 1λ˜ξ + 1λ˜ν ) + 1λ˜ξ (13 − R˜T ))ρˆuˆx
(1
3
( 1
λ˜ξ
− 1
λ˜ν
) + 1
λ˜ξ
(1
3
− R˜T ))ρˆuˆx

+
∂
∂yˆ
 (13( 1λ˜ξ − 1λ˜ν ) + 1λ˜ξ (13 − R˜T ))ρˆuˆy
(1
3
( 1
λ˜ξ
+ 1
λ˜ν
) + 1
λ˜ξ
(1
3
− R˜T ))ρˆuˆy
] (4.44)
Obviously if R˜T = 1
3
then some terms become zero, but if R˜T > 1
3
these terms dominate.
This is the source of the instability experienced at high values of R˜T , this will be seen
more clearly in Equation (4.48). The momentum flux density tensor is completed by the
third row of equation (4.43),
Π˜XY ≈ 2Maρˆuˆxuˆy − MaKn
1
λ˜ν
1
3
(
∂
∂xˆ
(ρˆuˆy) +
∂
∂yˆ
(ρˆuˆx)
)
(4.45)
4.3.3 The momentum equation
Consider the second and third moment of the evolution equation, (4.40).
MaKn
∂
∂tˆ
 Π˜X
Π˜Y
 + Kn ∂
∂xˆ
 Π˜XX
Π˜XY
 + Kn ∂
∂yˆ
 Π˜XY
Π˜Y Y
 =
 λ˜1(Π˜eqX − Π˜X)
λ˜2(Π˜eqY − Π˜Y )

(4.46)
The right hand side of this equation becomes zero in all cases. Substituting the definition
of the momentum flux density tensor, Equations (4.44) and (4.45) into Equation (4.46)
48
yields,
2MaKn
∂
∂tˆ
 ρˆuˆx
ρˆuˆy
+ Kn ∂
∂xˆ
 2Ma(pˆ+ ρˆuˆ2x)
0
+ Kn ∂
∂yˆ
 0
2Ma(pˆ+ ρˆuˆ
2
y)

+ Kn
∂
∂xˆ
−MaKn( ∂∂xˆ [(13( 1λ˜ξ + 1λ˜ν )+ 1λ˜ξ (13−R˜T ))ρˆuˆx]− ∂∂yˆ [(13( 1λ˜ξ− 1λ˜ν )+ 1λ˜ξ (13−R˜T ))ρˆuˆy])
2Maρˆuˆxuˆy − MaKn 13λ˜ν ( ∂∂xˆ(ρˆuˆy) + ∂∂yˆ (ρˆuˆx))

+ Kn
∂
∂yˆ
 2Maρˆuˆxuˆy − MaKn 13λ˜ν ( ∂∂xˆ(ρˆuˆy) + ∂∂yˆ (ρˆuˆx))
−MaKn
(
∂
∂xˆ
[(1
3
( 1
λ˜ξ
− 1
λ˜ν
)+ 1
λ˜ξ
(1
3
−R˜T ))ρˆuˆx]− ∂∂yˆ [(13( 1λ˜ξ + 1λ˜ν )+ 1λ˜ξ (13−R˜T ))ρˆuˆy]
)
= 0 (4.47)
These equations can be rearranged into the standard form for the Navier-Stokes equations,
2MaKn
∂
∂tˆ
 ρˆuˆx
ρˆuˆy
+ 2MaKn ∂∂xˆ
 ρˆuˆ2x
ρˆuˆxuˆy
+ 2MaKn ∂∂yˆ
 ρˆuˆxuˆy
ρˆuˆ2y
 =
− 2MaKn
 ∂∂xˆ(pˆ)
∂
∂yˆ
(pˆ)
+ Ma2Kn ∂∂xˆ
 ∂∂xˆ( 1λ˜ξ + 1λ˜ν )13 ρˆuˆx + ∂∂yˆ ( 1λ˜ξ − 1λ˜ν )13 ρˆuˆy
1
3λ˜ν
∂
∂xˆ
(ρˆuˆy) +
1
3λ˜ν
∂
∂yˆ
(ρˆuˆx)

+ Ma
2
Kn
∂
∂yˆ
 13λ˜ν ∂∂xˆ(ρˆuˆy) + 13λ˜ν ∂∂yˆ (ρˆuˆx)
∂
∂xˆ
( 1
λ˜ξ
− 1
λ˜ν
)1
3
ρˆuˆx +
∂
∂yˆ
( 1
λ˜ξ
+ 1
λ˜ν
)1
3
ρˆuˆy

+ Ma
2
Kn
∂
∂xˆ
 1λ˜ξ (13 − R˜T )(∂ρˆuˆx∂xˆ + ∂ρˆuˆy∂yˆ )
0
+ Ma2Kn ∂∂yˆ
 0
1
λ˜ξ
(1
3
− R˜T )(∂ρˆuˆx
∂xˆ
+ ∂ρˆuˆy
∂yˆ
)

(4.48)
The final terms are zero in the steady state, as can been seen when the mass conservation
equation (4.20) is considered, and when R˜T = 1
3
. If R˜T > 1
3
then these terms represent
a negative bulk viscosity, and it is no surprise that this results in unstable behaviour. At
49
this stage these terms are removed to reveal,
2MaKn
 ∂
∂xˆ
 ρˆuˆ2x
ρˆuˆxuˆy
+ ∂
∂yˆ
 ρˆuˆxuˆy
ρˆuˆ2y
 = −2MaKn
 ∂∂xˆ(pˆ)
∂
∂yˆ
(pˆ)

+ Ma
2
Kn
 ∂2
∂xˆ2
 ( 1λ˜ξ + 1λ˜ν )13 ρˆuˆx
1
λ˜ν
1
3
ρˆuˆy
+ ∂2
∂yˆ2
 1λ˜ν 13 ρˆuˆx
( 1
λ˜ξ
+ 1
λ˜ν
)1
3
ρˆuˆy
+ ∂2
∂xˆ∂yˆ
 1λ˜ξ 13 ρˆuˆy
1
λ˜ξ
1
3
ρˆuˆx

(4.49)
This equation has the same form as the compressible Navier-Stokes equation. The shear
viscosity terms can be controlled by a suitable choice for λν , while the bulk viscosity is
controlled by λξ.
The kinematic shear viscosity, ν, of the fluid is given by,
λ˜ν =
Kn
Ma
1
3νˆ
=
1
3ν˜
(4.50)
The bulk viscosity coefficient, ξ, with the same units as the kinematic viscosity, is given
by,
λ˜ξ =
Kn
Ma
1
3ξˆ
=
1
3ξ˜
(4.51)
There are several interpretations of bulk viscosity, see Dellar (2001) for details, and
this definition is not regarded as standard. The relationship between Dellar’s µ′ (which
is the same as µB used by Vincenti and Kruger (1976)) and ξ is ξρ = µ/3 + µ
′. For a
perfect monatomic gas µ′ = 0 and so ξ = ν/3. It is interesting to note that the standard
single relaxation time model is obtained by setting all the relaxation rates equal, including
λν = λξ meaning ξ = ν, and the Navier-Stokes equations recovered have bulk viscosity
terms.
50
Introducing the definition of the viscosities into Equation (4.49) yields,
∂
∂xˆ
 ρˆuˆ2x
ρˆuˆxuˆy
+ ∂
∂yˆ
 ρˆuˆxuˆy
ρˆuˆ2y
 = −
 ∂∂xˆ(pˆ)
∂
∂yˆ
(pˆ)
+ νˆ
 ∂2
∂xˆ2
 ρˆuˆx
ρˆuˆy
+ ∂2
∂yˆ2
 ρˆuˆx
ρˆuˆy

+ ξˆ
 ∂2
∂xˆ2
 ρˆuˆx
0
+ ∂2
∂yˆ2
 0
ρˆuˆy
+ ∂2
∂xˆ∂yˆ
 ρˆuˆy
ρˆuˆx
 (4.52)
4.4 Conclusions
There are 9 unique moments corresponding to the nine discrete velocities of the D2Q9
scheme. Three of the moments are the fluid density and velocity components. Another
three of the moments are the components of the momentum flux density tensor. The
third order moments, ΠXXY and ΠXY Y , must have specific equilibrium values that yield
the Navier-Stokes equations, but the non-equilibrium parts of these moments is meaning-
less. For the ΠXXY Y moment, the equilibrium value is selected for stability rather than
hydrodynamic reasons, and the non-equilibrium part is meaningless.
A model which recovers the steady-state isothermal compressible Navier-Stokes equa-
tions was described. The use of the ideal gas equation of state means that a time accurate
solution cannot be found. A restriction on the magnitude of R˜T , and thus on c, was
identified from the expression for the momentum flux density tensor. The shear and bulk
viscosities may be independently controlled by setting values for the relaxation frequencies
λν and λξ respectively.
51
Chapter 5
Boundary Conditions and Numerical
Implementation
In Chapter 4 a LB scheme was developed and analysed to show that it is equivalent to
the steady-state isothermal compressible Navier-Stokes equations. Although this scheme
was based on a reduced set of particle velocities, it is still a continuous scheme. To make
it suitable for computation, the derivatives must be replaced by discrete approximations.
In this chapter a commonly used method for an explicit second order accurate scheme
will be described. Before that, a new method of analysis for boundary conditions using
the moment basis, rather than the particle basis, will be introduced. This method allows
analysis of commonly used boundary conditions, and will be used when new boundary
conditions are developed in subsequent chapters.
All variables in Section 5.1 are in the computational scale, and for clarity the˜is dropped.
5.1 Boundary conditions for LB
5.1.1 Kinetic style boundaries
Reading the literature, two clear themes come to the fore. The first is the use of kinetic style
boundary conditions, such as specular, diffuse or bounce-back reflection. The motivation
52
for this approach is clear as these are the boundary conditions used when solving the
continuous Boltzmann equation. The second theme is the use of ‘half-grid-wall’ schemes,
where the location of the wall is assumed to be half way between the last and second from
last row of nodes. The motivation for this approach is that the kinetic style boundaries
do not work satisfactorily. Of course, in taking this approach the initial assumption that
boundary conditions for the continuous Boltzmann equation can be used for LB methods
is violated, but these schemes do work in some cases, as shown in Section 5.1.4.
Given the analysis of the previous chapters, it is clear that kinetic style boundaries
are not appropriate for the D2Q9 system. When dealing with the 9 moment truncated
system, rather than the full set of infinitely many moments, a key link with the continuous
Boltzmann equation has been lost. Principally, the length of the grid spacing is much larger
than the mean free path, and physical effects at a scale smaller than the hydrodynamic
scale are not recreated. Therefore, the conditions at the boundary should be viewed as
macroscopic, not microscopic.
5.1.2 Moment grouping
This section approaches the subject of boundary conditions from a fresh perspective, with-
out any pretensions to replicate the full Boltzmann equation. Instead, a relatively simple
problem is examined; at a straight boundary (corners are dealt with in Section 5.1.5) values
for the three unknown elements of f must be supplied. Therefore 3 independent equations
must be identified. As the problem is expressed in terms of the hydrodynamic moments,
the boundary conditions should be specified in these terms. To do this, the unknown f
values must be transformed into unknown moments. For clarity two general cases, a hor-
izontal wall and a vertical open boundary, will be considered. Of course any orientation
can be analysed by application of the principles of symmetry. The horizontal wall shown
in Figure 5.1a has unknown f values f2, f5 and f6. By setting the other values of f to zero
53
temporarily, it can be seen how these values appear in the moment basis.
M

0
0
f2
0
0
f5
0
0
f8

=

f2 + f5 + f6
f5 − f6
f2 + f5 + f6
f5 + f6
f2 + f5 + f6
f5 − f6
f5 − f6
f5 + f6
f5 + f6

where m =

Π0
ΠX
ΠY
ΠXX
ΠY Y
ΠXY
ΠXY Y
ΠXXY
ΠXXY Y

(5.1)
The unknown f values only occur in certain combinations, these combinations and their
associated moments are listed in Table 5.1. Each group contains moments with the same
power of x in their description, 1. For a vertical open boundary shown in Figure 5.1b with
unknown f values f1, f5 and f8 the process can be repeated and Table 5.2 is obtained.
Solid
Fluid
xy
FluidFluidoutside
domain
xy
a b
Figure 5.1: a: A node at a horizontal wall. b: A node at a vertical open boundary. The unknown
f values are shown with dotted lines. f values running along the boundary are streamed from
neighbouring cells.
1Interestingly, this agrees with the theory of Torrilhon and Struchtrup (2008). They found that from
their truncated 13 moment system the boundary conditions needed at the wall are all the moments with
an ‘odd’ power of y (ΠY ,ΠXY ,ΠXXY ).
54
Unknowns Moments
f2 + f5 + f6 Π0,ΠY ,ΠY Y
f5 − f6 ΠX ,ΠXY ,ΠXY Y
f5 + f6 ΠXX ,ΠXXY ,ΠXXY Y
Table 5.1: Moment groups for a horizontal
boundary
Unknowns Moments
f1 + f5 + f8 Π0,ΠX ,ΠXX
f5 − f8 ΠY ,ΠXY ,ΠXXY
f5 + f8 ΠY Y ,ΠXY Y ,ΠXXY Y
Table 5.2: Moment groups for a vertical bound-
ary
For each boundary, one macroscopic variable from each group must be specified to
uniquely specify the boundary condition. Specifying more than one from each group is not
possible.
5.1.3 Analysis of Zou and He boundary conditions
Zou and He (1997) used a very unusual concept, ‘non-equilibrium bounce-back’, to derive
a set of boundary conditions. Although the logic in the derivation is not made clear in
their paper, the conditions themselves have proved to be useful. The method introduced
here allows these conditions to be analysed in a more straightforward way.
The Zou and He condition for a wall boundary, when applied to a node on a wall, is
equivalent to setting,
f2 = f4 (5.2a)
f5 = f7 − 1
2
− (f1 − f3) + 1
2
ρux (5.2b)
f6 = f8 +
1
2
+ (f1 − f3)− 1
2
ρux (5.2c)
55
In terms of moments, the Zou and He boundary becomes,
m =

f0 + f1 + 2f4 + f3 + 2f7 + 2f8
ρux
0
f1 + f3 + 2f7 + 2f8
2f4 + 2f7 + 2f8
ρux + 2f7 − f1 + f3 − 2f8
ρux − f1 + f3
0
2f7 + 2f8

(5.3)
The only moments that do not contain references to f values are ΠX ,ΠY and ΠXXY . In
effect, the conditions that are being applied at the wall boundary are,
ΠX = ρux, ΠY = 0, ΠXXY =
1
3
ρuy = 0 (5.4)
The first two conditions contain hydrodynamic quantities and so are useful in defining
the problem. The third condition, however, seems arbitrary. In the preceding derivation,
setting a value for the quantity ρux was discussed, but for practical problems one may wish
to set just the velocity. However the density, ρ = Π0, is in the same group as ΠY . This
means that the value for ρ is already fixed as,
ρ = f0 + f1 + f3 + 2(f7 + f4 + f8) (5.5)
This value may be used in conjunction with the desired ux to give the desired value of ΠX .
56
The Zou and He condition for an open boundary is,
f1 = f3 +
2
3
ρux (5.6a)
f5 = f7 − 1
2
(f2 − f4) + 1
6
ρux (5.6b)
f8 = f6 +
1
2
(f2 − f4) + 1
6
ρux (5.6c)
Which is equivalent in moment terms to,
m =

f0 + 2f3 + ρux + f2 + f4 + 2f7 + 2f6
ρux
0
2f3 + ρux + 2f7 + 2f6
f2 + f4 + 2f7 +
1
3
ρux + 2f6
2f7 − f2 + f4 − 2f6
1
3
ρux
−f2 + f4
2f7 +
1
3
ρux + 2f6

(5.7)
The only moments that do not contain references to f values are ΠX ,ΠY and ΠXY Y . In
effect, the conditions that are being applied at the open boundary are,
ΠX = ρux, ΠY = 0, ΠXY Y =
1
3
ρux = Π
eq
XY Y (5.8)
where the final relation comes from the definition of the equilibrium function, Equa-
tion (4.35). The first two conditions are desirable, as these are hydrodynamic quantities,
but the final condition cannot easily be justified. The term ρux may be split according to
the rule
ρux = ρ− (f0 + f2 + f4 + 2(f3 + f6 + f7)) (5.9)
Such that either ρ or ρux may be set, but not both.
57
5.1.4 Analysis of ‘bounce-back’, specular and diffusive boundary
conditions
For completeness, the same analysis used in 5.1.3 is applied to the bounce-back boundary
condition, where the unknown f values on a wall are given by,
f2 = f4, f5 = f7, f6 = f8 (5.10)
Leading to the conditions in the moment basis,
ΠY = 0, ΠXY Y = 0, ΠXXY = 0 (5.11)
Two conditions are imposed on the higher non-hydrodynamic moments, while the only
hydrodynamic condition is that the flowrate through the wall is zero. The velocity along
the wall is undefined, and it is this velocity that is often reported as Knudsen slip. In fact
it is not, just a result of specifying the boundary conditions improperly! Nevertheless the
results obtained using bounce-back are close to the correct solution, and it is worthwhile
considering why this is the case. ΠXY Y , like the other moments, is composed of the equi-
librium value plus a non-equilibrium term, ΠXY Y =
1
3
ρux + Π
(neq)
XY Y . As the non-equilibrium
term is small, setting ΠXY Y = 0 is equivalent to setting ux to a small value. One could
derive an exact expression for this value, using a method similar to the one used to obtain
an expression for the non-equilibrium parts of the momentum flux density tensor. This
was done by Ginzburg and Adler (1994), who then set the higher relaxation frequency to
ensure zero velocity. Similarly setting ΠXXY = 0 implies
1
3
ρuy + Π
(neq)
XXY = 0. As uy = 0,
this is effectively specifying Π
(neq)
XXY = 0. Those who defend bounce-back would argue that
it is quick to implement. However, for many people the time saved using bounce-back
was subsequently wasted analysing the slip velocity it creates (see Verhaeghe et al. (2009)
again).
58
The same analysis applied to specular reflection, defined by;
f2 = f4, f5 = f8, f6 = f7 (5.12)
results in the following moment conditions,
ΠY = 0, ΠXY = 0, ΠXXY = 0 (5.13)
Again, these are a fairly arbitrary set of constraints. As there is no longer the constraint
that ΠXY Y = 0, the velocity at the wall is not set to zero, or even a value close to zero.
Instead, the conditions ΠXY = ΠY = 0 lead to the viscous shear stress being set to zero,
which is effectively setting ∂ux/∂y = 0. In channel flow this leads to the x velocity at the
wall being the same as the velocity along the centre of the channel. As with bounce-back,
setting ΠXXY = 0 is equivalent to Π
(neq)
XXY = 0, and it is not clear if this has any impact on
the solution.
A diffusive reflection boundary condition, a LB version of Maxwell’s boundary condition
for the continuum Boltzmann equation, is proposed by many authors, including Verhaeghe
et al. (2009) and Sofonea and Sekerka (2005a). The various schemes for implementing
diffusive reflection all start from arbitrary constraints in the particle basis, and when
viewed in the moment basis simply reveal different arbitrary constraints. The details
are therefore omitted from this work. Verhaeghe et al. (2009) and Guo et al. (2009)
suggest combining the diffusive reflection boundary condition with either the bounce-back
or specular reflection boundary condition. By specifying how much of one or the other
condition is used some control is possible over some of the moments, such as the slip velocity.
However, this seems an unnecessarily complicated and less effective way of achieving what
is described in this work.
59
Fluid
Fluid
outside
domain
Solid
xy
Figure 5.2: A concave corner node. The unknown f values are shown with dotted lines. The
dashed line is the open boundary, the solid line is the wall.
5.1.5 Treatment of corner nodes
At this stage it is worth mentioning the treatment of corners. The literature on LB hardly
touches on the issue, yet simulations show that the treatment is absolutely essential for a
well-behaved solution. The idea of studying boundaries in the moment basis can be used to
suggest an extension of the Zou and He boundary conditions to include corners. However,
this is not satisfactory and a better solution will be proposed later.
Corner nodes can be divided into two types, concave and convex corners. A concave
corner is such that a fluid node is surrounded by five solid nodes, such that there are five
unknown values of f . A convex corner has only one solid node as a neighbour, and as such
is, in most cases, a trivial problem.
As an example consider the concave corner shown in Figure 5.2 with unknown f values
f1, f2, f5, f7 and f8. Substituting zero for the known f values into the moment equation
60
gives:
M

0
f1
f2
0
0
f5
0
f7
f8

=

f1 + f2 + f5 + f7 + f8
f1 + f5 − f7 + f8
f2 + f5 − f7 − f8
f1 + f5 + f7 + f8
f2 + f5 + f7 + f8
f5 + f7 − f8
f5 − f7 + f8
f5 − f7 − f8
f5 + f7 + f8

m =

Π0
ΠX
ΠY
ΠXX
ΠY Y
ΠXY
ΠXY Y
ΠXXY
ΠXXY Y

(5.14)
The moments are not as neatly arranged in groups as they were for straight boundaries.
Obviously the moments are not fully independent, only five may be given values, however
there is some freedom as to which five are selected.
For example, for a corner which forms the intersection between a wall and an open
boundary enforced using the Zou and He method, Equations (5.2) and (5.6), one could
apply the conditions of both, giving,
Π0 = ρopen (5.15)
ΠX = ρopenux (5.16)
ΠY = 0 (5.17)
ΠXY Y =
1
3
ρux = Π
eq
XY Y (5.18)
ΠXXY = 0 (5.19)
61
These conditions in moment space lead to the following f values:
f1 = f3 +
2
3
ρopenux (5.20a)
f2 = f4 (5.20b)
f5 = f7 +
1
6
ρopenux (5.20c)
f6 = −1
2
f0 − f3 − f4 − f7 + 1
2
ρopen − 1
2
ρopenux (5.20d)
f8 = −1
2
f0 − f3 − f4 − f7 + 1
2
ρopen − 1
3
ρopenux (5.20e)
To summarise, this is not the recommended solution to the corner nodes, but a method
for corners that is consistent with the Zou and He open and wall boundary conditions.
To demonstrate the effect of changing the corner nodes, two different boundary condi-
tions were applied to the corner nodes of a simulation of Poiseuille flow in a channel. For
both simulations the boundary at the wall was enforced using bounce back, and the inlet
and outlet using the Zou and He conditions. Firstly, the Zou and He open boundary con-
dition (Equation (5.6)) and then the bounce-back rule (Equation (5.10)) were applied to
the corner nodes, and the result is shown in Figure 5.3. Secondly, the improved treatment
described here (Equation (5.20)) was applied, with zero slip velocity, and the results are
shown in Figure 5.4. It is clear that using the moment basis to apply consistent condi-
tions at the corner greatly improves the accuracy of the solution, the error is one thousand
times smaller in the corner. However, a better solution for Poiseuille flow is proposed in
Chapter 6.
5.2 Numerical implementation
To convert the PDEs described in Chapter 4 into a numerical scheme, a first-order accurate
Forward-Euler scheme is used. This scheme has the advantage of being explicit, but has
low accuracy. Following the approach of many authors, including Chen and Doolen (1998),
a modified collision term is introduced which makes the scheme second-order accurate.
62
y x
u
e
r
r
o
r
5
10
15
20
25
30
5
10
15
20
−0.2
−0.1
0
0.1
0.2
Figure 5.3: uerror is the difference between the velocity calculated and the analytical solution for
Poiseuille flow. The wall follows the lines y = 1 and y = 20. The corner nodes have the Zou and
He conditions for an open boundary, followed by the application of the bounce-back rule. The
corner nodes clearly cause a large error.
y x
u
e
r
r
o
r
5
10
15
20
25
30
5
10
15
20
×10−4
0
1
2
3
4
Figure 5.4: uerror is the difference between the velocity recovered and the analytical solution.
The wall follows the lines y = 1 and y = 20. The corner nodes have the conditions described in
Equation (5.20) imposed. The corner nodes cause a small disturbance to the whole domain
63
Introducing f˜+ to represent f˜(x˜+ 1, y˜ + 1, t˜+ 1), that is after collision and streaming,
the forward Euler scheme can be written as,
f˜+ = f˜ + M˜−1Λ˜′M˜(f˜ (eq) − f˜) (5.21)
Where the prime superscript (′) in Λ˜′ indicates that a correction for numerical errors has
been included. If Λ˜′ = Λ˜ this scheme would only be first-order accurate in space. Taking
moments of Equation (5.21), and picking the XX component of the momentum flux density
tensor evolution equation (as an example), gives,
Π˜+XX − Π˜XX =
λ˜ν′ + λ˜ξ′
2
(Π˜
(eq)
XX − Π˜XX) +
λ˜ν′ − λ˜ξ′
2
(Π˜
(eq)
Y Y − Π˜Y Y ) (5.22)
A similar equation can be derived using a second order accurate central-differencing method,
such as,
(Π˜+XX − Π˜XX) =
1
2
[
λ˜ν + λ˜ξ
2
(
(Π˜
(eq)
XX − Π˜XX) + (Π˜(eq)+XX − Π˜+XX)
)
+
λ˜ν − λ˜ξ
2
(
(Π˜
(eq)
Y Y − Π˜Y Y ) + (Π˜(eq)+Y Y − Π˜+Y Y )
)]
(5.23)
Now, assuming Π
(eq)+
αα = Π
(eq)
αα (where no summation on α is implied), as Π
(eq)+
αα is unknown
at the beginning of the time step, the central differencing method becomes,
(Π˜+XX − Π˜XX)
(
1 +
λ˜ν + λ˜ξ
4
)
=
λ˜ν + λ˜ξ
2
(Π˜
(eq)
XX − Π˜XX)
+
λ˜ν − λ˜ξ
2
(Π˜
(eq)
Y Y − Π˜Y Y )−
λ˜ν − λ˜ξ
4
(Π˜+Y Y − Π˜Y Y ) (5.24)
There is a similar equation which can be derived by repeating the above steps with the
Y Y component of the momentum flux density tensor, and this must be used to remove the
dependence on (Π˜+Y Y − Π˜Y Y ). After some algebra, these two equations may be written in
64
the form of Equation (5.22),
(Π˜+XX − Π˜XX) = 2
λ˜ν + λ˜ξ + λ˜νλ˜ξ
(λ˜ν + 2) + (λ˜ξ + 2)
(Π˜
(eq)
XX − Π˜XX)
− 2λ˜
ν − 2λ˜ξ
(λ˜ν + 2) + (λ˜ξ + 2)
(Π˜
(eq)
Y Y − Π˜Y Y ) (5.25)
For this equation to describe the same scheme as Equation (5.22) the modified collision
times must be,
λ˜ν′ =
λ˜ν
(1 + λ˜
ν
2
)
, λ˜ξ′ =
λ˜ξ
(1 + λ˜
ξ
2
)
(5.26)
65
Chapter 6
Solution of Poiseuille Flow
A LB model for single component flow was developed in Chapter 4, and a method for
analysing and creating boundary conditions was developed in Chapter 5. These will now
be used to model Poiseuille flow. The boundary conditions for Poiseuille flow will be
determined from the analytical solution. Non-dimensional results are included to show the
limits of the model, however the section of greatest interest is the solution to a real physical
problem. This is rarely (maybe never) done in the literature, especially for a gas obeying
the Ideal Gas law.
6.1 Theoretical solution
Before attempting to compute a solution to Poiseuille flow it is informative to recover the
well known analytical solution from the LB model. The Poiseuille problem is flow through
a channel driven by a pressure gradient in the x direction, with zero velocity at the wall
(although a slip velocity is carried through the algebra for generality). The following
conditions are met,
uˆy = 0
∂ρˆ
∂yˆ
= 0
∂pˆ
∂yˆ
= 0
∂
∂tˆ
= 0 (6.1)
66
Introducing these simplifications into the mass conservation equation, (4.20), yields
∂(ρˆuˆx)
∂xˆ
= −∂(ρˆuˆy)
∂yˆ
= 0 (6.2)
The terms are equal to zero because the y velocity is zero everywhere. Given these con-
ditions the momentum flux density tensor, Equation (4.44) and Equation (4.45) becomes:
Π˜XX = 
2
Ma(pˆ+ ρˆuˆ
2
x) (6.3a)
Π˜Y Y = 
2
Mapˆ (6.3b)
Π˜XY = −2Maνˆ(
∂ρˆuˆx
∂y
) (6.3c)
The momentum equation (4.46) becomes,
∂
∂xˆ
 Π˜XX
Π˜XY
 = − ∂
∂yˆ
 Π˜XY
Π˜Y Y
 (6.4)
Into which the values of the momentum flux density tensor given by Equations (6.3) can
be substituted to obtain,
∂
∂xˆ
 (pˆ+ ρˆuˆ2x)
−νˆ(∂ρˆuˆx
∂y
)
 = − ∂
∂yˆ
 −νˆ(∂ρˆuˆx∂y )
pˆ
 (6.5)
With the use of the identities,
∂ρˆuˆ2x
∂xˆ
= 2uˆx
∂ρˆuˆx
∂xˆ
− uˆ2x
∂ρˆ
∂xˆ
= 0− uˆ2x
∂ρˆ
∂xˆ
(6.6)
∂
∂xˆ
∂ρˆuˆx
∂yˆ
=
∂
∂yˆ
∂ρˆuˆx
∂xˆ
= 0 (6.7)
∂
∂yˆ
(pˆ) = 0 (6.8)
67
the second row of (6.5) becomes zero, while the first row becomes,
∂pˆ
∂xˆ
− uˆ2x
∂ρˆ
∂xˆ
=
∂pˆ
∂xˆ
(
1− uˆ
2
x
RˆT
)
= νˆ
(
∂2ρˆuˆx
∂yˆ2
)
(6.9)
where the second term was obtained by using the Ideal Gas equation. As the physical
problems which will be considered are low Mach number, uˆ2x  RˆT and (1− uˆ
2
x
RˆT
) ≈ 1, this
is easily confused with the incompressible limit, and it is worth clarifying what is being
limited. The assumption is that
√
RT is much greater than the mean speed of the flow,
ux. The density is not constant but changes in the density will be 1/RT times changes in
pressure, and therefore they will be very small. However, the relative change in pressure
and density are the same size, ∆p/p = ∆ρ/ρ.
Integrating equation (6.9) with respect to yˆ then setting boundary conditions of zero x
velocity at yˆ = ±Hˆ/2, where H is the channel width, gives,
uˆx =
1
2ρˆνˆ
∂pˆ
∂xˆ
(
yˆ2 − (Hˆ/2)2
)
(6.10)
This is in agreement with the usual theoretical solution of the Poiseuille flow.
6.1.1 Moment based boundary conditions for Poiseuille Flow
Using the analytical solution of the Poiseuille flow, one can see that the higher order
moments, ΠXXY , ΠXY Y and ΠXXY Y , are not needed in order to recover the correct velocity
profile. Therefore, using these moments in boundary conditions, as existing methods such
as bounce-back do, may introduce an error. An approach using only the lower moments,
as in the analytical solution, is now explored. Considering the grouping of moments at a
68
wall, Table 5.1, one possible combination is,
Π˜Y =0 (6.11)
Π˜X =ρ˜u˜slip (6.12)
Π˜XX =p˜+ ρ˜u˜
2
slip (6.13)
For the sake of generality a slip velocity along the wall is included. An expression for the
pressure, pˆ, is required, however. Remembering that the value of ρ˜ can be calculated from
the known f values, pˆ can be calculated using the equation of state, giving
Π˜XX = p˜+ ρ˜u˜
2
slip = ρ˜(R˜T + u˜
2
slip) (6.14)
At an open boundary the moments to be set are,
Π˜0 = ρ˜open (6.15)
Π˜Y = 0 (6.16)
Π˜Y Y = p˜open = R˜T ρ˜open (6.17)
At a corner all these conditions are imposed simultaneously,
Π˜0 = ρ˜open (6.18)
Π˜X = ρ˜openu˜slip (6.19)
Π˜Y = 0 (6.20)
Π˜XX = ρ˜open(R˜T + u˜
2
slip) (6.21)
Π˜Y Y = R˜T ρ˜open (6.22)
A summary of these conditions is shown in Figure 6.1. Using them the unknown f values
can be recovered. At a wall they are,
69
Figure 6.1: Boundary conditions for Poiseuille flow
f˜2 = f˜1 + f˜3 + f˜4 + 2f˜7 + 2f˜8 − (R˜T + u˜2x)ρ˜ (6.23a)
f˜5 = −f˜1 − f˜8 + u˜xρ˜
2
+
1
2
(R˜T + u˜2x)ρ˜ (6.23b)
f˜6 = −f˜3 − f˜7 − u˜xρ˜
2
+
1
2
(R˜T + u˜2x)ρ˜ (6.23c)
at an open boundary,
f˜1 = −f˜0 − f˜3 + ρ˜− R˜T ρ˜ (6.24a)
f˜5 = −f˜2 − f˜6 + 1
2
R˜T ρ˜ (6.24b)
f˜8 = −f˜4 − f˜7 + 1
2
R˜T ρ˜ (6.24c)
and at the intersection of a wall and open boundary,
f˜1 = −f˜0 − f˜3 + ρ˜− R˜T ρ˜ (6.25a)
f˜2 = −f˜0 − f˜4 + ρ˜− (R˜T + u˜2x)ρ˜ (6.25b)
f˜5 =
1
2
(−f˜1 − f˜2 + f˜3 + f˜4) + f˜7 + 1
2
ρ˜u˜x (6.25c)
f˜6 =
1
2
(−f˜0 − f˜2 − f˜4)− f˜3 − f˜7 + 1
2
ρ˜− u˜xρ˜ (6.25d)
f˜8 = −f˜0 − f˜1 − f˜2 − f˜3 − f˜4 − f˜5 − f˜6 − f˜7 + ρ˜ (6.25e)
70
  
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
ΠXXY Y
ΠXY Y /10
−4
ΠXXY /10
−7
ΠY YΠXY /10
−4
ΠXX
ΠY /10
−12
ΠX/10
−3
Π0
0.111
0.1111
0.1112
0
2
4
6
−5
0
5
0.333
0.3335
−5
0
5
0.333
0.3335
−1
0
1
0
1
2
0.999
1
1.001
Figure 6.2: Moments for a solution to the Poiseuille flow arranged according to their order in
Table 4.1. A 20×20 grid was used, the pressure drop was 1×10−3, R˜T = 1/3 and ν˜ = 1. The
moments near the bottom are higher order. Non-dimensional units.
6.2 Non-dimensional results for Poiseuille flow
In this section results from a simple non-dimensional Poiseuille flow are presented. Poiseuille
flow is particularly simple, and the stability and accuracy demonstrated here will be far
higher than could be expected in a more complicated geometry. Unless otherwise stated,
all results in this section are of a solution that has converged, defined by the change in total
momentum throughout the domain from one time step to the next being zero to machine
accuracy. The geometry is a 20×20 grid, with the pressure imposed on both ends using the
method from Section 6.1.1. The accuracy will be measured as the percentage error in the
71
  
 
 
 
 
 
 
 
 
 
 
Π
(neq)
XXY Y /10
−7
Π
(neq)
XY Y /10
−4
Π
(neq)
XXY /10
−7
Π
(neq)
Y Y /10
−12
Π
(neq)
XY /10
−4
Π
(neq)
XX /10
−7
−4
−2
0
−1.2281
−1.2281
−1.2281
−5
0
5
−1
0
1
−5
0
5
−4
−2
0
Figure 6.3: Non-equilibrium part of moments, or the difference between the equilibrium values
and the actual moments, shown in Figure 6.2. The non-equilibrium part of the conserved moments
is zero, and not shown. Non-dimensional units.
cross-channel average velocity half way along the channel compared with the theoretical
solution.
Figure 6.2 shows contours for all the moments, for a typical computation of Poiseuille
flow. It is informative to study carefully the shape of these contours, and compare them
with the equilibrium values; the difference between the two is shown in Figure 6.3. The
density and x momentum behave as the theoretical solution predicts. The y velocity is
interesting, it is small, but there is obviously a pattern to it. This is because of the
compressible nature of the solution. At lower values of x velocity, the y velocity becomes
too small to resolve and this pattern disappears.
The Y Y component of the momentum flux density tensor is very close to its equilibrium
values, however the equilibrium values of the XY momentum flux is effectively zero, and
so the non-equilibrium value is very substantial. As is to be expected, the XX viscous
stress is smaller that the XY viscous stress, and is constant along lines of constant y.
72
  
Solution after 5000 iterations
Converged solution
P
er
ce
n
ta
g
e
er
ro
r
ǫMa
10−15 10−10 10−5 100
10−10
10−5
100
105
Figure 6.4: Percentage error of average velocity mid-way along the channel compared with the-
oretical solution for Poiseuille flow. The pressure difference driving the flow was varied to give
different values of Ma
The first test conducted was to vary the overall pressure difference in order to change the
maximum velocity. The value of R˜T was 1/3, the average scaled pressure was 1/3 and all
relaxation frequencies were set to 1.5. Results are shown in Figure 6.4 for solutions which
were fully converged, and for solutions after 5,000 iterations. Looking at the converged
solutions, clearly the magnitude of Ma has an impact on the accuracy, if it is too high there
is an error due to compressibility effects, and if it is too low the limited numerical range
of the computer introduces an error. There is a range where the accuracy is acceptable.
Within this range the number of iterations to convergence tends to reduce with lower values
of Ma. After 5,000 iterations the main source of error is due to the time dependent nature
of the solution, and this error is roughly independent of Ma.
The effect of changing the value of R˜T was also investigated and the results are shown
in Figure 6.5. Again the average scaled pressure was 1/3 and all relaxation frequencies
were set to 1.5. The pressure gradient was varied to find the maximum stable value of
Ma. The definition of ‘stable’ is not exact, as the error is high with large Ma values, and
therefore these results are only approximate. Values of R˜T above 1/2 are not stable. The
73
R˜T
ǫ
M
a
10
−2
10
−1
10
0
0.1
0.2
0.3
0.4
0.5
Figure 6.5: Maximum stable value of Ma at different values of ˜(RT ). The dashed vertical line
shows ˜(RT ) = 1/3
results in Figure 6.5 are clear, there is maximum stability when R˜T = 1/3. This agrees
with the conclusions reached in Section 4.3.2. For constant values of Ma the accuracy
increases with small values of R˜T , as long as it remains stable, as the change in density is
smaller and thus the compressibility effect is less.
Figure 6.6 shows the result of changing all the relaxation times together, within the range
of stable values from 0.2 to 1.5, which corresponds to viscosities of 1.5 to 0.06. Everything
else was held constant. The accuracy improves with smaller values of λ˜, however Ma also
decreases, and it is not clear if this is the source of the improvement in accuracy.
Figure 6.7 shows the range of stability for λξ, at different values of all the other relaxation
frequencies (where λ8 = λ7 = λ6 = λν). The lowest value of λξ which is stable is 0.35
times λν , and therefore bulk viscosity effects will always be present in any solution. If
the Vincenti and Kruger (1976) definition of bulk viscosity is used, obtained by setting
λν = λξ/3 as discussed in the text following Equation (4.51), the model is not stable for
zero bulk viscosity.
74
P
er
ce
n
ta
g
e
er
ro
r
ν˜
0 0.5 1 1.5
10−6
10−5
10−4
10−3
10−2
Figure 6.6: Percentage error at different values of viscosity. Pressure difference driving the flow,
in non-dimensional units, is 10−4.
λ
ξ
λ
ν
0.5 1 1.5
0
0.5
1
1.5
2
Figure 6.7: The white region represents stable values of λ˜ξ at different values of λ˜ν . The pressure
difference driving the flow is 10−4 in non-dimensional units.
75
Channel Length (L) 1× 10−5 m
Channel Height (H) 1× 10−5 m
Average Pressure 1× 105 Pa
Total Pressure drop 10 Pa
Temperature (T) 300 K
Gas Constant (R) 4157 J kmol−1
Viscosity of Hydrogen (at 25◦C) (µ) 8.81× 10−6 kg m−1 s−1
Table 6.1: Parameters for dimensional Poiseuille problem
6.2.1 Higher moment relaxation frequencies
Changing λ6 and λ7 has very little effect on the results. In fact, it is only when these
relaxation frequencies are very small that they reduce the accuracy of the solution. There
is a stability limit on high and low values of λ6 and λ7 but, as these relaxation frequencies
play no part in the hydrodynamic equations, avoiding unstable values is easy.
Changes in λ8 produce an extremely small change in the accuracy of the solution. For
quite a wide range of values of λν , λ8 can be set to zero without affecting accuracy, for
values of λν outside this range λ8 has a stability boundary at a very small number. For
all values of λν , λ8 has a stability boundary at high values. Typically λ8 cannot be more
than 3/2 times λν .
In view of this situation, it seems that setting λ6 = λ7 = λ8 = λν is the best first
approximation. If different values are required, care should be taken, when working on
dimensional problems, that the value of the higher relaxation frequencies are set in the
fluid dynamic scale, not the computational scale, as discussed by Dellar (2003).
6.3 A dimensional problem
It is useful to examine a dimensional problem, to clarify the accuracy which might be
expected from this method. Poiseuille flow of hydrogen was simulated, with the parameters
shown in Table 6.3. The analytical solution is obtained by integrating Equation (6.10).
76
c
Nx
30 60 90
1300 -0.9427 -0.9434 -0.9436
2000 -0.9427 -0.9434 -0.9436
3000 -0.9427 -0.9434 -0.9436
4000 -0.9427 -0.9434 -0.9436
Table 6.2: Channel-averaged x velocity obtained from simulations [m/s]
Nx
30 60 90
0.9448 0.9456 0.9458
Table 6.3: Average x velocity (u¯x) obtained from trapezoidal integration of the analytical solution,
Equation (6.10), [m/s]
u¯x =
H2
µ12
dp
dx
= 0.94589 m/s (6.26)
Where u¯ is the cross-channel averaged velocity in the x direction. As ux 
√
RT this
problem is low Mach number. Solutions were obtained for different numbers of nodes (Nx)
and grid speed (c). The mean velocity was calculated from the numerical results using the
trapezium rule. For values of c below 1300 m/s the simulation was unstable. The results
at larger values of c are shown in table 6.2.
As can be seen, increasing c has no effect on the results to four significant figures.
The main cause of error is in the use of the trapezium rule, and the value of Kn. To
differentiate between these errors, trapezoidal integration of the analytical solution was
carried out, giving the results shown in table 6.3.
Generally, standard LBM (where R˜T = 1/3) is stable at lower values of c, and so
an increase in accuracy can be observed as c increases, consistent with reducing Ma. In
dimensional problems like this, where ˜(RT ) has a meaningful value, the requirement to have
R˜T ≤ 1/3 ensures that Ma is already so small as to make further increases in accuracy
insignificant.
Also recorded were the number of site updates, which is the number of iterations needed
77
c
Nx
30 60 90
1300 1.89×107 9.72×107 2.92×108
2000 1.80×107 1.37×108 4.29×108
3000 2.61×107 1.94×108 6.08×108
4000 3.33×107 2.45×108 8.26×108
Table 6.4: Number of site updates required to reach a converged solution
to reach convergence multiplied by the number of grid points. This is a good indication
of the computational time required, regardless of implementation and hardware. It can
be seen that increasing the value of c increases the computational time needed, without
significantly increasing accuracy, and therefore the lowest value of c possible should be
used.
The largest size of the grid spacing ∆x was 3 × 10−7 m and grid spacing much larger
than this was not stable. The size of the grid spacing is comparable to the mean free path
of the gas, but only the continuum equations of fluid dynamics are modelled. Therefore,
this method is not well suited for either microscopic or macroscopic modelling.
6.4 Conclusions
The moment based boundary condition method from Chapter 5 was applied to the Poiseuille
flow, using the analytical solution for guidance on the values of the moments at the bound-
aries. At the wall the velocity tangential and normal to the wall was specified, and the
viscous normal stress in the direction tangential to the wall was set to zero. At the open
boundary the density was specified and the velocity tangential to the boundary and the
viscous normal stress tangential to the boundary were set to zero. At the intersection
between a wall and an open boundary all these conditions were imposed simultaneously.
These boundary conditions are easy to implement and involve no significant increase in
computational cost compared with others, such as the bounce-back boundary condition.
Non-dimensional results were obtained for the Poiseuille flow, which highlighted some
78
features of the LB model. Significantly, the model is not stable for values of R˜T & 1/3
and the stability is highest when R˜T = 1/3. The accuracy improves as Ma is reduced
until the limit of the computer’s numerical range is reached. Changes to the bulk viscosity,
controlled by the relaxation time λξ, cause instability outside of a relatively small region
near λξ = λν . This region is so narrow that the theoretical solution for a simple monatomic
gas cannot be recovered.
Finally, the solution to a dimensional problem, Poiseuille flow of hydrogen, was obtained.
The Ideal Gas equation of state was included in the model. To ensure that the value of R˜T
is less than 1/3 the value of c must be large, thus the value of Ma is small and the solution
is sufficiently accurate for most purposes. Further increases in c increase the accuracy, but
there are few applications where that level of accuracy would be required. Finally, the
maximum physical grid spacing for which the scheme is stable is small, around 3 × 10−7
m.
79
Chapter 7
LB for Multi-Component Gases
The method developed in the previous chapters will now be extended to multi-component
gas mixtures. As discussed in Section 3.1.6, most authors seek to split the collision step
in the particle basis. Instead, the approach outlined here is to use the multiple relaxation
time model to split the collision step in the moment basis, permitting a different relaxation
frequency to be used for the diffusive effects and the viscous effects. The concept of a
viscosity of each species within a mixture is occasionally mentioned in the literature, for
example Asinari (2005). However, as it is poorly defined and not obviously useful it will
not be considered in this work. Instead the viscosity is treated as a property of the mixture.
This chapter consists of a brief introduction to the nomenclature used, and an analytical
solution for binary channel flow. Three lattice Boltzmann models are defined, analysed and
compared, and the most suitable model is selected. Finally the solution of a dimensional
binary channel flow problem is investigate. Unfortunately, this shows accuracy which is
only first order.
80
7.1 Multi-component nomenclature
In a multi-component mixture each species is identified with a subscript σ, and the mixture
pressure, density and velocity are defined as,
p =
∑
σ
pσ ρ =
∑
σ
ρσ u =
∑
σ ρσuσ
ρ
(7.1)
The term mixture velocity will always refer to this, the mass averaged rather than mole
averaged, velocity. The mole and mass fractions, Xσ and Yσ respectively, are defined as
Xσ =
pσ
p
Yσ =
ρσ
ρ
(7.2)
The equation of state for each species is,
pσ = ρσ
R0T
Mσ
= ρσ(RσT ) (7.3)
where R0 is the universal molar gas constant. The equation of state for the mixture is
p = ρ
R0T
M
= ρ(RT ) (7.4)
Where M is the mole-averaged molar mass.
7.2 Binary channel flow theory
The flow of two gases through a channel with fixed partial pressures at both ends is the
binary equivalent of Poiseuille flow and will now be analysed. The following simplifications
will be employed,
uˆσ,y = 0
∂ρˆσ
∂yˆ
= 0
∂pˆσ
∂yˆ
= 0
∂
∂tˆ
= 0 (7.5)
81
Each species obeys the steady-state mass continuity equation,
∂(ρˆσuˆσ,x)
∂xˆ
= −∂(ρˆσuˆσ,y)
∂yˆ
= 0 (7.6)
The general solution for channel flow of the mixture must be the same as that of a single-
component fluid, namely the solution of,
∂pˆ
∂xˆ
=
∂
∂y
(
νˆρˆ
∂uˆx
∂y
)
(7.7)
Where νˆ is a suitable scaled viscosity of the mixture. Equation (7.7) may be integrated
using a given slip velocity, uˆslip, at the wall. The value of the slip velocity is known from
diffusion slip theory, and will be discussed in Section 7.6. ρˆ and νˆ are uniform in the y
direction, and can be taken outside the derivative. The result is,
uˆx = uˆslip +
(yˆ2 − (Hˆ/2)2)
2ρˆνˆ
∂pˆ
∂xˆ
(7.8)
This equation shows that the mixture velocity is parabolic, shifted by the slip velocity at
the wall.
Chapman and Cowling (1970) give the following general equation for diffusion of a
binary mixture,
(uA,α − uB,α) =
−1
XAXBBA,B
[
∂XA
∂xˆα
+
nAnB(mB −mA)
nρ
∂(ln p)
∂xˆα
− ρAρB
pρ
(FA,α − FB,α) + kT ∂(lnT )
∂xˆα
]
(7.9)
where Fσ,α is a force acting on the species σ in direction of α and kT is the thermal diffusion
ratio. The last two terms are not relevant to this work, but are included for completeness.
mσ is the mass of a molecule of species σ and Bσ,ς is the diffusion resistance between
species σ and ς. The diffusion resistance is the inverse of the binary diffusion coefficient
(Bσ,ς = 1/Dσ,ς). Clearly the diffusion is governed by four physical effects, the concentration
82
gradient of the species, the pressure gradient of the mixture, the body forces on the different
species and the temperature gradient. The second term is known as ‘pressure diffusion’. A
good example of where it is important is in the separation of a gas mixture rotated around
an axis, such as in a centrifuge. The pressure gradient causes the species to separate until
the concentration gradients cause diffusion in the opposite direction and a steady state is
reached. In the applications considered in this work there are no body forces on the gases
apart from gravity, and where species of different mass are considered this means that the
third term is non-zero. However it will be small and may be safely ignored. The assumption
is made that the thermal term of Equation (7.9) is zero, as the model is isothermal.
If ‘pressure diffusion’ is ignored the well known Stefan-Maxwell equation is recovered,
which can be written as,
ρˆ
∂Yσ
∂xˆα
= Bˆσ,ς ρˆσ(uˆα − uˆσ,α) (7.10)
Equation (7.10) shows that in channel flow, where the species density is constant across
the channel, the species velocity takes a parabolic shape just like the mixture velocity
but shifted by a constant value. Using Equation (7.8) for the mixture velocity in Equa-
tion (7.10) results in an expression for the species velocity,
uˆσ,α = uˆslip +
(yˆ2 − (Hˆ/2)2)
2ρˆνˆ
∂pˆ
∂xˆα
− ρˆ
Bˆσ,ς ρˆσ
∂Yσ
∂xˆα
(7.11)
The species velocity is therefore determined by three factors; the slip velocity of the mixture
at the wall, the mixture pressure gradient, and the species concentration gradient. Only
the term involving the pressure gradient is a function of y. Figure 7.1 illustrates typical
velocity profiles for binary channel flow.
Constant total pressure
If the total pressure is constant everywhere, Equation (7.11) reduces to,
uˆσ,α = uˆslip − ρˆ
Bˆσ,ς ρˆσ
∂Yσ
∂xˆα
(7.12)
83
Slip Velocity
Velocity
Dis
tan
ce
 fro
m 
wa
ll
0xy
Figure 7.1: Typical velocity profiles for binary channel flow. The parabolic shape of the species
and mixture velocities is due to the total pressure gradient. The slip velocity and the diffusion
velocity are independent of y.
Using (7.12) and noting that if the total pressure is constant ∂YA/∂xα = −∂YB/∂xα, the
mixture velocity can be written as,
uˆα = uˆslip − ρˆ
ρˆBˆσ,ς
(
ρˆA
ρˆA
− ρˆB
ρˆB
)
∂YˆA
∂xˆα
= uˆslip (7.13)
The slip velocity is the only parameter that has an impact on the mixture velocity in
the case of zero total pressure drop. If both gases have the same molar mass, diffusion
slip theory predicts uˆslip = 0, and the mixture velocity uˆα is therefore zero. However,
each individual species will have a velocity at the wall according to the third term of
Equation (7.11).
7.3 Multi-component Lattice Boltzmann scheme
To define a Lattice Boltzmann scheme for multi-component gas mixtures, each species
must have its own distribution function, fσ. The lower moments of this function are the
84
species density and species velocity, in a similar way to the single-component equations.
Thus, 
Π˜σ,0
Π˜σ,X
Π˜σ,Y
 =

ρˆσ
Maρˆσuˆσ,x
Maρˆσuˆσ,y
 (7.14)
An equilibrium distribution function must also be defined. If the lower moments of the
equilibrium function are the same as the moments above, momentum will be conserved
by each species, and diffusion will not occur. If the lower moments of the equilibrium
function are centred on the mixture velocity, then the momentum of the mixture will be
conserved, but not that of each species. Therefore, the mixture velocity is used, and the
lower moments of the equilibrium function are defined as,
Π˜eqσ,0
Π˜eqσ,X
Π˜eqσ,Y
 =

ρˆσ
Maρˆσuˆx
Maρˆσuˆy
 (7.15)
It is not necessary that the equilibrium function takes this form, as long as the correct
diffusion equation can be recovered. Although no alterative forms of the lower moments of
the equilibrium function are considered in this work, that is not to say that such consider-
ations might not be profitable. Similar to (4.16) for a single-component gas, the evolution
equation for each species is,
MaKn
∂f˜σ,i
∂tˆ
+ Kne˜x,i
∂f˜σ,i
∂xˆ
+ Kne˜y,i
∂f˜σ,i
∂yˆ
= Ω˜(f˜σ,i) (7.16)
7.4 Mass conservation equation
The first moment of Equation (7.16) is,
KnMa
∂(ρˆσ)
∂tˆ
+ Kn
∂(Maρˆσuˆσ,x)
∂xˆ
+ Kn
∂(Maρˆσuˆσ,y)
∂yˆ
= 0 (7.17)
85
Thus, species mass conservation is ensured, and summing over the mixture it is trivial to
show that the mixture mass is also conserved.
7.5 Momentum equation
In Chapter 4 an approximate expression for the momentum flux density tensor was ob-
tained by assuming that the real distribution function was only slightly perturbed from
the equilibrium function. This is not strictly true anymore, as the mixture velocity can be
very different from the species velocity. In fact, this is an area which could certainly benefit
from further work. Lacking such work however, it is assumed the method developed in
Chapter 4 is still valid.
The right hand side of the momentum conservation equation is no longer zero, but
instead is given by,
λ˜δ(Π˜(eq)α − Π˜α) = 2MaKnλˆδρˆσ(uˆα − uˆσ,α) (7.18)
where λ˜δ = λ˜1 = λ˜2 is introduced as the relaxation frequency controlling diffusion. The
right hand side is obtained using the definitions of the lower moments from Section 7.3.
It is not possible to write the momentum equation for a single species in a mixture in the
form of the corresponding conservation equation for a simple gas, see Vincenti and Kruger
(1976). An approximation must be made but the direction to take is not clear. Hence,
three plausible possibilities for the equilibrium distribution functions are considered. The
first, and perhaps the most intuitive, is the Mixture Centred Equilibrium (MCE) where the
higher moments are centred on the mixture velocity. The second is the Species Centred
Equilibrium (SCE) where the higher moments are centred on the species velocity. The
third is the Corrected Species Centred Equilibrium (CSCE) and will be explained later.
86
7.5.1 Mixture Centred Equilibrium (MCE)
If all the moments of the equilibrium function are centred on the mixture velocity, the
higher moments are,
Π˜eqσ,XX
Π˜eqσ,Y Y
Π˜eqσ,XY
Π˜eqσ,XY Y
Π˜eqσ,XXY
Π˜eqσ,XXY Y

=

2Ma(pˆσ + ρˆσ(uˆx)
2)
2Ma(pˆσ + ρˆσ(uˆy)
2)
2Maρˆσuˆxuˆy
Maρˆσuˆx/3
Maρˆσuˆy/3
2Ma
1
3
(pˆσ + (ρˆσuˆx)
2 + (ρˆσuˆy)
2)

(7.19)
An approximate expression for the momentum flux density tensor is then obtained following
the approach of Section 4.3.2.
MaKn
∂
∂tˆ

2Ma(pˆσ + ρˆσuˆ
2
x)
2Ma(pˆσ + ρˆσuˆ
2
y)
2Ma(ρˆσuˆyuˆx)
+ Kn ∂∂xˆ

Maρˆσuˆx
Ma
1
3
ρˆσuˆx
Ma
1
3
ρˆσuˆy
+ Kn ∂∂yˆ

Ma
1
3
ρˆσuˆy
Maρˆσuˆy
Ma
1
3
ρˆσuˆx

≈

λ˜ν(Π˜eqσ,XX − Π˜σ,XX)
λ˜ν(Π˜eqσ,Y Y − Π˜σ,Y Y )
λ˜ν(Π˜eqσ,XY − Π˜σ,XY )
 (7.20)
There are three different velocities that enter this expression; the velocity from the first
order moment (like Π˜eqσ,X = ρ˜σu˜x), the velocity from the third order moment ( like Π˜
eq
σ,XXX =
Π˜eqσ,X = ρ˜σu˜x and Π˜
eq
σ,XY Y = ρ˜σu˜x/3), and the velocity from the continuity equation. In
this case these are the mixture, mixture and species velocities respectively. To reduce
the algebraic complexity λν = λξ is assumed or, equivalently, that the collision matrix is
diagonal. The result can be easily extended to the more general collision operator discussed
earlier.
Repeating the analysis of Chapter 4, the time derivative terms are ignored except those
87
which (via the equation of state and the species mass conservation Equation (7.17)) contain
˜(RT )σ. The following expressions for the components of the momentum flux density tensor
are obtained,

Π˜σ,XX
Π˜σ,Y Y
Π˜σ,XY
 ≈

2Ma(pˆσ + ρˆσuˆ
2
x))
2Ma(pˆσ + ρˆσuˆ
2
y))
0
+ KnMa 1λ˜ν ∂∂xˆ

ρˆσ(uˆx − ˜(RT )σuˆσ,x)
ρˆσ(
1
3
uˆx − ˜(RT )σuˆσ,x)
1
3
ρˆσuˆy

+ KnMa
1
λ˜ν
∂
∂yˆ

ρˆσ(
1
3
uˆy − ˜(RT )σuˆσ,y)
ρˆσ(uˆy − ˜(RT )σuˆσ,y)
1
3
ρˆσuˆx
 (7.21)
The comparison between Equation (7.21) and Equation (4.44) is particularly interesting,
due to the mix of species and mixture velocities in Equation (7.21). The necessary, but not
sufficient, criteria for stability from Equation (4.44) was that R˜T . 1/3 (to ensure that
the error terms in Equation (4.44) did not dominate). This condition is now replaced with
the condition ˜(RT )σu˜σ,α . u˜α/3. Because this contains the species and mixture velocities
a maximum stable value of ˜(RT )σ cannot be determined. However small values of
˜(RT )σ
are more likely to be stable than large values.
The non-equilibrium parts of the moments are defined as Π
(neq)
αβ = Παβ − Π(eq)αβ , and
will provide a useful shorthand in what follows. From Equation (7.21) it is clear that∑
σ Π
(neq)
σ,αβ = Π
(neq)
αβ , where Π
(neq)
αβ is the momentum flux density tensor for the mixture.
Substituting the momentum flux density tensor into the species momentum equation, which
is the multi-component version of Equation (4.46) gives,
∂
∂tˆ
 ρˆσuˆσ,x
ρˆσuˆσ,y
+ ∂
∂xˆ
 ρˆσuˆ2x
ρˆσuˆxuˆy
+ ∂
∂yˆ
 ρˆσuˆxuˆy
ρˆσuˆ
2
y
 = −
 ∂∂xˆ(pˆσ)
∂
∂yˆ
(pˆσ)

+
∂
∂xˆ
 Π(neq)σ,XX
Π
(neq)
σ,XY
+ ∂
∂yˆ
 Π(neq)σ,XY
Π
(neq)
σ,Y Y
+
 λˆδρˆσ(uˆx − uˆσ,x)
λˆδρˆσ(uˆy − uˆσ,y)
 (7.22)
88
Note the factor of 2MaKn from Equation (7.18) cancels with terms in the momentum
equation. Summing over the species results in the momentum equation for the mixture,
∂
∂tˆ
 ρˆuˆx
ρˆuˆy
+ ∂
∂xˆ
 ρˆuˆ2x
ρˆuˆxuˆy
+ ∂
∂yˆ
 ρˆuˆxuˆy
ρˆuˆ2y
 = −
 ∂∂xˆ(pˆ)
∂
∂yˆ
(pˆ)

+
∂
∂xˆ
 Π(neq)XX
Π
(neq)
XY
+ ∂
∂yˆ
 Π(neq)σ,XY
Π
(neq)
σ,Y Y
 (7.23)
This is the same equation as for a single-component fluid, Equation (4.52), and the momen-
tum flux density tensor is the same as Equations (4.44) and (4.45) but with the mixture
velocity and density. Multiplying Equation (7.23) by the mass fraction, Yσ, and subtracting
Equation (7.22) results in,
Yσ
∂
∂xˆα
(pˆ)− ∂
∂αˆ
(pˆσ) = (Yσ −Xσ) ∂
∂xˆα
(pˆ)− pˆ ∂
∂xˆα
(Xσ) = −λˆδρˆσ(uˆα − uˆσ,α) (7.24)
An expression linking the diffusion resistance to the relaxation frequency is now needed.
The expression given by Chapman and Cowling (1970) for the continuum BGK equation
is,
λˆδ =
ρˆkˆT
nˆmˆσmˆς
Bˆσς = (Xσ ˆ(RT )ς +Xς
ˆ(RT )σ)Bˆσς (7.25)
where k is the Boltzmann constant. Fortunately, this expression can be used without
modification in the D2Q9 Lattice Boltzmann scheme. Substituting Equation (7.25) into
Equation (7.24) results in,
(Yσ −Xσ) ∂
∂xˆα
(pˆ)− pˆ ∂
∂xˆα
(Xσ) =
kˆT nˆσnˆςBˆσς
nˆ
(uˆσ,α − uˆς,α) (7.26)
Comparison between Equation (7.9) and a manipulated form of Equation (7.26) shows that
the correct form for ordinary diffusion and pressure diffusion are exactly recovered.
The relaxation frequency related to viscosity is defined as it was in the single-component
89
case, Equation (4.50),
λ˜νσ =
1
3ν˜
(7.27)
7.5.2 Species Centred Equilibrium (SCE)
The second alternative equilibrium function is the Species Centred Equilibrium function,
for which the higher moments are defined as,

Π˜eqσ,XX
Π˜eqσ,Y Y
Π˜eqσ,XY
Π˜eqσ,XY Y
Π˜eqσ,XXY
Π˜eqσ,XXY Y

=

2Ma(pˆσ + ρˆσ(uˆσ,x)
2)
2Ma(pˆσ + ρˆσ(uˆσ,y)
2)
2Maρˆσuˆσ,xuˆσ,y
Maρˆσuˆσ,x/3
Maρˆσuˆσ,y/3
2Ma
1
3
(pˆσ + (ρˆuˆσ,x)
2 + (ρˆuˆσ,y)
2)

(7.28)
Using this equilibrium function the equivalent to Equation (7.20) is,
MaKn
∂
∂tˆ

2Ma(pˆσ + ρˆσuˆ
2
σ,x)
2Ma(pˆσ + ρˆσuˆ
2
σ,y)
2Ma(ρˆuˆyuˆx)
+ Kn ∂∂xˆ

Maρˆσuˆx
Ma
1
3
ρˆσuˆσ,x
Ma
1
3
ρˆσuˆσ,y
+ Kn ∂∂yˆ

Ma
1
3
ρˆσuˆσ,y
Maρˆσuˆy
Ma
1
3
ρˆσuˆσ,x

≈

λ˜ν(2Ma(pˆσ + ρˆσ(uˆσ,x)
2)− Π˜σ,XX)
λ˜ν(2Ma(pˆσ + ρˆσ(uˆσ,y)
2)− Π˜σ,Y Y )
λ˜ν(2Maρˆσuˆσ,xuˆσ,y − Π˜σ,XY )
 (7.29)
90
which can be rearranged as,

Π˜σ,XX
Π˜σ,Y Y
Π˜σ,XY
 ≈

2Ma(pˆσ + ρˆσuˆ
2
σ,x)
2Ma(pˆσ + ρˆσuˆ
2
σ,y)
0
+ KnMa 1λ˜ν ∂∂xˆ

ρˆσ(uˆx − ˜(RT )σuˆσ,x)
ρˆσ(
1
3
− ˜(RT )σ)uˆσ,x
1
3
ρˆσuˆσ,y

+ KnMa
1
λ˜ν
∂
∂yˆ

ρˆσ(
1
3
− ˜(RT )σ)uˆσ,y)
ρˆσ(uˆy − ˜(RT )σuˆσ,y)
1
3
ρˆσuˆσ,x
 (7.30)
The momentum flux density tensor still contains both mixture and species velocities, and
both the stability criteria R˜T σ . 1/3 and ˜(RT )σu˜σ,α . u˜α/3 now apply. It is not clear if
this is an advantage over the MCE or not. The species momentum equation is,
∂
∂tˆ
 ρˆσuˆσ,x
ρˆσuˆσ,y
+ ∂
∂xˆ
 ρˆσuˆ2σ,x
ρˆσuˆσ,xuˆσ,y
+ ∂
∂yˆ
 ρˆσuˆσ,xuˆσ,y
ρˆσuˆ
2
σ,y
 = −
 ∂∂xˆ(pˆσ)
∂
∂yˆ
(pˆσ)

+
∂
∂xˆ
 Π(neq)σ,XX
Π
(neq)
σ,XY
+ ∂
∂yˆ
 Π(neq)σ,XY
Π
(neq)
σ,Y Y
+
 λˆδρˆσ(uˆx − uˆσ,x)
λˆδρˆσ(uˆy − uˆσ,y)
 (7.31)
and the momentum equation for the mixture is,
∂
∂tˆ
 ρˆuˆx
ρˆuˆy
+ ∂
∂xˆ
 ∑σ ρˆσuˆ2σ,x∑
σ ρˆσuˆσ,xuˆσ,y
+ ∂
∂yˆ
 ∑σ ρˆσuˆσ,xuˆσ,y∑
σ ρˆσuˆ
2
σ,y
 =
−
 ∂∂xˆ(pˆ)
∂
∂yˆ
(pˆ)
+ ∂
∂xˆ
 Π(neq)XX
Π
(neq)
XY
+ ∂
∂yˆ
 Π(neq)XY
Π
(neq)
Y Y
 (7.32)
The advection terms in this equation are not correct. This is a serious problem and the
Corrected Species Centred Equilibrium (CSCE) is introduced to address this difficulty.
91
7.5.3 Corrected Species Centred Equilibrium (CSCE)
The Corrected Species Centred Equilibrium function consists of the SCE with a correc-
tion to the equilibrium momentum flux density tensor such that the equilibrium function
becomes, 
Π˜eqσ,XX
Π˜eqσ,Y Y
Π˜eqσ,XY
Π˜eqσ,XY Y
Π˜eqσ,XXY
Π˜eqσ,XXY Y

=

2Ma(pˆσ + ρˆσuˆσ,xuˆx)
2Ma(pˆσ + ρˆσuˆσ,yuˆy)
2Maρˆσ
1
2
(uˆσ,xuˆy + uˆσ,yuˆx)
Maρˆσuˆσ,x/3
Maρˆσuˆσ,y/3
2Ma
1
3
(pˆσ + (ρˆuˆσ,x)
2 + (ρˆuˆσ,y)
2)

(7.33)
The non-equilibrium momentum flux density tensor (the viscous stress tensor) is the same
as that recovered from the SCE, specifically

Π˜(neq)σ,XX
Π˜(neq)σ,Y Y
Π˜(neq)σ,XY
 ≈ KnMa 1λ˜ν ∂∂xˆ

ρˆσ(uˆx − ˜(RT )σuˆσ,x)
ρˆσ(
1
3
− ˜(RT )σ)uˆσ,x
1
3
ρˆσuˆσ,y

+ KnMa
1
λ˜ν
∂
∂yˆ

ρˆσ(
1
3
− ˜(RT )σ)uˆσ,y)
ρˆσ(uˆy − ˜(RT )σuˆσ,y)
1
3
ρˆσuˆσ,x
 (7.34)
This has the same stability criteria as the SCE. The species momentum equation is,
∂
∂tˆ
 ρˆσuˆσ,x
ρˆσuˆσ,y
+ ∂
∂xˆ
 ρˆσuˆσ,xuˆx
ρˆσ
1
2
(uˆσ,xuˆy + uˆσ,yuˆx)
+ ∂
∂yˆ
 ρˆσ 12(uˆσ,xuˆy + uˆσ,yuˆx)
ρˆσuˆσ,yuˆy
 =
−
 ∂∂xˆ(pˆσ)
∂
∂yˆ
(pˆσ)
+ ∂
∂xˆ
 Π(neq)σ,XX
Π
(neq)
σ,XY
+ ∂
∂yˆ
 Π(neq)σ,XY
Π
(neq)
σ,Y Y
+
 λˆδρˆσ(uˆx − uˆσ,x)
λˆδρˆσ(uˆy − uˆσ,y)
 (7.35)
The strange advection terms in the species momentum equation are the price paid for the
correct advection terms in the mixture momentum equation. Summing over the species
92
gives the momentum equation for the mixture,
∂
∂tˆ
 ρˆuˆx
ρˆuˆy
+ ∂
∂xˆ
 ρˆuˆ2x
ρˆuˆxuˆy
+ ∂
∂yˆ
 ρˆuˆxuˆy
ρˆuˆ2y
 =
−
 ∂∂xˆ(pˆ)
∂
∂yˆ
(pˆ)
+ ∂
∂xˆ
 Π(neq)XX
Π
(neq)
XY
+ ∂
∂yˆ
 Π(neq)XY
Π
(neq)
Y Y
 (7.36)
Multiplying Equation (7.36) by the mass fraction, Yσ, and subtracting Equation (7.35)
results in terms of the form ρσu(u−uσ). Thus, the Stefan-Maxwell equation is not recovered
at all times. However this is also true for continuum kinetic theory, where the Stefan-
Maxwell equation is only recovered when all velocity gradients giving rise to viscous stresses
are zero.
7.6 Expressions for diffusion slip
In Section 3.2.1, the Kramers and Kistemaker (1943) expression for diffusion slip was
derived from simple kinetic theory.
∑
σ
ρσuσ,x√
Mσ
= 0 (7.37)
This expression is valid along a wall where the x direction is parallel to the wall. If the
denominator of this expression is instead unity, a zero mixture velocity at the wall is
obtained. It is important to note that this is not the same as the zero species velocity
condition usually used in LB. In order to cover both possibilities a function of molecular
weight, S(Mσ), is introduced, where either S(Mσ) = 1 or S(Mσ) =
√
Mσ. Thus,
∑
σ
ρσuσ,x
S(Mσ)
= 0 (7.38)
This expression relates one species velocity to the other, but the absolute value of either
species velocity cannot be determined. A further condition is therefore required. As shown
93
in Equation (7.11) the species velocity is influenced by three terms. Across a channel the
term containing the slip velocity and the diffusion term are both constant and only the
viscous term varies. Therefore,
∂ux
∂y
=
∂uσ,x
∂y
(7.39)
Equation (7.38) and (7.39) prove to be sufficient conditions for the species velocities tan-
gential to the wall to be determined. The velocity of each species normal to the wall is
zero, uσ,y = 0.
Kramers and Kistemaker derived an expression based on these equations, and the
equipartition of energy, for a binary gas with MB > MA.
ux,slip =
√
MB −
√
MA
pA
√
MA + pB
√
MB
DAB
∂pA
∂x
(7.40)
Clearly, the slip velocity will take the same sign as the gradient of the lighter gas concen-
tration.
7.6.1 Implementing diffusion slip in the LB method
FluidFluidoutside
domain
xy
Solid
Fluid
xy
Fluid
Fluid
outside
domain
Solid
xy
a b c
Figure 7.2: Summary of the geometry of boundary nodes repeated from Chapter 5. a is an open
boundary, b is a wall, and c is a corner node.
In this section a diffusion slip condition will be implemented on a flat wall, with unknown
f values f2, f5 and f6 as shown in Figure 7.2b. Equation (7.39) must be transformed into
94
a condition for the XY component of the momentum flux density tensor of each species.
The viscous shear stress, in multicomponent channel flow, based on the SCE or CSCE
models (see Equation (7.30)) is,
Π˜σ,XY ≈ −2Maνˆρˆσ
∂uˆσ,x
∂yˆ
(7.41)
If the MCE (Equation (7.21)) is used the viscous shear stress is,
Π˜σ,XY ≈ −2Maνˆρˆσ
∂uˆx
∂yˆ
(7.42)
In fact the analysis is the same whichever equilibrium function is used. Equation (7.39)
becomes,
1
ρA
(ΠA,XY ) =
1
ρB
(ΠB,XY ) (7.43)
1
ρA
(f5,A − f6,A + f7,A − f8,A) = 1
ρB
(f5,B − f6,B + f7,B − f8,B) (7.44)
The species mass flux in the x direction is given by the equation,
fσ,5 − fσ,6 = ρσuσ,x − fσ,1 − fσ,8 + fσ,3 + fσ,7 (7.45)
To eliminate the unknown f values in Equation (7.44), substitute in Equation (7.45) to
give,
1
ρA
(ρAuA,x− f1,A− 2f8,A + f3,A + 2f7,A) = 1
ρB
(ρBuB,x− f1,B − 2f8,B + f3,B + 2f7,B) (7.46)
For convenience a new variable composed only of known f values is introduced, κσ =
−f1,σ − 2f8,σ + f3,σ + 2f7,σ. The density at the wall is already known, and is given by,
ρσ = fσ,0 + fσ,1 + fσ,3 + 2(fσ,2 + fσ,6 + fσ,5) (7.47)
95
Substituting in the other boundary condition necessary for diffusion slip, Equation (7.38),
gives,
1
ρA
(ρAuA,x + κA) =
1
ρB
(
−S(MB)
S(MA)
(ρAuA,x) + κB
)
(7.48)
which can be rearranged as
(ρAuA,x)
(
1
ρA
+
S(MB)
ρBS(MA)
)
=
1
ρB
(κB)− 1
ρA
(κA) (7.49)
This is an expression for the velocity to be set at the wall, and it can be used to set the
unknown f values. There is a similar expression for species B, obtained by swapping the
B and A subscripts.
The procedure described in Chapter 6 is then used for setting this velocity at the wall,
with the other conditions being u˜σ,y = 0 and Π˜σ,XX = p˜σ+ρ˜σu˜
2
x (MCE), Π˜σ,XX = p˜σ+ρ˜σu˜
2
σ,x
(SCE) or Π˜σ,XX = p˜σ + ρ˜σu˜σ,xu˜x (CSCE).
7.6.2 Intersections between walls and open boundaries
In a corner between an open boundary and a wall there are 5 unknowns, fσ,1, fσ,2, fσ,5,
fσ,6 and fσ,8, as shown in Figure 7.2c. The density is the same as at the open boundary
and the velocity and viscous shear stress are determined by the conditions for a solid wall.
The definitions of the moments, Πσ,0, Πσ,X , Πσ,Y = 0 and Πσ,Y Y = ρσ(RT )σ can be written
with the unknown fσ values on the left hand side,
fσ,1 + fσ,2 + fσ,5 + fσ,6 + fσ,8 = ρ− fσ,3 − fσ,4 − fσ,7 − fσ,0 (7.50)
fσ,1 + fσ,5 + fσ,8 − fσ,6 = ρσuσ,x + fσ,3 + fσ,7 (7.51)
fσ,2 + fσ,6 + fσ,5 − fσ,8 = fσ,7 + fσ,4 (7.52)
fσ,5 + fσ,6 + fσ,2 + fσ,8 = ρσ(RT )σ − fσ,7 − fσ,4 (7.53)
96
Combining these expressions allows the elimination of the unknowns in the XY component
of the species momentum flux density tensor, resulting in,
Πσ,XY = (4f7,σ + f0,σ + 2f3,σ + 2f4,σ − ρσ + ρσux,σ) (7.54)
For the sake of convenience, a new variable is introduced, Θσ = 4f7,σ + f0,σ + 2f3,σ +
2f4,σ, which is composed only of the known f values. Substituting Equation (7.54) into
Equation (7.43) yields,
(ΘA − ρA + ρAux,A) 1
ρA
= (ΘB − ρB + ρBux,B) 1
ρB
(7.55)
Finally using Equation (7.38),
(ΘA − ρA + ρAux,A) 1
ρA
=
(
ΘB − ρB − ρAux,AS(MB)
S(MA)
)
1
ρB
(7.56)
ux,A
(
1 +
ρA
ρB
S(MB)
S(MA)
)
=
1
ρB
(ΘB − ρB)− 1
ρA
(ΘA − ρA) (7.57)
This x velocity is to be set at the corner and can be imposed using the equations described
in Section 6.1.1.
xy
Figure 7.3: The boundary conditions at the wall, the open boundary and the intersection between
the two for binary channel flow.
Figure 7.3 shows a summary of all the boundary conditions in the moment basis for the
97
binary channel flow problem.
7.7 Modified relaxation times
Following the approach of Section 5.2, a modified diffusion resistance is defined to take
account of second order numerical effects. Taking the x direction as an example the
equation that will be implemented is,
Π˜+X − Π˜X = λ˜δ′(Π˜(eq)X − Π˜X) (7.58)
Compare Equation (7.58) above with a similar one derived using central differencing,
(Π˜+X − Π˜X) =
λ˜δ
2
(
(Π˜
(eq)
X − Π˜X) + (Π˜(eq)+X − Π˜+X)
)
(7.59)
Again assume Π
(eq)+
X = Π
(eq)
X , leading to,
(Π˜+X − Π˜X)
(
1 +
λ˜δ
2
)
= λ˜δ(Π˜
(eq)
X − Π˜X) (7.60)
For Equation (7.60) to be identical to Equation (7.58), the modified relaxation frequency
must be,
λ˜δ′ =
λ˜δ
(1 + λ˜
δ
2
)
(7.61)
In all the simulations performed the higher relaxation times were λ6′ = λ7′ = λδ′ and
λ8′ = λν′.
7.8 Self-diffusion test case
The simplest test for a multicomponent model is channel flow of a binary gas with both
components having the same molar mass. This is called self-diffusion. The two species
can be thought of as two colours or isotopes of the same gas, with a high concentration
98
of one colour at one end. Obviously, the mixture should behave in the same way as for a
single-component gas with the same molar mass. However, if one considers the values of
ma for each species, they may well be higher than the mixture value, and so the error may
be higher than for an equivalent single-component gas.
The wall boundary condition predicted by diffusion slip theory is zero mixture velocity.
However, each species is allowed to slip at the wall.
Simulations were performed with zero total pressure drop at different values of ν and
B, with the three different equilibrium functions. The average mixture velocity across
the channel was recorded, and also the sum of the diffusion velocities, which is defined
as
∑
σ ρσ(uσ,x − ux). By definition the sum of the diffusion velocities should always be
zero, but this was not observed exactly. The value of either ν or B should not affect
the result, which should be one of zero mixture velocity everywhere. Figure 7.4 shows the
mixture velocities recorded. As the value of B is increased, ma (if defined using the species
velocity) decreases, and the results become more accurate. The MCE is fairly accurate,
but does not consistently improve with values of B above a certain limit. The SCE shows
a comparable accuracy, but again does not behave in a predictable way. The CSCE is of
comparable accuracy to the other two, and exhibits equally unpredictable behaviour. It
should be noted that the range of stable values of B for the MCE is far greater than the
other two. In fact, the method is still stable with the value of ma in excess of 10! This is
because the species velocity does not appear in the equilibrium function, and so the species
velocity should not strictly be used in the definition of ma. The MCE is still stable at
higher values of B but the computational time increases dramatically.
Figure 7.5 shows the sum of the diffusion velocities. The results broadly reflect the
conclusions reached from the results for the mixture velocity; at low values of B the results
are less accurate. At this stage it is not clear which equilibrium function is superior,
although the MCE has a much greater range of stable values.
Figure 7.6 and 7.7 show the same metrics but this time on a flow with a total pressure
gradient. The value of ν is fixed, while B is varied. The fact that the SCE does not recover
the correct mixture equation is obvious, and the error in mixture velocity is large. The
99
CSCE
M
ix
tu
re
V
el
o
ci
ty
B˜
SCE
M
ix
tu
re
V
el
o
ci
ty
B˜
MCE
M
ix
tu
re
V
el
o
ci
ty
B˜
2 4 6 8 10 122 4 6 8 10 120.01 1 100
10−19
10−18
10−17
10−16
10−15
10−14
10−13
10−12
10−19
10−18
10−17
10−16
10−15
10−14
10−13
10−12
10−19
10−18
10−17
10−16
10−15
10−14
10−13
10−12
Figure 7.4: Mixture velocity recorded with different equilibrium functions and viscosities for self
diffusion with no total pressure gradient. The dotted line is for ν˜ = 25 and the solid line is for
ν˜ = 1.
100
CSCE
S
u
m
o
f
D
iff
u
s
io
n
V
e
lo
c
it
ie
s
B˜
SCE
S
u
m
o
f
D
iff
u
s
io
n
V
e
lo
c
it
ie
s
B˜
MCE
S
u
m
o
f
D
iff
u
s
io
n
V
e
lo
c
it
ie
s
B˜
2 4 6 8 10 122 4 6 8 10120.01 1 100
10−18
10−17
10−16
10−15
10−14
10−13
10−12
10−18
10−17
10−16
10−15
10−14
10−13
10−12
10−18
10−17
10−16
10−15
10−14
10−13
10−12
Figure 7.5: Sum of diffusion velocities recorded with different equilibrium functions and viscosities
for self diffusion with no total pressure gradient. The dotted line is for ν˜ = 25 and the solid line
is for ν˜ = 1.
101
corrections made in the CSCE function have helped slightly, and if the stability allowed
for greater values of B it is possible the accuracy would improve, In practice, however, it
does not. The MCE gives a considerably more accurate mixture velocity than the others.
In Figure 7.7 The sum of the diffusion velocities appears to show that the MCE is the
least accurate but this is actually due to the extended range of stable values of B. For
the same value of B all the equilibrium functions exhibit broadly comparable accuracy.
Worryingly, the accuracy of this metric decreases with increasing B for the MCE, while
the accuracy of the total mixture velocity increases.
CSCE
B˜
E
rr
o
r
in
M
ix
tu
re
V
el
o
ci
ty
SCE
B˜
E
rr
o
r
in
M
ix
tu
re
V
el
o
ci
ty
MCE
B˜
E
rr
o
r
in
M
ix
tu
re
V
el
o
ci
ty
2 4 6 8 10 122 4 6 8 10 120.01 1 100
10−8
10−6
10−4
10−2
100
10−8
10−6
10−4
10−2
100
10−8
10−6
10−4
10−2
100
Figure 7.6: Mixture velocity recorded with different equilibrium functions for self diffusion with
a total pressure drop of 10−4.
Real problems do not allow the luxury of changing the values of B and ν independently,
and it cannot be guaranteed that these values will fall in a regime of high accuracy. Also,
with species of different molar masses, it may be impossible to ensure that both species
will be in an accurate regime at the same time.
In conclusion, none of the equilibrium functions is entirely satisfactory. In the rest of this
102
CSCE
S
u
m
o
f
D
iff
u
s
io
n
V
e
lo
c
it
ie
s
B˜
SCE
B˜
S
u
m
o
f
D
iff
u
s
io
n
V
e
lo
c
it
ie
s
MCE
B˜
S
u
m
o
f
D
iff
u
s
io
n
V
e
lo
c
it
ie
s
2 4 6 8 10 122 4 6 8 10120.01 1 100
10−13
10−12
10−11
10−10
10−9
10−13
10−12
10−11
10−10
10−9
10−13
10−12
10−11
10−10
10−9
Figure 7.7: Sum of diffusion velocities recorded with different equilibrium functions for self diffu-
sion with a total pressure drop of 10−4.
103
work the MCE is used, because of the much greater stability range. There are also some
simplifications which offer an advantage when the advanced implementation of Chapter 8
is used.
7.9 A dimensional test case
To understand the effect of scaling, it is important to consider a dimensional test case.
Here the physical problem is fixed and then different scaling values are chosen to try to
recover the correct solution. The problem to be investigated is diffusion of hydrogen and
water vapour in a channel, involving both diffusion slip at the walls and a total pressure
gradient along the channel. The physical data is given in Table 7.1.
Channel Length 1× 10−5 m
Channel Height 1× 10−5 m
Average Pressure 1× 105 Pa
Total Pressure drop 10 Pa
Mole fraction Hydrogen at Inlet 0.55
Mole fraction Hydrogen at Outlet 0.45
Temperature 300 K
BH2,H2O 1.0789× 104 s m−2
Viscosity of Hydrogen (at 25◦C) 8.81× 10−6 kg m−1 s−1
Viscosity of Water Vapour (at 25◦C) 9.40× 10−6 kg m−1 s−1
Table 7.1: Physical data for the example problem
The viscosity of the mixture is calculated using the method of Wilke and Fuller as
described by Poling et al. (2000). Evaluating the size of the terms in Equation 7.9, the
ordinary diffusion term is more than a thousand times larger than the pressure diffusion
term for this problem. The Knudsen number is of the order of 0.01, which is small enough
that there is no viscous slip at the wall. Therefore, the results obtained from the LB
method can be compared with the numerical integration of the governing continuum fluid
flow equations, performed by the method described by Young and Todd (2005), which
neglects pressure diffusion. These results are shown in Table 7.2. Because both methods
104
are numerical in nature, exact agreement between the two is not expected.
Mixture Velocity -1.3935 m/s
Hydrogen Velocity 1.9404 m/s
Water Vapour Velocity -1.7672 m/s
Mixture Wall Slip Velocity -0.55625 m/s
Table 7.2: Solution obtained by the numerical integration method of Young and Todd (2005) for
the problem described in Table 7.1
The MCE equilibrium function was used for all calculations and the criteria for the
convergence was based on the change between time steps of the total momentum within
the domain becoming zero to machine accuracy. The values of R˜TH2 for different values
of c and Nx, the number of grid points forming one side of the square domain, are shown
in Table 7.3. All the values in the table are less than one third, and so it is likely that a
stable solution will be found for all cases considered.
c
Nx
30 60 90
2000 3.12×10−1 3.12×10−1 3.12×10−1
8000 1.95×10−2 1.95×10−2 1.95×10−2
32000 1.22×10−3 1.22×10−3 1.22×10−3
128000 7.61×10−5 7.61×10−5 7.61×10−5
512000 4.76×10−6 4.76×10−6 4.76×10−6
Table 7.3: ˜(RT )H2 calculated at different values of c and Nx
c
Nx
30 60 90
2000 3.714 7.428 1.114
8000 0.929 0.186 2.786
32000 0.232 0.464 0.696
128000 0.058 0.116 0.174
512000 0.015 0.029 0.044
Table 7.4: ν˜ calculated at different values of c and Nx
105
c
Nx
30 60 90
2000 0.072 0.036 0.024
8000 0.288 0.144 0.096
32000 1.151 0.575 0.384
128000 4.603 2.301 1.534
512000 18.412 9.206 6.137
Table 7.5: BH2,H2O calculated at different values of c and Nx
The dimensionless viscosity and diffusion resistance are shown in Tables 7.4 and 7.5.
A large value of B˜ was shown to improve the accuracy of the method in Section 7.8.
Increasing the number of grid points reduces the value of B˜, and increasing the value of c
increases it.
Mixture Velocity Hydrogen Velocity Water velocity
c
Nx Nx Nx
30 60 90 30 60 90 30 60 90
2000 -1.7258 -1.5537 -1.4964 3.6520 2.8026 2.5190 -2.3360 -2.0451 -1.9485
8000 -1.4728 -1.4274 -1.4122 2.3661 2.1587 2.0894 -1.9089 -1.8322 -1.8066
32000 -1.4094 -1.3958 -1.3911 2.0454 1.9979 1.9821 -1.8020 -1.7789 -1.7712
128000 -1.3935 -1.3879 -1.3859 1.9653 1.9577 1.9553 -1.7753 -1.7656 -1.7623
256000 NS -1.3865 -1.3850 NS 1.9510 1.9508 NS -1.7634 -1.7608
384000 NS -1.3861 -1.3847 NS 1.9488 1.9493 NS -1.7626 -1.7603
Table 7.6: Average velocities [m/s] of the mixture and species for the problem described in Table
7.1 with diffusive slip boundary condition. NS means the solution was not stable.
The results of the simulations are shown in Table 7.6. The agreement with Table 7.2 is
very poor with few grid points and low values of c but improves with higher values. The
rate of convergence is first order. The highest stable value of c is a function of the number
of grid points, and increases with number of grid points.
Also recorded was the number of site updates until the solution had converged. This
was defined as the number of grid points multiplied by the number of iterations, and this
is shown in Table 7.7. There is a disproportionate penalty to pay for a large number of
grid points as more iterations are needed as well as there being more points. An accurate
106
c
Nx
30 60 90
2000 3.690×107 2.844×108 9.477×108
8000 1.413×108 1.062×109 3.507×109
32000 6.030×108 4.136×109 1.308×1010
128000 2.428×109 1.756×1010 5.318×1010
256000 NS 3.585×1010 1.100×1011
384000 NS 5.350×1010 1.678×1011
Table 7.7: Total site updates to convergence
solution required at least 2.5× 109 site updates, taking approximately an hour on a stan-
dard desktop computer, and the accurate solutions with finer grids took 100 times longer.
Comparing this with the results for the single-component solution, which took a minimum
of 1.9 × 107 site updates to reach an accurate solution, it is clear that the extension to
multi-component mixtures has come at significant cost.
The maximum physical size of the grid spacing in these simulations was 300 nm, and this
seems close to the maximum physical size of the grid spacing for any problem. This would
certainly restrict the method to very small scale applications or very large computers.
7.9.1 Alternative boundary conditions
The test case was repeated with alternative boundary conditions of zero mixture slip ve-
locity and zero species slip velocity. To obtain zero mixture slip S(Mσ) was set to 1
in Equation (7.38). The results are shown in Table 7.8. The mixture velocity is nearly
constant with c, while the species velocity requires a much higher value of c before it is
consistent. This is because the boundary condition for the mixture does not depend on the
species velocity, as the diffusive slip condition does. Separating the mixture and species
velocity in this way makes it clear that the mixture velocity, and not the species velocity,
demonstrates a second order convergence rate.
The boundary condition that is often used for multi-component simulations in the LB
literature is that the species velocity at the wall is zero, leading to the mixture velocity
107
Mixture Velocity Hydrogen Velocity Water velocity
c
Nx Nx Nx
30 60 90 30 60 90 30 60 90
2000 -0.837 -0.835 -0.834 4.517 3.511 3.174 -1.447 -1.327 -1.287
8000 -0.837E -0.835 -0.834 2.991 2.745 2.663 -1.273 -1.241 -1.230
32000 -0.837E -0.835 -0.834 2.609 2.553 2.535 -1.229 -1.219 -1.215
128000 -0.8.37 -0.835 -0.834 2.514 2.505 2.503 -1.219 -1.213 -1.212
256000 NS -0.835 -0.834 NS 2.498 2.497 NS -1.213 -1.211
384000 NS -0.835 -0.834 NS 2.495 2.496 NS -1.212 -1.211
Table 7.8: Velocity [m/s] with zero mixture slip boundary condition
also being zero. This was obtained by setting the slip velocity to zero in the method
described in Section 6.1.1. Results for this boundary condition are shown in Table 7.9.
The mixture velocity is the same as the results from the zero mixture slip boundary, as
would be expected. The species velocities are slightly different from the zero mixture slip
case, the magnitude of this difference depending on the problem being considered. As with
the previous calculations the rate of convergence of the mixture velocity is second order,
while for the species velocity it is only first order. This means the source of the error is
not due to the expression for the diffusion slip velocity at the wall.
Mixture Velocity Hydrogen Velocity Water velocity
c
Nx Nx Nx
30 60 90 30 60 90 30 60 90
2000 -0.837 -0.835 -0.834 4.410 3.425 3.094 -1.435 -1.317 -1.278
8000 -0.837 -0.835 -0.834 2.914 2.674 2.593 -1.264 -1.232 -1.222
32000 -0.837 -0.835 -0.834 2.540 2.486 2.467 -1.222 -1.211 -1.207
128000 -0.837 -0.835 -0.834 2.447 2.439 2.436 -1.211 -1.206 -1.204
256000 NS -0.835 -0.834 NS 2.431 2.431 NS -1.205 -1.203
384000 NS -0.835 -0.834 NS 2.428 2.429 NS -1.205 -1.203
Table 7.9: Velocity [m/s] with no species slip boundary condition
108
7.10 Conclusions
Three different equilibrium functions were considered for multi-component gases. The
Mixture Centred Equilibrium (MCE) function had the widest stability range, and was
selected as the preferred model. The Navier-Stokes equations for the mixture are recovered
in terms of a mixture viscosity. The Stefan-Maxwell equation and the ‘pressure diffusion’
equation were recovered for the species velocities. An expression was developed for diffusive
slip boundaries, using the moment-based boundary method.
Although the MCE scheme with diffusive boundaries gives results consistent with the
work of Young and Todd (2005) for binary channel flow, the convergence of the species
velocity is only first order. This does not make the scheme a feasible tool for real applica-
tions. In fact, the source of this error was the assumption that Π
(eq)+
X = Π
(eq)
X made before
Equation (7.60). Because the species momentum is not conserved during the collision this
assumption is incorrect. The error term is proportional to the size of the time step, and
so as the time step reduces the error reduces also. In the next chapter a second order
numerical integration scheme will be used which eliminates this error. These results show
that when the LB method is applied to a problem where the momentum is not conserved
over a time step, the simple forward Euler integration is not sufficient to ensure second
order accuracy.
109
Chapter 8
Improvements to Multi-component
Methods
In Chapter 7 a continuous LB model was derived which recovered the steady-state macro-
scopic equations for a multi-component gas mixture. However, when this scheme was
numerically integrated using the forward Euler scheme, it was found to be only first order
accurate. In Chapter 6 the single-component scheme developed was second order accurate.
However, for large time steps the scheme was unstable meaning that the computational
time could not be reduced at the expense of spurious accuracy. In this chapter both these
issues are addressed. Firstly, a second order accurate scheme is used to ensure the species
velocity is second order convergent. Secondly, a pre-conditioning method is used to enable
calculations to be performed with much larger time steps.
In this chapter all variables are in the computation scale, and for clarity the˜is omitted.
8.1 Second order scheme
8.1.1 Crank-Nicolson integration
To avoid numerical errors, and obtain second order accuracy for the species velocity, He
et al. (1998), Asinari (2006) and Dellar (2002b) propose the use of the Crank-Nicolson rule
110
to replace Equation (5.21),
f+σ = fσ +
1
2
A(f (eq)σ − fσ) +
1
2
A(+)(f (eq)+σ − f+σ ) (8.1)
where the + superscript is the post collision and streaming state, and A = M−1ΛM .
Equation (8.1) is an implicit expression for f+σ , and threatens to ruin the explicit nature
of the LBM. Fortunately a transformation to a new variable gσ can be used,
gσ = fσ − 1
2
A(f (eq)σ − fσ) (8.2)
Substituting Equation (8.2) into Equation (8.1),
g+σ = gσ −A(I +
1
2
A)−1(f (eq)σ − gσ) (8.3)
where I is the identity matrix. This equation may be used directly, or alternatively in the
form,
g+σ = gσ −A
′
(f (eq)σ − gσ) (8.4)
where A
′
is a modified form of A, such that each value of λ is replaced with λ/(1 + λ/2).
Therefore, in effect, the algorithm used so far is unchanged, except that the vector fσ is
replaced by gσ. The following identities arise directly from Equation (8.2) when the MCE
model is used, for simplicity λξ = λν , so the matrix Λ is diagonal.
∑
i
gσ,i =
∑
i
fσ,i = ρσ (8.5)
The sum over i of gσ,i is equal to the sum over i of fσ,i as the species density is conserved
during the collision. However momentum is not, and a new velocity u¨σ,x, defined as ρσu¨σ,α =∑
i eα,igσ,i, is introduced. Using Equation (8.2) yields,
ρσu¨σ,α = ρσuσ,α − λ
δ
2
ρσ(uα − uσ,α) (8.6)
111
Summing over all the species, momentum is conserved, and
∑
σ ρσu¨σ,α = ρuα. Using Φσ,αβ
to represent the equivalent of the momentum flux density tensor calculated from gσ (i.e.∑
i eα,ieβ,igσ,i = Φσ,αβ),
Φσ,αβ = Πσ,αβ − λ
ν
2
(Π
(eq)
σ,αβ − Πσ,αβ) (8.7)
Only when Πσ,αβ = Π
(eq)
σ,αβ does Φσ,αβ = Πσ,αβ.
Once the calculation is complete ρσ can be calculated directly from the values of gσ,
but uσ must be calculated from u¨σ using Equation (8.6). Aside from this calculation,
which only needs to be performed once, there is no additional computational demand from
the variable substitution method. However, the boundary conditions must be rewritten in
terms of gσ, an exercise which is conceptually simple yet algebraically overwhelming.
8.1.2 Diffusion slip boundary
Unfortunately, the use of the variable transform method adds another layer of complexity
to the treatment of boundaries. In general (when Λ cannot be expressed as a scalar
multiplying the identity matrix, as in the single relaxation time model), Equation (8.2)
cannot be used to convert from gσ to fσ if there are unknown values of fσ. Therefore, the
boundary conditions must be reformulated in terms of gσ, instead of fσ.
For the wall boundary condition the approach of Section 7.6.1 is repeated in terms of
gσ. Equation (8.6) becomes, using ρu = ρAu¨A + ρBu¨B,
uA = u¨A + (
λδ
2
1 + λ
δ
2
)
ρB
ρA + ρB
(u¨B − u¨A) (8.8)
uB = u¨B + (
λδ
2
1 + λ
δ
2
)
ρA
ρA + ρB
(u¨A − u¨B) (8.9)
This may be substituted into the macroscopic boundary conditions from Section 7.6 ,∑
σ(ρσuσ,x/S(Mσ)) = 0, to obtain a condition for u¨σ,
u¨A = −
ρB
(
λδ (ρAS(MB) + ρBS(MA)) + 2(ρA + ρB)S(MA)
)
ρA (λδ (ρAS(MB) + ρBS(MA)) + 2(ρA + ρB)S(MB))
u¨B (8.10)
112
Equation (8.7) yields an expression for the XY momentum flux,
ΦXY = ΠXY − 1
2
λν(Π
(eq)
XY − ΠXY ) (8.11)
Which, when combined with the macroscopic boundary conditions from Section 7.6, ΠA,XY /ρA =
ΠB,XY /ρB, becomes
1
ρA
(ΦA,XY − 1
2
λνΠ
(eq)
A,XY ) =
1
ρB
(ΠB,XY − 1
2
λνΠ
(eq)
B,XY ) (8.12)
If either the SCE or CSCE model were used it would lead to yet more complexity, as uσ
appears in the expression for Πσ,XY . However in the MCE model Π
(eq)
σ,XY = ρσuxuy, leading
to a cancellation of the Π
(eq)
σ,XY terms and Equation (8.12) becomes,
1
ρA
(ΦA,XY ) =
1
ρB
(ΦB,XY ) (8.13)
which may be written in terms of the elements of gA and gB,
1
ρA
(gA,5 − gA,6 + gA,7 − gA,8) = 1
ρB
(gB,5 − gB,6 + gB,7 − gB,8) (8.14)
The u¨ values are defined by,
gσ,5 − gσ,6 = ρσu¨σ,x − gσ,1 − gσ,8 + gσ,3 + gσ,7 (8.15)
and substituting this into Equation (8.14) yields,
1
ρA
(ρAu¨A,x− gA,1− 2gA,8 + gA,3 + 2gA,7) = 1
ρB
(ρBu¨B,x− gB,1− 2gB,8 + gB,3 + 2gB,7) (8.16)
A new variable is introduced, κ¨σ = −gσ,1 − 2gσ,8 + gσ,3 + 2gσ,7 equivalent to the κσ defined
in Section 7.6 in terms of f . Using Equation (8.10) to eliminate u¨B,x yields an expression
113
for u¨A,x purely in terms of known values,
ρAu¨A,x
(
1
ρA
+
1
ρB
(λ (ρAS(MB) + ρBS(MA)) + 2(ρA + ρB)S(MB))
(λ (ρAS(MB) + ρBS(MA)) + 2(ρA + ρB)S(MA))
)
=
1
ρB
(κB)− 1
ρA
(κA)
(8.17)
To impose this value of u¨A,x on the wall, the method discussed in Section 6.1.1 is used.
Noting that if uσ,y = 0, u¨σ,y = 0 and if Πσ,XX = Π
(eq)
σ,XX then Φσ,XX = Π
(eq)
σ,XX , the references
to fσ in Equations (6.23) may be swapped for gσ.
8.1.3 Other boundaries
At an open boundary ρσ is set to a specific value, uσ,y = 0 and Πσ,Y Y = Π
(eq)
σ,Y Y , similar to the
conditions used in Chapter 6. In terms of g these conditions are the same, specifically u¨σ,y =
0 and Φσ,Y Y = Π
(eq)
σ,Y Y . Therefore the same conditions already developed in Section 6.1.1
may be used, substituting gσ for fσ in Equation (6.24).
It may be desirable to set the mass flux, instead of the density, as a boundary condition.
For given values of ρσu¨σ,x, the value of ρσ is not uniquely defined. If this value of ρσ can be
determined it can be used in conjunction with the method from the previous paragraph.
In a corner of the type shown in Figure 5.2, ρσu¨σ,x and ρσ are linked by the equations,
ρσ(1 + u¨σ,x) = gσ,0 + gσ,2 + gσ,4 + 2 (gσ,1 + gσ,5 + gσ,8) = Υσ (8.18)
Where Υσ is a shorthand involving only the known values of gσ. The values of u¨σ are given
by,
ρAu¨A,x = ρAuA,x − λ
δ
2
ρA(ux − uA,x) (8.19)
ρBu¨B,x = ρAuB,x − λ
δ
2
ρB(ux − uB,x) (8.20)
114
Substituting Equation (8.18) into Equations (8.19) and (8.20) yields,
ΥA − ρA = ρAuA,x − λ
δ
2
ρA(ux − uA,x) (8.21)
ΥB − ρB = ρBuB,x − λ
δ
2
ρB(ux − uB,x) (8.22)
which can be summed over the species to give
ΥA + ΥB − ρ = ρux (8.23)
If Υσ and ρux are known, the value of ρ can be determined, which allows the known value of
ρux to be split and ux to be recovered. The only unknowns remaining in Equations (8.21)
and (8.22) are the species densities. Rearranging Equation (8.23) gives,
ρA =
ΥA − ρAuA,x(1 + λδ2 )
1− λδ
2
ux
(8.24)
and a similar expression can be obtained for ρB. In fact, the situation is more complicated
than it appears, as λδ is itself a function of the species densities through Equation (7.25).
To circumvent this problem the species density from the same grid point at the previous
time step is used. Although a proper solution could be found this would add additional
complexity, and seems to give no corresponding benefit.
8.1.4 Intersection between a wall and another boundary
Open boundaries with known densities
Following the approach of Section 7.6.2, but in terms of gσ, the following conditions are
true for the gσ values at a corner node at the intersection of an open boundary and a wall,
115
if the density at the open boundary is known.
gσ,1 + gσ,2 + gσ,5 + gσ,6 + gσ,8 = ρσ − gσ,3 − gσ,4 − gσ,7 − gσ,0 (8.25)
gσ,1 + gσ,5 + gσ,8 − gσ,6 = ρσu¨σ,x + gσ,3 + gσ,7 (8.26)
gσ,2 + gσ,6 + gσ,5 − gσ,8 = gσ,7 + gσ,4 (8.27)
Introducing the new variable Θ¨σ = 4gσ,7 + gσ,0 + 2gσ,3 + 2gσ,4, then the equivalent of
Equation (7.55) is,
(Θ¨A − ρA + ρAu¨A,x) 1
ρA
= (Θ¨B − ρB + ρBu¨B,x) 1
ρB
(8.28)
Substituting Equation (8.10) leaves an expression for the value of u¨A,x to be imposed at
the corner node,
u¨A,x
(
1+
ρA
(
λδ (ρAS(MB) + ρBS(MA)) + 2(ρA+ρB)S(MB)
)
ρB (λδ (ρAS(MB) + ρBS(MA)) + 2(ρA+ρB)S(MA))
)
=
1
ρB
(Θ¨B−ρB)− 1
ρA
(Θ¨A−ρA)
(8.29)
To impose this, the approach from Section 6.1.1, with fσ,i values swapped for gσ,i values,
is used.
Open boundaries with known mass flux
If the open boundary that intersects the corner is specified by a mass flux rather than a den-
sity condition, the diffusion slip condition for these nodes cannot be applied. Therfore, the
values of ρσuσ,x are known, and uσ,y = 0, Πσ,Y Y = ρσ(RT )σ and Πσ,XX = ρσ(RT )σ + ρσu
2
x.
As there are five unknown g values there must be five constraints, the final condition being
Πσ,XY = 0. Certainly, along the wall ∂ρσuy/∂x = 0, however along the open boundary
∂ρσux/∂y is only zero if the mass flux is constant. If a non-constant mass flux profile is to
be used, a specific value of Πσ,XY = 0 must be supplied. At the corner, the density and
116
modified velocity are related by,
ρσ + ρσu¨σ,x = Θ¨σ (8.30)
Following the approach of Section 8.1.3 the density can be found as,
ρA =
Θ¨A − ρAuA,x(1 + λδ2 )
1− λδ
2
ux
(8.31)
For example the gσ values to be imposed at a corner with gσ,3 gσ,4 and gσ,7 known are,
gσ,1 = −gσ,0 − gσ,3 + (1− (RT )σ)ρσ (8.32)
gσ,2 = −gσ,0 − gσ,4 − ((RT )σ + u2x − 1)ρσ (8.33)
gσ,5 =
1
2
gσ,0 − gσ,7 +
(
(RT )σ +
1
2
(u2x − 1)
)
ρσ (8.34)
gσ,6 =
1
2
gσ,0 + gσ,4 + gσ,7 +
(
1
2
(RT )σ + u
2
x − 1
)
ρσ (8.35)
gσ,8 = −gσ,4 − gσ,7 + 1
2
(RT )σρσ (8.36)
8.1.5 Binary channel flow results
The variable substitution method described above was used to solve the problem from
Section 7.9. The impressive results are shown in Table 8.1. Once the value of c is high
enough such that R˜T < 1/3, Ma is so small that further reduction does not affect the
results. This behaviour is similar to that obtained from the dimensional single-component
Poiseuille flow. The use of the variable substitution method has led to a dramatic increase
in accuracy, leading to lower computational times.
8.2 Pre-conditioning
The method that has been developed so far is not time accurate, unless ˜(RT )σ = 1/3 which,
for species of different molecular weights, is impossible. Although errors are introduced
117
Mixture Velocity Hydrogen Velocity Water velocity
c
Nx Nx Nx
30 60 90 30 60 90 30 60 90
1300 -1.720 -1.385 -1.384 1.130 1.944 1.946 -2.051 -1.761 -1.759
2000 -1.388 -1.385 -1.384 1.939 1.944 1.946 -1.766 -1.761 -1.759
3000 -1.388 -1.385 -1.384 1.939 1.944 1.946 -1.766 -1.761 -1.759
4000 -1.388 -1.385 -1.384 1.939 1.944 1.946 -1.766 -1.761 -1.759
Table 8.1: Velocity [m/s] with second order boundary integration. Results using the method of
Young and Todd (2005) are Mixture velocity = -1.3935 m/s, Hydrogen Velocity = 1.9404 m/s
and Water Vapour velocity = -1.7672 m/s
into the transient solution, the computation will proceed at a pace determined by the
physics of the problem. It is desirable to completely sever this link with physical time,
and thus produce the quickest possible convergence to the steady state. Verberg and
Ladd (1999) suggest one method of doing this, another is a technique known as pre-
conditioning, suggested for LBM by Guo et al. (2004) and further refined by Izquierdo and
Fueyo (2009). The approach outlined here is a simplification of pre-conditioning, yet it is
effective. To accomplish this pre-conditioning all the terms except the time derivative term
in the species momentum equation for the MCE model, Equation (7.22), are multiplied by
the non-dimensional quantity γ. Thus,
∂
∂tˆ
 ρˆσuˆσ,x
ρˆσuˆσ,y
+ γ ∂
∂xˆ
 ρˆσuˆ2x
ρˆσuˆxuˆy
+ γ ∂
∂yˆ
 ρˆσuˆxuˆy
ρˆσuˆ
2
y
 = −γ
 ∂∂xˆ(pˆσ)
∂
∂yˆ
(pˆσ)

+ γ
∂
∂xˆ
 Π(neq)σ,XX
Π
(neq)
σ,XY
+ γ ∂
∂yˆ
 Π(neq)σ,XY
Π
(neq)
σ,Y Y
+ γ
 λˆδρˆσ(uˆx − uˆσ,x)
λˆδρˆσ(uˆy − uˆσ,y)
 (8.37)
In steady state conditions the time derivative term is zero and the correct momentum
equation is recovered. To implement pre-conditioning, modifications must be made to the
equilibrium function, (RT ) and ν. The modified quantities are indicated with a super-
118
c
γ
0.015 0.030 0.06 0.12 0.24
250 4.72×105 2.39×105 NS NS NS
500 9.06×106 4.65×105 2.38×105 1.21×105 NS
750 1.31×106 6.76×105 3.48×105 1.79×105 9.09E×104
1000 1.70×106 8.76×105 4.54×105 2.36×105 1.20E×105
Table 8.2: Number of iterations needed to reach converged solution to the problem from Section
7.9 with different values of c and γ. The accuracy of all solutions were within 0.01% of each other.
The diffusive slip boundary condition was used, and the variable transform method of Section
8.1. NS means no solution could be found.
script γ.
Πeq,γσ,XX = pσ + γρσu
2
x (8.38)
Πeq,γσ,Y Y = pσ + γρσu
2
y (8.39)
Πeq,γσ,XY = γρσuxuy (8.40)
(RT )γ = γ(RT ) (8.41)
νγ = γν (8.42)
Both Guo et al. (2004) and Izquierdo and Fueyo (2009) suggest modifications to the ΠXXY Y
moment. However, as this moment does not influence the momentum equation such mod-
ifications seem unjustified.
The results shown in Table 8.2 show that the pre-conditioning can have a major impact
on the convergence time of the results without reducing accuracy. Higher values of γ with
a constant value of c produce quicker convergence. To minimise the convergence time a
large value of c and the largest possible value of γ should be selected. However, it is not
clear why this is the case.
119
8.3 Conclusions
In conclusion both the variable transformation method and the pre-conditioning method
are highly successful in reducing the computational expense of the multi-component LB
model. As such they should both be considered vital components of the model.
Using the variable transformation method requires careful handling of the boundary
conditions. It is interesting to note that a bounce-back boundary condition is often used
on the values of the new variable g, and this is not necessarily the same condition as
bounce-back of the f values. The moment-based boundary conditions for the new variable
are no more computationally expensive than the conditions used on the old variables and,
crucially they are still local to each node. Implementing pre-conditioning is extremely
simple, requiring only minimal changes.
120
Chapter 9
Calculations of Binary Flow
In this chapter some calculations are performed using the method developed in the previous
chapters. These calculations will be loosely themed around fuel cell development and, as
such, will involve hydrogen and water vapour diffusion. Some types of hydrogen fuel cell
such as the SOFC consist of a reacting surface where oxygen ions, conducted through
this surface, react with the hydrogen in the gas mixture to produce water. The hydrogen
arrives at the surface by diffusing through a porous material and the water leaves through
the same porous material. Each mole of hydrogen that is absorbed will produce one mole
of H2O, therefore at the surface pH2OuH2O = −pH2uH2 , which may be written as,
ρH2uH2,x
ρH2OuH2O,x
= − MH2
MH2O
(9.1)
where MH2 and MH2O are the molar masses of hydrogen and water respectively. The quan-
tity ρH2uH2,x/ρH2OuH2O,x will be referred to as the species mass flux ratio, and −MH2/MH2O
is its stoichiometric value. Due to mass conservation of each species, the species mass flux
ratio must be maintained on average throughout the gas diffusion path. It was shown in
Section 7.6 that the diffusion slip boundary forced a different value of mass flux ratio, and
so the bulk of the flow must compensate to maintain the average value at the stoichiometric
ratio. For the rest of this section the species mass flux ratio is assumed to be averaged
across the channel, and is indicated by an overbar such as ρH2uH2,x/ρH2OuH2O,x.
121
The water is only considered in the gas phase, because SOFCs operate at high temper-
atures and pressures. Other types of fuel cell encounter condensation of water vapour, and
this effect is beyond the scope of the present work.
9.1 Fully developed binary channel flow
∆
p
[P
a]
XH2 at x = L
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.11
0.2
0.4
0.6
0.8
1
0.2 0.4 0.6 0.8
500
1000
1500
2000
2500
3000
Figure 9.1: Contours of the average species mass flux ratio (ρH2ux,H2/ρH2Oux,H2O) for large
values of ∆p. For large values of ∆p the species mass flux ratio is a function of concentration
gradient only. The region between ∆p = 200 and ∆p = −200 is not plotted here, and details of
this region may be found in Figure 9.2.
When the partial pressures of the species are set at each end of a channel, using the
boundary conditions from Section 7.9, the solution obtained is the binary equivalent of
fully-developed single-component flow. That is, the viscous boundary layer at the wall
has formed and the velocity profile is parabolic (or, with no total pressure gradient, flat).
There is no equivalent of the viscous boundary layer for diffusion.
122
∆
p
[P
a]
XH2 at x = L
-0.33
-0.33
-0.22
-0.22
-0.11
-0.11
-0.055
-0.055
0
0
0.055
0.055
0.11
0.11
0.22
0.22
0.33
0.33
0.2 0.4 0.6 0.8
−100
−80
−60
−40
−20
0
20
40
60
80
100
Figure 9.2: Contours of the average species mass flux ratio (ρH2ux,H2/ρH2Oux,H2O) for small
values of ∆p. The dashed line represents an infinite flux ratio , i.e. ρH2Oux,H2O = 0. At the
central point, when ∆p = 0 and XH2 = 0.5, there is a singularity in the flux ratio as both species
fluxes tend to zero. Note
√
MH2/MH2O = 0.3333 and MH2/MH2O = 0.1111.
123
A rectangular domain was considered, with a wall at y = 0 and y = H. The boundary
conditions at x = 0 are fixed, with p = 1 bar and XH2 = 0.5 unless otherwise stated. At
x = L different values of total pressure and XH2 were applied. The change in total pressure
is defined as ∆p = px=L − px=0. The physical size of the domain is 10 µm × 10 µm unless
otherwise stated, a grid size of Nx = 30 and a grid speed of c = 6000 m/s was used.
∆
p
[P
a]
XH2 at x = L
-1
-0.8
-0.6
-0.4
-0.2
00.20.4
0.2 0.4 0.6 0.8
−100
−80
−60
−40
−20
0
20
40
60
80
100
Figure 9.3: Contours of the average mass flux of hydrogen, ρH2ux,H2 [kg m
−2 s−1], at different
pressure and concentration gradients. Along the dashed line the mass flux ratio is stoichiometric.
If the magnitude of the change in total pressure is very large both species are carried
in the same direction by the pressure gradient from high pressure to low pressure, and
the species mass flux ratio is positive, as shown in Figure 9.1. However, for much smaller
values of pressure drop the species concentration gradient begins to play a more important
role, and it is possible to obtain negative values for the species mass flux ratio.
Figure 9.2 shows contours of the average species mass flux ratio for different concentra-
tion and pressure gradients. If there is zero total pressure gradient the species velocity is
governed by the diffusion slip wall boundary condition, and the species mass flux ratio is
124
∆
p
[P
a]
XH2 at x = L
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
0.2 0.4 0.6 0.8
−100
−80
−60
−40
−20
0
20
40
60
80
100
Figure 9.4: Contours of the average mass flux of water, ρH2Oux,H2O [kg m
−2 s−1], at different
pressure and concentration gradients. Along the dashed line the mass flux ratio is stoichiometric.
125
−√MH2/MH2O. If, instead, there is no concentration gradient then the species velocities
are equal and the value of the mass flux ratio is XH2MH2/XH2OMH2O. If there is neither
a species nor a concentration gradient then all velocities are zero and the species mass
flux ratio is undefined. Aside from these simple cases the lines of constant species mass
flux ratio form a radial pattern. When the hydrogen velocity is zero, as the convective
and diffusive fluxes cancel, the species mass flux ratio is zero and likewise when the water
velocity is zero the species mass flux ratio is infinite. Figure 9.3 and 9.4 show the mass
 
 
H = 30µm
H = 25µm
H = 20µm
H = 15µm
H = 10µm
H = 5µm
∆
p
[P
a]
XH2 at x = L
0.1 0.2 0.3 0.4
0
10
20
30
40
50
60
70
Figure 9.5: Pressure difference against XH2 at the outlet for different channel widths. The results
shown correspond to the mass flux ratio being stoichiometric. The pressure drop corresponding
to a specific value of XH2 depends on the channel width.
flux of hydrogen and water corresponding to Figure 9.2. For a positive flux of hydrogen
it is not sufficient to ensure the concentration gradient is negative, one must also consider
the pressure gradient. Likewise for a large enough total pressure drop it is possible for the
mass flux of hydrogen to be positive despite the concentration gradient also being positive.
However, if the species mass flux ratio is stoichiometric, a negative concentration gradient
is required for a positive hydrogen mass flux.
Figure 9.5 shows the relationship between total pressure change and XH2 at x = L for
126
  
H = 30µm
H = 25µm
H = 20µm
H = 15µm
H = 10µm
H = 5µm
ρ
u
H
2
[k
g/
m
2
s]
XH2 at x = L
0 0.1 0.2 0.3 0.4 0.5
0
0.005
0.01
0.015
Figure 9.6: Average hydrogen mass flux against XH2 at the outlet for different channel widths.
The results shown correspond to the mass flux ratio being stoichiometric.
different channel widths, when the species mass flux ratio is stoichiometric. The length of
the channel is constant at 5 µm. Figure 9.6 shows the average mass flux of hydrogen which
is nearly independent of channel width. If the hydrogen mole fractions at each end of the
channel are maintained constant, increasing the channel width increases the pressure drop
along the channel but does not affect the average mass flux though the channel. This is
because the hydrogen is driven by diffusive rather than viscous forces.
If the mole fraction of hydrogen at x = 0 is 0.25, instead of 0.5, the species mass flux
ratio for different pressure and concentration gradients are as shown in Figure 9.7. The
features discussed previously are still observable, with the flux ratio at zero concentration
gradient still being XH2MH2/XH2OMH2O.
Finally, a similar plot is obtained for the original problem (XH2 = 0.5 at x = 0) with
zero mixture slip at the walls. Figure 9.8 shows the results obtained. Note that for zero
pressure drop the average mass flux ratio is now -1. For zero concentration gradient the
flux ratio is the same as that obtained with the diffusive slip boundary condition as seen in
Figure 9.2. Comparison between these two figures illustrates the importance of diffusion
127
∆
p
[P
a]
XH2 at x = L
-0.33
-0.22
-0.22
-0.11
-0.11
-0.0371
-0.0371
0
0
0.0371
0.0371
0.11
0.11
0.22
0.22
0.33
0.33
0.2 0.4 0.6 0.8
−100
−80
−60
−40
−20
0
20
40
60
80
100
Figure 9.7: Contours of the species mass flux ratio for small values of ∆p, similar to Figure 9.2,
except that the mole fraction of hydrogen at x = 0 is 0.25. The dashed line represents an infinite
flux ratio.
128
∆
p
[P
a]
XH2 at x = L
-1 -1
-0.33
-0.33
-0.22
-0.22
-0.11
-0.11
-0.055
-0.055
0
0
0.055
0.055
0.11
0.22
0.22
0.33
0.33
1
1
0.2 0.4 0.6 0.8
−100
−80
−60
−40
−20
0
20
40
60
80
100
Figure 9.8: Contours of the species mass flux ratio for small values of ∆p, similar to Figure 9.2,
except that the wall boundary condition is zero mixture velocity. The dashed line represents an
infinite flux ratio.
129
slip for problems dominated by diffusion. Conversely, for small concentration gradients the
inclusion of diffusion slip at the walls has very little effect on the species mass flux ratio.
9.2 Binary flow with one inert component
9.2.1 The Stefan tube problem
The Stefan tube is a well known problem in the theory of gas diffusion. A tube is exposed
to air at one end, and water (as a liquid) at the other end. The water evaporates from the
liquid surface and diffuses through the tube into the air, while other gases in the air, mostly
nitrogen, cannot pass through the end of the tube blocked by the water. It is common to
reduce this to a binary problem, treating air as an inert ideal gas or sometimes as pure
nitrogen. This helps to reduce the complexity of the analysis.
The Stefan tube is traditionally treated as a one-dimensional problem, see, for example,
Eckert and Drake (1959) and Bird et al. (1960) who built their analysis on the 1874
paper by Stefan. The analysis hinges upon the assumption that the inert gas is stationary
everywhere, because the average mass flow rate is zero. After this assumption is made, the
analysis is straightforward and the result is that the mass flow rate of the water vapour,
labelled A, is given by,
ρAuA = − DAB
(RT )A
p
p− pA
dpA
dx
(9.2)
where x is the direction in which the diffusion is taking place. Equation (9.2) is known as
Stefan’s Law. The species partial pressures do not vary with y, subsequently the velocity of
gas A is a function of x only, and the solution is one-dimensional. In fact, this approach is
unsatisfactory because of the velocities it predicts at the wall, which agree with neither the
zero mixture slip approximation nor the diffusion slip theory (nor even a zero species slip
condition). Using the no-slip wall condition a two-dimensional solution, with circulation
of the inert gas, was predicted by Meyer and Kostin (1975).
Salcedo-Diaz et al. (2008) used a model that includes diffusion slip to solve the Stefan
tube problem. They found that the velocity of the inert gas was not zero, but at the wall it
130
flowed towards the end saturated with water, and in mid-channel it flowed in the opposite
direction. On average the mass flux of the inert gas was zero, as required by the principle
of conservation of mass. The flux of water vapour was close to that predicted by Equation
(9.2), for a wide channel, but for a narrow channel it was not. This was also predicted
by Whitaker (1991). Interestingly, it appears from the results of Salcedo-Diaz that in a
narrow channel the flux of the diffusing gas is almost independent of the mole fraction at
the bottom of the tube. However, the Salcedo-Diaz model includes non-continuum effects,
and it is likely that this result is due to these. The lack of agreement with the continuum
based Stefan’s Law is then hardly surprising.
The calculations in the subsequent section are not of the Stefan tube problem, but a
similar type of problem; binary flow through a channel with a reacting surface at one end.
The two gases will be labelled the inert and reacting gas, and will be water vapour and
hydrogen respectively. Therefore the water vapour is no longer the diffusing gas, and now
becomes the inert gas.
9.2.2 A LB model of a Stefan tube type problem
Instead of studying the classical Stefan tube problem, in this section a LB model of a
similar problem is investigated. This problem may be thought of as diffusion in a solid
oxide fuel cell with particularly low chemical activity, such that the products of the reaction
may be ignored. Two components, hydrogen and water vapour, are present in a channel
closed at one end by a reacting surface. The water vapour does not react at this surface
(at x = L), while the hydrogen is absorbed at a specific rate. At the other end of the
channel the gas has a known composition and total pressure.
To model this problem, the boundary conditions described in Section 8.1.3 are used. At
the boundary where x = 0 the conditions of 1 bar total pressure and XH2 = 0.5 are set. A
summary of this and the other boundary conditions is shown in Figure 9.9.
Unfortunately, the results obtained are not perfect. To illustrate this, an example
calculation was performed with the mass flux of hydrogen at x = L set to 0.1 kg m−2 s−1.
131
xy
Re
ac
tin
g s
urf
ac
e
Figure 9.9: The boundary conditions at the wall, the open boundary, the reacting surface and
the intersections between them for the Stefan tube type problem of Section 9.2.2.
The tube is 20 µm long and 10 µm wide and a 60×30 grid was used. Figures 9.10 and
9.11 show the x and y velocities of the species and the mixture. Figure 9.10 shows that
the hydrogen has a positive x velocity throughout the domain, which is larger along the
centreline than at the walls. The water vapour has a negative x velocity along the wall and
a positive x velocity near the centre of the channel. The mixture x velocity has a similar
pattern to that of the water, but the total mass flux in the x direction of the mixture is
not zero. In Figure 9.11 the y velocities of the species and mixture show that significant
circulation is taking place near the reacting surface at x = 60.
Figure 9.12 shows the total pressure of the mixture. Clearly, the corner nodes at the
intersection between the reacting surface and the wall have unrealistic total pressures. This
is because of the discontinuity between the mixture velocity along the wall, the diffusion
slip velocity, and the set mass flux imposed at the corners. Therefore, ∂ux/∂x is not zero.
The boundary conditions include the condition that ΠXX = Π
eq
XX , equivalently Π
neq
XX = 0
and ∂ux/∂x = 0. There is thus an inconsistency between the boundary conditions, not an
error with the implementation.
The error in the total pressure has a limited impact away from the affected corners.
Figure 9.13 shows the cross-channel-averaged total pressure along the channel. The effect
of the misbehaving corner nodes can be seen but only in the last ten grid points. Figure
132
  
 
 
 
 
Mixture
y
x
Water
y
x
Hydrogen
y
x
20 40 60
20 40 60
20 40 60
−0.2
0
0.2
0.4
10
20
30
−0.6
−0.4
−0.2
0
0.2
10
20
30
2
2.5
3
10
20
30
Figure 9.10: Contours of the x component of the species and mixture velocities [m/s] for the
Stefan tube type problem. The x velocity of the water is zero at x = 60.
133
  
 
 
 
 
Mixture
y
x
Water
y
x
Hydrogen
y
x
20 40 60
20 40 60
20 40 60
−0.2
0
0.2
10
20
30
−0.2
0
0.2
10
20
30
−0.2
0
0.2
10
20
30
Figure 9.11: Contours of the y component of the species and mixture velocities [m/s] for the
Stefan tube type problem. The y velocity of both species is zero at all the boundaries.
134
y x
T
ot
al
P
re
ss
u
re
(P
a)
m
in
u
s
1
b
ar
20
40
60
10
20
30
−150
−100
−50
0
50
Figure 9.12: Total pressure minus 1 bar, showing the unrealistic effect at the corners. The reacting
surface is located at x = 60.
C
ro
ss
-c
h
an
n
el
av
er
ag
ed
p
(P
a)
m
in
u
s
1
b
ar
x
20 40 60
−25
−20
−15
−10
−5
0
Figure 9.13: Cross-channel-averaged total pressure minus 1 bar at different distances along the
channel. The disturbance in the pressure field caused by the corners is only evident near the
reacting surface.
135
9.14 shows the mole fractions of the two gases along the centreline of the channel, and no
obvious influence of the corner nodes can be observed.
M
o
le
fr
a
c
t
io
n
x
20 40 60
0.35
0.4
0.45
0.5
0.55
0.6
0.65
0.7
Figure 9.14: Mole fraction of hydrogen (dashed line) and water vapour (solid line) recorded along
the centreline of the channel.
9.2.3 The effect of channel width
In this section, some results are presented for the problem described in Section 9.2.2 at
different channel widths and lengths. To alleviate the strange behaviour of the total pres-
sure at the corners where x = L, the pressures at x = L were linearly extrapolated from
the other half of the domain.
Figure 9.15 shows the relationship between the hydrogen pressure difference along the
channel and the average hydrogen mass flux for a channel 10 µm long and of varying
widths. Also plotted is the Stefan’s Law result, from Equation (9.2), which is independent
of channel width. Figure 9.16 shows the change in hydrogen mole fraction along the
channel. Figures 9.17 and 9.18 show similar graphs for a channel 20 µm long. Clearly, the
one-dimensional approximation of Stefan’s Law does not lead to a realistic solution.
These results may be explained by considering the domain as two regions; the region
136
  
Stefans Law
Width = 20µm
Width = 10µm
Width = 5µm
H
y
d
ro
ge
n
p
re
ss
u
re
d
iff
er
en
ce
(p
H
2
,x
=
0
−
p
H
2
,x
=
L
)
[P
a]
Mass flux of hydrogen [kg/m2 s]
0 0.1 0.2 0.3 0.4
102
103
104
105
Figure 9.15: The change in hydrogen pressure drop along the channel for the Stefan tube type
problem plotted against the average mass flux of hydrogen for different channel widths. Also
plotted is the Stefan’s Law one-dimensional solution.
137
  
Width = 20µm
Width = 10µm
Width = 5µm
H
y
d
ro
ge
n
m
ol
e
fr
ac
ti
on
d
iff
er
en
ce
(X
H
2
,x
=
0
−
X
H
2
,x
=
L
)
Mass flux of hydrogen [kg/m2 s]
0 0.1 0.2 0.3 0.4
10−3
10−2
10−1
100
Figure 9.16: The change in hydrogen mole fraction along the Stefan tube problem plotted against
the average mass flux of hydrogen at different channel widths.
138
close to the wall and the region in the middle of the channel. Close to the wall the ef-
fect of diffusion slip dominates so that, changes in the width of the channel have only a
small effect on the hydrogen mole fraction gradient and the diffusion slip velocity (through
Equation (7.40)). The mixture velocity is negative and acts to reduce the average mass
flux of hydrogen. In the region near the centre of the channel the total pressure gradient
dominates, causing the mixture velocity to be positive and the average mass flux of hy-
drogen to increase. Increasing the width of the channel increases the width of the central
region, and to maintain the same average hydrogen mass flux the total pressure gradient
must be reduced.
 
 
Stefans Law
Width = 40µm
Width = 20µm
Width = 10µm
H
y
d
ro
ge
n
p
re
ss
u
re
d
iff
er
en
ce
(p
H
2
,x
=
0
−
p
H
2
,x
=
L
)
[P
a]
Mass flux of hydrogen [kg/m2 s]
0 0.1 0.2 0.3 0.4
102
103
104
105
Figure 9.17: As Figure 9.15 but the channel length is 20µm
9.3 Conclusions
The results of the two problems in this chapter are very relevant in the context of solid
oxide fuel cells, although, on their own, they do not constitute a thorough examination of
139
  
Width = 40µm
Width = 20µm
Width = 10µm
H
y
d
ro
ge
n
m
ol
e
fr
ac
ti
on
d
iff
er
en
ce
(X
H
2
,x
=
0
−
X
H
2
,x
=
L
)
Mass flux of hydrogen [kg/m2 s]
0 0.1 0.2 0.3 0.4
10−3
10−2
10−1
100
Figure 9.18: As Figure 9.16 but the channel length is 20µm
the topic. The errors identified with the solution of the Stefan tube problem were related
to the formulation of the problem, rather than the Lattice Boltzmann scheme. Further
work is needed to identify exactly how a reacting surface behaves when intersected by a
solid boundary, using continuum kinetic theory.
The fully-developed flow problem is applicable to a fuel cell with only products and reac-
tants present, in which case the ratio of the species mass flux must take the stoichiometric
value. This constraint is sufficient to ensure that the total pressure gradient is small, and
the mass flux of hydrogen is driven by the concentration gradient and is independent of
channel width.
The Stefan tube type problem represents a fuel cell with products and an inert (heavy)
species present, while the reactants have been ignored. This problem turned out to be
two-dimensional in nature, with recirculation of the inert gas. If the fuel cell maintains a
constant hydrogen concentration at the reacting surface, then a wider channel will pass a
higher hydrogen mass flux per unit area.
140
Chapter 10
Conclusions
This chapter presents a summary of the work contained in this thesis along with sugges-
tions for future research. Overall, the work should be considered as proof of the viability
of using LB methods for multi-component gases, although it does not yet constitute a
comprehensive tool. Future work should follow the spirit of the current work in seeking a
fundamental approach to the solution of problems using the Lattice Boltzmann method.
10.1 Summary of results from this work
A key theme of this work has been the two viewpoints from which LB can be approached,
namely the ‘particle’ basis and the ‘moment’ basis. Exactly the same method can be
described from these entirely different starting points, and this makes the LB literature
needlessly complicated. This thesis hopefully adds to the growing number of voices which
claim that using the moment basis is the most effective way to view the LB method.
The single-component one-dimensional scheme in Chapter 2 demonstrated the impor-
tance of having unique second order moments in order to recover the correct momentum
equation. The D2Q9 scheme has this property but not all the third order moments are
unique, and the D2Q9 scheme can be thought of as representing the truncated nine moment
system as shown in Table 4.1. Upon further inspection, the second order moments are only
141
truly unique when a non-physical equation of state is used. Otherwise, a physical equa-
tion of state connects the density to the pressure, and hence the zero to the second order
moments. It is this link which destabilises the scheme, as demonstrated by the expression
for the momentum flux density tensor, Equation (4.44). A single-component model for the
steady-state compressible Navier-Stokes equations was developed in Chapter 4 which has
a physical equation of state.
Boundary conditions are commonly conceived and described in the particle basis, while
specifying some conditions in the moment basis and others in the particle basis leads to a
hybrid which is unsatisfactory. A new method of analysis, based entirely in the moment
basis, was introduced in Chapter 5. The key result is that, although three moments must be
given values at a straight boundary, the choice of which three moments to use is restricted.
This analysis may be used to explore existing schemes, such as the ‘bounce-back’ scheme
which introduces an unspecified slip velocity at the boundary.
The results for Poiseuille flow in Chapter 6, showed good agreement with the analytical
solution. However, the use of a physical equation of state increases computational time
by placing a lower limit on the size of the time step. The stability range for the bulk and
shear viscosities is such that they must have comparable values. Finally, it was observed
that the maximum physical grid spacing is small, for example for hydrogen it is around
3× 10−7m.
The best multi-component model examined in Chapter 7 was the Mixture Centred
Equilibrium (MCE), where all the moments of the equilibrium function are centred on the
mixture velocity. This provided the greatest range of stability compared with the other
equilibrium functions described. The species momentum flux density tensor is a func-
tion of both the species and mixture velocities, but the momentum flux density tensor for
the mixture is the same as that of a single component gas. The steady-state compressible
Navier-Stokes equations are recovered for the mixture, with a viscosity defined for the mix-
ture but not for the species. The diffusive behaviour is a result of the first order moments
of the collision integral and can be controlled separately from the viscous behaviour. The
Stefan-Maxwell diffusion equations are recovered, plus ‘pressure diffusion’ terms, however
142
pressure diffusion was small in all cases considered. Diffusive slip boundary conditions were
devised using the moment based method for boundary conditions, and validated against a
one-dimensional numerical solution of the traditional fluid dynamic equations.
For real applications the use of the variable transformation scheme of Chapter 8 is
essential for reducing computational time and, in the example considered, did so by a factor
of approximately fifty. This involved no changes to the model of Chapter 7, but required
significant changes to the implementation. For equilibrium functions other than the MCE,
the variable transformation would become even more complicated. Pre-conditioning shows
promise for further reducing the computational time by changing the time derivate terms of
the momentum equation. Further work is needed to optimise the pre-conditioning process.
Finally, in Chapter 9 the LB model was applied to two different diffusion problems,
demonstrating the versatility of the model. The LB solution identified inconsistencies with
the formulation of the Stefan tube type problem, and further work is needed in this area.
The results are relevant to fuel cell design but are far from extensive. For binary channel
flow, diffusion slip boundary conditions must be taken into account when the total pressure
gradient is small.
10.2 Suggestions for future research
The principle obstacle, which must be overcome before LB methods are considered to be
useful for multi-component gas modelling, is the limited grid size. It is with the greatest
regret that the author must concede that the current work has barely scratched the surface
of uncovering the causes of this limitation. Having a physically correct equation of state
introduced the limitation on grid size. Therefore, one must assume that the instability is
caused by the erroneous terms in the momentum flux density tensor which lead to changes
in the effective bulk viscosity. This is not the only possible explanation; Dellar (2002b)
identified an instability mode which was not visible in a long-wavelength analysis. The
instability with large grid spacing could be similarly invisible. However, concluding that
the instability is most likely caused by the bulk viscosity highlights the next failure of the
143
model, the inability to set the bulk viscosity to zero. For most engineering applications, the
bulk viscosity does not play a role (see Stewartson (1964)), and the ability to remove bulk
viscosity terms is desirable. The intention of the MRT collision operator used in this work
was to separate bulk and shear viscosity terms. However, the model was not stable for
values of bulk viscosity substantially different from those which correspond to a diagonal
matrix Λ. A better method of controlling bulk viscosity is therefore needed.
Perhaps the simplest solution would be to make the third order moments independent,
either by changing the structure of the D2Q9 model or by increasing the number of ve-
locities. Many alternative velocity sets exist, but it is not clear which is preferable and
additional velocities do increase computational requirements. It seems that nine velocities
could be sufficient if they are different from the usual D2Q9 configuration. For example, if
the ‘straight’ (f1, f2, f3 and f4) velocities have a magnitude of 2 while the ‘diagonal’ (f5,
f6, f7 and f8) velocities remain the same, then Π˜X is no longer equal to Π˜XXX but instead
Π˜XXX = 4Π˜X − 3Π˜XY Y . If an appropriate velocity set is found that has velocities which
do not stream to the nearest neighbour, boundary conditions must be applied to nodes
that are not directly on the boundary. Using the moment based method for boundary
conditions sensible conditions for the additional nodes may be devised.
If a new velocity set cannot be utilised, it may be possible to modify the existing D2Q9
scheme to remove the erroneous bulk viscosity terms using a finite difference scheme, such
as those proposed by Toelke (2006) and Prasianakis and Karlin (2008). The aim is to use
a finite difference approximation to the velocity gradients to calculate the bulk viscosity
terms and then to subtract them from the momentum flux density tensor during the
collision step. This author’s attempts to reproduce this work have not succeeded, but it
is a promising, if somewhat cumbersome, idea. However, this seems to be addressing the
symptoms of the problem not the cause.
If, by whatever means, a scheme were developed which was stable at larger grid spacings,
this would allow a re-examination of the multi-component scheme. The idea behind the
MCE, SCE and CSCE equilibrium functions was that the mixture or species velocity would
appear in the viscous terms. However, the expressions for the viscous stress tensor using
144
any of the equilibrium functions contained both mixture and species velocities, and all were
unsatisfactory. Furthermore, it is not clear that the method used to derive the expressions
for the species momentum flux density tensors was entirely appropriate. These are minor
points however, as the correct macroscopic equation for the mixture was recovered for the
MCE and CSCE.
Extending the model developed in this thesis to gases with more than two components
is feasible. Once again, the mixture Navier-Stokes equations could be used along with
the Stefan-Maxwell equation for the diffusion. All that is required is an extension to the
collision terms.
Aside from the development of the model, this thesis has highlighted the fact that there is
much that is not understood about gas mixtures. Chapter 9 was far from a comprehensive
review of gas mixtures. Fuel cells rarely have only two gas species present and further
inert species can influence the behaviour of the reacting species, as seen in the Stefan tube
problem. The behaviour of a gas mixture immediately adjacent to a reacting surface has
not been extensively studied, as was seen with the boundary conditions for the Stefan tube
type problem. Within a mean free path distance from the surface, the gases do not behave
as a continuum. Examining this layer would result in macroscopic conditions that could be
applied to the LB model, in a similar manner as the macroscopic conditions for diffusion
slip.
In summary, for those wishing to follow this work there is much for mathematicians,
physicists and engineers to do.
145
References
Ahrenholz, B., Toelke, J., and Krafczyk, M. (2006). Lattice-Boltzmann simulations in
reconstructed parametrized porous media. Int. J. Comput. Fluid Dyn., 20(6):369–377.
Andries, P., Aoki, K., and Perthame, B. (2002). A consistent BGK-type model for gas
mixtures. J. Statist. Phys., 106(5-6):993–1018.
Arcidiacono, S., Ansumali, S., Karlin, I. V., Mantzaras, J., and Boulouchos, K. B. (2006a).
Entropic lattice Boltzmann method for simulation of binary mixtures. Math. Comput.
Simulat., 72(2-6):79–83.
Arcidiacono, S., Karlin, I. V., Mantzaras, J., and Frouzakis, C. E. (2007). Lattice Boltz-
mann model for the simulation of multicomponent mixtures. Phys. Rev. E, 76(4, Part
2):046703.
Arcidiacono, S., Mantzaras, J., Ansumali, S., Karlin, I. V., Frouzakis, C., and Boulouchos,
K. B. (2006b). Simulation of binary mixtures with the lattice Boltzman method. Phys.
Rev. E, 74(5):056707.
Asinari, P. (2005). Viscous coupling based lattice Boltzmann model for binary mixtures.
Phys. Fluids, 17(6):067102.
Asinari, P. (2006). Semi-implicit-linearized multiple-relaxation-time formulation of lattice
Boltzmann schemes for mixture modeling. Phys. Rev. E, 73(5, Part 2):056705.
Asinari, P. (2008). Multiple-relaxation-time lattice Boltzmann scheme for homogeneous
mixture flows with external force. Phys. Rev. E, 77(5, Part 2):056706.
146
Asinari, P. (2009). Lattice Boltzmann scheme for mixture modeling: Analysis of the con-
tinuum diffusion regimes recovering Maxwell-Stefan model and incompressible Navier-
Stokes equations. Phys. Rev. E, 80(5, Part 2):056701.
Asinari, P. and Ohwada, T. (2009). Connection between kinetic methods for fluid-dynamic
equations and macroscopic finite-difference schemes. Comput. Math. Appl., 58(5):841–
861.
Asinari, P., Quaglia, M. C., von Spakovsky, M. R., and Kasula, B. V. (2007). Direct
numerical calculation of the kinematic tortuosity of reactive mixture flow in the an-
ode layer of solid oxide fuel cells by the lattice Boltzmann method. J. Power Sources,
170(2):359–375.
Baehr, H. D. and Stephan, K. (1998). Heat and Mass-transfer. Springer.
Bennett, C. O. and Myers, J. E. (1982). Momentum, Heat and Mass Transfer. Mcgraw
Hill Chemical Engineering Series.
Bhatnagar, P., Gross, E., and Krook, M. (1954). A model for collision processes in gases 1:
Small amplitude processes in charged and neutral one-component systems. Phys. Rev.,
94(3):511–525.
Bird, III, R. B., Stewart, W. E., and Lightfoot, E. N. (1960). Transport Phenomena. John
Wiley & Sons.
Brush, S. G. (1976). The Kind of Motion We Call Heat, volume 2. North Holland.
Cengel, Y. A. (1998). Introduction To Thermodynamics and Heat Transfer. McGraw-Hill
Higher Education.
Chapman, S. and Cowling, T. G. (1970). The Mathematical Theory of Non-Uniform Gases.
Cambridge University Press.
Chen, H. D., Chen, S. Y., and Matthaeus, W. H. (1992). Recovery of the Navier-Stokes
equations using a lattice-gas Boltzmann method. Phys. Rev. A, 45(8):5339–5342.
147
Chen, S. and Doolen, G. D. (1998). Lattice Boltzmann method for fluid flows. Annu. Rev.
Fluid Mech., 30:329–364.
Chen, S. Y., Martinez, D., and Mei, R. W. (1999). On boundary conditions in lattice
Boltzmann methods. Phys. Fluids, 8(9):2527–2536.
Dellar, P. (2002a). Nonhydrodynamic modes and a priori construction of shallow water
lattice Boltzmann equations. Phys. Rev. E, 65(3, Part 2B).
Dellar, P. J. (2001). Bulk and shear viscosities in lattice Boltzmann equations. Phys. Rev.
E, 64(3, Part 1):031203.
Dellar, P. J. (2002b). Nonhydrodynamic modes and a priori construction of shallow water
lattice Boltzmann equations. Phys. Rev. E, 65(3, Part 2B):036309.
Dellar, P. J. (2003). Incompressible limits of lattice Boltzmann equations using multiple
relaxation times. J. Comput. Phys., 190(2):351–370.
d’Humie`res, D. (1992). Generalized lattice-Boltzmann equations. In Shizgal, B. D. and
Weaver, D. P., editors, Rarefied Gas Dynamics: Theory and Simulations, volume 159 of
Prog. Astronaut. Aeronaut., pages 450–458, Washington, D.C. AIAA.
Dullien, F. A. L. and Scott, D. S. (1962). The flux-ratio for binary counterdiffusion of Ideal
Gases. Chem. Eng. Sci., 17(10):771–775.
Eckert, E. R. G. and Drake, Jr, R. M. (1959). Heat and Mass Transfer. McGraw-Hill Book
Company.
Evans, III, R. B., Truitt, J., and Watson, G. M. (1961a). Interdiffusion of helium and
argon in a large-pore graphite. J. Chem. Eng. Data, 6(4):522–525.
Evans, III, R. B., Watson, G., and Mason, E. (1962a). Gaseous diffusion in porous media
.2. effect of pressure gradients. J. Chem. Phys., 36(7):1894–1902.
148
Evans, III, R. B., Watson, G. M., and Mason, E. A. (1961b). Gaseous diffusion in porous
media at uniform pressure. J. Chem. Phys., 35(6):2076–2083.
Evans, III, R. B., Watson, G. M., and Truitt, J. (1962b). Interdiffusion of gases in a low
permeability graphite at uniform pressure. J. Appl. Phys., 33(9):2682–2688.
Evans, III, R. B., Watson, G. M., and Truitt, J. (1963). Interdiffusion of gases in a low-
permeability graphite .2. influence of pressure gradients. J. Appl. Phys., 34(7):2020–2026.
Facin, P. C., Philippi, P. C., and dos Santos, L. O. E. (2004). A non-linear lattice-
Boltzmann model for ideal miscible fluids. Future Gener. Comp. Sy., 20(6):945–949.
Geankoplis, C. J. (1993). Transport Processes and Unit Operations. PennWell Books.
Ginzburg, I. and Adler, P. M. (1994). Boundary flow condition analysis for the 3-
Dimensional lattice Boltzmann model. J. Phys. II, 4(2):191–214.
Ginzburg, I. and d’Humieres, D. (2003). Multireflection boundary conditions for lattice
Boltzmann models. Phys. Rev. E, 68(6, Part 2):066614.
Ginzburg, I., Verhaeghe, F., and d’Humieres, D. (2008). Study of simple hydrodynamic
solutions with the two-relaxation-times lattice Boltzmann scheme. Commun. Comput.
Phys., 3(3):519–581.
Grad, H. (1963). Asymptotic theory of the boltzmann equation. Phys. Fluids, 6(2):147–181.
Graham, T. (1876). Chemical and Physical Researches. [Printed by T. and A. Constable].
Guo, Z., Asinari, P., and Zheng, C. (2009). Lattice boltzmann equation for microscale gas
flows of binary mixtures. Phys. Rev. E, 79(2, Part 2):026702.
Guo, Z., Zhao, T. S., and Shi, Y. (2004). Preconditioned lattice-Boltzmann method for
steady flows. Phys. Rev. E, 70(6):066706.
Guo, Z. L. and Zhao, T. S. (2003). Discrete velocity and lattice Boltzmann models for
binary mixtures of nonideal fluids. Phys. Rev. E, 68(3, Part 2):035302.
149
Hamel, B. B. (1966). Two-fluid hydrodynamic equations for a neutral disparate-mass
binary mixture. Phys. Fluids, 9(1):12–22.
He, X., Chen, S., and Doolen, G. D. (1998). A novel thermal model of the lattice Boltzmann
method in incompressible limit. J. Comput. Phys., 146:282–300.
He, X. Y. and Luo, L. S. (1997a). Lattice Boltzmann model for the incompressible Navier-
Stokes equation. J. Statist. Phys., 88(3-4):927–944.
He, X. Y. and Luo, L. S. (1997b). A priori derivation of the lattice Boltzmann equation.
Phys. Rev. E, 55(6, Part A):6333–6336.
He, X. Y. and Luo, L. S. (1997c). Theory of the lattice Boltzmann method: From the
Boltzmann equation to the lattice Boltzmann equation. Phys. Rev. E, 56(6):6811–6817.
He, X. Y., Zou, Q. S., Luo, L. S., and Dembo, M. (1997). Analytic solutions of simple flows
and analysis of nonslip boundary conditions for the lattice Boltzmann BGK model. J.
Statist. Phys., 87(1-2):115–136.
Higuera, F. J., Succi, S., and Benzi, R. (1989). Lattice gas-dynamics with enhanced
collisions. Europhys. Lett., 9(4):345–349.
Hollis, A. P., Halliday, I., and Care, C. M. (2008). An accurate and versatile lattice closure
scheme for lattice Boltzmann equation fluids under external forces. J. Comput. Phys.,
227(17):8065–8082.
Holman, J. P. (1988). Heat Transfer. McGraw Hill Higher Education.
Hoogschagen, J. (1955). Diffusion in porous catalysts and adsorbents. Ind. Eng. Chem.,
47(5):906–913.
Ikenberry, E. and Truesdell, C. (1956). On the pressures and the flux of energy in a gas
according to Maxwell’s kinetic theory, I. J. Ratl. Mech. Anal., 5:1–54.
150
Inamuro, T., Yoshino, M., and Ogino, F. (1995). Non-slip boundary-condition for lattice
Boltzmann simulations. Phys. Fluids, 7(12):2928–2930.
Incropera, F. P., DeWitt, D. P., Bergman, T. L., and Lavine, A. S. (2001). Fundamentals
of Heat and Mass Transfer. John Wiley & Sons.
Ivchenko, I. N., Loyalka, S. K., and Tompson, R. V. (2002). Boundary slip phenomena in
a binary gas mixture. Z. Angew. Math. Phys., 53(1):58–72.
Izquierdo, S. and Fueyo, N. (2009). Optimal preconditioning of lattice Boltzmann methods.
J. Comput. Phys., 228(17):6479–6495.
Jackson, R. (1977). Transport in Porous Catalysts (Chemical engineering monographs).
Elsevier Science Ltd.
Joshi, A. S., Grew, K. N., Peracchio, A. A., and Chiu, W. K. S. (2007a). Lattice Boltz-
mann modeling of 2D gas transport in a solid oxide fuel cell anode. J. Power Sources,
164(2):631–638.
Joshi, A. S., Peracchio, A. A., Grew, K. N., and Chiu, W. K. S. (2007b). Lattice Boltzmann
method for multi-component, non-continuum mass diffusion. J. Phys. D: Appl. Phys.,
40(23):7593–7600.
Junk, M. and Yang, Z. X. (2005). One-point boundary condition for the lattice Boltzmann
method. Phys. Rev. E, 72(6, Part 2):066701.
Kim, S. H., Pitsch, H., and Boyd, I. D. (2009). Lattice Boltzmann modeling of multicom-
ponent diffusion in narrow channels. Phys. Rev. E, 79(1, Part 2):016702.
Knaff, G. and Schlunder, E. U. (1985). Experimental confirmation of graham law of
diffusion up to pore diameters of 2 µm. Chem. Eng. Process., 19(3):167–173.
Kramers, H. A. and Kistemaker, J. (1943). On the slip of a diffusing gas mixture along a
wall. Physica, 10:699–713.
151
Lallemand, P. and Luo, L. S. (2000). Theory of the lattice Boltzmann method: Dispersion,
dissipation, isotropy, Galilean invariance, and stability. Phys. Rev. E, 61(6, Part A):6546–
6562.
Landau, L. D. and Lifshitz, E. M. (1959). Fluid Mechanics. Pergamon Press.
Latt, J. and Chopard, B. (2006). Lattice Boltzmann method with regularized pre-collision
distribution functions. Math. Comput. Simulat., 72(2-6):165–168.
Latt, J., Chopard, B., Malaspinas, O., Deville, M., and Michler, A. (2008). Straight velocity
boundaries in the lattice Boltzmann method. Phys. Rev. E, 77(5, Part 2):056703.
Lee, J. H., Heo, J. W., Lee, D. S., Kim, J., Kim, G. H., Lee, H. W., Song, H. S., and Moon,
J. H. (2003). The impact of anode microstructure on the power generating characteristics
of SOFC. Solid State Ionics, 158(3-4):225–232.
Lee, J. H., Moon, H., Lee, H. W., Kim, J., Kim, J. D., and Yoon, K. H. (2002). Quantitative
analysis of microstructure and its related electrical property of SOFC anode, Ni-YSZ
cermet. Solid State Ionics, 148(1-2):15–26.
Lee, K. R., Choi, S. H., Kim, J., Lee, H. W., and Lee, J. H. (2005). Viable image analyzing
method to characterize the microstructure and the properties of the Ni/YSZ cermet
anode of SOFC. J. Power Sources, 140(2):226–234.
Loyalka, S. K. (1971). Velocity slip coefficient and the diffusion slip velocity for a multi-
component gas mixture. The Phys. Fluids, 14(12):2599–2604.
Luo, L. S. (2004). Comment on “Discrete Boltzmann equation for microfluidics”. Phys.
Rev. Letters, 92(13):139401.
Luo, L. S. and Girimaji, S. S. (2002). Lattice Boltzmann model for binary mixtures. Phys.
Rev. E, 66(3, Part 2A):035301.
Maier, R. S., Bernard, R. S., and Grunau, D. W. (1996). Boundary conditions for the
lattice Boltzmann method. Phys. Fluids, 8(7):1788–1801.
152
Mason, E. A. and Evans, III, R. B. (1969). Graham’s laws: Simple demonstrations of gases
in motion: Part I, Theory. J. Chem. Educ., 46(6):358–364.
Mason, E. A., Love, L. D., and Evans, III, R. B. (1969). Graham’s laws: Simple demon-
strations of gases in motion: Part II, Experiments. J. Chem. Educ., 46(7):423–427.
McCracken, M. E. and Abraham, J. (2005). Lattice Boltzmann methods for binary mix-
tures with different molecular weights. Phys. Rev. E, 71(4, Part 2):046704.
McNamara, G. R., Garcia, A. L., and Alder, B. J. (1995). Stabilization of thermal lattice
Boltzmann models. J. Statist. Phys., 81(1-2):395–408.
Meyer, J. P. and Kostin, M. D. (1975). Circulation phenomena in stefan diffusion. Int. J.
Heat Mass Transfer, 18(11):1293–1297.
Mills, A. F. (2003). On steady one-dimensional diffusion in binary ideal gas mixtures. Int.
J. Heat Mass Transfer, 46(13):2495–2497.
Mills, A. F. (2007). Diffusion creep. Int. J. Heat Mass Transfer, 50(25-26):5087–5098.
Mo, G. and Rosenberger, F. (1991). Molecular-Dynamics simulations of flow with bi-
nary diffusion in a 2-Dimensional channel with atomically rough walls. Phys. Rev. A,
44(8):4978–4985.
Noble, D. R., Chen, S. Y., Georgiadis, J. G., and Buckius, R. (1995a). A consistent
hydrodynamic boundary-condition for the Lattice Boltzmann method. Phys. Fluids,
7(1):203–209.
Noble, D. R., Georgiadis, J. G., and Buckius, R. O. (1995b). Direct assessment of lattice
Boltzmann hydrodynamics and boundary-conditions for recirculating-flows. J. Statist.
Phys., 81(1-2):17–33.
Noever, D. (1990a). A note on the no-slip condition applied to diffusing gases. Phys. Lett.
A, 144(4-5):253–255.
153
Noever, D. A. (1990b). The baroeffect and an appropriate momentum boundary-condition.
Phys. Fluids A-Fluid Dynamics, 2(5):858–865.
Noever, D. A. (1990c). Ternary baroeffect with a nondiffusing component. Phys. Rev.
Letters, 65(13):1587–1596.
Poling, B. E., Prausnitz, J. M., and O’Connell, J. (2000). The Properties of Gases and
Liquids. McGraw-Hill Professional, 5th edition.
Prasianakis, N. I. and Karlin, I. V. (2008). Lattice Boltzmann method for simulation of
compressible flows on standard lattices. Phys. Rev. E, 78(1, Part 2):016704.
Rama, P., Liu, Y., Chen, R., Ostadi, H., Jiang, K., Gao, Y., Zhang, X., Fisher, R., and
Jeschke, M. (2010). Multiscale modeling of single-phase multicomponent transport in the
cathode gas diffusion layer of a polymer electrolyte fuel cell. Energy Fuels, 24:3130–3143.
Reis, T. and Phillips, T. N. (2008). Alternative approach to the solution of the dispersion
relation for a generalized lattice Boltzmann equation. Phys. Rev. E, 77(2, Part 2):026702.
Remick, R. R. and Geankoplis, C. J. (1973). Binary diffusion of gases in capillaries in
transition region between Knudsen and molecular diffusion. Ind. Eng. Chem. Fundam.,
12(12):214–220.
Salcedo-Diaz, R., Ruiz-Femenia, R., Kerkhof, P. J. A. M., and Peters, E. A. J. F. (2008).
Velocity profiles and circulation in Stefan-diffusion. Chem. Eng. Sci., 63(19):4685–4693.
Sirovich, L. (1962). Kinetic modeling of gas mixtures. Phys. Fluids, 5(8):908–918.
Skordos, P. A. (1993). Initial and boundary-conditions for the lattice Boltzmann method.
Phys. Rev. E, 48(6):4823–4842.
Sofonea, V. and Sekerka, R. (2005a). Boundary conditions for the upwind finite difference
Lattice Boltzmann model: Evidence of slip velocity in micro-channel flow. J. Comput.
Phys., 207(2):639–659.
154
Sofonea, V. and Sekerka, R. F. (2001). BGK models for diffusion in isothermal binary fluid
systems. Physica A, 299(3-4):494–520.
Sofonea, V. and Sekerka, R. F. (2005b). Diffusivity of two-component isothermal finite
difference lattice Boltzmann models. Int. J. Mod. Phys. C, 16(7):1075–1090.
Stewartson, K. (1964). The Theory of Laminar Boundary Layers in Compressible Fluids.
Oxford University Press, Oxford.
Succi, S., Benzi, R., and Higuera, F. (1991). The Lattice Boltzmann-Equation - A New
Tool For Computational Fluid-Dynamics. Physica D, 47(1-2):219–230.
Succi, S., Karlin, I. V., and Chen, H. (2002). Role of the H theorem in lattice Boltzmann
hydrodynamic simulations. Rev. Mod. Phys., 74:1203–1220.
Toelke, J. (2006). A thermal model based on the lattice Boltzmann method for low Mach
number compressible flows. J. Comput. Theor. Nanosci., 3(4):579–587.
Torrilhon, M. and Struchtrup, H. (2008). Boundary conditions for regularized 13-moment-
equations for micro-channel-flows. J. Comput. Phys., 227(3):1982–2011.
VanKampen, N. G. (1987). Chapman-Enskog as an application of the method for elimi-
nating fast variables. J. Statist. Phys., 46:709–727.
Verberg, R. and Ladd, A. J. C. (1999). Simulation of low-Reynolds-number flow via a
time-independent lattice-Boltzmann method. Phys. Rev. E, 60:3366–3373.
Verhaeghe, F., Luo, L.-S., and Blanpain, B. (2009). Lattice Boltzmann modeling of mi-
crochannel flow in slip flow regime. J. Comput. Phys., 228(1):147–157.
Vincenti, W. G. and Kruger, C. H. (1976). Introduction to Physical Gas Dynamics. Krieger
Publishing Company.
Wang, L. and Afsharpoya, B. (2006). Modeling fluid flow in fuel cells using the lattice-
Boltzmann approach. Math. Comput. Simulat., 72(2-6):242–248.
155
Welty, J., Wicks, C. E., Rorrer, G. L., and Wilson, R. E. (1984). Fundamentals of Mo-
mentum, Heat and Mass Transfer. John Wiley & Sons.
Whitaker, S. (1991). Role of the species momentum equation in the analysis of the Stefan
diffusion tube. Ind. Eng. Chem. Res., 30(5):978–983.
White, F. M. (1988). Heat and Mass Transfer. Prentice Hall.
Young, J. B. and Todd, B. (2005). Modelling of multi-component gas flows in capillaries
and porous solids. Int. J. Heat Mass Transfer, 48(25-26):5338–5353.
Yu, D. Z., Mei, R. W., and Shyy, W. (2005). Improved treatment of the open boundary
in the method of lattice Boltzmann equation. Prog. Comput. Fluid Dy., 5(1-2):3–12.
Yu, H. and Zhao, K. (2000). Lattice Boltzmann method for compressible flows with high
Mach numbers. Phys. Rev. E, 61:3867–3870.
Ziegler, D. P. (1993). Boundary-conditions for lattice Boltzmann simulations. J. Statist.
Phys., 71(5-6):1171–1177.
Zou, Q. S. and He, X. Y. (1997). On pressure and velocity boundary conditions for the
lattice Boltzmann BGK model. Phys. Fluids, 9(6):1591–1598.
156