# Simulation et Visualisation de la Couche Limite Thermique au-Dessous d'un Disque Horizontal

Ahmed El motassadeq, Hassan Chehouani, Mohamed Waqif, Sauveur Benet*

LPEM, Laboratoire des Procédés Energétiques et Métrologie, BP 618 Faculté des Sciences et Techniques, Université Cadi Ayyad, Gueliz, Marrakech, Maroc.

**Résumé** – Ce travail concerne la simulation numérique et la visualisation de la couche limite thermique au voisinage d'un disque horizontal en convection libre dans l'air. Le modèle que nous avons développé est fondé sur la résolution des équations de Navier-Stokes couplées à l'équation de l'énergie pour décrire l'écoulement et le transfert de chaleur. Nous avons utilisé la technique des volumes finis pour résoudre numériquement les équations. On simule ainsi le champ de température et les lignes de courant dans le milieu. Les résultats mettent en évidence la présence de forts gradients de température entraînant l'établissement d'une mince couche limite thermique localisée près du disque. La visualisation de cette couche a été effectuée par la technique optique d'ombroscopie. Une approche numérique a été suggérée pour expliquer la caustique de la couche limite thermique qui joue le rôle d'une lentille thermique.

**Abstract** – This work deals with a numerical simulation and a visualisation of the free convection thermal boundary layer near a heated horizontal disk in ambient air. The developed model is based on the resolution of the coupled Navier-Stokes 's equations and the energy equation describing the flow patterns and the heat transfer. The model equations are solved numerically using a finite control volumes technique leading to a fine description of the temperature field and the streamlines near the heated disk. The results show strong gradients of temperature and velocity, which give rise to the creation of a thin thermal boundary layer localised in the neighbouring of the disk. The shadowgraph technique is used to visualise this layer. A numerical approach is suggested to explain the caustic of the thermal boundary layer, which plays the role of a thermal lens.

**Mots clés** : Transfert de chaleur – Couche limite thermique – Lentille thermique – Ombroscopie – Simulation numérique.

## 1. INTRODUCTION

Dans le domaine de la modélisation des transferts thermiques entre une paroi et un fluide, la connaissance des caractéristiques de la couche limite thermique et de la répartition de température est une étape fondamentale afin de tester la validation des modèles numériques. Cependant, une telle expérimentation n'est pas simple en général. En effet, il s'agit de mesurer des grandeurs locales dans un écoulement convectif sensible aux perturbations extérieures. Cette contrainte exclut les techniques de contact consistant à introduire un capteur thermique dans le fluide et favorise l'emploi des méthodes de contrôle non destructif (NDT). Dans ce contexte, les techniques optiques constituent les outils les plus utilisés en mécanique des fluides anisotropes [1, 2, 3]. Leur avantage majeur réside dans le fait qu'elles rendent certaines propriétés du champ d'écoulement directement accessibles à la perception visuelle sans aucune perturbation du milieu. Le principe de cette métrologie consiste à quantifier les modifications que subit une onde électromagnétique lors de sa traversée de l'écoulement soumis à l'observation.

<footnote>(*) Laboratoire de Physique Appliquée et Automatique, Université de Perpignan, Av. de Villeneuve, F66825 Perpignan Cedex, France.</footnote>
Les techniques d'interférométrie, qui mettent en évidence les variations de la phase de l'onde, ne sont pas directes car on doit effectuer tout un ensemble de traitements numériques complexes qui conditionnent la précision de la mesure [4,5,6]. Chehouani a utilisé l'interférométrie holographique pour étudier les transferts de chaleur autour d'un substrat cylindrique placé dans un réacteur CVD a paroi froide [6]. Les variations de température et donc d'indice de réfraction dans le milieu sont converties en un réseau de franges d'interférence qui permettent de visualiser la couche limite thermique. Cependant, l'auteur met en évidence la difficulté de restituer le champ de température à partir de l'interférogramme qui nécessite l'inversion de l'intégrale d'Abel dont la précision est extrêmement sensible aux erreurs expérimentales. Sabri et al. [7] ont également souligné la complexité du traitement numérique des interférogrammes réalisés pour mesurer les gradients thermique et massique au voisinage d'une plaque inclinée en transpiration.

En comparaison avec l'interférométrie, la technique d'ombroscopie qui exploite la déviation de la trajectoire des rayons lumineux par un gradient d'indice de réfraction, constitue le moyen le plus simple à mettre en œuvre pour visualiser la couche limite thermique [3]. Cette méthode considérée jusqu'à maintenant comme étant qualitative, peut être améliorée en exploitant quantitativement l'ensemble de l'information contenue dans l'image du milieu étudié.

L'objectif de ce travail est de présenter une première approche de la simulation de la trajectoire d'un faisceau laser traversant la couche limite thermique établie au-dessous d'un disque horizontal maintenu à température constante.

## 2. ETUDE EXPERIMENTALE

### 2.1. Dispositif expérimental

L'étude a été réalisée sur un disque à axe vertical chauffé dans l'air ambiant (figure 1).

<figure><img src="image_1.png" /><figcaption>Fig. 1. Schéma de principe du montage optique de la technique de l'ombroscopie</figcaption></figure>
C'est un cylindre en aluminium de rayon ($R_D$ = 19.75 mm), dans lequel est noyée une résistance électrique de telle manière à le porter à une température uniforme de 420 °C. Les courants de convection libre induisent l'établissement d'une couche d'air chaud au voisinage du cylindre. Les gradients de température dans cette couche entraînent des variations de la masse volumique et par conséquent de l'indice de réfraction de l'air. Un faisceau parallèle issu d'une source laser He-Ne de puissance 1 mW traverse le milieu et subit alors des modifications de son chemin optique par l'effet de la réfraction. Un écran placé en aval des rayons émergeants permet de visualiser et de déterminer les caractéristiques de la couche limite qui joue le rôle d'une lentille thermique [3].

## 2.2. Résultats expérimentaux

Les figures 2 et 3 représentent l'ombre portée du disque sur l'écran placé respectivement à deux distances différentes ($l$ = 1.71 m et $l$ = 2.48 m). Ces photographies illustrent la concentration de la lumière sous la forme d'une mince bande lumineuse entourant le disque. Ce phénomène optique est dû aux effets de la réfraction sur les rayons lumineux et permet de matérialiser la frontière de la couche limite thermique. Il faut signaler que cette bande n'apparaît d'une façon nette que pour une position de l'écran dans un intervalle déterminé. On peut faire une analogie de ce système avec une lentille astigmatique dont l'indice de réfraction est plus faible que le fluide environnant. L'enveloppe des rayons traversant cette couche limite forme la caustique de la lentille thermique.

<figure><img src="image_2.png" /><figcaption>Fig. 2. Photographie de l'objet porté à la température de 420 °C, $l$ = 1.71 m</figcaption></figure>

En outre, les photos des figures 2 et 3 révèlent la présence des franges d'interférence qui sont très serrées et localisées. Elles sont le résultat des intersections des rayons lumineux émergeants du milieu avec des différences de marche. En comparaison avec des travaux publiés dans ce domaine [2, 3], l'utilisation d'une source cohérente a rendu ce phénomène d'interférence visible avec un bon contraste.
<figure><img src="image_3.png" /><figcaption>Fig. 3. Photographie de l'objet porté à la température de 420 °C, $l = 2.48$ m</figcaption></figure>

## 3. SIMULATION NUMERIQUE

### 3.1. Equations du modèle

Pour comprendre les phénomènes optiques concernant cette lentille thermique et expliquer l'existence de cette caustique, il est nécessaire de connaître la distribution de l'indice de réfraction dans le milieu. Cette grandeur est liée à la température de l'air par la relation de Gladstone-Dale [1].

La première étape de cette approche consiste à simuler les écoulements et le champ de température dans l'air voisinant le disque chaud. Nous avons mis au point un modèle numérique fondé sur l'écriture des équations de Navier-Stokes et du transfert de chaleur avec les conditions aux limites imposées par la géométrie du système. Dans les hypothèses d'un écoulement axisymétrique, laminaire, incompressible en régime établi, ces équations s'écrivent en coordonnées cylindriques (r, z) de la manière suivante :

* **Equation de continuité**

$$ \frac{1}{r} \frac{\partial}{\partial r} (r \rho V_r) + \frac{\partial}{\partial z} (\rho V_z) = 0 \qquad (1) $$

* **Equation de mouvement**

• Composante selon r :

$$ \rho \left( V_r \frac{\partial V_r}{\partial r} + V_z \frac{\partial V_r}{\partial z} \right) = -\frac{\partial P}{\partial r} + \frac{1}{r} \frac{\partial}{\partial r} \left\{ r \mu \left[ 2 \frac{\partial V_r}{\partial r} - \frac{2}{3} \left( \frac{1}{r} \frac{\partial}{\partial r} (r V_r) + \frac{\partial V_z}{\partial z} \right) \right] \right\} \\ - \frac{\mu}{r} \left[ 2 \frac{V_r}{r} - \frac{2}{3} \left( \frac{1}{r} \frac{\partial}{\partial r} (r V_r) + \frac{\partial V_z}{\partial z} \right) \right] + \frac{\partial}{\partial z} \left[ \mu \left( \frac{\partial V_z}{\partial r} + \frac{\partial V_r}{\partial z} \right) \right] \qquad (2) $$
• Composante selon z :

$$ \rho \left( V_r \frac{\partial V_z}{\partial r} + V_z \frac{\partial V_r}{\partial z} \right) = - \frac{\partial P}{\partial z} + \frac{1}{r} \frac{\partial}{\partial r} \left[ r \mu \left( \frac{\partial V_z}{\partial r} + \frac{\partial V_r}{\partial z} \right) \right] + \frac{\partial}{\partial z} \left\{ \mu \left[ 2 \frac{\partial V_z}{\partial z} - \frac{2}{3} \left( \frac{1}{r} \frac{\partial}{\partial r} (rV_r) + \frac{\partial V_z}{\partial z} \right) \right] \right\} - \rho g \quad (3) $$

Les propriétés physiques de l'air sont considérées comme variables avec la température selon les expressions suivantes facilitant le calcul numérique :

Masse volumique : $\rho(T)=\frac{PM}{RT}$

Viscosité dynamique : $\mu(T)=1.875 \times 10^{-5} \left(\frac{T}{293}\right)^{0.6386}$

Conductivité thermique : $\lambda(T)=0.026 \left(\frac{T}{293}\right)^{0.7231}$

Capacité calorifique à pression constante : $C_p(T)=975.3+0.0368T+2 \times 10^{-4} T^2$ où T est la température absolue.

<figure><img src="image_4.png" /><figcaption>Fig. 4 : Maillage du domaine de calcul OABCDE comportant 83 x 109 nœuds respectivement selon les axes or et oz (OA=150mm, AB=300mm, BC=130.25mm, CD=30mm, DE=$R_D=19.75mm$, EO=270mm)</figcaption></figure>
Les conditions aux limites sont fondées sur la condition d'imperméabilité et de non glissement du fluide au niveau de la paroi du disque. L'hypothèse d'axisymétrie évoquée précédemment impose la parité des fonctions à calculer et permet de restreindre le domaine de calcul à la moitié du disque et son voisinage. Au niveau de la frontière du domaine de calcul, il est difficile de postuler les conditions aérodynamiques en convection libre. Pour remédier à cette indétermination, nous avons considéré un domaine de calcul suffisamment grand pour pouvoir supposer que l'air loin du disque est immobile.

La résolution numérique a été réalisée en utilisant un code de calcul aux volumes finis [6]. Le domaine de calcul (figure 4) est discrétisé selon un maillage rectangulaire à pas variable afin de s'affranchir des forts gradients de température [6, 8] en particulier au voisinage de l'échantillon (le pas est de l'ordre de 0.5 mm selon les deux directions).

Les équations du modèle sont mises selon une seule forme d'équation de conservation permettant d'effectuer une résolution commune des équations de transfert par intégration sur chaque volume de contrôle [9, 10]. Les systèmes d'équations algébriques résultant du schéma numérique sont résolus itérativement par la méthode de Gauss-Seidel [11]. Les calculs sont effectués sur un micro-ordinateur de type Pentium II et le critère de convergence est satisfait après 3000 à 5000 itérations (selon le nombre de nœuds composant le maillage) avec une durée de l'ordre de 20 minutes [8].

## 3.2. Résultats de la simulation

### 3.2.1 Ecoulement et champ de température

Le tracé des lignes de courant de la figure 5 montre la topologie de l'écoulement. La simulation a été réalisée dans les mêmes conditions expérimentales des figures 2 et 3 ($T_D = 420$ °C, $P = 1.013 \ 10^5$ Pa d'air et $R_D = 19.75$ mm). On constate que les lignes de courant se resserrent au voisinage de l'échantillon pour former une couche limite hydrodynamique mettant ainsi en évidence l'existence d'un fort gradient de vitesse dû aux effets de la convection libre.

<figure><img src="image_5.png" /><figcaption>Fig. 5 : Les lignes de courant de l'écoulement autour du disque à 420°C</figcaption></figure>
La figure 6 représente les isothermes de 40 °C à 400 °C avec un pas de 80°C. Ces valeurs sont localisées très près de l'échantillon illustrant la présence d'un fort gradient de température de l'ordre de 100 °C/mm. On peut considérer en première approximation que

<figure><img src="image_6.png" /><figcaption>Fig. 6 : Les isothermes autour du disque à 420°C de 40 à 400°C avec un pas de 80 °C</figcaption></figure>

l'isotherme 40 °C représente la frontière de la couche limite thermique dont l'épaisseur au centre est de l'ordre de 5 mm (figure7).

<figure><img src="image_7.png" /><figcaption>Fig. 7 : Schématisation de la couche limite thermique établie sous le disque</figcaption></figure>

Le tracé des distributions de température respectivement radiale et axiale (figures 8 et 9) permet d'évaluer le gradient thermique. Dans ces graphiques, l'origine de l'axe z est prise au centre de la face inférieure du disque. On constate que le profil radial est sensiblement constant sous une grande partie du disque (≈ 80%) et décroît rapidement en s'approchant du bord sur une distance de l'ordre de $\delta \approx 5$mm. Le profile axial T(r=0, z) illustre également une
forte décroissance de la température de 400 °C à 20 °C sur une distance voisine de 5mm. Ce résultat est en bon accord avec les observations et les mesures réalisées par Chehouani et al. en utilisant l'interférométrie holographique [12].

<figure><img src="image_8.png" /><figcaption>Fig. 8 : Profil radial de la température sous le disque pour plusieurs cotes (z=cst)</figcaption></figure>

Etant donnée l'importance de ces gradients de température et donc d'indice de réfraction engendrés par les effets de le convection naturelle, on peut s'attendre à une influence non négligeable de ce phénomène sur la trajectoire des rayons lumineux traversant la couche limite sous le disque.

<figure><img src="image_9.png" /><figcaption>Fig. 9 : Profil axial de la température sous le disque pour (r = 0)</figcaption></figure>
### 3.2.2 Trajectoire des rayons

A partir du champ de température ainsi simulé on peut déduire la distribution de l'indice de réfraction de l'air en appliquant l'équation de Gladstone-Dale (figure 10):

$$n(r, z) = 1 + \frac{KPM}{RT(r, z)} \qquad (4)$$

$$ \text{avec} \quad \begin{array}{l} K = 0.2256 \ 10^{-3} \text{ m}^3/\text{kg}, \text{ constante de Gladstone-Dale} \\ M = 28.810 \ 10^{-3} \text{ kg/mole}, \text{ masse molaire de l'air} \\ P = 1.013 \ 10^5 \text{ Pa} \qquad R = 8.314 \text{ J/mole/K} \end{array} $$

<figure><img src="image_10.png" /><figcaption>Fig. 10 : Distribution radiale de l'indice de réfraction sous le disque pour plusieurs cotes ($z = CST$)</figcaption></figure>

Les courbes de la distribution radiale de l'indice de réfraction $n(r, z)$, ont un comportement analogue à celui de la température $T(r, z)$. Les variations d'indice de réfraction sont très accentuées aux extrémités du disque. On peut s'attendre à ce que les effets de bord soient les plus significatifs sur la déviation des rayons lumineux par la réfraction.

Ainsi pour expliquer qualitativement les observations expérimentales représentées par les clichés des figures 2 et 3, nous avons tracé le profil simulé de la masse volumique et son influence qualitative sur le trajet des rayons lumineux (figure 11). Il faut noter que contrairement au profil axial de température (figure 9), la masse volumique présente un point d'inflexion pour $z = 2.5$ mm. L'éclairement sur l'écran présente alors un maximum donnant lieu à une bande lumineuse définissant l'image de la couche limite thermique.

Si l'on reprend le schéma de la figure 1, les rayons lumineux parallèles formés par la lentille $L_1$ traversent le milieu d'indice variable $n(r, z)$, représenté par la figure 10, et subissent une certaine déviation. La trajectoire de chaque rayon peut être déterminée par la résolution de l'équation de Snell [1, 3] :
$$ \frac{d}{ds} \left[ n \frac{d\bar{s}}{ds} \right] = \bar{\nabla} \cdot \mathbf{n} \qquad (5) $$

<figure><img src="image_11.png" /><figcaption>Fig. 11 : Interprétation qualitative du comportement optique de la couche limite thermique</figcaption></figure>

$\bar{s}$ est le vecteur position d'un point appartenant à un rayon et $s$ désigne son abscisse curviligne (figure 12).

<figure><img src="image_12.png" /><figcaption>Fig. 12 : Schéma de la déviation des rayons lumineux réfractés dans un milieu à gradient d'indice</figcaption></figure>

Cette équation vectorielle peut être décomposée en deux équations différentielles indépendantes [3] :
$$ \frac{\frac{d^2 r}{dx^2}}{1 + \left[ \frac{dr}{dx} \right]^2} = \frac{1}{n} \frac{\partial n(r, z)}{\partial r} \qquad (6) $$

$$ \frac{\frac{d^2 z}{dx^2}}{1 + \left[ \frac{dz}{dx} \right]^2} = \frac{1}{n} \frac{\partial n(r, z)}{\partial z} \qquad (7) $$

où x désigne la direction de propagation du faisceau incident.

Comme l'indice de réfraction n'est connu que d'une façon discrète sur le maillage ayant servi à la simulation des champs de température, les deux équations sont résolues numériquement.

Selon les travaux de Vest et de Lira [13, 14, 15], la solution de l'équation (6), nous donne la déviation radiale de chaque rayon lumineux traversant la couche limite thermique, sa trajectoire qualitative est illustrée sur la figure 13/b. En outre, la contribution du gradient axial de l'indice (équation 7), conduira à des chemins optiques analogues à ceux représentés par la figure 13/a.

<figure><img src="image_13.png" /><figcaption>Figure 13. Schéma de la trajectoire des rayons lumineux réfractés sous le disque a/ - gradient d'indice axial- b/ -gradient d'indice radial-</figcaption></figure>
## 4. CONCLUSION

Ce travail met l'accent sur les avantages qu'offre la technique d'ombroscopie pour visualiser la couche limite thermique et la simplicité de mise en oeuvre de cette méthode. La comparaison des résultats de la simulation avec les observations expérimentales a permis de valider notre code de calcul de transfert de chaleur. La simulation numérique a conduit au calcul de la distribution de l'indice de réfraction nécessaire à la résolution de l'équation qui régit la trajectoire des rayons lumineux.

Un travail est en cours pour simuler la caustique de la lentille thermique a partir de l'équation de Snell. Ceci nous permettra à priori d'interpréter d'une façon plus quantitative le phénomène d'interférence observé expérimentalement et de corréler l'interfrange au gradient de température.

## NOMENCLATURE

<table><tr><td>g</td><td>Accélération de la pesanteur, m.s<sup>-2</sup></td><td>R<sub>D</sub></td><td>Rayon du disque = 19.75 mm</td></tr><tr><td>K</td><td>Constante de Gladstone-Dale, m<sup>3</sup>/kg</td><td>s</td><td>Abscisse curviligne</td></tr><tr><td>l</td><td>Position de l'écran d'observation, m</td><td>s&#773;</td><td>Vecteur position</td></tr><tr><td>M</td><td>Masse molaire, kg.mole<sup>-1</sup></td><td>T</td><td>Température du gaz, K</td></tr><tr><td>n</td><td>Indice de réfraction du milieu</td><td>T<sub>D</sub></td><td>Température du disque = 420 °C</td></tr><tr><td>P</td><td>Pression, Pa</td><td>V<sub>r</sub></td><td>Vitesse radiale, m.s<sup>-1</sup></td></tr><tr><td>r, z</td><td>Coordonnées cylindriques, m</td><td>V<sub>z</sub></td><td>Vitesse radiale, m.s<sup>-1</sup></td></tr><tr><td>δ</td><td>Epaisseur de la couche limite, mm</td><td>μ</td><td>Viscosité dynamique, kg.m<sup>-1</sup>.s<sup>-1</sup></td></tr><tr><td>λ</td><td>Conductivité thermique, W.m<sup>-1</sup>.K<sup>-1</sup></td><td>ρ</td><td>Masse volumique, kg.m<sup>-3</sup></td></tr><tr><td>C<sub>P</sub></td><td>Chaleur spécifique à pression constante, J.kg<sup>-1</sup>.K<sup>-1</sup></td><td></td><td></td></tr><tr><td>R</td><td>Constante des gaz parfaits = 8.314 J.mole<sup>-1</sup>.K<sup>-1</sup></td><td></td><td></td></tr></table>

## REFERENCE

[1] C.H. Vest, *Holographic Interferometry*, J. WILEY & SON, New-York (1978) T.F. Irvine, Jr. Ed. Academic Press, New-York (1970)

[2] W. Merzkirch *Flow Visualization* Academic Press, New-York (1974)

[3] W. Hauf, U. Grigull, "Opticals Methods in Heat Transfer" in *Advances in Heat Transfer*, Vol 6, J.P Hartnett, (1970).

[4] T. Kreis, J. Opt. Soc. Am., Vol 3, N°6 (1986).

[5] D. A. Murio, *The Mollification Method and the Numerical Solution of Ill-Posed Problems*, John Wiley & Sons, New-York, (1992)

[6] H. Chehouani "Etude Théorique et Expérimentale des Phénomènes de Transport dans les Réacteurs CVD" Thèse de Doctorat d'état, Marrakech, Maroc (1999)

[7] A. Sabri, S. Benet, R. Berge, H. Chehouani.15ème Congrès Canadien de Mécanique Appliquée, Victoria, pp. 424-425, (1995).

[8] H. Chehouani, A. Belghit, B. Claudet, B. Armas, B. Granier, Entropie, 216, 14-20.(1999).

[9] A.D. Gosman, W.M. Pun, A.K. Runchal, D.B. Spalding and M. Wolfshtein, *Heat and Mass Transfer in Recirculating Flows*, Academic Press, London, (1973)

[10] S.V.PatankaR, *Numerical Heat Transfer and Fluid Flow* McGraw Hill New.York., (1980).
[11] A. Gourdin et M. Boumahrat; *Méthodes Numériques Appliquées*, Ed. Tec. Doc. Lavoisier, Paris (1983).

[12] H. Chehouani, S. Benet, S. Brunet, B Armas, J.Physique IV Coll. C3 , Volume 3, pp. 83-90, (1993).

[13] C.H. Vest, Appl. Opt. **14**,1601-1606. (1975).

[14] I.H. Lira and C. M. Vest, Appl. Opt. **18**, 3919-3928. (1987).

[15] I.H. Lira, "*Correcting for Refraction Effects in Holographic Interferometry of Transparent Objects*", Ph.D Thesis, U.Mchingan (1987).