Dans cette communication, nous présentons un schéma numérique hybride pour simuler l'interaction entre la convection naturelle et le rayonnement surfacique dans une cavité rectangulaire différentiellement chauffée. Les parois verticales de la cavité sont adiabatiques, tandis que ses parois horizontales supérieure et inférieure sont respectivement froide et chaude. Les vitesses sont déterminées par la méthode de Boltzmann sur réseau et l'équation de l'énergie est discrétisée à l'aide de la méthode des volumes finis. Les systèmes algébriques obtenus sont résolus par la méthode des gradients conjugués. Seul le rayonnement surfacique est pris en compte et les surfaces radiatives sont supposées grises et isotropes en émission/réflexion

MEZRHAB, Ahmed

NAJIB H.

I-Revues

18ème Congrès Français de Mécanique Grenoble, 27-31 août 2007 
 
Modélisation des transferts thermiques couplés à l’aide de la méthode de 
Boltzmann sur réseau. 
 
A. Mezrhaba,*, H. Najib 
 
a Faculté des Sciences, Département de Physique, Laboratoire de Mécanique & Energétique, 
60000 Oujda, Maroc. 
b
 Ustl/Polytech'lille/Lml UMR 8107 CNRS, Département de Mécanique, F-59655 Villeneuve 
d'Ascq Cedex, France. 
*(auteur correspondant:  mezrhab@fso.ump.ma ) 
 
Résumé : 
 
Dans cette communication, nous présentons un schéma numérique hybride pour simuler l’interaction entre la 
convection naturelle et le rayonnement surfacique dans une cavité rectangulaire différentiellement chauffée. Les 
parois verticales de la cavité sont adiabatiques, tandis que ses parois horizontales supérieure et inférieure sont 
respectivement froide et chaude. Les vitesses sont déterminées par la méthode de Boltzmann sur réseau et 
l’équation de l’énergie est discrétisée à l’aide de la méthode des volumes finis. Les systèmes algébriques obtenus 
sont résolus par la méthode des gradients conjugués. Seul le rayonnement surfacique est pris en compte et les 
surfaces radiatives sont supposées grises et isotropes en émission/réflexion.  
 
Abstract : 
 
In this communication, we present a hybrid numerical scheme to simulate the interaction between the natural 
convection and the surface radiation in a differentially heated square cavity. The vertical walls of the cavity are 
adiabatic, while its horizontal walls top and bottom are cold and hot, respectively. The air velocity is determined 
by the lattice Boltzmann method and the energy equation is discretized by using a finite volume method. The 
resulting systems of discretized equations were solved by an iterative procedure based on a preconditioned 
conjugate gradient method. Only the surface radiation is taken into account and the walls of the enclosure are 
assumed to be assumed to be diffuse-grey.  
 
Mots-clefs :  
 
Transfert thermique, volumes finis, méthode Boltzmann sur réseau. 
 
Nomenclature  
 
A rapport de forme de la cavité, L/H 
g accélération de la pesanteur, ms-2 
H hauteur de la cavité, m    
k conductivité thermique, Wm-1K-1  
L longueur de la cavité, m 
Nu nombre de Nusselt moyen 
Nr      nombre de rayonnement, σTc4/(k∆T/H) 
Pr nombre de Prandtl, ν /α   
qr      densité de flux radiatif, Wm-2 
Qr      densité de flux radiatif   
adimensionnelle, qr /σTc4 
Ra nombre de Rayleigh, gβ(Tc-Tf)H3/να 
T température, K 
T0 température moyenne, (Tc+Tf)/2, K 
Tf température de la paroi froide, K 
Tc température du cylindre chaud, K 
u, w composantes de la vitesse, ms-1 
U, V vitesses adimensionnelles selon x, y,  
  U=uH/α, W =wH/α  
x, y coordonnées cartésiennes, m 
X,Y coordonnées cartésiennes adimensionnelles,   
X = x/H, Y = y/H 
 
Symboles grecs 
α    diffusivité thermique, m2s-1 
β     constante de Boussinesq, K-1  
∆T  différence de température maximale 
                  , (Tc-Tf), K 
ν    viscosité cinématique, m2s-1 
ρ    densité du fluide, Kg.m-3     
θ     température adimensionnelle, 
       (T-T0)/(Tc-Tf) 
σ   constante de Stefan-Boltzmann,  
      Wm-2K-4 
 
  Indices 
  c     chaud      f    froid
18ème Congrès Français de Mécanique Grenoble, 27-31 août 2007 
 
 
1.  Introduction 
 
Les écoulements de convection naturelle en cavité bidimensionnelle ont été largement étudiés. Cet 
intérêt est dû à leurs applications industrielles dans divers problèmes d'ingénierie comme l'énergie 
solaire, la thermique de l’habitat et le refroidissement des composants électroniques. Habituellement, 
les méthodes numériques utilisées pour effectuer ces simulations reposent sur la discrétisation spatiale 
et temporelle des équations d’évolution macroscopique. Cependant, leurs solutions peuvent être très 
complexes lorsqu’il s’agit de traiter des problèmes où la géométrie du milieu est complexe ou 
lorsqu’on se trouve en présence de plusieurs phases d’un fluide ou de plusieurs fluides. C’est la raison 
pour laquelle une approche différente dite gaz sur réseau ‘‘Lattice Gas Automata’’ a été développée. 
Son but est de modéliser l’écoulement d’un fluide au niveau microscopique en termes d’interactions 
locales entre les particules. Cette méthode possède beaucoup d’avantages : il est plus facile et intuitif 
de traiter les conditions particulières comme la présence d’obstacles ou l’étude des écoulements 
multiphasiques. De plus, la répétition successive de collisions et de propagations permet une 
implémentation facile sur l’ordinateur. 
 Dans ce papier, nous avons analysé numériquement le transfert thermique par convection naturelle 
et par rayonnement dans une cavité carrée chauffée par le bas à l’aide de la méthode de Boltzmann sur 
réseau couplée à la méthode des volumes finis. 
 
 
2.  Formulation mathématique et modèle numérique 
 
 La cavité étudiée est présentée sur la figure 1. Les parois latérales sont adiabatiques, la paroi 
horizontale supérieure est maintenue à une température froide Tf, alors que la paroi horizontale 
inférieure est portée à une température chaude Tc. On suppose que la géométrie est bidimensionnelle, 
l’écoulement est laminaire et les propriétés physiques de l’air, à part sa densité, sont constantes et 
prises à la température moyenne To. L’approximation de Boussinesq est supposée vérifiée.          
La méthode numérique adoptée est celle des gaz sur réseaux : Gaz de Boltzmann sur réseau qui fait 
appel à des principes entièrement différents, et qui relève d'une description à échelle quasi-
microscopique du fluide. La méthode des "gaz sur réseau" suit le déplacement de milliers de 
"particules de fluide" qui sont contraintes à se déplacer sur un réseau  (qui peut être hexagonal, carré 
avec diagonales ou autre). Les vitesses sont déterminées par la méthode de Boltzmann sur réseau [1].     
 
                                                                                                    
0 1
2
3
4
56
7 8
e1
e2
e3
e4
e5e6
e7 e8
 
 
 
                                                                                         Figure 2 : Modèle D2Q9 
 
2.1. Champs de vitesses  
 
La discrétisation du domaine est faite suivant le modèle D2Q9 basé sur un réseau carré. Chaque 
motif du réseau est caractérisé par neuf vitesses discrètes (ci, (i = 1… , 9)) représentées sur la figure 2. 
Le système évolue par pas alternatifs de collision et de propagation suivant l’équation : 
 ),(),()1,( trftrftcrf i
rrrr Ω+=++   (1) 
paroi adiabatique 
 paroi chaude 
paroi adiabatique 
 y 
 
x 
L 
H 
0 
   
                                      g 
 
Figure 1. Géométrie étudiée 
  paroi froide 
18ème Congrès Français de Mécanique Grenoble, 27-31 août 2007 
 
où fi est une fonction de distribution d’une particule dans la direction i, rr  est la position de la particule 
à l’instant t et Ω est l’opérateur de collision. L’état du système peut être représenté par le vecteur 
Ttrftrftrf )),(),...,,((),( 91
rrr
=  défini dans un espace des phases (F = R9). Une autre idée, 
exprimée par  D’Humière [2],
 
consiste à étendre le vecteur f  en un vecteur m  appartenant à un 
espace des moments (M = R9). Certains de ces moments, comme la densité de masse ρ et les 
composantes de la quantité de mouvement ρu et ρw, se conservent après la collision. La présence de la 
force volumique gβ∆T ),( trrθ  modifie la conservation de la vitesse verticale dans l’étape de collision, 
ce qui se traduit par :  
 
 wac = wbc + gβ∆T ),( trrθ  (2) 
où wac et wbc représentent la composante verticale de la vitesse respectivement après et avant la 
collision.  
Les conditions aux limites dynamiques sur les parois de la cavité sont de types ‘rebonds’, ce qui 
permet d’obtenir u = w =0 au niveau de toutes les parois de la cavité. 
 
2.2.  Champs de température  
 
Les équations de bilan d’énergie  et de concentration sont données par :  
 








∂
∂
+
∂
∂
=
∂
∂
+
∂
∂
2
2
2
2
WU
YXYX
θ θθθ
 (3)  
Les conditions aux limites thermiques utilisées sont : 
 θ =  0.5 pour Y = 0 et  0 ≤ X ≤ A 
θ =  - 0.5 pour Y = 1 et  0 ≤ X ≤ A 
Au niveau des parois adiabatiques : X = 0 ou 1 et 0 ≤ Y ≤ 1 
                     0=−
∂
∂
rNrQY
θ
              (4) 
où Qr la densité adimensionnelle du flux radiatif le long de cette surface. 
 
 Une étude du maillage a été réalisée et le maillage optimal qui permet d’avoir un bon compromis 
(précision/temps de calcul) est de 101x101. 
 L’équation de l’énergie (Eq. 3) a  été discrétisée par une méthode de volumes finis avec un schéma 
de discrétisation centré pour les termes de transport et les systèmes obtenus ont été résolus par la 
méthode des gradients conjugués. 
 Afin de discrétiser l’équation (4), les surfaces radiatives du solide formant la cavité ont été 
discrétisées en un nombre N de surfaces radiatives. N est le nombre total des surfaces radiatives, qui 
est égal au nombre total des interfaces entre volumes de contrôles solide et fluide. Par conséquent, la 
densité adimensionnelle Qri de flux radiatif perdu par la surface Si est donnée par :  
∑
=
−−=
N
j
jijiir FRRQ
1
,           (5) 
et la radiosité adimensionnelle est obtenu en résolvant le système: 
4
1
))1(( iij
N
j
jiiij RF Θ=−−∑
=
− εεδ        (6) 
où Fi-j est le facteur de forme entre les surfaces Si et Sj. Fi-j est déterminé en utilisant une 
approximation par éléments de frontière, en effectuant l’intégration par méthodes de Monte Carlo [3].  
 Pour traiter le couplage entre la convection et le rayonnement, les températures de surfaces des 
parois adiabatiques sont calculées en résolvant l’équation non linéaire (Eq. 4). Cette résolution est 
effectuée en utilisant une procédure itérative interne avec sous relaxation à chaque pas de temps. A 
18ème Congrès Français de Mécanique Grenoble, 27-31 août 2007 
 
chaque itération interne (relative au calcul des températures de surfaces), le système linéaire des 
équations de radiosités (Eq. 6) est résolu par une méthode directe. 
 Enfin, pour caractériser les transferts de chaleur, le nombre de Nusselt moyen est donné par : 
 dXYXQN
YA
Num
A
wrr
YX w
∫ 






=+
∂
∂
−=
=
0
0,
)0,(1 θ  (7) 
 
3. Résultats et discussions 
 
Le code a été intensivement exercé sur des problèmes de référence pour vérifier sa validité. 
Nous rappelons dans le tableau ci-dessous quelques résultats obtenus, en convection naturelle, par 
notre code en comparaison avec ceux rapportés dans les références [4,5]  
 
Ra 
   103      104      105        106        107 
Présent  1.118      2.242     4.524      8.824       16.490 
Réf. [4,5]  1.118      2.243     4.523      8.826       16.510 
 
Quand l'échange de rayonnement est pris en considération, notre code a été validé avec les 
résultats numériques publiés antérieurement et rapportés par différents auteurs, et l'accord entre le 
présent et les résultats précédents était très bon dans la référence [6]. Pour cette raison, et pour la 
brièveté on ne le répète pas ici. Basé sur les études ci-dessus, on a conclu que le code pourrait être 
correctement appliqué au problème considéré.  
En convection naturelle pure, chaque configuration dépend au minimum de 3 paramètres 
adimensionnels (A, Pr, Ra).  La prise en compte du rayonnement thermique nécessite la donnée de 
quatre autres paramètres (T0, ∆Τ, L, ε). Pr, A et T0 sont fixés respectivement à 0.71, 1 et   300 K. Dans 
ces conditions, le terme gβ/να est de l’ordre de 9×107 m-3K-1 et le nombre de Rayleigh, variant de 103 
à 106, n’est fonction que de ∆Τ et L. L’émissivité ε des surfaces radiatives est prise égale à 0 en 
convection naturelle pure et égale à 1 en présence du rayonnement. Dans cette étude, deux cas ont été 
considérés : dans un premier temps, nous avons fixé ∆T à 20 K  et nous avons calculé L  pour chaque 
valeur de Ra et dans un second temps, nous avons fixé L à 3 cm et ∆T a été calculé pour chaque valeur 
de Ra. Le nombre de rayonnement Nr est calculé, à son tour, en fonction du nombre de Rayleigh selon 
que l’on fixe ∆T ou L. 
La figure 3 présente les isothermes et les lignes de courant en absence et en présence du 
rayonnement pour deux nombres de Rayleigh 104 et 105. 
Dans le cas de la cavité chauffée par la paroi horizontale basse, il a été démontré dans le passé, 
la présence de solutions multiples en régime permanent [7]. Les nombres de Rayleigh considérés dans 
cette étude produisent des écoulements monocellulaires avec des rotations possibles dans le sens des 
aiguilles d’une montre et dans le sens inverse. 
 En mode de convection naturelle combiné au rayonnement thermique (ε = 1), l’inclinaison des 
isothermes dans les régions se trouvant près des parois adiabatiques est due à l’importance des flux 
radiatifs. Les lignes de courant montrent que le rayonnement thermique réduit considérablement 
l’écoulement dans la cavité. En fait, le rayonnement réduit l’écart de températures entre les parois 
adiabatiques et donc, diminue la vitesse de l’air au voisinage de ces parois. 
 Le nombre de Nusselt moyen est présenté sur la figure 4. On observe que le nombre de Nusselt 
moyen augmente sous l’influence des forces de poussée (augmentation de Ra) et du rayonnement 
thermique. Cela s’explique par le fait que le nombre de Nusselt moyen (Eq. 7) est proportionnel au 
produit NrQr, qui augmente avec Ra, que ce soit lorsqu’on fixe ∆T ou L. 
 
 
 
18ème Congrès Français de Mécanique Grenoble, 27-31 août 2007 
 
4. Conclusion 
 
Le couplage convection naturelle rayonnement dans une cavité carrée différentiellement 
chauffée, dont les parois verticales sont adiabatiques, a été étudié numériquement. La méthodologie 
numérique utilisée est un schéma hybride basé sur la méthode de Boltzmann sur réseau combinée à la 
méthode des volumes finis. 
A la lumière de cette étude, les conclusions suivantes ont été formulées :   
1)  Le couplage de la méthode de Boltzmann sur réseau à la méthode des volumes finis donne de très 
bons résultats, 
2)  Le rayonnement thermique homogénéise la température à l’intérieur de la cavité en réduisant 
l’écart de température entre les parois adiabatiques. Il réduit la circulation d’air dans la cavité. 
3)  Le rayonnement thermique augmente le transfert de chaleur dans la cavité. 
 
Références  
 
 
[1]  A. Mezrhab, M. Bouzidi and P. Lallemand, “Hybrid lattice-Boltzmann finite-difference simulation 
of convective flows”, Computers & Fluids, vol. 33, pp. 623-641, 2004 
[2] D. d’Humières, Generalized lattice-Boltzmann equations, in : B.D. Shizgal, D.P. Weaver (Eds.), 
Rarefied Gas Dynamics : Theory and Simulations, Prog. Astronaut. Aeronaut., 159, AIAA, 
Washington, DC, (1992), 450-58. 
[3] A. Mezrhab and M. Bouzidi, “Computation of view factors for surfaces of complex shape including 
screening effects and using a boundary element ”, Engineering Computations: International Journal for 
Computer-Aided Engineering and software, vol. 22, n° 2, pp. 132-148, 2005 
[4] G. De Vahl Davis, “Natural convection in air in a square cavity: A Bench Mark Numerical 
Solution”, Int. J. Numer. Meth. Fluids, vol. 3, pp. 249-264, 1983 
[5] P. Le Quéré , “Accurate solutions to the square thermally driven cavity at high Rayleigh number”, 
Computers & Fluids, vol. 20, pp. 29-41, 1991 
[6] A. Mezrhab, M. Jami, M. Bouzidi, P. Lallemand, “Analysis of radiation–natural convection in a 
divided enclosure using the lattice Boltzmann method”, Computers & Fluids, Vol. 36, n°, pp. 423-434 
[7] M. Hasnaoui, E. Bilgen, P. Vasseur, “Natural convection Heat Transfer in rectangular cavities 
partially heated from below”, J. Thermophys. Heat Transfer, vol. 6, n°2, pp. 255-264, 
1992 
 
 
                Ra =  104                          Ra =   105         
 
 
        
                             ε = 0 
 
        
                           ε = 0 
 
         
                             ε = 1, L = 3 cm 
 
         
                             ε = 1, L = 3 cm 
  
  
          
                             ε = 1, ∆T = 20 K  
         
                             ε = 1, ∆T = 20 K  
                                                      
 
Figure 3: Isothermes et Lignes de courant en absence (ε = 0) et  en présence (ε = 1)  du rayonnement 
pour Ra = 104 et 105. 
 
103 104 105
2
4
6
8
10
Ra
 ε = 0
 ε = 1, L = 3 cm
 ε = 1, ∆T = 20 K
Nu
Ra
  
Figure 4 : Variation du nombre de Nusselt moyen en fonction de Ra.  


Modélisation des transferts thermiques couplés à l'aide de la méthode de Boltzmann sur réseau

Abstract

Similar works

Full text

Available Versions

I-Revues