Refroidissement de composants électroniques à l'aide du logiciel OpenFOAM

 

Refroidissement de composants électroniques

à l'aide du logiciel OpenFOAM


 

 

Hector Alvarado Florent Dizes Astrid Molinier

Elèves ingénieurs dernière année ENSEEIHT Mécanique des Fluides

 

Introduction

Ce Bureau d'Etudes Industriel (BEI), proposé par la société Epsilon et mené à l'ENSEEIHT dans le département Mécanique des Fluides, traite de l'étude du refroidissement de composants électroniques, principalement à l'aide du logiciel libre OpenFOAM.

Cette partie introductive à notre BEI a pour but d'exposer, tout d'abord, le contexte industriel qui lui est rattaché. Dans un second temps, les objectifs de celui-ci seront abordés avec en particulier un bref récapitulatif des moyens industriels utilisés pour refroidir les composants électroniques ainsi que les moyens mis en oeuvre dans ce BEI pour répondre à cette attente. L'entreprise partenaire Epsilon sera ensuite présentée avant de terminer par les aspects plus techniques concernant les logiciels utilisés pour notre étude.

1. Contexte industriel

Grâce aux progrès technologiques réalisés dans le domaine de l'électronique, les composants ont l'avantage de devenir de plus en plus performants tout en étant de plus en plus petits. Dans ce contexte, les densités surfacique et volumique de chaleur à évacuer deviennent plus importantes. Dans de multiples marchés, tels l'Aeronautique, le Spatial, le Naval ou les systèmes informatiques, ces composants électroniques sont très abondants et dégagent une chaleur non négligeable lorsqu'ils sont en fonctionnement, ce qui peut les endommager. En effet, les plages de température des composants ne peuvent excéder un intervalle de -20°C à 40°C en fonctionnement.
 
 
De plus, en ce qui concerne les composants dans les satellites, il y a bien plus de contraintes par rapport aux variations de température à cause de l'exposition de ceux-ci soit à l'ombre, soit en plein soleil. La figure ci-contre, par exemple, illustre le fort contraste de températures qu'il peut exister sur un satellite en orbite. En conséquence, il est nécessaire de refroidir ces composants soumis à une forte chaleur pour un fonctionnement optimal et durable, de même qu'il faut les protéger contre les températures trop froides : c'est le contrôle thermique des satellites.
 
Exemple d'une carte des températures d'un satellite en orbite

 

2. Présentation

Comment chauffe un composant électronique?

Lorsqu'un composant électronique est parcouru par un courant électrique créé par des porteurs de charges, ces derniers interagissent avec les atomes du composant, ce qui induit une résistance au déplacement. Ainsi, pour contrer cette résistance, il faut fournir une puissance suffisante qui est ensuite dissipée sous forme de chaleur lors des interactions entre atomes : c'est l'effet Joule. Cette chaleur doit être évacuée par des méthodes industrielles telles que présentées dans la prochaine section, sous peine de fonte du composant.

Loi de Joule pour un conducteur ohmique

 

Quelles sont les différentes façons de refroidir un composant électronique?

On trouve dans l'industrie trois façons de refroidir les composants électroniques : le refroidissement par liquide, par air ou par convection naturelle.

 

1. Refroidissement par liquide

  • le caloduc

Un caloduc (heat pipe) est un conducteur thermique, bien plus efficace que les métaux, servant à transporter la chaleur d'une source chaude vers une source froide, par le principe de changement de phase d'un fluide caloporteur.

Le principe du caloduc est basé sur la circulation d'un liquide en équilibre avec sa vapeur entre une région chauffée, l'évaporateur, et une région refroidie, le condenseur, le tout en l'absence d'air. Alors que la vapeur se déplace grâce à la différence de pression entre l'évaporateur et le condenseur, le condensat retourne alors vers l'évaporateur sous l'effet de forces. La circulation de ce liquide est induite soit par effet de la gravitation, soit par effet capillaire. Ainsi, selon ces effets, on distingue deux types de caloducs : les thermosiphons diphasiques et les caloducs capillaires.

Le caloduc est un dispositif statique permettant le transfert de flux thermiques très élevés avec un faible gradient thermique. Dès lors, il peut être très utilisé pour le refroidisssement de microprocesseurs d'ordinateurs ou dans le domaine spatial, où existent d'importants flux thermiques à évacuer pour éviter l'altération du satellite. En effet, le contrôle thermique des composants dans ce domaine a été une des premières applications des caloducs qui sont très efficaces et souples dans la gestion thermique.

 

Principe de fonctionnement d'un caloduc
  • la boucle fluide diphasique à pompage mécanique

Ce système est composé d'un ou plusieurs évaporateurs et condenseurs ainsi que d'une pompe permettant la circulation du liquide. Les composants électroniques à refroidir sont placés à proximité de l'évaporateur de la boucle puis le liquide saturé issu de la pompe soutire de la chaleur aux composants. Le liquide devient donc vapeur à partir de la température de saturation, ce qui augmente d'autant plus l'efficacité des transferts thermiques. La vapeur arrive ensuite au condenseur par le biais de la pompe, ce qui permet la liquéfaction de celle-ci.

 

  • la boucle fluide monophasique à pompage mécanique

La boucle fluide monophasique utilise la chaleur du fluide caloporteur. Ce dernier absorbe la puissance dissipée par les équipements puis la rejette en se refroidissant dans un radiateur sans changer de phase. Le fluide est mis en mouvement par un système de pompage mécanique mais qui peut présenter l'inconvénient d'une durée de vie courte, peu souhaitable lors d'une mission satellitaire.

 

  • la boucle fluide diphasique à pompage capillaire

Tel le caloduc, ce type passif de boucle fluide utilise les propriétés de changement de phase d'un fluide caloporteur, mais permet le transport de chaleur sur une plus grande distance pouvant aller jusqu'à une dizaine de mètres. Le fluide est mis en mouvement grâce aux forces capillaires s'exerçant au sein d'un milieu poreux. Il existe deux catégories de boucles fluides diphasiques à pompage capillaire : les CPL (Capillary Pumped Loop) et les LHP (Loop Heat Pipe).

____________________________________________________________________

Le comparatif des technologies de transport de chaleur ci-dessous permet de voir les performances des caloducs par rapport aux boucles fluides.

Classical fluid loops : boucles fluides classiques

Macro heat pipes : macro-caloducs

Known mini fluid loops : boucles fluides miniatures

Mini-heat pipes : caloducs miniatures

 

Comparatif des performances des boucles fluides par rapport aux caloducs - source EADS Astrium

 

2. Refroidissement par air

L'avantage du refroidissement de composants par un débit d'air incident par rapport aux caloducs est sa facilité de mise en place, sa fiabilité ainsi que son faible coût. En effet, le refroidissement par liquide peut endommager les composants si une fuite venait à se déclarer dans le système. De plus, une installation en plus pour contrôler et faire circuler le fluide de refroidissement est nécessaire, ce qui pose des problèmes d'encombrement. Le refroidissement par air est généralement mis en place avec un ventilateur.

 

3. Refroidissement par convection naturelle

L'étude du refroidissement par air, c'est-à-dire par convection forcée, étant déjà un cas référencé, Epsilon nous a plutôt proposé de traiter le refroidissement des composants par convection naturelle, beaucoup plus difficile à mettre en oeuvre. La convection naturelle est bien moins efficace que la convection forcée mais un avantage réside dans l'économie de moyens de refroidissement.

3. Objectifs

L'objet de ce Bureau d'Etudes Industriel (BEI), proposé par la société Epsilon-Alcen, est de procéder à l'étude du refroidissement de composants électroniques par convection naturelle en vue d'abaisser la température de ceux-ci pour éviter leur endommagement. Ce BEI est à la fois une mise en application des connaissances théoriques de transferts thermiques ainsi que des méthodes numériques de Dynamique des Fluides Numérique (CFD - Computational Fluid Dynamics); en effet, il y a un couplage entre Mécanique des Fluides et Thermique. Cette étude sera effectuée avec deux logiciels de CFD : tout d'abord le logiciel libre OpenFOAM puis une comparaison de ces résultats avec le logiciel commercial StarCCM+.

 

En quoi ce BEI peut répondre à ces objectifs?

Afin de satisfaire la demande d'Epsilon, le présent BEI s'appuira le programme suivant :

  • validation du logiciel OpenFOAM sur un cas simple, par exemple une plaque plane chauffée soumise à un débit d'air constant : couplage entre dynamique des fluides et transferts thermiques avec codage de l'équation de la chaleur dans le solveur incompresible icoFoam
  • simulations sur le premier cas-test proposé par Epsilon : la convection naturelle autour d'un câble électrique modélisé par un cylindre
  • simulation sur le deuxième cas-test proposé par Epsilon : la convection naturelle sur un composant électronique monté sur un PCB (Printed Circuit Board)
  • comparaison des deux cas précédents avec StarCCM+
  • synthèse des méthodes numériques utilisées
  • étude de la performance et de la précision du code OpenFOAM

 

Tout au long de ce BEI, les simulations seront effectuées avec OpenFOAM 2.1.0 et StarCCM+ v7.

4. L'entreprise partenaire : Epsilon-Alcen

 

Epsilon

Ce BEI est effectué en partenariat avec la société Epsilon, filiale du Groupe Alcen, fondée en 1992, spécialisée en ingénierie thermique. Il est à noter qu'Epsilon est la seule société privée européenne de recherche dont l'expertise est mise en valeur dans les hautes technologies civiles et militaires. Elle dispose de compétences dans les domaines de la simulation et la modélisation en électrothermique, thermofluidique et thermomécanique.

 

Source image

La société Epsilon collabore avec les équipementiers et les fournisseurs de composants des domaines de l'Aeronautique, du Spatial, des Systèmes Embarqués et de l'Energie. L’objectif principal d’Epsilon est d’optimiser le comportement thermique d'équipements et de systèmes complexes en tenant compte des phénomènes physiques associés à la thermique.

                                                        Source image

La société est référencée RANG 1 auprès des principaux représentants de ces secteurs d’activité. Ses docteurs et ingénieurs spécialisés en thermique et phénomènes physiques associés apportent une réponse innovante et fiable aux intégrateurs de systèmes, aux équipementiers et aux fournisseurs de composants.

Epsilon-Alcen

 

 

 

 

5. Les logiciels de CFD utilisés : OpenFOAM et StarCCM+

En vue d'accomplir les objectifs du BEI, les études thermiques seront effectuées sur deux logiciels, à la demande d'Epsilon.

 

Le logiciel libre OpenFOAM

Présentation

OpenFOAM (Open Field Operation and Manipulation) est un logiciel libre, open source, multiplateforme et multiphysique spécialement destiné à la Dynamique des Fluides Numérique (CFD) et développé par la société britannique OpenCFD Ltd mais initialement conçu à l'Imperial College de Londres dans les années 1980 par Henry Weller. OpenFOAM est codé en C++, utilise une approche orientée objet et est principalement dédié à la résolution des équations aux dérivées partielles par la méthode Volumes Finis. Ses principaux avantages sont qu'il est gratuit, performant, parallélisable et personnalisable, c'est-à-dire que des modules complémentaires peuvent être développés car sa source est accessible[2].

 

Solveurs

OpenFOAM, entre autres capable de résoudre des problèmes multiphysiques, comporte de nombreux solveurs :

  • écoulements incompressibles
  • écoulements compressibles
  • écoulements multiphasiques
  • transferts thermiques
  • écoulements particulaires
  • combustion
  • finance...

 

Modèles de turbulence

Outre sa performance en écoulement laminaire, OpenFOAM permet aussi de traiter des écoulements turbulents avec plusieurs modèles de turbulence :

  • Large Eddy Simulation (LES) : résout les grandes échelles de l'écoulement
  • Reynolds-averaged Navier-Stokes (RANS) : résout les équations de Navier-Stokes moyennées
  • Detached Eddy Simulation (DES) : modèle hybride entre la LES et le RANS
  • Direct Numerical Simulation (DNS) : résout d'équation de Navier-Stokes directement sans modèle de turbulence

Exemple de simulation OpenFOAM d'un champ de contraintes sur une plaque percée

 

Première prise en main d'OpenFOAM

La structure générale du logiciel OpenFOAM est présentée ci-dessous :

Chaque cas de calcul dans OpenFOAM est organisé selon une certaine structure. L'utilisateur donne un nom à son cas de calcul, par exemple "Validation". Ce nom devient alors le nom du dossier dans lequel tous les fichiers relatifs à ce cas de calcul seront stockés. Peu importe l'endroit où ce dossier est placé, mais il est préférable qu'il le soit dans le dossier run du dossier du projet de l'utilisateur $HOME/OpenFOAM/${USER}-2.1.0.

Capture du dossier {USER}-2.1.0, ici, login-2.1.0

 

Contenu du dossier "run" et du cas de calcul "Validation" avant simulation

 

Structure du dossier contenant le cas de calcul

 

Le logiciel commercial StarCCM+

Présentation

Le logiciel StarCCM+ (Simulation of Turbulent flow in Arbitrary Regions Computational Continuum Mechanics) est un outil commercial de CFD développé par la société CD-adapco. Il résout les équations de la Mécanique des Fluides grâce à la méthode Volumes Finis, tout comme le logiciel OpenFOAM présenté précédemment. StarCCM+ permet de simuler les écoulements diphasiques, incompressibles, compressibles ou encore la combustion dans différents domaines tels l'Aerospatial, l'Automobile ou l'Environnement. StarCCM+ comprend un outil de modélisation 3D de CAD (Computer-aided Design), une technologie de maillage automatique, des modélisations de la turbulence (RANS, Reynolds Stress Model, k-$\epsilon$, k-$\omega$, v2f, LES, Spalart Allmaras...) ainsi que des outils de post-traitement.

 

Exemple de simulation CFD sur StarCCM+ autour d'un navire de combat

 

 

 

I. Le phénomène de convection

La convection est un des trois modes de transfert de chaleur avec la conduction et le rayonnement. Le terme de convection fait référence aux transferts de chaleur se produisant entre une surface et un fluide en mouvement lorsque ceux-ci sont à des températures différentes. En plus du transfert d'énergie dû à la diffusion, il y a également transfert par le biais du mouvement du fluide. Ce dernier est associé au fait que de multiples molécules ont un mouvement collectif, ce qui implique un transfert de chaleur dans le cas où il existe un gradient thermique.

La contribution due au mouvement aléatoire des molécules, la diffusion, domine près de la surface où la vitesse du fluide est faible. En effet, à l'interface entre la surface et le fluide, étant donné que la vitesse du fluide est nulle, le seul mode de transfert est la diffusion. La contribution due au mouvement du fluide tient son origine du fait que la couche limite croît au fur et à mesure de l'avancée du fluide sur la surface.

Le transfert thermique par convection est divisé en deux parties suivant la nature de l'écoulement :

 

Les deux modes de transfert de chaleur par convection

 

Pour mesurer l'intensité du transfert thermique dans le fluide dû à ses mouvements et pour caractériser l'échange thermique entre le fluide et la paroi, on utilise le nombre de Nusselt.

Le nombre de Nusselt représente le gradient thermique adimensionné à la paroi :

$$Nu=\left (\frac {\partial T^+} {\partial y^+}\right )_{paroi}= \frac {h  x} {\lambda_{fluide}}$$

 

nombre de Reynolds :  $Re=\frac {U  L}{\nu}$

nombre de Prandtl :  $Pr=\frac {\nu}{\alpha}$

nombre de Grashof :  $Gr=\frac{g  \beta  \Delta T  L^3}{\nu^2}$

avec :

$U$ : vitesse caractéristique $(m/s)$

$L$ : longueur caractéristique $(m)$

$\nu$ : viscosité cinématique $(m^2/s)$

$\alpha$ : diffusivité thermique $(m^2/s)$

$g$ : accélération de la pesanteur $(m/s^2)$

$\beta$ : coefficient de dilatation thermique $(K^{-1})$

II. Données physiques

Dans nos problèmes de validation et des deux cas-test, le fluide considéré est de l'air pour lequel ses propriétés physiques varient avec la température, soit 293 K (20°C), on utilise donc les relations suivantes pour les calculer :

$\nu = -1,363.10^{-14} T^3 + 1,008.10^{-10} T^2+3,452.10^{-8} T-3,401.10^{-6}=1,495.10^{-5} m^2/s$

$\rho=1,293 \frac {273} {273+t(°C)}=1,205  kg.m^{-3}$

$\eta=\rho \times \nu = 1,8.10^{-5} Pl$

$\lambda=1,521.10^{-11} T^3-4,857.10^{-8}T^2+1,018.10^{-4}T-3,933.10^{-4}=0,026  W.m^{-1}.K^{-1}$

$C_p=1,933.10^{-10}T^4-7,999.10^{-7}T^3+1,141.10^{-3}T^2-4,489.10^{-1}T+1,057.10^{3}=1004,72  J.kg^{-1}.K^{-1}$

$\alpha=\frac {\lambda} {\rho C_p}=2,14.10^{-5} m^2/s$

III. Validation du logiciel OpenFOAM

Avant d'aborder les simulations numériques avec OpenFOAM sur les cas-tests proposés par Epsilon, il faut au préalable avoir au moins un point de validation du logiciel. Notre validation portera sur l'étude de la convection forcée entre une plaque plane chauffée et un écoulement d'air à température inférieure.

1. Analyse préliminaire

Cette validation concerne l'étude du transfert thermique entre une plaque plane chauffée et un courant d'air incident. Il s'agit par conséquent d'un cas de transfert thermique par convection forcée.

 

Présentation générale d'un écoulement sur plaque plane chauffée

Un écoulement incident sur une plaque plane chauffée sera principalement caractérisé par le développement d'une couche limite dynamique (en vitesse) et d'une couche limite thermique (en température).

  • L'origine du développement de la couche limite dynamique provient des effets de viscosité, c'est-à-dire de diffusion de quantité de mouvement, entre le fluide incident et la plaque. En effet, la vitesse du fluide est nulle sur la plaque, tandis qu'elle est non nulle à l'infini, ce qui implique qu'au voisinage de la plaque il y ait une zone de transition : la couche la plus lente freine la couche la plus rapide, qui, en retour, l'accélère. La couche limite dynamique est une zone où la variation de vitesse est la plus marquée, c'est-à-dire là où l'on trouve de forts gradients de vitesse. La vitesse à l'extrémité de la couche limite est égale à 99% de la vitesse à l'infini, on a donc :

$$U(\delta)=0,99  U_0$$

L'épaisseur de la couche limite varie selon la relation suivante[3] :

$$\delta(x)= 4,64 \sqrt \frac {\nu  x} {U_0}$$    

$\nu$ : viscosité cinématique du fluide incident en $m^2/s$

$x$ : abscisse de la plaque en $m$

$U_0$ : vitesse du fluide incident en $m/s$

 

  • De la même manière que la couche limite dynamique, la couche limite thermique est une zone proche de la plaque où la température varie de façon significative entre la température de la plaque et la température du fluide incident. La température à l'extrémité de la couche limite est telle que :

$$\frac{T(\delta_T)-T_P} {T_\infty - T_P} = 0,99$$

 

L'épaisseur de la couche limite thermique varie selon l'expression suivante[3] :

$$\delta_T= \frac { \delta} {1,026  Pr^{1/3}} =5,09 \sqrt \frac {\nu  x} {U_0}$$  

$\alpha$ : diffusivité thermique du fluide incident en $m^2/s$

Régimes d'écoulement sur plaque plane

 

Position du cas de validation

Notre cas de validation se présente tel que l'illustre la figure ci-dessous. De l'air incident à une vitesse de 50 cm/s et à une température de 293 K arrive sur une plaque plane chauffée à une température de 348 K.

Schéma du cas de validation du logiciel OpenFOAM

 

Avant de procéder à toute simulation, une analyse préliminaire est requise afin d'identifier les phénomènes et ordres de grandeur.

 

Hypothèses et mécanismes prédominants

  • écoulement stationnaire
  • écoulement établi
  • écoulement monophasique, pas de phénomène d'évaporation ou de condensation
  • problème à deux dimensions x et y
  • invariance par translation suivant l'axe z
  • convection prépondérante devant conduction et rayonnement
  • convection forcée
  • température considérée comme un scalaire passif (aucune influence sur l'hydrodynamique)

 

Échelles caractéristiques

  • longueur de la plaque : $L = 10  cm$
  • vitesse incidente : $U_0=50  cm/s$
  • iso-contour de température à l'épaisseur de la couche limite thermique : $T(\delta_T)=0,99  (T_{\infty}-T_P)+T_P=293,55  K$
  • iso-contour de vitesse à l'épaisseur de la couche limite dynamique :  $U(\delta)=0,495   m/s$

 

Détermination des nombres adimensionnels

  • nombre de Reynolds

Pour une plaque plane, le nombre de Reynolds de transition entre les régimes laminaire et turbulent est environ :

$${Re_{tr} \approx 5.10^5} $$

Dans notre cas de validation, le nombre de Reynolds à l'extrémité de la plaque s'écrit :

$${Re}=\frac{U_0 L} {\nu_{air}}=\frac {50.10^{-2} \times 10.10^{-2}} {1,495.10^{-5}}=3344,5$$

L'écoulement au-dessus de la plaque est donc laminaire en tout point de celle-ci.

  • nombre de Prandtl

$$Pr=\frac {\nu_{air}} {\alpha_{air}}=0,698$$

 

Corrélation théorique

Dans notre cas de convection forcée sur plaque plane, d'après l'ouvrage Fundamentals of Heat and Mass Transfers, Sixth Edition de Frank P. Incropera et David P. DeWitt, la corrélation qui sera utilisée pour vérifier nos résultats est la suivante, qui relie le nombre de Nusselt local aux nombres de Prandtl et de Reynolds locaux :

 

$$Nu_x = \frac{h_x x}{\lambda} = f(Pr) {Re_x}^{0,5}$$

avec :                                            $f(Pr)=0,332  Pr^{1/3}$    pour    $0,6<Pr<50$

 

2. Création du maillage

Le maillage sous OpenFOAM a été réalisé avec l'outil blockMesh. Cet outil permet de mailler des géométries relativement simples puisqu'il s'agit d'un mailleur minimaliste. Dans le cas d'une géométrie beaucoup plus compliquée qu'une plaque plane, on lui préfèrera des mailleurs externes tels que Salomé, I-deas, Gmsh... Cependant, OpenFOAM dispose d'autres outils de maillage.

 

Caractéristiques du maillage

Avec blockMesh, le maillage doit être entièrement réalisé sans interface graphique. Ainsi, il faut définir point par point les différents éléments constitutifs de la géométrie dans un fichier texte, puis spécifier les paramètres de maillage.

Pour notre cas de validation, nous devons définir 2 zones : celle dont la surface inférieure sera la plaque, et celle en amont de la plaque. Nous raffinons ensuite au niveau de plaque pour une meilleure résolution des couches limites dynamique et thermique en utilisant un maillage géométrique avec un rapport 100 entre la maille la plus grande (le plus loin de la paroi) et la maille la plus petite (proche de la paroi).

Schéma explicatif de la géométrie, division du domaine en deux blocs

 

Maillage du domaine obtenu avec blockMesh après visualisation avec paraFoam

 

Les caractéristiques du maillage sont éditées dans le fichier blockMeshDict.

 

Capture du fichier blockMeshDict correspondant à notre maillage de la plaque plane

3. Simulations et résultats

Mise en place des simulations

Dans ce cas de validation, nous allons faire l'hypothèse que la température n'a aucune influence sur l'hydrodynamique, c'est à dire que la température est considérée comme scalaire passif. L'écoulement est incompressible et les vitesses considérées en entrée sont telles que la transition vers un état turbulent ne sera jamais observée le long de la plaque.

Dans ces conditions, nous avons choisi d'utiliser le solveur incompressible icoFoam et d'y intégrer l'équation de la chaleur (avec terme de convection), le nouveau solveur porte ainsi le nom de my_icoFoam.

$$ \frac {\partial T} {\partial t} + Div(\phi T) = \frac {\lambda} {\rho c_p} \Delta T $$

Cette modification du solveur a déjà été étudiée : un tutoriel en anglais est disponible à cette adresse. L'idée pour nous était de voir comment il est possible d'intégrer simplement de nouvelles variables et équations à un modèle déjà existant. Ceci nous permet d'évaluer le potentiel d'OpenFOAM en terme de flexibilité, un des atouts majeurs de tout logiciel "open source".

De plus, avant toute simulation, il faut veiller à avoir un nombre de Courant inférieur à 1 :

$$Co=\frac {U  \Delta t} {\Delta x} < 1$$

$U$ : vitesse dans la direction $x$

$\Delta t$ : intervalle de temps à choisir pour avoir un nombre de Courant au plus égal à 1

$\Delta x$ : distance entre deux mailles

Dans notre cas, le pas de temps à choisir est au plus :

$$\Delta t=\frac {Co  \Delta x} {U} =\frac{1 \times  1.10^{-3}} {50.10^{-2}}=2  ms$$

 

Visualisation des résultats

La procédure pour effectuer la simulation et visualiser les résultats se fait par le biais du terminal en trois parties principales :

0. Se placer dans le répertoire de travail contenant le cas de calcul

         

Contenu du cas de calcul "Validation" avant simulation

 

1. Mailler le domaine en appelant l'outil blockMesh

 

2. Résoudre les équations mises en jeu en appelant le solveur approprié, ici my_icoFoam

3. Afficher les résultats dans le logiciel Paraview en appelant l'outil de post-traitement paraFoam

Champ de température sur la plaque avec l'iso-contour $T(\delta_T)=293,55  K$ avec une vitesse d'entrée $U_0=50  cm/s$ en régime établi (t > 0,4 s)

 

  

Champ de vitesse sur la plaque avec l'iso-contour $U(\delta)=0,495  m/s$ avec une vitesse d'entrée $U_0=50  cm/s$ en régime établi (t > 0,4 s)

 

 

Visualisation des couches limites

Afin de comparer les épaisseurs des couches limites thermique et dynamique, on peut les visualiser sous Paraview en utilisant la fonction isosurface en spécifiant les valeurs correspondantes en vitesse et température.

En exportant les résultats sur Matlab, il est alors possible de superposer les profils numérique et théorique :

 

Tracé des couches limites thermique et dynamique sur la plaque plane d'après les iso-contours OpenFOAM

 

On remarque d'une part que la couche limite thermique se développe plus vite que la couche limite dynamique. Ceci s'explique par le fait que le nombre de Prandtl est inférieur à 1.

D'autre part, on constate un certain écart avec la théorie : la couche limite thermique est sur-estimée tandis que la couche limite dynamique est quant à elle sous-estimée. Deux raisons peuvent expliquer cet écart :

  • Vitesse et température loin de la plaque sont mal évaluées (par exemple, la vitesse loin de la plaque n'est pas exactement égale à la vitesse imposée en entrée car la plaque modifie l'écoulement)
  • Les paramètres $\nu$ et $\alpha$ ne sont pas estimés à la bonne température

 

 

Evolution du Nusselt local

Afin de calculer le Nusselt local le long de la plaque, il nous faut connaître l'évolution du gradient de température sur la plaque, perpendiculairement à celle-ci.

Ceci a été réalisé en ajoutant un nouveau champ dans le fichier icoFoam_gradT.C et calculé seulement le long de la plaque. Le code à rajouter doit s'effectuer sur le fichier précédent placé dans le dossier nommé icoFoam_gradT localisé dans les dossiers applications puis solvers dont le contenu est le suivant :

Les lignes de code à ajouter sont les suivantes :

Capture du fichier texte "icoFoam_gradT.C"

 

La nouvelle variable ainsi créée apparaît désormais dans Paraview, nous permettant de connaître l'évolution de la dérivée normale à la paroi de la température. En adimensionnalisant cette grandeur (on la multiplie par $ \frac{x}{\Delta T} $), on obtient le profil du Nusselt local.

La figure suivante montre les résultats obtenus pour différentes vitesses en entrée :

vitesses nombre de Reynolds
$U_0=5  mm/s$ $33,44$
$U_0=5  cm/s$ $334,4$
$U_0=50  cm/s$ $3344$
$U_0=5  m/s$ $33440$

 

Le tracé du Nusselt montre d'excellents résultats vis-à-vis de la théorie, même dans le cas d'une vitesse d'entrée de 5 m/s qui donne un nombre de Reynolds s'approchant de la valeur critique de transition laminaire/turbulente. On peut cependant observer quelques différences au niveau du bord d'attaque de la paroi : cela s'explique par la résolution du maillage qui n'est pas assez élevée dans cette zone où les gradients sont les plus forts.

 

IV. Premier cas-test : convection naturelle autour d'un câble électrique

Ce premier cas-test traite de l'étude de la convection naturelle autour d'un câble électrique modélisé par un cylindre. Ce câble est soumis à un fort courant électrique et dissipe la chaleur qui est évacuée par convection naturelle entre le câble et l'air ambiant.

 

Cahier des charges

Les données utiles au traitement du problème, fournies par Epsilon, sont les suivantes :

 

 

1. Analyse préliminaire

 

Présentation générale d'un cas de convection naturelle autour d'un cylindre

En convection naturelle, lorsqu'un cylindre horizontal, de température $T_P=75 °C$ est immergé dans un fluide de température $T_{\infty}=20 °C$, une couche limite thermique se développe autour de celui-ci. D'après la figure ci-dessous[1], on peut constater que l'épaisseur de la couche limite thermique autour du cylindre est une fonction de l'angle $\theta$.

Panache de convection et couche limite thermique autour d'un cylindre chauffé

 

Analyse du cas-test

 

Hypothèses et mécanismes prédominants

  • écoulement stationnaire
  • écoulement établi
  • écoulement monophasique, pas de phénomène d'évaporation ou de condensation lors du transfert
  • problème à deux dimensions x et y
  • invariance par translation suivant l'axe z
  • symétrie selon le plan (y0z)
  • convection prépondérante devant conduction et rayonnement
  • convection naturelle
  • hypothèse de Boussinesq  $\rho=\rho_0 [1-\beta(T-T_0)]$   avec   $\beta= - \frac{1} {\rho} \left (\frac{\partial \rho} {\partial T} \right)_{P=P_0} $ le coefficient de dilatation thermique

 

Échelles caractéristiques

  • le diamètre du cylindre est $L = 1  cm$

 

Détermination des nombres adimensionnels

  • nombre de Prandtl

$$Pr=\frac {\nu_{air}} {\alpha_{air}}=0,698$$

 

  • nombre de Rayleigh basé sur le diamètre D du cylindre

$$Ra_D=\frac{g \beta (T_P-T_{\infty})D^3} { \nu \alpha}=5,26.10^3$$

avec, en considérant l'air comme un gaz parfait :

$\beta=\frac {1} { \overline T}=\frac {2} {T_P + T_{\infty}} = 3,12.10^{-3} K^{-1}$

 

  • nombre de Nusselt : déterminé par la corrélation théorique

 

Corrélation théorique

Dans le cas de la convection naturelle autour d'un cylindre, la corrélation à vérifier lors de nos simulations numériques est la suivante (Churchill et Chu - 1975) :

                              $\overline {Nu_D}=  \left \{0,60 + \Large{\frac {0,387 Ra_D^{1/6}} { \left [1+ \left (\frac{0,559} {Pr} \right)^{9/16} \right ]^{8/27}}} \right \}^2$       pour       $Ra_D<10^{12}$

Dans nos analyses, nous tracerons donc le nombre de Nusselt en fonction du nombre de Rayleigh afin de comparer les résultats numériques d'OpenFOAM et StarCCM+ avec cette corrélation.

 

2. Simulations OpenFOAM

Dans cette partie, nous allons simuler l'écoulement de convection naturelle autour du cylindre à l'aide du logiciel OpenFOAM. Comme nous allons le voir plus en détails par la suite, le phénomène de convection naturelle étudié ici va imposer certaines contraintes au niveau de la résolution numérique. Un des objectifs de l'étude est de comparer les résultats obtenus pour différents modèles de turbulence : standard k-$\epsilon$, realizable  k-$\epsilon$, k-$\omega$, k-$\omega$  SST (Shear Stress Transport) et Spalart Allmaras.

 

Considérations numériques

  • Une des difficultés rencontrées lorsque l'on souhaite simuler une configuration de type convection naturelle réside dans la définition de conditions aux limites entrée/sortie. En fait, le fluide, initialement au repos, est mis en mouvement par l'apparition de gradients de densité, eux-mêmes dus aux variations de température. L'air chauffé par le cylindre monte, puis se refroidit , faisant apparaître un courant descendant car plus dense que l'air chaud. Imaginons que l'on fixe une condition d'entrée au niveau de la frontière inférieure et de sortie au niveau de la frontière supérieure, en imposant seulement une valeur de pression. Du fait de l'apparition de rouleaux de convection dans le domaine, les 2 frontières considérées voient à la fois du fluide sortir et rentrer simultanément. D'un point de vue numérique, ceci n'est pas acceptable pour une grande partie des codes de calcul CFD, y compris OpenFOAM ou du moins nous n'avons pas réussi à trouver les conditions aux limites adéquates. Nous sommes donc obligés d'isoler complètement le domaine, en imposant une condition de type wall sur toutes les frontières.

 

  • Le fait d'isoler complètement le cylindre va avoir une conséquence sur le choix du solveur sous OpenFOAM. En effet, nous avons le choix pour ce problème thermique entre 4 solveurs :
    • buoyantBoussinesqPimpleFoam
    • buoyantBoussinesqSimpleFoam
    • buoyantPimpleFoam
    • buoyantSimpleFoam

Les deux premiers solveurs (Boussinesq) fonctionnent en écoulement incompressible sous l'hypothèse de Boussinesq, tandis que les deux derniers fonctionnent en écoulement compressible. Étant données les hypothèses formulées précédemment, nous devons choisir parmi les deux solveurs Boussinesq.

Pimple se réfère à l'étude du régime transitoire, Simple le régime permanent. Comme nous considérons une géométrie complètement fermée, le cylindre va chauffer continuellement la pièce jusqu'à ce qu'à atteindre un état permanent où la température du domaine sera uniforme et égale à celle du cylindre. Ce n'est pas ce que nous souhaitons observer a priori, nous écartons donc également le solveur en régime permanent.

Nous allons donc effectuer nos simulations avec le solveur buoyantBoussinesqPimpleFoam, le but étant d'étendre la géométrie au maximum (attention au temps de calcul cependant) et de constater l'établissement de la couche limite autour du cylindre ainsi que le panache qui s'en dégage, avant les que re-circulations ne perturbent l'ensemble.

 

Maillage de la géométrie

La géométrie étant assez simple à réaliser, nous avons utilisé l'outil blockMesh pour ce problème également. Afin de simuler correctement la couche limite thermique proche du cylindre, nous raffinons le maillage dans cette zone. La gravité est selon l'axe horizontal et orientée vers la gauche (le fluide chaud "monte" vers la droite). Nous avons ici, ainsi que pour le second cas-test pour OpenFOAM, seulement créé la moitié de la géométrie afin de simplifier l'écriture des coordonnées dans le fichier blockMeshDict.

 

Maillage du domaine avec raffinement autour de la paroi du cylindre

 

Exploitation des résultats

 

Visualisation des résultats

Comme énoncé précédemment, le solveur que nous utilisons simule le régime transitoire. Nous avons donc dû estimer, pour une simulation donnée, le temps à partir duquel le régime permanent est atteint. Cela est réalisé à l'aide de Paraview : en traçant l'évolution du nombre de Nusselt le long du cylindre et en faisant varier le temps de simulation, on peut observer que la courbe se stabilise pendant un certain temps, avant de subir des changements sous l'effet des re-circulations provoquées par le fermeture complète du domaine. Ce laps de temps pendant lequel la courbe se stabilise nous donne alors une assez bonne approximation du régime permanent au voisinage du cylindre (seule l'évolution du nombre de Nusselt le long du cylindre nous intéresse).

Nous avons utilisé plusieurs valeurs de gravité ($g=10  m/s^2$, $g=100  m/s^2$, $g=1000  m/s^2$) afin de faire varier le nombre de Rayleigh. Parallèlement, plusieurs modèles de turbulence sont utilisés.

 

Résultats obtenus avec le modèle de turbulence $k-\omega$ :

 

Champs de température

 

 

 

 

 

 

 

 

Résultat de simulation pour le modèle k-$\omega$ et pour g=10 m/s2

 

Résultat de simulation pour le modèle k-$\omega$ et pour g=100 m/s2

 

Résultat de simulation pour le modèle k-$\omega$ et pour g=1000 m/s2

 

Champs de vitesse

Résultat de simulation pour g=10 m/s2

 

Résultat de simulation pour g=100 m/s2

 

Résultat de simulation pour g=1000 m/s2

 

L'augmentation de la gravité, et donc du nombre de Rayleigh, s'accompagne :

  • d'une réduction de l'épaisseur de la couche limite thermique autour du cylindre : les gradients thermiques sont plus élevés
  • d'une augmentation de la vitesse du fluide particulièrement visible au coeur du panache (les échelles de vitesses ne sont pas les mêmes sur chaque image). Pour la dernière valeur de gravité, le fluide atteint une vitesse de plus d'un mètre par seconde, alors que le cylindre et le mur de droite (correspondant à la frontière supérieure en réalité) ne sont séparés que de 15 cm. Il faut donc être très précis quant au moment à choisir pour l'estimation du régime permanent.

Des simulations pour des plus grandes valeurs de gravité ($g=10  000  m/s^2$ et $g=100  000  m/s^2$) ont été réalisées afin d'étendre encore plus l'intervalle de Rayleigh. Seulement, et cela rejoint la remarque faite précédemment, le domaine n'est plus assez étendu, la solution étant vite contaminée par les re-circulations (cf l'évolution du Nusselt sur la figure qui suit). De ce fait, il n'a pas été possible d'étendre l'intervalle de nombre de Rayleigh, les calculs étant trop volumineux.

 

Calcul du nombre de Nusselt

Afin de pouvoir visualiser l'évolution du Nusselt local le long du cylindre, nous avons choisi d'utiliser l'outil de post-traitement wallHeatFlux. Cet outil permet, une fois la simulation effectuée, de calculer les flux thermiques et les Nusselt associés sur toutes les frontières de type wall. Cependant, il est pré-programmé pour ne fonctionner qu'avec les solveurs simulant le phénomène de combustion en écoulement compressible (par exemple : reactingFoam). Il a fallu donc l'adapter pour le rendre compatible au solveur buoyantBoussinesqPimpleFoam qui, rappelons le, fonctionne en incompressible.

Ce problème se résout notamment en changeant le modèle thermophysique de l'outil : celui-ci se doit d'être le même que celui utilisé par le solveur buoyantBoussinesqPimpleFoam. Une fois correctement re-compilé et exécuté, l'outil rajoute dans chaque répertoire de temps un nouveau champ, NusseltNumber, calculé le long de toutes les parois.

 

Évolution du nombre de Nusselt local

La figure qui suit montre l'évolution du nombre de Nusselt local autour du cylindre, pour différentes valeurs du nombre de Rayleigh associées aux valeurs de gravité décrites précédemment. Le modèle de turbulence utilisé est le modèle k-$\epsilon$.

Évolution du nombre de Nusselt en fonction de l'angle autour du cylindre en régime turbulent avec le modèle k-$\epsilon$

 

Les courbes correpondantes aux valeurs $Ra_D=4.10^3$, $4.10^4$ et $4.10^5$ montrent d'excellents résultats qualitativement parlant, au vu de la courbe théorique présentée dans la partie préliminaire. Comme prévu, les 2 dernières courbes ($Ra_D=4.10^6$ et $4.10^7$) s'écartent de la courbe théorique, mais peuvent nous donner un ordre de grandeur du Nusselt local sur la partie inférieure du cylindre ($ \theta = - \pi /2$), qui est moins perturbée par les re-circulations.

 

Calcul du nombre de Nusselt moyen

Nous nous intéressons maintenant au nombre de nusselt moyen, $Nu_D$, qui s'obtient en intégrant le nombre de Nusselt local le long du cylindre, tracé en fonction du nombre de Rayleigh pour les différents modèles de turbulence considérés. Ces résultats sont comparés à la corrélation de Churchill et Chu :

 

Évolution du nombre de Nusselt moyen en fonction du nombre de Rayleigh pour divers modèles de turbulence

On peut remarquer que pour les "petites" valeurs de $Ra_D$, c'est-à-dire pour $Ra_D<10^5$, les valeurs de $Nu_D$ sont très proches de la courbe théorique. Cela est particulièrement visible pour les modèles k-$\omega$ et k-$\omega$ SST. Néanmoins, il semble que globalement, pour ces valeurs de $Ra_D$, le nombre de Nusselt moyen soit légèrement sur-estimé, ce qui se traduit par une sur-estimation du flux dégagé par le cylindre.

Lorsque le nombre de Rayleigh augmente encore, la tendance inverse est observée (sous-estimation du Nusselt moyen). Les modèles standard  k-$\epsilon$ et Spalart Allmaras, qui donnent l'impression de globalement sur-estimer le nombre de Nusselt moyen pour les faibles valeurs du nombre de Rayleigh comparativement aux autres modèles, sont ceux qui donnent les meilleurs résultats à haut nombre de Rayleigh.

 

3. Simulations StarCCM+

Sous le logiciel commercial StarCCM+, les simulations sont beaucoup plus évidentes à mettre en place que sous OpenFOAM. Il faut par contre avoir les bonnes conditions aux limites et les bons modèles.

 

Géométrie du domaine et conditions aux limites

Afin de pouvoir observer le panache de convection, nous avons créé une géométrie rectangulaire beaucoup plus haute que large, en tenant compte des dimensions voulues par Epsilon pour le cylindre, celui-ci étant représenté par une condition limite de type wall à température fixée.

Schéma de la géométrie considérée et des conditions aux limites introduites dans StarCCM+

 

Maillage de la géométrie

Compte tenu du développement de la couche limite thermique sur le cylindre, il semble judicieux de raffiner le maillage d'autant plus que l'on se rapproche de la paroi du cylindre. Nous avons réalisé le maillage sous Ansys IcemCFD. Les détails du nombre de mailles et des lois de maillage sont présentées ci-dessous.

Maillage du domaine avec raffinement au niveau de la paroi du cylindre

 

Zoom sur la partie proche de la paroi du cylindre

 

Mise en place des simulations

Afin d'avoir des résultats les plus proches possible de ceux attendus, nous avons sélectionné plusieurs modèles dans StarCCM+ comme le montre la capture d'écran suivante. Nous avons entre autres sélectionné l'option segregated flow par opposition à coupled flow. La résolution de la vitesse avec segregated flow s'opère indépendamment de la pression, ce qui produit des simulations plus rapides, mais n'est pas adapté dans le cas d'écoulements supersoniques. Étant donné que nous avons des conditions de pressure outlet en entrée et en sortie du domaine, puisqu'il n'y a pas de vitesse en entrée, nous avons fixé la pression à 0 Pa. Nous avons imposé au cylindre une température de 348 K (75 °C).

Capture des modèles utilisés dans StarCCM+ pour la convection naturelle autour du cylindre

 

Résultats et vérification de la corrélation théorique

Nous avons ensuite réalisé les simulations pour différents modèles de turbulence qui seront présentés plus loin afin de voir l'influence qu'ils peuvent avoir sur la couche limite thermique et/ou sur la panache de convection.

 

Résultats en régime laminaire

Dans un premier temps, les simulations ont été effectuées en régime laminaire.

Résultat de simulation pour $T_P=348  K$  et  $T_\infty=293 K$ avec un modèle laminaire

 

L'objectif est de vérifier si les simulations apportent le résultat attendu. La figure ci-dessous montre l'évolution du nombre de Nusselt en fonction de l'angle $\theta$ autour du cylindre. Nous pouvons remarquer un maximum du nombre de Nusselt sur la partie inférieure du cylindre, c'est-à-dire pour un angle $\theta=-\frac {\pi} {2}$ car c'est sur cette partie que la couche limite est la plus fine, il y a donc des gradients thermiques beaucoup plus importants que sur la partie supérieure du cylindre, le nombre de Nusselt y est donc d'autant plus élevé.

Évolution du nombre de Nusselt en fonction de l'angle autour du cylindre en régime laminaire

 

Vérification de la corrélation

Maintenant que l'allure du nombre de Nusselt précédent est cohérent par rapport aux considérations physiques, nous allons dans un second temps vérifier la corrélation théorique de Churchill et Chu. Nous avons tracé les nombres de Nusselt théoriques puis ceux obtenus numériquement.

Rappelons la corrélation à vérifier ainsi que l'expression du nombre de Rayleigh :

                                            $\overline {Nu_D}=  \left \{0,60 + \Large{\frac {0,387 Ra_D^{1/6}} { \left [1+ \left (\frac{0,559} {Pr} \right)^{9/16} \right ]^{8/27}}} \right \}^2$       pour       $Ra_D<10^{12}$

 

Pour faire varier le nombre de Rayleigh, nous avons choisi de faire varier la gravité ainsi que la différence de température entre le cylindre et l'air environnant.

  • Pour la première figure ci-dessous, nous avons fait varier le nombre de Rayleigh pour les gravités suivantes : $g=9,81  m/s^2$, $g=98,1  m/s^2$, $g=981  m/s^2$ et $g=9810  m/s^2$.

Nombres de Nusselt numérique et théorique en fonction du nombre de Rayleigh calculé pour différentes gravités

 

  • Pour la seconde figure ci-dessous, nous avons fait varier le nombre de Rayleigh pour les différences de température : $\Delta T=55 °C$, $\Delta T=65 °C$, $\Delta T=75 °C$, $\Delta T=85 °C$, $\Delta T=95 °C$, $\Delta T=200 °C$  et  $\Delta T=1000 °C$.

Nombres de Nusselt numérique et théorique en fonction du nombre de Rayleigh calculé pour différents $\Delta T$

 

Les courbes numériques ci-dessous font apparaître des résultats relativement proches de la corrélation théorique avant de diverger de manière assez considérable à partir de nombres de Rayleigh relativement faibles. Dans le cas où nous avons fait varier le nombre de Rayleigh avec des gravités différentes (première figure), l'origine de cette divergence pourrait s'expliquer par l'aspect peu physique d'une gravité très importante, puisque nous somme allés jusqu'à une accélération de la pesanteur de $9810  m/s^2$.

 

Comparaison des modèles de turbulence

Nous avons également procédé à des simulations pour divers modèles de turbulence afin de voir leur influence sur la couche limite thermique ainsi que sur le panache de convection.

 

Résultats de simulation pour $T_P=348  K$  et  $T_\infty=293 K$ avec différents modèles de turbulence

Les simulations ci-dessous font apparaître des différences bien plus marquées au niveau du panache de convection qu'au niveau de la couche limite thermique. En effet, c'est dans ce panache que les effets de la turbulence peuvent se faire davantage ressentir du fait des instabilités qui peuvent se produire ainsi que les recirculations.

 

Le tracé du nombre de Nusselt ci-dessous permet aussi d'apprécier les différences qui apparaissent selon le modèle de turbulence utilisé. Ceux-ci sont quasi-confondus sur la partie supérieure du cylindre en $\theta=\frac {\pi} {2}$ tandis qu'ils diffèrent d'une unité pour les modèles Reynolds Stress (RSM) et realizable k-$\epsilon$ sur la partie inférieure de celui-ci. On peut donc en conclure que, quel que soit le modèle de turbulence utilisé, le panache qui s'élève du cylindre n'a aucune influence sur l'évolution du nombre de Nusselt.

Évolution du nombre de Nusselt numérique en fonction de l'angle autour du cylindre pour différents modèles de turbulence

 

V. Deuxième cas-test : convection naturelle autour d'un composant électronique

Ce second cas-test proposé par Epsilon traite de l'étude de la convection naturelle autour d'un composant électronique monté sur un circuit imprimé. Le composant électronique dissipe la chaleur et il faut donc contrôler sa température pour éviter un dysfonctionnement. Le problème est modélisé par un bloc dissipant un flux de chaleur monté sur une plaque plane horizontale. L'objectif est d'évaluer le champ de température dans le domaine fluide entourant le composant.

 

Cahier des charges

Les données numériques qui nous ont été fournies par Epsilon sont les suivantes :

 

1. Analyse préliminaire

Analyse du cas-test

 

Hypothèses et mécanismes prédominants

  • écoulement stationnaire
  • écoulement établi
  • écoulement monophasique, pas de phénomène d'évaporation ou de condensation lors du transfert
  • problème à deux dimensions x et y
  • invariance par translation suivant l'axe z
  • symétrie selon le plan (y0z)
  • convection prépondérante devant conduction et rayonnement
  • convection naturelle
  • hypothèse de Boussinesq  $\rho=\rho_0 [1-\beta(T-T_0)]$   avec   $\beta= - \frac{1} {\rho} \left (\frac{\partial \rho} {\partial T} \right)_{P=P_0} $

 

Échelles caractéristiques

  • la longueur du composant est $L=10  cm$
  • la largeur du composant est $l=10  cm$

 

Détermination des nombres adimensionnels

  • nombre de Prandtl

$$Pr=\frac {\nu_{air}} {\alpha_{air}}=0,698$$

2. Simulations OpenFOAM

Note : l'étude de ce cas-test sous OpenFOAM est incomplète, par manque de temps. Nous avons néanmoins réussi à simuler l'écoulement résultant du phénomène de convection naturelle

 

Maillage de la géométrie

 

Pour ce problème, nous allons considérer une géométrie 2D, avec symétrie axiale. L'outil blockMesh est encore une fois suffisant pour ce genre de géométrie. Le domaine rectangulaire est découpé en 3 blocs. Le maillage est raffiné au voisinage de la puce.

Comme pour le cas du cylindre, la puce électronique est entièrement isolée de l'extérieur : toutes les frontières sont de type wall, mis à part la condition de symétrie axiale sur la frontière de gauche. Nous avons également seulement considéré la moitié du domaine afin de simplifier la mise en place du maillage dans blockMesh.

 

Maillage du domaine englobant le composant électronique

 

Résultats de simulation OpenFOAM

 

Pour ce cas-test, nous avons choisi d'utiliser le solveur buoyantPimpleFoam. L'idée était de voir si ce solveur donne des résultats cohérents en comparaison avec le précédent solveur buoyantBoussinesqPimpleFoam.

Deux résultats de simulation ont pu être obtenus pour les modèles de turbulence laminar et k-$\epsilon$

 

                       laminar                                                                                       k-$\epsilon$

 

 

Les deux modèles fournissent un champ de température relativement comparable au niveau de la puce. Les différences entre ces deux modèles apparaissent sur le développement du panache de convection : dans le cas laminaire, celui-ci est fin et rectiligne, tandis que pour le modèle k-$\epsilon$ il se déforme rapidement sous l'effet de la turbulence.

3. Simulations StarCCM+

Géométrie du domaine et conditions aux limites

Schéma de la géométrie considérée et des conditions aux limites introduites dans StarCCM+

 

Maillage de la géométrie

De même que pour le cas-test précédent, nous avons maillé le domaine avec IcemCFD. Les détails du nombre de mailles et de la loi de maillage sont exposés sur la figure ci-dessous.

 

Maillage du domaine avec raffinement autour du composant électronique

 

Mise en place des simulations

Les modèles utilisés et toute la mise en place des simulations est identique à celle que nous avons pu faire dans le premier cas-test du cylindre.

Capture des modèles utilisés dans StarCCM+ pour la convection naturelle autour du composant électronique

 

Résultats

 

Résultats en régime laminaire

Résultat de simulation pour $T_P=348  K$ et $T_\infty = 293  K$  en régime laminaire

 

Tracé du nombre de Nusselt en fonction du nombre de Rayleigh

Étant donné que nous n'avons pas de corrélation pour ce cas, nous avons néanmoins jugé utile de tracer, comme pour le premier cas-test, le nombre de Nusselt en fonction du nombre de Rayleigh calculé pour différentes gravités : $g=9,81  m/s^2$, $g=98,1  m/s^2$, $g=981  m/s^2$ et $g=9810  m/s^2$.

 

Nombre de Nusselt numérique en fonction du nombre de Rayleigh calculé pour différentes gravités

 

Nous obtenons dans ce cas également une courbe du nombre de Nusselt croissante en fonction du nombre de Rayleigh.

 

Comparaison des modèles de turbulence

Tel que nous l'avons présenté dans le premier cas-test, nous avons procédé à des simulations pour divers modèles de turbulence afin de voir leur influence sur la couche limite thermique ainsi que sur le panache de convection. Les modèles k-$\omega$ et Reynolds Stress ne sont pas présentés ici car donnant des résultats peu exploitables.

 

Résultats de simulation pour $T_P=348  K$  et  $T_\infty=293 K$ avec différents modèles de turbulence

 

Nous pouvons constater des différences selon le modèles de turbulence utilisé, en particulier pour le modèle realizable k-$\epsilon$, contrairement au premier cas-test. En effet, les principales remarques concernent, ici également, le panache de convection dont la forme peut être plus évasée que pour le modèle laminaire, notamment pour le modèle standard k-$\epsilon$. Pour le modèle realizable k-$\epsilon$, le panache de convection oscille au cours du temps, ce qui diffère de manière significative par rapport aux autres modèles.

 

La figure ci-dessous présente les évolutions des nombres de Nusselt sur le composant pour différents modèles de turbulence. On peut constater l'adéquation entre les modèles tant du point de vue de la partie horizontale du composant que de ses parties verticales. Comme dans le premier cas-test, l'impact du choix du modèle de turbulence se voit principalement sur le panache de convection plutôt que sur la couche limite près du solide.

Évolution du nombre de Nusselt numérique en fonction de l'abscisse du composant pour différents modèles de turbulence

 

 

Synthèse

Dans cette partie, nous synthétisons les résultats que nous avons obtenus sur les différents cas d'étude.

 

Validation du logiciel OpenFOAM

 

 

Convection naturelle autour d'un cylindre horizontal

 

 

Conclusion et perspectives

A travers ce Bureau d'Etudes Industriel, nous avons pu nous familiariser avec le logiciel de CFD OpenFOAM qui nécessite une prise en main préalable afin d'être le plus efficace possible puisqu'il ne possède pas d'interface graphique. Il n'est en effet pas aisé de l'utiliser au premier abord, surtout lorsqu'on étudie des configurations, comme ce fut notre cas, qui nécessitent de modifier le code source et donc de comprendre son fonctionnement afin d'obtenir ce qui nous intéresse. Cependant, une fois cette difficulté appréhendée, on peut facilement concevoir les bénéfices qu'on peut tirer à utiliser OpenFOAM par la suite, que ce soit en terme de puissance de calcul ou de flexibilité.

 

Pour le cas du cylindre, une étude plus poussée permettrait de connaître réellement le comportement d'OpenFOAM vis-à-vis du phénomène de convection naturelle : extension de l'intervalle de Rayleigh ($Ra_D < 10^4$ et $Ra_D > 10^8$), agrandissement de la géométrie, ouverture du domaine en étudiant toutes les possibilités de conditions aux limites offertes par OpenFOAM.

Une autre amélioration consisterait à prendre en compte la conduction à l'intérieur du câble et du composant électronique, ce qui impliquerait de mailler ces solides et d'utiliser un autre solveur OpenFOAM (chtMultiRegionFoam par exemple). Il pourrait également être intéressant de tenir compte de la résistance de contact qui existe entre le composant électronique et son circuit imprimé dans le second cas-test. Par ailleurs, nous pourrions également procéder à des études où l'on considèrerait les transferts de chaleurs à travers plusieurs câbles et composants électroniques pour apprécier l'influence de la convection naturelle sur chacun d'eux.

 

 

 

Bibliographie

[1] Fundamentals of Heat and Mass Transfers - Frank P. Incropera, David P. DeWitt, Theodore L. Bergman et Adrienne S. Lavine - Sixth Edition

[2] Site web OpenFOAM

[3] Forced convection heat transfer -  Bill Garland, p 4.43 - p 4.48 : lien

 

Fundamentals of Heat and Mass Transfers