UNIVERSITE MOULAY ISMAIL

FACULTE DES SCIENCES MEKNES

DEPARTEMENT DE PHYSIQUE

Mémoire du Projet de Module

Rhéologie et écoulement polyphasique

FILIÈRE DE MASTER
Mécanique des fluides et Énergétique

Sous le thème :

Évaluation des modèles hybrides de Boltzmann sur

réseau pour la convection naturelle laminaire et
turbulente

Réalisé par les étudiants :

M. Agrimeh Iliass Encadré par :

Mme Dakkaki Khadija
Pr. Assoudi Redouane
Mme El Amri Souha
Mme Gara Douae

Année universitaire 2024–2025

Table des matières

Liste des figures 4

1 Fondements théoriques 6
1.1 Modes de tranfert thermique : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Types de Convection : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3 Domaine d’application de la convection naturelle . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4 Nombres adimensionnels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.5 Cadre Théorique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.6 Fonction de distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.7 Moments de la fonction de distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.8 Développement de la fonction de distribution à l’équilibre . . . . . . . . . . . . . . . . . . . . 15
1.9 Equation de Boltzmann . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.10 Théorème H de Boltzmann . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.11 l’opérateur BGK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.12 Équation de Boltzmann-BGK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.13 Développement de Chapman-Enskog . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2 Méthode de Lattice Boltzmann (LBM) 20

2.1 historique sur la LBM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.2 Les automates cellulaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3 Les gaz sur réseau . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4 Discrétisation des vitesses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.5 Discrétisation des moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.5.1 Méthode de quadrature de Gauss-Hermite . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.5.2 Application de la méthode de Gauss-Hermite pour discritiser les moments . . . . . . . 25
2.6 Discrétisation de la fonction de distribution à l’équilibre . . . . . . . . . . . . . . . . . . . . . 25
2.7 Equation de Boltzmann sur Réseau . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.7.1 Modèle dynamique de Boltzmann sur réseau (BGK) : . . . . . . . . . . . . . . . . . . 26
2.7.2 Modèle thermique de Boltzmann sur réseau : . . . . . . . . . . . . . . . . . . . . . . . 26
2.7.3 Modèle TRT (Two-Relaxation-Time) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.7.4 Modèle de Boltzmann sur réseau régularisé d’ordre élevé . . . . . . . . . . . . . . . . . 28
2.8 Schémas de relaxation multiples : MRT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.9 Implémentation d’une force externe dans la LBM . . . . . . . . . . . . . . . . . . . . . . . . . 33

3 Problèmatique et les Conditions aux limites 35

3.1 Description du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2 Hypothèses simplificatrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.3 Description mathématique du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.4 Modélisation sous forme adimensionnelle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.4.1 Mise en œuvre pratique de la méthode à double population de Boltzmann sur réseau
(TLBM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.4.2 Conditions aux limites adaptées à (LBM) . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.5 Mise en oeuvre informatique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

2
4 Résultas et discussion 41
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.2 Résultats et discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.3 Validation du modèle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.4 Choix du nombre de Rayleigh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.5 Evolution du Nombre de Nusselt en fonction du temps de relaxation . . . . . . . . . . . . . . 43
4.5.1 Évolution du nombre de Nusselt moyen pour différents nombre de Rayleigh . . . . . . 44
4.6 Résultats par régime laminaires 104 ≤ Ra ≤ 108 . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.6.1 Le champ de température (Isothermes) . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.6.2 Les lignes de courant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.6.3 Profil de température adimensionnelle et de la vitesse horizontale . . . . . . . . . . . 46
4.7 Résultats par régime Transition vers la turbulence (6, 4.108 < Ra < 1010 ) . . . . . . . . . . . 47
4.8 Résultats par régime fortement turbulents (Ra ≥ 1010 ) . . . . . . . . . . . . . . . . . . . . . . 48
4.9 Bilan des performances des modèles LBM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

3
Table des figures

1.1 Direction de la gravité et gradient de température. . . . . . . . . . . . . . . . . . . . . . . . . 8

1.2 Exemple de convection forcée. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Une particule se deplace librement dans un tube . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.4 Mouvement des particules du gaz dans un tube. . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.5 Mouvement d’un atome de gaz monoatomique. . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.6 Mouvement d’une molécule de gaz diatomique. . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.1 Quelques modèles de l’automate cellulaire "Vie" qui restent stables de génération en génération. 22
2.2 Le schéma D1Q3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3 Le schéma D2Q9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.4 Modèle D2Q5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.5 Réseau D2Q5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.1 Configuration physique et conditions aux limites. . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.2 les distributions manquantes sur la frontière inférieure pourle modèle D2Q9 à gauche et le
modèle D2Q5 à droite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.3 Organigramme de calcul de TLBM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

4.1 Comparaison entre le modèle d’article et le modèle obtenue par la simulation. . . . . . . . . . 42

4.2 Variation des nombres de Nusselt moyens pour Ra = 109 et P r = 0, 71 : (a) BGK-1 ; (b)
TRT-2 ; (c) RLB-4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.3 Variation des nombres de Nusselt moyens pour Pr = 0,71 : a) Ra = 6,4 × 108 , b) Ra = 2 × 109 . 44
4.4 Résultats du champ thermique pour les modèles BGK-1 , TRT-2 et RLB-4 à Ra = 106 .. . . . 45
4.5 Résultats du lignes de courant pour les modèles BGK-1 , TRT-2 et RLB-4à Ra = 106 .. . . . . 45
4.6 Profils de température dans une cavité chauffée latéralement pour Ra = 106 et P r = 0, 71
dans les modeles BGK-1 TRT-2 RLB-4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.7 Profils de température dans une cavité chauffée latéralement pour Ra = 106 et P r = 0, 71
dans les modeles BGK-1 TRT-2 RLB-4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.8 Champs thermiques moyennés (a, b, c) et instantanés (d, e, f), lignes de courant moyennées
dans le temps (g, h, j) ( a, d, g) BGK-1 ; (b, e, h) TRT-2 ; c, f, j) RLB-4 . . . . . . . . . . . . 48
4.9 Champs thermiques moyennés (a, b, c, d) et instantanés (e, f, g, h), lignes de courant moyen-
nées dans le temps (i, j, k, l) ( a, e, i) BGK-1 ; (b, f, g) TRT-2 ; c, g, k) RLB-4 . . . . . . . . 49
4.10 Champs thermiques et d’écoulement instantanés et moyennés temporellement pour les modèles
TRT-2 et RLB-4 avec Ra = 1011 et Pr = 0,71. . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.11 Analyse thermique et dynamique de l’écoulement. . . . . . . . . . . . . . . . . . . . . . . . . . 51

4
Introduction générale

La convection naturelle est un phénomène fondamental en mécanique des fluides et en transfert thermique,
jouant un rôle crucial dans de nombreuses applications industrielles et naturelles, telles que les systèmes de
refroidissement, les écoulements atmosphériques et les processus océaniques.
Comprendre et modéliser ce phénomène nécessite une compréhension fine des régimes d’écoulement (la-
minaire ou turbulent), qui dépendent essentiellement du nombre de Rayleigh, un paramètre adimensionnel
caractérisant l’intensité de la force de flottabilité par rapport aux forces visqueuses et conductives. Plus ce
nombre est élevé, plus le système tend vers un régime turbulent, rendant la modélisation des écoulements et
des transferts de chaleur plus complexe.[1]
Les méthodes traditionnelles de simulation numérique, comme la résolution directe des équations de
Navier-Stokes, peuvent s’avérer complexes et coûteuses en termes de calcul, surtout pour des écoulements
turbulents. Dans ce contexte, la méthode de Boltzmann sur réseau (Lattice Boltzmann Method, LBM)
s’impose comme une approche numérique innovante et efficace pour la simulation des écoulements fluides.
En se basant sur la théorie cinétique des gaz, la LBM permet de simuler de manière robuste aussi bien les
régimes laminaires que turbulents.[2]
Dans cette étude, nous nous concentrons sur trois modèles LBM : le modèle BGK (Bhatnagar-Gross-
Krook), le modèle TRT (Two-Relaxation-Time), et le modèle RLB-4 (Régularisé d’Ordre Élevé). Ces modèles
sont comparés pour leur capacité à reproduire les caractéristiques de la convection naturelle, en mettant
l’accent sur leur stabilité, leur précision et leur efficacité numérique.[3]
Les objectifs de ce rapport sont les suivants :
— Évaluer les modèles LBM : Analyser les performances des modèles BGK, TRT et RLB-4 pour des
écoulements laminaires et turbulents, en fonction de paramètres clés tels que le nombre de Rayleigh
(Ra) et le nombre de Prandtl (Pr).
— Valider les résultats : Comparer les simulations numériques avec des données théoriques et expé-
rimentales pour vérifier la fiabilité des modèles.
— Explorer les applications : Mettre en lumière les avantages et les limites de chaque modèle dans
des configurations pratiques, comme les cavités chauffées différentiellement.
Ce rapport est structuré en plusieurs chapitres :
— Le premier chapitre introduit les concepts théoriques de base, notamment les modes de transfert
thermique, les nombres adimensionnels, et les fondements de la théorie cinétique des gaz.
— Le deuxième chapitre présente les aspects numériques de la LBM, incluant la discrétisation des
vitesses et des moments, ainsi que les modèles de collision.
— Le troisième chapitre décrit la problématique étudiée, les hypothèses simplificatrices et les condi-
tions aux limites.
— Enfin, le quatrième chapitre expose les résultats des simulations et une analyse comparative des
modèles, suivie d’une discussion approfondie.
À travers cette étude, nous cherchons à contribuer à une meilleure compréhension des modèles LBM et
à leur optimisation pour des applications complexes en convection naturelle, tout en fournissant des pistes
pour de futures recherches dans ce domaine.

5
Chapitre 1

Fondements théoriques

Introduction
Ce premier chapitre pose les bases essentielles pour la compréhension des phénomènes de transfert ther-
mique étudiés dans le cadre de ce mémoire. Nous y examinerons les principales modalités de transfert de
chaleur la conduction, la convection et le rayonnement avec un focus particulier sur la convection naturelle,
un mécanisme fondamental dans de nombreux systèmes thermiques.
Une section sera dédiée à la classification des différents types de convection (forcée, naturelle, mixte),
en insistant sur les principes physiques qui les régissent et leur rôle dans les échanges thermiques. Cette
compréhension est primordiale pour l’analyse des écoulements thermoconvectifs présentée dans les chapitres
suivants.
Par ailleurs, nous introduirons les nombres adimensionnels incontournables en mécanique des fluides et
en transfert de chaleur. Les nombres de Rayleigh, Prandtl et Nusselt seront définis et interprétés, car ils
permettent de caractériser les régimes d’écoulement, d’évaluer les performances thermiques des systèmes, et
de faciliter la généralisation des résultats expérimentaux ou numériques.[4]
Enfin, une dernière partie de ce chapitre sera consacrée à la théorie cinétique des gaz, base fondamentale
pour le développement des modèles de réseau de Boltzmann. Nous y présenterons la fonction de distribution
des vitesses des particules, ainsi que ses moments, qui permettent de relier les descriptions microscopiques
aux grandeurs macroscopiques (densité, vitesse, énergie, etc.). Cette section constitue le socle théorique des
approches numériques utilisées dans ce travail, notamment la méthode Lattice Boltzmann.

1.1 Modes de tranfert thermique :

Le transfert thermique s’effectue spontanément dès qu’il existe une différence de température entre deux
points d’un même système ou entre deux systèmes différents, en l’absence de changement de phase. Il joue
un rôle crucial dans de nombreux domaines, notamment en mécanique des fluides et en génie thermique. Le
mécanisme de cette transmission de chaleur se produit selon trois modes

Conduction :
La conduction thermique est un mode de transfert thermique provoqué par une différence de température
entre deux milieux en contact, et se réalisant sans déplacement global de matière[5].
dT
Q = −k · A · (1.1)
dx
où :
Q flux thermique transféré par conduction (W),
k conductivité thermique du matériau (W m−1 K−1 ),
A surface à travers laquelle le transfert thermique a lieu (m2 ),
dT
gradient de température (K m−1 ), représentant la variation de température selon la direction x.
dx

6
Rayonnement thermique
Le rayonnement est une transmission d’énergie par le biais d’ondes électromagnétiques, sans déplacement
de matière et sans contact direct entre les éléments ou les milieux qui partagent l’énergie[6].

Q = ε · σ · A · Ts4 − Te4 (1.2)

où :
Q : la puissance thermique transférée par rayonnement (W),
ε : l’émissivité de la surface (sans unité),
σ : la constante de Stefan-Boltzmann (σ ≈ 5.67 × 10−8 W m−2 K−4 ),
A : la surface d’échange (m2 ),
Ts : la température de la surface émettrice (K),
Te : la température de l’environnement (K).

Convection thermique
La convection est un mécanisme de transfert de chaleur qui a lieu dans les fluides par le biais du dé-
placement des particules. Cette transmission découle du mélange thermique à l’intérieur du fluide et des
déplacements macroscopiques provoqués par les variations de température et de densité[6].

Q = h · A · (Ts − Tf ) (1.3)

où :
Q : puissance thermique transférée par convection (W),
2
h : coefficient de convection thermique (W/(m K)),
A : surface d’échange thermique (m2 ),
Ts : température de la surface (K),
Tf : température du fluide environnant (K).

1.2 Types de Convection :

Les transferts thermiques font partie des Sciences de base, peu de domaines industriels où la thermique
joue un rôle petit ou grand. Le chercheur, au minimum, doit donc connaître les principes fondamentaux, être
capable d’évaluer leur importance dans un problème particulier et dégager les ordres de grandeur par une
modélisation simple. Les transferts de chaleur qui s’effectuent simultanément avec des transferts de masse
sont des transferts de chaleur par convection[7].
Ce mode d’échange de chaleur existe au sein des milieux fluides dans lesquels il est généralement prépon-
dérant. La convection est un mode de transfert de chaleur où celui-ci est advecté (transporté, conduit) par
au moins un fluide.
Selon la nature du mécanisme qui provoque le mouvement du fluide on distingue :
— La convection libre ou naturelle
— La convection forcée
L’étude du transfert de chaleur par convection permet de déterminer les échanges de chaleur se produisant
entre un fluide et une paroi.

7
Convection naturelle ou libre
La convection naturelle se manifeste de manière spontanée. Elle a lieu dans un fluide ou un gradient de
température est présent. C’est le cas dans une salle où l’air chaud généré au niveau du sol par un radiateur
ou un convecteur va monter vers le plafond, tandis que l’air froid va redescendre. La poussée d’Archimède,
qui est due au fait que l’air chaud est moins dense que l’air froid, provoque le mouvement de l’air chaud vers
le haut[7].

N u = f (∆T, β, ρ) (1.4)

Figure 1.1 – Direction de la gravité et gradient de température.

Convection forcée
Ce second type de convection se manifeste lorsque le mouvement du fluide est forcé par une action
extérieure, sans lien avec la différence de température. Par exemple : une pompe, un ventilateur, un
agitateur ou même le vent. Dans le cas de la convection forcée, l’effet d’Archimède est insignifiant comparé
aux forces utilisées pour faire bouger le fluide[7].

N u = f (Re, P r) (1.5)

Figure 1.2 – Exemple de convection forcée.

Remarque : On identifie une troisième forme de convection : la convection mixte. Elle se manifeste
lorsqu’il existe un facteur externe au déplacement du fluide, mais ce dernier n’est pas suffisant pour que l’on
puisse ignorer la force d’Archimède.

8
1.3 Domaine d’application de la convection naturelle
Les applications de transfert thermique sont variées, dans lesquelles la convection naturelle est le phéno-
mène le plus dominant, la communauté scientifique autant qu’industrielle avaient un intérêt croissant aux
écoulements causés par la force de flottabilité, la meilleure compréhension de ce phénomène augmente le
nombre d’applications et mène à un certain nombre de conceptions industrielles et environnementales so-
phistiquées. Toutefois, les coûts de fonctionnement sont importants, les petites améliorations d’efficacité sont
essentielles et peuvent jouer un grand rôle dans la consommation d’énergie [8].
Le phénomène de la convection naturelle faire l’objet de maintes différente applications à savoir : les
problèmes océanographiques et atmosphériques tels que les effets de serre, les changements extrêmes de
climat, ainsi que les problèmes technologiques, à savoir les équipements électriques et les réacteurs nucléaires,
les capteurs solaires, le stockage des fluides, l’écoulement d’air dans les pièces d’habitation, les appareils
ménagers, les réfrigérateurs et les échangeurs de chaleur sont tous des problèmes, qui ont donné un intérêt
particulier à cette science [9].

1.4 Nombres adimensionnels

L’adimensionnalisation nous permet la mise en évidence de cinq paramètres de contrôle qui caracté-
risent le problème étudié, à savoir, respectivement, le nombre de Rayleigh Ra, le nombre de Prandtl
P r, et le nombre de Nusselt N u.

Nombre de Nusselt
Représente le rapport entre le transfert thermique total et le transfert thermique par conduction. Il s’agit
en fait du transfert de chaleur sans dimension. Si la conduction est le mode de transfert prédominant, alors le
nombre de Nusselt tendra vers 1. Si l’influence de la convection est importante, le nombre de Nusselt tendra
vers +∞. Le nombre de Nusselt s’écrit comme :
hL
Nu = (1.5)
λ
Avec h en [W·m−2 ·K−1 ] est le coefficient d’échange convectif. Ce coefficient d’échange de chaleur par convec-
tion dépend de l’écoulement du fluide, des propriétés thermiques Du milieu fluide et de la géométrie du
système[10]. [Grashof]. Au contact d’une paroi solide, la densité de flux de chaleur transmise au fluide est
donnée par :
q ′′ = h(Ts − Tf ) (1.6)
Avec q ′′ est la densité de flux de chaleur transmise au fluide ; Ts (K) est la température du fluide à la paroi ;
et Tf (K) est la température du fluide au loin de la paroi.

Nombre de Grashof
C’est un nombre sans dimension, utilisé en mécanique des fluides pour caractériser la convection naturelle
dans un fluide. Il correspond au rapport des forces de gravité sur les forces visqueuses. On le définit par :
∆T
Gr = βgL3 (1.7)
ν2
Avec L : La longueur caractéristique entre la paroi chaude et froide.

Nombre de Prandtl
Le nombre de Prandtl permet d’évaluer l’efficacité relative du transport de chaleur et de quantité de
mouvement en comparant les effets de diffusion de la chaleur et de quantité de mouvement, ce nombre est
défini par le rapport[10] :
ν λ
Pr = avec α=
α ρc
ν λ
Pr = avec α= (1.8)
α ρc

9
Nombre de Rayleigh
Dans la convection de Rayleigh-Bénard, la poussée résultant du gradient de température doit donc l’em-
porter sur la traînée viscous et la diffusion de la chaleur pour que la convection ait lieu.
Par conséquent, une relation entre ces trois paramètres s’exprime sous forme d’un rapport sans dimension :
la force de poussée divisée par le produit de la traînée visqueuse et du taux de diffusion de la chaleur est
appelée le nombre de Rayleigh (Ra) [10] :
 
Ts − Tf
Ra = βg L3 = Gr − P r (1.9)
ν

Où :
g : accélération due à la pesanteur
β : coefficient de dilatation thermique
Tc : température chaude
Tf : température froide
L : largeur
ν : viscosité cinématique
α : diffusivité thermique
Il est clair que le nombre de Rayleigh est un paramètre déterminant dans ce type de problème car il
exprime la force d’entraînement (flottabilité) dans la cavité, sans laquelle il n’y aurait aucune turbulence.
Un nombre de Rayleigh plus élevé engendre un effet de flottabilité plus grand et, par conséquent, plus de
turbulence. Un nombre de Rayleigh moins élevé signifie que la fluidité flottabilité s’applique sur l’écoulement,
le rendu laminaire. Plus le nombre de Rayleigh est élevé, plus la convection prend le pas sur la diffusion
de chaleur ou de quantité de mouvement. Ce nombre, une fois le fluide choisi pour une enceinte fermée, ne
dépend plus que de ∆T , qui est ainsi le paramètre de contrôle.

Influence du Nombre de Rayleigh sur les Régimes d’Écoulement

L’influence du nombre de Rayleigh (Ra) sur les régimes d’écoulement est un sujet clé dans le cadre de
la convection thermique. En effet, ce paramètre sans dimension détermine la transition entre les différents
types d’écoulement dans les fluides en mouvement sous l’effet d’un gradient de température.

Nombre de Rayleigh Régime d’Écoulement Description

(Ra)
Ra < 104 Laminaire L’écoulement reste laminaire. La convection ther-
mique est principalement due à la conduction, avec
des gradients de température relativement faibles.
104 < Ra < 106 Transition vers Turbu- La convection thermique commence à dominer.
lence Des cellules de Rayleigh-Bénard peuvent se former.
L’écoulement devient instable mais pas encore tur-
bulent.
Ra > 106 Turbulent L’écoulement est turbulent. Des mouvements chao-
tiques et irréguliers apparaissent. La convection ther-
mique est pleinement dominante et les structures de
convection sont complexes.

Table 1.1 – Classification des régimes d’écoulement en fonction du nombre de Rayleigh (Ra).

Le nombre de Rayleigh joue un rôle essentiel pour déterminer quel type d’écoulement prédomine dans un
fluide soumis à un gradient de température. Un faible nombre de Rayleigh entraîne un écoulement laminaire,
tandis qu’un nombre plus élevé peut aboutir à un écoulement turbulent.

10
1.5 Cadre Théorique
Théorie cinétique des gaz
Physiquement, un gaz est composé d’un grand ensemble d’atomes ou de molécules qui entrent en collision
les uns avec les autres et se déplacent de manière aléatoire. La description de la position d’une seule particule
(atome, molécule) peut être faite par trois coordonnées, mais le nombre de particules est très élevé, donc le
nombre total de variables est important. Par conséquent, l’étude de l’évolution des coordonnées de toutes
les particules, qui sont de l’ordre du nombre d’Avogadro, n’est pas envisageable.
La théorie cinétique des gaz, apparue au XVIIIe siècle, peut résoudre ce problème à l’aide de moyens
de description statistique [KineticTheory]. Cette théorie est une description des fluides qui est utilisée à
l’échelle microscopique, où on suit le mouvement des molécules individuelles, et à l’échelle macroscopique,
où on décrit le fluide en utilisant des lois statistiques. La théorie cinétique des gaz est principalement utilisée
pour décrire n’importe quel fluide, elle est la plus souvent appliquée au cas le plus simple d’un gaz dilué.
Dans ce cas, on peut supposer que les molécules constitutives passent très peu de temps à entrer en collision.
Cela équivaut à supposer que les molécules émergent presque toujours en collision en tête à tête, puisque les
particules ayant presque toujours impliqué simultanément une collision. Cette hypothèse ne tient pas aussi
bien pour les gaz denses ou les molécules plus proches les unes des autres et passant donc plus de temps à
entrer en collision, et elle ne s’applique pas du tout aux liquides où les molécules sont maintenues proches
les unes des autres par des forces d’attraction intermoléculaires et interagissent donc constamment[11] .

Pression cinétique
La cinétique du gaz idéal peut être comprise en déterminant la pression cinétique et l’énergie cinétique
moléculaire. Pour cela, nous supposons que les particules (molécules/atomes) du gaz sont des sphères rigides
dont le diamètre est négligeable par rapport aux distances intermoléculaires. Elles sont en agitation constante
et subissent des collisions élastiques sans interaction entre elles entre deux collisions. Les molécules suivent
des trajectoires rectilignes et sont réparties uniformément dans le récipient.
Considérons un gaz parfait constitué de N molécules dans un tube de longueur L, avec une particule de
masse m se déplaçant à la vitesse ⃗ci , entre deux parois séparées par une distance L. Lors des collisions, la
vitesse de la particule change de ⃗ci → ⃗cf . En effectuant le bilan des quantités de mouvement avant et après
la collision, nous obtenons :

Avant la collision :
P⃗i = P⃗p + P⃗m (1.11)
Après la collision :
P⃗f = P⃗p + P⃗m (1.12)
où :pi , pf , ppi et ppf sont les quantités de mouvement de la paroi et de la molécule, initiale et finale, respecti-
vement."

Figure 1.3 – Une particule se deplace librement dans un tube

11
 La conservation de la quantité de mouvement du système impose p⃗i = p⃗f . Pour un choc avec une particule,
la paroi reçoit la quantité de mouvement :

pp = p⃗f − p⃗i = p⃗f − p⃗m = p⃗f − p⃗m = 2mcx êx

∆⃗ (1.13)

Ainsi, la quantité de mouvement transférée à la paroi perpendiculaire à êx est ∆px = 2mcx à chaque collision
d’une molécule.
Le principe fondamental de la dynamique appliqué à la particule donne :

∆px mc2x
F = = (1.14)
∆t L
où F est la force exercée par la molécule sur la face perpendiculaire à x̂ et ∆t est le temps entre deux
collisions. L’intervalle de temps ∆t est égal à 2L/cx , représentant le temps nécessaire pour qu’une particule
se déplace d’une extrémité à l’autre et retourne au même endroit.
Pour obtenir la force totale, on somme toutes les molécules du gaz. Pour calculer la pression, il suffit de
diviser F par l’unité de surface :
N
1X 2
P = mcxi /L (1.15)
S i=1
et :
Nm 2
P = ⟨cx ⟩
V
où ⟨c2x ⟩ est la moyenne du carré des composantes des vitesses des molécules selon x̂ et V représente le volume
du tube (avec V = L × S).
En général, c2 = c2x + c2y + c2z , où cx , cy et cz sont les composantes de vitesses dans les directions x̂, ŷ
et ẑ, respectivement. Et comme les molécules ont des vitesses distribuées isotropiquement dans l’espace, par
conséquent ⟨c2x ⟩ = ⟨c2y ⟩ = ⟨c2z ⟩ = 13 c2 . Ainsi, l’expression finale de la pression cinétique est :

1 Nm 2
P = ⟨c ⟩ (1.16)
3 V
L’équation (1.8) montre que la pression exercée par le gaz sur les parois du tube est la manifestation ma-
croscopique des chocs des molécules sur les parois. Elle illustre l’esprit de la théorie cinétique des gaz en
mettant en relation une quantité macroscopique (pression) et une quantité microscopique (vitesse quadra-
tique moyenne).

Énergie cinétique moléculaire

Considérons un ensemble de particules d’un gaz parfait se déplaçant aléatoirement dans un tube (Figure
1.4). L’énergie possédée par le gaz est due uniquement aux mouvements des particules qui le composent
(énergie cinétique).
L’énergie cinétique moyenne (Ec ) d’une particule est :
1 2
Ec = mc (1.17)
2
Et l’énergie totale du système ET est donc :
1
ET = N ⟨Ec ⟩ = N mc2 (1,18)
2
On sait que l’équation d’état du gaz parfait s’écrit comme :

P V = nRT (1.19)

En introduisant la constante de Boltzmann (kB = NRA = 1.38 × 10−23 J/K) et le nombre de particules
(N = nNA ), l’équation d’état devient :
P V = N kB T (1.19)
Par identification avec l’équation précédente, il vient :

N mc2 = 3N kB T (1.20)

12
 Figure 1.4 – Mouvement des particules du gaz dans un tube.

On peut maintenant facilement déduire que l’énergie totale est :

3 3
ET = P V = N kB T (1.21)
2 2
Et l’énergie cinétique moyenne d’une particule est ainsi :
3
Ec = kB T (1.22)
2
ou encore :
3
Ec = kB T (1.23)
2
ou encore :
1 1 1 3
m(c2x ) + m(c2y ) + m(c2z ) = kB T (1.24)
2 2 2 2
It must be noted that the amount 21 kB T corresponds to the energy thermal of a particule and to a degree
of freedom. In effect, for an atome of gaz monoatomique, there are three degrees of freedom that
correspond to the three possible directions of travelvel (Figure 1.5).

Figure 1.5 – Mouvement d’un atome de gaz monoatomique.

Pour un gaz diatomique, il faut considérer deux degrés de liberté supplémentaires liés aux rotations de
la molécule (Figure 1.6). L’énergie associée à chaque molécule est alors :
5
Ec = kB T (1.24)
2

13
 Figure 1.6 – Mouvement d’une molécule de gaz diatomique.

On peut maintenant déduire que la température et la pression dans le monde macroscopique ne sont que
des échelles de l’énergie cinétique des molécules dans le monde microscopique. La section suivante traitera
d’un concept très intéressant dans le domaine de la physique statistique : il s’agit de la fonction de distribution.

1.6 Fonction de distribution

Définition
La fonction de distribution en physique statistique décrit la répartition des particules d’un système
en fonction d’une variable spécifique, telle que la vitesse. Elle permet de déterminer la probabilité qu’une
particule ait une vitesse dans un certain intervalle donné. Plus précisément, dans un système contenant un
grand nombre de particules, cette fonction permet de calculer la densité de particules ou la masse volumique
dans un espace donné, en fonction des propriétés du système étudié. Considérons un système constitué d’un
grand nombre N de particules [12].
On note F (⃗x, ⃗c, t) la fonction de distribution en nombre à une particule. Par définition, la grandeur

F (⃗x, ⃗c, t) d⃗x d⃗c (1.25)

représente la probabilité de trouver une particule ayant une vitesse comprise entre c et c + dc dans le
volume dx autour de x. Dans un espace de dimension D, la grandeur F s’exprime en L−D T −D . On a la
relation de fermeture :
Z
F (⃗x, ⃗c, t) d⃗x d⃗c = 1 (1.26)

La densité de particules, c’est-à-dire le nombre de particules par unité de volume, de dimension L−D , est
notée n(⃗x, t) et est définie par :
Z
n(⃗x, t) = F (⃗x, ⃗c, t) d⃗c (1.27)

On a évidemment :
Z
n(⃗x, t) d⃗x = 1 (1.28)

Si les particules ont toutes la même masse mp , alors on définit la fonction de distribution massique à une
particule f par :

f = N mp F (1.29)
En dimension D, la grandeur F s’exprime en M L−2 DT D . Autrement formulée, la grandeur :

f (⃗x, ⃗c, t).d⃗x.d⃗x (1.30)

représente la valeur espérée de la masse des particules ayant une vitesse comprise entre c et c + dc et se
trouvant dans le volume dx autour de x. Selon la dimension spatiale D du système étudié, on a :

 pour D = 1 : d⃗c ≡ dcx d⃗x ≡ dx



pour D = 2 : d⃗c ≡ dcx dcy d⃗x ≡ dx dy

pour D = 3 : d⃗c ≡ dcx dcy dcz d⃗x ≡ dx dy dz


14
1.7 Moments de la fonction de distribution
Dans le cadre de la théorie cinétique des gaz, la fonction de distribution f(x,c,t) représente la probabilité
de localiser une particule à la position x, avec une vitesse c, à un temps donné t. Les moments de cette
fonction facilitent le lien entre la description à l’échelle microscopique (niveau moléculaire) et les variables à
grande échelle utilisées dans la mécanique des [Link] moments de la fonction de distribution constituent
des intégrales pondérées de ladite fonction et facilitent le lien entre les quantités à l’échelle microscopique et
les grandeurs macroscopiques du fluide, telles que la densité, la vitesse moyenne et la pression[13].
Par exemple, la densité de masse macroscopique peut être déterminée en tant que moment lié à la
concentration de masse de toutes les particules à la position x et au temps t, c’est-à-dire l’agrégation de
toutes les fonctions de distribution à cette position et ce moment.
Z
ρ(x, t) = f (x, c, t)dc3 (1.31)

Avec, ρ représente la masse totale des particules présentes dans un petit volume autour du point x. On peut
aussi prendre en compte l’apport des particules f à la densité du moment macroscopique.
En révaluant toutes les vitesses possibles, on détermine que la densité du moment cinétique macroscopique
correspond à la somme vectorielle de la quantité de mouvement de l’ensemble des particules :
Z
ρ(x, t) ∗ u(x, t) = cf (x, c, t)dc3 , (1.32)

u(x,t) : vitesse moyenne du fluide en x.

On obtient la densité d’énergie cinétique macroscopique Ec en additionnant l’énergie cinétique de chaque
particule : Z
1
ρ(x, t) ∗ Ec (x, t) = |c|2 f (x, c, t) dc3 , (1.33)
2 R3
Ec (x, t) : énergie cinétique macroscopique par unité de masse.
e : énergie interne spécifique, qui représente l’énergie cinétique due aux fluctuations microscopiques des
particules autour de la vitesse moyenne u.

1.8 Développement de la fonction de distribution à l’équilibre

L’un des éléments clés pour l’application de la méthode LBM (Lattice Boltzmann Method) à divers pro-
blèmes physiques est la fonction de distribution à l’équilibre, notée f eq . En effet, différents phénomènes
physiques, tels que la propagation d’ondes ou le transfert de chaleur, peuvent être modélisés par la LBM, à
condition d’utiliser une expression appropriée de f eq et un temps de relaxation adapté [14].
La fonction de distribution à l’équilibre est généralement donnée par :

|⃗c − ⃗u|2
 
ρ
eq
f (⃗c) = √ exp − , (1.34)
D 2πRT 2RT

Cette expression exponentielle indique que les particules sont statistiquement concentrées autour de la
vitesse moyenne ⃗u du fluide.
Où :
T :température thermodynamique du fluide.
R :constante des gaz parfaits spécifique au fluide.
D :dimension spatiale du système.

Développement de l’expression de f eq
Nous cherchons à développer l’expression précédente de f eq en utilisant un développement limité. Pour
cela, nous nous intéressons à l’exposant de la fonction exponentielle :

(c − u)2
 
exp − , (1.35)
2RT
Développons l’expression de l’exposant :

15
 (c − u)2
 2
c + u2 − 2cu
  
2 2
(c − u) = c − 2cu + u 2
⇒ exp − = exp − , (1.36)
2RT 2RT
Grâce aux propriétés des exponentielles, cette expression peut être factorisée :

(c − u)2 c2
     2 
u − 2cu
exp − = exp − exp − , (1.37)
2RT 2RT 2RT
√
Sous l’hypothèse de faibles vitesses (i.e. u ≪ RT ), nous pouvons approximer le deuxième terme par un
développement limité d’ordre 2 :
 2
u2 (cu)2

u − 2cu cu
exp − ≈1+ − + (1.38)
2RT RT 2RT 2(RT )2
En combinant cette approximation avec le facteur gaussien initial, on obtient une forme développée de
la fonction de distribution à l’équilibre :

u2 (cu)2 c2
   
ρ cu
f eq (c) ≈ √ 1+ − + exp − (1.39)
2πRT RT 2RT 2(RT )2 2RT
Cette expression développée est particulièrement utile dans les implémentations de la méthode LBM
pour formuler les modèles simplifiés tels que BGK ou TRT, en facilitant le calcul explicite des moments et
l’ajustement des paramètres numériques.

1.9 Equation de Boltzmann

Bien que les particules d’un gaz puissent être considérées comme soumises aux équations de Newton, leur
nombre élevé rend impossible l’utilisation d’équations différentielles couplées pour suivre chacune d’entre
elles. Une description statistique simplifiée devient donc nécessaire. Cette description introduite il y a plus
d’un siècle par L. Boltzmann et J.C. Maxwell fait intervenir La fonction de distribution, ou densité, f =
f (t, x, c) du gaz constitué de particules ponctuelles repérées à l’instant t par leur position x ∈ RN et leur
vitesse c ∈ RN .
La fonction de distribution f (t, x, c) d’une particule simple donne la probabilité de trouver au temps t
une particule à une position donnée x et ayant une vitesse c. Le nombre probable de particules se trouvant
dans l’étendue x + dx et d’une vitesse entre c + dc est donné par f (t, x, c) dx dc.
Soit F une force externe (moyennement faible par rapport aux forces intermoléculaires). En incrémentant
le temps de t à t + dt, il existe des particules qui partent de (x, c) et arrivent à (x + dx, c + dc).
En raison du phénomène de collision durant le temps dt, il existe un nombre de particules qui partent de
(x, c) et n’arrivent pas à (x + dx, c + dc), soit Ω− dx dc, et un autre nombre de particules partant quelque
part autre que (x, c) mais y arrivent,soit Ω+ dxdc. En conséquence, la variation de la fonction de distribution
pendant le temps dt est régie par l’équation suivante [15] :

f (x + dx, c + dc, t + dt) dxdc − f (x, c, t) dxdc = (Ω+ − Ω− ) dx dc (1.40)

Nous appliquons d’abord un développement de Taylor au premier ordre sur le terme à gauche de cette
équation autour de f (x, c, t), ce qui nous donne :
   
∂f
f (x + dx, c + dc, t + dt) − f (x, c, t) = dx∇x f + dc∇c f + dt dx dc dt (1.41)
∂t
Ainsi, nous obtenons l’expression suivante :
   
∂f
dx∇x f + dc∇c f + dt dx dc dt = (Ω+ − Ω− ) dx dc dt (1.42)
∂t
En divisant ensuite cette équation par dt, nous arrivons à :
 
dc dx ∂f
∇x f + ∇c f + = Ω+ − Ω− (1.43)
dt dt ∂t
En réorganisant cette dernière expression, nous obtenons :
 
∂f
c∇x f + γ∇c f + = Ω+ − Ω− (1.44)
∂t

16
 dt représente la vitesse du système. En appliquant la loi de la dynamique de phase (PFD) et en
Ici, γ = dx
substituant la force F dans l’expression, nous arrivons à l’équation suivante :
 
F ∂f
c∇x f + ∇c f + = Ω+ − Ω− (1.45)
m ∂t
Cette équation est précisément l’équation de Boltzmann, où les termes à gauche décrivent la propagation
et l’effet des forces sur les particules, et le terme à droite représente les collisions entre les particules.

Absence de collision
Si aucune collision n’a lieu entre les particules, le nombre de molécules f (x, c, t) avant l’application de la
force extérieure est égal au nombre de molécules après la perturbation, f (x + cdt, c + m F
dt).
La dérivée de la fonction f par rapport au temps est donc nulle :
∂f ∂f F ∂f
+c + =0 (1.46)
∂t ∂x m ∂c
Cette équation représente l’équation de Boltzmann continue sans collision. Elle peut être généralisée comme
suit :

∂f F⃗
+ ⃗c · ∇x f + · ∇c f = 0 (1.1)
∂t m

Presence de collision
Si des collisions ont lieu entre les molécules, il y aura une différence nette entre les nombres de molécules
dans l’intervalle dx dc. Dans ce cas, le deuxième terme ne sera pas nul. L’équation de Boltzmann avec
collisions s’écrit alors [16] :
∂f F⃗
+ ⃗c · ∇⃗x f + · ∇⃗c f = Ω(f ) (1.47)
∂t m
où Ω(f ) = Ω+ − Ω− : représente le terme de collision, qui prend en compte les échanges de particules entre
les différents états de mouvement du système.

1.10 Théorème H de Boltzmann

Il est également possible de relier la fonction de distribution f à la propriété thermodynamique de
l’entropie. Par conséquent, l’entropie de la distribution de densité peut être exprimée sous la forme suivante
[17] :

s(t) = −kB H(t) + ϕ (1.49)

où :
— kB est la constante de Boltzmann ;
— ϕ est une constante arbitraire ;
— H(t) est une mesure introduite pour la première fois par Boltzmann, définie comme suit :
Z
H(t) = f (⃗r, ⃗c, t) ln f (⃗r, ⃗c, t) d⃗r d⃗c (1.50)

où :
— f (⃗r, ⃗c, t) est la fonction de distribution des particules en fonction de la position ⃗r, de la vitesse ⃗c, et
du temps t ;
— L’intégrale est effectuée sur tout l’espace des positions et des vitesses.
L’évolution de la quantité H par rapport au temps est :
Z Z
dH ∂f ∂ ∂f
= f ln f d⃗r d⃗c = (1 + ln f ) d⃗r d⃗c (1.51)
dt ∂t ∂f ∂t
En utilisant l’équation de transport de Boltzmann (Eq. (1.28)) pour évaluer ∂f /∂t, nous obtenons :
!
F⃗ ∂f
Z
dH ∂f
= (1 + ln f ) −⃗c · − · + Ω(f ) d⃗r d⃗c (1.52)
dt ∂⃗r m ∂⃗c

17
 Le développement mathématique de cette équation pour le cas des collisions binaires montre que la
quantité H est strictement décroissante dans le temps jusqu’à atteindre une valeur constante correspondant
à l’équilibre thermodynamique [18] :
dH
≤0 (1.53)
dt
par conséquent :
dS
≥0 (1.54)
dt
Ce théorème est connu sous le nom de théorème H de Boltzmann. Et d’après ce théorème, l’entropie
d’un gaz qui satisfait à l’équation de Boltzmann ne peut donc qu’augmenter avec le temps. Par conséquent,
l’équation de Boltzmann décrit la relaxation d’un gaz raréfié vers l’équilibre. Cela conduit à une forme simple
de l’opérateur de collision, comme nous le verrons par la suite, il s’agit de l’opérateur BGK.

1.11 l’opérateur BGK

Comme stipulé auparavant, la pertinence du théorème H découle de sa contribution à l’élaboration d’un
opérateur de collision très simple et distinct de celui suggéré par Boltzmann, qui se concentre uniquement
sur les collisions binaires. La manière la plus simple a été proposée par Bhatnagar, Gross et Krook en
1954 et porte le nom d’approximation BGK. Elle consiste à dire que chaque collision change la fonction
de distribution f d’une quantité proportionnelle à l’écart entre la distribution actuelle f et la distribution à
l’équilibre f eq [19]. On définit donc l’opérateur de collision BGK par :
L’opérateur de collision BGK est défini par :
1
ΩBGK = − (f − f eq ) (1.55)
τ
f est la fonction de distribution actuelle. f eq est la fonction de distribution d’équilibre obtenue à partir
de la distribution de Maxwell-Boltzmann. τ est le temps de relaxation, contrôlant la vitesse de retour à
1
l’équilibre. est la fréquence des collisions.
τ

1.12 Équation de Boltzmann-BGK

L’équation de Boltzmann intégrant l’opérateur BGK est :
∂f F ∂f 1
+ c · ∇f + · = − (f − f eq ) (1.55)
∂t m ∂c τ
Pour un système sans forces externes :
∂f 1
+ c · ∇f = − (f − f eq ) (1.56)
∂t τ
Le modèle suppose un retour progressif des particules vers l’équilibre à une vitesse contrôlée par τ .
Il permet de simplifier les calculs tout en maintenant les propriétés fondamentales des interactions par-
ticulaires. Il est une bonne approximation des équations de Navier-Stokes qui régissent les fluides.
La relation entre le temps de relaxation τ et la viscosité cinématique ν est donnée par :
 
1 2
ν= τ− c (1.57)
2 s

où
cs est la vitesse du son du modèle de réseau (pour D2Q9, c2s = 13 ).

18
1.13 Développement de Chapman-Enskog
Les équations macroscopiques obtenues dépendent alors des fonctions de distribution par l’intermédiaire
de Pij et qi . Il faut donc trouver un moyen d’exprimer ces grandeurs de façon explicite ne dépendant plus des
fonctions de distribution. Ce problème de fermeture des équations fut résolu par Chapman et Enskog dans
les années 1920. L’idée réside dans le développement des fonctions de distribution en fonction du nombre
de Knudsen. Ce nombre sans dimension, fondamental en physique statistique, traduit le rapport entre le
libre parcours moyen ℓ des particules, c’est-à-dire la distance sur laquelle ces particules ne subissent aucune
collision, et la longueur caractéristique L de l’écoulement [20] :
√
L λ0 rT0
ϵ= = (1.58)
L L
√ √
— rT0 : c’est la vitesse thermique moyenne. En effet, en physique des gaz, rT donne une échelle
typique des vitesses.
— L : échelle de longueur caractéristique (ex. taille du domaine, diamètre du tuyau. . .)
— λ0 : temps moyen entre deux collisions → permet de définir un temps caractéristique
En d’autres termes, l’analyse de Chapman–Enskog intervient en particulier pour déterminer la partie
de non-équilibre. Au cœur de cette analyse se trouve une expansion de perturbation de f autour de la
distribution d’équilibre feq (où feq = f (0) ), avec le nombre de Knudsen Kn comme paramètre d’expansion.
En utilisant l’étiquette εn pour indiquer les termes d’ordre Knn , l’expansion s’écrit :

f = f (0) + εf (1) + ε2 f (2) + O(ε3 ) (1.59)

le développement de Chapman-Enskog permet également de trouver l’expression mathématique du temps
de relaxation simple (τ ) du modèle BGK :

ν = ∆tc2s (τ − 0.5) (1.47)

√
En unité LBM, ∆t est égal à l’unité. Et pour le modèle D2Q9, la vitesse du son est fixée à 1/ 3, ainsi
l’équation (1.47) devient :
1
ν=(τ − 0.5) (1.60)
3
L’équation montre que τ > 1/2 est une condition obligatoire pour la stabilité du schéma BGK, puisque
τ < 1/2 provoquerait une situation macroscopiquement instable de viscosité négative.

Conclusion
Le premier chapitre a posé les bases théoriques essentielles pour comprendre les phénomènes de transfert
thermique, en mettant l’accent sur la convection naturelle. Nous avons exploré les différents modes de trans-
fert de chaleur (conduction, convection, rayonnement) et détaillé les mécanismes spécifiques à la convection
naturelle, forcée , Les nombres adimensionnels clés, tels que les nombres de Rayleigh, Prandtl et Nusselt,
ont été introduits pour caractériser les régimes d’écoulement et évaluer les performances thermiques des
systèmes.
La théorie cinétique des gaz a été présentée comme fondement des modèles de réseau de Boltzmann, avec
une attention particulière portée à la fonction de distribution et à ses moments. Ces concepts permettent
de relier les descriptions microscopiques aux grandeurs macroscopiques, telles que la densité, la vitesse et
l’énergie. L’équation de Boltzmann et ses simplifications, notamment via l’opérateur BGK, ont été discutées
pour montrer leur rôle dans la modélisation des écoulements fluides.
Enfin, le développement de Chapman-Enskog a été abordé pour expliquer comment les équations ma-
croscopiques, comme celles de Navier-Stokes, peuvent être dérivées à partir de l’équation de Boltzmann. Ce
cadre théorique est crucial pour la méthode de Boltzmann sur réseau (LBM), qui sera appliquée dans les
chapitres suivants pour simuler des écoulements laminaires et turbulents.
Ce chapitre a ainsi fourni les outils conceptuels nécessaires pour aborder les modèles numériques et les
simulations présentées dans la suite de l’étude, tout en soulignant l’importance des approches statistiques et
cinétiques dans la modélisation des systèmes thermo-fluidiques complexes.

19
Chapitre 2

Méthode de Lattice Boltzmann (LBM)

Introduction
Ce deuxième chapitre s’inscrit dans la continuité du cadre théorique établi précédemment. Après avoir
introduit la théorie cinétique des gaz et les principes fondamentaux de la fonction de distribution dans le
premier chapitre, nous nous intéressons ici aux aspects numériques de la méthode de Boltzmann sur réseau
(LBM). Nous commencerons par détailler la discrétisation des vitesses, des moments et des fonctions de
distribution , éléments centraux dans le passage du modèle cinétique continu au cadre du réseau. Nous
introduirons ensuite la fonction de distribution d’équilibre, qui joue un rôle fondamental dans les processus
de relaxation. La seconde partie du chapitre sera consacrée à la formulation discrète des modèles de collision,
essentiels dans la modélisation du comportement des fluides. Nous décrirons les trois modèles hybrides retenus
dans notre étude, le modèle BGK-1, TRT-2 et RLB-4.
Ces formulations discrètes poseront les fondations nécessaires à l’étude comparative que nous mènerons
dans le chapitre suivant, dans le contexte de la convection naturelle, pour différents régimes d’écoulement.

2.1 historique sur la LBM

La LBM est une méthode développée depuis une trentaine d’années pour simuler les écoulements fluides
à partir d’une approche basée sur la théorie cinétique des gaz. L’efficacité de la méthode et sa rapidité
ont permis d’élargir ses domaines d’applications. On peut citer par exemple : diffusion dans les milieux
poreux, écoulements avec changement de phase, écoulements de convection naturelle, . . .. et plus récemment
l’interaction fluide structure. D’une manière générale, pour n’importe quelle équation aux dérivées partielles
on peut associer un schéma de type LBM basé sur la modification de la fonction d’équilibre et l’opérateur
de collision. Pour les problèmes d’interaction fluide structure où cette dernière est rigide nous avons proposé
un schéma de discrétisation LBM s’appuyant sur une formulation de type pénalisation volumique. Cette
méthode consiste à résoudre les équations du fluide sur le domaine tout entier (fluide et solide), et à rajouter
à l’emplacement du domaine solide une force de pénalisation, consistant à assimiler la structure à un milieu
de perméabilité infinie. La construction de la méthode Boltzmann sur Réseau peut se résumer en deux
étapes paradoxalement indépendantes : le développement de la physique statistique d’une part et l’apparition
des automates cellulaires d’autre part [21]. La première étape nous propulse au milieu du XIXeme siècle
avec l’avènement de la physique statistique basée sur les travaux de l’Américain Josiah Willard Gibbs, des
Anglais Lord Kelvin et James Clerk Maxwell et de l’Autrichien Ludwig Eduard Boltzmann. La physique
statistique, dérivée de la théorie cinétique des gaz, s’intéresse non plus à l’évolution globale des quantités
macroscopiques (densité, vitesse, pression) mais au comportement microscopique de la matière et tente de
décrire l’évolution des particules qui la constitue. Ainsi en 1872, Boltzmann propose sa célèbre équation (1.1)
décrivant l’évolution spatio-temporelle d’une fonction f représentant la distribution des particules ayant une
vitesse donnée à un endroit et à un temps donnés. Cette fonction est souvent appelée fonction de distribution
et dépend donc de l’espace, de la vitesse et du temps f = f(x, c, t) [22] :
 
∂f ∂f Fi ∂f ∂f
+ ci + =
∂t ∂xi m ∂ci ∂t coll
Le membre de droite de l’équation de Boltzmann, appelé opérateur de collision, représente l’e et des
collisions entre particules. Si ce terme est nul, les particulessont alors simplement advectées et soumises
à l’action de la force F présente dans le membre de gauche. Si les particules se rencontrent, l’évolution du

20
système dépend alors de la forme de l’opérateur de collision. A la fin n du XIXeme siècle Boltzmann lui-même
n’avait qu’une idée approximative de cet opérateur et ne considérait que les collisions à deux particules (dites
collisions binaires). On sait aujourd’hui de puis les travaux de Cohen et Dorfman (1970) que l’équation de
Boltzmann avec un tel opérateur de collision ne peut pas s’appliquer aux liquides et aux gaz de haute densité.
Les travaux de Boltzmann, très controversés à l’époque ne furent validés expérimentalement qu’après sa mort
en 1906. Ce n’est que plus tard, dans les années 1920 que les travaux indépendants d’un astronome anglais,
Sydney Chapman et d’un mathématicien suisse, David Enskog s’intéressent à des collisions plus com plexes
mettant en jeu plus de deux particules. Le jeune Enskog eut l’idée d’utiliser la méthode de Hilbert (1912)
pour e ectuer un développement systématique de la fonction de distribution de l’équation de Boltzmann.
Une deuxième étape importante vers la LBM moderne a été franchie en 1988 par McNamara et Zanetti, qui
ont remplacé les particules individuelles des LGCA par une fonction de distribution moyenne, mais toujours
discrète sur le plan directionnel. Cela a permis d’éliminer complètement le bruit statistique de LGCA. Une
simplification majeure a été introduite par Qian et al. la matrice de collision de Higuera est remplacée par
un seul temps de relaxation, ce qui conduit au modèle BGK Par la suite, la LBM s’est développée très
rapidement et est devenue une technique populaire dans les simulations numériques. Plus tard, Lallemand et
Luo et Luo ont montré que la LBM peut être dérivée de l’équation continue de Boltzmann. Par conséquent,
elle peut être considérée comme une forme discrétisée spéciale de l’équation de Boltzmann . À partir de
l’expansion de Chapman Enskog, les équations de Navier-Stokes gouvernantes peuvent être récupérées à
partir de la LBM . En effet, l’utilisation de la LBM couvre aujourd’hui une grande variété de disciplines .
Et comme cette thèse porte sur la simulation de la propagation des ondes acoustiques, il est intéressant de
mentionner quelques applications de la LBM dans ce domaine, comme les ondes élastiques , les ondes de
choc, les ondes aéroacoustiques les ondes sonores et l’acoustic streaming [23].

2.2 Les automates cellulaires

L’histoire des automates cellulaires (cellular automata, CA, en anglais) remonte pratiquement aux années
1950, lorsqu’ils ont été introduits pour la première fois par Stanislas Ulam, John von Neumann et Konrad
Zuse. En général, ils peuvent être caractérisés comme suit :
• CA sont des arrangements réguliers de cellules individuelles de même nature.
• Chaque cellule contient un nombre fini d’états discrets.
• Les états sont mis à jour simultanément (de manière synchrone) à des niveaux temporels discrets.
• Les règles de mise à jour sont déterministes et uniformes dans l’espace et le temps.
• Les règles d’évolution d’une cellule ne dépendent que d’un voisinage local de cellules autour d’elle.
Un automate cellulaire est donc un concept mathématique constitué d’une grille qui est composée de
cellules. Ces cellules sont caractérisées par un nombre d’états précis. L’état d’une cellule à l’instant t est
déterminé par l’état de quatre cellules voisines à l’instant t − 1. À chaque nouvelle unité de temps, les mêmes
règles sont appliquées pour toutes les cellules de la grille, produisant une nouvelle génération de cellules
dépendant entièrement de la génération précédente.
Par exemple en deux dimensions, il y a beaucoup plus de liberté pour organiser les cellules et définir
les quartiers pour les règles de mise à jour. Il existe différentes configurations pour présenter les automates
cellulaires dans une grille bidimensionnelle, comme présenté dans la référence . Par exemple, au début des
années 70, Conway a introduit la notion de "Vie" comme un automate cellulaire synchrone bidimensionnel
qui simule l’évolution d’une société d’organismes vivants. La "vie" est définie par deux règles importantes
[24] :
• Chaque site vivant restera en vie à l’étape suivante s’il a deux ou trois voisins vivants, sinon il mourra.
• Sur un site mort, un nouveau site vivant ne naîtra que s’il y a exactement trois voisins vivants.
La "vie" contient de nombreux modèles qui restent stables d’une itération à l’autre lorsqu’ils ne sont pas
perturbés par d’autres objets (voir la figure 2.1 pour quelques exemples).

21
2.3 Les gaz sur réseau
La technique du gaz sur réseau fait partie de la classe des automates cellulaires. Son objectif principal est
de proposer des modèles capables de simuler les écoulements de fluides. Le principe de base est de simuler
de la manière la plus simple possible le comportement et l’interaction de nombreuses particules individuelles
se déplaçant sur des réseaux carrés ou triangulaires.
Les gaz sur réseau peuvent également être considérés comme des méthodes très simples de la dynamique
moléculaire. Des quantités macroscopiques telles que la densité et la vitesse des particules peuvent être récu-
pérées de ce mécanisme, ce qui permet d’étudier le comportement macroscopique d’un fluide dans différentes
géométries avec ce modèle. Le gaz est modélisé Comme une multitude de sphères dures se déplaçant le long
d’une grille régulière, avec un ensemble discret de vitesses possibles ci pour chaque particule en satisfaisant
les hypothèses suivantes :
• Principe d’exclusion : il ne sera pas possible pour plusieurs particules coexistant dans un même nœud
de se propager sur les éléments du réseau, sauf en raison d’une seule particule par direction.
• Les collisions conservent la masse (le nombre de particules) et la quantité de mouvement.
• Discrétisation de l’espace, du temps, des vitesses et du nombre de particules présentes à un instant
donné en un nœud donné du réseau.
Il existe plusieurs modèles de gaz sur réseau qui permettent de simuler les écoulements de fluides, en
particulier les écoulements de gaz. Les modèles les plus connus dans la littérature sont les modèles HPP et
FHP, qui sont abordés dans l’annexe A. Un autre modèle peut être utilisé. Il s’agit du modèle FCHC décrit
par Jami .

Figure 2.1 – Quelques modèles de l’automate cellulaire "Vie" qui restent stables de génération en génération.

Il est important de mentionner que la méthode de Boltzmann sur réseau est le résultat direct de l’adap-
tation de la méthode des gaz sur réseau à la fonction de distribution des particules f . En effet, au lieu de
travailler avec des particules discrètes du gaz, on utilise plutôt les flux de particules décrits par f dans un
processus d’advection et de collision, qui respectent les lois de la conservation. La technique LBM est vérifiée
par la récupération des équations gouvernantes de Navier-Stokes à l’aide de l’expansion de Chapman-Enskog
.

2.4 Discrétisation des vitesses

La discrétisation de l’espace de vitesse permet de définir le réseau LBM. Il est nécessaire de choisir
un ensemble de vitesses bien réduites pour optimiser le temps de calcul des simulations LBM. Cependant,
le nombre de vitesses doit être suffisant pour obtenir le comportement dynamique de l’écoulement et les
vitesses ne doivent pas être choisies de manière aléatoire. Le choix du réseau LBM est donc très important.
En général, il est nécessaire de choisir un réseau symétrique pour obtenir le comportement de l’écoulement
au niveau macroscopique. Typiquement, les réseaux D1Q3, D2Q9 et D3Q19 sont utilisés respectivement en
1D, 2D et 3D. Le choix du réseau dépend souvent du problème physique étudié. La notation DnQm est
utilisée pour indiquer qu’on a “n” dimensions et “m” vitesses.

22
 Les schémas D1Q3 et D2Q9 seront abordés dans les paragraphes suivants, et le modèle D3Q19 sera
discuté dans le cinquième chapitre.
Pour le réseau D1Q3, il existe trois vecteurs de vitesses c⃗0 , c⃗1 et c⃗2 correspondant aux fonctions f0 , f1 et
f2 , respectivement (voir la figure 2.2). La distance entre deux nœuds est donnée par ∆x et le temps nécessaire
à une particule pour passer d’un nœud à l’autre est ∆t. Ces deux paramètres sont généralement égaux à une
unité LBM (∆x = ∆t = 1). En d’autres termes, c1 = ∆x ∆t et c2 = − ∆t . Alors, pour l’arrangement D1Q3,
∆x

le nombre total de particules à un moment donné ne peut pas dépasser trois particules. Une particule de
vitesse nulle (c⃗0 = ⃗0) réside sur le site central. Les deux autres particules se déplacent vers la gauche ou la
droite dans le processus de transport .

Figure 2.2 – Le schéma D1Q3

schéma D2Q9
Le modèle D2Q9 est largement utilisé dans la littérature, en particulier pour résoudre les problèmes
d’écoulement des fluides . Il comporte neuf vecteurs de vitesses qui sont définis en fonction de leur position
dans le réseau (voir la figure 2.3). Ces vitesses sont c(0, 0), c(0, 1), c(1, 0), c(1, 1), c(0, −1), c(−1, 0), c(−1, −1), c(1, −1)
et c(−1, 1) pour f0 , f1 , f2 , f3 , f4 , f5 , f6 , f7 et f8 , respectivement.

Figure 2.3 – Le schéma D2Q9

La discrétisation des vitesses nous a permis d’avoir une idée sur la propagation (diffusion) des particules
dans chaqJue direction. La question évidente est maintenant de savoir quelle est la probabilité qu’une par-
ticule se propage dans une direction donnée ? Pour répondre à cette question, il est nécessaire de définir un
nouveau paramètre qui est le poids de discrétisation w [25] .
Les poids wi dépendent généralement du choix du réseau LBM et peuvent être trouvés via la quadrature
de Gauss-Hermite . Nous avons par exemple pour le schéma D2Q9 les neuf poids wi suivants :

schémaD2Q5
C’est le réseau le plus utilisé pour la simulation du transfert de chaleur dans la méthode à double
population. Le modèle D2Q5 est basé sur un réseau carré de pas ∆x = ∆y = 1, où chaque motif du réseau
est caractérisé par cinq vitesses discrètes ci (i = 1 . . . 5) (voir Fig. 2.6).

23
 Direction i 0 1 2 3 4 5 6 7 8
4 1 1 1 1 1 1 1 1
wi 9 9 9 9 9 36 36 36 36

Table 2.1 – Les poids de discrétisation pour le modèle D2Q9.

Figure 2.4 – Modèle D2Q5

Les particules fluides se déplacent d’un nœud de la grille vers le nœud voisin avec les vitesses discrètes
qui sont données par :

c1 = (0, 0), les composantes suivant x et y sont nulles.

c2 = (1, 0), la composante suivant x est égale à 1 et suivant y est nulle.
c3 = (0, 1), la composante suivant x nulle et suivant y est égale à 1.
c4 = (−1, 0), la composante suivant x est égale à − 1 et suivant y est nulle
c5 = (0, −1), la composante suivant x est nulle suivant y est égale à − 1.

Les facteurs de poids du réseau D2Q5 sont les suivants :

2 1
w1 = , w2 = w3 = w4 = w5 =
6 6
La diffusivité thermique est reliée au temps de relaxation par la relation suivante :
 
1 1 2
α= τT − c ∆t (2.1)
3 2

2.5 Discrétisation des moments

2.5.1 Méthode de quadrature de Gauss-Hermite
Pour faciliter la discrétisation des moments de la fonction de distribution f, nous utilisons la méthode de
Gauss-Hermite,elle est basé sur l’approximation des intégrales par une somme discrète [26].
Z X
φ(c)f (c) dc ≈ wk φ(ck )f (ck )
k

où ck sont les points de quadrature et wk sont les poids associés.

24
 Figure 2.5 – Réseau D2Q5

2.5.2 Application de la méthode de Gauss-Hermite pour discritiser les moments

.
Après avoir discrétisé l’intégrale en somme, on définit :

fk (x, t) = wk f (x, ck , t)
ainsi, l’équation devient : Z X
φ(c)f (c) dc ≈ φ(ck )fk
k

D’où il est possible de déterminer la densité, la quantité de mouvement et l’énergie cinétique évoquées dans
le premier chapitre à partir de la fonction de distribution discrète fk , respectivement :
X
ρ(⃗x, t)⃗u = fk (2.2)
k
X
ρ(⃗x, t)⃗u = fk⃗ck (2.3)
k

1X
ρ(⃗x, t)Ec (⃗x, t) = fk |⃗ck |2 (2.4)
2
k

2.6 Discrétisation de la fonction de distribution à l’équilibre

Comme illustré dans le chapitre 1, la fonction de distribution à l’équilibre s’écrit sous la forme suivante :

(u)2 (cu)2
 
ρ cu c2
fkeq = √ 1+ − + 2
e− 2RT
2πRT RT 2RT 2(RT )

avec wk = √2πRT
1
.
En remplaçant dans la fonction de distribution, on obtient :

(u)2 (cu)2
 
eq cu c2
− 2RT
fk = ρwk 1 + − + e
RT 2RT 2(RT )2

La vitesse du son dans l’air (cs ) est définie comme :

 
∂p
cs =
∂ρ
Pour le cas d’un gaz idéal, en thermodynamique, on a :

25
 c2s = RT
En substituant RT par c2s dans l’équation , on trouve :

⃗ck · ⃗u (⃗ck · ⃗u)2 |⃗u|2

 
fkeq (⃗ck ) = wk ρ 1 + 2 + −
cs 2c4s 2c2s

Dans le modèle Lattice Boltzmann D2Q9, la vitesse du son du fluide est donnée par cs = √13 , ce qui implique
que c2s = 13 .
Cette valeur est spécifique à la structure discrète du modèle D2Q9 et découle de sa symétrie ainsi que
de la manière dont les vitesses discrètes sont définies. En remplaçant donc c2s par 13 dans l’expression de la
fonction d’équilibre fkeq , on obtient :
 
9 3
fkeq (⃗ck ) = wk ρ 1 + 3(⃗ck · ⃗u) + (⃗ck · ⃗u)2 − |⃗u|2
2 2

2.7 Equation de Boltzmann sur Réseau

2.7.1 Modèle dynamique de Boltzmann sur réseau (BGK) :
L’équation continue de Boltzmann-BGK , peut être discrétisée en remplaçant la fonction continue f par
sa version discrète fk . On obtient :
∂fk 1
+ ⃗ck · ∇⃗r fk = (fkeq − fk )
∂t τ
La discrétisation temporelle mène à l’équation de Boltzmann sur réseau, permettant l’application d’un schéma
de différence finie du premier ordre.

∂fk dfk (⃗x, t) fk (⃗x + ⃗ck ∆t, t + ∆t) − fk (⃗x, t)

+ ⃗ck · ∇⃗x fk = ≈ (2.5)
∂t dt ∆t
En unités LBM, le pas de temps ∆t est égal à 1 (∆t = 1). Alors, la forme finale de la dérivée Lagrangienne
discrétisée est :
∂fk
+ ⃗ck · ∇⃗x fk ≈ fk (⃗x + ⃗ck , t + 1) − fk (⃗x, t) (2.6)
∂t
L’insertion de l’équation (5) dans l’équation (3) conduit à l’équation de Boltzmann sur réseau selon le
modèle BGK.
1 eq
fk (⃗x + ⃗ck , t + 1) − fk (⃗x, t) = (f − fk ) (2.7)
τ k

2.7.2 Modèle thermique de Boltzmann sur réseau :

Puisque les flux thermiques sont constamment présents, un modèle sans chaleur (IEDDF) a été suggéré,
permettant de résoudre le problème de déterminer le nombre de points à partir d’une nouvelle fonction de
distribution pour quantifier l’énergie interne exploitent l’idée que le processus de compression réalisé par la
pression et la dissipation de chaleur visqueuse peut être ignoré dans les flux incompressibles, ce qui permet
de simplifier l’IEDDF. Les nouvelles formes préservent la forme standard de LBM et ne contiennent pas de
termes dégradés.
L’équation d’évolution du nouveau modèle thermique est définie comme suit :
1
gk (x + ck ∆t, t + ∆t) − gk (x, t) = − [gk (x, t) − gkeq (x, t)]
τT
où τT est le temps de relaxation du transfert thermique.
La fonction de distribution d’équilibre correspondante est comme suit [33] :

gkeq = T ωk [1 + 3ϵ̃k · ⃗u]

Les facteurs de poids ωk correspondants sont définis comme suit :

26
 ω0 = 1/3
ω1,4 = 1/6
Les vitesses discrètes (⃗ck ) pour le modèle D2Q5 (Fig.2.5) et la diffusivité thermique (ω) liée au modèle
réseau D2Q5 voir (Table 2.1) :
1
τt = 3α +
2
Le nombre de Nusselt moyen est déterminé par :

⟨u · T ⟩
N umoy = 1 +
α · ∆T /L

2.7.3 Modèle TRT (Two-Relaxation-Time)

Le modèle Two-Relaxation-Time (TRT) est une variante améliorée de la méthode de Boltzmann sur
réseau (LBM), introduite par Ginzburg. Il est considéré comme un cas particulier de la méthode Multiple-
Relaxation-Time (MRT) et vise à améliorer la stabilité numérique par rapport au modèle BGK[16].
Dans le modèle TRT, la fonction de distribution fk est décomposée en deux parties :
— Une composante symétrique fks : elle influence les propriétés physiques de l’écoulement, notamment
la viscosité.
— Une composante antisymétrique fka : elle améliore la dissipation des erreurs numériques.
Les composantes symétrique et antisymétrique sont définies comme suit :

1
fks = (fk + f−k )
2
1
fka = (fk − f−k )
2

Équation de collision TRT

L’équation de collision TRT est donnée par :
1 s 1 a
fk (⃗x + c⃗k ∆t, t + ∆t) − fk (⃗x, t) = − (fk (⃗x, t) − fks,eq ) − (f (⃗x, t) − fka,eq )
τs τa k
On pose que :
1 1
ws = , wa =
τs τa
Donc, l’équation devient :
1 1
fk (⃗x + c⃗k ∆t, t + ∆t) − fk (⃗x, t) = − (ws + wa ) (fk − fkeq ) − (ws − wa ) f−k − f−k
eq

2 2

Calcul de ws :
Le coefficient de relaxation symétrique ws est lié à la viscosité cinématique ν par la relation classique en
LBM :
 
1 1
ν= τs −
3 2
Donc :
1
ws =
3ν + 0.5

27
Calcul de wa :
Pour minimiser la dépendance de la viscosité à la vitesse de glissement, wa est donné par :

8(2 − ws )
wa =
8 − ws

2.7.4 Modèle de Boltzmann sur réseau régularisé d’ordre élevé

La stabilité numérique du modèle LBM peut être considérablement améliorée à l’aide de techniques de
régularisation et de corrections d’ordre supérieur appliquées à l’opérateur de collision. En particulier, le
modèle de Boltzmann sur réseau régularisé d’ordre élevé proposé dans est très stable et permet de capturer
avec précision les motifs turbulents.
L’équation de base pour la méthode LBM est donnée par :

fk (x + ck ∆t, t + ∆t) − fk (x, t) = Ωk (f ) + ∆tFk (2.9)

Notations
— Ωk (f ) : Opérateur de collision (modèle BGK ou régularisé)
— Fk : Terme de force externe

Décomposition de l’Opérateur de Collision

L’opérateur de collision Ωk (f ) est défini comme suit :

Ωk (f ) = −ω(fk − fkeq ) + ψk (2.10)

où :
∆t
— ω= : Taux de relaxation
τ
— ψk : Terme de correction d’ordre élevé (régularisation)

Substitution dans l’Équation de Boltzmann sur Réseau

En substituant Ωk (f ) dans l’équation LBM, on obtient :

fk (x + ck ∆t, t + ∆t) = fk (x, t) − ω(fk − fkeq ) + ψk ∆t + ∆tFk (2.11)

En isolant fk (x + ck ∆t, t + ∆t) :

fk (x + ck ∆t, t + ∆t) = fkeq + (1 − ω)(fk − fkeq ) + ∆t(ψk + Fk ) (2.12)

On définit le terme de non-équilibre :

fkneq = fk − fkeq (2.13)
Ce qui donne la forme finale :

fk (x + ck ∆t, t + ∆t) = fkeq + (1 − ω)fkneq + ∆t(ψk + Fk ) (2.14)

Modèle RLB-4 et Hypothèse de Convergence

Dans le modèle RLB-4, le terme ψk est conçu pour éliminer les moments non-hydrodynamiques d’ordre
supérieur. Pour simplifier l’étude de stabilité, on suppose[27] :
ψk : Négligeable ou absorbé dans Fk
L’équation discrétisée devient alors :

fk (x + ck ∆t, t + ∆t) = fkeq + (1 − ω)fkneq + ∆tFk (2.15)

28
Résumé des Termes
— fk : Fonction de distribution associée à la vitesse discrète k
— fkeq : Fonction de distribution à l’équilibre
— fkneq : Partie non-équilibre de fk (utile dans les modèles régularisés)
— Fk : Terme de force externe
— ω : Paramètre de relaxation lié au temps de relaxation τ
Cette présentation synthétise les étapes clés de l’analyse de stabilité dans le modèle LBM et clarifie les
hypothèses utilisées.

Fonction d’Équilibre Étendue fkeq

La fonction de distribution à l’équilibre, notée fkeq , peut être développée sous forme d’une série de po-
lynômes de Hermite centrée sur la vitesse macroscopique. Cette expansion permet de capturer des effets
d’ordre élevé dans les modèles régularisés et améliore la précision pour les écoulements complexes.

Développement en Série de Hermite

L’expression générale est donnée par :
" N
#
1 X
fkeq = wk ρ 1 + (n)
Hα1 ...αn (ck )uα1 . . . uαn , (2.16)
n=1
n!c2n
s

où :
— wk : poids du réseau
— ρ : densité du fluide
— cs : vitesse du son du réseau
— H(n) : tenseur de Hermite d’ordre n
— u = (ux , uy ) : vitesse macroscopique

Développement Explicite jusqu’à l’Ordre 4 (Modèle D2Q9)

Ordres 1 et 2 (standard)
 
eq,(2) uα ck,α uα uβ 2
(2.17)


fk = wk ρ 1 + + c c
k,α k,β − c δ
s αβ
c2s 2c4s
Ces termes garantissent la conservation de la masse et de la quantité de mouvement.

Ordre 3 (corrections anisotropes)

u2x uy ck,y 2
u2y ux ck,x 2

eq,(3) 2
ck,y − c2s (2.18)


fk = ck,x − cs +
2c6s 2c6s
Ces termes améliorent la précision pour les écoulements asymétriques.

Ordre 4 (stabilité numérique)

eq,(4) u2x u2y 2

ck,x − c2s c2k,y − c2s (2.19)



fk = 8
4cs
Ce terme contribue à la stabilité du schéma, en particulier à hauts nombres de Mach.

Forme Finale Combinée

En combinant les termes jusqu’à l’ordre 4, on obtient l’expression étendue :
"
eq uα ck,α uα uβ
u2 uy ck,y 2
u2y ux ck,x 2
ck,α ck,β − c2s δαβ + x 6 ck,x − c2s + ck,y − c2s


fk = wk ρ 1 + 2
+ 4 6
cs 2cs 2cs 2cs
# (2.20)
u2x u2y 2
ck,x − c2s c2k,y − c2s



+ 8
4cs

29
 Cette formulation est particulièrement utile dans les modèles de régularisation avancée (RLB-4) pour
améliorer la stabilité et la fidélité des simulations.

Fonction de Non-Équilibre fkneq

La fonction de non-équilibre fkneq = fk − fkeq joue un rôle essentiel dans la méthode Lattice Boltzmann
(LBM), en permettant de capturer les effets dissipatifs et les gradients locaux de vitesse. Elle est souvent
développée à l’aide de termes d’ordre supérieur pour améliorer la précision dans les régimes complexes ou
fortement anisotropes.

Formulation générale :
fkneq = wk · [termes d’ordre 2 + ordre 3 + ordre 4] (2.21)

Décomposition des Termes de Non-Équilibre

Ordre 2 – Contraintes Visqueuses
(2) ταβ 2
(2.22)


fk = ck,α c k,β − cs δαβ
2c4s
Origine : Tenseur des contraintes visqueuses ταβ . Rôle : Représente la dissipation visqueuse (analogue à
ν∇2 u dans les équations de Navier–Stokes).

Ordre 3 – Couplages vitesse–gradients

(3) (uy ξxx + 2ux ξxy ) ck,y 2 (ux ξyy + 2uy ξxy ) ck,x 2
fk = (ck,x − c2s )+ (ck,y − c2s ) (2.23)
2c6s 2c6s
Origine : Couplages entre la vitesse macroscopique u et les moments non-hydrodynamiques ξαβ . Rôle :
Corrige les effets anisotropes dans les zones de cisaillement (ex. : couches limites).

Ordre 4 – Stabilisation numérique

(4) u2y ξxx + 4ux uy ξxy + u2x ξyy 2

fk = (ck,x − c2s )(c2k,y − c2s ) (2.24)
4c8s
Origine : Interaction entre les vitesses quadratiques et les moments d’ordre deux. Rôle : Améliore la
stabilité numérique à nombre de Reynolds élevé ou en présence de forts gradients.

Moments Non-Hydrodynamiques ξαβ

Les moments ξαβ mesurent les déviations locales par rapport à l’équilibre :
Q
X
ξαβ = ck,α ck,β (fk − fkeq ) , α, β ∈ {x, y} (2.25)
k=1

Ils constituent une base pour la régularisation des fonctions fk dans les modèles LBM avancés.

Récupération des Grandeurs Macroscopiques

Les variables d’intérêt (masse, vitesse, pression) s’obtiennent à partir des moments de fk :
— Densité massique :
Q−1
X
ρ(x, t) = fk (x, t)
k=0
— Quantité de mouvement :
Q−1
X ∆t
ρu(x, t) = ck fk (x, t) + F(x, t)
2
k=0

où F est une force externe (ex. : gravité).

30
 — Pression :
1 1
p(x, t) = c22 ρ(x, t), avec c2 = √ (ex. : pour D2Q9, c2 = √ )
q 3
Cette formulation non-équilibre enrichie est cruciale pour les simulations de convection naturelle, en
particulier pour des régimes instationnaires et fortement non-linéaires.

2.8 Schémas de relaxation multiples : MRT

Les schémas de relaxation multiples (MRT) constituent une version sophistiquée de la méthode Lattice
Boltzmann (LBM), dans laquelle plusieurs taux de relaxation sont assignés aux moments de distribution.
Une méthode MRT se caractérise par l’emploi de paramètres de relaxation distincts, ce qui améliore la
stabilité, l’exactitude et l’élimination des artefacts numériques comparativement aux méthodes tradition-
nelles[28].
De la même façon que le modèle LBM-BGK calcule les variables macroscopiques avec les fonctions de
distribution, le modèle LBM-MRT permet le calcul de chaque moment grâce aux fonctions de distribution
par la relation suivante :
m = Mf
où :
— M est la matrice de transformation de taille n × n.
— f est le vecteur des fonctions de distribution.
L’équation générale du modèle s’écrit alors :

fk (x̃ + c̃, t + 1) − fk (x̃, t) = −M−1 S (mk − meq )

où :
— r̃ représente la position dans le réseau, et c̃ est le vecteur de vitesse discrète.
— M−1 est l’inverse de la matrice de transformation.
— S est la matrice diagonale de relaxation.
— mk et meq sont les moments de distribution réel et équivalent, respectivement.
Dans le cas du réseau carré D2Q9, les vitesses discrètes C = ∆x ∆t sont définies comme suit :

pour k = 0

 (0, 0)
Ck = (1, 0), (0, 1), (−1, 0), (0, −1) pour k = 1, 2, 3, 4 (directions cardinales)
(1, 1), (−1, 1), (−1, −1), (1, −1) pour k = 5, 6, 7, 8 (directions diagonales)


Dans le modèle D2Q9, la matrice M est élaborée en sélectionnant des moments physiques significatifs. Chaque
ligne de M est liée à un instant précis, tandis que chaque colonne se réfère à l’une des 9 vitesses distinctes
Ck . Les 9 moments que nous voulons extraire à partir des fonctions de distribution fk sont :
8
X
m0 = fk
k=0

8
X
fk c2k,x + c2k,y


m1 =
k=0

8
X
fk c2k,x − c2k,y


m2 =
k=0

8
X
m3 = fk ck,x
k=0

8
X
m4 = fk ck,y
k=0

8
X
m5 = fk ck,x ck,y
k=0

31
 8
X
m6 = fk |ck,x |
k=0

8
X
m7 = fk |ck,y |
k=0

8
X
fk c2k,x c2k,y


m8 =
k=0

La matrice de transformation M est explicitement donnée par :

 
1 1 1 1 1 1 1 1 1
−4 −1 −1 −1 −1 2 2 2 2
 
4
 −2 −2 −2 −2 1 1 1 1 
0
 1 0 −1 0 1 −1 −1 1
0
M = 0 1 0 −1 1 1 −1 −1
0
 −2 2 −2 2 0 0 0 0 
0
 0 0 0 0 1 −1 1 −1
0 0 0 0 0 1 1 −1 −1
0 0 0 0 0 1 −1 −1 1

L’inverse de la matrice M est :

 
4 −4 4 0 0 0 0 0 0
4
 −1 −2 6 −6 0 0 9 0
4
 −1 −2 0 6 −6 −9 0 0
4 −1 −2 −6 6 0 0 9 0
M −1
 
4
= a −1 −2 0 −6 6 −9 0 0
4 2 1 6 3 6 3 0 9
 
4
 2 1 −6 −3 −6 −3 0 −9

4 2 1 −6 3 −6 3 0 9
4 2 1 6 −3 6 −3 0 −9

où le coefficient a = 36 .
1

S est une matrice diagonale composée de l’inverse de chacun des neuf temps de relaxation . Elle est
donnée par :

 
s0 0 0 0 0 0 0 0 0
0 s1 0 0 0 0 0 0 0
 
0 0 s2 0 0 0 0 0 0
 
0 0 0 s3 0 0 0 0 0
 
0
S= 0 0 0 s4 0 0 0 0 
0 0 0 0 0 s5 0 0 0
 
0 0 0 0 0 0 s6 0 0
 
0 0 0 0 0 0 0 s7 0
0 0 0 0 0 0 0 0 s8
Ces valeurs déterminent la rapidité avec laquelle chaque moment atteint son état d’équilibre.
— s0 = s3 = s5 = 1 (relaxation instantanée)
— s1 = s2 = 1.4
— s4 = s6 = 1.2
— s7 = s8 = 3ν+0.5
1
(lié à la viscosité cinématique)

32
Composition du vecteur des moments m
Le vecteur m est constitué de 9 composantes disposées dans l’ordre suivant :

m = (ρ, e, ϵ, jx , qx , jy , qy , pxx , pxy )T

Chaque élément représente une caractéristique physique distincte du fluide :

— ρ : Densité du fluide, qui représente la masse totale par unité de volume.
— e : Moment d’énergie, associé à l’énergie totale transportée par les particules.
— ϵ : Carré de l’énergie, lié aux fluctuations d’énergie dans le système.
— jx et jy : Composantes de l’impulsion dans les directions x et y. Elles sont définies par j = ρ(u, v),
où u et v sont les vitesses macroscopiques.
— qx et qy : Composantes du flux énergétique dans les directions x et y, représentant le transfert d’énergie.
— pxx : Composante diagonale du tenseur de contrainte, qui mesure les efforts normaux.
— pxy : Composante hors diagonale du tenseur de contrainte, qui mesure les efforts de cisaillement.

Moments d’équilibre meq

i
Les moments d’équilibre représentent les valeurs que les moments tendent à atteindre après chaque étape
de relaxation.

meq
0 =ρ

meq 2 2


1 = −2ρ + 3ρ jx + jy

meq 2 2


2 = ρ − 3ρ jx + jy

meq
3 = jx , meq
4 = −jx

meq
5 = jy , meq
6 = −jy

1 2
meq j − jy2


7 =
ρ x
1
meq
8 = jx jy
ρ
À l’opposé des techniques traditionnelles de la dynamique des fluides numérique (CFD), qui impliquent la
résolution directe d’équations différentielles partielles (telles que les équations de Navier-Stokes), la méthode
de Boltzmann sur réseau (LBM) adopte une approche discrète et statistique.
Par exemple, dans le modèle D2Q9, on peut déterminer la densité ρ ainsi que les composantes de vitesse
u et v de la manière suivante :
8
X 8
X 8
X
ρ= fk , ρu = fk ckx , ρv = fk cky
k=0 k=0 k=0

2.9 Implémentation d’une force externe dans la LBM

Pour aborder un problème physique soumis à une force externe comme la convection naturelle, il est
indispensable d’intégrer cette force dans le code LBM employé. Il y a plusieurs modèles disponibles pour
l’implémentation d’une force externe dans la technique LBM. Les modèles les plus couramment cités et
utilisés dans les publications sont celui de Shan et Chen ainsi que celui proposé par Luo[29].

33
Modèle de Shan et Chen
Ce modèle intègre la force de façon indirecte en modifiant l’impulsion à l’échelle locale, soit la quantité
de mouvement du fluide :
X ∆t Fext
j(x, y) = fk ck +
2
k

où Fext est la force externe appliquée.

Modèle de Luo
Le modèle de Luo intègre une force externe directement dans l’équation discrétisée de Boltzmann. Dans
ce modèle, la force externe Fext est discrétisée et intégrée de façon explicite dans l’équation de Boltzmann
via le terme ∆t Fk :

fk (x + ck ∆t, t + ∆t) − fk (x, t) = Ωk + ∆t Fk

La force discrétisée pour chaque direction k est donnée par :
ck · Fext
Fk = w k
c2s
où ck est le vecteur vitesse discret associé à la direction k, et cs est la vitesse du son.

Conclusion
Ce chapitre a présenté les fondements théoriques et numériques de la méthode de Boltzmann sur réseau
(LBM), en mettant l’accent sur les aspects clés de sa formulation discrète. Nous avons détaillé les étapes
essentielles de la discrétisation des vitesses, des moments et de la fonction de distribution, ainsi que les
différents modèles de collision utilisés pour simuler les écoulements fluides. Les modèles BGK, TRT et RLB-
4 ont été introduits, chacun offrant des avantages spécifiques en termes de simplicité, stabilité ou précision.
La LBM se distingue des approches traditionnelles de la dynamique des fluides numérique (CFD) par son
cadre théorique basé sur la physique statistique, permettant une modélisation efficace des phénomènes micro-
scopiques et macroscopiques. Les schémas de relaxation multiples (MRT) et les techniques de régularisation
avancée (comme le modèle RLB-4) améliorent significativement la stabilité et la précision des simulations,
en particulier pour les régimes turbulents.
En conclusion, ce chapitre a posé les bases méthodologiques nécessaires pour aborder les simulations numé-
riques de la convection naturelle dans les chapitres suivants Les formulations obtenues ici nous permettront
de mettre en œuvre les simulations numériques étudiées , où nous évaluerons le comportement de ces modèles
dans des conditions de convection naturelle, en régimes laminaire et turbulent.

34
Chapitre 3

Problèmatique et les Conditions aux

limites

Introduction
Ce chapitre marque la transition entre les fondements théoriques de la méthode de Boltzmann sur réseau
(LBM) et son application concrète à l’étude de la convection naturelle. Après avoir présenté les modèles
numériques (BGK, TRT et RLB-4) ainsi que leurs principes de discrétisation dans le chapitre précédent,
nous nous intéressons ici à leur mise en œuvre pour résoudre un problème physique bien défini : l’écoulement
thermo-convectif dans une cavité carrée différentiellement chauffée.
La convection naturelle, phénomène fondamental en mécanique des fluides et en transfert thermique,
joue un rôle crucial dans de nombreuses applications industrielles et naturelles, telles que les systèmes de
refroidissement, les écoulements atmosphériques ou les processus océaniques. Sa modélisation numérique
représente un défi majeur, notamment en raison de la complexité des interactions entre les transferts de
chaleur et les mouvements du fluide, particulièrement dans les régimes turbulents.
Dans ce chapitre, nous définissons précisément :
— la configuration géométrique et les conditions aux limites du problème étudié
— les hypothèses simplificatrices adoptées pour faciliter la résolution numérique
— les équations gouvernantes adimensionnelles et leur traduction dans le cadre de la LBM
— les spécificités liées à l’implémentation des modèles BGK, TRT et RLB-4
— les conditions aux limites hydrodynamiques et thermiques
Ce travail prépare le terrain pour le chapitre suivant, où les résultats des simulations seront analysés et
comparés, permettant d’évaluer la performance des différents modèles LBM dans la prédiction des régimes
laminaires et turbulents de la convection naturelle.

3.1 Description du problème

Le problème étudié dans ce chapitre porte sur la modélisation numérique de la convection naturelle dans
une cavité carrée remplie d’un fluide newtonien incompressible. La cavité, de côté H, est soumise à un
gradient de température imposé entre ses parois verticales : la paroi gauche est maintenue à une température
élevée Tc , tandis que la paroi droite est maintenue à une température plus basse Tf (Tc > Tf ). Les parois
horizontales (haut et bas) sont considérées comme adiabatiques, empêchant tout transfert de chaleur à travers
elles.

35
 Figure 3.1 – Configuration physique et conditions aux limites.

Il convient de signaler que les présentes simulations ont été réalisées avec les modelés utiliser dans les
mêmes conditions appliquées. Notamment le fluide emprisonné dans l’enceinte est d’un nombre de Prandtl
égal à 0.71 ; son écoulement est considéré comme étant bidimensionnel, incompressible et de régime laminaire
stationnaire, avec l’application d’un gradient de température orthogonal à la force de volume tandis que le
reste de l’enceinte est thermiquement isolé. Le nombre de Rayleigh est choisi entre 106 à 109 .

3.2 Hypothèses simplificatrices

Pour but de faciliter la résolution du problème considéré, les hypothèses simplificatrices suivantes sont
appliquées :
— Le problème est considéré comme bidimensionnel.
— Le fluide est newtonien et incompressible.
— L’écoulement est en régime laminaire permanent.
— L’approximation de Boussinesq est appliquée : la masse volumique du fluide ne dépend que de la
température.
— Les autres propriétés physiques du fluide sont supposées constantes et sont calculées à la température
moyenne du fluide.
— L’enceinte se trouve dans le champ gravitationnel de la Terre.
— Les parois de l’enceinte sont rigides.
— Le transfert de chaleur par rayonnement est considéré comme négligeable.

3.3 Description mathématique du problème

Modèle en deux dimensions
En utilisant les hypothèses simplificatrices, les Équations de continuité , Équation du bilan de quantité
de mouvement et Équation du bilan énergétique prennent respectivement les formes suivantes :
— Équation de Continuité :
∂u ∂v
+ =0 (3.1)
∂x ∂y

36
 — Équations de quantité de mouvement :
Suivant l’axe x :  2
∂ u ∂2u

∂u ∂u 1 ∂p
u +v =− +ν + 2 (3.2)
∂x ∂y ρ ∂x ∂x2 ∂y
Suivant l’axe y :
∂2v ∂2v
 
∂v ∂v 1 ∂p
u +v =− +ν + + gβ(T − Tf ) (3.3)
∂x ∂y ρ ∂y ∂x2 ∂y 2

— Équation d’énergie :

∂2T ∂2T
 
∂T ∂T
u +v =α 2
+ (3.4)
∂x ∂y ∂x ∂y 2

3.4 Modélisation sous forme adimensionnelle

Variables adimensionnelles
Ainsi s’écrivent les variables adimensionnelles :

x y uH vH T − Tf H2
X= , Y = , U= , V = , θ= , P = p (3.5)
H H α α Te − Tf ρα2
L’introduction de ces variables dans les équations (3.1), (3.2), (3.3) et (3.4) on aura :

Équation de continuité

∂U ∂V
+ =0 (3.6)
∂X ∂Y

Équation de quantité de mouvement

Suivant ox :
∂2U ∂2U
 
∂U ∂U ∂P
U +V =− + Pr + (3.7)
∂X ∂Y ∂X ∂X 2 ∂Y 2
Suivant oy :
∂2V ∂2V
 
∂V ∂V ∂P
U +V =− + Pr + + P r Ra θ (3.8)
∂X ∂Y ∂Y ∂X 2 ∂Y 2

Équation d’énergie

∂2θ ∂2θ
 
∂θ ∂θ
U +V = 2
+ (3.9)
∂X ∂Y ∂X ∂Y 2

conditions formulées mathématiquement

Les parois horizontales y = 0 et y = H sont considérées comme adiabatiques.
   
∂T ∂T
= =0 (3.10)
∂y y=0 ∂y y=H

La condition de non-glissement des parois permet d’écrire :

u = v = 0 en y = 0 et y = H

u = v = 0 en x = 0 et x = H

37
 La partie chauffée de la paroi gauche est maintenue à une température uniforme :
T = Tc en x = 0
Le reste de la paroi en x = 0 est thermiquement isolé.
 
∂T
=0
∂x x=0
La paroi droite est maintenue à une température uniforme : T = Tf en x = H

3.4.1 Mise en œuvre pratique de la méthode à double population de Boltzmann

sur réseau (TLBM)
Les conditions aux limites, tant hydrodynamiques que thermiques, associées au phénomène de convection
naturelle sont représentées sur la figure 4.11. La température est fixée à θc sur la portion chauffée de la paroi
ouest, et à θf sur la paroi est, avec θc > θf . Ce gradient thermique engendre une convection naturelle à
l’intérieur de la cavité.
Les parois nord, sud, ainsi que le reste de la paroi ouest, sont considérées comme adiabatiques, c’est-à-dire
thermiquement isolées. Pour modéliser les effets de la force de flottabilité, l’approximation de Boussinesq est
adoptée : le coefficient de dilatation thermique β et la viscosité cinématique ν sont supposés constants, et la
force de flottabilité est exprimée comme une fonction linéaire de la température.

T f = ρβ(Tc − To )⃗g (3.11)

où T0 est la température de référence ; To = (Tc + Tf )/2. Dans notre cas To = Tf .
La force de volume dans la direction oy apparaît dans les équations de Navier-Stokes. La similarité
dynamique dépend de deux paramètres adimensionnels : le nombre de Prandtl (1.10) et le nombre de Rayleigh
(1.11). Dans le présent travail, le nombre de Prandtl est égal à 0.71.
Dans la méthode à double population de Boltzmann sur réseau, on aura besoin de modéliser le champ
hydrodynamique et le champ thermique, donc on aura besoin de deux fonctions de distributions différentes
et deux modèles réseaux différents : D2Q9 et D2Q5.

3.4.2 Conditions aux limites adaptées à (LBM)

À partir du processus de propagation, les fonctions de distribution orientées vers l’extérieur du domaine
(Fig.3.2) sont connues, tandis que celles orientées vers l’intérieur du domaine sont inconnues et doivent être
déterminées. Il est nécessaire d’appliquer des conditions aux limites.

Figure 3.2 – les distributions manquantes sur la frontière inférieure pourle modèle D2Q9 à gauche et le
modèle D2Q5 à droite

38
 À titre d’illustration, la (Fig. 14) présente les fonctions de distribution absentes, représentées par des
lignes discontinues, dont la détermination s’effectue selon la procédure suivante :
Pour l’écoulement du fluide (Rebond)
Cette condition aux limites est appliquée dans un modèle dynamique de Boltzmann sur réseau lorsque une
particule atteignant une frontière de type paroi solide rebondit sur celle-ci et réintègre l’écoulement dans sa
direction initiale. Le processus de collision ne se produit pas à la frontière, mais uniquement dans la région
interne du fluide[30].
Les fonctions de distribution f5 , f8 et f9 sont connues par le processus de propagation. Lorsque ces
fonctions rencontrent le mur, le schéma de rebond est appliqué. Ainsi, on a :

f7 = f9
f3 = f5
f8 = f6

Le modèle de rebond garantit la conservation de la masse et de la quantité de mouvement. Il peut être

appliqué aux obstacles ainsi qu’aux parois, en tenant compte du fait que la vitesse est nulle en tout point
sur les surfaces solides.
Pour la température
La condition de Dirichlet en température imposée est utilisée pour modéliser des frontières soumises à des
températures différentes[31] :

θx=0 = θc = 1 (portion chauffée de la frontière ouest)

θx=H = θf = 0 (frontière est)

Les fonctions de distribution manquantes sont calculées selon l’expression suivante :

4
X
θ= gj (3.26)
j=0

Pour les frontières adiabatiques (nord, sud et le reste de la frontière ouest), un schéma d’extrapolation simple
est utilisé, conformément à [ref39], pour imposer une condition de Neumann (absence de flux thermique).
Concernant la frontière sud, la fonction de distribution est estimée selon la relation proposée dans [ref33] :

4θ1 − θ2
θ0 = (3.27)
3

où les indices 0, 1 et 2 désignent respectivement le nœud frontière, le premier et le deuxième voisin les plus
proches.

3.5 Mise en oeuvre informatique

La mise en œuvre de la méthode de Boltzmann sur réseau thermique nécessite un code informatique
(MATLAB). En suivant l’organigramme basé sur la méthode à double population de Boltzmann , un code
Matlab a été développé pour simuler le comportement du fluide ainsi que le transfert de chaleur à l’intérieur
de l’enceinte carrée.

39
 Figure 3.3 – Organigramme de calcul de TLBM.

Conclusion
Ce chapitre a permis de formaliser la problématique physique de la convection naturelle dans une ca-
vité carrée chauffée différemment, en définissant avec rigueur les paramètres géométriques, les hypothèses
simplificatrices, les équations gouvernantes adimensionnelles et les conditions aux limites pertinentes. Il a
également précisé la manière dont ces éléments sont traduits dans le cadre numérique de la méthode de
Boltzmann sur réseau, notamment à travers les modèles BGK, TRT et RLB-4. La mise en œuvre informa-
tique des conditions aux limites hydrodynamiques et thermiques, ainsi que l’architecture du code utilisé pour
les simulations, ont été détaillées pour assurer la reproductibilité et la robustesse de l’étude. Ces fondations
solides préparent l’analyse comparative qui sera menée dans le chapitre suivant, afin d’évaluer la performance
de chaque modèle dans la prédiction des régimes d’écoulement laminaire et turbulent.

40
Chapitre 4

Résultas et discussion

4.1 Introduction
Ce dernier chapitre clôture notre étude en synthétisant les principaux résultats obtenus et en mettant
en lumière les contributions majeures de ce travail. À travers une analyse comparative des modèles BGK-1,
TRT-2 et RLB-4, nous avons exploré leur performance dans la simulation de la convection naturelle, couvrant
des régimes allant du laminaire au turbulent.
Les résultats ont révélé des différences significatives en termes de stabilité, précision et capacité à capturer
les structures complexes des écoulements, notamment pour des nombres de Rayleigh élevés.
Dans cette partie, nous revenons sur les objectifs initiaux du projet, résumons les observations clés et
discutons des implications pratiques de nos résultats. Nous abordons également les limites des modèles
étudiés et proposons des pistes pour des recherches futures, visant à améliorer encore la modélisation des
phénomènes thermo-fluidiques complexes.
Enfin, cette synthèse souligne l’importance des méthodes de Boltzmann sur réseau (LBM) dans le domaine de
la mécanique des fluides et des transferts thermiques, tout en ouvrant des perspectives pour leur application
à des problèmes industriels et environnementaux plus vastes.

4.2 Résultats et discussions

Dans cette étude, nos simulations numériques visent à comparer systématiquement les performances des
différents modèles hybrides de Boltzmann sur réseau. Pour couvrir l’ensemble du spectre des comportements
convectifs, nous avons exploré une gamme exceptionnellement large de nombres de Rayleigh, s’étendant de
10 à 10¹². Cette approche nous permet d’observer la transition progressive depuis les écoulements laminaires
stables jusqu’aux régimes de turbulence pleinement développée. Notre analyse se structure en trois parties
distinctes, chacune ciblant une plage spécifique du nombre de Rayleigh et révélant des phénomènes physiques
différents. Pour chaque configuration, nous présentons des résultats détaillés incluant les champs thermiques
instantanés et moyennés dans le temps, les distributions des lignes de courant, ainsi que des indicateurs
turbulents avancés tels que l’énergie cinétique turbulente, la variance de température et la vitesse quadratique
moyenne. Une attention particulière est portée à l’évolution du nombre de Nusselt, paramètre clé quantifiant
l’efficacité du transfert thermique à travers la cavité. L’examen des statistiques de premier et second ordre
nous permet d’évaluer non seulement la précision globale des modèles, mais aussi leur capacité à capturer
les subtilités des interactions fluide-thermique à différentes échelles spatiales et temporelles.

41
4.3 Validation du modèle
Cette figure présente une validation du modèle numérique en comparant la courbe obtenue (en rouge
pointillé) à celle issue de l’article de référence (en noir). La forte concordance entre les deux courbes confirme
la précision et la fiabilité du modèle développé. Les écarts mineurs observés peuvent s’expliquer par des
différences de discrétisation ou de paramètres numériques, sans remettre en cause la qualité globale de la
simulation. Cette validation montre que le modèle reproduit fidèlement les résultats de la littérature et peut
être utilisé pour des analyses plus complexes.

Figure 4.1 – Comparaison entre le modèle d’article et le modèle obtenue par la simulation.

4.4 Choix du nombre de Rayleigh

Le choix des nombres de Rayleigh dans cette étude repose sur la volonté d’explorer l’ensemble des régimes
caractéristiques de la convection naturelle, du laminaire au turbulent. En sélectionnant des valeurs allant
de 106 à 1012 , l’objectif est de tester la robustesse et la précision des modèles BGK, TRT et RLB-4 face
à des écoulements de complexité croissante. Ce spectre permet non seulement de valider les modèles dans
des conditions bien établies, mais aussi d’analyser leur comportement face à l’apparition d’instabilités et
au développement de structures turbulentes. Ce choix réfléchi assure ainsi une évaluation complète des
performances numériques dans des contextes variés de la convection thermique. Ainsi que Le nombre de
Prandtl (Pr = 0, 71), caractéristique des gaz (comme l’air), est fixé pour toutes les simulations afin d’isoler
l’impact de Ra sur les dynamiques thermo-fluidiques

42
4.5 Evolution du Nombre de Nusselt en fonction du temps de re-
laxation
La (fig 4.2) montre l’évolution du nombre de Nusselt moyen (Nu) en fonction du temps (T), pour une
convection naturelle avec un nombre de Rayleigh Ra = 109 et un nombre de Prandtl Pr = 0,71. Les résultats
sont obtenus à l’aide des trois modèles numériques présentés dans le deuxième chapitre, utilisés dans la
méthode de Boltzmann sur réseau : BGK-1, TRT-2 et RLB-4.

Figure 4.2 – Variation des nombres de Nusselt moyens pour Ra = 109 et P r = 0, 71 : (a) BGK-1 ; (b)
TRT-2 ; (c) RLB-4.

Ces résultats mettent en évidence la sensibilité des modèles au raffinement du maillage dans le cadre
d’une simulation pseudo-directe (pDNS). Le modèle BGK-1 (figure a), basé sur un schéma de premier ordre,
présente une forte dépendance à la grille : les résultats ne convergent de manière satisfaisante qu’avec une grille
très fine (4012 ). En revanche, le modèle TRT-2 (figure b) offre une meilleure stabilité numérique, avec une
convergence plus rapide et une sensibilité réduite à la densité du maillage. Enfin, le modèle régularisé RLB-4
(figure c), couplé à un schéma de résolution d’ordre élevé, se distingue par une remarquable insensibilité à la
taille de la grille : les trois courbes se superposent presque parfaitement dès la grille 2012 . Cette robustesse
illustre l’efficacité des modèles d’ordre élevé dans une approche pDNS, qui vise à simuler fidèlement les
grandes structures turbulentes sans résoudre les plus petites échelles. Contrairement à la DNS classique,
la pDNS permet de s’affranchir des contraintes liées à la résolution des échelles de Kolmogorov, ce qui est
particulièrement adapté au cadre mésoscopique du LBM. Ainsi, les modèles hybrides de Boltzmann sur
réseau, notamment RLB-4, apparaissent comme des outils numériques prometteurs pour des simulations
efficaces et précises dans les régimes fortement convectifs.

43
4.5.1 Évolution du nombre de Nusselt moyen pour différents nombre de Ray-
leigh
Les figures montrent l’évolution du nombre de Nusselt moyen (Nu) en fonction du temps (T) pour
différentes valeurs du nombre de Rayleigh (Ra) et pour plusieurs modèles numériques : BGK-1, TRT-2 et
RLB-4.

Figure 4.3 – Variation des nombres de Nusselt moyens pour Pr = 0,71 : a) Ra = 6,4 × 108 , b) Ra = 2 × 109 .

Comme on pouvait s’y attendre, les modèles BGK-1 et TRT-2 montrent des limites dans la prédiction du
transfert thermique aux nombres de Rayleigh élevés. En particulier, le modèle BGK-1 tend à sous-estimer
significativement le taux de transfert thermique dès que le nombre de Rayleigh atteint ou dépasse 2 × 109 ,
tandis que cette sous-estimation apparaît pour le modèle TRT-2 lorsque Ra ≥ 1010 . Cette limitation est
probablement liée à leur formulation simplifiée, qui ne parvient pas à capturer avec précision les mécanismes
de transport thermique intensifiés par la turbulence à ces régimes.
En revanche, le modèle RLB-4 se distingue par sa capacité à reproduire de manière fidèle les valeurs
moyennes du nombre de Nusselt, et ce, même sans raffinement de la grille. Cela suggère une meilleure
aptitude à modéliser les effets de la turbulence thermique dans des conditions extrêmes.
Bien que le modèle TRT-2 parvienne à fournir des résultats raisonnables pour les statistiques de premier
et second ordre, ainsi que pour les valeurs globales de N u, tant que Ra ≤ 1010 , il présente néanmoins une
erreur non négligeable de 17,5 % sur le nombre moyen de Nusselt. Cette erreur, bien qu’acceptable dans
certains contextes, indique une perte de précision au niveau des échanges thermiques globaux.
À l’inverse, le modèle RLB-4 s’avère plus robuste pour décrire le comportement turbulent du fluide, en
particulier les structures globales de la turbulence thermique développée, ce qui en fait un choix plus fiable
pour les simulations à hauts nombres de Rayleigh.

44
4.6 Résultats par régime laminaires 104 ≤ Ra ≤ 108
4.6.1 Le champ de température (Isothermes)
Les résultats des simulations, affichés sous forme d’isothermes (Fig.4.4), comparent les champs de tem-
pérature prédits par les modèles BGK-1, TRT-2 et RLB-4 dans une cavité carrée soumise à un écoulement
de convection naturelle à Ra = 106 .

Figure 4.4 – Résultats du champ thermique pour les modèles BGK-1 , TRT-2 et RLB-4 à Ra = 106 ..

4.6.2 Les lignes de courant

Les (Figures 4.5 ) illustrent les lignes de courant obtenues avec les modèles BGK-1, TRT-2, et RLB-4
pour un écoulement en convection naturelle à Ra = 106 .

Figure 4.5 – Résultats du lignes de courant pour les modèles BGK-1 , TRT-2 et RLB-4à Ra = 106 ..

L’analyse des isothermes et des lignes de courant pour les modèles BGK1, TRT2 et RLB-4, dans le cas
d’un écoulement avec un nombre de Rayleigh (Ra) de 106 et un nombre de Prandtl (Pr) de 0.71, révèle
une similitude frappante entre les résultats obtenus par ces différentes approches. Les isothermes, qu’ils
soient générés par BGK1, TRT2 ou RLB-4, présentent des structures identiques en termes de distribution
spatiale et de valeurs numériques, avec des températures variant symétriquement autour de zéro. Cette
uniformité suggère que les trois modèles captent efficacement la même physique sous-jacente, caractérisée par
un écoulement laminaire stable et dépourvu de bifurcations complexes. De même, les lignes de courant, bien
que non détaillées visuellement dans les données fournies, semblent reproduire des motifs identiques pour
chaque modèle, ce qui confirme l’absence de différences significatives dans la dynamique de l’écoulement.
Cette cohérence entre les modèles s’explique probablement par la nature simple et bien ordonnée du régime
laminaire étudié, où les premières approximations suffisent à décrire toutes les échelles pertinentes sans
nécessiter de corrections supplémentaires. Ces résultats soulignent la robustesse des méthodes hybrides dans
des configurations où les effets turbulents sont négligeables.

45
4.6.3 Profil de température adimensionnelle et de la vitesse horizontale
La Figure 4.6 illustre le profil de température adimensionnelle dans la cavité pour Ra = 106 . On observe
une couche limite thermique prononcée près des parois, avec un gradient abrupt, tandis que la région centrale
présente une stratification régulière. Ce comportement, typique d’un écoulement laminaire, valide la précision
des schémas numériques utilisés (BGK, TRT, RLB) pour les faibles nombres de Rayleigh.

Figure 4.6 – Profils de température dans une cavité chauffée latéralement pour Ra = 106 et P r = 0, 71
dans les modeles BGK-1 TRT-2 RLB-4.

La Figure 4.7 montre le profil vertical de la vitesse horizontale pour Ra = 106 . On observe une distribution
parabolique caractéristique, avec des vitesses maximales près des parois et une annulation au centre. Ce profil
symétrique, associé à des nombres de Reynolds locaux faibles, confirme la nature laminaire de l’écoulement
pour ces paramètres. La comparaison avec les solutions de référence valide la précision du schéma numérique
employé.

Figure 4.7 – Profils de température dans une cavité chauffée latéralement pour Ra = 106 et P r = 0, 71
dans les modeles BGK-1 TRT-2 RLB-4.

Pour ces résultats, il ressort clairement que les modèles BGK1, TRT2 et RLB4 présentent les mêmes
valeurs de température et de vitesse. Aucun écart significatif n’est observée en écoulement laminaire.

46
4.7 Résultats par régime Transition vers la turbulence (6, 4.108 <
Ra < 1010 )
Pour Ra = 6,4 × 108
Les figures comparent trois méthodes de simulation thermique, BGK-1, TRT-2 et RLB-4, selon plusieurs
critères visuels, pour un cas à Ra = 6,4 × 108 .

47
Figure 4.8 – Champs thermiques moyennés (a, b, c) et instantanés (d, e, f), lignes de courant moyennées
dans le temps (g, h, j) ( a, d, g) BGK-1 ; (b, e, h) TRT-2 ; c, f, j) RLB-4 .

Comme on pouvait s’y attendre, les modèles BGK-1 et TRT-2 montrent des limites dans la prédiction du
transfert thermique aux nombres de Rayleigh élevés. En particulier, le modèle BGK-1 tend à sous-estimer
significativement le taux de transfert thermique dès que le nombre de Rayleigh atteint ou dépasse 2 × 109 ,
tandis que cette sous-estimation apparaît pour le modèle TRT-2 lorsque Ra ≥ 1010 . Cette limitation est
probablement liée à leur formulation simplifiée, qui ne parvient pas à capturer avec précision les mécanismes
de transport thermique intensifiés par la turbulence à ces régimes.
En revanche, le modèle RLB-4 se distingue par sa capacité à reproduire de manière fidèle les valeurs
moyennes du nombre de Nusselt, et ce, même sans raffinement de la grille. Cela suggère une meilleure
aptitude à modéliser les effets de la turbulence thermique dans des conditions extrêmes.
Bien que le modèle TRT-2 parvienne à fournir des résultats raisonnables pour les statistiques de premier
et second ordre, ainsi que pour les valeurs globales de N u, tant que Ra ≤ 1010 , il présente néanmoins une
erreur non négligeable de 17,5 % sur le nombre moyen de Nusselt. Cette erreur, bien qu’acceptable dans
certains contextes, indique une perte de précision au niveau des échanges thermiques globaux.
À l’inverse, le modèle RLB-4 s’avère plus robuste pour décrire le comportement turbulent du fluide, en
particulier les structures globales de la turbulence thermique développée, ce qui en fait un choix plus fiable
pour les simulations à hauts nombres de Rayleigh.

4.8 Résultats par régime fortement turbulents (Ra ≥ 1010 )

Pour Ra = 2 ×109
parmi les méthodes comparées, RLB-4 se distingue par sa haute précision et sa capacité à capturer les
phénomènes turbulents complexes. TRT-2 offre une bonne finesse dans les transitions thermiques, tandis que
BGK-1 reste plus simple mais cohérent. RLB-4 est la plus adaptée aux régimes fortement instables liés à ce
nombre de Rayleigh élevé.

48
Figure 4.9 – Champs thermiques moyennés (a, b, c, d) et instantanés (e, f, g, h), lignes de courant moyennées
dans le temps (i, j, k, l) ( a, e, i) BGK-1 ; (b, f, g) TRT-2 ; c, g, k) RLB-4

49
 Dans le cas de BGK-1, le profil thermique moyen met en évidence des couches quasi verticales avec de
légères ondulations, caractérisant une convection bien organisée. Sous TRT-2, une structure similaire est
observée, avec des oscillations latérales légèrement plus régulières, traduisant un régime convectif un peu
plus stabilisé. En RLB-4, les couches de température deviennent encore plus homogènes, ce qui suggère une
dissipation accrue des petites structures thermiques. Toujours en BGK-1, un motif ondulé serré et complexe
apparaît, révélateur d’une convection turbulente active, marquée par des structures fines de plumes ther-
miques. Avec TRT-2, les motifs restent denses mais sont visiblement plus ordonnés, avec des ondulations
mieux alignées. En RLB-4, la structure est comparable mais montre un léger élargissement des motifs ther-
miques, signe d’une turbulence en diminution. Dans la configuration BGK-1, on distingue de larges rouleaux
convectifs accompagnés de structures secondaires fines, indiquant une interaction complexe entre grandes et
petites échelles. En TRT-2, la topologie des rouleaux reste similaire, mais avec un léger lissage des struc-
tures secondaires. Enfin, en RLB-4, le motif devient plus dispersé, laissant émerger de nombreuses structures
secondaires, ce qui traduit une convection plus désorganisée.

Performance TRT-2 vs RLB-4 pour Ra = 1011 et P r = 0.71 :

Pour Ra = 1011
Le graphique présente la configuration d’un flux de convection naturelle turbulente pour Ra = 1011 et
P r = 0,71, réalisée à l’aide des modèles TRT-2 (en haut) et RLB-4 (en bas). Elle met en évidence les champs
thermiques instantanés et moyens, les trajectoires de courant moyennes, l’énergie cinétique turbulente ainsi
que la vitesse verticale RMS. On note une dynamique grandement instable, caractérisée par des panaches
thermiques puissants et un mélange turbulent bien établi, particulièrement dans le cas RLB. L’analyse dé-
montre l’impact du modèle sur la résolution des structures turbulentes et le transfert de chaleur.

Figure 4.10 – Champs thermiques et d’écoulement instantanés et moyennés temporellement pour les mo-
dèles TRT-2 et RLB-4 avec Ra = 1011 et Pr = 0,71.

La figure 10 met en évidence le comportement turbulent de la convection naturelle à très haut nombre
de Rayleigh (Ra = 1011 ) pour P r = 0,71, en comparant deux modèles LBM : TRT-2 (haut) et RLB-4
(bas). On observe que les champs de température instantanés (a, f) présentent de fins panaches thermiques
et des structures tourbillonnaires complexes, caractéristiques d’un régime pleinement turbulent. Les champs
thermiques moyennés (b, g) révèlent une stratification thermique influencée par les panaches proches des
parois. Les lignes de courant moyennées (c, h) indiquent une organisation cohérente de l’écoulement malgré
l’intensité du chaos instantané. L’énergie cinétique turbulente (d, i) est fortement concentrée au centre,
signe d’un brassage actif, tandis que la vitesse verticale RMS (e, j) montre des zones d’intenses fluctuations
verticales, en particulier près des parois chaudes et froides. Globalement, la figure confirme la présence d’un

50
régime de convection thermique très instable, avec une légère supériorité du modèle RLB-4 pour capturer
les détails fins des structures turbulentes.

Pour Ra = 1012
La figure 11 analyse un écoulement turbulent (Ra = 1012 ) avec stratification (P r = 0,71). Le champ
instantané (a) montre des structures thermiques complexes, typiques de la turbulence. Le champ moyen (b)
révèle des motifs thermiques stables. Les lignes de courant moyennes (c) indiquent des tourbillons influen-
çant les transferts. L’énergie cinétique turbulente (d) localise les zones de forte activité, tandis que la vitesse
verticale RMS (e) identifie les mouvements convectifs dominants. Ces résultats illustrent l’interaction entre
turbulence et stratification thermique.

Figure 4.11 – Analyse thermique et dynamique de l’écoulement.

La figure 11 révèle la dynamique d’un écoulement de convection turbulente fortement stratifiée (Ra =
1012 , P r = 0,71). Le champ thermique instantané (a) présente des structures fines et désorganisées, signature
d’une turbulence pleinement développée où les fluctuations thermiques sont importantes. En contraste, le
champ moyen (b) montre une organisation à grande échelle, avec des panaches thermiques caractéristiques.
Les lignes de courant moyennes (c) mettent en évidence des rouleaux convectifs dominants qui structurent
l’écoulement.
L’énergie cinétique turbulente (d) se concentre principalement aux interfaces entre les panaches, zones
de forts gradients où la dissipation est maximale. La vitesse verticale RMS (e) quantifie l’intensité des
mouvements convectifs, révélant comment la stratification modérée (P r = 0,71) filtre sélectivement les
fluctuations verticales tout en permettant une turbulence horizontale soutenue.

4.9 Bilan des performances des modèles LBM

Suite à une analyse rigoureuse des résultats numériques, ce bilan met en perspective les performances
des modèles Lattice Boltzmann (LBM) dans la modélisation des régimes de convection naturelle, couvrant
un large spectre de nombres de Rayleigh, de 104 à 1012 . L’évaluation porte à la fois sur la capacité physique
des modèles à restituer avec fidélité les structures thermofluidiques caractéristiques des régimes laminaires,
transitoires et turbulents, ainsi que sur leur stabilité et efficacité numérique face aux contraintes de calcul.
Cette étude souligne l’impact déterminant des paramètres numériques — notamment la résolution spatiale
de la grille et les temps de relaxation — sur la précision et la convergence des simulations. Elle met en lumière
les compromis entre robustesse numérique et fidélité physique propres à chaque modèle. En particulier,
l’analyse différenciée clairement la réponse des modèles BGK-1, TRT-2 et RLB-4, en fonction de leur capacité
à capturer les dynamiques complexes liées à l’apparition et au développement de la turbulence.
À travers cette synthèse, il devient possible d’apprécier précisément comment ces trois approches LBM
se comportent numériquement et physiquement selon le régime d’écoulement étudié. Voici les conclusions
détaillées issues de cette étude intégrée.

51
Régime laminaires (104 ≤ Ra ≤ 108 )
Dans ce régime caractérisé par un écoulement stable et régulier, les trois modèles BGK-1, TRT-2 et
RLB-4 affichent d’excellentes performances. Ils reproduisent fidèlement les distributions de température et
les champs de vitesse attendus, ainsi que les lignes de courant classiques d’un écoulement laminaire dans
une cavité chauffée. La robustesse et la précision sont maintenues même avec des maillages de résolution
modérée. Ce régime ne fait pas apparaître d’artefacts numériques majeurs, ce qui valide l’applicabilité des
trois modèles dans des conditions simples et bien contrôlées. Le nombre de Nusselt moyen calculé par chaque
modèle correspond bien aux références théoriques et expérimentales, attestant de leur fiabilité dans cette
plage.

Régime de transition vers la turbulence (6, 4.108 < Ra < 1010 )

Ce régime intermédiaire, où les instabilités commencent à se développer, met en lumière les différences
entre les modèles. Le modèle BGK-1 montre une certaine sensibilité à la résolution de la grille, ce qui peut
affecter la qualité des résultats lors de la capture des premières structures dynamiques instables. Le modèle
TRT-2 offre une meilleure stabilité numérique, permettant une représentation plus précise des structures
transitionnaires émergentes et une meilleure gestion des effets non-linéaires naissants. Le modèle RLB-4,
grâce à sa régularisation d’ordre élevé, se distingue par sa capacité à gérer ces écoulements complexes avec
une moindre sensibilité aux paramètres numériques. Il préserve une bonne stabilité et une représentation
fine des variations locales de température et de vitesse, ce qui est crucial dans ce régime où la dynamique
devient plus chaotique.

Régime turbulent pleinement développé (Ra ≥ 1010 )

Lorsque le flux devient turbulent, la complexité des structures thermo-fluidiques croît fortement. Le mo-
dèle RLB-4 montre alors ses meilleures performances. Il reste stable et précis même pour des nombres de
Rayleigh très élevés, jusqu’à 1012 dans certains cas. Sa formulation régularisée réduit les instabilités numé-
riques typiques des modèles classiques et permet de capturer les petites échelles turbulentes et les fluctuations
dynamiques avec finesse. Le modèle TRT-2 conserve une bonne capacité à gérer le régime turbulent, bien
que ses performances soient légèrement plus sensibles au maillage, nécessitant un raffinement plus important
pour garantir la précision des résultats. En revanche, le schéma RLB-4 résout avec succès les structures
turbulentes grâce à une approximation d’ordre supérieur. Il est très important de souligner que le modèle
RLB-4 est stable dans ce type de régimes. Le modèle BGK-1 présente ici ses limites, avec une stabilité plus
fragile et une plus grande dépendance aux paramètres numériques, ce qui peut conduire à des approximations
moins fiables dans les zones de forte turbulence.

Les résultats de cette étude démontrent clairement que la performance des modèles Lattice Boltzmann
dépend fortement du régime d’écoulement considéré, mettant en lumière les forces et limites spécifiques de
chaque approche.
• Le modèle BGK-1, par sa simplicité algorithmique, s’impose comme la solution la plus efficace pour
les écoulements laminaires. Il offre une bonne précision dans les configurations basiques, avec des temps
de calcul réduits et une implémentation directe. Cependant, cette simplicité devient un handicap dès
que le régime évolue vers des écoulements plus instables ou turbulents. Dans ces conditions, la stabilité
numérique du BGK-1 se dégrade, limitant sa capacité à capturer fidèlement les structures complexes
et les phénomènes dynamiques propres aux régimes transitoires et turbulents.
• Le modèle TRT-2 constitue un compromis équilibré entre stabilité numérique et précision physique,
particulièrement adapté aux régimes transitoires. Il améliore significativement la gestion des effets non
linéaires et permet une meilleure résolution des instabilités débutantes. Cette flexibilité le rend capable
d’aborder des écoulements modérément turbulents, avec une robustesse accrue par rapport au BGK-
1. Néanmoins, la sensibilité de ce modèle aux paramètres numériques, notamment la résolution du
maillage, exige des ajustements rigoureux pour garantir la qualité des simulations dans des régimes à
turbulence plus développée.

52
 • Le modèle RLB-4, grâce à sa formulation régularisée d’ordre élevé, établit une nouvelle référence
en termes de robustesse et de fidélité dans les simulations d’écoulements turbulents. Il excelle à
stabiliser les phénomènes non hydrodynamiques, ce qui se traduit par une meilleure tenue numérique
même pour des nombres de Rayleigh très élevés. Cette capacité permet au RLB-4 de reproduire avec
précision les structures fines et les interactions complexes qui caractérisent la convection naturelle
fortement turbulente. Par ailleurs, son insensibilité notable à la densité de la grille en fait un outil
particulièrement adapté aux simulations exigeantes en ressources, offrant un compromis optimal entre
qualité des résultats et coût de calcul.
Ce tableau résume de manière synthétique les performances et limites des modèles LBM selon les régimes
d’écoulement.

Modèle LBM Régime d’écoule- Points forts Limites / Contraintes Usage recommandé
ment adapté
BGK-1 Laminaire Cas simples, écoulements
• Simplicité algorith- • Stabilité et précision
laminaires uniquement
mique faibles dans les ré-
• Bonne précision pour gimes transitoires et
écoulements basiques turbulents
• Faible coût de calcul • Difficultés à capturer
les phénomènes com-
plexes

TRT-2 Transitoire à turbu- Études des régimes de

• Compromis entre sta- • Sensible aux para-
lence modérée transition et turbulences
bilité et précision mètres numériques
modérées
• Meilleure gestion des (maillage)
effets non linéaires • Nécessite réglages fins
• Robustesse accrue par pour régimes turbu-
rapport au BGK-1 lents plus forts

RLB-4 Turbulent pleine- Simulations avancées de

• Formulation régulari- • Complexité numérique
ment développé turbulence forte et com-
sée d’ordre élevé plus élevée
plexes
• Très bonne stabilité • Coût calcul plus im-
numérique portant
• Capture fine des struc-
tures turbulentes
• Insensible à la résolu-
tion de la grille

Table 4.1 – Comparaison des modèles LBM selon le régime d’écoulement

En bref, cette étude confirme que le choix du modèle LBM doit être rigoureusement adapté au régime
d’écoulement visé et aux objectifs de la simulation. La simplicité du BGK-1 convient parfaitement aux cas
laminaire basiques, tandis que le TRT-2 offre une extension fiable vers des régimes plus instables. Pour les
défis liés aux écoulements turbulents complexes, le RLB-4 se positionne comme le modèle incontournable,
combinant précision, stabilité et efficacité numérique. Cette hiérarchisation claire permet d’orienter les futures
études et applications vers la meilleure adéquation entre modélisation physique et contraintes pratiques.

Conclusion
Notre étude comparative des modèles BGK-1, TRT-2 et RLB-4 sur une large gamme de nombres de
Rayleigh (104 à 1012 ) démontre que leur performance dépend fortement du régime d’écoulement.
En régime laminaire (Ra ≤ 108 ), les trois modèles présentent des performances similaires et excellentes,
avec un avantage pour BGK-1 en termes de simplicité et rapidité de calcul. Dans le régime de transition
(108 < Ra < 1010 ), TRT-2 offre une meilleure stabilité numérique que BGK-1, tandis que RLB-4 se distingue
par sa robustesse face aux instabilités naissantes et sa faible sensibilité aux paramètres numériques.
En régime pleinement turbulent (Ra ≥ 1010 ), RLB-4 démontre une supériorité incontestable, maintenant
stabilité et précision même à Ra = 1012 . TRT-2 reste performant jusqu’à Ra ≈ 1010 mais sous-estime ensuite
le transfert thermique, tandis que BGK-1 atteint ses limites dès 2 × 109 . Cette hiérarchisation permet
d’optimiser le choix du modèle selon le régime visé : BGK-1 pour les écoulements laminaires, TRT-2 pour
les régimes transitoires, et RLB-4 pour les configurations turbulentes complexes.

53
Conclusion Générale

Cette étude approfondie des modèles hybrides de Boltzmann sur réseau pour la simulation de la convection
naturelle a permis d’évaluer les performances des formulations BGK-1, TRT-2 et RLB-4 à travers une
gamme exceptionnellement large de nombres de Rayleigh (104 à 1012 ). Notre analyse comparative révèle que
l’efficacité des modèles dépend fondamentalement du régime d’écoulement considéré.
En régime laminaire (Ra ≤ 108 ), les trois modèles présentent une précision comparable pour la prédiction
des champs thermiques et des structures d’écoulement. Le modèle BGK-1 se distingue par sa simplicité
algorithmique et son efficacité computationnelle, constituant ainsi le choix optimal pour ces configurations
stables où les phénomènes physiques restent relativement simples.
Dans la zone de transition (108 < Ra < 1010 ), l’émergence d’instabilités met en évidence la supériorité
du modèle TRT-2. Sa formulation à deux temps de relaxation offre un meilleur équilibre entre stabilité
numérique et précision physique, permettant une représentation plus fidèle des phénomènes non-linéaires
caractéristiques de ce régime. La distinction entre composantes symétriques et antisymétriques de la fonction
de distribution confère à ce modèle une robustesse accrue face aux oscillations numériques.
En régime pleinement turbulent (Ra ≥ 1010 ), le modèle RLB-4 démontre une supériorité incontestable.
Sa formulation régularisée d’ordre élevé maintient une stabilité et une précision remarquables même à des
nombres de Rayleigh extrêmement élevés (1012 ). L’analyse des champs instantanés et moyennés confirme
sa capacité exceptionnelle à capturer les structures multi-échelles de la turbulence thermique, avec une
prédiction plus exacte du nombre de Nusselt et une sensibilité réduite aux paramètres numériques.
Ces résultats établissent une hiérarchisation claire pour l’optimisation des simulations numériques : BGK-
1 pour les écoulements laminaires, TRT-2 pour les régimes transitoires, et RLB-4 pour les configurations
turbulentes complexes. Cette différenciation permet une allocation rationnelle des ressources de calcul selon
les exigences physiques du problème étudié.
Les performances démontrées ouvrent des perspectives prometteuses pour l’application de la méthode de
Boltzmann sur réseau à des problèmes industriels et environnementaux variés, depuis les systèmes de refroi-
dissement jusqu’aux phénomènes géophysiques. L’approche hybride développée constitue ainsi une contri-
bution significative à l’avancement des méthodes numériques en mécanique des fluides thermiques, offrant
aux chercheurs et ingénieurs un cadre méthodologique robuste adapté à la complexité croissante des défis
contemporains en thermohydraulique.

54
Bibliographie

[1] Convection naturelle - Aspects théoriques : dossier complet, Techniques de l’ingénierie.

[2] Bejan, A., Convection Heat Transfer (4e éd.), John Wiley & Sons, Inc., 2013.
[3] Succi, S., The Lattice Boltzmann Method for Fluid Dynamics and Beyond, Oxford University Press,
2018.
ISBN : 978-0198503989.
[4] Bejan, A. & Kraus, A.D., Heat Transfer Handbook, John Wiley & Sons, 2003.
ISBN : 978-0471390152.
[5] J.-F. Sacadura, Initiation aux transferts thermiques, Techniques et documentation, Éditions Lavoisier,
1993.
[6] José-Philippe Pérez et A. M. Romulus, Thermodynamique. Fondements et applications, Paris, Masson,
1993, p. 153.
[7] Kharbouch, Cours de convection (LF), Université Abdelmalek Essaâdi, 2019.
[8] Çengel, Y.A. & Ghajar, A.J., Heat and Mass Transfer : Fundamentals and Applications (6e éd.),
McGraw-Hill Education, 2020.
ISBN : 978-0073398181.
[9] Gebhart, B., Jaluria, Y., Mahajan, R.L., & Sammakia, B., Buoyancy-Induced Flows and Transport,
Routledge, 2018.
ISBN : 978-1138198098.
[10] Bergman, T.L., Lavine, A.S., Incropera, F.P., & DeWitt, D.P., Fundamentals of Heat and Mass Transfer
(8e éd.), John Wiley & Sons, 2017.
ISBN : 978-1119353881.
[11] D. Wolf-Gladrow, Lattice-gas cellular automata and Lattice Boltzmann models, Springer, Berlin, 2000.
[12] M. Neveu, Vérification de la méthode de Boltzmann sur réseau en vue de calculer la perméabilité de
milieux poreux, Thèse de doctorat, Université de Montréal, 2017.
[13] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, S. Goncalo, et E.M. Viggen, The Lattice Boltz-
mann Method : Principles and Practice, Springer, 2017.
[14] José-Philippe Pérez et A. M. Romulus, Thermodynamique. Fondements et applications, Paris, Masson,
1993, p. 153.
[15] M.C. Sukop et D.T. Thorne, Lattice Boltzmann Modeling – An Introduction for Geoscientists and En-
gineers, Springer, Heidelberg, Berlin, New York.
[16] A.A. Mohamad, Lattice Boltzmann Method : Fundamentals and Engineering Applications with Computer
Codes, Springer, London, 2011.
[17] C. Ngô et H. Ngô, Physique statistique (3e édition), Dunod, 2008.
[18] S. Marié, Étude de la Méthode Boltzmann sur Réseau pour les Simulations en Aéroacoustiques, Thèse
de doctorat, Université Pierre et Marie Curie - Paris VI, 2008.
[19] A.A. Mohamad, Lattice Boltzmann Method : Fundamentals and Engineering Applications with Computer
Codes, Springer, London, 2011.
[20] S. Chen et G.D. Doolen, “Lattice Boltzmann Method for Fluid Flows,” Annual Review of Fluid Mecha-
nics, vol. 30, pp. 329–364, 1998.
[21] H. Huang, M. Sukop, et X. Lu, Multiphase Lattice Boltzmann Methods : Theory and Application, John
Wiley & Sons, 2015.
[22] Y. H. Qian, D. d’Humières, et P. Lallemand, “Lattice BGK models for Navier-Stokes equation

55
[23] R. R. Nourgaliev, T. N. Dinh, T. G. Theofanous, et D. Joseph, “The lattice Boltzmann equation method :
theoretical interpretation, numerics and implications,” International Journal of Multiphase Flow, vol.
29, no 1, pp. 117–169, 2003.
[24] J.O. Hardy, O. de Pazzis, et Y. Pomeau, “Molecular dynamics of a classical lattice gas : transport
properties and time correlation functions,” Physical Review A, vol. 13, pp. 1949–1961, 1976.
[25] P. A. Skordos, “Initial and boundary conditions for lattice Boltzmann method,” Physical Review E, vol.
48, pp. 4824–4842, 1993.
[26] J. Pan et R. Thompson, “Gauss-Hermite Quadrature Approximation for Estimation in Generalised
Linear Mixed Models,” Journal Name, vol. 18, pp. 57–78, 2003.
[27] J. Latt et B. Chopard, “Lattice Boltzmann Method with Regularized Pre-collision Distribution Func-
tions,” Mathematics and Computers in Simulation, vol. 72, no. 2–6, pp. 165–168, 2006.
[28] J. Miranda Fuentes, Développement d’un modèle de Boltzmann sur gaz réseau pour l’étude du chan-
gement de phase en présence de convection naturelle et de rayonnement, Thèse de doctorat, INSA de
Lyon, 2013.
[29] A. A. Mohamad, Lattice Boltzmann Method : Fundamentals and Engineering Applications with Com-
puter Codes, 2e édition, Springer Nature, 2019, pp. 25, 31, 48–49.
[30] T. Inamuro, M. Yoshino et F. Ogino, “A non-slip boundary condition for lattice Boltzmann simulations,”
Physics of Fluids, vol. 7, pp. 2928–2930, 1995.
[31] H.-B. Huang, X.-Y. Lu et M.C. Sukop, “Numerical study of lattice Boltzmann methods for a convection-
diffusion equation coupled with Navier-Stokes equations,” Journal of Physics A : Mathematical and
Theoretical, vol. 44, pp. 1–8, 2011.

56