Accueil


 

Dimensionnement d'un évaporateur

Tubulaire pour le refroidissement

de composants électroniques d'un

satellite

 


Résumé

Les composants électroniques qui constituent un satellite sont sensibles aux fortes chaleurs ; leur température ne doit pas dépasser une température maximale. Il est donc nécessaire d'évacuer la chaleur qu'ils produisent et de les refroidir. La solution proposée par Thales Alenia Space consiste à faire circuler un fluide réfrigérant sous les composants électroniques à travers un évaporateur. La chaleur récupérée déclenche alors la vaporisation du fluide caloporteur sous le composant. L'efficacité du refroidissement provient du changement de phase. Cette méthode est plus efficace qu'un refroidissement classique en écoulement monophasique, car le flux de chaleur transféré est plus important. Plusieurs recherches approfondies ont été faites en ce qui concerne la modélisation de l'écoulement et des transferts de chaleur au sein de l'évaporateur. L'étude consiste à rechercher le modèle de simulation adéquat. La méthode consiste à tester tous les modèles empiriques susceptibles de répondre à nos besoins. A la suite de la programmation Matlab, il a été mis en évidence que le modèle thermique de Kandlikar était parfaitement adapté au problème en gravité normale et micro-gravité. Quant aux pertes de charge, il s'est avéré que le régime gravitationnel influençait leur modélisation. Ainsi, le modèle de Lockhart & Martinelli est le plus pertinent en gravité normale ; tandis qu'en micro-gravité, le modèle d'Awad est préférable. A la suite de ces résultats, les modèles choisis ont été implémentés dans un second programme Matlab afin de calculer les différents paramètres, comme les pertes de pression au sein de l'évaporateur, et la température des composants. Les résultats ont permis de vérifier que le cahier des charges était respecté dans les conditions d'utilisation spatiale.

 

Abstract

Two-phase thermal systems are broadly used in various industrial applications and engineering fields. These systems take advantage of latent heat transportation, which generally enables a good efficiency in heat exchanges. For that reason, two-phase thermal management systems are considered as extremely beneficial for space applications. Indeed, in satellites, the major thermal problem is currently to remove the vast heat amount generated by devices from the inside into the space, in order to ensure suitable environmental and working conditions. But boiling is a complex phenomenon which combines heat and mass transfers, hydrodynamics and interfacial phenomena. Furthermore, gravity consequently affects the fluid dynamics and may lead to unpredictable performances of thermal management systems. Information is needed to understand the thermal and hydraulic environment in the evaporating process. We need to predict the heat transfer coefficient and pressure drop. This informations can be obtained by simulations. Different models have been implemented in Matlab. First results show that the thermal model of Kandlikar is adapted to our problem both in normal gravity and microgravity. As for the pressure drops, gravity consequently affects the fluid dynamics. Then, they fit the Lockhart & Martinelli's correlation in normal gravity, whereas, in microgravity, the Awad's correlation is more suitable. Next step is to implement the selected model in a second Matlab program to calculate pressure drop and temperature of devices. Calculations show that spatial specifications are respected.


 

Dans le cadre de la troisième année de formation à l'ENSEEIHT, les élèves ont la chance de participer à un bureau d'étude industrielle de six semaines en partenariat avec de grandes compagnies. Ce projet permet à la fois de travailler en autonomie, mais aussi de découvrir une application industrielle ou un sujet de recherche particulier.

Pour notre part, nous avons choisi de travailler sur le système de refroidissement de composants électroniques d'un satellite. Il s'agit d'un sujet proposé par Thales Alenia Sapce et encadré par Catherine Colin de l'IMFT. Plus précisément, l'étude porte sur la modélisation de l'écoulement et des transferts de chaleur à l'intérieur de l'évaporateur du système de refroidissement.

 

                                  http://www.clubic.com/

L'Équipe

L'équipe de travail se constitue de deux élèves :

Les Encadrants

Ce projet est encadré par Mme Catherine COLIN, enseignant-chercheur à l'Institut de Mécanique des Fluides de Toulouse (IMFT) rattaché à l'ENSEEIHT. Tout au long du BEI, Mme COLIN aidera les étudiants et validera les résultats obtenus. De plus, Mr Julien HUGON ingénieur chargé du projet chez Thales Aliena Space, suivra également l'avancée du travail, fournira les renseignements utiles.

 

Introduction

 

Introduction

L’objectif de ce Bureau d’Études Industrielles est d'étudier les phénomènes thermo-hydrauliques au sein d'un évaporateur pour une boucle diphasique de refroidissement d'un satellite de télécommunication.

L'étude se décompose en 4 parties:

 

Thales Alenia Space

Thales Alenia Space est le numéro 1 européen des solutions par satellites et est un acteur majeur dans le domaine de l’infrastructure orbitale. Thales Alenia Space est aujourd’hui une référence mondiale dans le domaine des télécommunications, de la navigation, de la météorologie, de la gestion de l’environnement, de la défense et de la sécurité, de l’observation et de la science.

La branche de Thales Alenia Space présente un réseau de 11 sites industriels dans 4 pays européens et près de 7200 salariés dans le monde. Leader européen et numéro 3 mondial, Thales Alenia Space est à la pointe des technologies par satellites que ce soit dans le secteur civil ou militaire.

Aujourd’hui, Thales Alenia Space est le premier fournisseur de solutions par satellites en Europe dans le domaine de la défense et de la sécurité, avec de fortes positions dans les satellites et les systèmes sols. http://www.thalesgroup.com/space/

 

Contexte général

Les satellites sont constitués de systèmes électroniques qui dégagent de grandes quantités de flux de chaleur. Cette chaleur dégagée doit impérativement être évacuée, sous peine d'endommager les composants du satellite. La solution innovante proposée par Thales consiste à utiliser une boucle frigorifique pour refroidir les composants d'une étagère de satellite. La boucle est constituée d'un évaporateur, d'un surchauffeur de vapeur, d'un compresseur, d'un condenseur, et d'un détendeur. L'évaporateur est l'élément clés, c'est dans ce dernier que la chaleur dégagée par les composants électroniques est évacuée ; elle est transférée au fluide caloporteur. Le fluide caloporteur va donc s'échauffer et s'évaporer peu à peu. En résulte l'apparition d'un écoulement diphasique.

L'objectif de cette étude est de reprendre le projet réalisé l'année précédente et le compléter. Cette année, nous disposons de données expérimentales de l'IMFT. L'étude consiste à rechercher le modèle de simulation adéquat. La méthode consiste à tester tous les modèles empiriques susceptibles de répondre à nos besoins. Nous allons donc pouvoir comparer les données expérimentales aux simulations pour mettre en évidence lequel est le plus adapté. Après quoi, le modèle choisi sera implémenté dans un second programme Matlab afin de calculer les différents paramètres, comme les pertes de pression au sein de l'évaporateur, et la température des composants.

 

Schéma de la boucle de refroidissement.

Description du sujet

Les composants électroniques qui constituent un satellite sont sensibles aux fortes chaleurs ; leur température ne doit pas dépasser une température maximale. Il est donc nécessaire d'évacuer la chaleur qu'ils produisent et de les refroidir. Leur refroidissement n'est pourtant pas aisé : ceux situés aux extrémités du satellite évacuent leur chaleur avec l'air ambiant. Le satellite disposant de deux faces non exposées au soleil, les composants qui s'y situent évacuent leur chaleur par rayonnement. Le refroidissement du reste des composants électroniques est alors plus complexe.

La solution proposée par Thales Alenia Space consiste à faire circuler un fluide réfrigérant sous les composants électroniques à travers un évaporateur. La chaleur récupérée déclenche alors la vaporisation du fluide caloporteur sous le composant. L'efficacité du refroidissement provient du changement de phase. Cette méthode est plus efficace qu'un refroidissement classique en écoulement monophasique, car le flux de chaleur transféré est plus important.

Pour s'inscrire dans la continuité du travail de l'an passé, les modèles qui ont déjà été testés sont repris et complétés par de nouveaux de manière à avoir une liste exhaustive des modèles susceptibles de convenir au problème. De plus, nous reproduisons une géométrie exacte concernant la position des tubes, la localisation des imitateurs de composants électroniques et des coudes.

Pour satisfaire la requête de Thales Alenia Space, il faut premièrement déterminer le régime d'écoulement diphasique présent dans le circuit. Il en existe un éventail, mais dans notre cas, seul l'écoulement annulaire nous intéresse. De plus, il est nécessaire de connaître les pertes de charge au sein de l'écoulement, notamment les pertes de charge singulières dues aux nombreux coudes. Toute la difficulté est de modéliser ces pertes de charges car l'écoulement est diphasique. Pour cela, il faut se baser sur des études expérimentales déjà réalisées. Les frottements interfaciaux sont modélisés grâce au modèle de Wallis tandis que les frottements pariétaux avec les modèles de Lockhart & Martinelli, de Baroczy et d'Awad. En ce qui concerne les échanges thermiques convectifs dans le tube, il faut, là aussi, utiliser des corrélations empiriques diphasiques pour calculer le coefficient d'échange h qui varie au sein de l'écoulement. Les transferts de chaleur sous le régime d'ébullition saturée sont modélisés d'après les corrélations de Kandlikar, de Gunger & Winterton, de Schrock & Grossman et de Chen. Enfin, un modèle récemment développé par A. Cioncolini est aussi testé.

Dans un premier temps, tous les modèles cités plus haut sont comparés aux expériences de l'IMFT. Ces comparaisons sont effectuées en gravité et en micro-gravité, cela permet de sélectionner les plus performants en fonction du régime gravitationnel.

Dans un second temps, les modèles retenus sont intégrés dans un programme Matlab qui permet de calculer les pertes de pression au sein de l’évaporateur, la température des composants et de la paroi. Les résultats permettront de vérifier que le cahier des charges est respecté dans les conditions d'utilisation spatiales.

 

Prototype de l'évaporateur de la boucle de refroidissement.

Objectifs

L’objectif est de développer un programme MATLAB pour l’industriel Thales Alenia Space afin de calculer différents paramètres, comme les pertes de pression au sein de l’évaporateur, ou la température des composants avec les modèles les plus performants.

Après analyse des résultats, le programme pourra être corrigé, amélioré puis validé. Et les modèles pré-sélectionnés pour la simulations devront être triés en fonction de leurs performances.

Cahier des charges

Structure

Températures et débits

Géométrie

La géométrie est imposée, elle est fournie par l'industriel. Le circuit de refroidissement se compose de trois tubes parallèles fixés au support des composants.

Modélisation

Afin de modéliser au mieux les phénomènes thermo-hydrauliques, différents modèles diphasiques sont utilisés.

L'écoulement en entrée de l'évaporateur est supposé annulaire : la vapeur s'écoule au centre et le liquide est en paroi. La chaleur émise par les composants électroniques va initier un changement de phase : le liquide va peu à peu se vaporiser. Le titre massique de vapeur, noté x, va ainsi croître et la vitesse superficielle du gaz Jg va augmenter d'après l'équation :

$J_g = \frac{Gx}{\rho_g}$

G est le débit massique surfacique en kg/m2/s et ρg la masse volumique de la vapeur.

La vitesse du gaz étant supérieure à celle du liquide, il peut y avoir arrachement de gouttelettes de liquide à l'interface. Ce phénomène ne se produit qu'à partir d'une vitesse superficielle du gaz Jg critique. Pour se rendre compte de l'influence de l'arrachage sur la modélisation, deux modèles sont testés:

 

MODÈLES HYDRAULIQUES

Dans cette partie, l'écoulement est supposé annulaire sans arrachage. Les frottements interfaciaux τi sont modélisés grâce au modèle de Wallis tandis que les frottements pariétaux τp avec les modèles de Lockhart & Martinelli, de Baroczy et d'Awad. En ce qui concerne les pertes de charge singulières dans les coudes, un modèle issu d'expériences en laboratoire a été implémenté dans le programme.

Il faut toutefois rester vigilant à l'évolution du titre massique :

 

Modèle Annulaire sans arrachage

Lorsque le taux de vide est important, la phase gazeuse occupe le centre du tube tandis que la phase liquide se trouve sous la forme d'un film le long de la paroi du tube, c'est ce que l'on appelle un écoulement annulaire. En fonction des conditions imposées, il se peut également qu'une partie non négligeable du liquide se trouve sous forme de gouttelettes en suspension dans le cœur de vapeur ; l'écoulement est alors dit annulaire dispersé. Il s'agit du modèle annulaire avec arrachage qui est traité plus loin (paragraphe 5).

La connaissance de l’épaisseur de film, notée t, est importante pour l’estimation des transferts thermiques. Cette épaisseur du film diminue au fur et à mesure que la vitesse superficielle de la vapeur augmente.

Les équations qui régissent ce modèle s'écrivent comme suit :

$\frac{dR_g}{dz}G^2(\frac{R_l x^2}{\rho_g R_g ^2}+\frac{R_g (1-x)^2}{\rho_l R_l ^2}) = -\frac{\tau_{ig}4}{D}\sqrt{R_g}+R_g \frac{\tau_p 4}{D}-(\rho_l - \rho_g) R_g R_l g + G^2 \frac{dx}{dz}(\frac{2xR_l}{\rho_g R_g} + \frac{(1-x)(2 R_g -1)}{\rho_l R_l})$ [1]

 

$\frac{dp}{dz} = \frac{-d}{dz}(\frac{G^2 x^2}{\rho_g R_g})-\frac{d}{dz}(\frac{G^2 (1-x)^2}{\rho_l R_l}) + \frac{\tau_p 4}{D} - (\rho_l R_l + \rho_g R_g) g$ [2]

 

$\frac{dx}{dz} = \frac{4q}{G D h_lv}$ [3]

 

$t = \frac{D}{2}(1-\sqrt{R_g})$

L’équation [1] est obtenue en éliminant le gradient de pression entre les équations de quantité de mouvement pour le liquide et la vapeur. L’équation [2] est l’équation de quantité d emouvement pour le mélange et l’équation [3], le bilan d’enthalpie.

Le taux de vide, les pertes de charge et le titre sont ainsi déterminés le long de l'évaporateur, pour chaque pas ∆z.

NB : Ces équations sont utilisées lorsque le fluide réfrigérant est au niveau d’un composant. Lorsque ce n’est pas le cas, l’écoulement est adiabatique et il n'y aucun changement de phase. Le taux de vide Rg et le titre massique x sont constants et les termes faisant intervenir les dérivées de x et Rg sont nuls dans les équations [1] et [2].

 

Frottement interfacial

Le Modèle de Wallis [2] permet de modéliser les frottements interfaciaux τi entre la phase liquide et la phase gazeuse. Il est défini comme suit: $\tau_i = -\frac{1}{2} f_i \rho_g |U_g - U_l|(U_g -U_l)$ avec  $f_i = 0.005 (1+300 \frac{\delta}{D}) = 0.005 [1 + 150(1-\sqrt{R_g})]$.

Les vitesses du gaz Ug et du liquide Ul sont déterminées à partir des vitesses superficielles Jg et Jl calculées à chaque pas spatial de la façon :

$U_g = \frac{J_g}{R_g}$,  $U_l = \frac{J_l}{R_l}$ , avec  $R_g + R_l = 1$   et  $J_g = \frac{G x}{\rho_g}$ , $J_l = \frac{G (1-x)}{\rho_l}$.

Frottement pariétal

Plusieurs modélisations des frottements pariétaux sont envisageables:

Modèle de Lockhart & Martinelli

Le modèle de Lockhart et Martinelli [2] permet de déterminer le frottement pariétal τp. Ce modèle est l'un des plus vieux et des plus utilisés pour calculer les pertes de charges diphasiques par frottement. Cette corrélation a été établie pour des écoulements adiabatiques en tube de section cylindrique. Il utilise le paramètre de Martinelli χ, qui compare la dissipation des deux phases.

Étape 1 : On calcule les pertes de charge dans la phase liquide et la phase gazeuse telles que :

$(\frac{dp}{dz})_l = -\frac{S_p}{A}f_{pl}\frac{\rho_l J_l ^2}{2}$ , $(\frac{dp}{dz})_g = -\frac{S_p}{A}f_{pg}\frac{\rho_g J_g ^2}{2}$.

Les coefficients de frottement pariétal sont définis de la manière suivante :

$f_{pl} = K(\frac{J_l D}{\nu_l})^{-n}$ , $f_{pg} = K(\frac{J_g D}{\nu_g})^{-n}$

Avec :

Les nombres de Reynolds sont basés sur les vitesses superficielles : $Re_l = \frac{J_l D}{\nu_l}$ , $Re_g = \frac{J_g D}{\nu_g}$.

Étape 2 : On calcule le paramètre de Martinelli χ, tel que:

$X = \sqrt{\frac{(dp/dz)_l}{(dp/dz)_g}} = \frac{J_l}{J_g}\sqrt{\frac{\rho_l f_pl}{\rho_g f_pg}}$

Étape 3 : On calcule les coefficients multiplicatifs φl et φg tels que :

$\varphi_l ^2 = (1+\frac{C}{X}+\frac{1}{X^2})$ , $\varphi_g ^2 = (1+ CX+X^2)$.

 

Le paramètre C étant fonction du régime d'écoulement :

Paramètre de Lockhart & Martinelli.
Liquide Gaz C
Turbulent Turbulent 20
Laminaire Turbulent 12
Turbulent Laminaire 10
Laminaire Laminaire 5

Le régime laminaire est défini pour un Reynolds inférieur à 2000 et le régime turbulent pour un Reynolds supérieur à 3000. Dans la zone de transition, il n'y a donc pas de corrélation. Dans les simulations, le changement de régime occasionne des discontinuités dans l'expression de la contrainte de frottement diphasique. Le problème a été résolu en utilisant le lissage proposé par Erik de Malmazet, post-doctorant ayant travaillé sur le sujet à l'IMFT. Cette interpolation linéaire de C est une fonction de log(Rel) et/ou log(Reg). Lorsqu'une seule des deux phases est dans la zone de transition (cas le plus fréquemment rencontré), l'interpolation linéaire de C est de type C(Rel)=a*log(Rel)+b ou C(Reg)=a'*log(Reg)+b'.

A titre d'exemple, dans le cas où le gaz est en régime turbulent et le liquide en régime de transition, la corrélation utilisée est :

$C(Re_l)= \frac{20-12}{log (3000)-log(2000)}*(log( Re_l)-log(2000))+12$

Dans le cas où les deux phases sont en zone de transition, on utilise deux interpolations linéaires de C du type C(Rel, Reg)=a*log(Rel)+b*log(Reg)+c, une pour le cas Rel>Reg et l'autre pour le cas Rel<Reg.

Étape 4 : On peut alors calculer les pertes de charge par frottement :

$\left ( \frac{dp}{dz} \right )_{fr}=\frac{\tau_p S_p}{A}={\varphi_l}^2 \left ( \frac{dp}{dz} \right )_{l}={\varphi_g}^2 \left ( \frac{dp}{dz} \right )_{g}$.

Modèle de Baroczy

La corrélation de Baroczy [3] consiste à exprimer les pertes de charge par frottement diphasique en fonction des pertes de charge par frottement (dP/dz)l que l’on aurait si la phase liquide s’écoulait seule avec le même débit masse qu'en écoulement diphasique. Pour cela, il faut pouvoir exprimer le coefficient multiplicatif  défini de la manière suivante :

${\varphi_{l0}}^2=\frac{\left ({dp}/{dz} \right )_{fr}}{\left ( {dp}/{dz} \right )_l}$

Pour cela, on utilise aussi les pertes de charge par frottement (dP/dz)g que l’on aurait si la phase vapeur s’écoulait seule avec le même débit de masse que l’écoulement diphasique et on introduit la variable Y définie de la manière suivante :

$Y^2=\frac{\left ({dp}/{dz} \right )_g}{\left ( {dp}/{dz} \right )_l}= \frac{1}{X^2}$

Les pertes de charge (dP/dz)l et (dP/dz)g sont définies de la même manière que dans le modèle de Lockhart et Martinelli.

Les deux phases s’écoulant seules, avec le même débit masse qu'en écoulement diphasique, auraient respectivement les vitesses moyennes  $V_l=G/\rho_l$ et $V_g=G/\rho_g$, à partir desquelles on peut définir les nombres de Reynolds  $Re_{V l}$ et $Re_{V g}$ . On a alors :

$Y^2=\frac{\rho_g}{\rho_l} \frac{f_{pg}}{f_{pl}}\frac{{J_{g}}^2}{{J_{l}}^2}$

Baroczy donne une expression de ${\varphi_{l0}}$ à partir de nombreux points expérimentaux pour différents fluides. La corrélation qu’il obtient est graphique mais ses courbes ont par la suite été corrélées par Chisholm en 1973, qui a obtenu l’expression suivante :

${\varphi_{l0}}^2=1+(Y^2-1) [B x^{\frac{2-n}{2}} (1-x)^{\frac{2-n}{2}}+x^{2-n}]$

n est l’opposé de l’exposant de  que l’on a dans f($Re_{V l}$), soit :

B est une variable dépendant de Y et de G dont l’expression varie en fonction de la valeur de Y :

Paramètre de Baroczy.
Valeur de Y $0<Y<9,5$ $9,5 <Y<28$ $28<Y$
Valeur de B $\frac{55}{G^{1/2}}$ $\frac{520}{YG^{1/2}}$ $\frac{15000}{Y^2 G^{1/2}}$

Finalement, on peut calculer les pertes de charge par frottement :

$\left ( \frac{dp}{dz} \right )_{fr}= {\varphi_{l0}}^2 \left ( \frac{dp}{dz} \right )_{l}$

Modèle d'Awad

Le modèle d'Awad [4] est un autre modèle empirique qui permet de déterminer le frottement pariétal. On calcule le nombre de Reynolds de chaque phase tel que :

$Re_l=\frac{J_l D}{\nu_l} $ et $Re_g=\frac{J_g D}{\nu_g} $

Les coefficients de frottement pariétal liquide fpl et gazeux fpg sont définis de la même manière que dans le modèle de Lockhart et Martinelli.

On calcule alors le paramètre X, de la même façon que le paramètre de Martinelli :

$X=\frac{J_l}{J_g} \sqrt{\frac{\rho_l f_{pl}}{\rho_g F_{pg}}}$

Les définitions de $\varphi_l$ et $\varphi_g$ changent :

${\varphi_l}^2=\left ( 1+ \left (  \frac{1}{X^2} \right )^p \right )^{1/p}$  et  ${\varphi_l}^2=\left ( 1+ \left (  \frac{1}{X^2}\right )^p \right )^{1/p}$

 

La valeur du paramètre p est difficile à déterminer. Elle varie en fonction des cas et généralement, elle est choisie comme le minimum de l’erreur RMS. Pour assurer des conditions de robustesse, p = 2/7. On peut alors calculer le gradient de pression liquide et gazeux.

$\left ( \frac{dp}{dz} \right )_{l}=-\frac{S_p}{A}f_{pl}\frac{\rho_l {J_l}^2}{2}$ et  $\left ( \frac{dp}{dz} \right )_{d}=-\frac{S_p}{A}f_{pg}\frac{\rho_g {J_g}^2}{2}$

Le gradient de pression par frottement s'écrit alors :

$\left ( \frac{dp}{dz} \right )_{fr}=\frac{\tau_p S_p}{A}={\varphi_l}^2 \left ( \frac{dp}{dz} \right )_{l}={\varphi_g}^2 \left ( \frac{dp}{dz} \right )_{g}$

Finalement, le modèle d’Awad reprend les mêmes équations que le modèle de Lockhart et Martinelli. Les seuls paramètres qui changent sont les facteurs diphasiques $\varphi_l$ et $\varphi_g$.

 

Modèle Monophasique

Lorsque le titre de vapeur x atteint l'unité, c'est le modèle monophasique qui est utilisé. Le cahier des charges imposant en sortie un titre massique x proche de 1, il est normal de s'intéresser au cas d'un modèle à un fluide pour calculer les pertes de charge.

$\frac{dp}{dz}=\frac{4 \tau_p}{D}$ avec $\tau_p =-0,5f_c \rho_g {U_g}^2$.

Le coefficient de frottement pariétal fc dépend du régime d'écoulement dans lequel on se situe :$\left\{\begin{array}{l l} f_{c} = 16/Re_g & \quad \text{si $Re_g < 2000$}\\f_{c} = 0,079 Re_g ^{0,25} & \quad \text{si $Re_g > 2000$}\end{array} \right.$   

 

Perte de charge dans les coudes

La littérature [5] renseigne un modèle empirique pour les pertes de charge :

$\Delta p_{rb}= \Phi \Delta p_{sp}$

$\Delta p_{rb}$ est la perte de charge totale à travers un coude pour un écoulement diphasique et

$\Delta p_{sp}$ défini la perte de charge totale à travers un coude en écoulement monophasique. Elle s'écrit : $\Delta p_{sp}=K_{sp} \frac{G^2}{2 \rho_l}$.

Ksp est le coefficient de perte de charge local pour un écoulement monophasique liquide. Pour estimer sa valeur [Idelcik 1986] suggère l'expression suivante :

$K_{sp}=f_{pl} \frac{L}{D}+0,294 \left ( \frac{R}{D} \right )^{0,5}$

où R est le rayon de courbure et fpl le coefficient de frottement pariétal calculé par les modèles explicités précédemment.

Le coefficient multiplicateur pour l’écoulement diphasique ϕ est donné par : $\Phi=1+ \left ( \frac{\rho_l}{\rho_g} -1 \right ) x [b(1-x)+x]$ avec $b=1+ \frac{2,2}{K_{sp}(2+ \frac{R}{D})}$.

 

MODÈLES THERMIQUES

Du point de vue thermique, il s'agit de modéliser les transferts de chaleur. Le type de corrélation utilisé dépend tout d'abord du régime d'ébullition, lié au taux de vide et plus généralement aux conditions d'entrée.

Pour le régime d'ébullition saturée, le liquide en entrée est à la température de saturation et son taux de vide n'est pas nul. On observe une forte évolution longitudinale de l'écoulement, ce dernier pouvant être à bulles, poches bouchons et annulaire. Mais on ne s’intéresse qu'au régime annulaire.

L'ensemble des flux dégagés par les composants électroniques est connu au sein de la géométrie ; ils sont tous constants. Cette chaleur est dissipée grâce à la vaporisation du fluide circulant dans l'évaporateur. L'objectif est de déterminer la température des différents composants.

En régime d'ébullition nucléée saturée, on suppose que le fluide dans l'évaporateur est à sa température de saturation Tl= TSAT , et que l'intégralité du flux thermique sert à vaporiser le fluide réfrigérant. De plus, on suppose que les phases vapeur et liquide sont en équilibre thermodynamique lorsque un changement d'état se produit au niveau du fluide. On peut alors accéder à l'évolution du titre massique le long de l'évaporateur :

$\frac{dx}{dz}= \frac{qS_p}{AGh_{lg}D}$

Pour chaque pas d'espace Δz, on peut ainsi calculer le titre massique x puis en déduire le coefficient d'échange thermique le long du tube. Le coefficient d'échange thermique est déterminé d'après les corrélations de Kandlikar, de Gunger & Winterton, de Schrock & Grossman et de Chen ; toutes décrites ci-dessous.

Calcul du coefficient d'échange thermique h

Notations et définitions

On introduit les grandeurs suivantes :

 

Corrélations de Kandlikar

La corrélation de Kandlikar [2] exprime le coefficient d'échange thermique comme suit :

$h=h_l \left [ C_1 {C_0}^{C_2} (25Fr)^{C_5}+C_3 Bo^{C_4}F_k \right ]$

Pour déterminer les différentes constantes de l'équation, il est nécessaire de s'intéresser à l'écoulement dans la phase liquide. Pour cela, il faut calculer le nombre d'ébullition Bo, le nombre de Froude Fr et la constante C0 définie telle que   $C_0= \left ( \frac{1-x}{x} \right )^{0,8} \sqrt{\frac{\rho_g}{\rho_l}}$.

Les valeurs des constantes Ci sont déduites de C0 :

Paramètre de Kandlikar
  C0<0,65 C0>0,65
C1 1,1360 0,6683
C2 -0,9 -0,2
C3 667,2 1058
C4 0,7 0,7
C5 0,3 0,3

C5=0 pour les tubes verticaux et les tubes horizontaux quand Fr>0.04. Pour le fluide R245FA, Fk = 1,4.

 

Corrélation de Gunger & Winterton

Selon le modèle de Gunger et Winterton [2] $h=h_l \left [ 1+3000Bo^{0,86}+ \left ( \frac{x}{1-x} \right )^{3/4} \left ( \frac{\rho_l}{\rho_g} \right )^{0,41} \right ]$

 

Corrélation de Schrock & Grossman

Selon le modèle de Schrock et Grossman [2] $h=7,39.10^3 h_l \left [ Bo+0,00015 \left ( \frac{1}{X_{tt}} \right )^{0,66} \right ]$

 

Corrélation de Chen

Chen [2] a proposé la première corrélation pour de l'évaporation en tube vertical pour une large gamme de validité. Il a considéré le coefficient de transfert thermique de l'écoulement diphasique h comme la somme du coefficient de transfert thermique lié au phénomène de nucléation et au phénomène de convection.

Dans son modèle, il a supposé que le gradient de température imposé par les conditions de convection forcée supprimait une partie des sites de nucléation, réduisant ainsi le coefficient de transfert thermique associé à ce phénomène. Plus le titre augmente, plus la turbulence se développe favorisant ainsi les transferts de chaleur.

Le coefficient de transfert thermique peut ainsi être modélisé de la sorte : $h=Sh_n+Fh_l$ avec :

 

Modèle monophasique (x=1)

Lorsque le titre x tend vers un, l'écoulement est considéré comme monophasique vapeur. Le bilan thermique local s'écrit alors:GC_{pg} \frac{dT_g}{dz}=q \frac{S_p}{A}$. Tl correspondant à la température du liquide.

Pour un flux imposé: $T_g(z)-T_{ge}= \frac{q S_p}{AGC_{pg}}(z-z_e)$ et  $T_p(z)-T_{g}(z)= \frac{q}{h_g}$.

Le coefficient d'échange thermique est donné par : $h_g=0,023 \frac{\lambda_g}{D} [ \ \left (\frac{GD}{\mu_g} \right ) \ ] ^{0,8} Pr^{1/3}$

Calcul des températures

Maintenant que l'évolution du coefficient de transfert thermique est connue le long de l'évaporateur, la température en paroi et celle des différents composants peuvent être déterminées. Une schématisation simplifiée des résistances thermiques permet de relier ces deux températures. Les tubes cylindriques sont attachés à des ailettes. Ces ailettes sont-elles-mêmes en contact avec les composants en question. Les composants électroniques sont imités par des résistances qui dissipent de l'énergie de manière homogène. Les différentes résistances prises en compte sont les suivantes :

La figure ci-dessous représente le circuit de résistances associé au problème. Certains composants sont parcourus plusieurs fois par les tubes (n*3 tubes), dans ce cas, le sens de l'écoulement n'est pas le même: la gravité change de signe.

Par une analogie électrique, les températures de la paroi, des ailettes et des composants peuvent être déterminées à partir de celle du fluide.

La méthode est la suivante:

$T_p=T_{sat}+ \frac{q}{h}$ en diphasique; ou  $T_p=T_{sat}+ \frac{q S_p}{AGC_{pg} }\Delta z +\frac{q}{h_g} $ en monophasique;

$T_{ailette}=T_p + \frac{3 \pi n_{passages} DR_{cond}Q}{Ll^2}$

$T_{composant}=T_{ailette}+R_{contact} \frac{Q}{Ll}$

où  q est le flux de chaleur surfacique échangé à la paroi en W/m²,

Q est la puissance dissipée par l'élément en W,

l la largeur de l'ailette en m,

L la longueur du composant en m et

nPASSAGE le nombre de passages du tube.

 

L'ensemble des corrélations présentées jusqu'ici constitue la théorie pour un écoulement annulaire sans arrachage. Autrement dit, la phase gazeuse occupe le centre du tube tandis que la phase liquide se trouve sous la forme d'un film le long de la paroi. Dans le modèle de A. Cioncolini qui est présenté ci-dessous, il est fait l’hypothèse qu'une partie du liquide en paroi est arrachée sous forme de gouttelettes.

 

MODELE DE A. CIONCOLINI

Il s'agit d'un modèle annulaire avec arrachage [6]. Il est applicable aux écoulements annulaires et prédit l'évolution des valeurs du coefficient d'échange thermique , des pertes de charges, du taux d’entraînement, de l'épaisseur du film liquide et du taux de vide.

Le taux d'entraînement, noté e, est calculé comme une fonction du nombre de Weber associé à l'écoulement gazeux au cœur de la conduite : $e=(1+279,6{We_c}^{-0,8395)^{-2,209}} $ pour 10¹<WeC<10⁵. Avec $We_c = \frac{\rho_C {J_g}^2 d}{\sigma}$ où Jg est la vitesse superficielle du gaz Ug et $\rho_c$ la masse volumique du coeur définie comme $\rho_C = \frac{x+e(1-x)}{x/\rho_g +e(1-x)/\rho_l}$. [7]

 

Le taux de vide Rg lui est défini comme suit: $R_g =\frac{Kx^n}{1+(K-1)x^n}$ avec $\left\{\begin{array}{l l} K=a+(1-a)(\rho_g / \rho_l)^{\alpha} \\n=b+(1-b)(\rho_g /\rho_l)^{\beta} \end{array} \right.$ [8].

 

Le gradient de pression totale, c'est à dire qui fait intervenir l'accélération, la gravité et le frottement s'écrit:

$\frac{dp}{dz}= \frac{4 \tau_p}{D} - G^2\frac{d}{dz} \left [ \frac{x^2}{\rho_g R_g} \frac{e^2(1-x)^2}{\rho_l \gamma (1-R_g)} + \frac{(1-e)^2(1-x)^2}{\rho_l (1- \gamma) (1-R_g)} \right] - [\rho_l (1-R_g)+\rho_g R_g]g$

 

Le frottement pariétal  $\tau_p$  s'écrit : $\tau_p= \frac{1}{2} f_{tp} \rho_C V_C^2$ . Où le coefficient de frottement est déterminé par la corrélation de Fanning $f_{pl}=0,172 {We_c}^{-0,372}$ pour WeC>100.

Avec dC le diamètre interne du coeur:  $d_C=D \sqrt{R_g \frac{x \rho_l +e(1-x) \rho_g}{x \rho_l}}$ .

VC est la vitesse du coeur diphasique:  $V_C=\frac{4}{\pi} \frac{\left[ x+e(1-x) \right] \dot {m} }{\rho_C d_C^2}$ [9].

 

L'épaisseur du film liquide t est quant à elle définie comme   $t=y^{*} max \left(\sqrt{\frac{2 \Gamma_{lf}^{+}}{R^{+}}} ; 0,066 \frac{\Gamma_{lf}^{+}}{R^{+}} \right)$ [10].

Où $y^{*}=\frac{\mu_l}{\rho_l V^{"*"}}$ est l'échelle de longueur;

      $V^{*}=\sqrt{\frac{\tau_p}{\rho_l}}$ est l'échelle de vitesse;

      $\Gamma_{lf}^{+}=\frac{(1-e)(1-x) \dot {m}}{2 \pi \rho_l V^{*}{y^{*}}^2}$ est le débit dans le film liquide adimensionnalisé;

      $R^{+}=\frac{R}{y^{*}}$ es le rayon du tube adimensionnalisé.

 

Enfin, le coefficient d'échange thermique s'écrit:  $h=\frac{\lambda_l}{t} 0,0776 t^{0,9}{Pr_l^{0,52}}t$  avec  $t^{+}=\frac{t}{y^{*}}$.

 

Les résultats de ce modèle sont comparés à ceux fournis par le modèle annulaire sans arrachage et aux données expérimentales de l'IMFT afin de déterminer lequel est le plus adapté au problème.

Modèles hydrauliques

MODÈLES HYDRAULIQUES

Dans cette partie, l'écoulement est supposé annulaire sans arrachage. Les frottements interfaciaux τi sont modélisés grâce au modèle de Wallis tandis que les frottements pariétaux τp avec les modèles de Lockhart & Martinelli, de Baroczy et d'Awad. En ce qui concerne les pertes de charge singulières dans les coudes, un modèle issu d'expériences en laboratoire a été implémenté dans le programme.

Il faut toutefois rester vigilant à l'évolution du titre massique :

  • Lorsque la valeur de x devient importante (proche de 1), l'écoulement est alors considéré comme monophasique et il faut basculer vers des modèles monophasiques en transferts de chaleur et en hydraulique.

  • Inversement lorsque le titre x devient faible, les vitesses moyennes du liquide et du gaz sont du même ordre de grandeur Ul ~Ug, et le modèle homogène doit être appliqué. Ce modèle n'est cité qu'à titre informatif, car dans le cas présent, le titre est supérieur ou égale 0,3 ; il n'est pas question de l'utiliser.

 

Modèle Annulaire sans arrachage

Lorsque le taux de vide est important, la phase gazeuse occupe le centre du tube tandis que la phase liquide se trouve sous la forme d'un film le long de la paroi du tube, c'est ce que l'on appelle un écoulement annulaire. En fonction des conditions imposées, il se peut également qu'une partie non négligeable du liquide se trouve sous forme de gouttelettes en suspension dans le cœur de vapeur ; l'écoulement est alors dit annulaire dispersé. Il s'agit du modèle annulaire avec arrachage qui est traité plus loin (paragraphe 5).

La connaissance de l’épaisseur de film, notée t, est importante pour l’estimation des transferts thermiques. Cette épaisseur du film diminue au fur et à mesure que la vitesse superficielle de la vapeur augmente.

Les équations qui régissent ce modèle s'écrivent comme suit :

$\frac{dR_g}{dz}G^2(\frac{R_l x^2}{\rho_g R_g ^2}+\frac{R_g (1-x)^2}{\rho_l R_l ^2}) = -\frac{\tau_{ig}4}{D}\sqrt{R_g}+R_g \frac{\tau_p 4}{D}-(\rho_l - \rho_g) R_g R_l g + G^2 \frac{dx}{dz}(\frac{2xR_l}{\rho_g R_g} + \frac{(1-x)(2 R_g -1)}{\rho_l R_l})$ [1]

 

$\frac{dp}{dz} = \frac{-d}{dz}(\frac{G^2 x^2}{\rho_g R_g})-\frac{d}{dz}(\frac{G^2 (1-x)^2}{\rho_l R_l}) + \frac{\tau_p 4}{D} - (\rho_l R_l + \rho_g R_g) g$ [2]

 

$\frac{dx}{dz} = \frac{4q}{G D h_lv}$ [3]

 

$t = \frac{D}{2}(1-\sqrt{R_g})$

L’équation [1] est obtenue en éliminant le gradient de pression entre les équations de quantité de mouvement pour le liquide et la vapeur. L’équation [2] est l’équation de quantité d emouvement pour le mélange et l’équation [3], le bilan d’enthalpie.

Le taux de vide, les pertes de charge et le titre sont ainsi déterminés le long de l'évaporateur, pour chaque pas ∆z.

NB : Ces équations sont utilisées lorsque le fluide réfrigérant est au niveau d’un composant. Lorsque ce n’est pas le cas, l’écoulement est adiabatique et il n'y aucun changement de phase. Le taux de vide Rg et le titre massique x sont constants et les termes faisant intervenir les dérivées de x et Rg sont nuls dans les équations [1] et [2].

 

Frottement interfacial

Le Modèle de Wallis [2] permet de modéliser les frottements interfaciaux τi entre la phase liquide et la phase gazeuse. Il est défini comme suit: $\tau_i = -\frac{1}{2} f_i \rho_g |U_g - U_l|(U_g -U_l)$ avec  $f_i = 0.005 (1+300 \frac{\delta}{D}) = 0.005 [1 + 150(1-\sqrt{R_g})]$.

Les vitesses du gaz Ug et du liquide Ul sont déterminées à partir des vitesses superficielles Jg et Jl calculées à chaque pas spatial de la façon :

$U_g = \frac{J_g}{R_g}$,  $U_l = \frac{J_l}{R_l}$ , avec  $R_g + R_l = 1$   et  $J_g = \frac{G x}{\rho_g}$ , $J_l = \frac{G (1-x)}{\rho_l}$.

Frottement pariétal

Plusieurs modélisations des frottements pariétaux sont envisageables:

Modèle de Lockhart & Martinelli

Le modèle de Lockhart et Martinelli [2] permet de déterminer le frottement pariétal τp. Ce modèle est l'un des plus vieux et des plus utilisés pour calculer les pertes de charges diphasiques par frottement. Cette corrélation a été établie pour des écoulements adiabatiques en tube de section cylindrique. Il utilise le paramètre de Martinelli χ, qui compare la dissipation des deux phases.

Étape 1 : On calcule les pertes de charge dans la phase liquide et la phase gazeuse telles que :

$(\frac{dp}{dz})_l = -\frac{S_p}{A}f_{pl}\frac{\rho_l J_l ^2}{2}$ , $(\frac{dp}{dz})_g = -\frac{S_p}{A}f_{pg}\frac{\rho_g J_g ^2}{2}$.

Les coefficients de frottement pariétal sont définis de la manière suivante :

$f_{pl} = K(\frac{J_l D}{\nu_l})^{-n}$ , $f_{pg} = K(\frac{J_g D}{\nu_g})^{-n}$

Avec :

  • n = 1 et K=16 en écoulement laminaire ;

  • n = ¼ et K=0.079 en écoulement turbulent.

Les nombres de Reynolds sont basés sur les vitesses superficielles : $Re_l = \frac{J_l D}{\nu_l}$ , $Re_g = \frac{J_g D}{\nu_g}$.

Étape 2 : On calcule le paramètre de Martinelli χ, tel que:

$X = \sqrt{\frac{(dp/dz)_l}{(dp/dz)_g}} = \frac{J_l}{J_g}\sqrt{\frac{\rho_l f_pl}{\rho_g f_pg}}$

Étape 3 : On calcule les coefficients multiplicatifs φl et φg tels que :

$\varphi_l ^2 = (1+\frac{C}{X}+\frac{1}{X^2})$ , $\varphi_g ^2 = (1+ CX+X^2)$.

 

Le paramètre C étant fonction du régime d'écoulement :

Paramètre de Lockhart & Martinelli.
Liquide Gaz C
Turbulent Turbulent 20
Laminaire Turbulent 12
Turbulent Laminaire 10
Laminaire Laminaire 5

Le régime laminaire est défini pour un Reynolds inférieur à 2000 et le régime turbulent pour un Reynolds supérieur à 3000. Dans la zone de transition, il n'y a donc pas de corrélation. Dans les simulations, le changement de régime occasionne des discontinuités dans l'expression de la contrainte de frottement diphasique. Le problème a été résolu en utilisant le lissage proposé par Erik de Malmazet, post-doctorant ayant travaillé sur le sujet à l'IMFT. Cette interpolation linéaire de C est une fonction de log(Rel) et/ou log(Reg). Lorsqu'une seule des deux phases est dans la zone de transition (cas le plus fréquemment rencontré), l'interpolation linéaire de C est de type C(Rel)=a*log(Rel)+b ou C(Reg)=a'*log(Reg)+b'.

A titre d'exemple, dans le cas où le gaz est en régime turbulent et le liquide en régime de transition, la corrélation utilisée est :

$C(Re_l)= \frac{20-12}{log (3000)-log(2000)}*(log( Re_l)-log(2000))+12$

Dans le cas où les deux phases sont en zone de transition, on utilise deux interpolations linéaires de C du type C(Rel, Reg)=a*log(Rel)+b*log(Reg)+c, une pour le cas Rel>Reg et l'autre pour le cas Rel<Reg.

Étape 4 : On peut alors calculer les pertes de charge par frottement :

$\left ( \frac{dp}{dz} \right )_{fr}=\frac{\tau_p S_p}{A}={\varphi_l}^2 \left ( \frac{dp}{dz} \right )_{l}={\varphi_g}^2 \left ( \frac{dp}{dz} \right )_{g}$.

Modèle de Baroczy

La corrélation de Baroczy [3] consiste à exprimer les pertes de charge par frottement diphasique en fonction des pertes de charge par frottement (dP/dz)l que l’on aurait si la phase liquide s’écoulait seule avec le même débit masse qu'en écoulement diphasique. Pour cela, il faut pouvoir exprimer le coefficient multiplicatif  défini de la manière suivante :

${\varphi_{l0}}^2=\frac{\left ({dp}/{dz} \right )_{fr}}{\left ( {dp}/{dz} \right )_l}$

Pour cela, on utilise aussi les pertes de charge par frottement (dP/dz)g que l’on aurait si la phase vapeur s’écoulait seule avec le même débit de masse que l’écoulement diphasique et on introduit la variable Y définie de la manière suivante :

$Y^2=\frac{\left ({dp}/{dz} \right )_g}{\left ( {dp}/{dz} \right )_l}= \frac{1}{X^2}$

Les pertes de charge (dP/dz)l et (dP/dz)g sont définies de la même manière que dans le modèle de Lockhart et Martinelli.

Les deux phases s’écoulant seules, avec le même débit masse qu'en écoulement diphasique, auraient respectivement les vitesses moyennes  $V_l=G/\rho_l$ et $V_g=G/\rho_g$, à partir desquelles on peut définir les nombres de Reynolds  $Re_{V l}$ et $Re_{V g}$ . On a alors :

$Y^2=\frac{\rho_g}{\rho_l} \frac{f_{pg}}{f_{pl}}\frac{{J_{g}}^2}{{J_{l}}^2}$

Baroczy donne une expression de ${\varphi_{l0}}$ à partir de nombreux points expérimentaux pour différents fluides. La corrélation qu’il obtient est graphique mais ses courbes ont par la suite été corrélées par Chisholm en 1973, qui a obtenu l’expression suivante :

${\varphi_{l0}}^2=1+(Y^2-1) [B x^{\frac{2-n}{2}} (1-x)^{\frac{2-n}{2}}+x^{2-n}]$

n est l’opposé de l’exposant de  que l’on a dans f($Re_{V l}$), soit :

  • n = 1 si l’écoulement liquide seul est laminaire ( $Re_{V l}$< 2000)

  • n = 0.25 s’il est turbulent ( $Re_{V l}$> 2000).

B est une variable dépendant de Y et de G dont l’expression varie en fonction de la valeur de Y :

Paramètre de Baroczy.
Valeur de Y $0<Y<9,5$ $9,5 <Y<28$ $28<Y$
Valeur de B $\frac{55}{G^{1/2}}$ $\frac{520}{YG^{1/2}}$ $\frac{15000}{Y^2 G^{1/2}}$

Finalement, on peut calculer les pertes de charge par frottement :

$\left ( \frac{dp}{dz} \right )_{fr}= {\varphi_{l0}}^2 \left ( \frac{dp}{dz} \right )_{l}$

Modèle d'Awad

Le modèle d'Awad [4] est un autre modèle empirique qui permet de déterminer le frottement pariétal. On calcule le nombre de Reynolds de chaque phase tel que :

$Re_l=\frac{J_l D}{\nu_l} $ et $Re_g=\frac{J_g D}{\nu_g} $

Les coefficients de frottement pariétal liquide fpl et gazeux fpg sont définis de la même manière que dans le modèle de Lockhart et Martinelli.

On calcule alors le paramètre X, de la même façon que le paramètre de Martinelli :

$X=\frac{J_l}{J_g} \sqrt{\frac{\rho_l f_{pl}}{\rho_g F_{pg}}}$

Les définitions de $\varphi_l$ et $\varphi_g$ changent :

${\varphi_l}^2=\left ( 1+ \left (  \frac{1}{X^2} \right )^p \right )^{1/p}$  et  ${\varphi_l}^2=\left ( 1+ \left (  \frac{1}{X^2}\right )^p \right )^{1/p}$

 

La valeur du paramètre p est difficile à déterminer. Elle varie en fonction des cas et généralement, elle est choisie comme le minimum de l’erreur RMS. Pour assurer des conditions de robustesse, p = 2/7. On peut alors calculer le gradient de pression liquide et gazeux.

$\left ( \frac{dp}{dz} \right )_{l}=-\frac{S_p}{A}f_{pl}\frac{\rho_l {J_l}^2}{2}$ et  $\left ( \frac{dp}{dz} \right )_{d}=-\frac{S_p}{A}f_{pg}\frac{\rho_g {J_g}^2}{2}$

Le gradient de pression par frottement s'écrit alors :

$\left ( \frac{dp}{dz} \right )_{fr}=\frac{\tau_p S_p}{A}={\varphi_l}^2 \left ( \frac{dp}{dz} \right )_{l}={\varphi_g}^2 \left ( \frac{dp}{dz} \right )_{g}$

Finalement, le modèle d’Awad reprend les mêmes équations que le modèle de Lockhart et Martinelli. Les seuls paramètres qui changent sont les facteurs diphasiques $\varphi_l$ et $\varphi_g$.

 

Modèle Monophasique

Lorsque le titre de vapeur x atteint l'unité, c'est le modèle monophasique qui est utilisé. Le cahier des charges imposant en sortie un titre massique x proche de 1, il est normal de s'intéresser au cas d'un modèle à un fluide pour calculer les pertes de charge.

$\frac{dp}{dz}=\frac{4 \tau_p}{D}$ avec $\tau_p =-0,5f_c \rho_g {U_g}^2$.

Le coefficient de frottement pariétal fc dépend du régime d'écoulement dans lequel on se situe :$\left\{\begin{array}{l l} f_{c} = 16/Re_g & \quad \text{si $Re_g < 2000$}\\f_{c} = 0,079 Re_g ^{0,25} & \quad \text{si $Re_g > 2000$}\end{array} \right.$   

 

Perte de charge dans les coudes

La littérature [5] renseigne un modèle empirique pour les pertes de charge :

$\Delta p_{rb}= \Phi \Delta p_{sp}$

$\Delta p_{rb}$ est la perte de charge totale à travers un coude pour un écoulement diphasique et

$\Delta p_{sp}$ défini la perte de charge totale à travers un coude en écoulement monophasique. Elle s'écrit : $\Delta p_{sp}=K_{sp} \frac{G^2}{2 \rho_l}$.

Ksp est le coefficient de perte de charge local pour un écoulement monophasique liquide. Pour estimer sa valeur [Idelcik 1986] suggère l'expression suivante :

$K_{sp}=f_{pl} \frac{L}{D}+0,294 \left ( \frac{R}{D} \right )^{0,5}$

où R est le rayon de courbure et fpl le coefficient de frottement pariétal calculé par les modèles explicités précédemment.

Le coefficient multiplicateur pour l’écoulement diphasique ϕ est donné par : $\Phi=1+ \left ( \frac{\rho_l}{\rho_g} -1 \right ) x [b(1-x)+x]$ avec $b=1+ \frac{2,2}{K_{sp}(2+ \frac{R}{D})}$.

Modèles thermiques

 

MODÈLES THERMIQUES

Du point de vue thermique, il s'agit de modéliser les transferts de chaleur. Le type de corrélation utilisé dépend tout d'abord du régime d'ébullition, lié au taux de vide et plus généralement aux conditions d'entrée.

Pour le régime d'ébullition saturée, le liquide en entrée est à la température de saturation et son taux de vide n'est pas nul. On observe une forte évolution longitudinale de l'écoulement, ce dernier pouvant être à bulles, poches bouchons et annulaire. Mais on ne s’intéresse qu'au régime annulaire.

L'ensemble des flux dégagés par les composants électroniques est connu au sein de la géométrie ; ils sont tous constants. Cette chaleur est dissipée grâce à la vaporisation du fluide circulant dans l'évaporateur. L'objectif est de déterminer la température des différents composants.

En régime d'ébullition nucléée saturée, on suppose que le fluide dans l'évaporateur est à sa température de saturation Tl= TSAT , et que l'intégralité du flux thermique sert à vaporiser le fluide réfrigérant. De plus, on suppose que les phases vapeur et liquide sont en équilibre thermodynamique lorsque un changement d'état se produit au niveau du fluide. On peut alors accéder à l'évolution du titre massique le long de l'évaporateur :

$\frac{dx}{dz}= \frac{qS_p}{AGh_{lg}D}$

Pour chaque pas d'espace Δz, on peut ainsi calculer le titre massique x puis en déduire le coefficient d'échange thermique le long du tube. Le coefficient d'échange thermique est déterminé d'après les corrélations de Kandlikar, de Gunger & Winterton, de Schrock & Grossman et de Chen ; toutes décrites ci-dessous.

Calcul du coefficient d'échange thermique h

Notations et définitions

On introduit les grandeurs suivantes :

  • le nombre d'ébullition : $Bo=\frac{q}{Gh_{lg}}$
  • le nombre de Froude : $Fr=\frac{G^2}{{\rho_l}^2gD}$

  • le coefficient de transfert thermique de la phase liquide déterminé à partir de la corrélation de Dittus Boetler (1930) pour des écoulements turbulents: $h_l=0,023 \frac{\lambda_l}{D}{Re_l}^{0,8}Pr^{1/3}$

  • le Reynolds liquide basé sur le flux de masse et le taux de vide : $Re_l= \frac{G(1-x)D}{\mu_l}$

  • le Prandtl liquide : $Pr=\frac{\mu_l C_{pl}}{\lambda_l}$

  • le paramètre de Martinelli :  $X_{tt}= \frac{1-x}{x}  \sqrt{  \frac{\rho_g f_{pl}}{\rho_l f_{pg}} }$  et

 

Corrélations de Kandlikar

La corrélation de Kandlikar [2] exprime le coefficient d'échange thermique comme suit :

$h=h_l \left [ C_1 {C_0}^{C_2} (25Fr)^{C_5}+C_3 Bo^{C_4}F_k \right ]$

Pour déterminer les différentes constantes de l'équation, il est nécessaire de s'intéresser à l'écoulement dans la phase liquide. Pour cela, il faut calculer le nombre d'ébullition Bo, le nombre de Froude Fr et la constante C0 définie telle que   $C_0= \left ( \frac{1-x}{x} \right )^{0,8} \sqrt{\frac{\rho_g}{\rho_l}}$.

Les valeurs des constantes Ci sont déduites de C0 :

Paramètre de Kandlikar
  C0<0,65 C0>0,65
C1 1,1360 0,6683
C2 -0,9 -0,2
C3 667,2 1058
C4 0,7 0,7
C5 0,3 0,3

C5=0 pour les tubes verticaux et les tubes horizontaux quand Fr>0.04. Pour le fluide R245FA, Fk = 1,4.

 

Corrélation de Gunger & Winterton

Selon le modèle de Gunger et Winterton [2] $h=h_l \left [ 1+3000Bo^{0,86}+ \left ( \frac{x}{1-x} \right )^{3/4} \left ( \frac{\rho_l}{\rho_g} \right )^{0,41} \right ]$

 

Corrélation de Schrock & Grossman

Selon le modèle de Schrock et Grossman [2] $h=7,39.10^3 h_l \left [ Bo+0,00015 \left ( \frac{1}{X_{tt}} \right )^{0,66} \right ]$

 

Corrélation de Chen

Chen [2] a proposé la première corrélation pour de l'évaporation en tube vertical pour une large gamme de validité. Il a considéré le coefficient de transfert thermique de l'écoulement diphasique h comme la somme du coefficient de transfert thermique lié au phénomène de nucléation et au phénomène de convection.

Dans son modèle, il a supposé que le gradient de température imposé par les conditions de convection forcée supprimait une partie des sites de nucléation, réduisant ainsi le coefficient de transfert thermique associé à ce phénomène. Plus le titre augmente, plus la turbulence se développe favorisant ainsi les transferts de chaleur.

Le coefficient de transfert thermique peut ainsi être modélisé de la sorte : $h=Sh_n+Fh_l$ avec :

  • hn le coefficient de transfert thermique associé au phénomène de nucléation. On utilise pour cela la corrélation d'ébullition en cuve de Forster et Zuber (1955) :

    $h_n=0,00122 \left [ \frac{{k_l}^{0,79}{C_{pl}}^{0,45}{\rho_l}^{0,49}}{{\sigma}^{0,5}{h_{lg}}^{0,24}{\mu_l}^{0,29}{\rho_g}^{0,24}} \right ] (T_p-T_{sat})^{0,24}( P_{sat}(T_p)-P_l)^{0,75} $

  • S le facteur de suppression des sites de nucléation : $S=\frac{1}{1+2,53.10^{-6}(Re_lF(X_{tt})^{1,25})^{1,17}}$

  • F le coefficient multiplicateur qui représente l'augmentation du coefficient de transfert thermique par rapport à un écoulement monophasique : $F(X_{tt})=2,35 \left [ 0,213+ \frac{1}{X_{tt}} \right ]^{0,736}$  pour  ${X_{tt}}^{-1}> 0,1 $, sinon $F(X_{tt})=1$.

 

Modèle monophasique (x=1)

Lorsque le titre x tend vers un, l'écoulement est considéré comme monophasique vapeur. Le bilan thermique local s'écrit alors:GC_{pg} \frac{dT_g}{dz}=q \frac{S_p}{A}$. Tl correspondant à la température du liquide.

Pour un flux imposé: $T_g(z)-T_{ge}= \frac{q S_p}{AGC_{pg}}(z-z_e)$ et  $T_p(z)-T_{g}(z)= \frac{q}{h_g}$.

Le coefficient d'échange thermique est donné par : $h_g=0,023 \frac{\lambda_g}{D} [ \ \left (\frac{GD}{\mu_g} \right ) \ ] ^{0,8} Pr^{1/3}$

Calcul des températures

Maintenant que l'évolution du coefficient de transfert thermique est connue le long de l'évaporateur, la température en paroi et celle des différents composants peuvent être déterminées. Une schématisation simplifiée des résistances thermiques permet de relier ces deux températures. Les tubes cylindriques sont attachés à des ailettes. Ces ailettes sont-elles-mêmes en contact avec les composants en question. Les composants électroniques sont imités par des résistances qui dissipent de l'énergie de manière homogène. Les différentes résistances prises en compte sont les suivantes :

  • Une résistance de convection au niveau de l'échange fluide / paroi, déterminée par des modèles de coefficient d'échange thermique.

  • Une résistance de conduction entre la paroi des tubes et les ailettes, calculée en fonction du coefficient h et de l'ordre de 2.10-4 K.m2/W.

  • Une résistance de contact entre la semelle de l'évaporateur et l'imitateur qui vaut 2.10-4 K.m2/W.

La figure ci-dessous représente le circuit de résistances associé au problème. Certains composants sont parcourus plusieurs fois par les tubes (n*3 tubes), dans ce cas, le sens de l'écoulement n'est pas le même: la gravité change de signe.

Par une analogie électrique, les températures de la paroi, des ailettes et des composants peuvent être déterminées à partir de celle du fluide.

La méthode est la suivante:

$T_p=T_{sat}+ \frac{q}{h}$ en diphasique; ou  $T_p=T_{sat}+ \frac{q S_p}{AGC_{pg} }\Delta z +\frac{q}{h_g} $ en monophasique;

$T_{ailette}=T_p + \frac{3 \pi n_{passages} DR_{cond}Q}{Ll^2}$

$T_{composant}=T_{ailette}+R_{contact} \frac{Q}{Ll}$

où  q est le flux de chaleur surfacique échangé à la paroi en W/m²,

Q est la puissance dissipée par l'élément en W,

l la largeur de l'ailette en m,

L la longueur du composant en m et

nPASSAGE le nombre de passages du tube.

 

Modèle de A. Cioncolini

L'ensemble des corrélations présentées jusqu'ici constitue la théorie pour un écoulement annulaire sans arrachage. Autrement dit, la phase gazeuse occupe le centre du tube tandis que la phase liquide se trouve sous la forme d'un film le long de la paroi. Dans le modèle de A. Cioncolini qui est présenté ci-dessous, il est fait l’hypothèse qu'une partie du liquide en paroi est arrachée sous forme de gouttelettes.

 

MODELE DE A. CIONCOLINI

Il s'agit d'un modèle annulaire avec arrachage [6]. Il est applicable aux écoulements annulaires et prédit l'évolution des valeurs du coefficient d'échange thermique , des pertes de charges, du taux d’entraînement, de l'épaisseur du film liquide et du taux de vide.

Le taux d'entraînement, noté e, est calculé comme une fonction du nombre de Weber associé à l'écoulement gazeux au cœur de la conduite : $e=(1+279,6{We_c}^{-0,8395)^{-2,209}} $ pour 10¹<WeC<10⁵. Avec $We_c = \frac{\rho_C {J_g}^2 d}{\sigma}$ où Jg est la vitesse superficielle du gaz Ug et $\rho_c$ la masse volumique du coeur définie comme $\rho_C = \frac{x+e(1-x)}{x/\rho_g +e(1-x)/\rho_l}$. [7]

 

Le taux de vide Rg lui est défini comme suit: $R_g =\frac{Kx^n}{1+(K-1)x^n}$ avec $\left\{\begin{array}{l l} K=a+(1-a)(\rho_g / \rho_l)^{\alpha} \\n=b+(1-b)(\rho_g /\rho_l)^{\beta} \end{array} \right.$ [8].

 

Le gradient de pression totale, c'est à dire qui fait intervenir l'accélération, la gravité et le frottement s'écrit:

$\frac{dp}{dz}= \frac{4 \tau_p}{D} - G^2\frac{d}{dz} \left [ \frac{x^2}{\rho_g R_g} \frac{e^2(1-x)^2}{\rho_l \gamma (1-R_g)} + \frac{(1-e)^2(1-x)^2}{\rho_l (1- \gamma) (1-R_g)} \right] - [\rho_l (1-R_g)+\rho_g R_g]g$

 

Le frottement pariétal  $\tau_p$  s'écrit : $\tau_p= \frac{1}{2} f_{tp} \rho_C V_C^2$ . Où le coefficient de frottement est déterminé par la corrélation de Fanning $f_{pl}=0,172 {We_c}^{-0,372}$ pour WeC>100.

Avec dC le diamètre interne du coeur:  $d_C=D \sqrt{R_g \frac{x \rho_l +e(1-x) \rho_g}{x \rho_l}}$ .

VC est la vitesse du coeur diphasique:  $V_C=\frac{4}{\pi} \frac{\left[ x+e(1-x) \right] \dot {m} }{\rho_C d_C^2}$ [9].

 

L'épaisseur du film liquide t est quant à elle définie comme   $t=y^{*} max \left(\sqrt{\frac{2 \Gamma_{lf}^{+}}{R^{+}}} ; 0,066 \frac{\Gamma_{lf}^{+}}{R^{+}} \right)$ [10].

Où $y^{*}=\frac{\mu_l}{\rho_l V^{"*"}}$ est l'échelle de longueur;

      $V^{*}=\sqrt{\frac{\tau_p}{\rho_l}}$ est l'échelle de vitesse;

      $\Gamma_{lf}^{+}=\frac{(1-e)(1-x) \dot {m}}{2 \pi \rho_l V^{*}{y^{*}}^2}$ est le débit dans le film liquide adimensionnalisé;

      $R^{+}=\frac{R}{y^{*}}$ es le rayon du tube adimensionnalisé.

 

Enfin, le coefficient d'échange thermique s'écrit:  $h=\frac{\lambda_l}{t} 0,0776 t^{0,9}{Pr_l^{0,52}}t$  avec  $t^{+}=\frac{t}{y^{*}}$.

 

Les résultats de ce modèle sont comparés à ceux fournis par le modèle annulaire sans arrachage et aux données expérimentales de l'IMFT afin de déterminer lequel est le plus adapté au problème.

Démarche scientifique

Démarche scientifique

La démarche est la suivante : puisque nous disposons de données expérimentales pour le taux de vide, le titre massique et le coefficient d'échange thermique en micro-gravité et gravité normale, nous allons confronter ces données aux résultats fournis par les diverses corrélations recueillies. De cette étude préliminaire, nous pourrons alors sélectionner les modèles les plus performants. Cependant, avant de les implémenter au programme Matlab destiné au projet de Thales, il faut s'assurer que les conditions d'utilisation sont similaires. Il est donc nécessaire de faire une étude adimensionnelle. Si les nombres adimensionnels du projet Thales sont du même ordres de grandeurs que ceux des expériences de l'IMFT, alors les modèles présélectionnés devraient être applicables au projet Thales.

Le problème fait intervenir quatre dimensions (M, L, T, θ) et treize variables, à savoir :

Finalement, neuf nombres adimensionnels sont nécessaires pour définir le problème. Les 9 nombres adimensionnels en question sont les suivants :

$Re_l= \frac{G J_l}{\nu_l}$      $Re_g=\frac{G J_g}{\nu_g}$      $Fr_l=\frac{G^2}{{\rho_l}^2 g D}$      $Fr_g=\frac{G^2}{{\rho_g}^2 g D}$     $We_l=\frac{\rho_l D {J_l}^2}{\sigma}$     $We_g=\frac{\rho_g D {J_g}^2}{\sigma}$      $Bo=\frac{q}{G h_{lg}}$

Les valeurs de ces nombres sont résumés dans le tableau ci-dessous.

Étude adimensionnelle.
  IMFT 1-g THALES  1-g IMFT μ-g THALES  μ-g
ρgl 0,0061-0,0117 0,0117 0,0053-0,0095 0,0117
μgl  0,0337-0,0359 0,0329 0,0317-0,39 0,0329
λgl 0,1785-0,1829 0,185 0,1767-0,1848 0,185
Rel 301-6340 10-2260 146-6510 4-2260
Reg 6900-66000 29000-98000 7068-88100 29500-98100
Wel 0,13-60 6,5.10-6-3 0,032-60 10-5-3
Weg 13-1050 49-540 15-2050 50-545
Frl 0,076-1,44 0,04 / /
Frg 1530-36600 293 / /
Bo 9,39.10-5-0,003 0-2,94.10-4 1,12.10-4-0,0028 0-2,94.10-4

En général, les ordres de grandeurs des expériences sont plutôt proches de celles du projet. Il n'y a que le nombre de Froude qui fasse défaut, mais le débit, élevé au carré, diffère beaucoup d'une expérience à l'autre (G~88kg/m²/s pour Thales, G~[100;400] pour l'IMFT). En conséquence, les conclusions tirées des données de l'IMFT peuvent être appliquées au projet de Thales Alenia Space.

Choix des modèles

Comparaisons des modèles avec les données de l'IMFT

Dans cette partie, il s'agit de comparer les résultats fournis par les corrélations citées précédemment avec les données expérimentales de l'IMFT. L’ébullition nucléée en micro-gravité a été étudiée grâce à des expériences réalisées en vols paraboliques en avion. Des régimes d’ébullition saturée stationnaires à flux constant ont été obtenus pour un réfrigérant HFE7000 à différentes pressions. Les expériences ont permis de relever la pression et la température du liquide, le titre massique, le taux de vide le frottement pariétal et le coefficient d'échange h, pour plusieurs débits. Quelques résultats sont décrits dans les figures 1, 2 et 3.

En ce qui concerne l'épaisseur de film liquide en paroi, elle dépend fortement du titre massique puisqu'elle décroit significativement quand le titre augmente, quelque soit le régime gravitationnel. Il semble également que l'épaisseur augmente avec le débit mais cela n'est pas flagrant. Par contre, l'influence de la gravité est plus nette : le film liquide est plus fin en micro-gravité. Cela s'explique par le fait que la gravité accumule le liquide en bas du tube, en résulte un film plus épais en gravité.

Quant au coefficient d'échange, pour un débit et un flux de chaleur fixés, il augmente avec le titre massique. Le coefficient de transfert thermique augmente aussi avec le débit. Il semble qu'il soit moins important en micro-gravité qu'en gravité normal pour des titres faibles, x<0,2. Inversement, le coefficient h est plus important en micro-gravité qu'en gravité normale pour les titres plus élevés, x>0,2.

Dans le cas du projet Thales, le titre est supérieur ou égale 0,3, on en conclut que le transfert thermique est détérioré en gravité normale.

Figure 1: Épaisseur du film en fonction du titre massique en μ-g et 1-g.

 

 

Figure 2: Coefficient de transfert thermique en 1-g.

 

 

Figure 3: Coefficient de transfert thermique en μ-g.

 

Discutons maintenant des modélisations. Toutes les simulations ont été effectuées en micro-gravité et en gravité normale de manière à comparer les différents résultats et ne retenir que les modèles les plus performants en fonction du régime gravitationnel. Les différents modèles de frottement pariétal et de coefficient d'échange sont confrontés aux données expérimentales en micro-gravité et gravité.

 

Frottement pariétal

Dans un premier temps, ce sont les différents modèles de frottement pariétal qui sont confrontés aux expériences. Connaissant les débits et les titres associés, les corrélations de Lockhart & Martinelli, de Baroczy et d'Awad ont pu être appliquées. Ce sont les coefficients multiplicatifs φl puis les frottements interfaciaux τp de chacun des modèles qui sont comparés ; les résultats figurent dans les illustrations 4 à 7. Le taux d'erreur a également été calculé, cf. Tableaux 1 et 2.

L'erreur relative moyenne est définie comme suit : $\frac{1}{n}\sum_{i=1}^{n} \frac{\varphi_{i exp}-\varphi_{i th}}{\varphi_{i exp}}$  où ψ est la grandeur étudiée.

L'écart type est défini  $\sqrt{\frac{1}{n}\sum_{i=1}^{n}(\varphi_{i exp}-\varphi_{i th})^2}$  comme où ψ est la grandeur étudiée.

En gravité normale, on remarque que les expériences sont en très bon accord avec le modèle de Lockhart & Martinelli.

En micro-gravité, c'est le modèle d'Awad qui est le plus performant. Le modèle de Lockhart & Martinelli a de bons résultats mais ils sont dispersés. Quant au modèle de Baroczy, il sous-estime le coefficient multiplicatif, ils sont donc tous les deux donc à exclure de la modélisation.

En conclusion, ce sont les modèles de Lockhart & Martielli et d'Awad qui sont retenus pour calculer les pertes de charge diphasiques, respectivement en gravité normale et micro-gravité, dans le programme Matlab.

Tableaux 1: Erreur relative moyenne pour les modèles de pertes de charges.
  1-g μ-g
Awad 38,15% 8,26%
Baroczy 56,92% 62,65%
Lockhart & Martinelli 5,52% -17,73%
Tableaux 2: Écart-type pour les modèles de pertes de charges.
  1-g μ-g
Awad 3,49 2,39
Baroczy 2,83 8,71
Mockhart & Martinelli 2,13 3,36

Figure 4: Coefficient multiplicatif en 1-g.

 

 

Figure 5: Coefficient multiplicatif en μ-g.

 

 

Figure 6: Frottement pariétal en 1-g.

 

 

Figure 7: Frottement pariétal en μ-g.

 

 

Coefficient d'échange thermique

Ce sont maintenant les différents modèles de coefficient d'échange qui sont comparés aux expériences. De la connaissance des débits surfaciques, des titres associés et des flux surfaciques de chaleur, on déduit les valeurs de h selon les corrélations de Kandlikar, de Gunger & Winterton, de Schrock & Grossman et de Chen. Les résultats sont résumés dans les figures 8 et 9. Le taux d'erreur relative et l'écart moyen ont également été calculés, cf. Tableau 4 et 5.

En gravité normale, c'est le modèle de Kandlikar qui approxime au mieux les expériences ; le modèle de Gunger & Winterton n'est pas mauvais, mais les résultats sont plus dispersés relativement à ceux de Kandlikar. Les modèles de Chen et de Schrock & Grossman ont tendance à sous-estimer le coefficient d'échange.

En micro-gravité, ce sont les modèles de Kandlikar, de Gunger & Winterton et de Chen les plus performants. Mais d'après les taux d'erreurs, celui de Kandlikar reste le meilleur. Le modèle de Schrock & Grossman sous-estime toujours le coefficient h.

En conclusion, ce sera le modèle de Kandlikar qui modélisera le coefficient d'échange h pour les deux régimes gravitationnels.

Tableau 4: Erreur relative moyenne pour les corrélations du coefficient d'échange.
  1-g μ-g
Chen 23,57% 7,33%
Gunger & Winterton 14,03% 6,13%
Kandlikar 9,35% -4,1%
Schrock & Grossman 39,35% 38,37%
Tableau 5: Écart-type pour les corrélations du coefficient d'échange.
  1-g μ-g
Chen 974 1044
Gunger & Winterton 758 936
Kandlikar 547 786
Schrock & Grossman 1562 1728

 

Figure 8: Coefficient de transfert thermique comparé aux expériences en 1-g.

 

 

Figure 9: Coefficient de transfert thermique comparé aux expériences en μ-g.

 

 

Bilan

On est finalement en mesure de sélectionner les modèles les plus performants en fonction du régime gravitationnel. Le tableau 6 résume la situation.

Tableau 6: Sélection des modèles en fonction de leurs performances et du régime gravitationnel.
  Frottement interfacial Frottement pariétal Coefficient d'échange Thermique
μ-g Wallis Lockhart & Martinelli Kandlikar
1-g Wallis Awad Kandlikar

 

 

Comparaisons du modèle de A. Cioncolini avec les données de l'IMFT

Il s'agit maintenant de comparer les résultats fournis par le modèle de A. Cioncolini [3] avec les expériences de l'IMFT. En effet, le Dr A. Cioncolini a étudié les écoulements annulaires avec arrachage et a développé avec son équipe un modèle. Il l'avait implémenter dans Matlab pour différents débits de fluide HFE-7000, le même que celui utilisé par l'IMFT. L'objectif est donc de comparer ses résultats avec les expériences. On pourra alors estimer l'influence de l'arrachage. Il est important de remarquer que le modèle de A. Cioncolini ne dépend ni de la gravité, ni du flux de chaleur.

 

Taux d'entraînement

La figure 10 représente le taux d'entraînement. Il augmente de manière significative avec le débit surfacique G.

Cela s'explique par le fait que la vitesse superficielle du gaz Jg augmente avec ce dernier selon la relation :

$J_g=\frac{Gx}{\rho_g}$

La vitesse du gaz étant supérieure à celle du liquide, il y a alors arrachement de gouttelettes de liquide à l'interface.

Le taux d'entraînement n'a pas été mesuré dans l'étude de l'IMFT, c'est la raison pour laquelle aucune donnée expérimentale ne figure sur le graphique.

Figure 10: Taux d'entraînement selon A. Cioncolini.

 

Taux de vide

Le taux de vide selon le modèle de A. Cioncolini ne dépend pas du débit. On remarque sur les figures 11 et 12, qu'il est plus adapté aux expériences en micro-gravité. En effet, les écarts sont élevés en gravité normale. Il est important de noter le taux de vide intervient dans le calcul du frottement pariétal, l'épaisseur du film liquide et le coefficient d'échange h ; ces erreurs vont donc se répercuter sur ces derniers.

Figure 11: Taux de vide comparé au modèle de A. Cioncolini, en 1-g

 

 

Figure 12: Taux de vide comparé au modèle de A. Cioncolini, en μ-g.

 

Frottement pariétal

Selon les figures 13 et 14, le frottement pariétal augmente avec le débit et le titre massique. De plus, il semble que les données expérimentales soient plus proches du modèle de A. Cioncolini pour des débits faibles, donc pour de faibles taux d'entraînement. On remarque aussi que pour de faibles titres, le frottement pariétal est sous estimé avec le modèle de A. Cioncolini.

 

Figure 13: Frottement pariétal comparé au modèle de A. Cioncolini, en 1-g.

 

 

Figure 14: Frottement pariétal comparé au modèle de A. Cioncolini, en μ-g.

 

Épaisseur du film liquide

A débit fixé, l'épaisseur du film liquide diminue avec le titre massique. Il diminue également avec le débit, c'est en accord avec le fait que le taux d'entraînement croit avec le débit.

Comme remarqué plus haut, les écarts sont, plus faibles en micro-gravité. Le modèle de A. Cioncolini sous-estime les épaisseurs en gravité normale, notamment pour des titres élevés.

Il est important de remarquer que le modèle de A. Cioncolini prédit une diminution de l'épaisseur avec l'augmentation du débit ; tandis que l'IMFT prédit le contraire. Cependant, le modèle de A. cioncolini prédit également une augmentation de l'arrachage avec le débit ; il est donc normal de constater une diminution de l'épaisseur avec une augmentation du débit.

Figure 15: Épaisseur du film liquide comparée au modèle de A. Cioncolini, en 1-g.

 

 

Figure 16: Épaisseur du film liquide comparée au modèle de A. Cioncolini, en μ-g.

 

Coefficient d'échange thermique

Le modèle de A. Cioncolini prédit une augmentation du coefficient d'échange h avec le titre et le débit. Il ne privilégie pas les expériences en micro-gravité : les deux configurations gravitationnelles présentent des écarts du même ordre de grandeur. Cependant, d'après les figures 17 et 18, ces écarts sont d'autant plus grands que les débits sont importants.

Figure 17: Coefficient de transfert thermique comparé au modèle de A. Cioncolini, en 1-g.

 

 

Figure 18: Coefficient de transfert thermique comparé au modèle de A. Cioncolini, en μ-g.

 

 

Prolongement

Pour approfondir le modèle de A. Cioncolini, nous l'avons aussi étudier en supposant qu'il n'y avait aucun arrachage donc pour un taux d’entraînement nul, e=0. Les écarts entre le modèle avec et sans entraînement sont illustrés en figures 19 à 26; ils sont discutés dans ce qui suit.

On remarque que les écarts sont d'autant plus importants que le débit est élevé. Pour un débit de 100 kg/m²/s, il est nul. Pour des débits faibles, G<100, l'arrachage est donc négligeable.

Cependant, on peut se poser la question à savoir lequel des deux modèles, avec ou sans arrachage, est le plus adéquat. Le sujet est traité dans le paragraphe qui suit.

Remarque : Dans le cas de l'évaporateur, le débit est de l'ordre 88 kg/m²/s, on conclut que l'arrachage y est négligeable. En conclusion, l'hypothèse qui annule le taux d'entraînement dans le modèle de A. Cioncolini n'est pas nécessaire pour obtenir des résultats pertinents vis à vis de l'évaporateur.

Figure 19: Résultats du modèle de A. Cioncolini sans arrachage.

 

 

Figure 20: Taux de vide comparée au modèle de A. Cioncolini sans arrachage, en 1-g.

 

 

Figure 21: Taux de vide comparée au modèle de A. Cioncolini sans arrachage, en μ-g.

 

 

Figure 22: Frottement pariétal comparé au modèle de A. Cioncolini sans arrachage, en 1-g.

 

 

Figure 23: Frottement pariétal comparé au modèle de A. Cioncolini sans arrachage, en μ-g.

 

 

Figure 24: Épaisseur du film liquide comparée au modèle de A. Cioncolini sans arrachage, en 1-g.

 

 

Figure 25: Épaisseur du film liquide comparée au modèle de A. Cioncolini sans arrachage, en μ-g.

 

 

Figure 26: Coefficient de transfert thermique comparé au modèle de A. Cioncolini sans arrachage, en 1-g.

 

 

Figure 26: Coefficient de transfert thermique comparé au modèle de A. Cioncolini sans arrachage, en μ-g.

 

Conclusion: Avec ou sans arrachage

Maintenant que les deux régimes d'écoulement, avec et sans arrachage, ont été étudiés, il s'agit d'identifier lequel est le plus adapté à la problématique. En gravité normale, c'est le modèle de Lockhart & Martinelli et celui de Kandlikar qui avaient été sélectionnés, respectivement pour les pertes de charges et le coefficient d'échange. Tandis qu'en micro-gravité, ce sont les modèles d'Awad et de Kandlikar qui ont été choisis. Nous allons donc les comparer au modèle de A. Cioncolini (e≠0 et e=0).

Les résultats sont résumés dans les figures 27 à 30 et les erreurs dans le tableau 7.

 

Tableau 7: Erreur relative moyenne et écart-type pour les modèles avec et sans arrachage.
  1-g μ-g
Pertes de charge

Lockhart & Martinelli:   5,52%   2,13

Cioncolini e≠0:            23,32%  3,04

Cioncolini e=0:               31%      3,9

Awad:                     8,26%     2,39

Cioncolini e≠0:   -31,74%    2,40

Cioncolini e=0:   -19,58%    3,40

Coefficient h

Kandlikar:                     9,35%    547

Cioncolini e≠0:           23,25%   977

Cioncolini e=0:          27,78%   1174

Kandlikar:               -4,1%      782

Cioncolini e≠0:      9,94%      994

Cioncolini e=0:      -1,77%   1102

 

Au vue des résultats, nous sommes en mesure de recommander l'utilisation du modèle de Kandlikar pour le coefficient d'échange quelque soit le régime gravitationnel. En ce qui concerne les pertes de charge, le modèle de Lockart & Martinelli est préférable en gravité normale ; en micro-gravité, ce sera le modèle d'Awad.

Le modèle de A. Cioncolini est finalement inutilisable dans notre cas puisque l'arrachage est négligeable. En effet, il avait été remarqué que ses résultats étaient plus précis uniquement pour des taux d’entraînement faibles.

Le tableau 8 résume les conclusions.

Tableau 8: Bilan des modèles sélectionnés.
  1-g μ-g
Pertes de charge Lockhart & Martinelli Awad
Coefficient h Kandlikar Kandlikar

 

Figure 27: Comparaison des modèles de Lockhart & Martinelli et de Cioncolini, en 1-g.

 

 

Figure 28: Comparaison des modèles de Awad  et de Cioncolini, en μ-g.

 

 

Figure 29: Comparaison des modèles de Kandlikar et de Cioncolini, en 1-g.

 

 

Figure 30: Comparaison des modèles de Kandlikar et de Cioncolini, en μ-g.

Résultats des simulations

Résultats des simulations

A la suite de ces conclusions, un second programme Matlab a été élaboré afin de calculer la température des composants et les pertes de charges le long de l'évaporateur. En réalité, il s'agit du programme développé par le groupe de l'an dernier qui a été amélioré et corrigé.

L'évolution de ces paramètres est illustrée dans les figures 31 à 36. Les simulations ont été réalisées en micro-gravité et gravité normale. La première partie du tube est ascendante.

On remarque, en figure 31, que le titre augmente bien le long de l'évaporateur au fur et à mesure des refroidissements des composants, ; il atteint l'unité à la fin du tube. Dans le programme, la condition x=1 a été remplacée par x=0,99 pour éviter les divisions par zéro. Au delà de x=0,99, le titre est considéré constant, l'écoulement est supposé monophasique.

En ce qui concerne le taux de vide, il est lui aussi limité à 0,99 pour des raisons pratiques. La figure 32 décrit son évolution. En gravité normale, dans les tubes descendants, le gaz est ralenti (Ug diminue) donc le taux de vide augmente. Inversement, dans les tubes ascendants, le liquide est accumulé au fond, le gaz s'écoule sans gène, sa vitesse augmente et le taux de vide diminue. Par contre, en micro-gravité la courbe est plus lisse, il n'y a que la présence d'un composant électronique qui fait varier le taux de vide.

L'évolution des vitesses, illustrée en figure 33, le démontre parfaitement : en gravité normale, dans les parties ascendantes la vitesse du gaz augmente, celle du liquide diminue. Puis, lorsque le titre tend vers zéro, le vitesse du liquide s'annule.

L'évolution des pertes de charges (cf. illustration 34) comprend les pertes par frottement pariétal, interfacial, et dans les coudes. Dans le cas monophasique, les pertes de charges dues au gaz sont également prises en considération. Le gradient de pression totale est plus important en gravité normale qu'en micro gravité, cela s'explique par le fait qu'en gravité, le liquide frotte d'avantage la paroi.

Quant à la prédiction de la température des composants (cf. figures 35 et 36), ils n'atteignent pas leur température maximale. L'augmentation à la fin de l'évaporateur est due à la limite x=0,99 imposée, elle entraîne l'utilisation des modèles en écoulement monophasique pour lesquels la paroi du tube s'échauffe considérablement. De plus, à titres élevés, les modèles perdent de leurs performances pour prédire les échanges thermiques.

Finalement, les simulations prédisent un bon fonctionnement de l'évaporateur, le cahier des charges est respecté.

Figure 31: Évolution du titre massique le long de l'évaporateur, en 1-g et μ-g.

 

 

Figure 32: Évolution du taux de vide le long de l'évaporateur, en 1-g et μ-g.

 

 

Figure 33: Évolution des vitesses  le long de l'évaporateur, en 1-g et μ-g.

 

 

Figure 34: Évolution des pressions  le long de l'évaporateur, en 1-g et μ-g.

 

 

Figure 35: Température des composants, en 1-g.

 

 

Figure 36: Température des composants, en μ-g.

 

Conclusion & Perspectives

Conclusion & Perspectives

Ce travail de recherche appliquée n'est pas le premier en la matière. A notre connaissance, plusieurs recherches approfondies ont été faites en ce qui concerne la modélisation de l'écoulement et des transferts de chaleur au sein de l'évaporateur. L'étude consiste à rechercher le modèle de simulation adéquat. La méthode consiste à tester tous les modèles empiriques susceptibles de répondre à nos besoins.

La première phase du travail a consisté à faire une étude bibliographique. Elle a permis de cerner le sujet et sa problématique. A la suite de ces recherches, un modèle analytique a pu être déterminé et a servi à la programmation puis à réfuter ou valider les modèles testés.

La seconde phase du travail, qui constitue le cœur du sujet, a traiter de l'étude de différents modèles pour les pertes de charge par frottement et pour le coefficient d'échange thermique.

A la suite de la programmation Matlab, il a été mis en évidence que le modèle thermique de Kandlikar était parfaitement adapté au problème en gravité normale et micro-gravité. Quant aux pertes de charge, il s'est avéré que le régime gravitationnel influençait leur modélisation. Ainsi, le modèle de Lockhart & Martinelli est le plus pertinent en gravité normale ; tandis qu'en micro-gravité, le modèle d'Awad est préférable.

L'utilisation du modèle de A. Cioncolini, qui fait intervenir l'arrachage de gouttelettes à l'interface, ne permet pas de modéliser l'écoulement et les transferts de chaleur correctement. Il semble que l'arrachage soit négligeable au sein de l'évaporateur.

A la suite de ces résultats, les modèles choisis ont été implémentés dans un second programme Matlab afin de calculer les différents paramètres, comme les pertes de pression au sein de l'évaporateur, et la température des composants. Les résultats ont permis de vérifier que le cahier des charges était respecté dans les conditions d'utilisation spatiales.

A l’issue de ces recherches, des perspectives de travail sont envisageables. Certaines sont en cours de réalisation : un prototype de la boucle frigorifique est sur le point de délivrer les premiers résultats. Une fois obtenus, ces résultats pourront être comparés à ceux du programme Matlab.

Après quoi, si la modélisation s'avère exploitable, il sera également utile de mener des études de sensibilité. Elles permettront de l'estimation de l'influence des paramètres hydrauliques et thermiques sur le refroidissement des composants. Afin de quantifier l'influence ce ces différents paramètres, une étude paramétrique pourra être menée. Il s'agit de faire varier le débit.

Enfin, à moyen terme, il sera possible d'utiliser cet outil, dès la conception de satellites, en phase de développement. L'enjeu sera donc de prévoir la température des composants.

Nomenclature

 

Nom                       Description

A                                    Section du tube (m²)

Cp                                  Capacité calorifique (J/kg/K)

D                                    Diamètre du tube (m)

Dh                                  Diamètre hydraulique (m)

e                                    Épaisseur du tube (m)

E                                   Taux d’entraînement (-)

f                                     Coefficient de frottement (-)

G                                   Débit surfacique (kg/m²/s)

h                                    Coefficient d'échange thermique (J/kg)

hlg                                  Chaleur latente de vaporisation ((J/kg)

J                                    Vitesse superficielle (m/s)

m                                   Débit massique (kg/s)

L                                    Longueur de l'évaporateur (m)

P                                    Pression (Pa)

q                                    Flux de chaleur surfacique échangé à la paroi (W/m²)

R                                   Résistance thermique (K.m²/W)

Rg,Rl                              Taux de vide gazeux ou liquide (-)

Sp                                  Périmètre mouillé du tube (m)

 t                                   Epaisseur du film liquide (m)

T                                   Température (°C)

U                                   Vitesse moyenne (m/s)

x                                    Titre massique de vapeur (-)

z                                    Abscisse du tube (m)

 

Symbole grecs

λ                                   Conductivité thermique (W/m/K)

μ                                   Viscosité dynamique (Pa.s)

ν                                   Viscosité cinématique (m²/s)

ρ                                   Masse volumique (kg/m³)

σ                                  Tension de surface (N/m)

 

Notations

e                                   Entrée

g                                   Gaz

i                                    Interface

l                                    Liquide

m                                  Mélange

p                                   Paroi

s                                   Sortie

sat                                Saturation

 

Références

Références

[1] B. Boujida, D. Cicoria et A. Voinis, Dimensionnement d'un évaporateur pour une boucle de contrôle thermique d'un satellite de communication, Rapport Bureau d’Étude Industriel, 2012.

[2] C. Colin, Changement de phase : Ébullition et condensation convective, Cours ENSEEIHT-INPT, 2012- 2013

[3] P. Vassallo and K. Keller, Two-Phase Frictional Pressure Drop Multipliers for SUVA R-134a Flowing in a Rectangular Duct,US Government, 2004.

[4] Awad and Muzychka, Review and modelling of two-phase frictional pressure with gradient at microgravity conditions, 2010.

[5] Padilla and Revellin, Prediction and simulation of two-phase pressure drop in return bends, Thèse INSA de Lyon, 2009.

[6] N. Antonsen and J. R. Thome, Predictions of annular flow parameters for the MANBO project, Réunion Projet MANBO, ESTEC Noordwijk, 6-7 février 2013.

[7] A. Cioncolini and J. R. Thome, Entrained liquid fraction prediction in adiabatic and evaporating annular two-phase flow, Nuclear Engineering and Design , 2012.

[8] A. Cioncolini and J. R. Thome, Void fraction prediction in annular two-phase flow, Int. J. of Multiphase Flow , 2012.

[9] A. Cioncolini and J. R. Thome and C. Lombardi, Unified macro-to-microscale model to predict two- phase flow, Int. J. of Multiphase Flow, 2009.

[10] A. Cioncolini and J. R. Thome, Algebraic turbulence modeling in adiabatic and evaporating annular two- phase flow, Int. J. of Heat and Fluid Flow , 2011.