Abstract

Faculte des sciences

Professeur Bastien Chopard
Professeur Michel Droz

`
Universite de Geneve
Dpartement dinformatique
e
Dpartement de physique
e

Hydrodynamic Limit of Lattice

Boltzmann Equations

`
THESE
prsente ` la Facult des sciences de lUniversit de Gen`ve
e
e a
e
e
e
pour obtenir le grade de Docteur `s sciences, mention interdisciplinaire
e

par

Jonas Ltt
a
de Mhledorf (SO)
u

Th`se No 3846
e

Gen`ve
e
Atelier de reproduction de la Section de physique
2007


FACULTE DES SCIENCES
Doctorat `s Sciences
e
mention interdisciplinaire

Th`se de Monsieur Jonas LATT

e
intitule:
e

Hydrodynamic Limit of Lattice Boltzmann Equations

La facult des sciences, sur le pravis de Messieurs B. CHOPARD, professeur adjoint
e
e
et codirecteur de th`se (Dpartement dinformatique), M. DROZ, professeur titulaire
e
e
et codirecteur de th`se (Dpartement de physique thorique), S. SUCCI, professeur
e
e
e
(Istituto per le Applicazioni del Calcolo Mauro Picone, Rome, Italie), autorise
limpression de la prsente th`se, sans exprimer dopinion sur les propositions qui y
e
e
sont mises.
e

Gen`ve, le 23 avril 2007

e

Th`se -3846e
Le Doyen, Pierre SPIERER

This dissertation was presented on March 9, 2007, at the Faculty of Science of the
University of Geneva.
The members of the comittee were the following:
Prof. Bastien Chopard, Supervisor
Prof. Michel Droz, Supervisor
Prof. Sauro Succi, Istituto Applicazioni Calcolo, Rome

Remerciements
En prenant mod`le sur la dynamique des uides mergeante dune microscopie riche
e
e
et abondante, cette th`se est issue dun tissu dinteractions complexe avec mon
e
entourage qui ma soutenu ` un niveau scientque, technique et motionnel. Je
a
e
suis reconnaissant ` tous les gens qui mont suivi durant ma th`se et esp`re quils
a
e
e
trouveront plaisir, en lisant ce document, ` y reconnaitre une contribution, une ide
a
e
ou un biais qui leur est d.
u
Je pense en particulier ` mes directeurs de th`se Bastien Chopard et Michel Droz
a
e
qui mont suivi au long du travail, et ` Sauro Succi qui a t conseiller et jur de
a
ee
e
th`se. Des conversations stimulantes ont permis de situer plusieurs sujets dans un
e
contexte plus large, et jaimerais faire valoir entre autres les contributions de Peter
Wittwer, de Vincent Heuveline et de Michel Deville.
Jai prot dun environnement de travail inspirant grce ` Jean-Luc et Bernhard
e
a a
qui mont aid ` aronter le monde hostile de linformatique, Fokko qui a soigneuseea
ment relu ma th`se, Orestis qui a prt une oreille critique ` mes ides farfelues et
e
ee
a
e
mes autres coll`gues qui mont tous donn un soutien remarquable.
e
e
La contribution la plus prcieuse vient sans aucun doute de mon entourage
e
proche, et je remercie Roxane, ma m`re, les parents de Roxane et toute ma famille
e
pour leur patience et leur appui.

Contents
Abstract

Rsum en Franais
e
e
c

vii

Introduction

1 Computational uid dynamics and lattice Boltzmann

1.1 Vector and tensor notation . . . . . . . . . . . . . . . .
1.2 Macroscopic equations . . . . . . . . . . . . . . . . . .
1.2.1 Governing equations . . . . . . . . . . . . . . .
1.2.2 Isothermal and incompressible ows . . . . . . .
1.2.3 Dimensionless formulation . . . . . . . . . . . .
1.2.4 Further reading . . . . . . . . . . . . . . . . . .
1.3 Lattice Boltzmann method . . . . . . . . . . . . . . . .
1.3.1 Discrete variables . . . . . . . . . . . . . . . . .
1.3.2 Lattice Boltzmann model . . . . . . . . . . . .
1.3.3 Lattice structures . . . . . . . . . . . . . . . . .
1.3.4 Further reading . . . . . . . . . . . . . . . . . .
2 Multiscale Chapman-Enskog analysis
2.1 Series expansion with scale separation . . . . .
2.1.1 Lattice symmetries . . . . . . . . . . .
2.1.2 Multi-scale expansion . . . . . . . . . .
2.1.3 Conservation laws . . . . . . . . . . . .
2.2 Advection-diusion: Chapman-Enskog ansatz
2.3 Fluid ows: Chapman-Enskog ansatz . . . . .
2.4 Dimensionless formulation and accuracy . . .
2.4.1 Dimensionless formulation . . . . . . .
2.4.2 Accuracy of LB methods . . . . . . . .
2.4.3 Accuracy of specic models . . . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

3 Corrections to the BGK dynamics

3.1 Advection-diusion revisited . . . . . . . . . . . .
3.1.1 Second-order time stepping . . . . . . . .
3.1.2 Numerical verication . . . . . . . . . . .
3.2 Bulk viscosity for uid ows . . . . . . . . . . . .
3.2.1 Numerical error in incompressible ows . .
3.2.2 Correction term for the bulk viscosity . . .
3.3 Alternative LB models . . . . . . . . . . . . . . .
3.3.1 LB models with multiple relaxation times
3.3.2 Entropic LB models . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

5
5
6
6
7
8
8
9
9
10
11
13

.
.
.
.
.
.
.
.
.
.

15
15
15
17
18
21
23
26
27
28
28

.
.
.
.
.
.
.
.
.

31
31
31
33
34
34
36
39
39
41

ii

Contents

4 Regularized lattice Boltzmann

4.1 Introduction . . . . . . . . . . . . . . .
4.2 Regularization procedure . . . . . . . .
4.2.1 Regularized particle distribution
4.2.2 Dimensionless formulation . . .
4.2.3 Algorithm . . . . . . . . . . . .
4.2.4 Related models . . . . . . . . .
4.3 RLB and MRT . . . . . . . . . . . . .
4.4 Numerical verication . . . . . . . . .

. . . . . .
. . . . . .
functions
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .

5 Boundary and initial conditions

5.1 Introduction . . . . . . . . . . . . . . . . . .
5.2 Boundary condition . . . . . . . . . . . . . .
5.2.1 Overview . . . . . . . . . . . . . . .
5.2.2 Implementation of the pressure . . .
5.2.3 Regularization of on-wall distribution
5.2.4 Numerical verication . . . . . . . .
5.3 Initial condition . . . . . . . . . . . . . . . .
5.3.1 Introduction . . . . . . . . . . . . . .
5.3.2 The benchmark . . . . . . . . . . . .
5.4 Numerical results . . . . . . . . . . . . . . .
5.4.1 Time evolution of the ow . . . . . .
5.4.2 Benchmark values . . . . . . . . . . .
5.5 Setup of the initial condition . . . . . . . . .
5.6 Discussion: choice of the time step . . . . .
6 Adaptive space and time steps
6.1 Introduction . . . . . . . . . . . . . . . . .
6.2 Grid renement . . . . . . . . . . . . . . .
6.2.1 Introduction . . . . . . . . . . . . .
6.2.2 Analytical prole on the boundaries
6.2.3 Numerical implementation . . . . .
6.2.4 Simulation results . . . . . . . . . .
6.3 Adaptive time . . . . . . . . . . . . . . . .
6.3.1 Introduction . . . . . . . . . . . . .
6.3.2 Numerical experiment . . . . . . .
7 Coupling with other tools of CFD
7.1 Introduction . . . . . . . . . . . . . . .
7.2 FD model and computational grids . .
7.3 Choice of units . . . . . . . . . . . . .
7.4 The FD-LB interface . . . . . . . . . .
7.5 The coupling algorithm . . . . . . . . .
7.5.1 Coupling the pressure . . . . .
7.5.2 Coupling the velocity . . . . . .
7.5.3 Coupling the velocity gradients
7.5.4 Summary . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

. . . . . .
. . . . . .
. . . . . .
. . . . . .
functions
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

43
43
45
45
46
46
47
48
48

.
.
.
.
.
.
.
.
.
.
.
.
.
.

53
53
54
54
55
55
57
59
59
59
61
61
61
63
64

.
.
.
.
.
.
.
.
.

65
65
67
67
68
69
70
71
71
72

.
.
.
.
.
.
.
.
.

75
75
76
78
78
79
79
79
80
80

Contents

iii

7.6 Numerical validation . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

7.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
8 A model for uid turbulence based on kinetic variables
8.1 Turbulence modeling . . . . . . . . . . . . . . . . . . . . .
8.1.1 Overview . . . . . . . . . . . . . . . . . . . . . . .
8.1.2 Large eddy simulations in 3D incompressible uids
8.2 Averaged LB equation . . . . . . . . . . . . . . . . . . . .
8.3 Direct numerical simulation . . . . . . . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

85
85
85
86
88
91

Publications

95

Bibliography

99

Abstract
An analysis of the lattice Boltzmann (LB) method is conducted, and various conclusions are drawn on how to exploit this method for the numerical resolution of
the Navier-Stokes equation. This focus contrasts with the traditional approach
to the LB method, which usually considers complex uids with extended physical
properties. The specic scope of the present thesis is motivated by current practice in the domain of computational uid dynamics. It is indeed observed that
an increasing number of scientists and engineers simply use the LB method as an
alternative to conventional numerical solvers for the Navier-Stokes equation. The
conceptual framework presently available for an application of the LB method in
this area is however still weak. This leads to a number of limitations, to which some
workarounds are presented in the thesis.
In the rst part, the LB method is reformulated such as to put the main emphasis
on macroscopically relevant variables, instead of microscopical quantities originated
from the kinetic theory of gases. From this perspective, the accuracy and the limitations of the lattice Boltzmann method are discussed in simple terms.
The second part introduces novel LB models, rst for the simulation of advectiondiusion problems, and then for the resolution of the Navier-Stokes equation. The
latter model, presented under the name of regularized lattice Boltzmann is shown
to substantially increase the stability and accuracy of LB models in numerical simulation. The model is furthermore found to be innovative due to the fact that it
possesses an accurate dimensionless formulation in terms of macroscopic variables.
The technical implications of this observation are investigated in the third part:
the regularized LB method is used to solve various problems with varying time and
space scales, ranging from the use of spatially rened grids to the coupling of the
LB method with a nite dierence method. The point of view adopted throughout
this discussion is that practically all implementation issues in the LB method can be
addressed properly by systematically applying the results from the corresponding
Chapman-Enskog analysis.
The last part of the thesis deliberately changes the focus and utilizes the original,
non-regularized LB method for the simulation of a complex uid. Based on the microscopic LB variables, a model for the simulation of uid turbulence is formulated.
This document presents a selective choice of the topics investigated during the
thesis, aiming at a succinct overview of the main ndings. It focuses on theoretical
aspects and does not reect another important part of the thesis, related to implementation issues of the numerical models. In can be pointed out in particular that
the thesis led to the creation of three independent software projects. They introduce
new programming paradigms for the implementation of parallel programs and LB
models and provide freely available and well documented open source codes. This
part of the work is summarized at the end of the present document, along with a
listing of achieved publications.

Rsum en Franais
e
e
c
Aperu de la problmatique
c
e
Introduction
La prsente th`se de doctorat propose une tude de la mthode de Boltzmann sur
e
e
e
e
rseau, dnote par lacronyme LB pour son nom anglais de lattice Boltzmann.
e
e e
Diverses conclusions sont tires sur les possibilits dexploiter cette mthode pour
e
e
e
la rsolution numrique de lquation de Navier-Stokes. Cette mani`re restrictive
e
e
e
e
dutiliser la mthode LB doit tre vue par opposition ` son utilisation traditione
e
a
nelle, se situant dans le contexte de ltude de uides complexes avec des proprits
e
ee
physiques tendues. Le choix adopt dans la th`se est motiv par ltat de lart
e
e
e
e
e
dans le domaine de la simulation de uides (CFD pour lexpression anglaise computational uid dynamics). En eet, on observe quun nombre croissant de scientiques et ingnieurs utilisent la mthode de LB comme outil alternatif aux moyens
e
e
de rsolution numrique conventionnels pour lquation de Navier-Stokes. Le cadre
e
e
e
conceptuel ` disposition actuellement pour une application du LB dans ce domaine
a
est par contre toujours faible. Il en ensuit un nombre de limitations, auxquelles
certains rem`des sont proposs dans la th`se.
e
e
e
Dans la premi`re partie, la mthode LB est reformule dune telle mani`re `
e
e
e
e a
mettre en vidence les variables de nature macroscopique, au lieu des quantits
e
e
microscopiques issues de la thorie cintique des gaz. De cette nouvelle perspective,
e
e
la prcision ainsi que les limitations de la mthode de Boltzmann sur rseau sont
e
e
e
tudis en termes simples.
e
e
La deuxi`me partie introduit de nouveaux mod`les LB, dabord pour la simulae
e
tion de probl`mes ` convection-diusion, et ensuite pour la rsolution de lquation
e
a
e
e
de Navier-Stokes. Il est dmontr que ce dernier mod`le, prsent sous le nom de
e
e
e
e
e
lattice Boltzmann rgularis, augmente de mani`re substantielle la stabilit et la
e
e
e
e
prcision des mod`les LB au cours de simulations. Un autre aspect innovant de
e
e
ce mod`le se rf`re au fait quil scrit de mani`re exacte en terme de variables
e
ee
e
e
macroscopiques.
Les implications techniques de cette observation sont tudies dans la troisi`me
e
e
e
partie: la mthode LB rgularise est utilise pour rsoudre divers probl`mes pose
e
e
e
e
e
sdant des chelles temporelles et spatiales variables. Il sagit l` par exemple de
e
e
a
grilles numriques localement ranes ou de couplages entre une mthode LB et
e
e
e
une mthode de dirences nies. Le point de vue adopt ` travers cette tude
e
e
e a
e
sappuie sur lide que pratiquement toute question dimplmentation pratique en
e
e
LB se rsoud ` travers une analyse de Chapman-Enskog correspondante.
e
a
Cette optique est dlibrment modie lors de la derni`re partie de la th`se,
e ee
e
e
e
dans laquelle la mthode LB originale, non-rgularise est utilise pour la simulation
e
e
e
e
dun uide complexe. Sur la base de variables LB microscopiques, un mod`le pour
e
la simulation de uides turbulents est formul.
e

viii

Rsum en Franais
e
e
c

Ce document contient un choix slectif des sujets tudis au cours de la th`se,

e
e
e
e
an de rsumer bri`vement les dcouvertes majeures. Il se concentre sur les aspects
e
e
e
thoriques et omet de dtailler une autre partie importante de la th`se, relative aux
e
e
e
questions dimplmentation numrique. On peut en particulier mettre faire valoir
e
e
le travail investi dans la conception de trois projets de logiciel indpendants, ine
troduisant de nouveaux paradigmes de programmation pour limplantation de programmes parall`les et de mod`les LB. En outre, ces projets proposent des codes
e
e
source libres et pleinement documents. Cette partie du travail est rsume ` la n
e
e
e a
de la th`se, accompagne dune liste des articles publis au cours de la th`se.
e
e
e
e

Simulations en dynamique des uides

La recherche en simulation de dynamique des uides (CFD) se concentre sur une
quation a dirences partielles (PDE pour langlais partial dierential equation)
e
e
connue sous le nom dquation de Navier-Stokes incompressible qui exprime une loi
e
de conservation locale pour la quantit de mouvement du syst`me. Elle reprsente
e
e
e
un mod`le largement simpli pour un uie inviscide ` une seule composante dans un
e
e
a
environnment sans variations de densit ou de temprature, et sans autre ingrdient
e
e
e
physique. Bien que cette quation nadresse que partiellement la complexit de
e
e
la plupart des uides considrs en ingnierie, elle est utilise avec beaucoup de
ee
e
e
succ`s dans dirents domaines an de prdire qualitativement le comportement de
e
e
e
syst`mes immergs dans un uide. Comme exemple, on peut nommer les probl`mes
e
e
e
doptimisation en arodynamique, tel que la conception dune carrosserie de voiture
e
ou la cration dailes davion. Lquation de Navier-Stokes incompressible narrive
e
e
bien videmment ` dcrire la nature de lair compressible que de mani`re tr`s
e
a e
e
e
approximative. Nanmoins, cette quation m`ne ` des prdiction de tr`s bonne
e
e
e a
e
e
qualit quon utilise frquemment, conjointement aux rsultats de la simulation
e
e
e
exprimentale. On peut par ailleurs noter que la rsolution de cette PDE reprsente
e
e
e
une tche de calcul intensive, malgr la simplicit apparente de lquation. On ina
e
e
e
vestit donc beaucoup deorts dans le dveloppement de techniques de rsolution
e
e
numrique, plutt que de rajouter davantage de complexit en ajoutant des termes
e
o
e
physiques additionnels.
Dun autre ct, les domaines de recherche se concentrant sur des uides dits
oe
complexes ont galement une grande importance dans le domaine. Les uides
e
en question possdent une physique enrichie: pour ceux-ci, lapproximation de
e
lquation de Navier-Stokes est donc insusante. Lapproche classique adopte en
e
e
CFD pour traiter ce type de uides consiste en la description de nouvelles proprits
ee
physiques en terme de phnom`nes de transport relies ` une proprit macroe
e
e a
ee
scopique. Une nouvelle PDE est formule pour la dynamique de cette proprit,
e
ee
quon rsout ensuite par une technique numrique approprie. Linteraction entre
e
e
e
les dirents phnom`nes de transport est traite ` travers de termes ajouts dans
e
e
e
e a
e
les PDEs respectives. Dans un uide avec dimportantes variations de temprature
e
par exemple, la temprature est introduite comme une nouvelle proprit, et sa
e
ee
dynamique est dcrite par une quation de transport pour la chaleur. En plus,
e
e
lquation de conservation du moment est tendue pour tenir compte de la come
e
pressibilit du uide. La temprature et la vitesse sont ensuite lis ` travers un
e
e
e a
terme convectif dans lquation de la chaleur et un terme de pousse verticale dans
e
e

Rsum en Franais
e
e
c

ix

lquation du moment. Une quation dtat additionnelle couple la temprature

e
e
e
e
avec la densit. On peut dire que le domaine traditionnel de la CFD adopte une
e
approche de haut niveau ` la modlisation dun uide, en commenant par une loi
a
e
c
de conservation simple, et en rajoutant de nouvelles quations sur demande, bases
e
e
sur des considrations thoriques et phnomnologiques.
e
e
e
e

Mthode de Boltzmann sur rseau

e
e
Le sujet principal de la th`se se traduit par une approche numrique du nom de
e
e
Boltzmann sur rseau (LB pour le terme anglais lattice Boltzmann). Il sagit
e
dun mthode realtivement rcente qui se distingue mthodes traditionelle en adope
e
e
tant une approche de bas niveau ` la modlisation de uides. A cet eet, elle dcrit
a
e
e
le uide ` un niveau molculaire et labore des mod`les pour linteraction entre
a
e
e
e
molcules. La physique compl`te du uide ` un niveau du continu se trouve ime
e
a
plicitement dcrite par le mod`le, et le travail conceptuel quon doit eectuer se
e
e
trouve ` loppos de ce quon ferait dans lapproche de haut niveau classique. En
a
e
eet, les ingrdients physiques du mod`le doivent tre identis un par un et exe
e
e
e
plicits propement, an de permettre une sparation entre proprits relevantes et
e
e
ee
ngligeables. En se basant sur ces considrations, le mod`le de collision est largee
e
e
ment simpli pour les besoins du traitement numrique. En tant que tel, le LB
e
e
est une mthode bien adapte ` la simulation de uides complexes, tant donn
e
e a
e
e
que les ingrdients physiques relevants sont pris en charge tr`s naturellement par la
e
e
mthode. Il est mis en vidence ` quel point une formulation LB dun probl`me de
e
e
a
e
dynamique des uides peut tre lgante lorsquon on simule par exemple des ue
ee
ides ` deux phases. En CFD traditionnelle, une PDE est formule pour chacune des
a
e
deux phases, ainsi que pour la dynamique dinteraction sur linterface entre cellesci. Limplmentation numrique dun tel mod`le requiert lutilisation de techniques
e
e
e
avances en gnie logiciel, an de correctement tracer la position de linterface et
e
e
implmenter la dynamique correspondante. Ces dicults napparaissent pas dans
e
e
la mthode LB, car les deux phases peuvent tre modlises par un mod`le de uide
e
e
e e
e
unique, et linterface entre les phases est traite automatiquement, en ingrdient
e
e
naturel du mod`le.
e

La mthode LB applique au domaine de la CFD

e
e
Lorsque lapproche LB ` la modlisation physique est importe dans le domaine de
a
e
e
CFD classique, deux mondes tr`s dirents, possdant leurs propres convictions et
e
e
e
traditions sont confronts. Ceci m`ne ` de nombreux malentendus, dun ct de
e
e a
oe
la part de la communaut de la CFD qui manque de reconnaitre le potentiel de la
e
mthode LB par sa juste valeur, et dun autre ct de la communaut LB qui prouve
e
oe
e
e
des dicults ` reconnaitres les enjeux rels de la CFD, ou ` comprendre quels sont
e a
e
a
les domaines dans lesquels le LB peut tre comptitifs avec dautres approches.
e
e
Un obstacle majeur ` un change mutuel entre ces deux mthodes vient du fait
a
e
e
quelles dcrivent le uide ` travers un choix de variables dirent. Les variables
e
a
e
centrales en CFD sont les quantits macroscopiques, observables. Dans lquation
e
e
de Navier-Stokes incompressible, les degrs de libert sont dnis par la vitesse et
e
e
e
la pression, dpendants de lespace et du temps. Les outils de rsolution numrique
e
e
e

Rsum en Franais
e
e
c

en CFD classique mod`lisent souvent directement ces variables-ci, cest-`-dire, ils

e
a
reprsentent la valeur adopte par ces variables sur les noeuds dune grille numrique
e
e
e
donne. La mthode LB dun autre ct simule la dynamique dune quantit appele
e
e
oe
e
e
fonction de distribution de particules, drive de la thorie cintique des gaz. Sil
e e
e
e
est bien connu comment cette quantits peut tre convertie en variables dintrt
e
e
ee
macroscopique, la procdure inverse par contre est plus complique. De fait, il
e
e
nest pas toujours clair de quelle mani`re une conguration donne des variables
e
e
macroscopiques est utilise pour initialiser la fonction de distribution. Ce probl`me
e
e
est rarement mis en vidence par une communaut habitue ` penser en terme
e
e
e a
de quantits cintiques uniquement. Il est nanmoins ` la source de complications
e
e
e
a
substantielles lorsque des variables macroscopiques sont voques explicitement dasn
e
e
la simulation, ou lorsque les chelles despace et de temps macroscopiques du mod`les
e
e
sont ajustes.
e

Comment traiter des chelles temporelles et spatiales varie

ables?
Les probl`mes lis ` la conversion entre variables macroscopiques et fonctions de
e
e a
distribution apparaissent par exemple lorsque lon applique une technique de ranement de grille, au cour de laquelle le domaine de la simulation est partitionn en souse
domaines, dont chacun est reprsent par une grille numrique possdant un degr
e
e
e
e
e
de rsolution variable. Le but de cet exercice est dadapter la structure de la grille `
e
a
la gomtrie du probl`me, et par consquent, daugmenter lecacit du calcul. Lors
e e
e
e
e
de lexcution dune simulation, il sav`re ncessaire de communiquer les valeurs cale
e
e
cules au fur et ` mesure dune grille ` lautre. Ce procd sapplique sans dicult
e
a
a
e e
e
dans le contexte des outils de CFD classiques. En eet, ces mthodes gardent en
e
mmoire la valeur sans dimensions des variables macroscopiques, indpendante donc
e
e
de la rsolution de la grille ou du pas temporel, schangeant entre dirente grilles
e
e
e
sans besoin de conversion. La fonction de distribution de particules utilise en LB
e
par contre nest pas adimensionnelle par rapport au choix dunits macroscopiques.
e
Elle doit donc tre remise ` lchelle lors du transfer dune grille ` une autre, puisque
e
a e
a
sa valeur dpend de lespacement de la grille. La r`gle qui dtermine cette remise `
e
e
e
a
lchelle nest aucunement vidente, et doit tre applique avec beaucoup de soins.
e
e
e
e
Par consquent, une certaine confusion r`gne ` ce sujet dans la litrature, et de
e
e
a
e
nombreuses approches direntes ont t proposes, dont les ides de base di`rent
e
ee
e
e
e
substantiellement.
Le ranement nest quun exemple particulier, introduit ici dans un souci dillustration. La gure 1 ache dautres probl`mes ` chelles despace et de temps varie
ae
ables, ncessitant une transformation entre les units macroscopiques et celles du LB:
e
e
limplantation de conditions aux bords et de conditions initiales, et limplantation
de mod`les hybrides, faisant appel en mme temps ` une mthode LB et ` un outil
e
e
a
e
a
classique de la CFD.

Rsum en Franais
e
e
c

xi

Conditions aux bords

Conditions initiales

Couplage avec dautres mthodes

e

Applications multi-grille

Figure 1: Divers probl`mes se posent lors de limplantation dun mod`le LB, pour
e
e
lesquels une conversion entre variables macroscopiques et fonction de distribution
est requise.

Mthode LB et dveloppement multi-chelle

e
e
e
Lquation de Navier-Stokes incompressible
e
Dans ce qui suit, les probl`mes de la dynamique des uides sont rduits ` un point
e
e
a
de vue purement macroscopique. Dans ce contexte, les questions concernant la
simulation dun uide quivalent ` des questions de rsolution numrique dune
e
a
e
e
quation direntielle ` drives partielles. Un uide incompressible et isotherme
e
e
a e e
de viscosit est par exemple dcrit par lquation suivante, connue sous le nom
e
e
e
dquation de Navier-Stokes:
e

p
u
+ u u + I 2S
t
0

= 0.

Elle exprime la conservation de la quantit de mouvement dun uide caractris par

e
e e
sa vitesse u et sa pression p. La PDE scrit sous forme dune divergence dun terme
e
tensoriel. Il est sous-entendu que la divergence est applique au deuxi`me indice
e
e
du tenseur. En outre, on dsigne par le symbole le produit tensoriel entre deux
e
vecteurs. Le tenseur I dsigne lidentit, et S le tenseur des taux de dformation,
e
e
e
dni comme la composante symtrique du tenseur des drives premi`res du champ
e
e
e e
e
de vitesse:
S=

1
2

u+

xii

Rsum en Franais
e
e
c

Figure 2: Le rseau D2Q9 pour la simulation de uides en deux dimensions.

e

Mthode de Boltzmann sur rseau

e
e
La mthode de Boltzmann sur rseau est prsente dans le chapitre 1.3 de la th`se,
e
e
e
e
e
ainsi que dans les rfrences indiques ` la n de ce chapitre. Un aperu rapide de
ee
e a
c
la mthode est reproduit ici, an dintroduire la notation.
e
La mthode LB consid`re une version discr`te de la fonction de distribution
e
e
e
f (x, c, t), dnie sur lespace des positions x et des vitesses c. Lors de la discrtisation,
e
e
lespace des vitesses c est omis en faveur dune augmentation du nombre de fonctions
de distribution, dont un nouvel exemplaire est introduit pour chaque individu dune
population de vecteurs de vitesse discr`te. Chacune de ces fonctions de distribution
e
est identie par un indice k: fk = fk (x, t). Les vecteurs vitesse utiliss dans le
e
e
cas discret reprsentent la relation spatiale entre une cellule du rseau et ses plus
e
e
proches voisins. En deux dimensions, on utilise par exemple frquemment le rseau
e
e
D2Q9 qui contient neuf vecteurs, un pour chacun des huit voisins, et un pour le
vecteur nul, reprsentant des particules au repos ne quittant pas la cellule. Ce
e
rseau est illustr sur la gure 2.
e
e
La dynamique discr`te dun mod`le LB exprime la propagation de ux de partice
e
ules dun noeud xij vers ses voisins. Lquation est crite dans un syst`me dunits
e
e
e
e
du rseau, dans lequel lespacement spatial entre deux noeuds voisins de la grille,
e
ainsi que lespacement temporel entre deux pas ditration, sont unitaires. La dye
namique LB sexprime donc par la formule
fk (xij + ck , t + 1) fk (xij , t) = k (xij , t),
o` le terme dsigne loprateur de collision entre particules. Dans le mod`le le
u
e
e
e
plus couramment utilis, connu sous le nom de BGK, la collision scrit comme une
e
e
relaxation vers un quilibre local:
e
eq
k = (fk fk ) ,

o` on a dni lquilibre, une version discr`te de la distribution de Maxwell-Boltzmann,

u
e
e
e
ne dpendant que des variables macroscopiques, de la mani`re suivante:
e
e
eq
eq
fk = fk (, u) = tk 1 +

1
1
1
c u + 4 (ck u)2 2 |u|2 .
2 k
cs
2cs
2cs

Le param`tre de relaxation peut tre formellement reli ` la viscosit du uide.

e
e
ea
e

Rsum en Franais
e
e
c

xiii

Les quantits macroscopiques peuvent tre calcules ` partir de fonctions de

e
e
e a
distribution en calculent les moments dordre zro, et du premier ordre. On trouve
e
ainsi
La densit: =
e

La vitesse: u = 1/

fk .
k ck fk .

Bien quen principe la mthode LB dcrive un uide compressible, on retrouve

e
e
lquation dun uide incompressible dans la limite des faibles nombres de Mach.
e
Dans ce cas, la relation entre la densit et la pression est dcrite par lquation pour
e
e
e
un gaz parfait:
p = c2 .
s

Dveloppement multi-chelle
e
e
Bien quon ait montr dans la section prcdente comment calculer les variables
e
e e
macroscopiques ` partir des fonctions de distribution, la procdure inverse na
a
e
jusqualors pas t aborde. An de connecter en dtail la dynamique LB avec
ee
e
e
une dynamique macroscopique, la limite du continu doit dabord tre value par
e
e
e
un dveloppement de Taylor du deuxi`me ordre, et une technique de sparation des
e
e
e
chelles temporelles et spatiales est utilise pour valuer la limite hydrodynamique
e
e
e
du mod`le. Le dveloppement de Taylor se rsume par lapproximation suivante:
e
e
e
fk (xij + ck , t + 1) fk (xij , t) (t +

ck ) fk +

1
2

2
t + 2t

ck +

ck

fk .

Dans le processus de sparation des chelles, les fonctions fk , ainsi que les drives
e
e
e e
spatiales et temporelles, sont dvelopps en fonction dun param`tre formel petit,
e
e
e
cest-`-dire
a
1:
(0)

(1)

fk = fk + fk + O( 2 ),
t = t1 + 2 t2 + O( 3 ) et

(1)
=
+ O( 2 ).
En rintroduisant ces approximations dans le mod`le BGK, on retrouve bien
e
e
lquation de Navier-Stokes, et on conclut que ce mod`le simule proprement la dye
e
namique dun uide. Ce rsultat est bien connu est reprsente un des piliers sur
e
e
lesquels repose la mthode. Ce nest par contre pas ce rsultat en lui-mme qui
e
e
e
nous intresse prsentement, mais plutt un rsultat intermdiaire, reliant les fonce
e
o
e
e
tions de distribution aux champs macroscopiques.
Le terme constant en est en eet quivalent ` la fonction dquilibre, et ne
e
a
e
dpend donc que de et de u:
e
(0)

eq
fk = fk (, u).

Le terme dordre 1 est li aux drives de la vitesse, ainsi que du carr de la

e
e e
e
vitesse:

tk
(1)
fk = 2
Qk : u
cs

1
ck : u u + 2 (ck )(Qk : u u) ,
2cs

xiv

Rsum en Franais
e
e
c

o` on a dni le tenseur symtrique

u
e
e
Qk = ck ck c2 I.
s
Il est important de remarquer que le dtail des fonctions de distribution nest
e
pas relevant dans le calcul de la limite hydrodynamique du mod`le. Limportance
e
revient plutt au moment du deuxi`me ordre, (1) , calcul selon la formule suivante:
o
e
e
(1)

(1) =

ck ck fk .
k

Au cours du dveloppement multi-chelle, on montre que ce tenseur est proportionnel

e
e
au tenseur des taux de dformation:
e
(1) =

2c2
s
S.

Lors du calcul du deuxi`me moment des fonctions de distribution, tous les termes
e
contenus dans la fonction de distribution nentrent pas en compte. En eet, les
termes dnis au carr de la vitesse sliminent pour des raisons de symtrie. Dans
e
e
e
e
le but dobtenir le bon tenseur (1) produisant une hydrodynamique approprie,
e
il est donc susant de considrer une version rduite des fonctions fk :
e
e

tk
(1)
fk = 2
Qk : u .
cs
Cette expression relie de mani`re simple le champs macroscopique aux valeurs des
e
fonctions de distribution.

Mod`le de Boltzmann sur rseau rgularis

e
e
e
e
Procd de rgularisation
e e
e
Dans le chapitre prcdent, on a introduit une version simplie des fonctions de dise e
e

tribution fk . Elle se dsigne par le terme fk et reprsente proprement les composants

e
e
hydrodynamiques du mod`le. Ces variables ont t relies par une relation simple
e
ee
e
aux drives du champ de vitesse. Pour lutilisation courante de cette expression, il
e e
serait par contre plus utile que les fonctions de distribution soient directement lies
e
aux fonctions macroscopiques simules, donc la densit , la vitesse u et le tenseur
e
e
des taux de dformation S. En eet, il a t montr que ces trois grandeurs peuvent
e
ee
e
tre values localement en chaque noeud du rseau LB, en calculant les moments
e
e
e
e
des fonctions de distribution.
Une telle relation se trouve en exploitant la symtrie du tenseur Q:
e

(1) = tk Q : u
fk
k
2
cs

tk
1
= 2 Qk :
( u) + ( u)T
cs
2
tk
= 2 Qk : S
cs
tk
Q : (1)
=
2c4 k
s

Rsum en Franais
e
e
c

xv

Cette relation est centrale dans pratiquement tous les sujets dvelopps dans la
e
e
th`se. Elle relie les champs macroscopiques aux valeurs des fonctions de distribution,
e
et peut donc tre exploite systmatiquement, ` chaque fois quune telle transfore
e
e
a
mation est ncessite. Pour souligner cette importance, la formule de rgularisation
e
e
e
est reproduite, en omettant cette fois-ci les calculs qui ont permis de la trouver:
tk
(1)
fk = 4 Qk : (1)
2cs

En suivant lenchainement suivant, la fonction de rgularisation est utilise pour

e
e
valuer les valeurs des fonctions de distribution ` partir des variables macroscopiques:
e
a

u
S

(1)
eq
fk

r
e
eq
(1)
fk egularis = fk + fk

Mod`le LB rgularis
e
e
e
La formule de rgularisation est utilise dans plusieurs domaines qui sont prsents
e
e
e
e
dans la th`se. A titre dexemple, lapplication probablement la plus importante est
e
explique dans ce qui suit: le mod`le de Boltzmann sur rseau modi, dont les
e
e
e
e
fonctions de distribution sont rgularises.
e
e
Lalgorithme pour la collision entre particules dans ce mod`le se dcompose en
e
e
trois tapes:
e
1. Calcul de , u et neq =

k ck

eq
ck (fk fk ).

r
e
2. Rgularisation fk fk egularis
e

3. Excution de la collision BGK sur les fonctions de distribution rgularises.

e
e
e
An dapprcier les qualits de ce nouveau mod`le, un ux bidimensionnel connu
e
e
e
sous le nom de ux de Kovasznay a t implant numriquement, une fois en utilisant
ee
e
e
la mthode BGK, et une fois en sappuyant sur le mod`le rgularis. Une solution
e
e e
e
analytique pour ce ux est connue, et on peut donc vrier la prcision de la solution
e
e
numrique, en la comparant ` lexpression analytique. Le rsultat de ce test a montr
e
a
e
e
que le mod`le rgularis est nettement plus prcis que le mod`le BGK, et que le gain
e e
e
e
e
en prcision augmente au fur et ` mesure que la grille est rane. Sur des grilles
e
a
e
de taille commune, on peut gagner un ordre de prcision en utilisant le mod`le
e
e
rgularis.
e
e
Un autre mod`le bidimensionnel a t implant, dsign par le nom de ux de
e
ee
e e
e
cavit. Il reprsente un uide contenu dans une boite quadratique, dont la paroi
e
e
suprieure est entraine ` une vitesse constante. Cette fois-ci, la stabilit numrique
e
e a
e
e
du mod`le a t teste, cest-`-dire, la valeur maximale du nombre de Reynolds
e
ee
e
a
pouvant tre atteinte avant que la simulation soit numriquement instable. Dans
e
e
les deux cas, BGK et rgularis, une relation linaire entre le nombre de Reynolds
e
e
e
maximal et le param`tre de rsolution N dune grille N N est tablie. La croissance
e
e
e

xvi

Rsum en Franais
e
e
c

linaire est par contre sept fois plus rapide dans le cas du mod`le BGK. Dans les
e
e
dtails, on trouve
e
Pour le mod`le BGK:
e
Remax = 0.4 + 1.0 N
Pour le mod`le rgularis: Remax = 16 + 7.4 N
e e
e
Il est donc tabli que pour ce type de probl`mes, le mod`le rgularis est nete
e
e e
e
tement plus prcis et plus stable que le mod`le BGK. Cet avantage sexplique par
e
e
lorigine de la mthode rgularise. Elle a t construite en imposant une sparation
e
e
e
ee
e
des chelles au mod`le original, le mod`le BGK. Uniquement les termes hydroe
e
e
dynamiques, cest-`-dire, le termes contribuant ` la dynamique de lquation de
a
a
e
Navier-Stokes, ont t prservs, et les termes dordre suprieur sont ngligs. Il
ee e
e
e
e
e
est vrai qu` cause de cette simplication, la mthode rgularise sest apauvrie.
a
e
e
e
On sattend en particulier ` ce quelle soit moins apte ` prsenter spontanment des
a
a e
e
ingrdients physiques dun uide complexe. Dans un but particulier par contre, celui
e
de la rsolution de lquation de Navier-Stokes, elle est plus concise, plus pratique
e
e
a
` manier, plus prcise et plus stable.
e

Introduction
Computational uid dynamics
The present thesis investigates a methodology in the domain of computational uid
dynamics (CFD), that is, the numerical simulation of uid ows. The core of this
research area is represented by a partial dierential equation (PDE) widely known
as the incompressible Navier-Stokes equation, which expresses a local conservation
law for the momentum in the system. It constitutes a largely simplied model for
a viscid, single component uid in an environment without density or temperature
variations nor any other specic physical ingredients. Although this equation only
partially addresses the complexity of most uids of interest in engineering applications, it is successfully applied in dierent areas for qualitative predictions of uid
ows. A case in point are optimization problems in aerodynamics, as for example
the design of a car body, or the creation of airfoils. It is obvious that the incompressible Navier-Stokes equation only roughly approximates the nature of compressible
air. Nevertheless, this equation leads to predictions of high quality that are successfully exploited, jointly with the results of physical experiments. Solving the
Navier-Stokes equation is actually a numerically very challenging task, in spite of
the apparent simplicity of the PDE. Much eort is therefore invested in developing
numerical approaches to the resolution of this equation, rather than adding even
more complexity by additional physical terms.
On the other hand, some research branches investigate so-called complex uids that exhibit extended physical properties, and to which the basic Navier-Stokes
equation is an insucient approximation. The classical approach in CFD to the
treatment of such uids is to describe the new physical properties in terms of transport phenomena related to a new observable, macroscopic property. A new PDE is
written down for the dynamics of this property, which then is solved by an appropriate numerical technique. The interaction between the dierent transport phenomena is handled via interaction terms between the respective PDEs. In a uid with
important temperature variations for example, a new observable property, the temperature, is introduced and its dynamics is described by a heat transport equation.
Furthermore, the momentum equation is extended to comprehend compressible uids. The temperature and the velocity are then linked through a convective term in
the heat equation and a buoyancy term in the momentum equation. An additional
equation of state links the temperature to the density. One can view the traditional
eld of CFD as a top-down approach to uid modeling, starting with a simple conservation law, and adding new equations on demand, based on both theoretical and
phenomenological considerations.

Introduction

Lattice Boltzmann method

The main topic of the thesis is a numerical approach known under the name of lattice
Boltzmann method. This method is relatively new and contrasts with the traditional
approach to CFD by adopting a bottom-up approach to uid modeling. To achieve
this, it describes the uid at a molecular level and proposes models for the collision
between molecules. The full continuum-level physics of the uid is implicitly contained in this model, and the conceptual work to be eectuated is opposite to the
one in the top-down approach of classical CFD. Indeed, the various physical ingredients contained in the model need to be identied one by one and properly explicited
in order to allow for a segregation between relevant and negligible properties. Based
on those considerations, the collision model is substantially simplied for the needs
of the numerical treatment. As such, the lattice Boltzmann (LB) method is well
adapted to the simulation of complex uids, as the relevant physical ingredients are
taken in charge very naturally by the method. How elegant a LB formulation of
a uid problem can be becomes apparent for example in the context of two-phase
ows. In traditional CFD, a PDE is written down for each of the two phases, as
well as for the interaction of the phases on their common interface. The implementation of such a model requires advanced software engineering techniques to track
the position of this interface and to implement its dynamics. Those problems dont
appear in the LB method, where the two phases can be represented by the same
uid model, and the interface between the phases is handled automatically, as a
natural ingredient of the model.
In spite of this apparent adequacy of the LB model for complex uids, the recent
historical development shows that much eort is invested in solving the raw NavierStokes equation with this method. This is reected by numerous publications and
technical reports, in which the LB method is used to solve benchmark applications
of classical CFD, and its relative eciency is compared to the one of other methods.
A reason for the success of the LB method in this domain can possibly be found
in the relative simplicity of its numerical implementation. As a matter of fact, the
physical principles that stand behind the method can be intimately connected with
classical programming paradigms of high performance computing. For example, the
fact that collisions between particles are local in their nature encourages the use of
parallel architectures in which spatial subdomains of the problem are solved locally
on dierent computational nodes.

The LB method applied to CFD

When the LB approach to physical modeling is imported in the domain of classical
CFD, two very dierent worlds with dierent traditions and convictions are confronted. This leads to numerous misunderstandings, on the one hand from the CFD
community that fails to recognize the potential of the LB method by its right value,
and on the other hand from the LB community that has diculties to recognize the
real challenges of CFD, or the areas in which the LB method can prove competitive
with other approaches. One main obstacle to a mutual exchange between those two
methods stems from the fact that they describe the uid via a dierent choice of vari-

Introduction

ables. The variables of concern in CFD are the macroscopic, observable quantities.
In the incompressible Navier-Stokes equation, the degrees of freedom are dened
by the space and time dependent velocity u = u(r, t) and the pressure p = p(r, t).
Classical CFD solvers most often calculate directly a discrete representation of these
variables, that is, they store the value adopted by the variables on the nodes of a
given numerical grid. The LB method on the other hand simulates the dynamics of
a so-called particle distribution function, a quantity derived from the kinetic theory
of gases that represents a probability density for the presence of particles in phase
space. Although it is well known how the macroscopic variables are computed from
a given state of the particle distribution function, the reverse procedure is much
more contrived. It is actually not always clear how a given conguration of the
macroscopic variables can be used to initialize the distribution function. This fact
is rarely recognized as a problem by a community that is used to thinking exclusively in terms of kinetic quantities. It does however generate many complications
whenever macroscopic variables are explicitly evoked in the simulation, or when the
macroscopic length and time scales of the model need to be adapted.

How to treat varying space and time scales?

As an example, the diculty becomes apparent when a grid renement technique
is used, in which the computational domain is partitioned into subdomains that
are represented by grids with a varying degree of resolution. The purpose of this
procedure is to adapt the structure of the grid to the geometry of the problem,
and thus to gain computational eciency. During the simulation, it is necessary to
communicate successively computed values from one grid to another. This is easily
performed in classical uid solvers. Indeed, they store the value of dimensionless
macroscopic variables that do not depend on a given resolution of the grid and can
therefore communicate those values without the need for conversions. On the other
hand, the particle distribution function used in LB methods is not dimensionless
with respect to a choice of macroscopic units. It must therefore be rescaled during
the transfer from a grid to another, because its value depends on the specic grid
spacing. The rules that dictate this rescaling are in no way straightforward, and
much care must be used to get them right. A certain confusion reigns consequently in
the literature on this point, and numerous dierent approaches have been proposed,
based on substantially dierent points of view. Grid renement is just a case in
point to illustrate the encountered diculties, but similar problems are encountered
with essentially all practical implementation issues of LB, such as the denition of
initial conditions or boundary conditions, the implementation of an adaptive time
interval, the coupling of LB methods with other techniques, and many more.

Multiscale expansion of lattice Boltzmann

The value of the distribution function can be formally related to the macroscopic
variables through an innite series via the so-called Chapman-Enskog multi-scale
expansion. This expansion is commonly used to verify that a given LB method
eectively solves the macroscopic PDEs asymptotically. To do this, the innite

Introduction

series is truncated. The remaining lower order terms, sucient to recovering the
Navier-Stokes equation, are called hydrodynamic terms. This technique is however
not so frequently used to actually develop new methods in the domain of LB. The
point of view adopted in the present thesis is that practically all implementation
issues in the LB method can be addressed properly by systematically taking into
account the results of the Chapman-Enskog analysis. In this approach, intuitive
reasoning is underlined by formal proofs, and a uniform strategy is devised to solve
dierent problems in a common framework.

Content of the thesis

In the rst part of the thesis, various theoretical results from the literature are
gathered to show that a LB implementation eectively produces results equivalent to
those of a conventional CFD solver. The necessity to include some less known terms
into the model is highlighted, and the limitations of the approach are emphasized.
This part includes also a novel contribution that shows how to properly implement
the application of an external force term in a compressible uid model with adaptable
bulk viscosity. Given that the theoretical developments can be somewhat hard to
follow, they are preceded, for illustration purposes, by similar-looking but simpler
discussions of a LB model for advection-diusion problems. On the way, it is shown
that this advection-diusion model is treated erroneously in the literature, and a
workaround is suggested.
In the second part, a modied LB model is introduced under the name of regularized lattice Boltzmann model. In this model, the original LB distribution function
is replaced by a regularized function, inspired by the hydrodynamic terms found
in the Chapman-Enskog expansion. The specicity of the regularized distribution
function is that it can be fully expressed in terms of macroscopic, observable variables. It can therefore be decomposed into components that are dimensionless in
a system of macroscopic variables, and is rescaled in a controllable way when the
parameters of the space or time discretization vary. Furthermore, as it is shown on
selected CFD applications, the resulting numerical model is more stable and accurate than the original one. The ability to rescale the variables of the regularized
model in a controllable way is exploited in the subsequent chapters. Through a unied approach, it is explained how to implement grid renement, an adaptive time
interval, boundary conditions, initial conditions, and a coupling of the LB model
with another CFD tool. In the last chapter nally, a new point of view is adopted,
and the focus turns back to the original, non-regularized LB method. It is discussed
how the non-hydrodynamical terms can be purposely exploited to create a model
of uid turbulence. The theoretical considerations are accompanied by results of
direct numerical simulations of a turbulent uid.

Chapter 1
Computational uid dynamics and lattice
Boltzmann

1.1

Vector and tensor notation

Throughout this thesis, vectors in the 2D or 3D Euclidean space are noted with
an arrow on top of them as in a, and higher order tensors by a bold face notation
as in A. Their indexes are labeled by Greek letters, as in a and in A . Two
types of notation are used to express algebraic operations on vectors and tensors.
In the index notations, all indexes are explicited, but the sums are skipped: it
is implicitly understood that they must be taken whenever an index is repeated
twice in the same term. With this notation, a scalar product between a and b is
written as = a b . In a full vector notation, the indexes are skipped, and the
scalar product is written as = a b. The shorter and more readable vector/tensor
notation is generally preferred. The more explicit index notation is however still
used in complicated expressions involving higher order tensors, to make sure they
are properly interpreted. The following list shows through several examples how
operations between vectors and tensors are denoted:
Vector/tensor notation Index notation
Scalar product
=ab
= a b
Tensor contraction
=A:B
= A B
Dyadic vector product
A = ab
A = a b

Gradient
a=
a =

Divergence of vector
= a
= a

Divergence of tensor (order 2) a = A

a = A

Divergence of tensor (order 3) A = T

A = T
Another vector space, the space Rq of the particle distribution functions, sometimes appears in the thesis. In order to clearly distinguish between the two types
of vectors, the indexes of vectors in Rq are labeled by a Latin index, and none of
the shorthand notations introduced above is used for them. Occasionally, a scalar
product in Rq is executed, for which a bracket notation < | > is used. Furthermore, Section 3.3.1 makes a short exception and uses a full vector notations for the
vectors in Rq , in order to underline the origin of the computations as linear algebra

Chapter 1

expressions.

1.2
1.2.1

Macroscopic equations
Governing equations

Computational Fluid Dynamics (CFD) provides a qualitative, and in some cases

even quantitative, prediction of uid ows by means of mathematical and numerical
tools. The chain leading to a full numerical model includes the following steps:
The mathematical model consists most often of a set of partial dierential equations (PDEs). They describe the evolution of the observable, physical quantities in the ow and are based on phenomenological and/or theoretical considerations. As an alternative, a simplied microscopic model can be directly
simulated on a computer, without referring explicitly to a PDE.
The Numerical method provides a way to nd approximate solutions of the
PDE, or to solve the dynamics of a microscopic model, with the help of computer simulations.
Some Software tools , consisting of specialized algorithms, are used to implement
parts of the numerical model. Pre- and postprocessing utilities are required to
handle the setup of a simulation and the analysis of the simulated data.
The present section introduces the most commonly used mathematical model,
which consists of a set of PDEs. It describes the uid at a continuum level, at
which its particulate nature can be neglected. Instead, one takes interest in the
dynamics of a reduced set of macroscopic variables that describe the state of the
uid suciently well. Whenever the uid exhibits a physical behavior that cannot
be described by the macroscopic model, new macroscopic variables and PDEs are
introduced to adapt the model. This technique is opposed to an approach in which
one tries to obtain an understanding of the microscopic nature of the uid. In that
case, the observed macroscopic physics of the uid is an emergent property that
arises naturally from the microscopic model.
The state of a simple uid is described by the following macroscopic variables:
the
the
the
the
the

uid density
ow velocity
pressure
energy
temperature

,
u,
p,
E and
T.

The governing equations for those variables can be derived by methods of statistical
physics from the equations of motion of the microscopic model. Alternatively, they
can be interpreted as conservation laws for the concerned variables. To reect this
fact, they are written in a divergence form:

U
+ F = Q.
t

(1.1)

Computational uid dynamics and lattice Boltzmann

The most commonly used equations express the conservation of

mass:

momentum: t (u)
energy:

t E

+
+
+

(u)

= 0,

(1.2)

(uu + pI )

= g

and (1.3)

(E + p)u T (u )T

= (q + g u).
(1.4)

The parameter stands for the thermal conductivity, and q for internal heat sources.
The tensor represents the deviatoric stresses in the uid. For a Newtonian uid
it has, by denition, the following form:

= ( u) I + 2 S,

(1.5)

where the strain rate tensor S is dened as

1
S = ( u + ( u)T ).
2

(1.6)

The parameters and stand for the dynamic shear and bulk viscosity respectively.
Based on a microscopic model, in which the molecules of a gas are represented by
a single rigid sphere, one can derive the following theoretical relationship between
those parameters (see e.g. [1]):
2
(1.7)
= .
3
In real gases and liquids however, the actual value of can vary substantially from
this prediction. Equation (1.3) is commonly known as the compressible NavierStokes equation.
Eqs. (1.2,1.3,1.4) represent a system of d + 2 equations for d + 4 unknowns, in a
d-dimensional space. Is must be closed by additional equations of state that relate
the pressure with the temperature and the temperature with the energy.

1.2.2

Isothermal and incompressible ows

In an isothermal ow, eects of the temperature are neglected, and Eq. 1.4 is not
taken into account. In that case, an equation of state must be devised that relates
the pressure and the density in the uid. In the following, an equation for ideal
gases is used, in which the pressure scales linearly with the density:
p = c2 ,
s

(1.8)

where cs is the speed of sound.

In an incompressible uid, the density takes a constant value = 0 which does
not vary in time and space. The conservation law for the mass is simplied, and
states that the velocity eld is divergence-less (solenoidal):

u = 0.

(1.9)

Chapter 1

Eq. (1.9) is often used as a denition for uid incompressibility. The conservation law
for the momentum is also simplied and leads to the incompressible Navier-Stokes
equation, often called Navier-Stokes equation for short:
t u + (u

)u =

1
p+
0

u.

(1.10)

In order to nd the value of p, one takes the divergence of Eq. (1.10) and makes use
of Eq. (1.9). This yields the Poisson equation:
2

p = 0 ( u) : ( u)T .

(1.11)

This time-independent equation replaces Eq. (1.9). When the evolution of an incompressible ow is computed numerically, Eq. (1.11) needs to be solved by an iterative
procedure at every discrete time step. It adjusts the value of the pressure in such a
way that the velocity eld remains divergence-less during its time evolution.

1.2.3

Dimensionless formulation

Before the governing equations for uid motion can be solved on a computer, one
needs to get rid of the physical units of the macroscopic variables. This leads to a
set of PDEs that act on dimensionless variables, and whose properties are tuned
via generic, dimensionless parameters. For the purpose of illustration, Eqs. (1.9)
and (1.10) are now cast in a dimensionless form. For this, a length scale l0 and a
time scale t0 are introduced that are representative for the ow conguration. The
length l0 could for example stand for the size of an obstacle which is immersed in
the uid, and t0 could be the time needed by a passive scalar in the uid to travel
a distance l0 . The physical variables for time and position, t and r, are replaced
by their dimensionless counterpart t = t/t0 and r = r/l0 . In the same manner, a

2
= 1/l0 and p = 0 l0 /t2 p is
change of variables for u = l0 /t0 u , t = 1/t0 t ,
0
performed. This leads to the following dimensionless Navier-Stokes equation:

t u + (u

1
)u = p +
Re

u = 0,

and

(1.12)
(1.13)

2
where Re = l0 /(t0 ) is the dimensionless Reynolds number. Two ows that obey
the same Navier-Stokes equation and which have the same Reynolds number are
equivalent. In this way, it is possible to solve the dimensionless equations on a
computer, and to relate the results with a physical ow that possesses the same
Reynolds number and the same ow geometry.

1.2.4

Further reading

A large number of books discuss the numerical resolution of PDEs and the topic of
computational uid dynamics. Only a few are mentioned here to guide the interested
reader.

Computational uid dynamics and lattice Boltzmann

A well written and theoretically sound introduction to the numerical treatment

of PDEs can be found in Ref. [2]. A general introduction to the problems of computational uid dynamics is available on the web presented as a one-semester introduction course [3]. An undergraduate-level, practically oriented introduction
to computational uid dynamics with a specic emphasis on the nite dierence
method is presented in Ref. [4]. Reference [5] nally gives an extended overview
on various methods of computational uid dynamics in general, and Ref. [6] on the
nite element method in particular.

1.3
1.3.1

Lattice Boltzmann method

Discrete variables

The problems of uid dynamics possess an innite number of degrees of freedom,

because the considered quantities (velocity, pressure, etc.) are not just single values,
but spatially extended elds. Computers can however only handle a nite amount
of variables, and represent them with a nite accuracy. The spatially extended elds
must therefore be replaced by a nite set of scalar values that are appropriate for the
numerical investigation of the problem. This step is called the space discretization
of the problem. To keep the discussion simple, it is restricted to problems whose
domain of denition is conned within a regular box of nite extent. A possible
way to discretize this space is to partition it into cells. In that case, a numerical
representation of the uid consists of a set of numbers, of which each stands for the
average value of a uid variable over a cell. Another approach to space discretization consists in replacing spatially extended elds by their value on a chosen nite
population of space points. The lattice Boltzmann method, as it is presented here,
adopts the latter point of view. The space is hence discretized on a regular lattice
with xed grid spacing. The present discussion concentrates on a cubic subsample
of a three-dimensional system whose size is l0 l0 l0 (l0 is a characteristic length of
the system, introduced in Section 1.2.3). The coordinate system is dened in such
a way that one corner of the cube is situated at the origin, and the opposite corner
at the position l0 e0 + l0 e1 + l0 e2 . The cube is represented numerically by N 3 points
rijk , with i, j, k = 0 N 1, situated at the position rijk = ie0 x + je1 x + ke2 x .
Those positions are called the lattice sites. The lattice spacing x = 1/(N 1) is
equal to the distance between two neighboring lattice sites. A scalar eld (r, t) is
represented by N 3 values ijk (t) that are numerical approximations of (rijk , t).
The time axis is discretized in the same manner by a set of equidistributed nite time steps tn = n t . In the numerical representation that is adopted here, a
time-dependent variable f (t) is replaced by a numerical approximation fn f (ti ).
Finally, to combine the notation for space and time discretization, the somewhat
obscure index notation is again skipped for the more readable notation with parentheses. That is, the expression f (rijk , tn ) is used to represent the numerical approximation (and not the exact value) of f at the position rijk and at the time step
tn .

10

1.3.2

Chapter 1

Lattice Boltzmann model

Each node of a lattice Boltzmann simulation holds a set of q variables fi , i =

0 q 1, called the particle distribution functions. Each of those functions is
responsible for carrying information from a node to one of its neighbors. The position
of this neighbor is dened by a vector ci , which points from a lattice node to the
corresponding lattice node. It is characteristic for the lattice and does not depend
on space or time. The following general rule for the evolution of a LB model explicits
the role of the vector ci :
fi (r + ci , t + 1) fi (r, t) = i .

(1.14)

This equation has been written in lattice units, in which the spacing between two
adjacent lattice nodes, as well as the time interval from one iteration to the next,
are unitary. In this way, a close relationship between the theory and the application
is kept, as it is common to write numerical implementations of LB models in lattice
units. This approach is opposed to the one of other numerical models that are most
often implemented in a system of dimensionless units introduced in Section 1.2.3).
The term i on the right hand side of Eq. 1.14 is the collision operator that describes how the q values fi dened on the same node at given time step interact.
The macroscopic variables, the density and the velocity u are dened locally as
moments of the distribution functions:
q1

fi

and

(1.15)

i=0
q1

u=

ci fi .

(1.16)

i=0

In the following, whenever a sum over all distribution functions is taken, the range
of the variable i is not explicited any more to keep the notation short: the sum over
the q elements of a variable Ei , q1 Ei , is simply written as i Ei .
i=0
The most commonly used collision operator, known under the name of BGK
collision and discussed extensively in Ref. [7], implements a relaxation dynamics
towards a local equilibrium with a relaxation parameter :
BGK = (fi fieq ).
i

(1.17)

The local equilibrium is dened as

(eq)

fi

= ti 1 +

1
1
c u + 4 Qi : uu ,
2 i
cs
2cs

(1.18)

and the tensor Q is dened as follows:

Qi = ci ci c2 I.
s

(1.19)

The constant cs is the speed of sound of the model. This parameter, as well as the
q parameters ti are lattice constants. Conceptually speaking, there is some freedom
in the choice of the constant cs . All that needs to be done is to modify the value of

Computational uid dynamics and lattice Boltzmann

11

the rest particle weight t0 correspondingly, in order to recover the lattice symmetries
described in Section 2.1. In practice, a value of c2 = 1/3 is found to be numerically
s
most stable, and this choice is therefore most commonly adopted. In conclusion,
there are two ingredients to the denition of a lattice Boltzmann model: on the one
hand the details of the collision operator i , and on the other hand the structure
of the lattice, dened by the constants q, ci , cs and ti . Some of the most commonly
used lattice structures are summarized in Section 1.3.3.
The numerical model dened by Eqs.(1.14,1.18) solves asymptotically the equations of motion for a compressible uid, when the discrete steps t and x , as well
as the relative velocity u/cs are small enough. This fact is proved with some reservations in Section 2.3. It is also shown that the relaxation parameter is directly
related to the dynamic shear viscosity of the uid via Eq. 2.73:
= c2
s

1 1

(1.20)

The formulas of the lattice Boltzmann method can be derived theoretically via
dierent approaches. One approach, advocated for example in Ref. [8], views the
LB method as a continuous version of a microscopic model known as a Cellular Automata. In another approach, described for example in Ref. [9], the dynamics of the
LB method is derived from the continuum Boltzmann equation. This equation considers the movement of molecules in a gas and describes their behavior statistically
at a continuum level (see e.g. Refs. [10, 1]). For this reason, the theory behind the
Boltzmann equation is often called kinetic theory. By extension, the LB method is
sometimes called a lattice kinetic scheme, and the quantities fi used in the method
are referred to as kinetic variables.
This concludes the presentation of the lattice Boltzmann method. An overview
of lattice structures that can be used in common with the BGK model is presented
in the next section. A discussion of how to choose the system of units and relate the
lattice variables to macroscopic ones is found in Section 2.4. For a more detailed
discussion of implementation issues related to the LB method, the reader is referred
to Refs. [11, 12].

1.3.3

Lattice structures

A lattice structure with q lattice directions, dened on a d-dimensional space, is commonly identied by the name DdQq lattice. A few 2D and 3D lattice structures
that are commonly used for isothermal ows are listed in this section. For those
ows, it is sucient to consider structures in which the lattice vectors ci point to
the neighbors included in the immediate neighborhood, so-called nearest neighbors.
The lattice weights ti are used to account for the fact that not all lattice vectors
have the same length. Each lattice structure has three types of lattice weights, the
weight t0 corresponding to the zero velocity c0 = 0, the weight ts corresponding to
the short velocities and the weight tl corresponding to the long velocities. In 2D for
example, the short velocities, parallel to the lattice edges are of length 1, and the

long velocities following diagonal directions are of length 2.

12

Chapter 1

D2Q9 lattice
c2 =
s
t0 =

1
3
4
9

ts =

1
9

tl =

1
36

c0 = (0, 0)
c1 = (1, 1) c2 = (1, 0) c3 = (1, 1) c4 = (0, 1)
c5 = (1, 1) c6 = (1, 0)
c7 = (1, 1)
c8 = (0, 1)

D3Q15 lattice
c2 =
s
t0 =

1
3
2
9

ts =

1
9

tl =

c0 = (0, 0, 0)
c1 = (1, 0, 0)
c4 = (1, 1, 1)
c8 = (1, 0, 0)
c11 = (1, 1, 1)

1
72

c2 = (0, 1, 0)
c5 = (1, 1, 1)
c9 = (0, 1, 0)
c12 = (1, 1, 1)

c3 = (0, 0, 1)
c6 = (1, 1, 1) c7 = (1, 1, 1)
c10 = (0, 0, 1)
c13 = (1, 1, 1) c14 = (1, 1, 1)

D3Q19 lattice
c2 =
s
t0 =

1
3
1
3

ts =

1
18

c0 = (0, 0, 0)
c1 = (1, 0, 0)
c4 = (1, 1, 0)
c7 = (1, 0, 1)
c10 = (1, 0, 0)
c13 = (1, 1, 0)
c16 = (1, 0, 1)

tl =

1
36

c2 = (0, 1, 0)
c5 = (1, 1, 0)
c8 = (0, 1, 1)
c11 = (0, 1, 0)
c14 = (1, 1, 0)
c17 = (0, 1, 1)

c3 = (0, 0, 1)
c6 = (1, 0, 1)
c9 = (0, 1, 1)
c12 = (0, 0, 1)
c15 = (1, 0, 1)
c18 = (0, 1, 1)

D3Q27 lattice
The D3Q27 lattice uses lattice vectors of three dierent lengths: the short ones,
such as (1, 0, 0) with weight ts , the medium ones like (1, 1, 0) with weight tm and
the long ones like (1, 1, 1) with weight tl .
c2 =
s
t0 =

1
3
8
27

ts =

2
27

tm =

1
54

tl =

1
216

Computational uid dynamics and lattice Boltzmann

c0 = (0, 0, 0)
c1 = (1, 0, 0)
c4 = (1, 1, 0)
c7 = (1, 0, 1)
c10 = (1, 1, 1)
c14 = (1, 0, 0)
c17 = (1, 1, 0)
c20 = (1, 0, 1)
c23 = (1, 1, 1)

1.3.4

c2 = (0, 1, 0)
c5 = (1, 1, 0)
c8 = (0, 1, 1)
c11 = (1, 1, 1)
c15 = (0, 1, 0)
c18 = (1, 1, 0)
c21 = (0, 1, 1)
c24 = (1, 1, 1)

13

c3 = (0, 0, 1)
c6 = (1, 0, 1)
c9 = (0, 1, 1)
c12 = (1, 1, 1) c13 = (1, 1, 1)
c16 = (0, 0, 1)
c19 = (1, 0, 1)
c22 = (0, 1, 1)
c25 = (1, 1, 1)
c26 = (1, 1, 1)

Further reading

It is more dicult to nd well written textbooks on the lattice Boltzmann method

than on other topics of computational uid dynamics. A few are however worth
mentioning and are listed in the following.
Some general information on the LB method, links to further literature and
an open source LB simulation code are found on the web page of the OpenLB
project [13]. As interesting introductory textbooks with some theoretical backgrounds Refs. [8, 14, 1] should be mentioned. The Refs. [11, 12] complete the
list with a more practically oriented approach and a discussion of implementation
details. An overview of LB methods, including multi relaxation time models, is
contained in the review paper [9], of which a reprint can be downloaded from the
internet.

Chapter 2
Multiscale Chapman-Enskog analysis

2.1

Series expansion with scale separation

In this section, the generic evolution equation of LB methods is inspected by means

of a truncated Taylor expansion and a multi-scale analysis. The results of this study
can be applied to show that a specic LB method solves asymptotically the dynamics
of the desired macroscopic PDE. In the present thesis, they are furthermore used
systematically, whenever a LB model needs to be adapted to a specic situation.
All developments contained in this section are largely inspired from corresponding
sections in Ref. [8].

2.1.1

Lattice symmetries

Not all lattice types are adequate for the execution of a lattice Boltzmann simulation.
In order for the simulation to yield the desired asymptotic PDE, the lattice must
verify a set of symmetry conditions. The BKG model for a uid for example requires
that there exists a constant cs and a set of weights ti for the lattice velocities ci such
that the following relations are veried:
(a)

ti = 1

ti ci ci = c2
s

(c)

(e)

ci ci ci ci =

(2.1)

c4 ( + + )
s
(b)

ti ci = 0
i

(d)

ti ci ci ci = 0
i

(f )

ti ci ci ci ci ci = 0.
i

The value of cs may vary from one lattice to another. In the following, this parameter
will be shown to be equal to the speed of sound in a simulation. The symmetries
displayed in Eq. (2.1) are present in all types of lattices listed in Section 1.3.3.
Note 2.1: Extended lattice structures

The lattice symmetries shown in

Eq. (2.1) are milestones to the theoretical developments that follow and are fundamental to the presented BGK lattice Boltzmann models. There exist however
specic lattice Boltzmann models that work on lattices with weaker symmetries.
A case in point is the model presented in Ref. [15] that works ne on a 13-velocity

16

Chapter 2

3D lattice (D3Q13), although this lattice lacks sucient symmetries to be used

with BKG dynamics. On the other hand, there exist models that require additional symmetry properties to those of Eq. (2.1). So-called thermal LB models
for example use an extended set of lattice velocities that include interaction with
next-to-nearest neighbors.
Equations (2.1) are sometimes viewed as orthogonality properties between lattice
vectors dened on the vector space Rq . The rst of those lattice vectors is called C 0
and is dened through
Ci0 = 1 i = 0 q.
(2.2)
1
It follows a set of d vectors C dened as
1
Ci, = ci

for = 0 d 1 and i = 0 q.

(2.3)

2
Finally, a set of d2 lattice vectors C is introduced that obey the relation
2
Ci, = Qi

for , = 0 d 1 and i = 0 q,

(2.4)

where the tensor Q is dened in Eq. (1.19). It is pointed out that by symmetry of
2
2
2
the tensor Q, the vectors Ci verify the relation Ci = Ci . Equation (2.4) denes
therefore only d(d + 1)/2 independent vectors, instead of d2 . A scalar product
between two vectors a and b in Rq is nally dened as follows:
a|b =

ti ai bi .

(2.5)

With this, it is easily concluded from Eq. (2.1) that the lattice vectors dened above
are orthogonal to each other, but not necessarily unitary:
1
C 0 |C = 0,

(2.6)

2
C 0 |C = 0,

(2.7)

1 2
C |C
0 0

= 0,

(2.8)

= 1,

(2.9)

C |C

C 1 |C 1 = c2 I
s
2
2
C |C

c4
s

and

( + ).

(2.10)
(2.11)

The fact that Eq. (2.11) expresses an orthogonality relation between all vectors of
the family C 2 becomes clear when the tensor resulting from a scalar product between
two of those vectors is contracted with an arbitrary tensor T :
2
2
C |C T = c4 (T + T )
s

(2.12)

The q-dimensional vector space introduced in the previous paragraph is of both

theoretical and technical interest. To see this, the unweighted particle distribution
functions
gi = fi /ti
(2.13)

Multiscale Chapman-Enskog analysis

17

are introduced. They are interpreted as vectors in Rq , and the index i is consequently skipped. All hydrodynamic variables in a LB model, which were previously
introduced as moments of the distribution functions, are now reinterpreted as projections of g on the lattice vectors. The mass density for example is computed as
= C 0 |g , and the momentum as u = C 1 |g . Most of the algebra shown in the
next section can be viewed as the computation of similar projections. This point of
view is useful, because vanishing terms in the computation are at once identied as
projections between orthogonal vectors.
Note 2.2: Hermite polynomials

The vectors C 0 , C 1 and C 2 introduced in this

section can be viewed as tensor Hermite polynomials. In some theoretical approaches to the LB method, the equilibrium distribution function in Eq. (1.18)
is derived by means of a truncated expansion of the velocity in tensor Hermite
polynomials. The key ideas of this approach are developed in Ref. [16].

2.1.2

Multi-scale expansion

LB models with a generic collision operator are dened by the evolution equation (1.14) on page 10. The left-hand side of this equation can be expanded in a
second-order Taylor series as follows:
i = fi (r + ci , t + 1) fi (r, t)

1 2
(t + ci ) fi + (t + 2t ci +
: ci ci ) fi .
(2.14)
2
In order to relate the LB equation with a macroscopic PDE, it is necessary to
formally separate dierent time scales. In this way, physical phenomena occurring
at dierent scales are discussed separately and contribute individually to the nal
equations of motion. To obtain this, the time derivative is expanded in terms of a
formal, small parameter :
t = t1 +

t2 + O( 3 ).

(2.15)

The idea behind this expansion is that all involved terms (t1 and t2 ) are of the
same order of magnitude, and the smallness is introduced via the parameter . The
space derivative is not expanded beyond its rst-order term:

(1)
=
+ O( 2 ).
(2.16)
The particle distribution functions are likewise expanded, starting with a zerothorder contribution f (0) :
(0)

fi = fi

(1)

+ fi

2 (2)
fi

+ O( 3 ).

(2.17)

It is concluded from Eq. (2.14) that the collision term i does not contain constant
(0)
contributions with respect to the parameter : i = 0. It is therefore expanded as
follows, starting with the O( ) term:
(1)

i = i +

(2)

i + O( 3 ).

(2.18)

18

Chapter 2

This expansion, up to a second order accuracy with respect to , appears to be sucient to nd the Navier-Stokes equation. However, higher order terms, known under
the name of Burnett terms are present in the LB dynamics and can be computed.
In the following, all developments are written down separately for the O( ) and the
O( 2 ) scales, and all terms of the order O( 3 ) or higher are neglected.
The scale separated version of Eq. (2.14) reads:
(1)

(2)

i + 2 i = ( t1 + 2 t2 +

1 ci +

1 2 2 2
+ t1
2 t2

1 ci +

1 2
(0)
(1)
1 1 : ci ci )(fi + fi ).
2
(2.19)

Note 2.3: Knudsen number The parameter is often identied as the Knudsen number, that is, the ratio f /l0 between the mean free path of a gas molecule
and a macroscopic length scale. This terminology is motivated by a dimensional
analysis of the Boltzmann equation. Indeed, when all terms of this equation are
made dimensionless with respect to microscopic scales, the collision term is left
with a trailing factor f /l0 (see e.g. Ref. [14]). As a result of the truncated
multi-scale expansion, the LB equation (as well as the Navier-Stokes equation)
is valid only for low Knudsen number ows. It does not apply well for example
to dilute gases or microuids.

2.1.3

Conservation laws

The macroscopic variables of the ow are dened as moments of the particle distribution functions, that is, as the projection of the unweighted distribution functions
gi on the base vectors C 0 , C 1 and C 2 . This leads to the denition of the zerothorder scalar moment , the rst-order vector moment j and the second order tensor
moment :
=

fi ,

(2.20)

j=

ci fi

and

(2.21)

Qi fi .

(2.22)

The dynamics of a ow can be expressed by means of conservation laws, applied to

some of those moments. A global conservation of a macroscopic quantity is expressed
locally by a collision invariant. As an example, conservation of mass in a uid is
enforced by a local conservation of mass during the collision between particles. At
the level of the LB equations, a collision invariance means that the projection of the
unweighted collision operator i /ti on the corresponding base vector vanishes.
Zeroth order moment balance. Conservation of the zeroth order moment is
ensured by the collision invariant
i = 0,
i

(2.23)

Multiscale Chapman-Enskog analysis

19

and by the requirement that rst-order contributions to vanish:

=

(0)

fi =

fi .

(2.24)

This latter condition can be understood intuitively as follows: the conserved variables are moments of f (0) only, whereas the collision operator acts on f (1) , which
leaves the conserved variables unchanged during the collision. Alternatively, Eqs. (2.23)
and (2.24) can simply be taken as technical requirements that serve the nal aim,
namely to nd a LB scheme leading to the Navier-Stokes equation. In Sections (2.2)
and (2.3), LB collision operators are explicitly such as to respect the conservation
laws.
Expanding Eq. (2.23) on two dierent scales and 2 , and using Eq. (2.19), yields
(1)

(0)

i = t1
i

fi

(0)

ci fi

= t1 +

j (0) = 0

(2.25)

and
(1)

(2)

fi

i = t1

(0)

fi

+ t2

+ t1

(1)

ci fi
i

(0)

ci fi

= t2 +

1
1 1 :
2

1 2
+ t1
2

(0)

fi
i

(0)

(Qi + c2 I)fi
s

(2.26)

1 2
1
1
(0)
(1)
+ t1 + t1 1 j (0) +
+ c2
1j
1 1 :
2
2
2 s

2
1 .

Equations (2.25) and (2.26) are combined to eliminate the second-order time deriva2
tive t1 :
t2 +

1
1
1
(0)
+ c2
j (1) + t1 1 j (0) +
1 1 :
2
2
2 s

2
1

= 0.

(2.27)

First order moment conservation. At the rst order moment, a source term
F = F (1)

(2.28)

is added. It will be identied in the following as an external force acting on a uid.

Depending on the specic requirements for the numerical model, a corresponding
source term can of course also be added to the zeroth order moment. This could for
example lead to an advection-diusion equation with source terms. The rst order
collision invariant reads
ci i = F .
(2.29)
i

In analogy with the zeroth-order moment, the rst-order moment is required to

depend only on zeroth-order particle distribution functions. For reasons that become
clear shortly, a correction term due to the source F must be added:
j=

(0)

ci fi =
i

ci fi
i

F
.
2

(2.30)

20

Chapter 2

The O( ) scale contribution to Eq. (2.29) reads

(1)
(0)
ci i = t1
ci fi + 1
i

= t1 j

(0)

(0)

(Qi + c2 I)fi
s

(2.31)

(0)

c2 1
s

=F

(1)

Two O( 2 ) terms appearing in Eq. (2.27) can now be evaluated and inserted into
this equation:

(2.32)
t1 1 j (0) = 1 F (1) 1 1 : (0) c1 2
by (2.31), and
s 1

(1)

1
= 1F
(2.33)
by (2.30).
1j
2
Eqs. (2.27), (2.32) and (2.33) are combined to obtain
t2 = 0,

(2.34)

which completes Eq. (2.25) to form the continuity equation:

t + j (0) = 0.

(2.35)

It is now clear that the particular form of Eq. (2.30), including a force term, was
required to cancel erroneous force contributions to the continuity equation.
Contributions of the order O( 2 ) to Eq. (2.29) yield:

(1)
(0)
(1)
(2)
(Qi + c2 I)fi
ci fi + 1
ci fi + t2
ci i = t1
s
1 2
+ t1
2

(0)

ci fi

+ t1

(0)

(Qi + c2 I)fi
s

1
(0)
1 1 : R
2

(2.36)

1 2 (0)
= t1 j (1) + t2 j (0) + 1 (1) + t1 j + t1 1 (0)
2

1
(0)
+ c2 t1 1 +
= 0,
1 1 : R
s
2
where the tensor R with components R = i ci ci ci fi has been introduced.
The second order time derivative is canceled with the help of Eq. (2.31),

2
(2.37)
t1 j (0) = t1 F (1) 1 (0) c2 1 ,
s
i

which simplies Eq. (2.36) as follows:

1
1
1
(0)
(1) + t1 1 (0) + c2 t1 1 +
= 0.
(2.38)
1 1 : R
s
2
2
2
Finally, the O( ) term in Eq. (2.31) and the O( 2 ) term in Eq. (2.36) are combined
and form the full momentum conservation equation:

t j (0) +
+ c2 I +
t1 (0) + c2 I + 1 R(0)
= F . (2.39)
s
s
2
This equation is written in a divergence form, just like Eq. (1.1) on page 6. In
Section 2.3, a collision operator is developed in such a way that an exact match
between the terms in Eq. (2.39) and the macroscopic conservation equation (1.3) is
obtained. More precisely, the O( ) term (0) + c2 I is identied with uu + p I,
s
and the remaining O( 2 ) terms with .
t2 j (0) +

Multiscale Chapman-Enskog analysis

21

Note 2.4: Conservation laws A striking feature of LB models is the fact that
they implement conservation laws without error. Indeed, the conservation laws
prescribed by Eqs. (2.23) and (2.29) can be interpreted as follows: the conserved
variables and j yield the same value, independently on whether they are computed from the incoming distribution functions fi or the outgoing distribution
functions fi +i . This implies that the space average of those variables does not
change during the time evolution, unless a source term is introduced by the domain boundaries or by the external momentum source F . Of course, a numerical
implementation of the LB model exhibits errors in the conservation laws anyway,
due to the limited precision of oating point representations in a computer. By
exact conservation it is meant that none of the typical error terms that scale
3
3
at an order O(x ), O(t ), O( 3 ) or O(Ma 3 ) is present in the conservation laws.

2.2

Advection-diusion: Chapman-Enskog ansatz

The advection-diusion equation describes the evolution of a scalar quantity that

is diused and exposed to the advective inuence of an external term u:
t +

(u) = d

(2.40)

where d is the diusivity constant. The scalar represents for example the density
of a component in a multi-component uid, or the temperature distribution is a
uid. It is possible to solve this equation by means of a lattice Boltzmann method.
This is actually a relatively simple task, because the governing equation is linear,
and only one quantity, , is conserved. For this reason, this problem is used as an
introductory example to explain the ideas behind the Chapman-Enskog expansion.
The collision operator for the advection-diusion LB model can be written as a
BGK term, a relaxation towards a local equilibrium distribution:
i =

1
(eq)
fi fi
+ ti ci j.

(2.41)

Here, is the relaxation time, and j = j (1) an external source term for the
non-conserved rst order moment. It is introduced formally and will be used later
to correct errors in the numerical model. The equilibrium distribution reads
fieq = ti (1 +

1
ci u).
c2
s

(2.42)

It is clear that this collision operator veries the conservation requirements of

Eqs. (2.23) and (2.24).
The governing equation for this model is given by the conservation laws for the
zeroth-order moment, Eqs. (2.25) and (2.27), that are combined as follows:
t +

j (0) =

2
1

1
1
1
(0)
j (1) + t1 j (0) +
+ c2 1 .
1
s
2
2
2

(2.43)

22

Chapter 2

In order to nd the unknown term in this equation, a scale separation for the particle distribution function fi must be explicitly dened. In the Chapman-Enskog
approximation, the model is expanded around the local equilibrium term. That is,
the equilibrium term is identied with f (0) :
(0)

fi

= fieq .

(2.44)

It is immediately concluded, using the symmetry properties (a)-(d) in Eq. (2.1), that
j (0) = u and
(0)

(2.45)

= 0.

(2.46)

First-order terms with respect to are computed by rewriting Eq. (2.41) as follows
(the equation is divided by , by which this parameter is skipped):
(1)

fi

= i + ti ci j.

(2.47)

The collision term i is expanded in a Taylor series as in Eq. (2.19). Only rst-order
derivatives need to be retained at the O( ) scale:
(1)

fi

1
ci u + ti ci j
c2
s

= 2 ti Qi : 1 (u) ti ci 1 2 ti t1 ci (u) + ti ci j.
cs
cs
= ti (t1 +

ci ) 1 +

(2.48)

Projection on the rst-order base vectors C 1 yields

(1)

j (1) =

ci fi

1
t1
c2
s

ti ci ci u +

= t1 (u) +

ti ci ci

c2 1
s

c2
s

ci ci j
i

j.
(2.49)

Eqs. (2.49), (2.45) and (2.46) are inserted into Eq. (2.43), which nally yields the
O( 2 ) scale of the dynamics:
t +

(u) = d

+ d/c2 t
s

(u) c2
s

unwanted term

j ,

(2.50)

correction term

where the diusivity constant d has been dened as

d = c2
s

1
2

(2.51)

Equation (2.50) is almost equivalent to the desired PDE, Eq. (2.40), except for a
correction term on the right hand side. To the knowledge of the author, this term
has been overlooked in the literature so far. Reference [17] for example presents a
full Chapman-Enskog development of the temperature model (2.41) but somehow
fails to nd the error term. In Section 2.4, numerical evidence for the negative

Multiscale Chapman-Enskog analysis

23

inuence of error term on the accuracy of a simulation is presented. It is shown

in Section 3.1 that the free variable j can be adjusted to eliminate the unwanted
term.
This section concludes with a discussion of lattice structures that are appropriate
for the implementation of the LB advection-diusion model. It is pointed out that
among the symmetry relations in Eq. (2.1), only the properties (a)-(d) were used
in the previous theoretic development, and the property (e) in particular is not
required. Additionally to the lattice structures listed in Section 1.3.2, one may
use settings with weaker symmetry properties. A case in point are topologies that
include only velocity vectors parallel to the space axes. In 2D, this includes the D2Q4
model with two horizontal vectors ch = (1, 0) and vertical vectors cv = (0, 1).
The lattice weights are all equal and take the value ti = 1/4, and the lattice constant
c2 is dened as c2 = 1/2. Another possible choice is the D2Q5 model that includes
s
s
the zero-speed velocity c0 . In this model, the constant c2 is a free parameter. The
s
lattice weights take the value t0 = 1 2 c2 for the zero-speed velocity, and ti = c2 /2
s
s
for the four other ones. It is common to choose c2 = 1/3, with corresponding lattice
s
weights t0 = 1/3 and ti = 1/6 for i = 1 4. The same arguments lead to the 3D
lattice structures D3Q6 (ti = 1/6, c2 = 1/3) and D3Q7 (t0 = 1 3 c2 , ti = c2 /2 for
s
s
s
i > 0).

2.3

Fluid ows: Chapman-Enskog ansatz

The arguments leading to a LB model for uid ows are similar to those developed in the previous section for the advection-diusion model. The related algebra
is however more involved, because both the zeroth order moment (mass) and the
second-order moments (momentum) are conserved. The collision operator is again
represented by a BGK relaxation term with corrections:
(eq)

i = (fi fi

) + FT i + CT i .

(2.52)

The force term FT i and the correction term CT i are developed on the rst and
second order vectors C 1 and C 2 :
a
b
c F + ti 4 Qi : (F u + uF ) and
2 i
cs
2 cs
c
CT i = ti 4 Qi : .
2 cs
FT i = ti

(2.53)
(2.54)

The term multiplied by a factor a adds to the total momentum of the system and
represents the external force contribution in the momentum conservation law. The
term preceded by a factor b acts on the higher order moment and introduces a correction to numerical errors of the force term, due to the second order time derivative.
This approach has rst been described in Ref. [18]. Finally, a formal term is
introduced. It is used in Section 3.2 to correct a numerical deciency of the model.
This term is introduced, just like the correction term j in Eq. (2.41), at the level
of the rst non-conserved moment. In this way, it does not interfere with the exact conservation laws of the numerical model, which are explained in Note 2.4 on
page 21.

24

Chapter 2

The equilibrium distribution function fieq is constructed in such a way that the
collision term respects all the conservation laws listed in Section 2.1.2, and the
moments of the particle distribution functions are such that the momentum conservation equation (2.39) is equivalent to the Navier-Stokes equation. For this, the
moments of the distribution functions are expanded on the base vectors C 0 , C 1 and
C 2 , which leads to the expression dened in Eq. 1.18. Actually, the dynamics described by Eqs. (2.52), (2.54) and (1.18) can only be used to solve the Navier-Stokes
equation when the relative velocity u/cs is small. More precisely, the quantity u/cs is
said to be of the order of magnitude of the Mach number (Ma), and the model contains a numerical error that scales at the third order of the Mach number, O(Ma 3 ).
In the following development, all terms that scale at the third order of the Mach
number are therefore neglected. An additional assumption is used, stating that
the time-variations of the density are of the same order of magnitude as the Mach
number, or smaller:
t = O(Ma).
(2.55)
This assumption is compatible with results from the kinetic theory of gases close to
the isothermal limit [1]. Close to the limit of uid incompressibility, the density variations are actually even found to scale like the square Mach number (t = O(Ma 2 )
and = O(Ma 2 )). The assumption of Eq. (2.55) is less restrictive and leaves
space for the model to be pushed further away from its limit of incompressibility .
Note 2.5: Discrete Maxwellian The lattice Boltzmann equation is sometimes
derived by discretizing the continuous Boltzmann equation. In this derivation, the
equilibrium distribution corresponds to a Maxwellian distribution of the velocity,
and it is obtained by expanding this distribution for small Mach numbers, up to
a second order. This point of view is for example developed in Ref. [14].
It is obvious that the collision term in Eq. (2.52) conserves the zeroth order
moment (mass). Conservation of the rst order moment (momentum) reads
ci i = a +
i

F
2

by Eqs. (2.52) and (2.30).

For Eq. (2.29) to be respected, the value of the constant a must be set to

a=1 .
2

(2.56)

(2.57)

The governing equations for the LB uid model are Eqs. (2.35) and (2.39). To
obtain the missing terms in these equations, the dynamics is developed around the
equilibrium distribution function, as it was done in the previous section for the
advection-diusion model. This time, all symmetry properties listed in Eq. (2.1) are
used. The zeroth order terms read
j (0) = u,
(0) = uu and

2
T
1 R = cs
1 (u) + ( 1 (u)) +
1 (u) I .

(2.58)
(2.59)
(2.60)

Multiscale Chapman-Enskog analysis

25

At the rst order with respect to , it follows from Eq. (2.52) that

ti
1
1
1
= (t1 + 1 ci ) 1 + 2 ci u + 4 Qi : uu + (FT i + CT i )

cs
2cs

ti
1
1
= t1 + 2 t1 (ci u) + 4 t1 (Qi : uu)

cs
2cs

1
1
+ci 1 + 2 ci ci : 1 (u) + (ci 1 )(Qi : uu) + (FT i + CT i ).
cs

(2.61)

(1)

fi

The time-derivatives in this equation are replaced by space derivatives through the
following substitutions:

t1 = 1 u
by (2.25,2.58),

(2.62)

1
1
1
t1 (ci u) = 2 ci F (1) 2 ci 1 : uu ci 1
c2
cs
cs
s

by (2.31,2.59),
(2.63)

1
1
Qi : t1 ( uu) = 4 Qi : t1 (u)u + u t1 (u) (t1 )uu
4
2 cs
2 cs
O(Ma 3 )

1
Q : t1 (u)u
(because Q is symmetric)
c4 i
s

1
= 4 Qi : F (1) 1 (uu) c2 1 u
by (2.31)
s
cs

O(Ma 3 )

1
Q : F (1) c2 1 u.
s
4 i
cs

(2.64)

Inserting all these terms into Eq. (2.61) yields

(1)

fi

ti
1
Qi : 1 u ci 1 : uu + 2 (ci 1 )(Qi : uu)

2cs
1 ti
(b 1)ti
c ti
2 ci F +
Qi : (F u + uF ) + 4 Qi : .
4
2 cs
2 cs
2 cs

c2
s

(2.65)

When this equation is projected on the vectors of C 2 , using Eq. (2.12), only two
terms are left:
(1)

(1) =

Qi fi
i

2 c2
c
b1
s
S +
(F u + uF ) + ,

(2.66)

where the symmetric tensor S is dened by Eq. (1.6) on page 7. In order to nd the
last missing term, the time derivative of the (0) tensor, the same approximations
are used as those that lead to Eq. (2.64):

t1 (0) = t1 (uu) F (1) u + uF (1) c2 ( 1 )u + u( 1 ) .

s

(2.67)

26

Chapter 2

Equations ((2.58)-(2.60)), (2.62), (2.66), and (2.67) are now collected and inserted
into Eqs. (2.35) and (2.39). The resulting momentum equation contains among others an unwanted term that depends on the external force F , introduced by Eqs. (2.66)
and (2.67):
1
(1) + t1 (0) =
2

b1 1
+

(F u + uF )

2 c2
1
s
S + c2 (u)I.
s

(2.68)

This force term is eliminated by setting the constant b to

b=1

.
2

(2.69)

For the present discussion, the correction term is left out, and thus the constant
c is set to 0. This choice of the constants leads to the following equations:

(2.70)
t + (u) = 0 and

t (u) + c2 I + uu = F ,
(2.71)
s
where the deviatoric stress tensor reads
= 2 S,

(2.72)

and the dynamic shear viscosity is related to the relaxation parameter as follows:
= c2
s

1 1

(2.73)

This set of PDEs is identical to the equations of uid motion, introduced in Section 1.2.1. The bulk viscosity can however not be chosen as a free parameter in
the LB method, and it takes the xed value = 0. In order to vary this parameter,
one needs to exploit the degree of freedom introduced by the correction term .
Section 3.2 explains how this is done.
For convenience, the force term FT i in Eq. (2.53) is nally rewritten using the
choice of the constants a and b derived in the present section:
FT i = 1

ti
2

(ci u) ci u
+ 4 ci F .
c2
cs
s

(2.74)

This form of the force term has been suggested in Ref. [18] and shown to be superior
to other commonly used terms. Section 3.2 shows that one more correction needs
to be added to the force when uids with adaptable bulk viscosity are simulated.

2.4

Dimensionless formulation and accuracy

In the multi-scale analysis, the equivalence between the LB method and its associated macroscopic PDE is explicited via a nite development of the dierence
fi (r + ci , t + 1) fi (r, t) into a Taylor series. Because this series is truncated after
the second order term, it is common to say that lattice Boltzmann is a second order
accurate method in space and time. The present section discusses in some detail
what is meant by this statement, and why it should be appreciated with care.

Multiscale Chapman-Enskog analysis

2.4.1

27

Dimensionless formulation

During the implementation of a simulation, LB users are not always aware of the time
and space step t and x they are actually using. This question is therefore claried
to begin with. The approach adopted here is to relate lattice units to dimensionless
units (see Section 1.2.3) through those parameters. Another approach would be
to relate the lattice units directly with the units of a given physical system, but
this is less general and does not serve the purpose of the present discussion. It is
supposed that a reference length in the simulation has been resolved by means of
N lattice sites, and a reference laps of time is simulated by T iteration steps. The
corresponding scales are, by denition, unitary in the dimensionless system of units:
T = 1 and L = 1. This leads to the denition of the discrete time and space steps:
T
1
=
and
T
T
L
1
x =
=
.
N 1
N 1
t =

(2.75)
(2.76)

The term N 1 in the denition of x stems from a vertex-centered interpretation

of the location of lattice sites, as opposed to a cell-centered interpretation. The
lattice sites are in that case situated at the intersection between two cell edges, and
the number of cells spanned by N lattice sites amounts to N 1. It is common
to specify a reference speed U instead of a reference time T to set up a simulation.
The reference speed relates to the other parameters as follows:
U=

N 1
t
=
T
x

(2.77)

Furthermore, the uid viscosity = c2 ( 1/2), expressed in lattice units, is related

s
to the Reynolds number. To simplify the discussion, only the incompressible NavierStokes equation (1.10) is considered. By casting this equation into its dimensionless
form, given by Eq. (1.12), one obtains:
=

Note 2.6: Dimensionless or not?

1 t
.
2
Re x

(2.78)

The expression U = t /x in Eq. (2.77) is

somewhat puzzling, as the right hand side of the equation does not seem to have
the dimensions of a velocity. To be sure it is understood properly, let us go
through this once more. The velocity U can be viewed in two dierent systems
of units. One of them is the dimensionless system: it might be more appropriate presently to call it physical system, to emphasize that it is independent
of the lattice paramters. The other one is the system of lattice units. The lattice parameters t and x are used to restore the units of the physical velocity
Uphysical from the lattice velocity Ulattice : Uphysical = x /t Ulattice . Equation (2.77)
is deduced when the physical system of units is chosen such that Uphysical = 1.

28

2.4.2

Chapter 2

Accuracy of LB methods

The signication of the expression space accuracy of a numerical method is best

understood when a stationary (time independent) system is simulated, and the
obtained numerical solution solution u(r) is compared with an exact solution v(r)
of the problem. A common approach to dening the numerical accuracy is to say
that a method is of order k in space when there exists a constant for which the
following relation holds:
r (|u(r)

v(r)|)2

Nd

< k ,

(2.79)

where the sum is taken over all sites of a grid of size N d . The time accuracy can
be dened similarly for ows that are time periodic with period T . The accuracy is
then said to be of the order k in time when there exists a constant by which the
error between the numerical solution u(x) and the exact solution v(x) is constrained
as follows:
2
1
r |u(r, t) v(r, t)|
< k,
(2.80)
T t
Nd
provided the space discretization error is small enough to be negligible.
Additionally to the LB approach, there exist many numerical methods for the
resolution of PDEs in which the space derivatives are replaced by a nite Taylor
series up to the order kx , and the time derivatives by a series up to the order kt . It
can be formally shown that those methods have a space accuracy of the order kx
and a time accuracy of the order kt , in the sense of Eqs. (2.79) and (2.80). Although
no such formal proof exists for the LB methods, it is often argued on intuitive
arguments that a similar relation must hold for this method too. This is further
underlined by numerical evidence, such as the one presented in Section 3.2.

2.4.3

Accuracy of specic models

Some care must however be taken when the LB method is used to simulate incompressible ows. In those cases, the Mach number is chosen to be very small, and it
is concluded that the time derivative of the density scales then like the square Mach
number (t = O(Ma 2 )), by which the density is asymptotically described as a time
independent constant. In the LB method, the speed of sound is a lattice constant,
and the Mach number of the system is therefore equivalent to the reference speed
U expressed in lattice units: Ma = U/cs . Based on Eq. (2.77), it is concluded that
2
2
the error due to the compressibility eects scales like EMa = O(t /x ). There are
therefore in total three terms that contribute to the overall error E:
2
2
2
2
E = Ex + Et + EMa = O(x ) + O(t ) + O(t /x ).

(2.81)

The additional error term EMa introduces a coupling between the space and the
time discretization error. When the space step x is decreased, the time step must
2
be readapted for the overall error to decrease at a rate of O(x ). That is, the EMa
2
error term must scale as EMa = O(x ), and the time and space step are therefore
2
related through t x . This in turn means that the time discretization error really

Multiscale Chapman-Enskog analysis

29

behaves as in numerical methods with rst order accurate time stepping. Several
numerical case studies are conducted in Sections 5.3 and 6.2 that show that the
error term EMa has severe drawbacks on the eciency of the LB model.
Note 2.7: Keeping xed

When benchmarks are executed on several simulations with dierent grid size, a decision needs to be taken on how to adapt the
discrete time step to a given value of the space step. A choice often observed
in the literature consists in changing the scales in such a way that the Reynolds
number Re on one hand and the relaxation parameter on the other hand remain identical from one grid to another. This choice is a priori surprising,
because it is based on a microscopic property whose signicance is not simple
to interpret from a macroscopic point of view, as would be for example the time
step t , or the Mach number. With the help of the results presented in this Section, it is however possible to better understand this choice. Indeed, Eq. (2.78)
2
shows that the parameter is proportional to t /x . Keeping constant thus
2
amounts to varying t at a rate proportional to x . It has been argued in the previous paragraph that this is an appropriate way to scale the lattice parameters.
It guarantees in particular that the discretization error of a simulation decreases
at a second order rate when the grid resolution is increased.
The full power of LB methods is thus best exploited when compressible ows are
simulated at low Mach numbers. But even then, the LB method does not possess all
features one might expect from a good numerical method. First, the bulk viscosity
can not be adjusted in the basic BGK uid model. This drawback can luckily be
xed by taking into account a correction term discussed in Section 3.2. Another
problem however still subsists: the choice of the Mach number is intimately coupled
with the choice of the discrete time and space step, from which it is concluded
that those three parameters cannot be chosen independently. Indeed, in the LB
model that was introduced in this text, and which is systematically used in the
community, the speed of sound is a lattice constant. This means that it takes
a constant value (independent of x and t ) when it is expressed in lattice units.
The Mach number is thus identied with the reference speed expressed in lattice
units. From Eq. (2.77), the following relationship between the Mach number and
the discretization parameters is deduced:
Ma

t
.
x

(2.82)

For consistency, it should be remarked that it is principally possible to vary the

speed of sound in LB models. This is done by adapting the rest-particle weight t0 .
Choosing a speed of sound dierent than 1/3 has however been observed to induce
numerical instabilities and is therefore very uncommon.
It is interesting to note that the LB equation for advection-diusion problems
introduced by Eqs. (2.41) and (2.42) contains an error term that is very similar to
the error term EMa in the context of uid dynamics. To see this, it suces to turn
the PDE (2.50) to a dimensionless system, just as was done in Section 1.2.3 for the

30

Chapter 2

Navier-Stokes equations (the correction term j is excluded from the discussion):

(u ) = d

t 1
+ 2 2 d t
(u ),
x cs

(2.83)

2
where d = x /t d is the dimensionless diusion parameter. Given that c2 is an
s
immutable lattice constant, this equation contains an error term ELE due to lattice
2
2
eects that scales like t /x . As in Eq (2.81), the overall error contains therefore
three contributing terms:
2
2
2
2
E = Ex + Et + ELE = O(x ) + O(t ) + O(t /x ).

(2.84)

The eect of those error terms on actual numerical simulations is discussed in Section 3.1, and a strategy is devised to eliminate the error term ELE .
Note 2.8: Summary

All points discussed in this section can be summarized

by answering the following
Question: What is the time and space accuracy of a LB method?
Answer: The LB method is generally thought of as a second order accurate
numerical scheme, both in space and time, for the simulation of compressible, isothermal ows at small Mach numbers. This conjecture has however
never been proved formally, and it makes sense only when the widely unknown correction term in Section 3.2 is used. Furthermore, the discrete
time and space steps cannot be chosen freely in compressible ows, because
they are coupled with a physically signicant value, the Mach number. This
problem does not appear when incompressible ows are simulated, but in
that context, the LB model behaves like a numerical scheme with rst order
accurate time stepping.

Chapter 3
Corrections to the BGK dynamics

In Chapter 2, two lattice Boltzmann models are analyzed, one that describes advectiondiusion phenomena, and one that simulates uid motion. In both models, a correction term is formally introduced, which oers an additional degree of freedom
for handling deciencies of the numerical models. Those degrees of freedom are
exploited in this chapter. In the advection-diusion model, the correction term j
is adjusted so as to eliminate the lattice error ELE (cf. Section 2.4). In the uid
model, an appropriate choice of the correction term leads to a numerical scheme
in which the bulk viscosity is a free parameter (it was shown in Section 2.3 that
is immutable in the basic BGK model). A strong analogy in the way those two
problems are treated becomes apparent. In particular, it is shown that the deciency
of both numerical models can be eliminated by explicitly referring to the value of
the rest-particle distribution function f0 .

3.1
3.1.1

Advection-diusion revisited
Correction term for second order accurate time stepping

In Section 2.2 a BGK model for the modeling of advection-diusion processes is

proposed. Equation (2.50) on page 22 shows the macroscopic PDE related to this
model, obtained by a multiscale Chapman-Enskog analysis. This equation contains
an unwanted error term, as well as a correction term, depending on a not yet specied
correction vector j. In the present section, this vector is chosen in such a way as
to cancel out the unwanted term:

1
2

t1

(u) c2
s

j = 0,

(3.1)

which leads to the condition

j =

1
c2
s

1
2

t1 (u) =

1
c2
s

1
2

( t1 u u t1 ).
j1

(3.2)

j2

The correction term j contains two contributions. The rst, which will be named
j1 , as indicated in Eq. 3.2, depends on the time derivative of the velocity u, and the

32

Chapter 3

second, j2 , on the time derivative of the density. It is emphasized that a rst order
accurate estimate of those derivatives is sucient for the present purpose, because
the term t1 acts at a O( ) scale. In Eq. (2.19), this scale is related to the rst order
terms of the nite Taylor series.
The value of j1 can for example be found by a nite dierence approximation
to the time derivative. The dierence of the velocity between two successive time
steps is expanded in a Taylor series and scale separated as follows:
u(t + 1) u(t) = t1 u +

1 2
(t2 + t1 )u + O( 3 ),
2

(3.3)

which leads to an estimate of the O( ) scale of the time derivative:

t1 u(r, t) = u(r, t + 1) u(r, t) + O( 2 ).

(3.4)

Given that the correction term j scales at an order O( ), the rst-order nite
dierence suggested by Eq. (3.4) is sucient to approximate the time derivative
t1 u.
Whether this expression can easily be evaluated or not depends on the origin of
the external convective term u. A typical application area for the advection-diusion
equation is the analysis of thermal ows using the Boussinesq approximation. In that
case, the external vector eld u is obtained from a numerical solver for the NavierStokes equation. It is common for such solvers to keep in memory two successive
time steps of the velocity eld, and the nite dierence in Eq. (3.4) can be evaluated
eciently.
(1)
The value of j2 is directly found from the o-equilibrium part f0 of the restparticle distribution function. As a matter of fact, Eq. (2.48) leads to the following
relation:
1 (1)
f .
(3.5)
t1 =
t0 0
Equations (3.4) and (3.5) are nally combined with Eq. (3.2) as follows:
j =

1
c2
s

1
2

(u(t + 1) u(t))

1 (1)
f u .
t0 0

(3.6)

This correction term can of course only be implemented on lattice structures that
possess a rest particle distribution function, which rules out the D2Q4 lattice, but
works just ne with the D2Q5 lattice.
Note 3.1: Eciency

When the correction term in Eq. (3.6) is taken into account, the LB method for advection-diusion problems is unconditionally second
order accurate, both in space and time. Furthermore, it belongs to the family
of so-called explicit schemes, because the data at a time step t + 1 is obtained
from the data at a time t by an explicit calculation, and no expensive iterative
procedure needs to be undertaken. It is quite unusual for numerical methods with
second order accurate time stepping to be explicit, and the LB method is exceptional in this sense. Compared to other methods, the LB method can be said to

Corrections to the BGK dynamics

33

be very memory expensive, as it requires for example ve variables per lattice

site on a D2Q5 lattice, but very fast, as it steps through time at a second order
accuracy and with an explicit scheme.

3.1.2

Numerical verication

This section presents the results of a 2D numerical experiment on which the LB

model for advection-diusion (2.41,2.42) is tested, both in its original BGK formulation with j = 0, and with the modication term suggested in Eq. (3.6). The
theoretical discussion of the error terms in Eq. (2.84) is numerically veried, and
the benet of the correction term j on the accuracy of the model is underlined.
In the numerical example, a xed, divergence free and time independent velocity
eld u is imposed, and the advection-diusion equation for the scalar is solved
under the advective inuence of u. The velocity eld is dened on a square domain
whose side length serves as reference for the denition of a dimensionless length
scale. It corresponds to a Poiseuille prole whose maximum velocity umax denes
the reference scale for the velocity:
u (x ) = (u (x , y ), v (x , y )),
u (x , y ) = 4 y (1 y ) and
v (x , y ) = 0.

where

(3.7)

The domain is periodically extended in both the x- and the y-direction. The diused
scalar is initially dened through a uniform Gaussian distribution:

ini (r ) =

|r c|2
1

e 2 2 ,
2

(3.8)

where the distribution, in dimensionless variables, is centered at c = (0.5, 0.4) and

has a width of = 0.1. All simulations are run on a D2Q5 lattice.
The numerical results are compared to the results of a high accuracy reference
simulation. To be sure the time and space steps of the reference simulation are
suciently small, they are successively decreased until time and space convergence
is reached, that is, until a further renement does not inuence any more the conclusions drawn from the simulations. The reference simulation uses the BGK model
without correction term, so that no prejudice on the correctness of this approach is
made. The simulations are executed on a time interval 0 t < 0.2. The error of
the numerical result at a given time step, as compared to the reference simulation,
is dened through a l2 space norm, just as in Eq. (2.79). Similarly, the error of
a simulation on a full time interval is dened on ground of a l2 time norm, as in
Eq. (2.80). This procedure cannot be used to evaluate quantitatively the accuracy
of the time evolution, as it is dened in Section 2.4, because the ow in not periodic. It leads however to a qualitative estimation of the error, and an interesting
comparison between the BGK approach and the corrected model.
Figure 3.1 displays the numerical results obtained on this problem by means of
the original BGK and the corrected model. In Figure 3.1 (a), the space resolution

34

Chapter 3
(a)

(b)

10

10
BGK model
Corrected model
Slope 2

Relative error

Relative error

10

BGK model
Corrected model
Slope 1

10

10

10

10

10

10
1/t

10

10

10
Grid resolution N

Figure 3.1: Accuracy of the BGK and the corrected advection diusion models on
a chosen 2D problem. (a) Fixed diusivity d = 0.01 and grid resolution N = 100,
varying time step t. (b) Fixed diusivity d = 0.05 and time step t = 0.002,
varying grid resolution N .

is kept constant, and the time step is progressively rened. It is appropriate at

this point to recapitulate the conclusions drawn in Section 2.4 on the accuracy of
the numerical model. In Eq. (2.84) on page 30, three error terms were explicited,
2
2
the space error Ex = O(x ), the time error Et = O(t ) and an additional lattice
2
2
eect ELE = O(t /x ) that is present only in the uncorrected BGK model. The
space error Ex is not visible on Figure 3.1 (a), because the lattice resolution has
been chosen deliberately high enough. The remaining error terms Et and ELE scale
2
both at a same rate O(t ). This is conrmed on the graph, on which the error
visibly decreases at the same rate for both numerical methods. However, the oset
between the two curves shows that the term ELE has a more severe inuence on the
overall error than Et . On Figure 3.1(b), the time step t is kept constant, at a value
suciently small for the time error Et to be negligible, and the lattice resolution is
progressively increased. This time, the BGK simulation contains a decreasing error
2
2
Ex = O(x ), and an increasing error ELE = O(1/x ). This is clearly reected on
the graph: the error curve corresponding to BGK reaches a minimum after which it
increases again, whereas the error in the corrected model continues to decrease. As a
2
side remark, the rate of decrease is proportional to x and not x , as one might have
expected. It must be again underlined that those measurements of the accuracy
may not be interpreted quantitatively, because the ow is not periodic. It would
be interesting to devise a new numerical experiment that is time periodic, and on
which the rate of decrease of the numerical error can be veried quantitatively.

3.2
3.2.1

Bulk viscosity for uid ows

Numerical error in incompressible ows

Equation (2.81) on page 28 shows that there are three contributions to the overall
error of a simulation when an incompressible uid is being simulated: a space error

Corrections to the BGK dynamics

35

2
2
2
2
Ex = O(x ), a time error Et = O(t ) and a compressibility error EMa = O(t /x ).
In the present section, those theoretical considerations are rst underlined by a
numerical experiment. Then, a modied BGK model is introduced in which the
error term cancels.
The experiment, known under the name of Womersley ow, is now shortly introduced and presented in a system of dimensionless variables. Just like the Poiseuille
ow used in Section 3.1, the Womersley ow takes place in a channel whose width
determines the length scale for the dimensionless variables. A Poiseuille ow is a
ow parallel to the channel boundaries with parabolic velocity prole that is driven
either by an constant body force or by a constant pressure gradient DP . A pressure
driven Poiseuille velocity prole has the following form:

u (y ) =
x

Re
DP (y 2 y ).
2

(3.9)

The maximum value of the velocity in the middle of the channel is often taken as a
reference value to choose a time scale for the dimensionless variables. In that case,
the pressure gradient DP takes the value
8
p
= DP = .

x
Re

(3.10)

In a Womersley ow, the pressure gradient is varied periodically with a frequency

as follows:
p
= A cos( t ).
(3.11)
x
The resulting velocity prole is time dependent and substantially more complex than
the Poiseuille prole in Eq. (3.9). When the system is iterated long enough, it enters
a time periodic state for which an analytical solution is known. To choose a time
scale for this ow, the Poiseuille ow is taken as a reference, and the amplitude A in
Eq. 3.11 is identied with the pressure gradient DP in Eq. 3.10: A = 8/Re. Furthermore, a parameter is introduced that is known under the name of Womersley
number and dened as
Re
=
.
(3.12)
4
With those denitions, the analytical Womersley prole reads

cosh 2 y 1 ( + i)
DP it
2
,
u = Real
e
1
(3.13)

x
2
i
cosh
( + i)
2

where i is the imaginary number, and the term Real indicates that the real component of the right hand side expression must be extracted.
All simulations are executed with an unmodied BGK model, dened by Eqs. (2.52,
2.54, 1.18), with = 0 on a D2Q9 lattice. The parameters are Re = 1 and = 5.
The analytical solution of the Womersley ow is actually one-dimensional, given
that Eq. (3.13) depends only on the y-direction and not the x-direction. Therefore,
the x-direction, parallel to the channel walls, is resolved by only two lattice sites (one
lattice site would be insucient for the implementation of the pressure gradient).

36

Chapter 3

0.1
0.08

t*/T = 0.2

Velocity

0.06

t /T = 0.3

t*/T = 0.1

0.04
0.02
t*/T = 0

0
0.02
0

0.2

0.4
0.6
0.8
Position y* in the channel

Figure 3.2: Womersley prole, predicted by Eq. (3.13), at four chosen time steps
within a period T = 2/ , for Re = 1 and = 5.

Figure 3.2 shows the analytical solution of the Womersley problem, obtained from
Eq. (3.13) with the mentioned parameters. It is interesting to note that the prole
is seriously attened, compared to a Poiseuille prole whose maximum value would
be 1.
The simulations are executed until a time periodic state is reached, and their
relative error to the analytical solution in Eq. (3.13) is calculated by a time average,
as in Eq. 2.80. It is argued in Section 2.4 that the time step t in the BGK method
2
must be varied like x so as to cancel compressibility eects. The results of applying
2
this procedure, with t = 2/5 x , is shown in Fig. 3.3(a). As expected, the error
decreases exactly at a second order rate with the space step x . This could however
be achieved only via a dramatic adaptation of the time step. For that reason, as it is
argued in Section 2.4, the BGK method behaves like a numerical scheme with rst
order accurate time discretization, when used for incompressible ows. Figure 3.3(b)
shows what happens when the time step is xed to a value t = 2/5 1002 = 0.004.
When the grid resolution is rened, the error, dominated by the space error Ex ,
decreases with a second order rate. At some point however, the compressibility error
2
EMa , which scales like 1/x , becomes dominant. From then on, the error increases
at a second order rate.

3.2.2

Correction term for the bulk viscosity

The error ELE in the advection diusion model could be eliminated by an appropriate choice of the correction term CT i in Section 3.1. It would be tempting to take
a similar approach for the uid model and to look for a correction for which
the model is fully incompressible: EMa = 0. Things are however not that easy, and
to the knowledge of the author, there exists currently no lattice Boltzmann model
that simulates incompressible ows with an adequate accuracy. Though a seemingly
incompressible uid model has been presented in Refs. [19, 20], this model has been

Corrections to the BGK dynamics

37

(a)

(b)
2

10

10

Relative error

Relative Error

10

on

ve
r

ge

10

nc

10

N 2

rge

nc

10

nc

Co
nv
e

e
erg
Div

10

10
Grid resolution N

10

10

20

40
60
Grid resolution N

80

100

Figure 3.3: Accuracy of the BGK scheme compared to the analytical solution of an
2
incompressible Womersley ow. (a) t = 2/5 x (b) t = const. = 2/5 1002

shown to be erroneous (see e.g. Ref. [21]), as it eliminates compressibility eects

only in stationary ows. For that reason, the LB method is most eciently used to
simulate weakly compressible ows. In this regime at least, it should therefore be
able to provide the expected results. This is not the case with the plain BGK model
discussed in Section 2.3, because it does not oer the possibility to choose the bulk
viscosity in an arbitrary way. In the present section, this problem is circumvented
by an appropriate choice of the correction term CT i . This correction is proposed
in the literature in Ref. [22], and it is extended in the present context by the inclusion of a body force term. Indeed, the following developments show that a coupling
between the body force and the bulk viscosity correction occurs and requires special care. The section ends with the suggestion of an alternative correction term,
which achieves the same eect as the one proposed in [22] but is more ecient to
implement numerically.
In order to obtain a LB model with adjustable bulk viscosity, the correction
is implemented by means of the following ansatz:
1
= (Tr ((1) + k F u) I,
d

(3.14)

where Tr (T ) yields the trace of a tensor T and k is a constant to be determined.

The constant d corresponds to the number of space dimensions, and it enters the
calculations through the identity Tr (I) = d. Taking the trace of Eq. (3.14) yields
Tr () = Tr ((1) ) + k F u.

(3.15)

From Eqs. (2.66) and (2.69) on page 26 the trace of (1) is evaluated as follows:
Tr ((1) ) =

c
2 c2
Tr () s 1 u F u.,

(3.16)

Inserting Eq. (3.15) leads to

Tr ((1) ) =

1
1

c
2 c2
k 1 F u s 1 u .

(3.17)

38

Chapter 3

Equation (2.66) shows that the correction to the tensor (1) due to (1) is of
the magnitude c/ (1) . With the help of Eqs. (3.14) and (3.17) this expression
evaluates to
c
1
c
=

1 c/ d

1+

c
2 c2
k 1 F u s 1 u I.

(3.18)

The force term is eliminated by choosing the constant k to be

k=

1
c .
1+

(3.19)

Furthermore, by setting
c=

1
1

2 c2
s
d

(3.20)

Eq. (3.18) becomes

c
(3.21)
= 1 u I.

The momentum conservation law derived from the modied model is then the same
as in Eq. (2.71), but it includes the full deviatoric stress , as in Eq. (1.5).
To explicit the nal form of the correction term, Eq. (3.14) is inserted into
Eq. (2.54). Using the identity Qi : I = |ci |2 c2 d, one obtains
s
CT i =

c ti
|ci |2 c2 d
s
2 d c4
s

Tr ((1) ) + k F u ,

(3.22)

where the constants c and k are dened by Eqs. (3.20) and (3.19). The expression
Tr ((1) ) is evaluated using the denition of in Eq. (2.22), and the value of (0)
for the BGK model in Eq. (2.59). One nds
(1) =

ci ci fi uu c2 I, and consequently,
s

(3.23)

|ci |2 fi |u|2 c2 d.
s

(3.24)

i
(1)

Tr ( ) =
i

It is interesting to note that the calculations required in Eq. (3.24) can be circumvented by referring explicitly to the rest particle distribution function f0 , as was
done in Section 3.1 for the advection-diusion model. Indeed, the following relation
is directly deduced from Eq. (2.65):
(1)

f0 = f0 t0 =

t0
Tr ((1) ).
2 c2
s

(3.25)

The correction term CT i can therefore alternatively be written as follows:

CT i =

c ti
|ci |2 c2 d
s
2 d c4
s

2 c2 (1)
s
f +kF u .
t0 0

(3.26)

It is reminded that the adjustable bulk viscosity enters the equation via the
constant c, dened in Eq. (3.20), and that the constant d corresponds to the number
of space dimensions (2 or 3) of the system.

Corrections to the BGK dynamics

3.3
3.3.1

39

Alternative LB models
LB models with multiple relaxation times

The BGK model introduced in Section 1.3.2 can be extended to a more general
class of models, in which the collision operator , dened in Eq. (1.14) on page 10
takes the form of a diagonalizable matrix ij that acts as a linear operator on the
elements of the o-equilibrium particle distribution functions:
ij (fj fjeq ).

fi (x + ci , t + 1) fi (x, t) =

(3.27)

The restriction to the o-equilibrium part of the distribution functions is motivated theoretically by the multi-scale analysis in Section (2.1), in which the scale
of the collision operator is formally identied with the scale of the o-equilibrium
distribution functions. Using the denition of the vector space Rq introduced in Section 2.1.1, with the scalar product of Eq. (2.5) and the denition of the unweighted
particle distribution functions g in Eq. (2.13), Eq. (3.27) is interpreted as a MatrixVector product between the matrix and the vector g. Given that the matrix
is diagonalizable, there exists an invertible matrix M , whose lines represent the
eigenvectors of , and a diagonal matrix D, for which the following relation holds:
= M 1 D M .

(3.28)

Equation (3.27) is therefore rewritten in matrix notation as follows:

g out = g + (g g eq ) = g + M 1 D M (g g eq ),

(3.29)

out
where the short-hand notation gi g(x + ci , t + 1) for outgoing particle distribution functions has been used.

Note 3.2: Matrix notation In this section, a matrix notation is exceptionally

used on the space Rq of the particle distribution functions. The same notation
is used in the original papers that introduce multiple relaxation time models, and
it underlines the fact that linear algebra tools are used for developing new LB
models. One must be careful though not to mistake space vectors for the vectors
in Rq , and thus, keep in mind that M and D are matrices of size q q, whereas
the tensors Q and are represented by matrices of size d d.
The space representation of the distribution functions is now replaced by what
is called a moment space representation for reasons that become clear shortly, by
the change of variables G = M g. The moment space representation of Eq. (3.29)
reads
G out = D (G G eq ).
(3.30)
The last step consists in the formulation of a base formed by eigenvectors of the
desired collision operator. The choice presented here is taken from Ref. [21] and
encourages a physical interpretation of the algebraic operations. In this approach

40

Chapter 3

the base is rst populated by the 1 + d + d(d + 1)/2 vectors dened by C 0 , C 1 and C 2
in Eqs. (2.2), (2.3) and (2.4) respectively. It is now clear that the rst components
of the the vector G are constituted by the moments of the distribution functions.
With this new knowledge, Eq. (3.30) can be inspected more closely. It is obvious
that in the BGK model, both and D are dened as = D = I. Equation (3.30) shows that in this model, all moments are relaxed by a factor . An
exception is made for the conserved moments and u, because they are not represented in the vector G. Indeed, the projection of g g eq on the vectors C 0 and C 1
vanishes by virtue of Eqs.(2.24, 2.30, 2.44).
The idea behind multiple relaxation time (MRT) models is that each moment can
be relaxed at a dierent rate, by modifying the corresponding value on the diagonal
of D (see Ref. [23]). This approach can be used either to modify the physics of a
model by adapting the relaxation parameter of a hydrodynamically relevant moment,
or to increase the stability of the model by adapting the relaxation time of physically
irrelevant moments. The latter aim is for example achieved in Ref. [24] via a linear
stability analysis of the model.
The moments , u and related to the base vectors introduced so far are all
relevant to the hydrodynamics of the model. To form a complete base of Rq , they
must be completed by additional, physically irrelevant base vectors known as ghost
vectors. Their corresponding variables are named ghost variables. This point is
illustrated with the simple D2Q9 lattice, in which the space of particle distribution
functions, and thus also the space of the moments, is nine-dimensional. This moment
space can be specied via 6 physically relevant base vectors (1 for the density, 2 for
the velocity and 3 for the -tensor) and 3 ghost vectors. Guided by Ref. [22], a new
set of lattice constants hi is introduced as follows:
hi = (1, 2, 2, 2, 2, 4, 4, 4, 4) for i = 0 8.

(3.31)

It is used to dene a new lattice vector C 3 by

Ci3 = hi
4
and two lattice vectors C by

(3.32)

4
Ci, = hi ci .

(3.33)

A distribution function f can be projected on this base by purely algebraic means,

without inclusion of any physical ingredient. Some additional factors appear during
this projection, because the base vectors are not unitary:
fi = ti

1
1
ci u + 4 : Q + hi
2
cs
2 cs

3
1
N + ci J
4
8

(3.34)

where the projections on the ghost vectors have been named C 3 |f = N and
C 4 |f = J respectively. It is illuminating to split Eq. (3.34) into an equilibrium
and an o-equilibrium contribution:
fieq = ti
fineq = ti

1
1
ci u + 4 eq : Q
and
2
cs
2 cs
1
3
1 neq
: Q + hi
N + ci J
.
2 c4
4
8
s

(3.35)
(3.36)

Corrections to the BGK dynamics

41

With this notation, the BGK collision is written as follows:

i = ti

1 neq
: Q + hi
2 c4
s

1
3
N + ci J
4
8

(3.37)

A MRT model has exactly the same form, except for the fact that all non-conserved
moments possess their own relaxation parameter, numbered through from 3 to 8 :
1
3 neq Qxx + 4 neq Qxy + 5 neq Qyy
xx
xy
yy
c4
s
1
3
+ hi 6 N + (7 cix Jx + 8 ciy Jy ) .
4
8

i = ti

(3.38)

Equation (3.38) is particularly illuminating, because it displays both the space and
moment representations of the dynamics at the same time. In Section 4, a new
LB model is introduced whose collision operator is linear and diagonalizable. The
MRT interpretation of this model, presented in Section 4.3, will show that the
hydrodynamic parameters 3 to 5 are the same as in the BGK model, whereas 6
to 8 are modied.
Note 3.3: Terminology

Some confusion in the terminology occurs sometimes

when the term MRT model is used for all LB models with a linear, diagonalizable collision operator. This puts the emphasis on the wrong place. Such a terminology would actually imply that practically all commonly used LB models are
MRT models, including the BGK model, and several other models that were developed even before the term MRT existed. In particular, specic LB models based
on a linearized collision operator have been described in Refs. [10, 25, 26, 27].
In the literature, the name MRT is rather used for models that where explicitly
constructed by switching to a moment representation and ne-tuning the relaxation parameters.

3.3.2

Entropic LB models

The family of entropic LB models has emerged over the past decade and is successively gaining importance. Two approaches were developed independently by
Karlin, Ansumali e.a., as documented in Refs. [28, 29], and by Boghosian e.a., documented in Ref. [30]. Those models exploit a discrete-space analogue of Boltzmanns
H-theorem to achieve unconditional numerical stability. Numerical stability means
that the simulated values always remain within a xed range and never exceed the
storage capacity of oating point representations in a computer. It does however
not imply that the simulated values systematically stay close to the exact solution of
the problem, nor that they even yield a physically meaningful result. Unconditional
numerical stability can be very useful when simulations are run close to a regime in
which they are underresolved. It might happen that the simulation is locally underresolved during a short time interval, which should only aect the overall quality
of the simulation to a small extent. In such a situation, numerical instabilities can
however occur in relation with the BGK model, but not with entropic models.

42

Chapter 3

Boltzmanns H-theorem states that each uid possesses a property H that never
increases during the time-evolution of a closed system (a system without energy
input). A local version of this theorem additionally states that a locally dened
property H = H(f (x, t)) never increases during particle collisions. Because of numerical errors, this theorem is only approximately true in BGK simulations. In
entropic models however, the collision operator is modied so as to obtain an exact
discrete version of the H theorem. This can be compared to the local conservation laws for mass and momentum, which are exactly respected in the discrete LB
equation.
The unconditional numerical stability is achieved by means of the three following
steps:
1. A function H = H(fi (x)) is constructed, which is dened over the space of
strictly positive particle distribution functions fi , and which is convex. The
minima of this function is given by the local equilibrium fieq .
2. The H value for the incoming distribution functions H in = H(fiin ) is evaluated,
and the set of distribution functions fi for which H(fi ) H in is evaluated.
3. The collision is modied so as to ensure that the outgoing distribution functions fiout are contained within the set of distribution functions evaluated previously. This in its turn ensures that all distribution functions remain strictly
positive. It has actually been observed empirically that numerical instabilities
only occur in the LB method when negative distribution functions appear.
This eect is technically prevented in entropic methods by the exact H theorem.

Chapter 4
Regularized lattice Boltzmann

4.1

Introduction

In Section 1.3, a lattice Boltzmann model is introduced that is expressed in a system

of lattice units. In this system, the spacing between two adjacent lattice cells, as
well as the time span between two successive iteration steps, is unitary. To motivate
why this system of units is a natural choice in LB simulation, it is useful to adopt the
interpretation of the equilibrium distribution function in Eq. (1.18) on page 10 as a
small Mach number expansion of a Maxwellian velocity distribution (see Note 2.5).
The Mach number is a dimensionless number computed from the ratio between a
characteristic velocity of the system and the speed of sound. As it was shown in
Section 2.3, the speed of sound depends on the choice of the discrete space and time
steps, and it is constant in a system of lattice units. On a D2Q9 lattice for example,
and in a system of lattice units, the speed of sound is dened by c2 = 1/3. This
s
shows that apart from this factor 1/3, the characteristic ow velocity u is equivalent
to the Mach number. The equilibrium distribution function is thus most naturally
written in lattice units, in which it does not explicitly depend on the parameters x
and t . This approach contrasts with one of other common CFD tools, where the
system of dimensionless units introduced in Section 1.2.3 is the most natural choice
for the formulation of a dynamics independent of discretization parameters.
It is of course possible to switch from one system of units to another. To interpret the results of a LB simulation in terms of dimensionless macroscopic variables,
the moments of the distribution functions are computed and then rescaled in an
appropriate way. The dimensionless velocity is for example obtained through the
expression
u =

x
x 1
u=
t
t

fi ci .

(4.1)

It is much more contrived however to initialize the particle distribution functions of a

LB simulation from a given prescription of the macroscopic ow elds. The results
from the Chapman-Enskog multiscale analysis are used to explicit this relation.
From Eqs. (1.18) and (2.65), sparing out the force term and the correction term to

44

Chapter 4

keep the discussion simple, one obtains:

(1)

fi = fieq + fi

+ O( 2 )

1
1
c u + 4 Qi : uu
2 i
cs
2cs

1
Qi : u ci : uu + 2 (ci )(Qi : uu) + O( 2 ).
2cs

= ti 1 +

ti
c2
s

(4.2)

This equation shows that in order to relate a macroscopic ow eld to a LB eld,

various types of velocity gradients need to be computed. Furthermore, the relatively
complicated expression of Eq. (4.2) needs to be evaluated. This is technically feasible, but tends to be very expensive if it has to be done at every iteration step.
This conversion is to be executed not only when values from an external ow eld
are imported into a LB simulations, but also when two LB simulations, executed
on dierent grids, are related to each other. This becomes clear when Eq. (4.2) is
rewritten in terms of dimensionless variables:
t 1
2 1
ci u + t 4 Qi : u u
2
x c2
x 2cs
s

2
1
t Qi : u t ci : u u + 2 (ci )(Qi : u u )
x
2cs

fi = ti 1 +

ti
c2
s

.
(4.3)

It is concluded from this equation that dierent contributions to the particle distribution functions scale at a dierent rate under a change of x and t . To compare
the results from two LB simulations that are run on dierent grids, it is therefore
necessary to split the fi into their components, and to rescale each component individually. It is possible, using the values contained locally on a lattice cell, to compute
and u, and thus to compute the various components of f eq , and to separate the
equilibrium contributions from the o-equilibrium part. It is however not possible to
split the o-equilibrium part into its components, because not all velocity gradients
are known locally on a lattice cell: only the symmetric part of the Jacobian matrix
can be computed by means of Eq. (2.66). The skew symmetric part needs to be
evaluated by some other, non-local method.
It is possible to rescale the particle distribution functions by purely local means
if the O(u2 ) terms in the o-equilibrium distribution functions are neglected. In
that case, all components scale in a known way: the density is invariant, the
(1)
velocity u scales like t /x , and the o-equilibrium distributions fi scale like t .
This approach is chosen in Ref. [31] to formulate a grid-renement algorithm for
the LB method. It must be pointed out that this is one of the rare papers on a
related topic that explicitly tackles the problems of scale invariance between two
grids during a grid renement procedure.
A dierent approach is presented here that is not based on an approximation
of the distribution function. The idea behind this method is to replace the BGK
model by a new model, the regularized LB (RLB) model, which has less degrees
of freedom. This new model depends only on the value of , u and (1) . Each of
them can be rescaled individually, as the tensor (1) is related to the tensor S of
Eq. (1.6) through Eq. (2.66). The details of the regularized method are explained

Regularized lattice Boltzmann

45

in Section 4.2. It is found that additionally to the advantages mentioned previously,

the regularized LB model is exceptionally accurate and numerically stable. This
point is veried in Section 4.4 on chosen numerical examples.
In the subsequent chapters, various typical problems of computational uid dynamics are listed in which the variables of a lattice Boltzmann simulation need to be
rescaled. The regularized LB model is proposed as a solution in each of those cases,
and the change of scales is handled in a straightforward way via a dimensionless formulation of the variables , u and (1) . The described procedures are actually also
applicable to the BGK model, in which the tensor (1) can be computed locally.
The obtained results are however slightly less accurate, because by doing so, the
O(u2 ) terms in the o-equilibrium parts of the distribution functions are not taken
into account.

4.2

Regularization procedure

4.2.1

Regularized particle distribution functions

The particle distribution functions fi are expanded by means of a multi-scale analysis

(1)
in Section 2.3. In particular, the value of the rst order contribution fi is explicited
in Eq. (2.65) in terms of the macroscopic variables and u. This term is however
of limited interest, because it does not explicitly appear in the macroscopic PDEs
that describe the hydrodynamic behavior of the model, Eqs. (2.35,2.71). Instead,
(1)
Eq. (2.71) refers to the tensor (1) , which is the projection of fi on the second
2
order base vectors C shown in Eq. (2.66). During this projection, many terms
(1)
2
cancel out, and only those components of fi that are expanded on C are kept.
(1)
The idea of the regularized LB method is to simplify the value of fi accordingly,
and to keep only those terms that are relevant to the computation of (1) .
(1)
It is useful at this point to take up again the expression for fi from Eq. (2.65)
and to split it into contributions that are relevant (Ri ) resp. irrelevant (Ii ) to the
computation of (1) :
(1)

fi

= Ri + Ii ,

where

1
1
c
and
Ri = ti Qi : 2 1 u 4 F u + uF + 2
cs
4 cs
2 cs

1
1 ti
ti
ci 1 : uu 2 ci 1 (Qi : uu) 2 ci F .
Ii = 2
cs
2 cs
2 cs

(4.4)
(4.5)
(4.6)

The term Ri is of the form Ri = ti Qi : T , where the tensor T stands for the whole
content of the parentheses on the right hand side of Eq. (4.5). The tensor (1) was
2
calculated in Eq. (2.66) by projecting Ri on the vectors C as follows:
(1)

ti Qi Qi T = c4 (T + T ) by Eq. (2.12).
s

(4.7)

The key to the regularization procedure is the observation that the evaluation of
(1) in Eq. (4.7) can be reversed. Once the tensor (1) is known, the tensor T can

46

Chapter 4

be evaluated in its turn. This is done by contracting (1) with the tensor Q:
Qi : (1) = c4 Qi : (T + T T ) = 2 c4 Qi : T
s
s

by symmetry of the tensor Q. (4.8)

(1)
(1)
From this, a regularized substitute fi for fi is derived:

ti
(1)
fi = Ri = 4 Qi : (1) .
2 cs

4.2.2

(4.9)

Dimensionless formulation

It is interesting to remind that up to the accuracy of the hydrodynamic terms, the

tensor (1) can be related to the strain rate tensor S. Indeed, using Eq. (2.66) on
page 25 yields, when the force and the correction term are not considered:
ti
(1)
fi 2 Qi : S.
cs

(4.10)

The tensor S, dened by Eq. (1.6), is a macroscopic quantity that depends only on
the gradients of the velocity eld. To conclude this discussion, it is pointed out that
the set of lattice kinetic variables dened by the regularized particle distribution

(1)
functions fi = fieq + fi is equivalent to a set of macroscopic variables dened by
the scalar , the vector u and the tensor (1) . Those variables, and thus the state
of the simulation, can be represented in a dimensionless system of units:
= ,
t
u =
u and
x
(1) = t (1) .

4.2.3

(4.11)

Algorithm

To summarize, the regularized LB method replaces the usual BGK method in

Eqs. (1.14,2.52) by the following procedure:
(eq)

1. Compute the equilibrium distribution fi

tions by means of Eq. (1.18).

from the particle distribution func(neq)

2. Determine the o-equilibrium distribution functions fi

(neq)
related tensor (neq) = i Qi fi
.

(eq)

= fi fi

and the

(1)
3. Derive the rst order regularized distribution functions fi from (neq) through
the approximation (1) (neq) in Eq. (4.9).

4. Replace the distribution functions fi by their regularized counterpart fi =

(eq)
(1)

fi + fi .
5. Execute the BGK collision and the streaming step.

Regularized lattice Boltzmann

47

In order to clarify this algorithm, the original BGK LB model (Eqs. (1.14,2.52))
is rst reformulated as follows:
(eq)

fi (r + ci , t + 1) = fi (r, t) (fi (r, t) fi

(eq)

= fi

(neq)

(, u) + (1 ) fi

(, u)) + FT i + CT i

(4.12)

(r, t) + FT i + CT i .

By analogy to this equation, the regularized model can be summarized as follows:

(eq)

(1)
(, u) + (1 ) fi (r, t) + FT i + CT i
(1 ) ti
(eq)
= fi (, u) +
Qi : (1) (r, t) + FT i + CT i .
2 c4
s

fi (r + ci , t + 1) = fi

4.2.4

(4.13)

Related models

It is emphasized that the regularization procedure described in this section principally applies to all LB models of the BGK family. The idea is throughout the same:
(1)
2
the rst-order term fi is projected on the vectors C corresponding to the rst
(1)
non-conserved moment, and the resulting regularized term fi is directly computed
from the value of this non-conserved moment. For the purpose of illustration, the
advection-diusion model described in Section 2.2 is now regularized. When the error and correction terms of this model cancel, for example by means of the procedure
described in Section 3.1, the rst-order contribution to the model reads
(1)

fi

t Qi : 1 (u) ti ci 1 .
2 i
cs

(4.14)

The rst non-conserved moment of this model is the momentum j, and the regularized version of Eq. (4.14) takes therefore the following form:

(1)
fi = ti ci 1 .
Indeed, the momentum computed from those two expressions is the same:

j = c2 1 ,
s

(4.15)

(4.16)

(1)
from which the regularized term fi is computed as follows:
ti
(1)
fi = 2 ci j.
cs

(4.17)

The LB model introduced here under the name of RLB has previously been
described in the literature, where it was derived from dierent principles and applied
in dierent settings. In Ref. [32], this model is used for the simulation of particle
simulations, whereas in Ref. [33], it is introduced as a model for the simulation of
high Knudsen number ows. The derivation of RLB from the Chapman-Enskog
expansion of a BGK model however is novel. Furthermore, the particular properties
of RLB, such as the enhanced stability and accuracy, as well as the adequacy of the
model to situations with varying time and space scales, have not been pointed out
before.

48

Chapter 4

Note 4.1: Memory savings

The state of a RLB simulation is fully described

by the three quantities , u and (1) . It is therefore possible to write a program
that holds only those variables in memory instead of the particle distribution
functions. This leads to memory savings because it requires a number of 1 +
d + d(d + 1)/2 variables to be stored, less than the usual q particle distribution
functions. On the other hand, the resulting code is quite dierent from a typical
LB code and may be more dicult to work with. Furthermore, it depends on the
particular software and hardware platform whether this alternative formulation
of RLB leads to performance benets or not.

4.3

RLB and MRT

The collision operator of the RLB model is linear and acts only on the o-equilibrium
parts of the particle distribution functions. This fact appears clearly when Eq. (4.13)
is formulated as follows:
(eq)

fi (r + ci , t + 1) = fi

(, u) + (1 )
j

ti
Q : cj cj fj fjeq ,
4
2 cs

(4.18)

or, to emphasize the analogy with a general collision operator as in Eq. (3.27) on
page 39:
1
fiout =
Qi : cj cj ij fj fjeq .
(4.19)
4
2 cs
j
The RLB model can thus be reinterpreted as a LB model with multiple relaxation
times, by using the theory introduced in Section 3.3.1. As an example, the D2Q9
lattice is now considered. Using the same notation as in Eq. (3.37), and after some
algebra, the collision operator of Eq. (4.19) is written as follows:
i = ti

1 neq
: Q + hi
2 c4
s

1
3
N + ci J
4
8

(4.20)

This means that all modes of neq are relaxed with a parameter , as in the BGK
model, and all other modes are relaxed with a parameter 1.

4.4

Numerical verication

We now turn to numerical verications of the regularized model on two 2D ows

using a D2Q9 lattice. The rst test case, known under the name of Kovasznay ow,
approximates a stationary 2D ow behind a regular grid. An analytical solution for
this ow, proposed in [34], is written as follows in dimensionless units:
u = 1 exp( x ) cos(2y ) and
x

u =
exp( x ) sin(2y ),
y
2

(4.21)

Regularized lattice Boltzmann

49

Relative l2 norm error

10

Slope 2

10

10

10

Slope 3

BGK, Inamuro: slope 2.42

BGK, Skordos: slope 2.11
RLB, Inamuro: slope 2.94
RLB, Skordos: slope 2.87
20
30
Grid resolution N

40

50

60

Figure 4.1: Relative error of the numerical result on a Kovasznay ow in the wake
region. Both traditional BGK and the regularized model are tested with two commonly used boundary conditions.

with
= Re/2

4 2 + Re2 /4.

(4.22)

The Reynolds number Re = u L/ is dened as a function of the asymptotic

velocity u of the velocity in the wake, far from the grid, the spacing L between
two grid nodes and the dynamic viscosity of the uid. The simulations are executed
in the wake of the grid, in the intervals x [L/2, 2L] and y [L/2, 3L/2], with
Re = 10. The space step x is varied linearly, while the ratio between the discrete
time and space step is xed to a value t /x = 0.01. As it follows from the discussion
in Section 2.4, keeping the ratio t /x constant amounts to xing the Mach number
M a = u/cs . This value is chosen small enough to mimic an incompressible ow.
Given that the ow is periodic in the y-direction, the upper and the lower boundary
of the simulation can be chosen to be periodic, whereas the Kovasznay solution
in Eq. (4.21) is imposed through Dirichlet boundary conditions on the left and
right boundary. After the simulation has reached a time independent state, the
numerical result is compared with the solution in Eq. (4.21) through an l2 norm on
each grid point, and then averaged over space. The result is shown in Fig. 4.1, on two
commonly used implementations of the boundary conditions. One of them has been
described by T. Inamuro e.a. in Ref. [35] and the other one by P. Skordos in Ref. [36].
The accuracy of the simulation with respect to the grid resolution is of order 2 to
2.5 when the BGK model is used, whereas the regularized model is almost thirdorder accurate. On the BGK simulations with the Inamuro boundary condition,
data points for small grids are missing because numerical instabilities make them
impossible, whereas the regularized model has no such stability deciencies.
The second test case implements a ow in a 2D square cavity whose top-wall
moves with a uniform velocity. Reference values for the denition of the Reynolds
number are dened by the side length L of the box and the top-wall velocity u0 .
Both standard BGK and the regularized model are rst compared with the reference

50

Chapter 4

solution of Ghia e.a. [37], on a lattice size of N N with N = 129, at Re = 100 and
a Mach number, in lattice units, dened as u0 = t /x = 0.02. This time, another
commonly used boundary condition is applied, which is described by Q. Zou and
X. He in Ref. [38]. The reference solution [37] proposes a set of accurate numerical
values for some x- and some y-components of the velocity on chosen space points.
A l1 norm error with respect to these reference points is averaged over all available
points and normalized with respect to u0 . For the BGK model, this yields an error of
= 3.71103 , and for the regularized method, of = 2.40103 . Thus, both methods
solve the problem with satisfying accuracy. This fact is underlined graphically on
Fig. 4.3, on which the numerical results of the RLB simulation and the solution by
Ghia e.a. are compared. The regularized model is however found to be substantially
more stable. To make this statement more quantitative, a series of simulations is
run, on which the velocity is again xed to a value of u0 = 0.02. For several chosen
grid sizes N , the maximal Reynolds number Remax at which the simulation remains
stable (i.e. delivers nite numerical values) is determined. Figure 4.2 shows that,
although both methods exhibit a linear relationship between Remax and N , the
observed increase rate is 7.7 times higher for the regularized method than for BGK.

Regularized lattice Boltzmann

51

Maximal Reynolds number Remax

1000
Regularized: Remax = 16.3+7.36*N
800

600

400

200
100
0

BGK: Remax = 0.391+0.955*N

50

60

70

80 90 100 110 120 130

Grid resolution N

Figure 4.2: Simulation of 2D cavity ow for xed Mach number. , : maximal

stable Reynolds number, numerically determined; solid line: least-square linear t
of the data points (parameters of the t are indicated on the graph).

(a)

(b)
0.2

u* on horizontal line
y

u* on vertical line
x

0.5

0.1
0
0.1
0.2

0.5
0

32
64
96
Position [lattice units]

128

0.3
0

32
64
96
Position [lattice units]

128

Figure 4.3: Comparison of computed RLB results with reference values from the
literature for the 2D lid-driven cavity ow at Re = 100. The solid line represents
the reference values, which have been interpolated to guide the eye, and the dots
stem from numerical results. (a) x-component of the velocity along a vertical line
passing through the center of the cavity. (b) y-component of the velocity along a
horizontal line passing through the center of the cavity.

Chapter 5
Boundary and initial conditions

5.1

Introduction

The Navier-Stokes equation possesses a unique solution when adequate initial and
boundary values for a given problem are specied. The problems of uid dynamics
are therefore initial and boundary value problems. It is common to specify the initial
value for a problem by imposing the velocity eld at the initial time t = 0. Several
approaches exist for specifying the boundaries of a ow. It is possible to implement
a so-called Dirichlet boundary condition by specifying a velocity prole along the
domain boundaries of the problem. In other commonly adopted approaches, the
value of the velocity gradient in the direction of the boundary normal is specied,
or the value of the pressure on the boundary.
Whatever approach is adopted, one major diculty needs to be addressed in
LB simulations: the conversion from macroscopic variables, which are taken from
the initial or boundary condition, to the particle distribution functions of a LB
simulation. When the regularized LB model is used, it is natural to implement this
conversion by means of the approach described in Chapter 4. In this approach,
the tensor (1) is computed additionally to the pressure p and the velocity u. By
virtue of Eq. (4.9) on page 46, this set of variables fully species the state of the LB
simulation.
During the implementation of an initial condition, the tensor (1) can be related
to the strain rate tensor S, as dened in Eq. (1.6) on page 7. The velocity gradients
contained in this tensor can in their turn be computed either analytically, when
the initial velocity prole is dened through an analytical expression, or otherwise
numerically, by applying for example a nite dierence scheme. It is actually shown
in Section 5.3 that it is more dicult to properly initialize the pressure eld in the
initial condition. In this section, a time-dependent benchmark problem is introduced
whose results depend heavily, among others, on the initial condition of the problem.
The initialization procedure described here, inspired from a macroscopic vision of
the problems of CFD, is tested numerically. It is compared with another approach
taken from the literature that rather adopts a microscopic view and is implemented
fully in terms of LB variables.
When the boundary condition is implemented, the velocity gradients cannot be
evaluated analytically, even when the velocity prole along the boundary is known
in terms of an analytical formula. Indeed, the gradients normal to the boundary

54

Chapter 5

c1
c2

c3

c4

c7

c6

c8

Boundary

c5

e1

e0

Figure 5.1: Orientation of a 2D boundary for the example developed in the text.
Dashed arrows label directions that originate from locations outside the lattice.

depend on the results of the simulation and are a priori unknown. It is possible
to evaluate those gradients by means of a nite dierence scheme, and close the
system of equations on the boundary in this way. Although this approach is feasible,
it violates partly the spirit of the LB methods, in which all operations, except for
the streaming step, tend to be purely local. In Section 5.2, another approach is
therefore shown, in which the elements of S are computed from the known particle
distribution functions on the wall.

5.2
5.2.1

Boundary condition
Overview

This section focuses on the implementation of Dirichlet boundary conditions for the
velocity eld. Following the approach discussed in the introduction, implementing
the boundary condition amounts to computing the values for the missing variables
and (1) . In this discussion, only straight boundaries are considered that are parallel
to the lattice edges and pass through lattice nodes. Other types of boundaries, such
as corner nodes, can be implemented via the same approach. The discussion of this
topic is however omitted here because it is quite technical. The interested reader
is invited to directly refer to the publicly available and well documented code of
the OpenLB project in Ref. [13]. The orientation of the boundary is taken to be
orthogonal to the unitary vector e1 . This means that for a lattice location r that
lies on the boundary, the location r + x e1 lies outside the computational domain,
and the location r x e1 lies inside, as shown on Fig. 5.1. As a consequence, all
neighbors of a boundary location r, placed at a position r ck , are unreachable when
ck1 = 1. However, particle streams parallel to the boundary, along directions cl
for which cl1 = 0 participate in the LB dynamics of the simulation. This approach
is based on a point of view advocated in Ref. [39], according to which boundary
nodes are located inside the uid, but innitesimally close to the wall. The problem
with the computation of and (1) resides in the fact that among all distribution
functions fi on the boundary those with an index i {k | ck1 > 0} are unknown and
must be evaluated on ground of some additional closure relations.

Boundary and initial conditions

5.2.2

55

Implementation of the pressure

The rst step consists in computing the variable on the wall. In an incompressible
uid, Eq. (1.8) on page 7 shows that this is equivalent to computing the pressure.
A fairly standard procedure for doing so is described for example in Refs. [38, 35].
It is based on the denition of three separate contributions to the density:
1 =

fk ,

(5.1a)

{k | ck1 =1}

0 =

fk ,

and

(5.1b)

{k | ck1 =0}

1 =

fk ,

(5.1c)

{k | ck1 =1}

of which the component 1 is unknown on the boundary. From Eqs. (1.15) and (1.16)
on page 10, one nds that
= 1 + 0 + 1

and u1 = 1 1 ,

(5.2)

and can be formulated in terms of known quantities:

=

1
(21 + 0 ).
1 u1

(5.3)

Note 5.1: Pressure boundaries

Although only Dirichlet velocity boundaries

are discussed in this text, other boundary types are dened in a very similar
way. Equation (5.3) can for example be inverted to compute the normal velocity component u1 from the density . This leads to the denition of a so-called
pressure boundary, on which the density and the tangential velocity u0 are
specied, and the normal velocity u1 is free.

5.2.3

Regularization of on-wall distribution functions

Now that the density and the velocity u are known, a last closure relation for the
evaluation of the stress tensor (1) is required, before the regularized distribution
functions can be calculated via Eq. (4.9) on page 46. An alternative approach
consists in computing the macroscopic strain rate tensor S instead of (1) and by
evaluating the regularized particle distribution functions via Eq. (4.10) on page 46.
In that case, the velocity gradients of the strain rate tensor are found with an
asymmetric second-order accurate nite dierence scheme on the boundary. This
computation involves values residing on nearest and a next-to-nearest neighbors of
the boundary node. The resulting BC is similar to the one presented in Ref. [36] in
which the values of the distribution functions on the boundary are tted according to
their zeroth- and rst-order terms of the Chapman-Enskog expansion. It is however
not typical for simple LB methods to make use of non-local nite dierence schemes,
as one prefers to implement a collision step depending only on values dened locally
on a lattice site, both for conceptual and technical simplicity. A purely local means

56

Chapter 5

to computing (1) is therefore presented, which is based on knowledge of the particle

distribution functions fi that are known on the boundary. The dynamics of the wall
is implemented in two steps:
Compuation of (1) . To all missing distribution functions fi on the wall, attribute
a value fi fieq (, u) + fjneq , where j is the index of the velocity opposite to ci :
ci = cj . This procedure is commonly applied and has been called bounceback of o-equilibrium parts in Ref. [38]. It is motivated by symmetries in
the regularization formula Eq. (4.9), from which it is easily deduced that the

o-equilibrium part for two distribution functions fi attributed to opposite

directions is equal. In conclusion, the obtained distribution functions possess
values that are compatible with the order O( ) development of the BGK dynamics. They are however not yet fully satisfying, as they lead to slightly
erroneous values for the density and the velocity u. For this reason, they are
only used to compute the value of the stress tensor (1) = i fi c2 I uu
s
on the wall.
Regularization. Replace all distribution functions on the wall by their value obtained from the variables , u and (1) via Eq. (4.9).
Note 5.2: Conservation laws

In the two-step implementation of a boundary

condition presented in this paragraph, it can seem appropriate to simply skip step
(2). Indeed, the missing particle distribution functions on the wall are initialized
in step (1) to a reasonable value that respects the conclusions of the ChapmanEnskog analysis. The problem is that the exact conservation laws for the mass
and the momentum (see Note 2.4 on page 21) would not be respected by this
procedure. The distribution functions computed in step (1) do for example not
respect the relation i fineq = 0, which is required for exact mass conservation.
It must be noted that exact conservation laws are important both for theoretical
and practical reasons, as they ensure for example a good numerical stability of
LB models. The regularization step (2) restores symmetries of the distribution
functions and leads to exact conservation laws. The boundary condition described
in Ref. [38] chooses to apply the bounce-back of o-equilibrium parts-rule of
step (1) only to one distribution function which is orthogonal to the wall, and it
uses the remaining degrees of freedom to explicitly enforce the conservation laws.
One particular strength of the described model is that it can be implemented
on just any lattice, independently on the number of dimensions and the number of
distribution functions. This is a striking feature that distinguishes the regularized
boundary condition from other approaches. It is for example reected in the LB
source code of Ref. [13] in which a generic implementation for the lattice boundaries
exists. An application programmer can for example chose to dene a new lattice
structure and immediately apply it to a given problem, as both the LB dynamics
and the boundary conditions are generic and adapt to the new lattice. Nevertheless,
to help with the implementation of typical codes, explicit formulas for the tensor
(1) are listed for the D2Q9 and the D3Q19 lattice structures (those two structures

Boundary and initial conditions

57

are specied in Section 1.3.3 on page 11). On a D2Q9 lattice, the following value
are obtained:
(1)

00

= f2 + f6 + 2(f1 + f7 )
u1
1

+ u2 +
,
0
3
3

(1)

= 2 (f1 + f7 + f8 ) u1 + u2 +
1

(1)

= 2 (f7 f1 )

11
01

(5.4a)
1
3

u0
+ u0 u1 .
3

(5.4b)
(5.4c)

The value of (1) on a D3Q19 lattice is:

(1)

00

(1)

11

(1)

22

(1)

01

(1)

12

(1)

02

5.2.4

= f1 + f2 + f3 + f5 + f6 + f7 + f8 +
1
u1
+ u2 +
,
2(f10 + f12 )
0
3
3
= 2(f9 + f10 + f11 + f12 + f13 )
1
,
u1 + u2 +
1
3
= f2 + f4 + f5 + f6 + f7 + f8 + 2(f11 + f13 )
u1
1
+ u2 +

,
2
3
3
u0
= 2(f10 f12 )
+ u1 u0 ,
3
u2
= 2(f11 f13 )
+ u1 u2 ,
3
= f5 f6 + f7 f8 u0 u2 .

(5.5a)

(5.5b)

(5.5c)
(5.5d)
(5.5e)
(5.5f)

Numerical verication

The quality of the regularized boundary condition is tested on a numerical implementation of the 2D Kovasznay ow described in Section 4.4. Four dierent types
of boundary conditions are implemented for this ow:
Local regularized: The regularized boundary condition, in which the tensor (1)
is computed from local values on the boundary node.
Non-local regularized: The regularized boundary condition computed from the
strain rate S, for which a nite dierence scheme on neighboring nodes is
applied.
Zou-He condition: A commonly used local boundary condition described in Ref. [38].
Inamuro condition: A commonly used local boundary condition described in Ref. [35].
The Reynolds number in all simulations is xed to a value Re = 20, and the grid
resolution is progressively rened to observe how the accuracy scales with the discrete time step x . Taking into account the observations of Section 2.4, the time step

58

Chapter 5

10

Relative error

10

10

Regularized local
Regularized nonlocal
ZouHe
Inamuro

10

10

10

10
Grid resolution N

Figure 5.2: Accuracy of a 2D Kovasznay ow, implemented with four dierent

boundary conditions. The curves corresponding to the Zou-He and the Inamuro
methods are incomplete, because those methods produce numerical instabilities on
coarse grids.

2
t is rescaled so as to keep constant the ratio t /x = 0.42. According to Note 2.7
on page 29, this amounts to xing the relaxation parameter to = 1.7762. The
results of the benchmark are presented on Figure 5.2 on simulations with a varying
grid size. It is good to remember at this point that the value N = 1/x does not
correspond to the resolution of the complete simulated domain, but rather to the
resolution of a reference distance, corresponding to the spacing between to points
of the grid (the physical grid represented in the Kovasznay ow, not the numerical
grid). The simulations are implemented with a BGK dynamics instead of RLB in
order to respect the spirit in which the Zou-He and the Inamuro boundary conditions were initially developed. As it is expected from the discussion of the accuracy
of LB methods in Section 2.4, the accuracy of all simulations scales at a second order
rate with the grid resolution. This conrms on the one hand the choice of how to
rescale the time step t properly, and on the other hand the fact that all boundary
conditions exhibit an accuracy that is compatible with the accuracy of the LB dynamics. Although the convergence rate is the same, the Zou-He and the Inamuro
boundary conditions are somewhat more accurate than the regularized one. This
dierence can be understood by the fact that the benchmark simulations implement
a BGK dynamics, whereas the regularized boundary condition cuts o all higher
order modes of the dynamics. This dierence between the boundary conditions
vanishes when they are applied to a regularized LB simulation. It is furthermore
remarked that the Zou-He and Inamuro boundary conditions are substantially less
stable than the regularized one. Simulations using those boundary conditions are
subject to numerical instabilities, whenever the grid is too coarse or the Reynolds
number too high.

Boundary and initial conditions

59

The regularized boundary condition is hence found to be an excellent alternative

to boundary conditions traditionally used in the LB method. It exhibits the desired
second order accuracy, it is numerically more stable and, a point that must be
particularly lined out, it is generic and works with all types of lattice structures.
The accuracy of the simulation is somewhat better when a non-local version of the
regularized boundary is used. The local version delivers however good results too,
and it is often preferred for conceptual reasons.

5.3
5.3.1

Initial condition
Introduction

The impact of initial conditions on the evolution of a time-dependent ow can be

just as crucial as the choice of an appropriate boundary condition. This issue is
illustrated in the present section with the help of an interesting 2D benchmark
problem, which is taken from the recent literature and describes the turbulent interaction between two vortexes and a no-slip wall. This problem was originally
devised as a benchmark case to test the quality of boundary conditions in CFD
software. It is however shown that when using the LB method, the real diculty
stems from an appropriate initialization of the particle distribution functions. It is
especially crucial to compute a correct initial value for the pressure and to initialize
the equilibrium contribution to the distribution functions properly. The numerical
experiment is explained in some detail, and the proposed macroscopic approach to
the implementation of the initial condition is compared with an approach described
in the literature, based on lattice kinetic considerations. The benchmark application
is also used to show the limits of the LB method for incompressible ows due to
constraints on the discrete time step.
In uid turbulence, the interaction between vortexes and no-slip walls is known to
have an important inuence on the evolution of a ow. Numerous numerical studies
have shown that the wall acts like a source for small-scale structures, which on their
hand interact with the uid and aect its overall characteristics. In the numerical
experiment used in this section, a complex of two counter rotating vortexes is set up
and used as an initial condition for the simulation. The ow is incompressible, and
it is conned within a quadratic box with no-slip walls. The rotation of the vortexes
induces a resulting forward momentum, and the dipole is self-propelled towards
one of the walls. During the rebound that follows, various additional vortexes are
generated and detach from the wall. They interact on their turn with the original
vortexes and form, among others, systems of secondary dipoles.

5.3.2

The benchmark

The benchmark is based on an incompressible uid, described by the Navier-Stokes

equation in Eq. (1.10) on page 8, together with a continuity condition in Eq.(1.9).
The ow is conned in a box of dimension [1, 1] [1, 1] with no-slip walls. The
initial condition is chosen so as to dene two counter-rotating monopoles, one with
positive core vorticity at the position (x1 , y1 ) and one with negative core vorticity

60

Chapter 5

Figure 5.3: Vorticity in the initial conguration. The system is self-propelling and
is going to move towards the right wall, as it is suggested by the arrow.

at (x2 , y2 ). This is achieved with the following velocity eld u0 = (u0 , v0 ):

1
1
u0 = |e | (y y1 ) exp (r1 /r0 )2 + |e | (y y2 ) exp (r2 /r0 )2 , (5.6)
2
2
1
1
2
v0 = + |e | (x x1 ) exp (r1 /r0 ) |e | (x x2 ) exp (r2 /r0 )2 ,(5.7)
2
2
where ri = (x xi )2 + (y yi )2 , the parameter r0 labels the width of a monopole
and e its core vorticity.
The average kinetic energy of this system at a given time is dened by the
expression
1 1 1 2

E(t) =
|u| (x, t)d2 x,
(5.8)
2 1 1
and the average enstrophy by
1

(t) =
2

2 (x, t)d2 x,
1

(5.9)

where = x v y u is the ow vorticity.

In all presented simulations, the initial kinetic energy is xed to E = 2, by

requiring e = 299.5286. Furthermore, the Reynolds number and the monopole
radius are set to Re = 625 and r0 = 0.1. The monopoles are aligned symmetrically
with the box, in such a way that the dipole-wall collision is frontal and takes place
in the middle of a wall. The position of the monopole centers is (x1 , y1 ) = (0, 0.1)
and (x2 , y2 ) = (0, 0.1).
The enstrophy eld for the initial conguration is shown on Fig. 5.3. The upper
monopole has a positive, and the lower one a negative core vorticity.
The ow is simulated by means of the RLB method, implemented on a D2Q9
lattice.

Boundary and initial conditions

(a)

61
(b)

2000

1800
1600

enstrophy

kinetic energy

1400

1.5

1200
1000
800
600
400
200

0.5
0

0.5

1
t

1.5

0
0

0.5

1
t

1.5

Figure 5.4: Time evolution of the kinetic energy (a) and the enstrophy (b)

5.4
5.4.1

Numerical results
Time evolution of the ow

Figure 5.4 displays the time evolution of the energy and the enstrophy of the system,
with the numerical parameters dened in the previous section. The rst rebound
of the dipole with the wall happens approximately at a time t = 0.37. This event
is represented by a peak in the enstrophy curve, due to the generation of numerous
small-scale vortexes. At the same time, the slope of the energy curve becomes
steeper, as the energy dissipation rate is directly dependent on the enstrophy. Two
smaller peaks in the enstrophy curve indicate successive rebounds of the dipole with
the wall. The value of those maxima, and the time at which they occur, are used as
benchmark values to test the accuracy of the simulation. A chosen number of results
have been simulated by a nite dierence and a spectral method and reported in
Ref. [40].
In order to illustrate the physical processes occurring during the rebound of the
dipole on the wall, Figure 5.5 shows a contour plot of the vorticity at four chosen
time steps. At t = 0.30, the lower monopole is observed to be approaching the wall.
At t = 0.34, small scale structures are generated in the boundary layer close to the
wall. At t = 0.38, a vortex with positive core vorticity detaches from the wall and
creates a secondary dipole with the monopole that initially collided with the wall.
At t = 0.38, this secondary dipole turns over and prepares for a second collision
with the wall.

5.4.2

Benchmark values

In order to complete the discussion of this benchmark problem, the obtained numerical results are shortly presented. The chosen approach for the implementation
of the LB simulation is summarized in the next two sections.
The choice of a suciently small grid spacing depends on the width of the boundary layer, as the smallest relevant structures are located within this layer. The
boundary layer scales like the inverse square root of the Reynolds number, and so
does consequently the size of a grid cell. The resolution required to obtain grid

62

Chapter 5

Figure 5.5: Vorticity contours at four chosen time steps for the collision of the lower
monopole with the right wall. The shown system is an extract of the full simulated
space.

Boundary and initial conditions

63

convergence with the LB method is practically identical with the one described in
Ref. [40] on a FD scheme. At the given number of Re = 625, the system size is
for example 700 700. This is surprising, because the LB model is much more
straightforward than the FD model, which makes use of a multi-grid technique. Additionally, the LB model is based on a simple bounce-back model for the simulation
of the no-slip wall, which contrasts with the sophisticated boundary condition used
in Ref. [40]. The following table shows the computed position and value of the rst
enstrophy peak with the LB, FD and spectral methods:
Position (time)
Value (enstrophy)

5.5

Lattice Boltzmann
0.371
933.8

Finite Dierence Spectral Method

0.371
0.371
932.8
933.6

Setup of the initial condition

The initial value for the velocity eld is described by Eqs. (5.6) and (5.7). It is not
straightforward to set up such an initial condition in a LB simulation. This is because
the simulated degrees of freedom in this method are not the macroscopic variables p
and u, but rather the so-called particle distribution functions fi , with i = 0..8. The
particle distribution functions can however be related to the macroscopic elds and
their space derivatives throught the formula of the regularized LB method, Eq. 4.10:
(0)

(1)

fi = fi (p, u ) + fi (p, u ).

(5.10)

Given that the ow is considered in its incompressible limit, the pressure p and the
density are identied through the ideal gas law, Eq. (1.8) on page 7.
Two approaches to setting up an appropriate initial condition are listed in the
following. The rst approach has been cited in the literature [41] and shares the
philosophy of LB methods. It proposes to obtain the initial condition iteratively
by implementing a common LB dynamics. During those iterations, the local uid
velocity u is however not recomputed on ground of the simulated particle distribution
functions, but it is rather kept at the value one wants it to be in the initial condition.
In this way, only the pressure p and the rst-order contributions f (1) to the particle
distribution functions are free variables. They correspondingly converge to their
appropriate value for the initial condition.
The second approach, inspired from a macroscopic approach to uid dynamics,
consists in computing the pressure p and the strain rate S directly. The velocity
gradients are e.g. computed using a nite dierence scheme, and the pressure by
solving an iterative Gauss-Seidel scheme for the Poisson equation, Eq. (1.11).
Both approaches to setting up the initial condition require an iterative algorithm
to be implemented. Their eciency is compared in Fig. 5.6 (a). In this comparison, the LB approach was implemented in terms of BGK iterations, and the
macroscopic approach in terms of a successive over-relaxation (SOR) scheme. It
is obvious that the SOR scheme converges at a much faster rate, with a gain in
eciency that implies several orders of magnitude. For the sake of completeness, it
must be pointed out that the scheme suggested in the Ref. [41] is based on a MRT
approach to LB, in which some relaxation parameters are ne-tuned so as to speed
up the convergence. It is nevertheless clear that a dedicated method to solving the

64

Chapter 5
(a)

(b)
2

10

Relative error on peak enstrophy

Relative L2 error of the pressure

10

10

10

10

10

BGK convergence
Overrelaxed GaussSeidel

10

10

2000

4000
6000
Iteration steps

8000

10000

10

10

10
Time step t

10

Figure 5.6: (a) Comparison of two iterative schemes for the setup of the initial
condition: BGK iterations and SOR scheme (log-log plot). The rate of convergence
of the SOR scheme exceeds the one of the BGK approach by several orders of
magnitude. (b) Relative error on the value of the enstrophy at the rst enstrophy
peak. The simulated values are compared with the numerical reference data obtained
with a spectral method, and plotted for various choices of the discrete time step.

pressure equation is qualitatively and quantitatively faster than the resolution of a

lattice-kinetic scheme.

5.6

Discussion: choice of the time step

It is argued in Section 2.4 that the time step t of a LB simulation cannot be chosen
freely because it is coupled to the Mach number. This is a serious drawback of the
LB method, whose inuence is clearly visible in the present benchmark application.
In order for the LB uid to mimic an incompressible uid, the Mach number, and
consequently the time step, must be chosen very small. The actual time step required
in the LB method is typically one order of magnitude smaller than the one used in
the FD method.
Figure 5.6 (b) shows how the error on the position of the rst enstrophy-peak
decreases as a function of the discrete time step. To achieve the accuracy of three
signicant digits, a time step as small as t = 106 must be chosen. For comparison,
it can be noted that the time step adopted in Ref. [40] in a FD approach takes the
value t = 6.25 105 . One is therefore tempted to conclude that the error plotted
in Fig. 5.6 (b) is dominated by contributions due to uid compressibility.

Chapter 6
Adaptive space and time steps

6.1

Introduction

In the LB method, the discrete space and time steps x and t are constants that
can vary from one simulation to another but are xed within a simulation. The
grid spacing is for example the same in all space directions, and it does not vary
from one position of the discrete space to another. The same can be said for the
time step t , which does not change during the time evolution of the system. Those
assumption were implicitly taken for granted in the theoretical analysis of the LB
method in Chapter 2, and they must hold in order for the simulation to asymptotically solve the associated PDE. Some workarounds to this have been proposed in
the literature, where a mapping between a regular and an irregular mesh is obtained
via interpolation and extrapolation schemes. These techniques will however not be
considered further, as they represent hybrid constructs rather than pure LB models
according to the line of thought of this thesis.
Other types of CFD methods are not subject to such restrictions on the regularity
of the time and space discretization. Finite dierence methods, presented e.g. in
Ref. [4], are principally executed on regular grids with xed spacing as well. The grid
can however be anisotropic, that is, the value of the grid spacing can dier between
from a space direction to another. Furthermore, the discrete time step can change
during the evolution of the system. Finite volume methods [5] and nite element
methods [6, 5] are even more general and can be implemented on unstructured grids:
the grid points can be situated on arbitrary space positions, and it is not necessarily
possible to represent the relative position of grid points to each other by a matrix
data structure.
The use of inhomogeneous grids makes sense when a problem with an inhomogeneous geometry is being simulated. Some parts of the simulated domain require
a higher grid resolution than others in order to reach the required level of accuracy.
The grid resolution needs for example to be increased close to an obstacle with complicated shape to ensure that the discretized version of the obstacle resembles the
original one suciently well. In other cases, the resolution needs to be increased
because the uid ow exhibits small-scale patterns, such as the small vortexes generated close to the wall in the numerical experiment of Section 5.3. In simulations
with xed-sized grids, the overall resolution needs to be adapted to the maximal
required value, by which much computational eort is wasted. It should be pointed

66

Chapter 6

out that approaches to calculate a local a priori estimate of the numerical error
exist in many numerical methods, as those described in Ref. [6]. With this estimate,
the grid can be dynamically adapted during a simulation to reach the desired tradeo between accuracy and eciency. The LB approach to CFD has however been
developed quite recently, and no formal framework exists currently that would lead
to such an error estimate in LB simulations.

An adaptive time interval is likewise useful to adjust the time resolution, depending on how rapidly the ow structures vary in the simulation. In some methods
such as the one presented in Ref. [4], the time step must furthermore be chosen on
ground of well known stability criteria. Again, experimental evidence shows that LB
simulations are also numerically unstable when the discrete time step is too large,
but a theoretical framework is lacking from which exact stability criteria could be
deduced.

In order to overcome the limitation of a xed grid spacing in the LB method

and to adapt the simulation at least roughly to the geometry of a problem, a socalled multi block technique is often used. The idea is to partition the domain of
interest into rectangular subdomains, for each of which a dierent grid is used with
its own parameters x and t . It is obviously necessary to transfer data from one of
those grids to another during the execution of a simulation. Two diculties need
to be overcome in order to do so. First, the data on the interface between two grids
needs to be interpolated, because the dierent grids overlap only partially, as well in
space as in time. Second, the data needs to be rescaled to account for the changing
value of x and t . As a matter of fact, it was shown in Section 2.4 that the LB
variables are not dimensionless with respect to a macroscopic system of units and
thus depend on the particular value of the discretization parameters. When the
RLB model is used, it is however easy to rescale the variables of interest. They can
be reduced to a set of variables , u and (1) which are cast into a dimensionless
representation according to Eq. (4.11) on page 46. An application of this approach
to multi block LB modeling is presented in Section 6.2. The numerical example
presented in this section illustrates how one can compute the drag force acting on
an obstacle immersed in a uid with innite extent. A system of successively rened,
nested grids is used to account for the need of a high resolution close to the obstacle.
By combining this technique with a novel approach for the implementation of open
boundaries, which has recently been described in the literature, benchmark results
of an exceptionally good quality are obtained.

In Section 6.3, the RLB model with adaptive time stepping is introduced. It is
used to simulate the values of an incompressible stationary ow. A large time step is
chosen initially, so as to rapidly converge to the stationary state, and a smaller time
step in a later stage, so as to decrease the numerical error due to the compressible
nature of the uid model.

Adaptive space and time steps

6.2
6.2.1

67

Grid renement
Introduction

The multi block grid renement technique of the regularized LB model discussed in
the previous section is now tested on a numerically challenging task: simulating a
stationary incompressible uid ow past a rigid obstacle. The uid is assumed to
ll the whole 2D space, the obstacle is placed at the origin of the coordinate system,
and the uid velocity is asymptotically constant far from the obstacle. Thus, the
problem consists in the resolution of the stationary Navier-Stokes equation for the
velocity eld u(r) with boundary condition
lim u(r) = u0 .

|r|

(6.1)

Apart from the grid renement procedure, a major diculty stems from the necessity to truncate the innite domain for numerical purposes, and to nd an articial
boundary condition for the boundaries of the truncated domain. A straightforward
approach consists in using the asymptotic condition u = u0 on the numerical domain
boundaries. Although this method is easy to implement, it appears to be quite inappropriate for the needs of numerical modeling, as it requires the use of excessively
large domains. Indeed, it will be shown in the present section that the structure of
the ow is strongly inuenced by the shape of the obstacle even far from the center.
Other approaches to this problem use extrapolation schemes on the boundaries so as
to ensure a vanishing gradient perpendicular to the boundaries, for the velocity or
other physically relevant quantities. The drawback of those approaches is that they
are unable to properly tune the asymptotic velocity u0 in the uid and can therefore not be used on all boundaries. Furthermore, they make it dicult to ensure
conservation of mass and momentum across domain boundaries.
For those reasons, an alternative technique is introduced here that has been
recently described in the literature [42, 43, 44]. In this method, an explicit vector
eld is proposed that can be used to implement a Dirichlet boundary condition for
the uid velocity in a region reasonably far from the domain center. The expression
for this vector eld is obtained from a truncated asymptotic expansion of a solution
to the stationary Navier-Stokes equation and approximates the structure of the
ow with considerably higher precision than the constant approximation u = u0 .
The drawback of this method is that it depends on the drag and lift coecients of
the obstacle which are a priori unknown. Therefore, the solution process involves
a series of iteration steps during which the formula of the boundary condition is
updated on ground of the drag coecients of the obstacle measured at this state
of the simulation. A brief overview of the method is found in the next section,
Section 6.2.2.
The RLB method, executed on a set of progressively rened grids, is used to
solve this problem. The simulations are run at very low Mach numbers in order to
approach the nature of an incompressible uid with sucient accuracy. The following sections present a case study for the numerical evaluation of a drag coecient,
and serve three main purposes. First, they present an introduction to the boundary
condition of Refs. [42, 43, 44] and demonstrate its eciency and simplicity in the
context of LB simulations. Second, it is shown that this method can be coupled

68

Chapter 6
(a)

(b)
0.02

0.025

0.018
0.016

0.02

0.014
0.012

0.015

0.01
0.008

0.01

0.006
0.004

0.005

0.002

Figure 6.1: (a) Flow prole close to the obstacle when a zeroth-order boundary
condition u = u0 is used. (b) Velocity prole with the rst-order boundary condition
introduced in this section.

with a numerical technique based on iterative grid renement. The grid renement
approach of RLB, in which all variables are cast into a dimensionless representation
according to Eq. (4.11) is observed to produce excellent results. Finally, it is argued
that although the boundary condition of Refs. [42, 43, 44] has been developed for
incompressible ows, it also proves useful for simulations of compressible ows at
low Mach numbers. For further precision, the inuence of the uid compressibility
on the computation of a drag force is analyzed.

6.2.2

Analytical prole on the boundaries

In Ref. [42], the solution to the 2D incompressible Navier-Stokes equation is expanded in a nite series, as a function of formal parameters depending on the drag
and lift coecients of the obstacle. The corresponding theory for the 3D case is
presented in Ref. [43]. It is recognized that at a certain distance from the center,
the structure of the ow doesnt depend for the specic details of the obstacle geometry, but only on 2 (2D) or 3 (3D) distinct parameters. These considerations
result in the prescription of an explicit vector eld that can be used as a boundary
condition on the numerical domain boundaries. The zeroth- and rst-order terms of
the expansion for a 2D ow on the velocity eld u = (u, v) at a position r = (x, y)
read
d 1
y2
d x
b
y
+l 2
(x) l e 4lx
x2 + y 2
x + y2
x
d y
y2
y
b
x
d
l 2
(x) l 3/2 e 4lx
l 2
x + y2
x + y2
2 x

u(x, y) = u 1 + l
v(x, y) = u

,
(6.2)

where d = Fx /(2 l u2 ) and b = Fy /(2 l u2 ) are the drag and the lift coecient,

and l = /u is the viscous length, dependent on the dynamic uid viscosity .

The formula also uses the Heaviside function (x), dened as (x) = 1 for x > 0 and
(x) = 0 for x < 0. Without loss of generality, the asymptotic velocity u = (u , 0)
has been taken to be parallel to the x-axis. This boundary condition is implicit in

Adaptive space and time steps

(a)

69
(b)

Interface

Figure 6.2: (a) Structure of the numerical grid close to a rectangular obstacle. One
dot on the gure represents a square of 88 grid nodes. (b) Schematic representation
of the interface between two adjacent grids. In the actual application, the grids
overlap by additionally one node to allow for accurate interpolations.

the sense that it depends on two constants, d and b that are in general unknown
before the execution of the simulation.
The computed velocity prole close to the obstacle, depending on whether the
zeroth-order boundary condition u = u0 or the boundary condition dened by
Eq. (6.2) is used, is displayed on Fig. 6.1. The zeroth-order boundary condition
on Fig. 6.1 (a) produces quite unexpected results, as it imposes an unphysical zeroux condition through the domain boundaries. The rst-order boundary condition
on Fig. 6.1 (b) is free from such unphysical eects and produces, at least graphically,
the illusion of an innitely extended domain.

6.2.3

Numerical implementation

The RLB simulation uses a grid with high resolution close to the center, given that
in this region, the uid is subject to sharp pressure and velocity variations. A grid
renement technique is therefore applied with a hierarchy of nested grids that have
a successively ner resolution as they approach the system center. This hierarchy
of nested grids is schematically represented on Fig. 6.2 (a). At each level of grid
renement, both the space step x and the time step t are divided by two. This
choice was taken for a particular reason, because one purpose of this simulation
was to analyze the eect of uid compressibility on the accuracy of the result. For
this to be possible, the Mach number, which is proportional to the ratio t /x ,
needs to be the same from one grid to another. It should be mentioned however
that the most ecient choice for the simulation of incompressible uids consists in
2
keeping the ratio x /t constant. This topic is discussed extensively in Section 2.4.
During the streaming steps, values are transferred from one grid to another on the
common interface. To do this, they are rst cast into a dimensionless representation
according to Eq. 4.11 in order to be independent of the parameters x and t . Then,
the values are interpolated in space (because one grid has twice as many nodes as

70

Chapter 6

the other) and in time (because one grid iterates twice while the other one executes
only one iteration). The procedure adopted for those interpolations follows closely
the method proposed in Ref. [31]. In this reference however, the particle distribution
functions fi are directly interpolated, while in the present work, the interpolation is
applied to the dimensionless macroscopic variables , u and (1) . As an additional
dierence, it is remarked that Ref. [31] makes use of a rst order accurate time
interpolation scheme. Although the accuracy of the time interpolation is not crucial
here, as a stationary ow is being simulated, it is generally recommended to use
second order accurate time interpolations for consistency with the accuracy of the
LB method discussed in Section 2.4. Numerous other grid renement techniques
have been suggested in the literature that mostly adopt a microscopic point of view
to justify their approach to rescaling the particle distribution functions. Some of
those techniques can be found in Refs. [45, 46, 47, 48, 49].
The nested grids are deployed progressively in such a way as to speed up the convergence of the simulation towards a stationary state. In a rst stage, the simulation
is run on a small domain close to the obstacle. Then, at chosen time intervals, the
size of the physical domain represented in the simulation is enlarged by implementing a new, coarser grid. This is a convenient way of doing to eciently converge
towards a stationary state. Another approach to reach a rapid convergence would
consist in rening the time step t as shown in Section 6.3.
As it has been presented so far, the utilized numerical approach contains two
separate iterative processes, one for the implementation of the boundary condition,
as explained in Section 6.2.2, and one for the convergence to the stationary state,
as explained in the current section. It is observed that both processes take place at
comparable time scales, and can therefore be coupled in a simple manner. In the
present simulations, the drag and lift coecients needed for the implementation of
the boundary condition are updated each time the grid is enlarged.
The no-slip boundaries of the obstacle are implemented via a so-called bounceback boundary condition, introduced e.g. in Ref. [11]. The value of the force acting
on the obstacle is evaluated by computing the momentum which is exchanged on its
surface. It is easy to evaluate this momentum exchange on bounce-back boundaries,
on which particle distribution functions entering a node from the bulk are simply
reected. The exchanged momentum amounts to twice the value of the momentum
carried by the reected distribution functions.

6.2.4

Simulation results

For illustration purposes, the simulation results of a 2D ow across a rectangular

obstacle are presented. The ratio between the width and the height of the obstacle
is 5 : 1, as it is shown in Fig. 6.2. The simulations are run on quadratic domains of
varying size, up to three orders of magnitude larger than the obstacle. The Reynolds
number Re = A/l, dened with respect to the height A of the obstacle, is xed at

Re = 1, and the Mach number at M a = 0.02 3. On the domain boundaries, both

the constant boundary condition u = u0 , and the boundary condition described in
Section 6.2.2 are tested. The measured drag coecient d is plotted on Fig. 6.3 (a)
as a function of the domain size, compared to a reference solution d = 5.0268 that
was computed on a substantially larger domain. This gure shows that the accuracy

Adaptive space and time steps

(a)

71
(b)

10

10

Compressibility error

Relative error of drag

10

10

10

10

10

10

10

10

10
10
10
Diameter of computational domain

10

10
0.02

0.04
t / x

0.08

Figure 6.3: (a) Drag coecient as a function of the system size. The x-axis represents the ratio between the system size and the height of the obstacle. Solid line:
asymptotic boundary condition u = u0 . Dotted line: rst-order accurate boundary
condition described by Eq. (6.2). (b) Solid line: error on the drag coecient as a
function of the Mach number. Dotted line: reference curve, representing a power
law with slope 2.

of the drag coecient converges approximately at the same rate in the presence of
either boundary condition. However, using the boundary condition described by
Eq. (6.2) gains at least an order of magnitude on the accuracy of the drag coecient,
independently of the system size.
As a conclusion, the combination of a nested grid technique with an appropriate
boundary condition results in an ecient numerical scheme. Indeed, the computation of drag forces as those shown on Fig. (6.3) (a) were performed within one day
on a common desktop computer. In comparison, the same simulation would take
about one month to be executed on a homogeneous, small sized grid. It must also be
noted that at the high level of accuracy that is applied here, compressibility eects
play an important role. This statement is quantied by a plot on Fig. 6.3 (b), which
shows the dierence between the drag force exerted by a compressible uid and a
reference solution for an incompressible uid found in Ref. [44]. As it is expected,
this dierence decreases at roughly a second order rate with respect to the Mach
number. It is however not possible to implement uids at an arbitrary small Mach
number in the LB BGK model used here, as the choice of this parameter is limited
by CPU time constraints.

6.3
6.3.1

Adaptive time
Introduction

It is relatively simple to implement adaptive time stepping in the regularized LB

method. At a time step tn , the variables can be cast into their dimensionless representation using the previous time interval t1 = tn tn1 . After this, lattice units are
recovered again, based on the new time interval t2 = tn+1 tn . For the velocity, one
computes for example the dimensionless value u = x /t1 u and the rescaled value

72

Chapter 6

u = t2 /x u . Combining those two expressions yields u = t2 /t1 u. An adaptive

time interval can be used only to simulate incompressible ows, for which the Mach
number is phyiscally irrelevant. Indeed, the Mach number of a simulation is related
to the lattice parameters through the expression M a t /x and thus varies with
the time step.
Compared to other numerical techniques, the LB method is not very sensitive,
from the point of view of numerical stability, to the choice of t . Numerically stable
simulations are observed at impressively large values of the discrete time step. This
parameter is therefore rather chosen as a criteria related to numerical accuracy and
the limitation of compressibility eects. In this section, an adaptive time interval is
used to nd the solution to stationary ow problems as fast as possible. The initial
time step is chosen quite large in order to rapidly drive the system from its initial
condition to the stationary state. During the evolution of the system, the time step is
however progressively decreased so as to obtain accurate results. A similar approach
to accelerating LB simulations has been previously described in the literature under
the name of Mach number annealing. This technique is for example described in
Ref. [50]. Compared to this method, the novelty of the approach described in the
following sections consists in a proper rescaling of the particle distribution functions
at each modication of the time step. This is obtained through a regularization
procedure explained in Section 4 and ensures an accurate representation of the
time-dependent dynamics.

6.3.2

Numerical experiment

The stationary ow in a lid-driven cavity introduced in Section 4.4 is used again as a

benchmark application to verify the eciency of an approach based on an adaptive
time interval. Again, the accuracy of the numerical simulations is evaluated through
a comparison with high resolution reference results reported in the literature in
Ref. [37]. For this, the value of the velocity eld is compared with the reference
value at some chosen space positions along a horizontal and a vertical line passing
through the cavity center, as it is shown on Fig. 4.3 on page 51.

Note 6.1: Stationary ows and LB

In the present section, as well as in Section 5.2, incompressible stationary ow problems are considered. Those problems
are described by a time-independent version of the Navier-Stokes equation (1.10)
on page 8, in which the term containing a time-derivative is simply skipped. To
solve such a problem with the LB method, a time-dependent system is simulated
which progressively approaches a steady state. When the macroscopic variables
of the system do not change any more from one time step to the next, the system
is considered to have reached the time-independent solution to the Navier-Stokes
equation. This is usually not an ecient way of doing. Other iterative procedures
exist that nd solutions to stationary ow problems substantially faster. In those
methods, the intermediate steps of the iteration do not necessarily represent a
physically meaningful ow, by which shortcuts can be taken to reach the desired
nal result. This limitation of the LB method can be partially overcome with

Adaptive space and time steps

(a)

73
(b)

Figure 6.4: Convergence of a 2D lid-driven cavity ow towards the steady state,

with a lattice Boltzmann RLB method and a nite dierence method. (a) Fixed
2
2
time step t = 1/2 x in FD, xed time step t = 5/2 x in RLB. (b) Fixed time step
2
2
2
t = 1/2 x in FD; variable time step from t = 25 x to t = 5/2 x in RLB.

the help of the techniques described in this thesis. In Section 5.2, the problem is
rst solved on a domain with limited extent, surrounding the region of interest.
This gives a rst rough approximation to the solution, which in the following is
improved by progressively enhancing the extent of the simulation. In the present
section, a rst rough approximation is obtained via a simulation with a large
time step, and the accuracy is increased by a progressive decrease of the time
step.
In order to appreciate the eciency of the LB model, the ow problem is also
solved with the help of a nite dierence model which is taken from Ref. [4] and
described in more detail in Chapter 7. The FD solver has similar properties as the
LB model, as it also makes use of a xed-width regular grid. With both the LB
and the FD model, the accuracy of the solution is monitored during the evolution
of the simulation. Given that the emphasis is put on rapidly obtaining a stationary
solution, rather than analyzing the time evolution of the ow, the real wall-clock
time instead of the simulated physical time is measured. Figure 6.4 displays results
of a simulation run on a common desktop computer, with a Reynolds number of
Re = 100 and on a grid of size 129 129. When the time step t is xed, as on
Fig. 6.4 (a), both simulations converge at a comparable speed. The time step of the
FD method is xed by stability criteria, whereas in the LB method, it is chosen so
as to obtain results of equivalent accuracy in both methods. On Fig. 6.4 (b), the
RLB simulation makes use of an adaptive time interval, starting out at a value of
2
t = 25 x . This means that in the beginning of the simulation, the velocity of the
lid adopts a value of u0 = t /x 0.2. In the following, t is decreased at each
2
iteration step, until the nal value t = 5/2 x is reached. With this approach, the
convergence is accelerated by roughly an order of magnitude. Those comparisons
must of course be appreciated with some care, because the speed of a simulation
depends among others on the quality of its numerical implementation. One can

74

Chapter 6

however consider the FD and LB implementations used here to be of equivalent

quality, as they were both written by the same programmer in the C language, and
comparable eorts were made to optimize the execution speed of either code. As
a conclusion, it can be said that the speed up obtained by applying an adaptive
time interval is considerable, and this technique is obviously of high interest in the
practical work with the LB method.

Chapter 7
Coupling with other tools of CFD

7.1

Introduction

The topic of computational uid dynamics is introduced in Section 1.2. The overview
presented there shows that many dierent models for the numerical analysis of uid
ows exist. On the one hand there are solvers based on discrete representations of
the Navier-Stokes equation, most notably nite dierence, nite element or nite
volume methods. On the other hand, a new category of solvers has emerged over
the past decades which describe the uid at a molecular level, such as Cellular
Automata, or at an intermediate level, such as LB methods.
In this section, the possibility of solving separate spatial regions of a simulation
with a dierent solver is investigated. In this approach, a nite dierence (FD)
solver is coupled with a LB method. The motivation for developing such a type of
coupling is that, depending on the geometry of the ow, one technique can be more
ecient, less memory consuming, or physically more appropriate than the other
in some regions, whereas the converse is true for other parts of the same system.
One can also imagine that a given system solved, say by FD, can be augmented
in some spatial regions with a new physical process that is better treated by a
LB model. With the approach shown here, only the concerned region is modied,
without altering the rest of the computation.
In order to couple a LB with a FD model, it is crucial to understand how the LB
set of variables is related to the FD set and vice versa. As has often been the case in
previous chapters, the transformation from LB variables to macroscopic variables is
easy to perform, whereas the reverse transformation from macroscopic variables to
the particle distribution functions requires special attention. Also, as many times
before, it seems most natural to solve this problem via the regularized LB method,
in which the state of the LB system is expressed in terms of macroscopic variables
and their derivatives. In order to obtain the desired coupling, it only remains to
nd appropriate interpolations for the concerned variables on the interface between
a LB and a FD region.
The procedure is of course also applicable when the BGK method is used. In
that case, the formula of the regularized LB method, Eq. (4.10) on page 46, is used
to convert from macroscopic variables to the particle distribution functions, after
which the BGK collision step is executed. The numerical application presented in
Section 7.6 makes actually use of the BGK method. At the time this simulation was

76

Chapter 7
.

.
fi1,j
j

vi,j

fi,j

ui1,j

pi,j

ui,j

fi1,j1 vi,j1 fi,j1

FD variables
LB variables
i
.

Figure 7.1: Choice of indexes for FD and LB variables on a chosen grid cell (i, j).

prepared, the benets of using the regularized LB on the full domain was not yet
understood, and this method was used only on boundary nodes.

7.2

FD model and computational grids

The FD model used in the present section is taken from Ref. [4]. Its scheme is explicit
for the velocities and implicit for the pressure. In order to be suciently stable, the
method makes use of a staggered grid, that is, a grid on which all macroscopic
variables are dened on dierent space locations. The relative arrangement of those
variables, as well as their space position with respect to the LB set of variables, is
discussed in the following.
The full computational domain is considered a two-dimensional, rectangular region
= [0, l] [0, h] R2 ,
on which a regular grid is dened. This grid is divided into imax cells of equal size
in the x-direction and jmax cells in the y-direction, resulting in grid lines spaced at
a distance
x = l/imax and y = h/jmax .
Although the FD model is potentially anisotropic and can handle cases in which
x is dierent from y , this is not true for the LB method. Therefore, those two
parameters are always taken equal in the following and are simply labelled x .
The FD model is based on three quantities that are dened at each cell position:
the pressure (p), the x-component (u) and the y-component (v) of the velocity. They
are however placed on a staggered grid. A given index (i, j) of the cell is assigned
to the pressure at the cell center, to the x-component of the velocity at the right
edge and the y-component at the upper edge (cf. Figure 7.1). The reason for this
staggered arrangement is that it prevents possible pressure oscillations which could
occur had all three variables u, v and p be evaluated at the same grid points.

Coupling with other tools of CFD

77

.
LB

FD
u
v
u p u
v
v
u p u
v
v
u p u
v
u

u
v
p u
v
p u
v
p u
v
u

v
p u
v
v
p u
v
v
p u
v

j=4
j=3
j=2
j=1
j=0

i=0 i=1 i=2 i=3 i=4

i=0 i=1 i=2 i=3 i=4

FD bulk node
FD boundary node

LB bulk node
LB boundary node

Figure 7.2: Computational grid representing a domain of size (imax x )(jmax x )

with imax = jmax = 3. The left hand side depicts the staggered arrangement of the
variables over the grid when the domain is resolved by a FD scheme. In the case of
a LB solver, all variables are located on cell edges, as shown on the right hand side.
The location of the boundary strip is indicated by a dashed line.

The LB model uses 9 variables fk , k = 0 8 which are all evaluated at the same
location of a cell. In the present coupling model, the LB variables are situated in
the upper right corner. This choice is of course arbitrary. It is however reasonable
to situate the f s on cell corners, to ensure that the LB and the FD model have
a compatible interpretation of the location of the domain boundary . Indeed,
this boundary is dened on a cell edge in the FD model. Considering that most
implementations of LB boundary conditions set the domain boundary on top of a
LB node, this leads to placing the LB node on the intersection of two cell edges.
The situation is depicted on Figure 7.2 for a system of extent imax = jmax = 3.
It shows also that as a result of the staggered arrangement, not all extremal grid
points of the FD set of variables come to lie on the domain boundary. For this
reason, an extra boundary strip of grid cells is introduced, so that the boundary
conditions may be found by linear interpolations between the nearest grid points on
either side.
The FD model is based on a discretization of the incompressible Navier-Stokes
equation (1.10) and the continuity equation (1.9). The computation of successive
iterations u(t) , v (t) , p(t) u(t+1) , v (t+1) , p(t+1) contains two distinct steps:
1. Resolution of the Poisson equation (Eq. (1.11) on page 8) to obtain the new
pressure eld. This computation utilizes the values of the pressure and the
velocity at the time t: u(t) , v (t) , p(t) p(t+1) . In presence of Dirichlet
boundary conditions, this procedure has a unique solution (except for an integration constant). In particular, there is no need to know the value of the
pressure on the boundary.

78

Chapter 7

2. Computation of the new velocity eld according to a nite dierence scheme.

It uses the pressure eld at step t + 1: u(t) , v (t) , p(t+1) u(t+1) , v (t+1) .

7.3

Choice of units

The LB method as presented so far makes use of a system of lattice units, in which
the distance between two adjacent grid points, as well as the time interval between
two successive iteration steps, is unitary. The FD method on the other hand uses the
dimensionless system introduced in Section 1.2.3, in which the simulated variables
are independent of x and t . In order to easily formulate the coupling between the
two methods, the variables of the LB method are also cast into dimensionless units.
The spacing between adjacent grid points thus takes the value x , and the interval
between successive iterations the value t . The lattice constants ci are replaced by
the lattice velocities vi = x /t ci , and the LB dynamics in Eq. (1.14) on page 10 is
reformulated as follows:
fi (r + t vi , t + t ) fi (r, t) = i .

(7.1)

Whenever some data is transferred from one grid to another, the dimensionless form
, u and (1) of the variables in the RLB model is used. The stars are however
omitted to simplify the notation.

7.4

The FD-LB interface

The full domain is now split into two subdomains 1 and 2 such that =
1 2 . In domain 1 the FD method is active and in 2 the LB method is used.
For simplicity, the subdomains are taken to be rectangular, 1 occupying the left
hand side and 2 the right hand side of the domain. This procedure can however
be extended with few changes to a general boundary.
Figure 7.2 displays the position of extremal grid points, drawn as white circles
and squares, on which a boundary condition needs to be implemented to achieve a
coupling between the two grids. On the interface between 1 and 2 , those boundary
values are taken from the boundaries of the opposite domain. As a consequence of
the staggered arrangement of the LB values with respect to the FD values, there is
need for an overlap between 1 and 2 . There are several ways the coupling can
be implemented. In the approach chosen here, the two grids overlap by roughly one
lattice site. The position of this site is indexed by i = iint (see Figure 7.3).
As a conclusion, the FD domain requires knowledge of the values of ui,j for i = iint
and j = 1 jmax , and the values of vi,j for i = iint + 1 and j = 0 jmax . Those
values are easily obtained from the LB eld, on which the macroscopic variables are
computed by means of Eqs. (1.15) and (1.16) on page 10.
The LB domain on the other hand requires the knowledge of the 9 values fk;i,j
at i = iint 1 and j = 0 jmax . To obtain those, the regularization formula from
Eq. (4.10) on page 46 is used, which connects the pressure, the velocity and the the
velocity derivatives to the value of the particle distribution functions.

Coupling with other tools of CFD

79

.
i = iint

FD
v
v
v
v
u p u p u p u
v
v
v
v
p u p u p u
u
v
v
v
v

FD bulk node
FD boundary node

LB

LB bulk node
LB boundary node.

Figure 7.3: Subdivision of the computational domain into a FD subdomain (1 )

and a LB subdomain (2 ). One lattice cell on the interface between 1 and 2 ,
at the position i = iint , is computed by both methods. The boundary nodes of
the subdomains are represented by white symbols. They must be implemented by
means of a coupling term between the two methods.

7.5
7.5.1

The coupling algorithm

Coupling the pressure

Before turning to the velocity eld, the treatment of the pressure in the FD simulation and the density in the LB simulation is discussed. It is useful to remember that
the LB simulation is presently considered in an incompressible regime, in which the
density is connected to the pressure via the ideal gas law in Eq. (1.8) on page 7.
The pressure eld of an incompressible ow is not uniquely dened. It contains
a constant oset which can be chosen arbitrarily. In the FD method used here,
the value of the oset depends on how fast the numerical method converges, and it
varies from one time step to another. It is therefore dicult to couple the pressure
eld of the FD and LB simulations together. To do so would require a calibration of
the pressure elds by computing the average pressure on each side of the simulation.
It is fortunately possible to consistently compute the value of the pressure on each
computational subdomain, without transferring values from 1 to 2 . As a matter
of fact, the value of the pressure in an incompressible ow can be directly deduced
from the value of the velocity eld through the Poisson equation (1.11). Practically
speaking, on the side of the FD simulation, the FD-LB interface is simply considered
as a domain boundary, as far as the pressure is concerned, and the Poisson equation
is solved on the subdomain 1 . Similarly, the pressure of the LB simulation is
computed by applying a usual boundary condition on the FD-LB interface as it is
described in Section 5.2.2.

7.5.2

Coupling the velocity

Obtaining values for the velocity in the FD eld from the LB eld is a simple
exercise. The two components of the velocity are computed in a straightforward

80

Chapter 7

manner from Eqs. (1.15) and (1.16) on page 10. Because the FD and the LB variables
are not dened at the same point of a lattice cell, we need to adjust the values by
an interpolation which, for consistency with the accuracy of the overall numerical
scheme, is second order accurate. This leads to the following relations:
uFD,j = 1/2 uLB ,j + uLB ,j1
iint
iint
iint
FD
LB
LB
viint +1,j = 1/2 viint +1,j + viint ,j

(7.2)
(7.3)

A similar procedure is used to translate the values of the velocity from the FD
eld to the LB eld:
uLB 1,j = 1/2(uFD1,j + uFD1,j+1 )
iint
iint
iint
LB
FD
FD
viint 1,j = 1/2(viint 1,j + viint i,j )

7.5.3

(7.4)
(7.5)
(7.6)

Coupling the velocity gradients

The regularization formula (4.9) refers to the tensor S which, in Eq. (1.6) on page 7,
is dened in terms of velocity gradients. Figure 7.2 shows that two of those gradients
can be approximated by a centered dierence of half the mesh width:
y uLB 1,j = 1/x (uFD1,j+1 uFD1,j )
iint
iint
iint

(7.7)

LB
x viint 1,j

(7.8)

FD
1/x (viint ,j

FD
viint ,j1 )

One gradient needs to be calculated as a centered dierence of mesh width, based

on interpolated values of the velocity:
LB
y viint 1,j =

1 FD
FD
FD
FD
(v
+ viint 1,j+1 (viint ,j1 + viint 1,j1 )).
4x iint ,j+1

(7.9)

There are not enough grid points at hand for computing the fourth gradient at the
same level of precision. Luckily enough, this value can be computed with the help
of the continuity equation (1.9) on page 7:
LB
x uLB 1,j = y viint 1,j
iint

7.5.4

(7.10)

Summary

Now that all missing variables have been computed, it is useful to take a step back
and discuss the overall algorithm of a LB iteration step. For the purpose of this
discussion, the BGK dynamics is split into two steps. The rst, the collision step,
handles the computation of the equilibrium distribution and maps the incoming
(in)
(out)
particle stream fi onto the outgoing particle stream fi . It is followed by
a streaming step that transports the particles by a value of t in time and vk t in
space. The details of an iteration step are as follows:
1. On bulk nodes, (t) and u(t) are computed from the incoming particle densities
(in)
(in)
fk (t). On boundary nodes, all the values (t), u(t) and fk (t) are obtained
from the variables of the FD eld at time t.

Coupling with other tools of CFD

81
(out)

2. All nodes, bulk and boundaries, perform the collision step: fk

(in)
(eq)
)fk (r, t) + fk (r, t).
(in)

(r, t) = (1

(out)

3. The bulk nodes perform the streaming step: fk (r + vk t , t + t ) = fk

for all r such that r + vk t lies on a bulk node.

(r, t),

Alternatively, it is possible to extend the streaming step to boundary nodes for those
(in)
values of fk that are incoming from the bulk of the LB simulation. They are kept
(in)
unchanged, unlike the remaining set of fk and the macroscopic variables and u
that are provided by the FD eld. In the simulations, this procedure appears to
produce results equivalent to those of the proposed algorithm.

7.6

Numerical validation

The coupling algorithm is validated by the simulation of a Poiseuille ow. This

is a stationary ow in a channel of innite length with no-slip boundaries. An
analytical solution for this ow is shown in Eq. (3.7) on page 33. In lattice units,
the boundaries extend horizontally at a height y = 0 and y = L. The uid velocity
is strictly horizontal and does not depend on the x position: v = 0; x u = 0. In
the present example, the uid is driven by a body force. The left and the right
borders implement periodic boundary conditions in order to simulate a channel of
innite length. Special care must be taken when implementing this boundary in the
FD model. Indeed, during the simulation it tends to build up a pressure gradient
that must be eliminated by imposing a constant value of the pressure on the left
and right boundary.
The simulation is performed on a grid of size imax = 3 and jmax = 50. The
physical channel width is set to L = 1 and the body force has the value Fx = 0.01.
The numerical values usim are compared to the analytical solution by means of the
overall error:
jmax
2
j=0 (usim (, j) uana (, j))
=
(7.11)
jmax
uana (, j)2
j=0
Three simulations have been conducted that can be compared to each other.
The rst simulation implements a pure LB scheme with the BGK model, the second
simulation a pure FD model, and the third simulation is a FDLB hybrid. In the
rst case, the no-slip property of the walls is implemented by a boundary condition
known under the name of bounce-back (see e.g. [8]). In the third case, the top and
bottom strips of size jmax = 3 are computed by the FD model and the bulk domain
by the LB model (see Figure 7.4).
A remarkable result of the simulations is that the FD-only model reaches the
analytical solution at the machine level of precision (1015 ). Although there exist
LB boundary conditions which obtain the same result, as for example the one in
Ref. ([35]), their implementation is less natural and straightforward than the one
of the FD model. It is further remarked that the LB model converges faster (in
terms of iteration steps) than the FD model. The stationary velocity eld of the
LB simulation is however distinct from the analytical prediction, due to the limited

82

Chapter 7
.

3 (FD)
u

2 (LB)

(3)
jmax = 3

(2)
jmax = 46

L=1
jmax = 50
y

1 (FD)

(1)
jmax = 3

x.

Figure 7.4: The computational grid for the simulation of a Poiseuille ow is partitioned into three subdomains 1 , 2 and 3 . The FD scheme is used on the
boundary domains 1 and 3 , and the LB scheme on the bulk domain 2 .

precision of the boundary condition. The hybrid simulation shows the expected
convergence speed of the LB model and an error due to the limited precision of
the coupling. However, the error is two orders of magnitude smaller than the one
due to the bounce-back boundaries of the LB-only simulation. The results of the
simulations are shown on Figure 7.5.
The order of precision of the coupling can be estimated by varying the grid
resolution (jmax ) while keeping the physical quantities(L, Fx ) constant. Figure 7.6
plots the error of the stationary velocity eld as a function of the grid resolution.
It appears clearly that the coupling acts like a second-order boundary condition for
the velocity eld. No conclusion can be taken concerning the implementation of the
pressure eld, because the latter is constant in a Poiseuille ow.

7.7

Conclusion

In this chapter, a LB scheme for 2D incompressible uid ows is spatially coupled to

a FD scheme on a computational domain partitioned in two regions. The methods of
the regularized LB scheme are used to relate the LB distribution functions fi with the
classical physical quantities and their derivatives. A specic FD schemes has been
introduced, for which a complete coupling algorithm has been explicited. At the
interface, the LB and FD are connected so as to preserve continuity of the physical
quantities. The connection between the fi variables and the standard macroscopic
physical quantities is obtained through a dimensionless representation of the LB
variables via the regularized LB method: the equilibrium part of fk is related to
the macroscopic quantities and the non-equilibrium part to the gradients thereof. A
validation was performed by simulating a Poiseuille ow with FD boundary strips
and LB bulk and comparing it with an analytical solution. The simulation shows
that in this case, the coupling of the velocity eld is of second order in the grid
resolution.
It would be interesting to push this work further and test the results of the
coupling on more interesting geometries. It would also be interesting to extend the

Coupling with other tools of CFD

83

Figure 7.5: Simulation of a body-force driven Poiseuille ow with (1) a LB model,

(2) a FD model and (3) a FD-LB hybrid. The curves show the time-evolution of
the error, compared to the analytical solution of the Poiseuille ow.
0

10

LB with bounceback
FD
FDLB coupled
1

error

10

10

10

10

20

40

60
80
100
(iteration step)* t

120

140

160

Figure 7.6: Error of a FD-LB hybrid Poiseuille ow simulation as a function of the

grid resolution. A log-log plot shows the coupling of the velocity eld to be of second
order.
1

10

error

10

10

10

Numerical result
Slope 2

10

10

10

grid resolution (jmax)

10

10

84

Chapter 7

work presented here to other methods in CFD, such as the nite element method.
One key issue in this procedure would be to relate the values from the irregular
grid of typical nite element methods to the regular grid of the LB method. Such
questions are at the heart of the OpenLB project in Ref. [13], in which this topic is
approached from a software engineering point of view. The developed LB code is
encapsulated in object-oriented data structures that are very similar to those used
in a related nite element package, which encourages the use of those two tools on
common problems.

Chapter 8
A model for uid turbulence based on
kinetic variables

8.1
8.1.1

Turbulence modeling
Overview

In all numerical experiments conducted in the previous chapters, the simulated ows
have a laminar structure. At all time steps, it is possible to describe the velocity
eld exhaustively via a restricted number of streamlines. Those lines are dened
as integral curves to the ow eld: on every point of the curve, the ow velocity is
tangential to the curve. In a laminar ow, the streamlines vary smoothly in space and
time, and they are often used to illustrate the ow visually. In high Reynolds number
ows, a turbulent regime is observed to take place in which such a description is not
useful any more. Turbulent ows are characterized by the formation of structures
at many length scales. A description that focuses on a specic length scale, such as
a ow visualization through streamlines, captures only a partial aspect of its nature
and neglects other, physically relevant components to an exhaustive description of
its state.
The number of degrees of freedom required for the numerical simulation of ows is
substantially higher in presence of turbulent ows than with laminar ows. In order
to obtain an accurate representation of the ow evolution, all structures down to
the smallest ones must be representable on the computational grid. With currently
available computing resources, it is actually impossible to reach any of the Reynolds
numbers that are of interest in most engineering applications. To face this problem,
simulations are often simply run on an underresolved grid. They are known as large
eddy simulations, because they simulate the eect of the large eddies in the ow but
do not take into account the small ones. The eect of the small eddies is however
not neglected. Instead, a new, modied dynamics is executed on the computational
grid, which in the following is called the coarse grid. This dynamics replaces the
eect of the small eddies, had the simulation been executed on a fully resolved ne
grid, by an analytical or numerical prediction. This prediction is made on ground of
variables simulated on the coarse grid and is of course only an approximation to the
exact, fully resolved dynamics. The most common of those subgrid models include
the notion of an eective viscosity: they predict that the adapted dynamics on the
coarse grid implements the Navier-Stokes equation with a time and space dependent

86

Chapter 8

viscosity. A rough approximation, the so-called Smagorinsky model, relates the

eective viscosity to the value of the local strain rate tensor S. In more advanced
models, a new family of observables is introduced, for which a PDE is solved, and
to which the eective viscosity is related. It must be emphasized that all those
eorts are required only for the simulation of turbulent ows. When laminar ows
are simulated on an underresolved grid, an approximation to the exact solution is
obtained without the need for subgrid modeling. This is for example visible in the
numerical experiments conducted in Sections 2.4, 3.1 and 4.4, in which the accuracy
of the numerical result is observed to scale with the grid resolution. If on the other
hand subgrid modeling is omitted in turbulent ows, the numerical results dier
fundamentally from the ones of a fully resolved simulation, and, as it is most often
the case, numerical instabilities appear.
The aim of this chapter is to analyze whether the closure equations required for
the computation of the subgrid viscosity can be simulated at low cost in the LB
method, by exploiting the high number of degrees of freedom in typical LB lattice
structures. This new focus contrasts fundamentally with the point of view taken so
far in this document with respect to the LB method. In the regularized LB model
it is argued that, in order to solve the Navier-Stokes equation, it is sucient to
simulate the dynamics of some macroscopic variables, the density , the velocity
u and the strain rate tensor S. The set of q particle distribution functions fi is
thus shown to be overdetermined, and it is consequently reduced in the regularized
model. This is no longer the case in the present chapter, in which the original BGK
model is used, and the enhanced number of degrees of freedom is explicitly exploited
for the needs of subgrid modeling.
The next section contains a brief overview of large eddy models in classical
CFD tools. This is followed by the introduction of a new turbulence model, which is
directly based on the simulated variables of a LB simulation, the particle distribution
functions. The validity of this model is then tested with the help of numerical data
obtained from a fully resolved simulation of a turbulent ow.

8.1.2

Large eddy simulations in 3D incompressible uids

At high Reynold numbers, a uid enters a turbulent regime in which it is unsteady

and non-periodic. Turbulent ows are characterized by the occurrence of eddies
whose size may vary over a large range. The larger eddies contain the main portion
of the ow energy. This energy is successively transferred to the smaller eddies, and
is eliminated by viscous dissipation in the smallest ones. This process is described
by the theory of Kolmogorov [51]. It predicts that the size of the smallest eddies is
proportional to Re3/4 .
In a numerical simulation of a turbulent ow, the smallest eddies must be resolved
by the grid. Given three space dimensions, this requires N = O(Re9/4 ) grid points in
the discretization. This straightforward approach, involving the discretization of the
grid suciently ne for resolving all occurring eddies, is known as direct numerical
simulation (DNS). In industrial applications such as aerodynamic investigations
of automobiles, typical Reynolds numbers are found at a value of 106 and above.
Hence, solving these types of problems properly would require over 1013 grid points, a
number which is currently inaccessible in term of storage space or CPU performance.

A model for uid turbulence based on kinetic variables

87

For this reason, the emphasis in turbulent simulations is put on the macroscopically observed mean values rather than the microscopic detail. The physics of the
smaller, subgrid scales is no longer simulated, but replaced by an appropriate model.
Along this line of thought the velocity eld u is split into a mean part U = u and
the small variations u known as turbulent uctuations, such that
u=U +u.

(8.1)

The mean part is, for example, a componentwise statistical, temporal or spatial
average. The operator is often called a lter.
The lter can be applied to the Navier-Stokes equation, Eq. (1.10) on page 8.
The corresponding averaged equation is known under as the Reynolds equation for
the averaged velocity U :

U+ U
U =
P +
t
0

U+

uu ,

(8.2)

where P = p is the averaged pressure. The continuity equation, Eq. (1.9) becomes

U = 0.

(8.3)

Equation (8.2) is almost equivalent to the Navier-Stokes equation, except for an

additional term depending on the Reynolds stress tensor (u ) dened as
(u ) = u u .

(8.4)

A turbulence model expresses the Reynolds stress tensor in terms of the resolved
mean quantities and potentially some new variables described by a set of additional
equations. Widely used and easy to apply is the Smagorinsky approach [52], in
which the Reynolds stress tensor is dependent only on the local strain rate. This
model is extremely convenient in numerical simulations as it does not introduce
new degrees of freedom and does not require the resolution of additional PDEs. It
shows however good results only in two dimensions and for very ne mesh widths.
An improved version is the very popular k model [53, 54] that describes the
Reynolds stress in terms of the strain rate, the turbulent kinetic energy k and the
dissipation rate :
1
|u |2
and
2

=
| u + u |2 .
2

k=

(8.5)
(8.6)

The value of those additional variables is determined by two additional transport

equations. The Reynolds stresses are approximated through the following term:
(u ) R(S, k, )

2
k2
k I 2 c S.
3

(8.7)

This expression for the Reynolds stresses is substituted into Eq. (8.2), and a renormalized viscosity is introduced:
= + T = + 2c

k2

(8.8)

88

Chapter 8

as well as a renormalized pressure:

2
P = P + PT = P + 0 k.
3

(8.9)

This leads to the following averaged momentum equation:

1
U+ U
U =
P + ( S) .
t
0

(8.10)

This equation, as well as the continuity equation and the closure equations for
k and can now be simulated on a coarse grid with reasonable requirements for
computational resources. The constant c , as well as some additional constants
appearing in the closure equations, cannot be found from theoretical arguments and
must be calibrated on numerical experiments. The k model has also successfully
been used together with the LB scheme, as has been reported in Ref. [55].
Appropriate turbulence models are dicult to implement in practice. A major
diculty stems from the treatment of the domain boundaries, along which boundary
layers are observed, which dont obey the universal laws of the Kolmogorov theory
and need to be treated separately. An overview on the basics of turbulence modeling,
and the key issues for the simulation of turbulent ows, is found in Ref. [56].

8.2

Averaged LB equation

The basic idea developed in this section consists in applying the averaging procedure
described in Section 8.1.2 directly to the variables of a lattice Boltzmann simulation, the particle distribution functions fi . The following calculations are restricted
to the BGK model and the the 3-dimensional D3Q19 lattice (see Section 1.3.3).
The concepts are extended to other lattices in a straightforward manner. The LB
dynamics can be viewed as a map that applies the family of particle distribution
functions dened on the computational grid at the initial stage t0 on their state at
a time t = t0 + T :
DT : f (t0 ) f (t0 + T ).
(8.11)
The dynamics of the BGK model is governed by a relaxation towards an equilibrium
term, described by Eqs. (1.14) and (1.17) in Section 1.3.2. It is described as a
relaxation
fi (r + vci , t + t) fi (r, t) = (fi (r, t) fieq (u(r, t), (r, t)))

(8.12)

The equilibrium term is explicited in Eq. (1.18), and it is useful to repeat it here:
fieq (u(r, t), (r, t)) = ti (1 + 3 ci u +

9
Q : uu).
2 i

(8.13)

In this equation, the generic lattice constants have been replaced by their value on
the D3Q19 lattice. The velocity u, specied by Eq. (1.16), depends linearly on the
particle ow densities fi . Furthermore, the ow is considered incompressible, which
is the case at low Mach numbers (see Section 2.4). In that case, the non-linearity
of dynamics is entirely captured in the last term of equation 8.13. In this term, the

A model for uid turbulence based on kinetic variables

89

2
velocities are projected on the second order base vectors C , dened by means of
the tensors Qi specied in Eq. (1.19).
In order to formulate a turbulence model for the LB dynamics, a lter
is
introduced which acts as a linear map on the distribution functions f :

: f (t) F (t).

(8.14)

The dynamics of the ltered population

DT : F (r, t) F (r, t + T )

(8.15)

is dened as follows, using the functions composition operator :

DT

DT .

(8.16)

Equation (8.16) has already been used in Ref. [57] with the hope that a turbulence
model can be formulated from the obtained expression. An interpretation of this
procedure from the point of view of kinetic theory is presented in [58], and previous
work on subgrid models in the framework of kinetic theory is found in Refs. [59, 60].
The main inconvenience of an approach in which the lter is applied individually
to each of the particle distribution functions stems from the fact that the obtained
model depends on the index i and introduces an unwanted coupling between the
distribution functions. A novel approach is therefore used here, in which the ltered
2
equations are projected on the vectors C . The resulting quantities are dened as
tensor in the physical space, and they replace the q-dimensional vectors in the space
of particle distribution functions. This point of view is more favorable to a physical
interpretation of the obtained laws. For this, a map S is introduced which, from a
given conguration of the lattice kinetic variables f , yields the strain rate tensor,
dened in Eq. (1.6) on page 7, at each space position:
S : f (t) S(t).

(8.17)

With this denition, a weaker formulation of Eq. (8.16) is obtained as follows:

S DT

=S

DT ,

(8.18)

Calling U the velocity associated to F , (8.18) becomes

S[u(r, t)] = S[U (r, t)].

(8.19)

This expression is especially useful in the context of LB, as in this case, the rate of
strain can be expressed to the rst order in terms of local variables on every lattice
site.

Turbulence modeling consists in approximating DT by a term depending only on

the ltered variables (in some models, an additional set of variables is introduced).
A very straightforward model introduces the idea of a subgrid viscosity, in which
the original BGK model is implemented, with a renormalized relaxation parameter
:
Fi (r + vci , t + t) Fi (r, t) = (r, t)(Fi (r, t) fieq (U (r, t), (r, t))).

(8.20)

90

Chapter 8

Although this model seems natural, there is no a priori reason to support it. There
are other models, such as the k model presented in the previous section, that make
use of an additional pressure term. For the sake of completeness, it would further
be useful to view turbulence models under the light of the multiple relaxation time
approaches introduced in Section 3.3.1. It is indeed natural to expect that the
value of the renormalized relaxation parameter is dierent for each of the relaxed
models. For the time being, the discussion is however restricted to the case of a
single relaxation time, and a numerical verication of assumption (8.20) is proposed
in the next section.
It is useful to remember that the strain rate tensor is formally related to the
o-equilibrium part of the tensor . The ndings of Eq. (2.66) on page 25 are
reproduced by the following equation:
S

3 neq
;

neq [f ] =

ci ci (fi fieq [f ]) .

(8.21)

With assumption (8.20), equation 8.19 becomes

neq [f ] = neq [F ],

(8.22)

which gives an estimate for the value of the renormalized relaxation parameter:
neq [f ]

= neq

[F ]

, 0, 1, 2.

(8.23)

The numerator of the right hand side is further expanded:

ci ci (Fi fieq [F ] fieq ) = neq [F ]

neq [f ] =

ci ci fieq ,

(8.24)

where the correction fieq to the equilibrium term has been dened as
fieq

fieq [f ] fieq [F ]
9
9
ti Qi : uu U U = ti Qi : ( u u )
=
2
2
9
= ti Qi : .
2

(8.25)
(8.26)
(8.27)

Using the properties of the D3Q19 grid introduced in Section 2.1 ( i ti Qi = 0

and i ti Qi Qi T = 1 (T + T ) for any second order tensor T ), one writes
9
ti Qi Qi =
i

2
,
9

(8.28)

and thus, using the denition of Qi ,

ci ci fieq =
i

1
Qi fieq +
3

fieq =

(8.29)

This yields
neq [f ] = neq [F ] + ,

(8.30)

A model for uid turbulence based on kinetic variables

91

Figure 8.1: 2D cut through the vorticity eld of the homogeneous isotropic turbulent
ow.

and, using the denition of the tensor :

= 1 + neq
= 1 3
,

[F ]
S

(8.31)

In Eq. 8.31, the eective viscosity is expressed in terms of a ratio between the
Reynolds and the strain rate tensor. It is assumed that the various components
of the stress tensor share the same eective relaxation time. This is only true for
homogeneous isotropic ows, but is nonetheless consistent with the point of view
taken by turbulent viscosity models (e.g. Smagorinsky and k ).

8.3

Direct numerical simulation

In this section, Eq. (8.31) is veried numerically with the data of a DNS. The simulation represents a statistically homogeneous and isotropic ow. The system consists
of a cubic lattice of system size 128 128 128 with periodic boundary conditions.
The Reynolds number, dened with respect to the RMS value of the velocity and
the side length of the full cubic system, is Re = 1513. From this, the Kolmogorov
length is found to be = 0.5 in lattice units. Although the Reynolds number is quite
low, it is sucient for obtaining a weakly turbulent ow. The type of forcing needs
however to be chosen very carefully. The method used here is known under the name
of spectral forcing and has been introduced in Ref. [61]. In this method, the force
eld is constructed in the wavenumber space. The force is constructed in such a way
that it is non-zero only at the two lowest wavenumbers, in the low frequency limit.
In this way, the energy is introduced at large wavelengths, and the Kolmogorov
energy cascade can freely develop on the intermediate wavelengths. To obtain a
statistically homogeneous and isotropic eld, all components of the force eld in the
wavenumber space are subject to a random shift of the phase. A pseudo-random
number generator of good quality is used to generate fully decorrelated values for
those shifts. To end with, the force eld is transformed into real space by a Fourier

92

Chapter 8

transform. The method described in Ref. [61] guarantees that the obtained force
eld is divergenceless. In such a way, the simulation stays close to the limit of uid
incompressibility at each iteration step.
Figure 8.1 shows a two-dimensional cut through the uid vorticity during a simulation. This picture is typical for the data observed in turbulent ows. In particular,
eddies are observed at all length scales, from the size of the system down to the size
of a grid cell. The average kinetic energy contained at dierent wavenumbers is
displayed in Fig. 8.2. The curve is irregular at low wavenumbers, where the energy
is injected by the method described in the previous paragraph. After that, the energy scales down along the universal range, and it is nally dissipated in the low
wavelengths, close to the Kolmogorov length. As it is seen on the gure, the intermediate range agrees only partially with the prediction of the Kolmogorov theory.
This is however usual for simulations at such low Reynolds numbers, in which the
turbulence is not fully developed.
The model proposed in Eq. 8.31 is nally cross-checked with the data of the
DNS. The lter introduced in (8.14) is interpreted as a local spatial average over
a cube of size h3 . On every lattice site, is computed from (8.14). Figure 8.3
plots the measured average value for / = / on o-diagonal parts of the strain
rate over all lattice nodes at a given time step. The variance of the value, which
is also plotted on the graph, is unfortunately larger than the average. Therefore,
no denite conclusions can be taken from the data. This is most probably due to
a lack of statistics, and better results can be expected to be obtained at higher
Reynolds numbers. Nevertheless, the graph shows a clear increasing trend of the
average value, and shows that the theory developed in the previous section nds at
least a partial verication in the numerical data.

A model for uid turbulence based on kinetic variables

93

Figure 8.2: Energy spectrum of the turbulent ow.

Figure 8.3: Spatial distribution of the eective relaxation time / based on nondiagonal components of and neq . The marker () labels the value of the eective
relaxation time, averaged over space at a given time step. The standard deviation
over this average is also indicated. In spite of the important standard deviation, a
trend towards an increasing eective relaxation time is visible, which increases with
the depth h of the average. Although the statistics is rather poor, as a low Reynolds
number is used, turbulence modeling via an eective viscosity seems to be supported
by numerical evidence.

Publications
Software projects
An important part of the work achieved during the time of the thesis does not appear
in the present document. Most notably, important eorts were invested in developing
concepts for the numerical implementation of lattice Boltzmann models. Software
paradigms for the object-oriented implementation of advanced data structures, as
well as programming concepts for the development of massively parallel programs
are at the heart of this research topic. During the years of the thesis, three software
projects were formulated, each of which is published on an online web page. On
those web pages the source code of the project is available for free download, under
an open source software license. All projects are furthermore accompanied by a solid
documentation for software users and developers. The OpenLB project is particularly
pointed out, as at the time of the presentation of the thesis, this project is actively
maintainted and developed by several authors. The three projects are listed in the
following:

Vladymir
Vladymir [62] is a library for the C++ programming language that denes a new,
multi-functional array type. The library is specically designed for the purpose
of array-based programming - a programming style in which the use of loops over
array indices is replaced by simple expressions involving the array globally. Based
on this seemingly restrictive programming paradigm, it oers the user two powerful
features: straightforward programming style and automatic parallelization of the
code on both shared and distributed memory environments. Although it certainly
proves useful in many other contexts, the library has mainly been used and tested on
the implementation of scientic models such as Cellular Automata and LB schemes,
thus implementing roughly nearest-neighbor dynamics.

MPTL
MPTL [63] stands for MultiProcessing Template Library, a name that hints both
at the OpenMP standard and the STL. The MPTL proposes a programming model for
writing thread-level parallel programs in C++. Within this model, the programmer
can parallelize algorithms like those of the Standard Template Library (STL) in such
a way that they are executed concurrently by dierent threads. The parallelization
is transparent to the programmer, which means that he doesnt have to tackle any
explicit thread-related code, nor even, in most cases, think about a strategy for how
to parallelize a given algorithm. The syntax of a parallel code remains very similar
to the one used in typical C++/STL programs.

96

Publications

Figure 8.4: Flow behind a 3D backward facing step simulated by the OpenLB software.

OpenLB
The OpenLB project [13] provides a C++ package for the implementation of lattice
Boltzmann simulations that is general enough to address a vast range of problems
in computational uid dynamics. The package is mainly intended as a programming
support for researchers and engineers who simulate uid ows by means of a lattice
Boltzmann method. The source code is publicly available and constructed in a
well readable, modular way. This enables for a fast implementation of both simple
applications and advanced CFD problems. It is also easily extensible to include new
physical content. Pre- and postprocessing tools are also available to congure the
domain geometry and analyze the numerical results. Figure 8.4 shows the output
produced by an example code which is shipped with the standard OpenLB package.

Published Articles
Several topics of this thesis have been published in scientic journals, as it can be
seen from the following list:
[1] J. Latt, S. Succi B. Chopard, and F. Toschi. Numerical analysis of the averaged
ow eld in a turbulent lattice Boltzmann simulation. Physica A, 362:610,
2006.
[2] J. Latt and B. Chopard. Lattice Boltzmann method with regularized nonequilibrium distribution functions. Math. Comp. Sim., 72:165168, 2006.
[3] J. Latt, Y. Grillet, B. Chopard, and P. Wittwer. Simulating an exterior
domain for drag force computations in the lattice Boltzmann method. Math.
Comp. Sim., 72:169172, 2006.
[4] S. Orszag, H. Chen, S. Succi, J. Latt, and B. Chopard. Turbulence eects on
kinetic equations. Journal of Scientic Computing, 28:459466, 2006.

Publications

97

[5] J. Latt, G. Courbebaisse, B. Chopard, and J.-L. Falcone. Lattice Boltzmann

modeling of injection moulding process. In Cellular Automata: 6th International Conference on Cellular Automata for Research and Industry, page 345,
Amsterdam, The Netherlands, October 2004. ACRI 2004.
[6] J. Latt and B. Chopard. Vladymir a c++ matrix library for data-parallel
applications. Future Generation Computer Systems, 20:10231039, 2004.
[7] J. Latt and B. Chopard. An implicitly parallel object-oriented matrix library
and its application to medical physics. IPDPS, 00:254a, 2003.
[8] S. Marconi, B. Chopard, and J. Latt. Reducing the compressibility of a lattice
Boltzmann uid using a repulsive force. Int. J. Mod. Phys. C, 14:10151026,
2003.
[9] L3 Collaboration. The e+ e z q q reaction at lep and constraints

on anomalous quartic gauge boson couplings. Phys. Lett. B, 540:4351, 2002.

Articles to appear
Some submitted articles have been accepted for publication, but have presently not
yet appeared:
[10] J. Latt and B. Chopard. A benchmark case for lattice Boltzmann: turbulent
dipole-wall collision. To appear in Int. J. Mod. Phys. C.
[11] L. Abrahamyan, J. Latt, A. G. Hoekstra, B. Chopard, and P. M. A. Sloot.
Simulating time harmonic ows with the regularized L-BGK method. To
appear in Int. J. Mod. Phys. C.
[12] V. Heuveline and J. Latt. The OpenLB project: an open source and object
oriented implementation To appear in Int. J. Mod. Phys. C.

Bibliography
[1] D. Hnel. Molekulare Gasdynamik. Springer, 2004.
a
[2] A. Quarteroni and A. Valli. Numerical approximation of partial dierential
equations. Springer, 1997.
[3] [Link]
[4] M. Griebel, T. Dornseifer, and T. Neunhoer. Numerical simulation in uid
dynamics. SIAM, 1998.
[5] T. J. Chung. Computational uid dynamics. Cambridge University Press, 2002.
[6] O. C. Zienkiewicz, R. L. Taylor, and P. Nithiarasu. Finite element method for
uid dynamics. Elsevier, 2005.
[7] S. Chen and G.D. Doolen. Lattice Boltzmann method for uid ows. Annu.
Rev. Fluid Mech., 30:329364, 1998.
[8] B. Chopard and M. Droz. Cellular Automata Modeling of Physical Systems.
Cambridge University Press, 1998.
[9] D. Yu, R. Mei, L.-S. Luo, and W. Shyy. Viscous ow computations with
the method of lattice Boltzmann equation. Progr. Aerosp. Sciences, 39:329
367, 2003. [Link] [Link].
[10] C. Cercignani. Theory and application of the Boltzmann equation. Scottisch
Academic Press Ltd., 1974.
[11] S. Succi. The lattice Boltzmann equation, for uid dynamics and beyond. Oxford
University Press, 2001.
[12] M. C. Sukop and D. T. Thorne. lattice Boltzmann modeling; an introduction
for geoscientists and engineers. Springer, 2006.
[13] [Link]
[14] D. A. Wolf-Gladrow. Lattice-Gas Cellular Automata and Lattice Boltzmann
models: an introduction. Lecture Notes in Mathematics, 1725. Springer, 2000.
[15] D. dHumi`res, M. Bouzidi, and P. Lallemand. Thirteen-velocity threee
dimensional lattice Boltzmann model. Phys. Rev. E, 63:066702, 2001.
[16] X. He and L.-S. Luo. Theory of the lattice Boltzmann method: From the
Boltzmann equation to the lattice Boltzmann equation. Phys. Rev. E, 56:6811
6817, 1997.

100

Bibliography

[17] Z. Guo, B. Shi, and C. Zheng. A coupled lattice BGK model for the Boussinesq
equations. Internat. J. Numer. Methods Fluids, 39:325342, 2002.
[18] Z. Guo, C. Zheng, and B. Shi. Discrete lattice eects on the forcing term in
the lattice Boltzmann method. Phys. Rev. E, 65:046308, 2002.
[19] Q. S. Zou, S. L. Hou, S. Y. Chen, and G. D. Doolen. An improved incompressible
lattice Boltzmann model for time-independent ows. J. Stat. Phys., 81:3548,
1995.
[20] X. He and L.-S. Luo. Lattice Boltzmann model for the incompressible NavierStokes equation. J. Stat. Phys., 88:927944, 1997.
[21] P. J. Dellar. Incompressible limits of lattice Boltzmann equations using multiple
relaxation times. J. Comput. Phys., 190:351370, 2003.
[22] P. J. Dellar. Bulk and shear viscosities in lattice Boltzmann equations. Phys.
Rev. E, 64:031203, 2001.
[23] D. dHumi`res. Generalized lattice-Boltzmann equations. Prog. Astronaut.
e
Aeronaut., 159:450458, 1992.
[24] P. Lallemand and L.-S. Luo. Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, galilean invariance, and stability. Phys. Rev. E,
61:65466562, 2000.
[25] R. Benzi, S. Succi, and M. Vergassola. The lattice Boltzmann equation: Theory
and applications. Phys. Reports, 222:145197, 1992.
[26] Y. Chen, H. Ohashi, and M. Akiyama. Prandtl number of lattice bhatnagargross-krook uid. Phys. Fluids, 7:22802282, 1995.
[27] H. Chen, C. Teixeira, and K. Molvig. Digital physics approach to computational
uid dynamics: some basic theoretical features. Int. J. Mod. Phys. C, 8:675
684, 1997.
[28] S. Ansumali and I. V. Karlin. Stabilization of the lattice Boltzmann method
by the H theorem: a numerical test. Phys. Rev. E, 62:79998003, 2000.
[29] S. S. Chikatamarla, S. Ansumali, and I. V. Karlin. Entropic lattice Boltzmann
models for hydrodynamics in three dimensions. Phys. Rev. Lett., 97:010201,
2006.
[30] B. M. Boghosian, P. J. Love, P. V. Coveney, I. V. Karlin, S. Succi, and J. Yepez.
Galilean-invariant lattice-Boltzmann models with H theorem. Phys. Rev. E,
68:025103, 2003.
[31] A. Dupuis and B. Chopard. Theory and applications of an alternative lattice
Boltzmann grid renement algorithm. Phys. Rev. E, page 066707, 2003.
[32] A. J. C. Ladd and R. Verberg. Lattice-Boltzmann simulations of particle-uid
suspensions. J. Stat. Phys., 104:11911251, 2001.

Bibliography

101

[33] H. Chen, R. Zhang, I. Staroselsky, and M. Jhon. Recovery of full rotational

invariance in lattice Boltzmann formulations for high Knudsen number ows.
Physica A, 362:125131, 2006.
[34] L. Kovasznay. Laminar ow behind a two-dimensional grid. Proc. Cambridge
Philos. Soc., 44:5862, 1948.
[35] T. Inamuro, M. Yoshino, and F. Ogino. A non-slip boundary condition for
lattice Boltzmann simulations. Phys. Fluids, 7:29282930, 1995.
[36] P. A. Skordos. Initial and boundary conditions for the lattice Boltzmann
method. Phys. Rev. E, 48:48234942, 1993.
[37] U. Ghia, N. Ghia, and C. T. Shin. High-Re solutions for incompressible ow
using the Navier-Stokes equations and a multigrid method. J. Comp. Phys.,
48:387411, 1982.
[38] Q. Zou and X. He. On pressure and velocity boundary conditions for the lattice
Boltzmann BGK model. Phys. Fluids, 9:15911598, 1997.
[39] I. Halliday, L. Hammond, and C. Care. Enhanced closure scheme for lattice
Boltzmann equation hydrodynamics. J. Phys. A: Math. Gen., 35:L157L166,
2002.
[40] H.J.H. Clercx and C.-H. Bruneau. The normal and oblique collision of a dipole
with a no-slip boundary. Comp. uids, 35:245279, 2006.
[41] R. Mei, L.-S. Luo, P. Lallemand, and D. DHumi`res. Consistent initial condie
tions for lattice Boltzmann simulations. Comp. uids, 35:855862, 2006.
[42] S. Bnisch, V. Heuveline, and P. Wittwer. Adaptive boundary conditions for
o
exterior ow problems. J. Math. Fluid Mech., 7:87107, 2005.
[43] S. Bnisch, V. Heuveline, and P. Wittwer. Leading order down-stream asympo
totics of stationary navier-stokes ows in three dimensions. J. Math. Fluid
Mech., 7:147186, 2005.
[44] P. Wittwer. Second order adaptive boundary conditions for exterior ow problems: non-symmetric stationary ows in two dimensions. J. Math. Fluid Mech.,
8:126, 2006.
[45] X. He, Q. Zou, L. S. Luo, and M. Dembo. Some progress in the lattice Boltzmann method. part I, non-uniform mesh grids. J. Comput. Phys., 129:357363,
1996.
[46] O. Filippova and D. Hnel. Grid renement for lattice-bgk models. J. Comput.
a
Phys., 147:219228, 1998.
[47] J. Tlke, M. Krafczyk, M. Shulz, and E. Rank. Implicit discretization and
o
nonuniform mesh renement approaches for FD discretization of LBGK models.
Int. J. Mod. Phys. C, 9:11431157, 1998.

102

Bibliography

[48] D. Kandhai, W. Soll, S. Chen, A. Hoekstra, and P. Sloot. Finite-dierence

lattice-BGK methods on nested grids. Comput. Phys. Commun., 129:100109,
2000.
[49] M. Bouzidi, D. dHumi`res, P. Lallemand, and L. S. Luo. Lattice Boltzmann
e
equation on a two-dimensional rectangular grid. J. Comput. Phys., 172:704
717, 2001.
[50] A. M. Artoli, A. G. Hoekstra, and P. M. A. Sloot. Accelerated lattice bgk
method for unsteady simulations through mach number annealing. Internat. J.
Modern Phys. C, 14:835845, 2003.
[51] A. Kolmogorov. The equations of turbulent motion in an incompressible uid.
Izv. Sci. USSR Phys., 6:5658, 1942.
[52] J. Smagorinsky. General circulation model of the athmosphere. Mon. Wheather
Rev., 91:99164, 1963.
[53] Launder and Spalding. Mathematical models of turbulence. Academic Press,
New York, 1975.
[54] B. Mohammadi and O. Pironneau. Analysis of the k-epsilon turbulence model.
Wiley, 1993.
[55] C. Teixeira. Incorporating turbulence models into the lattice Boltzmann
method. Int. J. Mod. Phys. C, 9:11591175, 1998.
[56] S. Pope. Turbulent ows. Cambridge University Press, 2000.
[57] S. Succi, O. Filippova, H. Chen, and S. Orszag. Towards a renormalized lattice
Boltzmann equation for uid turbulence. Journ. Stat. Phys., 107:261278, 2002.
[58] S. Ansumali, I. Karlin, and S. Succi. Kinetic theory of turbulence modeling:
smallness parameter, scaling and derivation of Smagorinsky model. Physica A,
338:379394, 2004.
[59] S. Chen, S. Hou, J. Sterling, and G. D. Doolen. A lattice subgrid model for
high reynolds number ows. Fields Inst. Comm., 6:151166, 1996.
[60] H. Chen, S. Succi, and S. Orszag. Analysis of subgrid scale turbulence using
the Boltzmann BGK kinetic equation. Phys. Rev. E, 59:25272530, 1999.
[61] K. Alvelius. Random forcing of three-dimensional homogeneous turbulence.
Phys. Fluids, 11:18801889, 1999.
[62] [Link]
[63] [Link]