Thesis PDF

Thesis
Hydrodynamic limit of lattice Boltzmann equations
LATT, Jonas
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. A novel LB model is introduced, first for the simulation of advection-diffusion
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 simulaiton. The model is
furthermore found to be innovative due to the fact that it possesses an accurate
dimensionless formulation in terms of macroscopic variables. Based on this observation, the
regularized model is investigated to address various challenges of LB, such as the
implementaiton of boundary conditions and spatially refined grids.
Reference
LATT, Jonas. Hydrodynamic limit of lattice Boltzmann equations. Thèse de doctorat :
Univ. Genève, 2007, no. Sc. 3846
URN : urn:nbn:ch:unige-4641
DOI : 10.13097/archive-ouverte/unige:464
Available at:
http://archive-ouverte.unige.ch/unige:464
Disclaimer: layout of this document may differ from the published version.
Université de Genève Faculté des sciences
Département d’informatique Professeur Bastien Chopard
Département de physique Professeur Michel Droz
Hydrodynamic Limit of Lattice

Boltzmann Equations
THÈSE
présentée à la Faculté des sciences de l’Université de Genève
pour obtenir le grade de Docteur ès sciences, mention interdisciplinaire
par
Jonas Lätt
de Mühledorf (SO)
Thèse No 3846
Genève
Atelier de reproduction de la Section de physique
2007
FACULTÉ DES SCIENCES
Doctorat ès Sciences

mention interdisciplinaire
Thèse de Monsieur Jonas LÄTT
intitulée:
Hydrodynamic Limit of Lattice Boltzmann Equations

La faculté des sciences, sur le préavis de Messieurs B. CHOPARD, professeur adjoint
et codirecteur de thèse (Département d’informatique), M. DROZ, professeur titulaire
et codirecteur de thèse (Département de physique théorique), S. SUCCI, professeur
(Istituto per le Applicazioni del Calcolo “Mauro Picone”, Rome, Italie), autorise
l’impression de la présente thèse, sans exprimer d’opinion sur les propositions qui y
sont émises.
Genève, le 23 avril 2007
Thèse -3846-
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 modèle sur la dynamique des fluides émergeante d’une microscopie riche
et abondante, cette thèse est issue d’un tissu d’interactions complexe avec mon
entourage qui m’a soutenu à un niveau scientfique, technique et émotionnel. Je
suis reconnaissant à tous les gens qui m’ont suivi durant ma thèse et espère qu’ils
trouveront plaisir, en lisant ce document, à y reconnaitre une contribution, une idée
ou un biais qui leur est dû.
Je pense en particulier à mes directeurs de thèse Bastien Chopard et Michel Droz
qui m’ont suivi au long du travail, et à Sauro Succi qui a été conseiller et juré de
thèse. Des conversations stimulantes ont permis de situer plusieurs sujets dans un
contexte plus large, et j’aimerais faire valoir entre autres les contributions de Peter
Wittwer, de Vincent Heuveline et de Michel Deville.
J’ai profité d’un environnement de travail inspirant grâce à Jean-Luc et Bernhard
qui m’ont aidé à affronter le monde hostile de l’informatique, Fokko qui a soigneuse-
ment relu ma thèse, Orestis qui a prêté une oreille critique à mes idées farfelues et
mes autres collègues qui m’ont tous donné un soutien remarquable.
La contribution la plus précieuse vient sans aucun doute de mon entourage
proche, et je remercie Roxane, ma mère, les parents de Roxane et toute ma famille
pour leur patience et leur appui.
Contents
Abstract v
Résumé en Français vii
Introduction 1
1 Computational fluid dynamics and lattice Boltzmann 5

1.1 Vector and tensor notation . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Macroscopic equations . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . 6
1.2.2 Isothermal and incompressible flows . . . . . . . . . . . . . . . 7
1.2.3 Dimensionless formulation . . . . . . . . . . . . . . . . . . . . 8
1.2.4 Further reading . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Lattice Boltzmann method . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3.1 Discrete variables . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3.2 Lattice Boltzmann model . . . . . . . . . . . . . . . . . . . . 10
1.3.3 Lattice structures . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3.4 Further reading . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2 Multiscale Chapman-Enskog analysis 15

2.1 Series expansion with scale separation . . . . . . . . . . . . . . . . . . 15
2.1.1 Lattice symmetries . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.2 Multi-scale expansion . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.3 Conservation laws . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2 Advection-diffusion: Chapman-Enskog ansatz . . . . . . . . . . . . . 21
2.3 Fluid flows: Chapman-Enskog ansatz . . . . . . . . . . . . . . . . . . 23
2.4 Dimensionless formulation and accuracy . . . . . . . . . . . . . . . . 26
2.4.2 Accuracy of LB methods . . . . . . . . . . . . . . . . . . . . . 28
2.4.3 Accuracy of specific models . . . . . . . . . . . . . . . . . . . 28
3 Corrections to the BGK dynamics 31

3.1 Advection-diffusion revisited . . . . . . . . . . . . . . . . . . . . . . . 31
3.1.1 Second-order time stepping . . . . . . . . . . . . . . . . . . . 31
3.1.2 Numerical verification . . . . . . . . . . . . . . . . . . . . . . 33
3.2 Bulk viscosity for fluid flows . . . . . . . . . . . . . . . . . . . . . . . 34
3.2.1 Numerical error in incompressible flows . . . . . . . . . . . . . 34
3.2.2 Correction term for the bulk viscosity . . . . . . . . . . . . . . 36
3.3 Alternative LB models . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.3.1 LB models with multiple relaxation times . . . . . . . . . . . 39
3.3.2 Entropic LB models . . . . . . . . . . . . . . . . . . . . . . . 41
ii Contents
4 Regularized lattice Boltzmann 43

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.2 Regularization procedure . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.2.1 Regularized particle distribution functions . . . . . . . . . . . 45
4.2.3 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.2.4 Related models . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.3 RLB and MRT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.4 Numerical verification . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5 Boundary and initial conditions 53

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.2 Boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.2.2 Implementation of the pressure . . . . . . . . . . . . . . . . . 55
5.2.3 Regularization of on-wall distribution functions . . . . . . . . 55
5.2.4 Numerical verification . . . . . . . . . . . . . . . . . . . . . . 57
5.3 Initial condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.3.2 The benchmark . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.4 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.4.1 Time evolution of the flow . . . . . . . . . . . . . . . . . . . . 61
5.4.2 Benchmark values . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.5 Setup of the initial condition . . . . . . . . . . . . . . . . . . . . . . . 63
5.6 Discussion: choice of the time step . . . . . . . . . . . . . . . . . . . 64
6 Adaptive space and time steps 65

6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
6.2 Grid refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
6.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
6.2.2 Analytical profile on the boundaries . . . . . . . . . . . . . . . 68
6.2.3 Numerical implementation . . . . . . . . . . . . . . . . . . . . 69
6.2.4 Simulation results . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.3 Adaptive time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
6.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
6.3.2 Numerical experiment . . . . . . . . . . . . . . . . . . . . . . 72
7 Coupling with other tools of CFD 75

7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
7.2 FD model and computational grids . . . . . . . . . . . . . . . . . . . 76
7.3 Choice of units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
7.4 The FD-LB interface . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
7.5 The coupling algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 79
7.5.1 Coupling the pressure . . . . . . . . . . . . . . . . . . . . . . 79
7.5.2 Coupling the velocity . . . . . . . . . . . . . . . . . . . . . . . 79
7.5.3 Coupling the velocity gradients . . . . . . . . . . . . . . . . . 80
7.5.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
Contents iii
7.6 Numerical validation . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

7.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
8 A model for fluid turbulence based on kinetic variables 85

8.1 Turbulence modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
8.1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
8.1.2 Large eddy simulations in 3D incompressible fluids . . . . . . 86
8.2 Averaged LB equation . . . . . . . . . . . . . . . . . . . . . . . . . . 88
8.3 Direct numerical simulation . . . . . . . . . . . . . . . . . . . . . . . 91
Publications 95
Bibliography 99
Abstract
An analysis of the lattice Boltzmann (LB) method is conducted, and various con-
clusions 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 fluids with extended physical
properties. The specific scope of the present thesis is motivated by current prac-
tice in the domain of computational fluid 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 first 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 limi-
tations of the lattice Boltzmann method are discussed in simple terms.
The second part introduces novel LB models, first for the simulation of advection-
diffusion 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 sim-
ulation. 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 refined grids to the coupling of the
LB method with a finite difference 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 fluid. Based on the mi-
croscopic LB variables, a model for the simulation of fluid turbulence is formulated.
This document presents a selective choice of the topics investigated during the
thesis, aiming at a succinct overview of the main findings. It focuses on theoretical
aspects and does not reflect another important part of the thesis, related to imple-
mentation 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.
Résumé en Français
Aperçu de la problématique
Introduction
La présente thèse de doctorat propose une étude de la méthode de Boltzmann sur
réseau, dénotée par l’acronyme LB pour son nom anglais de “lattice Boltzmann”.
Diverses conclusions sont tirées sur les possibilités d’exploiter cette méthode pour
la résolution numérique de l’équation de Navier-Stokes. Cette manière restrictive
d’utiliser la méthode LB doit être vue par opposition à son utilisation tradition-
nelle, se situant dans le contexte de l’étude de fluides complexes avec des propriétés
physiques étendues. Le choix adopté dans la thèse est motivé par l’état de l’art
dans le domaine de la simulation de fluides (CFD pour l’expression anglaise “com-
putational fluid dynamics”). En effet, on observe qu’un nombre croissant de scien-
tifiques et ingénieurs utilisent la méthode de LB comme outil alternatif aux moyens
de résolution numérique conventionnels pour l’équation de Navier-Stokes. Le cadre
conceptuel à disposition actuellement pour une application du LB dans ce domaine
est par contre toujours faible. Il en ensuit un nombre de limitations, auxquelles
certains remèdes sont proposés dans la thèse.
Dans la première partie, la méthode LB est reformulée d’une telle manière à
mettre en évidence les variables de nature macroscopique, au lieu des quantités
microscopiques issues de la théorie cinétique des gaz. De cette nouvelle perspective,
la précision ainsi que les limitations de la méthode de Boltzmann sur réseau sont
étudiés en termes simples.
La deuxième partie introduit de nouveaux modèles LB, d’abord pour la simula-
tion de problèmes à convection-diffusion, et ensuite pour la résolution de l’équation
de Navier-Stokes. Il est démontré que ce dernier modèle, présenté sous le nom de
“lattice Boltzmann régularisé”, augmente de manière substantielle la stabilité et la
précision des modèles LB au cours de simulations. Un autre aspect innovant de
ce modèle se réfère au fait qu’il s’écrit de manière exacte en terme de variables
macroscopiques.
Les implications techniques de cette observation sont étudiées dans la troisième
partie: la méthode LB régularisée est utilisée pour résoudre divers problèmes pos-
sédant des échelles temporelles et spatiales variables. Il s’agit là par exemple de
grilles numériques localement raffinées ou de couplages entre une méthode LB et
une méthode de différences finies. Le point de vue adopté à travers cette étude
s’appuie sur l’idée que pratiquement toute question d’implémentation pratique en
LB se résoud à travers une analyse de Chapman-Enskog correspondante.
Cette optique est délibérément modifiée lors de la dernière partie de la thèse,
dans laquelle la méthode LB originale, non-régularisée est utilisée pour la simulation
d’un fluide complexe. Sur la base de variables LB microscopiques, un modèle pour
la simulation de fluides turbulents est formulé.
viii Résumé en Français
Ce document contient un choix sélectif des sujets étudiés au cours de la thèse,

afin de résumer brièvement les découvertes majeures. Il se concentre sur les aspects
théoriques et omet de détailler une autre partie importante de la thèse, relative aux
questions d’implémentation numérique. On peut en particulier mettre faire valoir
le travail investi dans la conception de trois projets de logiciel indépendants, in-
troduisant de nouveaux paradigmes de programmation pour l’implantation de pro-
grammes parallèles et de modèles LB. En outre, ces projets proposent des codes
source libres et pleinement documentés. Cette partie du travail est résumée à la fin
de la thèse, accompagnée d’une liste des articles publiés au cours de la thèse.
Simulations en dynamique des fluides

La recherche en simulation de dynamique des fluides (CFD) se concentre sur une
équation a différences partielles (PDE pour l’anglais “partial differential equation”)
connue sous le nom d’équation de Navier-Stokes incompressible qui exprime une loi
de conservation locale pour la quantité de mouvement du système. Elle représente
un modèle largement simplifié pour un fluie inviscide à une seule composante dans un
environnment sans variations de densité ou de température, et sans autre ingrédient
physique. Bien que cette équation n’adresse que partiellement la complexité de
la plupart des fluides considérés en ingénierie, elle est utilisée avec beaucoup de
succès dans différents domaines afin de prédire qualitativement le comportement de
systèmes immergés dans un fluide. Comme exemple, on peut nommer les problèmes
d’optimisation en aérodynamique, tel que la conception d’une carrosserie de voiture
ou la création d’ailes d’avion. L’équation de Navier-Stokes incompressible n’arrive
bien évidemment à décrire la nature de l’air compressible que de manière très
approximative. Néanmoins, cette équation mène à des prédiction de très bonne
qualité qu’on utilise fréquemment, conjointement aux résultats de la simulation
expérimentale. On peut par ailleurs noter que la résolution de cette PDE représente
une tâche de calcul intensive, malgré la simplicité apparente de l’équation. On in-
vestit donc beaucoup d’efforts dans le développement de techniques de résolution
numérique, plutôt que de rajouter davantage de complexité en ajoutant des termes
physiques additionnels.
D’un autre côté, les domaines de recherche se concentrant sur des fluides dits
complexes ont également une grande importance dans le domaine. Les fluides
en question possédent une physique enrichie: pour ceux-ci, l’approximation de
l’équation de Navier-Stokes est donc insuffisante. L’approche classique adoptée en
CFD pour traiter ce type de fluides consiste en la description de nouvelles propriétés
physiques en terme de phénomènes de transport reliées à une propriété macro-
scopique. Une nouvelle PDE est formulée pour la dynamique de cette propriété,
qu’on résout ensuite par une technique numérique appropriée. L’interaction entre
les différents phénomènes de transport est traitée à travers de termes ajoutés dans
les PDE’s respectives. Dans un fluide avec d’importantes variations de température
par exemple, la température est introduite comme une nouvelle propriété, et sa
dynamique est décrite par une équation de transport pour la chaleur. En plus,
l’équation de conservation du moment est étendue pour tenir compte de la com-
pressibilité du fluide. La température et la vitesse sont ensuite liés à travers un
terme convectif dans l’équation de la chaleur et un terme de poussée verticale dans
Résumé en Français ix
l’équation du moment. Une équation d’état additionnelle couple la température

avec la densité. On peut dire que le domaine traditionnel de la CFD adopte une
approche de haut niveau à la modélisation d’un fluide, en commençant par une loi
de conservation simple, et en rajoutant de nouvelles équations sur demande, basées
sur des considérations théoriques et phénoménologiques.
Méthode de Boltzmann sur réseau

Le sujet principal de la thèse se traduit par une approche numérique du nom de
Boltzmann sur réseau (LB pour le terme anglais “lattice Boltzmann”). Il s’agit
d’un méthode realtivement récente qui se distingue méthodes traditionelle en adop-
tant une approche de bas niveau à la modélisation de fluides. A cet effet, elle décrit
le fluide à un niveau moléculaire et élabore des modèles pour l’interaction entre
molécules. La physique complète du fluide à un niveau du continu se trouve im-
plicitement décrite par le modèle, et le travail conceptuel qu’on doit effectuer se
trouve à l’opposé de ce qu’on ferait dans l’approche de haut niveau classique. En
effet, les ingrédients physiques du modèle doivent être identifiés un par un et ex-
plicités propement, afin de permettre une séparation entre propriétés relevantes et
négligeables. En se basant sur ces considérations, le modèle de collision est large-
ment simplifié pour les besoins du traitement numérique. En tant que tel, le LB
est une méthode bien adaptée à la simulation de fluides complexes, étant donné
que les ingrédients physiques relevants sont pris en charge très naturellement par la
méthode. Il est mis en évidence à quel point une formulation LB d’un problème de
dynamique des fluides peut être élégante lorsqu’on on simule par exemple des flu-
ides à deux phases. En CFD traditionnelle, une PDE est formulée pour chacune des
deux phases, ainsi que pour la dynamique d’interaction sur l’interface entre celles-
ci. L’implémentation numérique d’un tel modèle requiert l’utilisation de techniques
avancées en génie logiciel, afin de correctement tracer la position de l’interface et
implémenter la dynamique correspondante. Ces difficultés n’apparaissent pas dans
la méthode LB, car les deux phases peuvent être modélisées par un modèle de fluide
unique, et l’interface entre les phases est traitée automatiquement, en ingrédient
naturel du modèle.
La méthode LB appliquée au domaine de la CFD

Lorsque l’approche LB à la modélisation physique est importée dans le domaine de
CFD classique, deux mondes très différents, possédant leurs propres convictions et
traditions sont confrontés. Ceci mène à de nombreux malentendus, d’un côté de
la part de la communauté de la CFD qui manque de reconnaitre le potentiel de la
méthode LB par sa juste valeur, et d’un autre côté de la communauté LB qui éprouve
des difficultés à reconnaitres les enjeux réels de la CFD, ou à comprendre quels sont
les domaines dans lesquels le LB peut être compétitifs avec d’autres approches.
Un obstacle majeur à un échange mutuel entre ces deux méthodes vient du fait
qu’elles décrivent le fluide à travers un choix de variables différent. Les variables
centrales en CFD sont les quantités macroscopiques, observables. Dans l’équation
de Navier-Stokes incompressible, les degrés de liberté sont définis par la vitesse et
la pression, dépendants de l’espace et du temps. Les outils de résolution numérique
x Résumé en Français
en CFD classique modèlisent souvent directement ces variables-ci, c’est-à-dire, ils

représentent la valeur adoptée par ces variables sur les noeuds d’une grille numérique
donnée. La méthode LB d’un autre côté simule la dynamique d’une quantité appelée
fonction de distribution de particules, dérivée de la théorie cinétique des gaz. S’il
est bien connu comment cette quantités peut être convertie en variables d’intérêt
macroscopique, la procédure inverse par contre est plus compliquée. De fait, il
n’est pas toujours clair de quelle manière une configuration donnée des variables
macroscopiques est utilisée pour initialiser la fonction de distribution. Ce problème
est rarement mis en évidence par une communauté habituée à penser en terme
de quantités cinétiques uniquement. Il est néanmoins à la source de complications
substantielles lorsque des variables macroscopiques sont évoquées explicitement dasn
la simulation, ou lorsque les échelles d’espace et de temps macroscopiques du modèles
sont ajustées.
Comment traiter des échelles temporelles et spatiales vari-

ables?
Les problèmes liés à la conversion entre variables macroscopiques et fonctions de

distribution apparaissent par exemple lorsque l’on applique une technique de raffine-
ment de grille, au cour de laquelle le domaine de la simulation est partitionné en sous-
domaines, dont chacun est représenté par une grille numérique possédant un degré
de résolution variable. Le but de cet exercice est d’adapter la structure de la grille à
la géométrie du problème, et par conséquent, d’augmenter l’efficacité du calcul. Lors
de l’exécution d’une simulation, il s’avère nécessaire de communiquer les valeurs cal-
culées au fur et à mesure d’une grille à l’autre. Ce procédé s’applique sans difficulté
dans le contexte des outils de CFD classiques. En effet, ces méthodes gardent en
mémoire la valeur sans dimensions des variables macroscopiques, indépendante donc
de la résolution de la grille ou du pas temporel, s’échangeant entre différente grilles
sans besoin de conversion. La fonction de distribution de particules utilisée en LB
par contre n’est pas adimensionnelle par rapport au choix d’unités macroscopiques.
Elle doit donc être remise à l’échelle lors du transfer d’une grille à une autre, puisque
sa valeur dépend de l’espacement de la grille. La règle qui détermine cette remise à
l’échelle n’est aucunement évidente, et doit être appliquée avec beaucoup de soins.
Par conséquent, une certaine confusion règne à ce sujet dans la litérature, et de
nombreuses approches différentes ont été proposées, dont les idées de base diffèrent
substantiellement.
Le raffinement n’est qu’un exemple particulier, introduit ici dans un souci d’illus-
tration. La figure 1 affiche d’autres problèmes à échelles d’espace et de temps vari-
ables, nécessitant une transformation entre les unités macroscopiques et celles du LB:
l’implantation de conditions aux bords et de conditions initiales, et l’implantation
de modèles hybrides, faisant appel en même temps à une méthode LB et à un outil
classique de la CFD.
Résumé en Français xi
Conditions aux bords Conditions initiales
Couplage avec d’autres méthodes Applications multi-grille
Figure 1: Divers problèmes se posent lors de l’implantation d’un modèle LB, pour
lesquels une conversion entre variables macroscopiques et fonction de distribution
est requise.
Méthode LB et développement multi-échelle
L’équation de Navier-Stokes incompressible

Dans ce qui suit, les problèmes de la dynamique des fluides sont réduits à un point
de vue purement macroscopique. Dans ce contexte, les questions concernant la
simulation d’un fluide équivalent à des questions de résolution numérique d’une
équation différentielle à dérivées partielles. Un fluide incompressible et isotherme
de viscosité ν est par exemple décrit par l’équation suivante, connue sous le nom
d’équation de Navier-Stokes:

∂~u −→ p
+ ∇ · ~u ⊗ ~u + I − 2νS = 0.
∂t ρ0
Elle exprime la conservation de la quantité de mouvement d’un fluide caractérisé par

sa vitesse ~u et sa pression p. La PDE s’écrit sous forme d’une divergence d’un terme
tensoriel. Il est sous-entendu que la divergence est appliquée au deuxième indice
du tenseur. En outre, on désigne par le symbole ⊗ le produit tensoriel entre deux
vecteurs. Le tenseur I désigne l’identité, et S le tenseur des taux de déformation,
défini comme la composante symétrique du tenseur des dérivées premières du champ
de vitesse:
T
1 − → −→
S= ∇ ⊗ ~u + ∇ ⊗ ~u
2
xii Résumé en Français
Figure 2: Le réseau D2Q9 pour la simulation de fluides en deux dimensions.
Méthode de Boltzmann sur réseau

La méthode de Boltzmann sur réseau est présentée dans le chapitre 1.3 de la thèse,
ainsi que dans les références indiquées à la fin de ce chapitre. Un aperçu rapide de
la méthode est reproduit ici, afin d’introduire la notation.
La méthode LB considère une version discrète de la fonction de distribution
f (~x, ~c, t), définie sur l’espace des positions ~x et des vitesses ~c. Lors de la discrétisation,
l’espace des vitesses ~c est omis en faveur d’une augmentation du nombre de fonctions
de distribution, dont un nouvel exemplaire est introduit pour chaque individu d’une
population de vecteurs de vitesse discrète. Chacune de ces fonctions de distribution
est identifiée par un indice k: fk = fk (~x, t). Les vecteurs vitesse utilisés dans le
cas discret représentent la relation spatiale entre une cellule du réseau et ses plus
proches voisins. En deux dimensions, on utilise par exemple fréquemment le réseau
D2Q9 qui contient neuf vecteurs, un pour chacun des huit voisins, et un pour le
vecteur nul, représentant des “particules au repos” ne quittant pas la cellule. Ce
réseau est illustré sur la figure 2.
La dynamique discrète d’un modèle LB exprime la propagation de flux de partic-
ules d’un noeud ~xij vers ses voisins. L’équation est écrite dans un système d’unités
du réseau, dans lequel l’espacement spatial entre deux noeuds voisins de la grille,
ainsi que l’espacement temporel entre deux pas d’itération, sont unitaires. La dy-
namique LB s’exprime donc par la formule
fk (~xij + ~ck , t + 1) − fk (~xij , t) = Ωk (~xij , t),
où le terme Ω désigne l’opérateur de collision entre particules. Dans le modèle le

plus couramment utilisé, connu sous le nom de BGK, la collision s’écrit comme une
relaxation vers un équilibre local:
Ωk = −ω (fk − fkeq ) ,
où on a défini l’équilibre, une version discrète de la distribution de Maxwell-Boltzmann,

ne dépendant que des variables macroscopiques, de la manière suivante:

eq eq 1 1 2 1 2
fk = fk (ρ, ~u) = ρ tk 1 + 2 ~ck · ~u + 4 (~ck · ~u) − 2 |~u| .
cs 2cs 2cs
Le paramètre de relaxation ω peut être formellement relié à la viscosité du fluide.
Résumé en Français xiii
Les quantités macroscopiques peuvent être calculées à partir de fonctions de

distribution en calculent les moments d’ordre zéro, et du premier ordre. On trouve
ainsi
P
• La densité: ρ = k fk .
P
• La vitesse: ~u = 1/ρ k ~ck fk .
Bien qu’en principe la méthode LB décrive un fluide compressible, on retrouve
l’équation d’un fluide incompressible dans la limite des faibles nombres de Mach.
Dans ce cas, la relation entre la densité et la pression est décrite par l’équation pour
un gaz parfait:
p = c2s ρ.
Développement multi-échelle
Bien qu’on ait montré dans la section précédente comment calculer les variables
macroscopiques à partir des fonctions de distribution, la procédure inverse n’a
jusqu’alors pas été abordée. Afin de connecter en détail la dynamique LB avec
une dynamique macroscopique, la limite du continu doit d’abord être évaluée par
un développement de Taylor du deuxième ordre, et une technique de séparation des
échelles temporelles et spatiales est utilisée pour évaluer la limite hydrodynamique
du modèle. Le développement de Taylor se résume par l’approximation suivante:
2
→
− 1 2
−
→
fk (~xij + ~ck , t + 1) − fk (~xij , t) ≈ (∂t + ∇ · ~ck ) fk + ∂t + 2∂t ∇ · ~ck + ∇ · ~ck fk .
2
Dans le processus de séparation des échelles, les fonctions fk , ainsi que les dérivées
spatiales et temporelles, sont développés en fonction d’un paramètre formel petit,
c’est-à-dire 1:
(0) (1)
fk = fk + fk + O(2 ),
∂t = ∂t1 + 2 ∂t2 + O(3 ) et
→
− →
−
∇ = ∇ (1) + O(2 ).
En réintroduisant ces approximations dans le modèle BGK, on retrouve bien
l’équation de Navier-Stokes, et on conclut que ce modèle simule proprement la dy-
namique d’un fluide. Ce résultat est bien connu est représente un des piliers sur
lesquels repose la méthode. Ce n’est par contre pas ce résultat en lui-même qui
nous intéresse présentement, mais plutôt un résultat intermédiaire, reliant les fonc-
tions de distribution aux champs macroscopiques.
Le terme constant en est en effet équivalent à la fonction d’équilibre, et ne
dépend donc que de ρ et de ~u:
(0)
fk = fkeq (ρ, ~u).
Le terme d’ordre 1 est lié aux dérivées de la vitesse, ainsi que du carré de la
vitesse:
(1) tk →
−
fk = − 2 Qk : ρ ∇ ⊗ ~u −
cs ω
→
− 1 →
−
~ck ⊗ ∇ : ρ~u ⊗ ~u + 2 (~ck · ∇)(Qk : ρ~u ⊗ ~u) ,
2cs
xiv Résumé en Français
où on a défini le tenseur symétrique

Qk = ~ck ⊗ ~ck − c2s I.
Il est important de remarquer que le détail des fonctions de distribution n’est
pas relevant dans le calcul de la limite hydrodynamique du modèle. L’importance
revient plutôt au moment du deuxième ordre, Π(1) , calculé selon la formule suivante:
X (1)
Π(1) = ~ck ⊗ ~ck fk .
k
Au cours du développement multi-échelle, on montre que ce tenseur est proportionnel

au tenseur des taux de déformation:
2c2s
Π(1) = − S.
ρ
Lors du calcul du deuxième moment des fonctions de distribution, tous les termes
contenus dans la fonction de distribution n’entrent pas en compte. En effet, les
termes définis au carré de la vitesse s’éliminent pour des raisons de symétrie. Dans
le but d’obtenir le “bon” tenseur Π(1) produisant une hydrodynamique appropriée,
il est donc suffisant de considérer une version réduite des fonctions fk :
(1) tk →
−
f¯k = − 2 Qk : ρ ∇ ⊗ ~u .
cs ω
Cette expression relie de manière simple le champs macroscopique aux valeurs des
fonctions de distribution.
Modèle de Boltzmann sur réseau régularisé

Procédé de régularisation
Dans le chapitre précédent, on a introduit une version simplifiée des fonctions de dis-
tribution fk . Elle se désigne par le terme f¯k et représente proprement les composants
hydrodynamiques du modèle. Ces variables ont été reliées par une relation simple
aux dérivées du champ de vitesse. Pour l’utilisation courante de cette expression, il
serait par contre plus utile que les fonctions de distribution soient directement liées
aux fonctions macroscopiques simulées, donc la densité ρ, la vitesse ~u et le tenseur
des taux de déformation S. En effet, il a été montré que ces trois grandeurs peuvent
être évaluées localement en chaque noeud du réseau LB, en calculant les moments
des fonctions de distribution.
Une telle relation se trouve en exploitant la symétrie du tenseur Q:
(1) tk ρ →
−
f¯k = − 2 Qk : ∇ ⊗ ~u
cs ω
tk ρ 1h − → →
− T
i
= − 2 Qk : ( ∇ ⊗ ~u) + ( ∇ ⊗ ~u)
cs ω 2
tk ρ
= − 2 Qk : S
cs ω
tk
= Q : Π(1)
2c4s ω k
Résumé en Français xv
Cette relation est centrale dans pratiquement tous les sujets développés dans la
thèse. Elle relie les champs macroscopiques aux valeurs des fonctions de distribution,
et peut donc être exploitée systématiquement, à chaque fois qu’une telle transfor-
mation est nécessitée. Pour souligner cette importance, la formule de régularisation
est reproduite, en omettant cette fois-ci les calculs qui ont permis de la trouver:
(1) tk
f¯k = 4 Qk : Π(1)
2cs ω
En suivant l’enchainement suivant, la fonction de régularisation est utilisée pour

évaluer les valeurs des fonctions de distribution à partir des variables macroscopiques:
ρ
Π(1) (1)
~u fkrégularisé = fkeq + f¯k
S fkeq
Modèle LB régularisé
La formule de régularisation est utilisée dans plusieurs domaines qui sont présentés
dans la thèse. A titre d’exemple, l’application probablement la plus importante est
expliquée dans ce qui suit: le modèle de Boltzmann sur réseau modifié, dont les
fonctions de distribution sont régularisées.
L’algorithme pour la collision entre particules dans ce modèle se décompose en
trois étapes:
1. Calcul de ρ, ~u et Πneq = k ~ck ⊗ ~ck (fk − fkeq ).

P
2. Régularisation fk ← fkrégularisé
3. Exécution de la collision BGK sur les fonctions de distribution régularisées.
Afin d’apprécier les qualités de ce nouveau modèle, un flux bidimensionnel connu

sous le nom de flux de Kovasznay a été implanté numériquement, une fois en utilisant
la méthode BGK, et une fois en s’appuyant sur le modèle régularisé. Une solution
analytique pour ce flux est connue, et on peut donc vérifier la précision de la solution
numérique, en la comparant à l’expression analytique. Le résultat de ce test a montré
que le modèle régularisé est nettement plus précis que le modèle BGK, et que le gain
en précision augmente au fur et à mesure que la grille est raffinée. Sur des grilles
de taille commune, on peut gagner un ordre de précision en utilisant le modèle
régularisé.
Un autre modèle bidimensionnel a été implanté, désigné par le nom de flux de
cavité. Il représente un fluide contenu dans une boite quadratique, dont la paroi
supérieure est entrainée à une vitesse constante. Cette fois-ci, la stabilité numérique
du modèle a été testée, c’est-à-dire, la valeur maximale du nombre de Reynolds
pouvant être atteinte avant que la simulation soit numériquement instable. Dans
les deux cas, BGK et régularisé, une relation linéaire entre le nombre de Reynolds
maximal et le paramètre de résolution N d’une grille N ×N est établie. La croissance
xvi Résumé en Français
linéaire est par contre sept fois plus rapide dans le cas du modèle BGK. Dans les
détails, on trouve
Pour le modèle BGK: Remax = 0.4 + 1.0 ∗ N
Pour le modèle régularisé: Remax = 16 + 7.4 ∗ N
Il est donc établi que pour ce type de problèmes, le modèle régularisé est net-
tement plus précis et plus stable que le modèle BGK. Cet avantage s’explique par
l’origine de la méthode régularisée. Elle a été construite en imposant une séparation
des échelles au modèle original, le modèle BGK. Uniquement les termes hydro-
dynamiques, c’est-à-dire, le termes contribuant à la dynamique de l’équation de
Navier-Stokes, ont été préservés, et les termes d’ordre supérieur sont négligés. Il
est vrai qu’à cause de cette simplification, la méthode régularisée s’est apauvrie.
On s’attend en particulier à ce qu’elle soit moins apte à présenter spontanément des
ingrédients physiques d’un fluide complexe. Dans un but particulier par contre, celui
de la résolution de l’équation de Navier-Stokes, elle est plus concise, plus pratique
à manier, plus précise et plus stable.
Introduction
Computational fluid dynamics
The present thesis investigates a methodology in the domain of computational fluid

dynamics (CFD), that is, the numerical simulation of fluid flows. The core of this
research area is represented by a partial differential 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 simplified model for
a viscid, single component fluid in an environment without density or temperature
variations nor any other specific physical ingredients. Although this equation only
partially addresses the complexity of most fluids of interest in engineering applica-
tions, it is successfully applied in different areas for qualitative predictions of fluid
flows. 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 incompress-
ible Navier-Stokes equation only roughly approximates the nature of compressible
air. Nevertheless, this equation leads to predictions of high quality that are suc-
cessfully 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 effort 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 flu-
ids that exhibit extended physical properties, and to which the basic Navier-Stokes
equation is an insufficient approximation. The classical approach in CFD to the
treatment of such fluids is to describe the new physical properties in terms of trans-
port 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 appropri-
ate numerical technique. The interaction between the different transport phenom-
ena is handled via interaction terms between the respective PDE’s. In a fluid with
important temperature variations for example, a new observable property, the tem-
perature, is introduced and its dynamics is described by a heat transport equation.
Furthermore, the momentum equation is extended to comprehend compressible flu-
ids. 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
field of CFD as a top-down approach to fluid modeling, starting with a simple con-
servation law, and adding new equations on demand, based on both theoretical and
phenomenological considerations.
2 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 fluid modeling. To achieve
this, it describes the fluid at a molecular level and proposes models for the collision
between molecules. The full continuum-level physics of the fluid is implicitly con-
tained in this model, and the conceptual work to be effectuated is opposite to the
one in the top-down approach of classical CFD. Indeed, the various physical ingredi-
ents contained in the model need to be identified 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 simplified for the needs
of the numerical treatment. As such, the lattice Boltzmann (LB) method is well
adapted to the simulation of complex fluids, as the relevant physical ingredients are
taken in charge very naturally by the method. How elegant a LB formulation of
a fluid problem can be becomes apparent for example in the context of two-phase
flows. 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 implemen-
tation of such a model requires advanced software engineering techniques to track
the position of this interface and to implement its dynamics. Those problems don’t
appear in the LB method, where the two phases can be represented by the same
fluid 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 fluids, the recent
historical development shows that much effort is invested in solving the raw Navier-
Stokes equation with this method. This is reflected by numerous publications and
technical reports, in which the LB method is used to solve benchmark applications
of classical CFD, and its relative efficiency 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 different 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 different worlds with different traditions and convictions are con-
fronted. 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 difficulties 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 fluid via a different choice of vari-
Introduction 3
ables. The variables of concern in CFD are the macroscopic, observable quantities.
In the incompressible Navier-Stokes equation, the degrees of freedom are defined
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 configuration 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 exclu-
sively 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 difficulty becomes apparent when a grid refinement 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 efficiency. During the simulation, it is necessary to
communicate successively computed values from one grid to another. This is easily
performed in classical fluid 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 specific 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 different approaches have been proposed,
based on substantially different points of view. Grid refinement is just a case in
point to illustrate the encountered difficulties, but similar problems are encountered
with essentially all practical implementation issues of LB, such as the definition 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 infinite series via the so-called Chapman-Enskog multi-scale
expansion. This expansion is commonly used to verify that a given LB method
effectively solves the macroscopic PDE’s asymptotically. To do this, the infinite
4 Introduction
series is truncated. The remaining lower order terms, sufficient 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
different problems in a common framework.
Content of the thesis

In the first part of the thesis, various theoretical results from the literature are
gathered to show that a LB implementation effectively 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 fluid 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-diffusion problems. On the way, it is shown
that this advection-diffusion model is treated erroneously in the literature, and a
workaround is suggested.
In the second part, a modified LB model is introduced under the name of regular-
ized 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 specificity of the regularized distribution
function is that it can be fully expressed in terms of macroscopic, observable vari-
ables. 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 accu-
rate 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 uni-
fied approach, it is explained how to implement grid refinement, an adaptive time
interval, boundary conditions, initial conditions, and a coupling of the LB model
with another CFD tool. In the last chapter finally, 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 fluid turbulence. The theoretical considerations are accompanied by results of
direct numerical simulations of a turbulent fluid.
Chapter 1
Computational fluid 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 λ = ~a · ~b λ = 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, some-
times 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. Further-
more, 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
6 Chapter 1
expressions.
1.2 Macroscopic equations

1.2.1 Governing equations
Computational Fluid Dynamics (CFD) provides a qualitative, and in some cases
even quantitative, prediction of fluid flows 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 differential equa-
tions (PDE’s). They describe the evolution of the observable, physical quan-
tities in the flow and are based on phenomenological and/or theoretical con-
siderations. As an alternative, a simplified microscopic model can be directly
simulated on a computer, without referring explicitly to a PDE.
The Numerical method provides a way to find approximate solutions of the

PDE, or to solve the dynamics of a microscopic model, with the help of com-
puter 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 PDE’s. It describes the fluid 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
fluid sufficiently well. Whenever the fluid exhibits a physical behavior that cannot
be described by the macroscopic model, new macroscopic variables and PDE’s 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 fluid. In that
case, the observed macroscopic physics of the fluid is an emergent property that
arises naturally from the microscopic model.
The state of a simple fluid is described by the following macroscopic variables:
the fluid density ρ,

the flow velocity ~u,
the pressure p,
the energy E and
the temperature 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 reflect this
fact, they are written in a divergence form:
∂U − →
+ ∇ · F = Q. (1.1)
∂t
Computational fluid dynamics and lattice Boltzmann 7
The most commonly used equations express the conservation of

→
−
mass: ∂t ρ + ∇ · (ρ~u) = 0, (1.2)
→
−
momentum: ∂t (ρ~u) + ∇ · (ρ~u~u + pI − τ ) = ρ~g and (1.3)
→
− →
−
energy: ∂t E + ∇ · (ρ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 fluid. For a Newtonian fluid
it has, by definition, the following form:
→
−
τ = −ν 0 ρ( ∇ · u) I + 2 ρ ν S, (1.5)
where the strain rate tensor S is defined as

1 −
→ →
−
S = ( ∇~u + ( ∇~u)T ). (1.6)
2
The parameters ν and ν 0 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
ν 0 = ν. (1.7)
3
In real gases and liquids however, the actual value of ν 0 can vary substantially from
this prediction. Equation (1.3) is commonly known as the compressible Navier-
Stokes 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 flows

In an isothermal flow, effects 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 fluid. In the following, an equation for ideal
gases is used, in which the pressure scales linearly with the density:
p = c2s ρ, (1.8)
where cs is the speed of sound.

In an incompressible fluid, the density takes a constant value ρ = ρ0 which does
not vary in time and space. The conservation law for the mass is simplified, and
states that the velocity field is divergence-less (solenoidal):
−
→
∇ · ~u = 0. (1.9)
8 Chapter 1
Eq. (1.9) is often used as a definition for fluid incompressibility. The conservation law
for the momentum is also simplified and leads to the incompressible Navier-Stokes
equation, often called “Navier-Stokes equation” for short:
~ u=−1−
∂t~u + (~u · ∇)~
→
∇p + ν∇2~u. (1.10)
ρ0
In order to find 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 incom-
pressible flow 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 field remains divergence-less during its time evolution.
1.2.3 Dimensionless formulation

Before the governing equations for fluid 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 PDE’s 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 flow configuration. The
length l0 could for example stand for the size of an obstacle which is immersed in
the fluid, and t0 could be the time needed by a passive scalar in the fluid 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
→
− →
−
change of variables for ~u = l0 /t0~u∗ , ∂t = 1/t0 ∂t∗ , ∇ = 1/l0 ∇ ∗ and p = ρ0 l02 /t20 p∗ is
performed. This leads to the following dimensionless Navier-Stokes equation:
→
− →
− 1 2∗ ∗
∂t∗~u∗ + (~u∗ · ∇ ∗ )~u∗ = − ∇ ∗ p∗ + ∇ ~u and (1.12)
Re
→∗ ∗
−
∇ · ~u = 0, (1.13)
where Re = l02 /(t0 ν) is the dimensionless Reynolds number. Two flows 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 flow that possesses the same
Reynolds number and the same flow geometry.
1.2.4 Further reading

A large number of books discuss the numerical resolution of PDE’s and the topic of
computational fluid dynamics. Only a few are mentioned here to guide the interested
reader.
A well written and theoretically sound introduction to the numerical treatment

of PDE’s can be found in Ref. [2]. A general introduction to the problems of com-
putational fluid dynamics is available on the web presented as a one-semester in-
troduction course [3]. An undergraduate-level, practically oriented introduction
to computational fluid dynamics with a specific emphasis on the finite difference
method is presented in Ref. [4]. Reference [5] finally gives an extended overview
on various methods of computational fluid dynamics in general, and Ref. [6] on the
finite element method in particular.
1.3 Lattice Boltzmann method
1.3.1 Discrete variables
The problems of fluid dynamics possess an infinite number of degrees of freedom,

because the considered quantities (velocity, pressure, etc.) are not just single values,
but spatially extended fields. Computers can however only handle a finite amount
of variables, and represent them with a finite accuracy. The spatially extended fields
must therefore be replaced by a finite 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 definition is confined within a regular box of finite extent. A possible
way to discretize this space is to partition it into cells. In that case, a numerical
representation of the fluid consists of a set of numbers, of which each stands for the
average value of a fluid variable over a cell. Another approach to space discretiza-
tion consists in replacing spatially extended fields by their value on a chosen finite
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 fixed 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 defined 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 = i~e0 δx + j~e1 δx + k~e2 δ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 field µ(~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 fi-
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 paren-
theses. That is, the expression f (~rijk , tn ) is used to represent the numerical approx-
imation (and not the exact value) of f at the position ~rijk and at the time step
tn .
10 Chapter 1
1.3.2 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 defined 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 de-
scribes how the q values fi defined on the same node at given time step interact.
The macroscopic variables, the density ρ and the velocity ~u are defined locally as
moments of the distribution functions:
q−1
X
ρ= fi and (1.15)
i=0
q−1
1X
~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 explicitedP
any more to keep the notationPshort: the sum over
the q elements of a variable Ei , q−1
i=0 Ei , is simply written as i Ei .
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
i = −ω(fi − fieq ). (1.17)
The local equilibrium is defined as

(eq) 1 1
fi = ρ ti 1 + 2 ~ci · ~u + 4 Qi : ~u~u , (1.18)
cs 2cs
and the tensor Q is defined as follows:
Qi = ~ci~ci − c2s I. (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
the rest particle weight t0 correspondingly, in order to recover the lattice symmetries
described in Section 2.1. In practice, a value of c2s = 1/3 is found to be numerically
most stable, and this choice is therefore most commonly adopted. In conclusion,
there are two ingredients to the definition 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, defined 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 defined by Eqs.(1.14,1.18) solves asymptotically the equa-
tions of motion for a compressible fluid, when the discrete steps δt and δx , as well
as the relative velocity ~u/cs are small enough. This fact is proved with some reser-
vations in Section 2.3. It is also shown that the relaxation parameter ω is directly
related to the dynamic shear viscosity ν of the fluid via Eq. 2.73:

1 1
ν= c2s − . (1.20)
ω 2
The formulas of the lattice Boltzmann method can be derived theoretically via
different 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 Au-
tomata. In another approach, described for example in Ref. [9], the dynamics of the
LB method is derived from the continuum Boltzmann equation. This equation con-
siders 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, defined on a d-dimensional space, is com-
monly identified by the name “DdQq lattice”. A few 2D and 3D lattice structures
that are commonly used for isothermal flows are listed in this section. For those
flows, it is sufficient 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
1
c2s = 3
4 1 1
t0 = 9
ts = 9
tl = 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
1
c2s = 3
2 1 1
t0 = 9
ts = 9
tl = 72
~c0 = (0, 0, 0)
~c1 = (−1, 0, 0) ~c2 = (0, −1, 0) ~c3 = (0, 0, −1)
~c4 = (−1, −1, −1) ~c5 = (−1, −1, 1) ~c6 = (−1, 1, −1) ~c7 = (−1, 1, 1)
~c8 = (1, 0, 0) ~c9 = (0, 1, 0) ~c10 = (0, 0, 1)
~c11 = (1, 1, 1) ~c12 = (1, 1, −1) ~c13 = (1, −1, 1) ~c14 = (1, −1, −1)
D3Q19 lattice
1
c2s = 3
1 1 1
t0 = 3
ts = 18
tl = 36
~c0 = (0, 0, 0)
~c1 = (−1, 0, 0) ~c2 = (0, −1, 0) ~c3 = (0, 0, −1)
~c4 = (−1, −1, 0) ~c5 = (−1, 1, 0) ~c6 = (−1, 0, −1)
~c7 = (−1, 0, 1) ~c8 = (0, −1, −1) ~c9 = (0, −1, 1)
~c10 = (1, 0, 0) ~c11 = (0, 1, 0) ~c12 = (0, 0, 1)
~c13 = (1, 1, 0) ~c14 = (1, −1, 0) ~c15 = (1, 0, 1)
~c16 = (1, 0, −1) ~c17 = (0, 1, 1) ~c18 = (0, 1, −1)
D3Q27 lattice
The D3Q27 lattice uses lattice vectors of three different 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 .
1
c2s = 3
8 2 1 1
t0 = 27
ts = 27
tm = 54
tl = 216
~c0 = (0, 0, 0)
~c1 = (−1, 0, 0) ~c2 = (0, −1, 0) ~c3 = (0, 0, −1)
~c4 = (−1, −1, 0) ~c5 = (−1, 1, 0) ~c6 = (−1, 0, −1)
~c7 = (−1, 0, 1) ~c8 = (0, −1, −1) ~c9 = (0, −1, 1)
~c10 = (−1, −1, −1) ~c11 = (−1, −1, 1) ~c12 = (−1, 1, −1) ~c13 = (−1, 1, 1)
~c14 = (1, 0, 0) ~c15 = (0, 1, 0) ~c16 = (0, 0, 1)
~c17 = (1, 1, 0) ~c18 = (1, −1, 0) ~c19 = (1, 0, 1)
~c20 = (1, 0, −1) ~c21 = (0, 1, 1) ~c22 = (0, 1, −1)
~c23 = (1, 1, 1) ~c24 = (1, 1, −1) ~c25 = (1, −1, 1) ~c26 = (1, −1, −1)
1.3.4 Further reading

It is more difficult to find well written textbooks on the lattice Boltzmann method
than on other topics of computational fluid 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 back-
grounds 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 specific 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 specific 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 fluid 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 verified:
X X X
(a) ti = 1 (c) ti ciα ciβ = c2s δαβ (e) ciα ciβ ciγ ciδ = (2.1)
i i i
c4s (δαβ δγδ + δαγ δβδ + δαδ δβγ )
X X X
(b) ti ciα = 0 (d) ti ciα ciβ ciγ = 0 (f ) ti ciα ciβ ciγ ciδ ci = 0.
i i 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 fun-
damental to the presented BGK lattice Boltzmann models. There exist however
specific lattice Boltzmann models that work on lattices with weaker symmetries.
A case in point is the model presented in Ref. [15] that works fine on a 13-velocity
16 Chapter 2
3D lattice (D3Q13), although this lattice lacks sufficient symmetries to be used

with BKG dynamics. On the other hand, there exist models that require addi-
tional 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 defined on the vector space Rq . The first of those lattice vectors is called C 0
and is defined through
Ci0 = 1 ∀i = 0 · · · q. (2.2)
It follows a set of d vectors Cα1 defined as
1
Ci,α = ciα for α = 0 · · · d − 1 and i = 0 · · · q. (2.3)
Finally, a set of d2 lattice vectors Cαβ

2
is introduced that obey the relation
2
Ci,αβ = Qiαβ for α, β = 0 · · · d − 1 and i = 0 · · · q, (2.4)
where the tensor Q is defined 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) defines
therefore only d(d + 1)/2 independent vectors, instead of d2 . A scalar product
between two vectors a and b in Rq is finally defined as follows:
X
ha|bi = ti ai bi . (2.5)
i
With this, it is easily concluded from Eq. (2.1) that the lattice vectors defined above
are orthogonal to each other, but not necessarily unitary:

0 1
C |C = 0, (2.6)

0 2α
C |Cαβ = 0, (2.7)

1 2
Cα |Cβγ = 0, (2.8)

0 0
C |C = 1, (2.9)

1 1
C |C = c2s I and (2.10)

2 2 4
Cαβ |Cγδ = cs (δαγ δβδ + δαδ δβγ ). (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γδ = c4s (Tαβ + Tβα ) (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 conse-
quently skipped. All hydrodynamic variables in a LB model, which were previously
introduced as moments of the distribution functions, are now reinterpreted as pro-
jections of g on the lattice vectors. The mass density for example is computed as
ρ = hC 0 |gi, and the momentum as ρ~u = hC 1 |gi. 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 identified 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 ap-
proaches 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 defined by the evolution equa-
tion (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)
≈ (∂t + ∇ ~ · c~i + −
~ · c~i ) fi + 1 (∂t2 + 2∂t ∇ →−→
∇ ∇ : ~ci~ci ) fi . (2.14)
2
In order to relate the LB equation with a macroscopic PDE, it is necessary to
formally separate different time scales. In this way, physical phenomena occurring
at different scales are discussed separately and contribute individually to the final
equations of motion. To obtain this, the time derivative is expanded in terms of a
formal, “small” parameter :
∂t = ∂t1 + 2 ∂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 first-order term:
→
− →
−
∇ = ∇ (1) + O(2 ). (2.16)
The particle distribution functions are likewise expanded, starting with a zeroth-
order contribution f (0) :
(0) (1) (2)
fi = fi + fi + 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) (2)
Ω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 suffi-

cient to find 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) ~ 1 ·~ 1 2 2 ~ 1 ·~ 1 −→− → (0) (1)
Ωi +2 Ωi = ( ∂t1 +2 ∂t2 +∇ ci + 2 ∂t2 + ∂t1 ∇ ci + 2 ∇1 ∇1 : ~ci~ci )(fi + fi ).
2 2
(2.19)
Note 2.3: Knudsen number The parameter is often identified as the Knud-
sen 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 flows. It does not apply well for example
to dilute gases or microfluids.
2.1.3 Conservation laws

The macroscopic variables of the flow are defined as moments of the particle distri-
bution 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 definition of the zeroth-
order scalar moment ρ, the first-order vector moment ~j and the second order tensor
moment Π:
X
ρ= fi , (2.20)
i
X
~j = ~ci fi and (2.21)
i
X
Π= Qi fi . (2.22)
i
The dynamics of a flow 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 fluid 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
X
Ωi = 0, (2.23)
i
and by the requirement that first-order contributions to ρ vanish:

X X (0)
ρ= fi = fi . (2.24)
i i
This latter condition can be understood intuitively as follows: the conserved vari-
ables 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 final aim,
namely to find 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 different scales and 2 , and using Eq. (2.19), yields
X (1) X (0) − → X (0) →
−
Ωi = ∂t1 fi + ∇1 · ~ci fi = ∂t1 ρ + ∇1 · ~j (0) = 0 (2.25)
i i i
and
X (2)
X (1)
X → X (1) 1 2 X (0)
− (0)
Ωi = ∂t1 fi + ∂t2 + ∇1 ·
fi ~ci fi + ∂t1 fi
i i i i
2 i
→ X (0) 1 −
− →−→ X (0)
+ ∂t1 ∇1 · ~ci fi + ∇1 ∇1 : (Qi + c2s I)fi (2.26)
i
2 i
→
− 1 2 →
− 1−→−
→ 1
= ∂t2 ρ + ∇ 1 · ~j (1) + ∂t1 ρ + ∂t1 ∇1 · ~j (0) + ∇1 ∇1 : Π(0) + c2s ∇12 ρ.
2 2 2
Equations (2.25) and (2.26) are combined to eliminate the second-order time deriva-
2
tive ∂t1 ρ:
→
− 1 − → 1−→−
→ 1
∂t2 ρ + ∇1 · ~j (1) + ∂t1 ∇1 · ~j (0) + ∇1 ∇1 : Π(0) + c2s ∇12 ρ = 0. (2.27)
2 2 2
First order moment conservation. At the first order moment, a source term
F~ = F~ (1) (2.28)
is added. It will be identified in the following as an external force acting on a fluid.

Depending on the specific 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-diffusion equation with source terms. The first order
collision invariant reads X
~ci Ωi = F~ . (2.29)
i
In analogy with the zeroth-order moment, the first-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 =
X X (0) F~
~ci fi = ~ci fi − . (2.30)
i i
2
20 Chapter 2
The O() scale contribution to Eq. (2.29) reads

X (1)
X (0) − → X (0)
~ci Ωi = ∂t1 ~ci fi + ∇1 · (Qi + c2s I)fi
i i i (2.31)
→
− →
−
= ∂t1~j (0) + ∇1 · Π(0) + c2s ∇1 ρ = F~ (1) .
Two O(2 ) terms appearing in Eq. (2.27) can now be evaluated and inserted into
this equation:
→
− →
− →−
− →
∂t1 ∇1~j (0) = ∇1 · F~ (1) − ∇1 ∇1 : Π(0) − c1s ∇21 ρ by (2.31), and (2.32)
−
→ ~ (1) 1−→
∇1 · j = − ∇1 · F~ by (2.30). (2.33)
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:
X (2)
X (1) X (0) − →X (1)
~ci Ωi = ∂t1 ~ci fi + ∂t2 ~ci fi + ∇1 (Qi + c2s I)fi
i i i i
1 2 X (0) → X
− (0) 1−→−
→
+ ∂t1 ~ci fi + ∂t1 ∇1 · (Qi + c2s I)fi + ∇1 ∇1 : R(0)
2 2
i i (2.36)
−
→ 1 →
−
= ∂t1~j (1) + ∂t2~j (0) + ∇1 · Π(1) + ∂t1
2 ~ (0)
j + ∂t1 ∇1 · Π(0)
2
→
− 1−→−→
+ c2s ∂t1 ∇12 ρ + ∇1 ∇1 : R(0) = 0,
2
P
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),
→
− →
−
∂t1 j = ∂t1 F~ (1) − ∇1 · Π(0) − c2s ∇1 ρ ,
2 ~ (0)
(2.37)
which simplifies Eq. (2.36) as follows:

→
− 1 − → 1 →
− 1−→−→
∂t2~j (0) + ∇1 · Π(1) + ∂t1 ∇1 · Π(0) + c2s ∂t1 ∇1 ρ + ∇1 ∇1 : R(0) = 0. (2.38)
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:
→
− (0) −→
∂t~j (0) + ∇ · Π + c2s ρ I + ∂t1 Π + c2s ρ I + ∇1 · R(0) = F~ . (2.39)
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) + c2s ρ I is identified with ρ ~u~u + p I,
and the remaining O(2 ) terms with −τ .
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 com-
puted 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 do-
main 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 floating point representations in a computer. By
“exact conservation” it is meant that none of the typical error terms that scale
at an order O(δx3 ), O(δt3 ), O(3 ) or O(Ma 3 ) is present in the conservation laws.
2.2 Advection-diffusion: Chapman-Enskog ansatz

The advection-diffusion equation describes the evolution of a scalar quantity ρ that
is diffused and exposed to the advective influence of an external term ~u:
→
−
∂t ρ + ∇ · (ρ~u) = d ∇2 ρ, (2.40)
where d is the diffusivity constant. The scalar ρ represents for example the density
of a component in a multi-component fluid, or the temperature distribution is a
fluid. 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-diffusion LB model can be written as a
BGK term, a relaxation towards a local equilibrium distribution:
1 (eq)

Ωi = − fi − fi + ti ~ci · δ~j. (2.41)
λ
Here, λ is the relaxation time, and δ~j = δ~j (1) an external source term for the
non-conserved first order moment. It is introduced formally and will be used later
to correct errors in the numerical model. The equilibrium distribution reads
1
fieq = ρ ti (1 + ~ci · ~u). (2.42)
c2s
It is clear that this collision operator verifies 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:

→ ~ (0)
− 2−
→ ~ (1) 1 ~ (0) 1 − → (0) 1 2−→
∂t ρ + ∇1 · j = − ∇1 · j + ∂t1 j + ∇1 · Π + cs ∇1 ρ . (2.43)
2 2 2
22 Chapter 2
In order to find the unknown term in this equation, a scale separation for the par-
ticle distribution function fi must be explicitly defined. In the Chapman-Enskog
approximation, the model is expanded around the local equilibrium term. That is,
the equilibrium term is identified 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 (2.45)
(0)
Π = 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 first-order
derivatives need to be retained at the O() scale:

(1) →
− 1
fi = −ti λ (∂t1 + ∇1 · ~ci ) ρ 1 + 2 ~ci · ~u + λ ti~ci · δ~j
cs
(2.48)
λ →
− →
− λ ~
= − 2 ti Qi : ∇1 (ρ~u) − λ ti ~ci · ∇1 ρ − 2 ti ∂t1~ci · (ρ~u) + λ ti ~ci · δ j.
cs cs
Projection on the first-order base vectors C 1 yields

! !! !
~j (1) =
X (1) 1 X →
− X X
~ci fi = −λ 2 ∂t1 ti~ci~ci ρ~u + ∇1 · ρ ti~ci~ci +λ ~ci~ci δ~j
i
cs i i i
→
−
= −λ ∂t1 (ρ~u) + c2s ∇1 ρ + λ c2s δ~j.
(2.49)
Eqs. (2.49), (2.45) and (2.46) are inserted into Eq. (2.43), which finally yields the
O(2 ) scale of the dynamics:
→
− →
− →
− →
−
∂t ρ + ∇ · (ρ~u) = d ∇ 2 ρ + d/c2s ∂t ∇ · (ρ~u) − λ c2s ∇ · δ~j , (2.50)
| {z } | {z }
unwanted term correction term
where the diffusivity constant d has been defined as

2 1
d = cs λ − . (2.51)
2
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 find the error term. In Section 2.4, numerical evidence for the negative
influence 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-diffusion 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
c2s is defined as c2s = 1/2. Another possible choice is the D2Q5 model that includes
the zero-speed velocity ~c0 . In this model, the constant c2s is a free parameter. The
lattice weights take the value t0 = 1 − 2 c2s for the zero-speed velocity, and ti = c2s /2
for the four other ones. It is common to choose c2s = 1/3, with corresponding lattice
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, c2s = 1/3) and D3Q7 (t0 = 1 − 3 c2s , ti = c2s /2 for
i > 0).
2.3 Fluid flows: Chapman-Enskog ansatz

The arguments leading to a LB model for fluid flows are similar to those devel-
oped in the previous section for the advection-diffusion 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 first and
second order vectors C 1 and C 2 :
a b
FT i = ti 2
~ci · F~ + ti 4 Qi : (F~ ~u + ~uF~ ) and (2.53)
cs 2 cs
c
CT i = ti 4 Qi : δΠ. (2.54)
2 cs
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 cor-
rection to numerical errors of the force term, due to the second order time derivative.
This approach has first been described in Ref. [18]. Finally, a formal term δΠ is
introduced. It is used in Section 3.2 to correct a numerical deficiency of the model.
This term is introduced, just like the correction term δ~j in Eq. (2.41), at the level
of the first non-conserved moment. In this way, it does not interfere with the ex-
act 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 con-
servation 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 defined in Eq. 1.18. Actually, the dynamics de-
scribed 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 con-
tains 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 fluid incompressibility, the density vari-
ations 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 first order moment (momentum) reads
X ω ~
~ci Ωi = a + F by Eqs. (2.52) and (2.30). (2.56)
i
2
For Eq. (2.29) to be respected, the value of the constant a must be set to
ω
a=1− . (2.57)
2
The governing equations for the LB fluid 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-diffusion model. This time, all symmetry properties listed in Eq. (2.1) are
used. The zeroth order terms read
~j (0) = ρ~u, (2.58)
Π(0) = ρ ~u~u and (2.59)
→
− −→ →
− →
−
∇1 · R = c2s ∇1 (ρ~u) + ( ∇1 (ρ~u))T + ∇1 · (ρ~u) I . (2.60)
At the first order with respect to , it follows from Eq. (2.52) that

(1) ti →
− 1 1 1
fi = − (∂t1 + ∇1 · ~ci ) ρ 1 + 2 ~ci · ~u + 4 Qi : ~u~u + (FT i + CT i )
ω cs 2cs ω
ti 1 1
= − ∂t1 ρ + 2 ∂t1 (~ci · ρ~u) + 4 ∂t1 (Qi : ρ~u~u)
ω cs 2cs
→
− 1 →
− →
− 1
+~ci · ∇1 ρ + 2 ~ci~ci : ∇1 (ρ~u) + (~ci · ∇1 )(Qi : ρ~u~u) + (FT i + CT i ).
cs ω
(2.61)
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
∂t1 (~
c i · ρ~
u ) =
1
~
c i · ~ (1) − 1 ~ci · −
F
→
∇ 1 : ρ~
u ~
u − ~
c i ·
→
−
∇1 ρ by (2.31,2.59),
c2s c2s c2s
(2.63)
1 1
4
Qi : ∂t1 (ρ ~u~u) = 4 Qi : ∂t1 (ρ~u)~u + ~u ∂t1 (ρ~u) − (∂t1 ρ)~u~u
2 cs 2 cs | {z }
O(Ma 3 )
1
= Q : ∂t1 (ρ~u)~u (because Q is symmetric)
c4s i
1 →
− →
−
= 4 Qi : F~ (1) − ∇1 · (ρ~u~u) −c2s ∇1 ρ ~u by (2.31)
cs | {z }
O(Ma 3 )
1 2−
→

~ (1)
= 4 Qi : F − cs ∇1 ρ ~u. (2.64)
cs
Inserting all these terms into Eq. (2.61) yields

(1) ti →
− →
− 1 →
−
fi = − 2 Qi : ρ ∇1~u − ~ci ∇1 : ρ ~u~u + 2 (~ci · ∇1 )(Qi : ρ ~u~u)
cs ω 2cs
(2.65)
1 ti (b − 1)ti c ti
− 2 ~ci · F~ + 4
Qi : (F~ ~u + ~uF~ ) + 4 Qi : δΠ.
2 cs 2 cs ω 2 cs ω
When this equation is projected on the vectors of C 2 , using Eq. (2.12), only two
terms are left:
X (1) 2 c2 b−1 ~ c
Π(1) = Qi fi = − s ρ S + (F ~u + ~uF~ ) + δΠ, (2.66)
i
ω ω ω
where the symmetric tensor S is defined by Eq. (1.6) on page 7. In order to find the
last missing term, the time derivative of the Π(0) tensor, the same approximations
are used as those that lead to Eq. (2.64):
−→ →
−
(0) ~ (1) ~ (1) 2
∂t1 Π = ∂t1 (ρ~u~u) − F ~u + ~uF − cs ( ∇1 ρ)~u + ~u( ∇1 ρ) . (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 oth-
ers an unwanted term that depends on the external force F~ , introduced by Eqs. (2.66)
and (2.67):
2 c2

1 b−1 1 1 − →
(1)
Π + ∂t1 Π = (0)
+ (F~ ~u + ~uF~ ) − s ρ S + c2s ∇ · (ρ~u)I. (2.68)
2 ω 2 ω 2
This force term is eliminated by setting the constant b to
ω
b=1− . (2.69)
2
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:
→
−
∂t ρ + ∇ · (ρ~u) = 0 and (2.70)
→
−
∂t (ρ~u) + ∇ · c2s ρI + ρ~u~u − τ = F~ ,

(2.71)
where the deviatoric stress tensor τ reads
τ = 2 νρ S, (2.72)
and the dynamic shear viscosity ν is related to the relaxation parameter ω as follows:

2 1 1
ν = cs − . (2.73)
ω 2
This set of PDE’s is identical to the equations of fluid motion, introduced in Sec-
tion 1.2.1. The bulk viscosity ν 0 can however not be chosen as a free parameter in
the LB method, and it takes the fixed value ν 0 = 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 finally rewritten using the
choice of the constants a and b derived in the present section:

ω (~ci − ~u) ~ci · ~u
FT i = 1 − ti + 4 ~ci · F~ . (2.74)
2 c2s cs
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 fluids 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 as-
sociated macroscopic PDE is explicited via a finite development of the difference
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.

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 clarified
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 definition, unitary in the dimensionless system of units:
T ∗ = 1 and L∗ = 1. This leads to the definition of the discrete time and space steps:
T∗ 1
δt = = and (2.75)
T T
L∗ 1
δx = = . (2.76)
N −1 N −1
The term N − 1 in the definition 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:
N −1 δt
U= = (2.77)
T δx
Furthermore, the fluid viscosity ν = c2s (τ − 1/2), expressed in lattice units, is related
to the Reynolds number. To simplify the discussion, only the incompressible Navier-
Stokes equation (1.10) is considered. By casting this equation into its dimensionless
form, given by Eq. (1.12), one obtains:
1 δt
ν= . (2.78)
Re δx2
Note 2.6: Dimensionless or not? 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 different systems
of units. One of them is the “dimensionless” system: it might be more appro-
priate 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 lat-
tice 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 Chapter 2
2.4.2 Accuracy of LB methods

The signification 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 defining 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:
rP
u(~r) − ~v (~r)|)2
r (|~
~
< κk , (2.79)
Nd
where the sum is taken over all sites of a grid of size N d . The time accuracy can
be defined similarly for flows 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 κ0 by which the
error between the numerical solution ~u(~x) and the exact solution ~v (~x) is constrained
as follows: s
u(~r, t) − ~v (~r, t)|2
P
1 r |~
X
~
< κ0k , (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 PDE’s in which the space derivatives are replaced by a finite 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 specific models

Some care must however be taken when the LB method is used to simulate incom-
pressible flows. 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
the error due to the compressibility effects scales like EMa = O(δt2 /δx2 ). There are
therefore in total three terms that contribute to the overall error E:
E = Eδx + Eδt + EMa = O(δx2 ) + O(δt2 ) + O(δt2 /δx2 ). (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
be readapted for the overall error to decrease at a rate of O(δx2 ). That is, the EMa
error term must scale as EMa = O(δx2 ), and the time and space step are therefore
related through δt ∼ δx2 . This in turn means that the time discretization error really
behaves as in numerical methods with first 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 efficiency of the LB model.
Note 2.7: Keeping ω fixed

When benchmarks are executed on several simu-
lations with different 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 re-
main identical from one grid to another. This choice is a priori surprising,
because it is based on a microscopic property ω whose significance 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 Sec-
tion, it is however possible to better understand this choice. Indeed, Eq. (2.78)
shows that the parameter ω is proportional to δt /δx2 . Keeping ω constant thus
amounts to varying δt at a rate proportional to δx2 . It has been argued in the pre-
vious 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 flows 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 fluid model. This drawback can luckily be
fixed 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 identified 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:
δt
Ma ∼ . (2.82)
δx
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 different 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-diffusion 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 fluid dynamics. To see this, it suffices 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):
−∗
→ ∗−
→∗2 δt2 1 ∗ ∗ −
→
∂t∗ ρ ∗
+ ∇ · (ρ~u ) = d ∇ ρ + 2 2 d ∂t ∇ ∗ · (ρ~u∗ ), (2.83)
δx cs
where d∗ = δx2 /δt d is the dimensionless diffusion parameter. Given that c2s is an
immutable lattice constant, this equation contains an error term ELE due to lattice
effects that scales like δt2 /δx2 . As in Eq (2.81), the overall error contains therefore
three contributing terms:
E = Eδx + Eδt + ELE = O(δx2 ) + O(δt2 ) + O(δt2 /δx2 ). (2.84)
The effect of those error terms on actual numerical simulations is discussed in Sec-
tion 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 compress-
ible, isothermal flows at small Mach numbers. This conjecture has however
never been proved formally, and it makes sense only when the widely un-
known correction term in Section 3.2 is used. Furthermore, the discrete
time and space steps cannot be chosen freely in compressible flows, because
they are coupled with a physically significant value, the Mach number. This
problem does not appear when incompressible flows are simulated, but in
that context, the LB model behaves like a numerical scheme with first order
accurate time stepping.
Chapter 3
Corrections to the BGK dynamics
In Chapter 2, two lattice Boltzmann models are analyzed, one that describes advection-
diffusion phenomena, and one that simulates fluid motion. In both models, a cor-
rection term is formally introduced, which offers an additional degree of freedom
for handling deficiencies of the numerical models. Those degrees of freedom are
exploited in this chapter. In the advection-diffusion model, the correction term δ~j
is adjusted so as to eliminate the lattice error ELE (cf. Section 2.4). In the fluid
model, an appropriate choice of the correction term δΠ leads to a numerical scheme
in which the bulk viscosity ν 0 is a free parameter (it was shown in Section 2.3 that
ν 0 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 deficiency
of both numerical models can be eliminated by explicitly referring to the value of
the “rest-particle” distribution function f0 .
3.1 Advection-diffusion revisited

3.1.1 Correction term for second order accurate time step-
ping
In Section 2.2 a BGK model for the modeling of advection-diffusion 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 specified
correction vector δ~j. In the present section, this vector is chosen in such a way as
to cancel out the unwanted term:

1 →
− →
−
λ− ∂t1 ∇1 · (ρ~u) − λ c2s ∇1 · δ~j = 0, (3.1)
2
which leads to the condition

1 1 1 1
δ~j = 2 1 − ∂t1 (ρ~u) = 2 1 − (ρ ∂t1~u − ~u ∂t1 ρ). (3.2)
cs 2λ cs 2λ | {z } | {z }
~j1 ~j2
The correction term δ~j contains two contributions. The first, 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 first order
accurate estimate of those derivatives is sufficient for the present purpose, because
the term ∂t1 acts at a O() scale. In Eq. (2.19), this scale is related to the first order
terms of the finite Taylor series.
The value of ~j1 can for example be found by a finite difference approximation
to the time derivative. The difference of the velocity between two successive time
steps is expanded in a Taylor series and scale separated as follows:
1 2
~u(t + 1) − ~u(t) = ∂t1~u + 2 (∂t2 + ∂t1 )~u + O(3 ), (3.3)
2
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 first-order finite
difference suggested by Eq. (3.4) is sufficient 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-diffusion
equation is the analysis of thermal flows using the Boussinesq approximation. In that
case, the external vector field ~u is obtained from a numerical solver for the Navier-
Stokes equation. It is common for such solvers to keep in memory two successive
time steps of the velocity field, and the finite difference in Eq. (3.4) can be evaluated
efficiently.
(1)
The value of ~j2 is directly found from the off-equilibrium part f0 of the rest-
particle distribution function. As a matter of fact, Eq. (2.48) leads to the following
relation:
1 (1)
∂t1 ρ = − f . (3.5)
λt0 0
Equations (3.4) and (3.5) are finally combined with Eq. (3.2) as follows:

~ 1 1 1 (1)
δj = 2 1 − ρ (~u(t + 1) − ~u(t)) − f ~u . (3.6)
cs 2λ λ t0 0
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 fine with the D2Q5 lattice.
Note 3.1: Efficiency

When the correction term in Eq. (3.6) is taken into ac-
count, the LB method for advection-diffusion 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 excep-
tional 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 five 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 verification

This section presents the results of a 2D numerical experiment on which the LB
model for advection-diffusion (2.41,2.42) is tested, both in its original BGK formu-
lation with δ~j = 0, and with the modification term suggested in Eq. (3.6). The
theoretical discussion of the error terms in Eq. (2.84) is numerically verified, and
the benefit of the correction term δ~j on the accuracy of the model is underlined.
In the numerical example, a fixed, divergence free and time independent velocity
field ~u is imposed, and the advection-diffusion equation for the scalar ρ is solved
under the advective influence of ~u. The velocity field is defined on a square domain
whose side length serves as reference for the definition of a dimensionless length
scale. It corresponds to a Poiseuille profile whose maximum velocity umax defines
the reference scale for the velocity:
u~∗ (x~∗ ) = (u∗ (x∗ , y ∗ ), v ∗ (x∗ , y ∗ )), where (3.7)

u∗ (x∗ , y ∗ ) = 4 y ∗ (1 − y ∗ ) and
v ∗ (x∗ , y ∗ ) = 0.
The domain is periodically extended in both the x- and the y-direction. The diffused
scalar ρ is initially defined through a uniform Gaussian distribution:
1 |~r∗ −~c|2
∗ −
ρini (~r ) = √ e 2 σ2 , (3.8)
2π σ
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
sufficiently small, they are successively decreased until time and space convergence
is reached, that is, until a further refinement does not influence any more the con-
clusions 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 defined through a l2 space norm, just as in Eq. (2.79). Similarly, the error of
a simulation on a full time interval is defined 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 defined in Section 2.4, because the flow in not peri-
odic. 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)
0 −1
10 10
BGK model BGK model
Corrected model Corrected model
−2 Slope −2 Slope −1
10
Relative error
Relative error
−4 −2
10 10
−6
10
−8 −3
10 2 3 4
10 2
10 10 10 10
1/δt Grid resolution N
Figure 3.1: Accuracy of the BGK and the corrected advection diffusion models on
a chosen 2D problem. (a) Fixed diffusivity d∗ = 0.01 and grid resolution N = 100,
varying time step δt. (b) Fixed diffusivity d∗ = 0.05 and time step δt = 0.002,
varying grid resolution N .
is kept constant, and the time step is progressively refined. 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,
the space error Eδx = O(δx2 ), the time error Eδt = O(δt2 ) and an additional lattice
effect ELE = O(δt2 /δx2 ) that is present only in the uncorrected BGK model. The
space error Eδx is not visible on Figure 3.1 (a), because the lattice resolution has
been chosen deliberately high enough. The remaining error terms Eδt and ELE scale
both at a same rate O(δt2 ). This is confirmed on the graph, on which the error
visibly decreases at the same rate for both numerical methods. However, the offset
between the two curves shows that the term ELE has a more severe influence on the
overall error than Eδt . On Figure 3.1(b), the time step δt is kept constant, at a value
sufficiently small for the time error Eδt to be negligible, and the lattice resolution is
progressively increased. This time, the BGK simulation contains a decreasing error
Eδx = O(δx2 ), and an increasing error ELE = O(1/δx2 ). This is clearly reflected 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
side remark, the rate of decrease is proportional to δx and not δx2 , as one might have
expected. It must be again underlined that those measurements of the accuracy
may not be interpreted quantitatively, because the flow 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 verified quantitatively.
3.2 Bulk viscosity for fluid flows

3.2.1 Numerical error in incompressible flows
Equation (2.81) on page 28 shows that there are three contributions to the overall
error of a simulation when an incompressible fluid is being simulated: a space error
Eδx = O(δx2 ), a time error Eδt = O(δt2 ) and a compressibility error EMa = O(δt2 /δx2 ).
In the present section, those theoretical considerations are first underlined by a
numerical experiment. Then, a modified BGK model is introduced in which the
error term cancels.
The experiment, known under the name of Womersley flow, is now shortly intro-
duced and presented in a system of dimensionless variables. Just like the Poiseuille
flow used in Section 3.1, the Womersley flow takes place in a channel whose width
determines the length scale for the dimensionless variables. A Poiseuille flow is a
flow parallel to the channel boundaries with parabolic velocity profile that is driven
either by an constant body force or by a constant pressure gradient DP . A pressure
driven Poiseuille velocity profile has the following form:
Re
u∗x (y ∗ ) = DP (y ∗2 − y ∗ ). (3.9)
2
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
∂p 8
∗
= DP = − . (3.10)
∂x Re
In a Womersley flow, the pressure gradient is varied periodically with a frequency
ω ∗ as follows:
∂p
= A cos(ω ∗ t∗ ). (3.11)
∂x∗
The resulting velocity profile is time dependent and substantially more complex than
the Poiseuille profile 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 flow, the Poiseuille flow is taken as a reference, and the amplitude A in
Eq. 3.11 is identified with the pressure gradient DP in Eq. 3.10: A = −8/Re. Fur-
thermore, a parameter α is introduced that is known under the name of Womersley
number and defined as r
Re ω ∗
α= . (3.12)
4
With those definitions, the analytical Womersley profile reads
  √ ∗ 1

DP iωt  cosh 2 y − (α + iα)
u∗x = Real  e 1− √ 2  , (3.13)
iω cosh 2
(α + iα)
2
where i is the imaginary number, and the term Real indicates that the real compo-
nent of the right hand side expression must be extracted.
All simulations are executed with an unmodified BGK model, defined 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 flow 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 insufficient for the implementation of the pressure gradient).
36 Chapter 3
0.1
t*/T = 0.2 *
t /T = 0.3
0.08
0.06
Velocity t*/T = 0.1
0.04
0.02
t*/T = 0
0
−0.02
0 0.2 0.4 0.6 0.8 1
Position y* in the channel
Figure 3.2: Womersley profile, 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 profile
is seriously flattened, compared to a Poiseuille profile 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
must be varied like δx2 so as to cancel compressibility effects. The results of applying
this procedure, with δt = 2/5 δx2 , 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 first
order accurate time discretization, when used for incompressible flows. Figure 3.3(b)
shows what happens when the time step is fixed to a value δt = 2/5 100−2 = 0.004.
When the grid resolution is refined, the error, dominated by the space error Eδx ,
decreases with a second order rate. At some point however, the compressibility error
EMa , which scales like 1/δx2 , 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 diffusion model could be eliminated by an appropri-
ate choice of the correction term CT i in Section 3.1. It would be tempting to take
a similar approach for the fluid 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 flows with an adequate accuracy. Though a seemingly
incompressible fluid model has been presented in Refs. [19, 20], this model has been
(a) (b)
−1 −2
10 10
−2
Relative Error
10
Relative error
N −2
−3
C 10 N −2
on e
ve Co
nv enc
−3 rge e erg
10
nc
rge
nc 2 Div
e e N
−4 −4
10 1 2 3 10
10 10 10 20 40 60 80 100
Grid resolution N Grid resolution N
Figure 3.3: Accuracy of the BGK scheme compared to the analytical solution of an
incompressible Womersley flow. (a) δt = 2/5 δx2 (b) δt = const. = 2/5 · 100−2
shown to be erroneous (see e.g. Ref. [21]), as it eliminates compressibility effects

only in stationary flows. For that reason, the LB method is most efficiently used to
simulate weakly compressible flows. 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 offer 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 inclu-
sion 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 spe-
cial care. The section ends with the suggestion of an alternative correction term,
which achieves the same effect as the one proposed in [22] but is more efficient 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, (3.14)
d
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:
c 2 c2 −
→
Tr (Π(1) ) = Tr (δΠ) − s ρ ∇1 · ~u − F~ · ~u., (3.16)
ω ω
Inserting Eq. (3.15) leads to
2

1 c 2 c →
−
Tr (Π ) = (1)
c k − 1 F~ · ~u − s
ρ ∇1 · ~u . (3.17)
1− ω
ω ω
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
2

c 1 c c 2 c →
−
δΠ = 1+ k − 1 F~ · ~u − s
ρ ∇1 · ~u I. (3.18)
ω 1 − c/ω ω d ω ω
The force term is eliminated by choosing the constant k to be
1
k= . (3.19)
1 + ωc
Furthermore, by setting
1
c=ω 2 c2s
, (3.20)
1− d ω ν0
Eq. (3.18) becomes
c →
−
δΠ = ν 0 ρ ∇1 · ~u I. (3.21)
ω
The momentum conservation law derived from the modified 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 final form of the correction term, Eq. (3.14) is inserted into
Eq. (2.54). Using the identity Qi : I = |~ci |2 − c2s d, one obtains
c ti 2 2
(1) ~

CT i = |~
c i | − c s d Tr (Π ) + k F · ~
u , (3.22)
2 d c4s
where the constants c and k are defined by Eqs. (3.20) and (3.19). The expression
Tr (Π(1) ) is evaluated using the definition of Π in Eq. (2.22), and the value of Π(0)
for the BGK model in Eq. (2.59). One finds
X
Π(1) = ~ci~ci fi − ρ ~u~u − c2s I, and consequently, (3.23)
i
X
(1)
Tr (Π ) = |~ci |2 fi − ρ |~u|2 − c2s d. (3.24)
i
It is interesting to note that the calculations required in Eq. (3.24) can be cir-
cumvented by referring explicitly to the rest particle distribution function f0 , as was
done in Section 3.1 for the advection-diffusion model. Indeed, the following relation
is directly deduced from Eq. (2.65):
(1) t0
f0 = f0 − ρ t0 = − Tr (Π(1) ). (3.25)
2 c2s
The correction term CT i can therefore alternatively be written as follows:
2 c2s (1)

c ti 2 2

~
CT i = |~ci | − cs d − f + k F · ~u . (3.26)
2 d c4s t0 0
It is reminded that the adjustable bulk viscosity ν 0 enters the equation via the
constant c, defined in Eq. (3.20), and that the constant d corresponds to the number
of space dimensions (2 or 3) of the system.
3.3 Alternative LB models

3.3.1 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 Ω, defined 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 off-equilibrium particle distribution functions:
X
fi (~x + ~ci , t + 1) − fi (~x, t) = Ωij (fj − fjeq ). (3.27)
j
The restriction to the off-equilibrium part of the distribution functions is moti-

vated theoretically by the multi-scale analysis in Section (2.1), in which the scale
of the collision operator is formally identified with the scale of the off-equilibrium
distribution functions. Using the definition of the vector space Rq introduced in Sec-
tion 2.1.1, with the scalar product of Eq. (2.5) and the definition of the unweighted
particle distribution functions g in Eq. (2.13), Eq. (3.27) is interpreted as a Matrix-
Vector 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)
where the short-hand notation giout ≡ g(~x + ~ci , t + 1) for “outgoing” particle distri-
bution 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 first populated by the 1 + d + d(d + 1)/2 vectors defined by C 0 , C 1 and C 2
in Eqs. (2.2), (2.3) and (2.4) respectively. It is now clear that the first 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 defined as Ω = D = −ω I. Equa-
tion (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 repre-
sented 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 different 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 specified 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 define a new lattice vector C 3 by
Ci3 = hi (3.32)
and two lattice vectors Cα4 by
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:

1 1 1 3
fi = ti ρ + 2 ρ ~ci · ~u + 4 Π : Q + hi N + ~ci · J~ , (3.34)
cs 2 cs 4 8
where the projections on the ghost vectors have been named hC 3 |f i = N and
hC 4 |f i = J~ respectively. It is illuminating to split Eq. (3.34) into an equilibrium
and an off-equilibrium contribution:

eq 1 1 eq
fi = ti ρ + 2 ρ ~ci · ~u + 4 Π : Q and (3.35)
cs 2 cs

1 neq 1 3
neq
fi = ti Π : Q + hi N + ~ci · J~ . (3.36)
2 c4s 4 8
With this notation, the BGK collision is written as follows:

1 neq 1 3
Ωi = −ω ti Π : Q + hi N + ~ci · J~ . (3.37)
2 c4s 4 8
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
Ωi = −ti 4 ω3 Πneq neq neq

xx Q xx + ω4 Π xy Q xy + ω5 Π yy Q yy
c
s (3.38)
1 3
+ hi ω6 N + (ω7 cix Jx + ω8 ciy Jy ) .
4 8
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 modified.
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, diagonaliz-
able collision operator. This puts the emphasis on the wrong place. Such a ter-
minology would actually imply that practically all commonly used LB models are
MRT models, including the BGK model, and several other models that were devel-
oped even before the term “MRT” existed. In particular, specific 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 fine-tuning the relax-
ation parameters.
3.3.2 Entropic LB models

The family of entropic LB models has emerged over the past decade and is suc-
cessively gaining importance. Two approaches were developed independently by
Karlin, Ansumali e.a., as documented in Refs. [28, 29], and by Boghosian e.a., doc-
umented in Ref. [30]. Those models exploit a discrete-space analogue of Boltzmann’s
H-theorem to achieve unconditional numerical stability. Numerical stability means
that the simulated values always remain within a fixed range and never exceed the
storage capacity of floating 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 under-
resolved during a short time interval, which should only affect 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
Boltzmann’s H-theorem states that each fluid 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 defined
property H = H(f (~x, t)) never increases during particle collisions. Because of nu-
merical errors, this theorem is only approximately true in BGK simulations. In
entropic models however, the collision operator is modified so as to obtain an exact
discrete version of the H theorem. This can be compared to the local conserva-
tion 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 defined 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 modified so as to ensure that the outgoing distribution func-

tions fiout are contained within the set of distribution functions evaluated pre-
viously. 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 effect is technically prevented in entropic methods by the exact H theo-
rem.
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 defined by c2s = 1/3. This
shows that apart from this factor 1/3, the characteristic flow 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 inter-
pret 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
δx δx 1 X
~u∗ = ~u = fi ~ci . (4.1)
δt δt ρ i
It is much more contrived however to initialize the particle distribution functions of a

LB simulation from a given prescription of the macroscopic flow fields. 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
= ρti 1 + 2 ~ci · ~u + 4 Qi : ~u~u
c 2cs (4.2)
s
ti →
− →
− 1 →
−
− 2 Qi : ρ ∇~u − ~ci ∇ : ρ~u~u + 2 (~ci · ∇)(Qi : ρ~u~u) + O(2 ).
cs ω 2cs
This equation shows that in order to relate a macroscopic flow field to a LB field,
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 fea-
sible, 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 flow field
are imported into a LB simulations, but also when two LB simulations, executed
on different grids, are related to each other. This becomes clear when Eq. (4.2) is
rewritten in terms of dimensionless variables:
δt2 1

δt 1 ∗ ∗ ∗
fi = ρti 1 + ~ci · ~u + 2 4 Qi : ~u ~u
δx c2s δx 2cs
→∗ ∗ δt2

ti − →∗
− ∗ ∗ 1 →∗
− ∗ ∗
− 2 δt Qi : ρ ∇ ~u − ~ci ∇ : ρ~u ~u + 2 (~ci · ∇ )(Qi : ρ~u ~u ) .
cs ω δx 2cs
(4.3)
It is concluded from this equation that different contributions to the particle distri-
bution functions scale at a different rate under a change of δx and δt . To compare
the results from two LB simulations that are run on different grids, it is therefore
necessary to split the fi into their components, and to rescale each component indi-
vidually. 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 off-equilibrium part. It is however not possible to
split the off-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 off-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 off-equilibrium distributions fi scale like δt .
This approach is chosen in Ref. [31] to formulate a grid-refinement 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 refinement procedure.
A different 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 verified in Section 4.4 on chosen numerical examples.
In the subsequent chapters, various typical problems of computational fluid dy-
namics 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 for-
mulation 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 off-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 first 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 PDE’s
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 (4.4)

1 →
− 1 ~ c
Ri = ti Qi : − 2 ρ ∇1~u − 4 F ~u + ~uF~ + 2 δΠ and (4.5)
cs ω 4 cs 2 cs ω

ti →
− 1 − → 1 ti
Ii = 2 ~ci ∇1 : ρ~u~u − 2 ~ci · ∇1 (Qi : ρ~u~u) − 2 ~ci · F~ . (4.6)
cs ω 2 cs 2 cs
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)
X
Παβ = ti Qiαβ Qiγδ Tγδ = c4s (Tαβ + Tβα ) by Eq. (2.12). (4.7)
i
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) = c4s Qi : (T + T T ) = 2 c4s Qi : T by symmetry of the tensor Q. (4.8)

(1) (1)
From this, a regularized substitute f¯i for fi is derived:
(1) ti
f¯i = Ri = 4 Qi : Π(1) . (4.9)
2 cs

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:
(1) ρ ti
f¯i ≈ − 2 Qi : S. (4.10)
cs ω
The tensor S, defined by Eq. (1.6), is a macroscopic quantity that depends only on
the gradients of the velocity field. To conclude this discussion, it is pointed out that
the set of lattice kinetic variables defined by the regularized particle distribution
(1)
functions f¯i = fieq + f¯i is equivalent to a set of macroscopic variables defined 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 (4.11)
δx
Π(1)∗ = δt Π(1) .
4.2.3 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 from the particle distribution func-
tions by means of Eq. (1.18).
(neq) (eq)
2. Determine the off-equilibrium distribution functions fi = fi − fi and the
(neq)
related tensor Π(neq) = i Qi fi
P
.
(1)
3. Derive the first order regularized distribution functions f¯i from Π(neq) through
the approximation Π(1) ≈ Π(neq) in Eq. (4.9).
4. Replace the distribution functions fi by their regularized counterpart f¯i =

(eq) (1)
fi + f¯i .
5. Execute the BGK collision and the streaming step.

In order to clarify this algorithm, the original BGK LB model (Eqs. (1.14,2.52))
is first reformulated as follows:
(eq)
fi (~r + ~ci , t + 1) = fi (~r, t) − ω (fi (~r, t) − fi (ρ, ~u)) + FT i + CT i
(eq) (neq)
(4.12)
= fi (ρ, ~u) + (1 − ω) fi (~r, t) + FT i + CT i .
By analogy to this equation, the regularized model can be summarized as follows:
(eq) (1)
fi (~r + ~ci , t + 1) = fi (ρ, ~u) + (1 − ω) f¯i (~r, t) + FT i + CT i
(eq) (1 − ω) ti (4.13)
= fi (ρ, ~u) + Qi : Π(1) (~r, t) + FT i + CT i .
2 c4s
4.2.4 Related models

It is emphasized that the regularization procedure described in this section princi-
pally applies to all LB models of the BGK family. The idea is throughout the same:
(1) 2
the first-order term fi is projected on the vectors Cαβ corresponding to the first
(1)
non-conserved moment, and the resulting regularized term f¯i is directly computed
from the value of this non-conserved moment. For the purpose of illustration, the
advection-diffusion model described in Section 2.2 is now regularized. When the er-
ror and correction terms of this model cancel, for example by means of the procedure
described in Section 3.1, the first-order contribution to the model reads
(1) λ −
→ →
−
fi =− 2
ti Qi : ∇1 (ρ~u) − λ ti ~ci · ∇1 ρ. (4.14)
cs
The first non-conserved moment of this model is the momentum ~j, and the regular-
ized version of Eq. (4.14) takes therefore the following form:
(1) →
−
f¯i = −λ ti ~ci · ∇1 ρ. (4.15)
Indeed, the momentum computed from those two expressions is the same:
~j = −λ c2s −
→
∇1 ρ, (4.16)
(1)
from which the regularized term f¯i is computed as follows:
(1) ti
f¯i = 2 ~ci · ~j. (4.17)
cs
The LB model introduced here under the name of RLB has previously been
described in the literature, where it was derived from different principles and applied
in different 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 flows. 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 different from a typical
LB code and may be more difficult to work with. Furthermore, it depends on the
particular software and hardware platform whether this alternative formulation
of RLB leads to performance benefits or not.
4.3 RLB and MRT

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

fi (~r + ~ci , t + 1) = fi (ρ, ~u) + (1 − ω) 4
Q : ~
c j ~
c j fj − fj , (4.18)
j
2 c s
or, to emphasize the analogy with a general collision operator as in Eq. (3.27) on
page 39: !
X1−ω eq
out

fi = Qi : ~
c j ~
c j − δij fj − fj . (4.19)
j
2 c4s
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:

1 neq 1 3 ~
Ωi = −ti ω 4 Π : Q + hi N + ~ci · J . (4.20)
2 cs 4 8
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 verification

We now turn to numerical verifications of the regularized model on two 2D flows
using a D2Q9 lattice. The first test case, known under the name of Kovasznay flow,
approximates a stationary 2D flow behind a regular grid. An analytical solution for
this flow, proposed in [34], is written as follows in dimensionless units:
u∗x = 1 − exp(λ x∗ ) · cos(2πy ∗ ) and

λ (4.21)
u∗y = exp(λ x∗ ) · sin(2πy ∗ ),
2π
−3
10
Slope −2
Relative l2 norm error ε

−4
10 Slope −3
10
−5 BGK, Inamuro: slope −2.42
BGK, Skordos: slope −2.11
RLB, Inamuro: slope −2.94
RLB, Skordos: slope −2.87
10 20 30 40 50 60
Grid resolution N
Figure 4.1: Relative error of the numerical result on a Kovasznay flow in the wake
region. Both traditional BGK and the regularized model are tested with two com-
monly used boundary conditions.
with p
λ = Re/2 − 4π 2 + Re2 /4. (4.22)
The Reynolds number Re = u∞ L/ν is defined 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 fluid. 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 fixed 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 fixing the Mach number
M a = u/cs . This value is chosen small enough to mimic an incompressible flow.
Given that the flow 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 third-
order 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 deficiencies.
The second test case implements a flow in a 2D square cavity whose top-wall
moves with a uniform velocity. Reference values for the definition of the Reynolds
number are defined by the side length L of the box and the top-wall velocity u~0 .
Both standard BGK and the regularized model are first 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, defined 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.71·10−3 , and for the regularized method, of = 2.40·10−3 . 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 fixed 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 finite 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.
1000
Regularized: Remax = 16.3+7.36*N

Maximal Reynolds number Remax
800
600
400
200
BGK: Remax = 0.391+0.955*N
100
0
50 60 70 80 90 100 110 120 130
Grid resolution N
Figure 4.2: Simulation of 2D cavity flow for fixed Mach number. ◦, ∗: maximal
stable Reynolds number, numerically determined; solid line: least-square linear fit
of the data points (parameters of the fit are indicated on the graph).
(a) (b)
1 0.2
0.1
u*y on horizontal line
u*x on vertical line
0.5
0
−0.1
0
−0.2
−0.5 −0.3
0 32 64 96 128 0 32 64 96 128
Position [lattice units] Position [lattice units]
Figure 4.3: Comparison of computed RLB results with reference values from the
literature for the 2D lid-driven cavity flow 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 specified. The problems of fluid dynamics
are therefore initial and boundary value problems. It is common to specify the initial
value for a problem by imposing the velocity field at the initial time t = 0. Several
approaches exist for specifying the boundaries of a flow. It is possible to implement
a so-called Dirichlet boundary condition by specifying a velocity profile 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 specified,
or the value of the pressure on the boundary.
Whatever approach is adopted, one major difficulty 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 specifies 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 defined 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 profile is defined through an analytical expression, or otherwise
numerically, by applying for example a finite difference scheme. It is actually shown
in Section 5.3 that it is more difficult to properly initialize the pressure field 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 profile along the boundary is known
in terms of an analytical formula. Indeed, the gradients normal to the boundary
54 Chapter 5
. .
~c1 ~c8 ~c7
~c2 ~c6 Boundary
~e1
~c3 ~c4 ~c5
~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 finite difference 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 Boundary condition

5.2.1 Overview
This section focuses on the implementation of Dirichlet boundary conditions for the
velocity field. 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 fluid, but infinitesimally 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 55
5.2.2 Implementation of the pressure

The first step consists in computing the variable ρ on the wall. In an incompressible
fluid, 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 definition of three separate contributions to the density:
X
ρ1 = fk , (5.1a)
{k | ck1 =1}
X
ρ0 = fk , and (5.1b)
{k | ck1 =0}
X
ρ−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 finds that
ρ = ρ1 + ρ0 + ρ−1 and ρu1 = ρ1 − ρ−1 , (5.2)
and ρ can be formulated in terms of known quantities:

1
ρ= (2ρ−1 + ρ0 ). (5.3)
1 − u1
Note 5.1: Pressure boundaries Although only Dirichlet velocity boundaries

are discussed in this text, other boundary types are defined in a very similar
way. Equation (5.3) can for example be inverted to compute the normal veloc-
ity component u1 from the density ρ. This leads to the definition of a so-called
pressure boundary, on which the density ρ and the tangential velocity u0 are
specified, 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 finite difference 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 fitted according to
their zeroth- and first-order terms of the Chapman-Enskog expansion. It is however
not typical for simple LB methods to make use of non-local finite difference schemes,
as one prefers to implement a collision step depending only on values defined 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 “bounce-
back of off-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
off-equilibrium part for two distribution functions f¯i 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 dy-
namics. They are however not yet fully satisfying, as they lead to slightly
erroneous values for the density ρ and the velocity ~u. For this P reason, they are
(1) 2
only used to compute the value of the stress tensor Π = i fi − cs ρ I − ρ ~u~u
on the wall.
Regularization. Replace all distribution functions on the wall by their value ob-
tained 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 Chapman-
Enskog 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.
P
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 off-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 reflected 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 define 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
are specified 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 2 1
ρ + u0 + , (5.4a)
3 3

(1) 2 1
Π11 = 2 (f1 + f7 + f8 ) − ρ u1 + u1 + , (5.4b)
3
u
(1) 0
Π01 = 2 (f7 − f1 ) − ρ + u0 u1 . (5.4c)
3
The value of Π(1) on a D3Q19 lattice is:
(1)
Π00 = f1 + f2 + f3 + f5 + f6 + f7 + f8 +

u1 2 1
2(f10 + f12 ) − ρ + u0 + , (5.5a)
3 3
(1)
Π11 = 2(f9 + f10 + f11 + f12 + f13 ) −

2 1
ρ u1 + u1 + , (5.5b)
3
(1)
Π22 = f2 + f4 + f5 + f6 + f7 + f8 + 2(f11 + f13 ) −

u1 2 1
ρ + u2 + , (5.5c)
3 3
u
(1) 0
Π01 = 2(f10 − f12 ) − ρ + u1 u0 , (5.5d)
u3
(1) 2
Π12 = 2(f11 − f13 ) − ρ + u1 u2 , (5.5e)
3
(1)
Π02 = f5 − f6 + f7 − f8 − ρ u0 u2 . (5.5f)
5.2.4 Numerical verification

The quality of the regularized boundary condition is tested on a numerical imple-
mentation of the 2D Kovasznay flow described in Section 4.4. Four different types
of boundary conditions are implemented for this flow:
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 finite difference 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 fixed to a value Re = 20, and the grid
resolution is progressively refined to observe how the accuracy scales with the dis-
crete time step δx . Taking into account the observations of Section 2.4, the time step
58 Chapter 5
−2
10
−3
10
Relative error
−4
10
−5 Regularized local
10
Regularized non−local
Zou−He
−6
Inamuro
10 1 2
10 10
Grid resolution N
Figure 5.2: Accuracy of a 2D Kovasznay flow, implemented with four different

boundary conditions. The curves corresponding to the Zou-He and the Inamuro
methods are incomplete, because those methods produce numerical instabilities on
coarse grids.
δt is rescaled so as to keep constant the ratio δt /δx2 = 0.42. According to Note 2.7
on page 29, this amounts to fixing 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 flow, 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 condi-
tions 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 confirms 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 dy-
namics. Although the convergence rate is the same, the Zou-He and the Inamuro
boundary conditions are somewhat more accurate than the regularized one. This
difference can be understood by the fact that the benchmark simulations implement
a BGK dynamics, whereas the regularized boundary condition cuts off all higher
order modes of the dynamics. This difference 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.
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 Initial condition

5.3.1 Introduction
The impact of initial conditions on the evolution of a time-dependent flow 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 in-
teraction 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 difficulty
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 flows due to
constraints on the discrete time step.
In fluid turbulence, the interaction between vortexes and no-slip walls is known to
have an important influence on the evolution of a flow. Numerous numerical studies
have shown that the wall acts like a source for small-scale structures, which on their
hand interact with the fluid and affect 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 flow is incompressible, and
it is confined 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 fluid, described by the Navier-Stokes
equation in Eq. (1.10) on page 8, together with a continuity condition in Eq.(1.9).
The flow is confined in a box of dimension [−1, 1] × [−1, 1] with no-slip walls. The
initial condition is chosen so as to define 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 configuration. 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 field ~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
v0 = + |ωe | (x − x1 ) exp −(r1 /r0 ) − |ωe | (x − x2 ) exp −(r2 /r0 )2 ,(5.7)
2

2 2
p
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 defined by the
expression
1 1 1 2
Z Z
Ē(t) = |~u| (~x, t)d2 x, (5.8)
2 −1 −1
and the average enstrophy by
Z 1 Z 1
1
Ω̄(t) = ω 2 (~x, t)d2 x, (5.9)
2 −1 −1
where ω = ∂x v − ∂y u is the flow vorticity.

In all presented simulations, the initial kinetic energy is fixed 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 field for the initial configuration is shown on Fig. 5.3. The upper
monopole has a positive, and the lower one a negative core vorticity.
The flow is simulated by means of the RLB method, implemented on a D2Q9
lattice.
(a) (b)
2 2000
1800
1600
1400
1.5
kinetic energy
1200
enstrophy
1000
800
1 600
400
200
0
0.5 0 0.5 1 1.5 2
0 0.5 1 1.5 2 t
t
Figure 5.4: Time evolution of the kinetic energy (a) and the enstrophy (b)
5.4 Numerical results

5.4.1 Time evolution of the flow
Figure 5.4 displays the time evolution of the energy and the enstrophy of the system,
with the numerical parameters defined in the previous section. The first 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 finite difference 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 nu-
merical 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 sufficiently small grid spacing depends on the width of the bound-
ary 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.
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. Ad-
ditionally, 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 first
enstrophy peak with the LB, FD and spectral methods:
Lattice Boltzmann Finite Difference Spectral Method
Position (time) 0.371 0.371 0.371
Value (enstrophy) 933.8 932.8 933.6
5.5 Setup of the initial condition

The initial value for the velocity field 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 fields 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 flow is considered in its incompressible limit, the pressure p and the
density ρ are identified 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 first 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 fluid
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 first-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 fluid dynamics,
consists in computing the pressure p and the strain rate S directly. The velocity
gradients are e.g. computed using a finite difference 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 efficiency is compared in Fig. 5.6 (a). In this compar-
ison, 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
efficiency 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 fine-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
0 10
10
Relative L2 error of the pressure
Relative error on peak enstrophy

−1
10
−2
10
−3
10
−3
10
−4
10
BGK convergence
−5
Overrelaxed Gauss−Seidel
10 −4
0 2000 4000 6000 8000 10000 10 −6 −5 −4
Iteration steps 10 10 10
Time step ∆ t
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 first 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 influence is clearly visible in the present benchmark application.
In order for the LB fluid to mimic an incompressible fluid, 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 first enstrophy-peak
decreases as a function of the discrete time step. To achieve the accuracy of three
significant digits, a time step as small as δt = 10−6 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 · 10−5 . One is therefore tempted to conclude that the error plotted
in Fig. 5.6 (b) is dominated by contributions due to fluid 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 fixed 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 asymptot-
ically 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 difference methods, presented e.g. in
Ref. [4], are principally executed on regular grids with fixed spacing as well. The grid
can however be anisotropic, that is, the value of the grid spacing can differ 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 finite 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 inhomo-
geneous 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 com-
plicated shape to ensure that the discretized version of the obstacle resembles the
original one sufficiently well. In other cases, the resolution needs to be increased
because the fluid flow exhibits small-scale patterns, such as the small vortexes gen-
erated close to the wall in the numerical experiment of Section 5.3. In simulations
with fixed-sized grids, the overall resolution needs to be adapted to the maximal
required value, by which much computational effort 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 trade-
off between accuracy and efficiency. 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, de-
pending on how rapidly the flow 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 fixed grid spacing in the LB method

and to adapt the simulation at least roughly to the geometry of a problem, a so-
called multi block technique is often used. The idea is to partition the domain of
interest into rectangular subdomains, for each of which a different 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 difficulties need
to be overcome in order to do so. First, the data on the interface between two grids
needs to be interpolated, because the different 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 fluid with infinite extent. A system of successively refined,
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 flow. 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 fluid model.
Adaptive space and time steps 67
6.2 Grid refinement

6.2.1 Introduction
The multi block grid refinement technique of the regularized LB model discussed in
the previous section is now tested on a numerically challenging task: simulating a
stationary incompressible fluid flow past a rigid obstacle. The fluid is assumed to
fill the whole 2D space, the obstacle is placed at the origin of the coordinate system,
and the fluid velocity is asymptotically constant far from the obstacle. Thus, the
problem consists in the resolution of the stationary Navier-Stokes equation for the
velocity field ~u(~r) with boundary condition
lim ~u(~r) = ~u0 . (6.1)
|~
r|→∞
Apart from the grid refinement procedure, a major difficulty stems from the neces-
sity to truncate the infinite domain for numerical purposes, and to find an artificial
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 in-
appropriate 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 flow is strongly influenced 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 fluid and can there-
fore not be used on all boundaries. Furthermore, they make it difficult 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
field is proposed that can be used to implement a Dirichlet boundary condition for
the fluid velocity in a region reasonably far from the domain center. The expression
for this vector field is obtained from a truncated asymptotic expansion of a solution
to the stationary Navier-Stokes equation and approximates the structure of the
flow with considerably higher precision than the constant approximation ~u = ~u0 .
The drawback of this method is that it depends on the drag and lift coefficients 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 coefficients 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 refined 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 fluid with sufficient accuracy. The follow-
ing sections present a case study for the numerical evaluation of a drag coefficient,
and serve three main purposes. First, they present an introduction to the boundary
condition of Refs. [42, 43, 44] and demonstrate its efficiency 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.01 0.008
0.006
0.005 0.004
0.002
Figure 6.1: (a) Flow profile close to the obstacle when a zeroth-order boundary
condition ~u = ~u0 is used. (b) Velocity profile with the first-order boundary condition
introduced in this section.
with a numerical technique based on iterative grid refinement. The grid refinement
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 flows, it also proves useful for simulations of compressible flows at
low Mach numbers. For further precision, the influence of the fluid compressibility
on the computation of a drag force is analyzed.
6.2.2 Analytical profile on the boundaries

In Ref. [42], the solution to the 2D incompressible Navier-Stokes equation is ex-
panded in a finite series, as a function of formal parameters depending on the drag
and lift coefficients 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 flow doesn’t depend for the specific details of the obstacle ge-
ometry, but only on 2 (2D) or 3 (3D) distinct parameters. These considerations
result in the prescription of an explicit vector field that can be used as a boundary
condition on the numerical domain boundaries. The zeroth- and first-order terms of
the expansion for a 2D flow on the velocity field ~u = (u, v) at a position ~r = (x, y)
read
√ d 1 − y2

d x b y
u(x, y) = u∞ 1 + l 2 +l 2 − θ(x) l √ √ e 4lx ,
π x + y2 π x + y2 π x

√ d y − y2
(6.2)
d y b x
v(x, y) = u∞ l 2 −l 2 − θ(x) l √ 3/2 e 4lx
π x + y2 π x + y2 2 πx
where d = Fx /(2 ρ l u2∞ ) and b = Fy /(2 ρ l u2∞ ) are the drag and the lift coefficient,
and l = ν/u∞ is the viscous length, dependent on the dynamic fluid viscosity ν.
The formula also uses the Heaviside function θ(x), defined 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
(a) (b)
. .
Interface
. .
Figure 6.2: (a) Structure of the numerical grid close to a rectangular obstacle. One
dot on the figure represents a square of 8×8 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 profile close to the obstacle, depending on whether the
zeroth-order boundary condition ~u = ~u0 or the boundary condition defined 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 zero-
flux condition through the domain boundaries. The first-order boundary condition
on Fig. 6.1 (b) is free from such unphysical effects and produces, at least graphically,
the illusion of an infinitely 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 fluid is subject to sharp pressure and velocity variations. A grid
refinement technique is therefore applied with a hierarchy of nested grids that have
a successively finer 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
refinement, 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 effect of fluid 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 efficient choice for the simulation of incompressible fluids consists in
keeping the ratio δx2 /δ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 first 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
difference, it is remarked that Ref. [31] makes use of a first order accurate time
interpolation scheme. Although the accuracy of the time interpolation is not crucial
here, as a stationary flow 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 refinement 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 con-
vergence of the simulation towards a stationary state. In a first 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 implement-
ing a new, coarser grid. This is a convenient way of doing to efficiently converge
towards a stationary state. Another approach to reach a rapid convergence would
consist in refining 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 coefficients 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 bounce-
back 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
reflected. The exchanged momentum amounts to twice the value of the momentum
carried by the reflected distribution functions.
6.2.4 Simulation results

For illustration purposes, the simulation results of a 2D flow 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, defined with respect to the√height A of the obstacle, is fixed 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 coefficient 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 figure shows that the accuracy
(a) (b)
1 0
10 10
0
10
Relative error of drag
Compressibility error
−1
−1
10
10
−2
10
−2
10
−3
10
−4 −3
10 0 1 2 3 4 10
10 10 10 10 10 0.02 0.04 0.08
Diameter of computational domain δt / δx
Figure 6.3: (a) Drag coefficient as a function of the system size. The x-axis repre-
sents the ratio between the system size and the height of the obstacle. Solid line:
asymptotic boundary condition ~u = ~u0 . Dotted line: first-order accurate boundary
condition described by Eq. (6.2). (b) Solid line: error on the drag coefficient as a
function of the Mach number. Dotted line: reference curve, representing a power
law with slope −2.
of the drag coefficient 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 coefficient,
independently of the system size.
As a conclusion, the combination of a nested grid technique with an appropriate
boundary condition results in an efficient numerical scheme. Indeed, the computa-
tion 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 effects
play an important role. This statement is quantified by a plot on Fig. 6.3 (b), which
shows the difference between the drag force exerted by a compressible fluid and a
reference solution for an incompressible fluid found in Ref. [44]. As it is expected,
this difference decreases at roughly a second order rate with respect to the Mach
number. It is however not possible to implement fluids 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 Adaptive time

6.3.1 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 repre-
sentation using the previous time interval δt1 = tn −tn−1 . 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
~u0 = δt2 /δx~u∗ . Combining those two expressions yields ~u0 = δt2 /δt1~u. An adaptive
time interval can be used only to simulate incompressible flows, 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 effects. In this section, an adaptive time interval is
used to find the solution to stationary flow 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 modification 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 flow in a lid-driven cavity introduced in Section 4.4 is used again as a
benchmark application to verify the efficiency 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 field 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 flows and LB

In the present section, as well as in Sec-
tion 5.2, incompressible stationary flow 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 efficient way of doing. Other iterative procedures
exist that find solutions to stationary flow problems substantially faster. In those
methods, the intermediate steps of the iteration do not necessarily represent a
physically meaningful flow, by which shortcuts can be taken to reach the desired
final result. This limitation of the LB method can be partially overcome with
(a) (b)
Figure 6.4: Convergence of a 2D lid-driven cavity flow towards the steady state,
with a lattice Boltzmann RLB method and a finite difference method. (a) Fixed
time step δt = 1/2 δx2 in FD, fixed time step δt = 5/2 δx2 in RLB. (b) Fixed time step
δt = 1/2 δx2 in FD; variable time step from δt = 25 δx2 to δt = 5/2 δx2 in RLB.
the help of the techniques described in this thesis. In Section 5.2, the problem is
first solved on a domain with limited extent, surrounding the region of interest.
This gives a first rough approximation to the solution, which in the following is
improved by progressively enhancing the extent of the simulation. In the present
section, a first 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 efficiency of the LB model, the flow problem is also
solved with the help of a finite difference 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 fixed-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 flow, 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 fixed, as on
Fig. 6.4 (a), both simulations converge at a comparable speed. The time step of the
FD method is fixed 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
δt = 25 δx2 . 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
iteration step, until the final value δt = 5/2 δx2 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 efforts 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 fluid dynamics is introduced in Section 1.2. The overview
presented there shows that many different models for the numerical analysis of fluid
flows exist. On the one hand there are solvers based on discrete representations of
the Navier-Stokes equation, most notably finite difference, finite element or finite
volume methods. On the other hand, a new category of solvers has emerged over
the past decades which describe the fluid 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 different solver is investigated. In this approach, a finite difference (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 flow, one technique can be more
efficient, 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 modified,
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
find 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
. .
fi−1,j
r v
b i,j rfi,j
j ui−1,j
b
p
b i,j
u
b i,j
fi−1,j−1
r v f
b i,j−1 r i,j−1
b
FD variables
r LB variables
i
. .
Figure 7.1: Choice of indexes for FD and LB variables on a chosen grid cell (i, j).
prepared, the benefits 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 sufficiently stable, the
method makes use of a staggered grid, that is, a grid on which all macroscopic
variables are defined on different 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 re-
gion
Ω = [0, l] × [0, h] ∈ R2 ,
on which a regular grid is defined. 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 different 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 defined 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
. .
FD LB
cu
b cu
b
c
bv c
bv bv
c
j=4 rs rs rs rs
cu bp bu
b bp bu bp b
cu
cv
b bv bv bv bcv
j=3 rs r r rs
c bp bu
bu bp bu p
b bcu
cv
b bv bv bv bcv j=2 rs r r rs
cu bp bu
b bp bu bp b
cu
cv
b cv
b cv
b j=1 rs rs rs rs
cu
b cu
b
j=0
i=0 i=1 i=2 i=3 i=4 i=0 i=1 i=2 i=3 i=4
b FD bulk node r LB bulk node
c FD boundary node
b rs 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 defined 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
(t) (t) (t) (t+1) (t+1) (t+1)
iterations u , v , p ⇒ u ,v ,p contains two distinct steps:
1. Resolution of the Poisson equation (Eq. (1.11) on page 8) to obtain the new
pressure field. 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 in-
tegration 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 field according to a finite

difference scheme.

(t) (t) (t+1)
It uses the pressure field at step t + 1: u , v , p ⇒ 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

TheS 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 field, 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.
. .
i = iint
FD LB
bv bv r
s bv r bcv r
bu bp bu bp bu bp bcu
Ω1 bv bv r
s bv r bcv r Ω2
p p
bu b bu b bu b
p bcu
bv bv r
s bv r bcv r
b FD bulk node r LB bulk node

. b
c FD boundary node rs 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 The coupling algorithm

7.5.1 Coupling the pressure
Before turning to the velocity field, the treatment of the pressure in the FD simula-
tion 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 field of an incompressible flow is not uniquely defined. It contains
a constant offset which can be chosen arbitrarily. In the FD method used here,
the value of the offset depends on how fast the numerical method converges, and it
varies from one time step to another. It is therefore difficult to couple the pressure
field of the FD and LB simulations together. To do so would require a calibration of
the pressure fields 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 flow can be directly deduced
from the value of the velocity field 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 field from the LB field 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 defined 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 LB LB

iint ,j = 1/2 uiint ,j + uiint ,j−1 (7.2)
viFD = 1/2 viLB + viLB

int +1,j int +1,j int ,j
(7.3)
A similar procedure is used to translate the values of the velocity from the FD
field to the LB field:
uLB FD FD
iint −1,j = 1/2(uiint −1,j + uiint −1,j+1 ) (7.4)
viLB
int −1,j
= 1/2(viFD
int −1,j
+ viFD
int i,j
) (7.5)
(7.6)
7.5.3 Coupling the velocity gradients

The regularization formula (4.9) refers to the tensor S which, in Eq. (1.6) on page 7,
is defined in terms of velocity gradients. Figure 7.2 shows that two of those gradients
can be approximated by a centered difference of half the mesh width:
∂y uLB FD FD
iint −1,j = 1/δx (uiint −1,j+1 − uiint −1,j ) (7.7)
∂x viLB
int −1,j
= 1/δx (viFD
int ,j
− viFD
int ,j−1
) (7.8)
One gradient needs to be calculated as a centered difference of mesh width, based

on interpolated values of the velocity:
1 FD
∂y viLB
int −1,j
= (v + viFD
int −1,j+1
− (viFD + viFD
int −1,j−1
)). (7.9)
4δx iint ,j+1 int ,j−1
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:
∂x uLB LB
iint −1,j = −∂y viint −1,j (7.10)
7.5.4 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 first, 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 field at time t.
(out)
2. All nodes, bulk and boundaries, perform the collision step: fk (~r, t) = (1 −
(in) (eq)
ω)fk (~r, t) + ωfk (~r, t).
(in) (out)
3. The bulk nodes perform the streaming step: fk (~r +~vk δt , t + δt ) = fk (~r, t),
for all ~r such that ~r + vk δt lies on a bulk node.
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 field. 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 flow. This
is a stationary flow in a channel of infinite length with no-slip boundaries. An
analytical solution for this flow 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 fluid velocity
is strictly horizontal and does not depend on the x position: v = 0; ∂x u = 0. In
the present example, the fluid is driven by a body force. The left and the right
borders implement periodic boundary conditions in order to simulate a channel of
infinite 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: qP
jmax 2
j=0 (usim (·, j) − uana (·, j))
= qP (7.11)
jmax 2
u
j=0 ana (·, j)
Three simulations have been conducted that can be compared to each other.
The first simulation implements a pure LB scheme with the BGK model, the second
simulation a pure FD model, and the third simulation is a FD–LB hybrid. In the
first 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
0
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 (10−15 ). 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 field of the
LB simulation is however distinct from the analytical prediction, due to the limited
82 Chapter 7
. .
(3)
Ω3 (FD) jmax =3
L=1
~u (2)
jmax = 50
Ω2 (LB) jmax = 46
y
(1)
. Ω1 (FD) jmax =3 x.
Figure 7.4: The computational grid for the simulation of a Poiseuille flow is par-
titioned 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 field as a function of the grid resolution.
It appears clearly that the coupling acts like a second-order boundary condition for
the velocity field. No conclusion can be taken concerning the implementation of the
pressure field, because the latter is constant in a Poiseuille flow.
7.7 Conclusion
In this chapter, a LB scheme for 2D incompressible fluid flows 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 specific 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 flow 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 field 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
Figure 7.5: Simulation of a body-force driven Poiseuille flow 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 flow.
0
10
LB with bounceback
FD
FD−LB coupled
−1
10
error
−2
10
−3
10
−4
10
0 20 40 60 80 100 120 140 160
(iteration step)*δ t
Figure 7.6: Error of a FD-LB hybrid Poiseuille flow simulation as a function of the
grid resolution. A log-log plot shows the coupling of the velocity field to be of second
order.
−1
10
−2
10
error
−3
10
−4
10
Numerical result
−5
Slope −2
10 0 1 2 3
10 10 10 10
grid resolution (jmax)
84 Chapter 7
work presented here to other methods in CFD, such as the finite element method.
One key issue in this procedure would be to relate the values from the irregular
grid of typical finite 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 finite element package, which encourages the use of those two tools on
common problems.
Chapter 8
A model for fluid turbulence based on

kinetic variables
8.1 Turbulence modeling

8.1.1 Overview
In all numerical experiments conducted in the previous chapters, the simulated flows
have a laminar structure. At all time steps, it is possible to describe the velocity
field exhaustively via a restricted number of streamlines. Those lines are defined
as integral curves to the flow field: on every point of the curve, the flow velocity is
tangential to the curve. In a laminar flow, the streamlines vary smoothly in space and
time, and they are often used to illustrate the flow visually. In high Reynolds number
flows, a turbulent regime is observed to take place in which such a description is not
useful any more. Turbulent flows are characterized by the formation of structures
at many length scales. A description that focuses on a specific length scale, such as
a flow 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 flows is
substantially higher in presence of turbulent flows than with laminar flows. In order
to obtain an accurate representation of the flow 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 effect of the large eddies in the flow but
do not take into account the small ones. The effect of the small eddies is however
not neglected. Instead, a new, modified dynamics is executed on the computational
grid, which in the following is called the coarse grid. This dynamics replaces the
effect of the small eddies, had the simulation been executed on a fully resolved fine
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 effective 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

effective 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 effective viscosity is related. It must be emphasized that all those
efforts are required only for the simulation of turbulent flows. When laminar flows
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 flows, the numerical results differ
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 sufficient 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 flow.
8.1.2 Large eddy simulations in 3D incompressible fluids

At high Reynold numbers, a fluid enters a turbulent regime in which it is unsteady
and non-periodic. Turbulent flows are characterized by the occurrence of eddies
whose size may vary over a large range. The larger eddies contain the main portion
of the flow 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 Re−3/4 .
In a numerical simulation of a turbulent flow, 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 sufficiently fine 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 fluid turbulence based on kinetic variables 87
For this reason, the emphasis in turbulent simulations is put on the macroscop-
ically 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.
~ = h~ui and
Along this line of thought the velocity field ~u is split into a mean part U
the small variations ~u0 known as turbulent fluctuations, such that
~ + ~u0 .
~u = U (8.1)
The mean part is, for example, a componentwise statistical, temporal or spatial
average. The operator hi is often called a filter.
The filter hi 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~:
∂ ~ ~ − → ~ 1 −
→ ~ +−→
U+ U·∇ U = − ∇P + ν ∇2 U ∇ · h~u0~u0 i , (8.2)
∂t ρ0
where P = hpi 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 σ(~u0 ) defined as
σ(~u0 ) = − h~u0~u0 i . (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 PDE’s. It
shows however good results only in two dimensions and for very fine 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
0 2
k= |~u | and (8.5)
2
ν
|∂α u0β + ∂β u0α |2 .

= (8.6)
2
The value of those additional variables is determined by two additional transport
equations. The Reynolds stresses are approximated through the following term:
2 k2
σ(~u0 ) ≈ R(S, k, ) ≡ k I − 2 cµ S. (8.7)
3
This expression for the Reynolds stresses is substituted into Eq. (8.2), and a renor-
malized viscosity ν 0 is introduced:
k2
ν 0 = ν + νT = ν + 2cµ , (8.8)

88 Chapter 8
as well as a renormalized pressure:

2
P 0 = P + PT = P + ρ0 k. (8.9)
3
This leads to the following averaged momentum equation:
∂ ~ ~ →− ~ 1 −→ →
−
U+ U·∇ U = − ∇P 0 + ∇ · (ν 0 S) . (8.10)
∂t ρ0
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 difficult to implement in practice. A major
difficulty stems from the treatment of the domain boundaries, along which boundary
layers are observed, which don’t 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 flows, 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 simula-
tion, 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 defined 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 + v~ci , 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:
9
fieq (~u(~r, t), ρ(~r, t)) = ρ ti (1 + 3 ~ci · ~u + Q : ~u~u). (8.13)
2 i
In this equation, the generic lattice constants have been replaced by their value on
the D3Q19 lattice. The velocity ~u, specified by Eq. (1.16), depends linearly on the
particle flow densities fi . Furthermore, the flow 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
2
velocities are projected on the second order base vectors Cαβ , defined by means of
the tensors Qiαβ specified in Eq. (1.19).
In order to formulate a turbulence model for the LB dynamics, a filter hi is
introduced which acts as a linear map on the distribution functions f :
h i : f (t) −→ F (t). (8.14)
The dynamics of the filtered population
DT∗ : F (~r, t) → F (~r, t + T ) (8.15)
is defined as follows, using the functions composition operator ◦:
DT∗ ◦ hi = hi ◦ 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 filter 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 filtered
2
equations are projected on the vectors Cαβ . The resulting quantities are defined 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 configuration of the lattice kinetic variables f , yields the strain rate tensor,
defined in Eq. (1.6) on page 7, at each space position:
S : f (t) −→ S(t). (8.17)
With this definition, a weaker formulation of Eq. (8.16) is obtained as follows:
S ◦ DT∗ ◦ h i = S ◦ h i ◦ DT , (8.18)
~ the velocity associated to F , (8.18) becomes
Calling U
~ (~r, t)].
hS[~u(~r, t)]i = S[U (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 first order in terms of local variables on every lattice
site.
Turbulence modeling consists in approximating DT∗ by a term depending only on
the filtered 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
ω∗:
~ (~r, t), hρi (~r, t))).
Fi (~r + v~ci , t + δt) − Fi (~r, t) = −ω ∗ (~r, t)(Fi (~r, t) − fieq (U (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 different 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 verification 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
off-equilibrium part of the tensor Π. The findings of Eq. (2.66) on page 25 are
reproduced by the following equation:
3ω neq X
Sαβ ≈ − Π ; Πneq
αβ [f ] = ciα ciβ (fi − fieq [f ]) . (8.21)
ρ αβ i
With assumption (8.20), equation 8.19 becomes
ω hΠneq [f ]i = ω ∗ Πneq [F ], (8.22)
which gives an estimate for the value of the renormalized relaxation parameter:

neq
ω∗ Παβ [f ]
= neq ∀ α, β ∈ 0, 1, 2. (8.23)
ω Παβ [F ]
The numerator of the right hand side is further expanded:

neq X X
Παβ [f ] = ciα ciβ (Fi − fieq [F ] − δfieq ) = Πneq
αβ [F ] − ciα ciβ δfieq , (8.24)
i i
where the correction δfieq to the equilibrium term has been defined as
δfieq ≡ hfieq [f ]i − fieq [F ] (8.25)

9
~ = 9 ρ ti Qi : (h~u0~u0 i)
~U

= ti Qi : h~u~ui − U (8.26)
2 2
9
= − ρ ti Qi : σ. (8.27)
2
P
Using
Pthe properties of the D3Q19 grid introduced in Section 2.1 ( i ti Qiαβ = 0
and i ti Qiαβ Qiγδ Tγδ = 19 (Tαβ + Tβα ) for any second order tensor T ), one writes
X 2
ti Qiαβ Qiγδ σγδ = σαβ , (8.28)
i
9
and thus, using the definition of Qiαβ ,

X X 1 X eq
ciα ciβ δfieq = Qiαβ δfieq + δαβ δfi = −ρ σαβ (8.29)
i i
3 i
This yields
Πneq neq

αβ [f ] = Παβ [F ] + ρ σαβ , (8.30)
Figure 8.1: 2D cut through the vorticity field of the homogeneous isotropic turbulent
flow.
and, using the definition of the tensor Π:
ω∗ ρ σαβ σαβ
= 1 + neq = 1 − 3ω , ∀αβ. (8.31)
ω Παβ [F ] Sαβ
In Eq. 8.31, the effective 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 effective relaxation time. This is only true for
homogeneous isotropic flows, 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 verified numerically with the data of a DNS. The simu-
lation represents a statistically homogeneous and isotropic flow. The system consists
of a cubic lattice of system size 128 × 128 × 128 with periodic boundary conditions.
The Reynolds number, defined 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 sufficient for obtaining a weakly turbulent flow. 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
field 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 field, all components of the force field 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 field is transformed into real space by a Fourier
92 Chapter 8
transform. The method described in Ref. [61] guarantees that the obtained force
field is divergenceless. In such a way, the simulation stays close to the limit of fluid
incompressibility at each iteration step.
Figure 8.1 shows a two-dimensional cut through the fluid vorticity during a simu-
lation. This picture is typical for the data observed in turbulent flows. 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 different 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 en-
ergy scales down along the universal range, and it is finally dissipated in the low
wavelengths, close to the Kolmogorov length. As it is seen on the figure, the inter-
mediate 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 finally cross-checked with the data of the
DNS. The filter 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 off-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 definite 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 finds at
least a partial verification in the numerical data.
Figure 8.2: Energy spectrum of the turbulent flow.
Figure 8.3: Spatial distribution of the effective relaxation time τ ∗ /τ based on non-
diagonal components of σ and Πneq . The marker (∗) labels the value of the effective
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 effective 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 effective 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 efforts 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 defines a new,
multi-functional array type. The library is specifically 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 offers 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 scientific 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 different threads. The parallelization
is transparent to the programmer, which means that he doesn’t 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 soft-
ware.
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 fluid dynamics. The package is mainly intended as a programming
support for researchers and engineers who simulate fluid flows 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 configure 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 scientific 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
flow field in a turbulent lattice Boltzmann simulation. Physica A, 362:6–10,
2006.
[2] J. Latt and B. Chopard. Lattice Boltzmann method with regularized non-
equilibrium distribution functions. Math. Comp. Sim., 72:165–168, 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:169–172, 2006.
[4] S. Orszag, H. Chen, S. Succi, J. Latt, and B. Chopard. Turbulence effects on

kinetic equations. Journal of Scientific Computing, 28:459–466, 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 Interna-
tional 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:1023–1039, 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 fluid using a repulsive force. Int. J. Mod. Phys. C, 14:1015–1026,
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:43–51, 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 flows 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. Hänel. Molekulare Gasdynamik. Springer, 2004.
[2] A. Quarteroni and A. Valli. Numerical approximation of partial differential

equations. Springer, 1997.
[3] http://www.mathematik.uni-dortmund.de/∼kuzmin/cfdintro/cfd.html.
[4] M. Griebel, T. Dornseifer, and T. Neunhoffer. Numerical simulation in fluid

dynamics. SIAM, 1998.
[5] T. J. Chung. Computational fluid dynamics. Cambridge University Press, 2002.
[6] O. C. Zienkiewicz, R. L. Taylor, and P. Nithiarasu. Finite element method for

fluid dynamics. Elsevier, 2005.
[7] S. Chen and G.D. Doolen. Lattice Boltzmann method for fluid flows. Annu.
Rev. Fluid Mech., 30:329–364, 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 flow computations with
the method of lattice Boltzmann equation. Progr. Aerosp. Sciences, 39:329–
367, 2003. http://research.nianet.org/∼luo/Reprints-luo/2003/YuDZ PAS-
2003.pdf.
[10] C. Cercignani. Theory and application of the Boltzmann equation. Scottisch

Academic Press Ltd., 1974.
[11] S. Succi. The lattice Boltzmann equation, for fluid 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] http://www.openlb.org.
[14] D. A. Wolf-Gladrow. Lattice-Gas Cellular Automata and Lattice Boltzmann

models: an introduction. Lecture Notes in Mathematics, 1725. Springer, 2000.
[15] D. d’Humières, M. Bouzidi, and P. Lallemand. Thirteen-velocity three-

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:325–342, 2002.
[18] Z. Guo, C. Zheng, and B. Shi. Discrete lattice effects 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 flows. J. Stat. Phys., 81:35–48,
1995.
[20] X. He and L.-S. Luo. Lattice Boltzmann model for the incompressible Navier-
Stokes equation. J. Stat. Phys., 88:927–944, 1997.
[21] P. J. Dellar. Incompressible limits of lattice Boltzmann equations using multiple

relaxation times. J. Comput. Phys., 190:351–370, 2003.
[22] P. J. Dellar. Bulk and shear viscosities in lattice Boltzmann equations. Phys.
Rev. E, 64:031203, 2001.
[23] D. d’Humières. Generalized lattice-Boltzmann equations. Prog. Astronaut.

Aeronaut., 159:450–458, 1992.
[24] P. Lallemand and L.-S. Luo. Theory of the lattice Boltzmann method: Dis-
persion, dissipation, isotropy, galilean invariance, and stability. Phys. Rev. E,
61:6546–6562, 2000.
[25] R. Benzi, S. Succi, and M. Vergassola. The lattice Boltzmann equation: Theory
and applications. Phys. Reports, 222:145–197, 1992.
[26] Y. Chen, H. Ohashi, and M. Akiyama. Prandtl number of lattice bhatnagar-

gross-krook fluid. Phys. Fluids, 7:2280–2282, 1995.
[27] H. Chen, C. Teixeira, and K. Molvig. Digital physics approach to computational

fluid 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:7999–8003, 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 refinement algorithm. Phys. Rev. E, page 066707, 2003.
[32] A. J. C. Ladd and R. Verberg. Lattice-Boltzmann simulations of particle-fluid

suspensions. J. Stat. Phys., 104:1191–1251, 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 flows.
Physica A, 362:125–131, 2006.
[34] L. Kovasznay. Laminar flow behind a two-dimensional grid. Proc. Cambridge

Philos. Soc., 44:58–62, 1948.
[35] T. Inamuro, M. Yoshino, and F. Ogino. A non-slip boundary condition for

lattice Boltzmann simulations. Phys. Fluids, 7:2928–2930, 1995.
[36] P. A. Skordos. Initial and boundary conditions for the lattice Boltzmann
method. Phys. Rev. E, 48:4823–4942, 1993.
[37] U. Ghia, N. Ghia, and C. T. Shin. High-Re solutions for incompressible flow
using the Navier-Stokes equations and a multigrid method. J. Comp. Phys.,
48:387–411, 1982.
[38] Q. Zou and X. He. On pressure and velocity boundary conditions for the lattice
Boltzmann BGK model. Phys. Fluids, 9:1591–1598, 1997.
[39] I. Halliday, L. Hammond, and C. Care. Enhanced closure scheme for lattice
Boltzmann equation hydrodynamics. J. Phys. A: Math. Gen., 35:L157–L166,
2002.
[40] H.J.H. Clercx and C.-H. Bruneau. The normal and oblique collision of a dipole
with a no-slip boundary. Comp. fluids, 35:245–279, 2006.
[41] R. Mei, L.-S. Luo, P. Lallemand, and D. D’Humières. Consistent initial condi-
tions for lattice Boltzmann simulations. Comp. fluids, 35:855–862, 2006.
[42] S. Bönisch, V. Heuveline, and P. Wittwer. Adaptive boundary conditions for

exterior flow problems. J. Math. Fluid Mech., 7:87–107, 2005.
[43] S. Bönisch, V. Heuveline, and P. Wittwer. Leading order down-stream asymp-

totics of stationary navier-stokes flows in three dimensions. J. Math. Fluid
Mech., 7:147–186, 2005.
[44] P. Wittwer. Second order adaptive boundary conditions for exterior flow prob-
lems: non-symmetric stationary flows in two dimensions. J. Math. Fluid Mech.,
8:1–26, 2006.
[45] X. He, Q. Zou, L. S. Luo, and M. Dembo. Some progress in the lattice Boltz-
mann method. part I, non-uniform mesh grids. J. Comput. Phys., 129:357–363,
1996.
[46] O. Filippova and D. Hänel. Grid refinement for lattice-bgk models. J. Comput.
Phys., 147:219–228, 1998.
[47] J. Tölke, M. Krafczyk, M. Shulz, and E. Rank. Implicit discretization and

nonuniform mesh refinement approaches for FD discretization of LBGK models.
Int. J. Mod. Phys. C, 9:1143–1157, 1998.
102 Bibliography
[48] D. Kandhai, W. Soll, S. Chen, A. Hoekstra, and P. Sloot. Finite-difference

lattice-BGK methods on nested grids. Comput. Phys. Commun., 129:100–109,
2000.
[49] M. Bouzidi, D. d’Humières, P. Lallemand, and L. S. Luo. Lattice Boltzmann

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:835–845, 2003.
[51] A. Kolmogorov. The equations of turbulent motion in an incompressible fluid.

Izv. Sci. USSR Phys., 6:56–58, 1942.
[52] J. Smagorinsky. General circulation model of the athmosphere. Mon. Wheather

Rev., 91:99–164, 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:1159–1175, 1998.
[56] S. Pope. Turbulent flows. Cambridge University Press, 2000.
[57] S. Succi, O. Filippova, H. Chen, and S. Orszag. Towards a renormalized lattice

Boltzmann equation for fluid turbulence. Journ. Stat. Phys., 107:261–278, 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:379–394, 2004.
[59] S. Chen, S. Hou, J. Sterling, and G. D. Doolen. A lattice subgrid model for
high reynolds number flows. Fields Inst. Comm., 6:151–166, 1996.
[60] H. Chen, S. Succi, and S. Orszag. Analysis of subgrid scale turbulence using
the Boltzmann BGK kinetic equation. Phys. Rev. E, 59:2527–2530, 1999.
[61] K. Alvelius. Random forcing of three-dimensional homogeneous turbulence.

Phys. Fluids, 11:1880–1889, 1999.
[62] http://cui.unige.ch/∼latt/vladymir/index.html.
[63] http://spc.unige.ch/mptl.

Thesis PDF

Uploaded by

Copyright:

Available Formats

Thesis PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Thesis PDF

Uploaded by

Copyright:

Available Formats

Thesis

Hydrodynamic limit of lattice Boltzmann equations

Hydrodynamic Limit of Lattice

Doctorat ès Sciences

Thèse de Monsieur Jonas LÄTT

Hydrodynamic Limit of Lattice Boltzmann Equations

Genève, le 23 avril 2007

The members of the comittee were the following:

• Prof. Bastien Chopard, Supervisor

• Prof. Michel Droz, Supervisor

• Prof. Sauro Succi, Istituto Applicazioni Calcolo, Rome

Résumé en Français vii

1 Computational fluid dynamics and lattice Boltzmann 5

2 Multiscale Chapman-Enskog analysis 15

3 Corrections to the BGK dynamics 31

4 Regularized lattice Boltzmann 43

5 Boundary and initial conditions 53

6 Adaptive space and time steps 65

7 Coupling with other tools of CFD 75

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

8 A model for fluid turbulence based on kinetic variables 85

Ce document contient un choix sélectif des sujets étudiés au cours de la thèse,

Simulations en dynamique des fluides

l’équation du moment. Une équation d’état additionnelle couple la température

Méthode de Boltzmann sur réseau

La méthode LB appliquée au domaine de la CFD

en CFD classique modèlisent souvent directement ces variables-ci, c’est-à-dire, ils

Comment traiter des échelles temporelles et spatiales vari-

Les problèmes liés à la conversion entre variables macroscopiques et fonctions de

Conditions aux bords Conditions initiales

Couplage avec d’autres méthodes Applications multi-grille

Méthode LB et développement multi-échelle

L’équation de Navier-Stokes incompressible

Elle exprime la conservation de la quantité de mouvement d’un fluide caractérisé par

Figure 2: Le réseau D2Q9 pour la simulation de fluides en deux dimensions.

Méthode de Boltzmann sur réseau

fk (~xij + ~ck , t + 1) − fk (~xij , t) = Ωk (~xij , t),

où le terme Ω désigne l’opérateur de collision entre particules. Dans le modèle le

où on a défini l’équilibre, une version discrète de la distribution de Maxwell-Boltzmann,

Les quantités macroscopiques peuvent être calculées à partir de fonctions de

où on a défini le tenseur symétrique

Au cours du développement multi-échelle, on montre que ce tenseur est proportionnel

Modèle de Boltzmann sur réseau régularisé

En suivant l’enchainement suivant, la fonction de régularisation est utilisée pour

1. Calcul de ρ, ~u et Πneq = k ~ck ⊗ ~ck (fk − fkeq ).

3. Exécution de la collision BGK sur les fonctions de distribution régularisées.

Afin d’apprécier les qualités de ce nouveau modèle, un flux bidimensionnel connu

Computational fluid dynamics

The present thesis investigates a methodology in the domain of computational fluid

Lattice Boltzmann method

The LB method applied to CFD

How to treat varying space and time scales?

Multiscale expansion of lattice Boltzmann

Content of the thesis

Computational fluid dynamics and lattice

1.1 Vector and tensor notation

Vector/tensor notation Index notation

1.2 Macroscopic equations

The Numerical method provides a way to find approximate solutions of the

This expansion, up to a second order accuracy with respect to , appears to be suffi-

F~ = F~ (1) (2.28)

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

∂t1~u(~r, t) = ~u(~r, t + 1) − ~u(~r, t) + O(2 ). (3.4)