Modélisation du flux de chaleur en paroi d'un réacteur à eau pressurisée en vue de prédire la crise d'ébullition

Avant-propos

 

Le bureau d'études industriel (BEI) est un projet qui s'inscrit dans l'enseignement de L'ENSEEIHT. En effet, en fin de formation les étudiant sont mise en situation d'ingénieur pour résoudre une problématique réelle proposée par des industriels. Ainsi, ils s'investissent à plein temps sur  le sujet choisi par binôme ou trinôme pendant six semaines, avec au préalable une période d'impréniation du sujet qui se résume à sept séances de deux heures qui ont lieu entre novembre et janvier. Dans ce sens, des réunions avec l encadrant, et éventuellement le contact industriel ont lieu pour préciser le cadre de l'étude et la définition des objectifs.

L'intiltulé du bureau d'étude sur lequel nous allons travaillé est :

"Modélisation du flux de chaleur en paroi d'un réacteur à eau pressurisée en vue de prédire la crise d'ébullition"

Ce sujet a été proposé par EDF. Dès lors, ce projet sera réalisé avec la collaboration de monsieur Michaël Montout, ingénieur recherche chez EDF, et sera encadré par madame Catherine Colin, chercheuse à l'IMFT, au sein du groupe Interface et professeur à L'INPT. En ce qui concerne, les étudiants chargés de l'étude, il s'agit de mademoiselle Aminata Collé Diagne et de monsieur Yahia Sadik.

Le groupe EDF par son statut de premier électricien nucléaire mondial, se doit d'assurer la pérennité de ses ouvrages. Le parc nucléaire français est le plus fourni d'Europe et il réside des risques inhérents au bon fonctionnement de la centrale, parmi lesquels la crise d'ébullition dans le circuit primaire. En effet, ce phénomène déplorable détériore la surface des barreaux de combustible. Le but de cette étude serait donc d'améliorer les modèles locaux actuels de prédiction de ce phénomène.

Contexte

Le contexte

 

 

     Les réacteurs à eau pressurisé constituent l'essentiel du parc nucléaire français. Ainsi, l'énergie issue de la réaction sert à vaporiser de l'eau qui entraînera les turbines qui produiront l'électricité. La vaporisation de l'eau se fait au sein des générateurs de vapeur. Dans le circuit primaire, l'eau pressurisée n'est pas vaporisée dans des conditions normales de fonctionnement. Cependant, en cas de dépression brutale, ou d'augmentation ponctuelle de température une ébullition pourrait se déclencher. Ainsi, cette étude s'inscrit dans une problématique sécuritaire au niveau du fonctionnement du cycle primaire.

 

Schéma de principe du fonctionnement d'un REP

 

        Source

       En effet, une ébullition nuclée peut évoluer en une ébullition en film. Dans ces conditions, les coefficients d'échange chutent brutalement. L'énergie contenue dans les barreaux n'est plus dissipée par le fluide; la température des parois augmente donc de façon drastique. Cet échauffement ponctuel risque d'endommager lesdites parois et de compromettre la longévité du matériel. Dans ce sens, des travaux antérieures ont été réalisés dans le développement d'un modèle prédictif local de la crise d'ébullition, notamment une thèse intitulée "Contribution au développement d'une Approche Prédictive Locale de la crise d'ébullition", dont le titulaire est monsieur Montout.

 

Érosion d'un barreau par la crise d'ébullition

 

       Dès lors, le but de cette étude serait de consigner l'extension possible des corrélations de référence, et de modéliser la génération de la bulle; un nouveau modèle de prédiction de l'ébullition à l'échelle locale des bulles, inspiré du modèle de Kurul et Podowsky serait alors mis sur pied. Dans ce sens, un modèle local requière de renseigner le diamètre de détachement des bulles, leur fréquence de formation et la densité de sites activés. Ainsi, les différentes corrélations seront comparées entre elles pour déterminer une extension de leurs domaines de validité.

Choix d'une corrélation de référence

Choix d'une corrélation de référence

 

La littérature fournit de nombreuses corrélations pour prédire le flux de chaleur pariètal. Étant donné qu'il n y existe pas de banque de données expérimentales pour toutes les gammes de surchauffe, de pression, et de sous-refroidissement, la comparaison entre elles des différentes corrélations disponibles permettra d'extraire au moins une corrélation  qui servira de référence par la suite pour déterminer la fiabilité des modèles développés.

Corrélations

Corrélations

 

Corrélation de Jens et Lottes [1]

Cette corrélation a été établie en 1951 à partir de banque de données du Massachusetts Institut of technologie et de l'Université California Los Angeles. Dans celle-ci, la pression est exprimée en bar et le flux obtenu est en MW/m². Elle ne prend guère en compte, ni la température du liquide loin de la paroi, ni certaines propriétés intrinsèques au fluides.

$$\Phi_w=(\frac{(T_w-T_{sat})}{25}exp(\frac{P}{62}))^{4}$$

 

Corrélation de Thom et al [2]

A l'instar de la corrélation précédente dont se sont inspirés les auteurs de celle-ci, la corrélation de Thom et al. a été développée en 1965, à partir de données expérimentales pour des plages de flux et de débit bien définis, répertoriés dans le tableau de synthèse. La pression s'exprime dans le cas présent toujours en bar et le flux en MW/m²

$$\Phi_w=(\frac{(T_w-T_{sat})}{22,65}exp(\frac{P}{87}))^{2}$$

 

Corrélation de Chen [3]

La corrélation de Chen établi en 1966, était destinée de prime abord à la prédiction d'ébullition nucléée saturée, puisque ayant été concçue sur la base de données collectées dans ce régime. Ce n'est qu'en 1970 qu'elle a été étendu à l'ébullition sous saturée par Butterworth.

Le principe de cette corrélation repose sur la conception de flux tributaire de deux contributions: un flux de convection forcée et un flux d'ébullition nucléée. Deux paramètres S et F permettent de pondérer ces deux flux, ils sont fonctions du nombre de Reynolds et des autres caractéristiques de l'écoulement.

Le coefficient d'échange thermique par convection forcée est modélisé par la corrélation de Dittus-Boelter établi en 1930.

$$h_{fc}=0,023\frac{k_l}{D_{hyd}}{Re_l}^{0,8}{Pr_l}^{0,4}$$

Le coefficient d'échange thermique par ébullition est modélisé par la corrélation de Forster et Zuber établi en 1955 pour de l'ébullition en vase/

$$h_{nb}=0,00122\frac{{k_l}^{0,79}{c_pl}^{0,45}{\rho_l}^{0,49}}{\sigma^1,5{\mu_l}^{0,29}{h_lg}^{0,24}{\rho_g}^{0,24}}(T_w-T_{sat})^0,25\Delta P_s^0,75$$

$\Delta P_s$ désigne l'écart en pascal entre la pression à saturation à la température de paroi Tw et la pression de saturation à la température de saturation Tsat.

Les paramètres de pondération sont définis ci-après, le symbole Xtt désigne le paramètre de Martinelli.

$$S=\frac{1}{1+2,53*10^{-6}{Re_l}^{1,17}}$$

$$\frac{1}{X_{tt}}={(\frac{x}{1-x})}^{0,9}{(\frac{\rho_l}{\rho_g})}^{0,5}{(\frac{\mu_l}{\mu_g})}^{0,1}$$

F=1 si $\frac{1}{X_{tt}}$

$F=2,35*(\frac{1}{X_{tt}}+0,213)^0,736$ si $\frac{1}{X_{tt}}$

Dans le cas présent, x désigne le taux de vide; vu que l'ébullition ne survient qu'en fonctionnement accidentel et incidentiel, le taux de vide est supposé infiniment petit. Dans le cadre des simulations, il sera supposé nul.

Le flux thermique total s'exprime donc comme suit.

$$\Phi_w=\Phi_{fc}+\Phi_{nb}=h_{fc}(T_w-T_l)F+h_{nb}(T_w-T_{sat})S$$

 

Corrélation de Gungor et Winterton [4]

Tout comme la corrélation de Chen, elle est basée sur la superposition des flux thermiques induits par la convection forcée et par le nucléation, et elle a été développée en 1991.

Le coefficient d'échange par convection forcée reste celui proposé par Dittus-Boelter et celui par ébullition nucléée sera dorénavant, celui proposé par Cooper et qui se décline comme suit.

$$h_{nb}=55p^{0,12}(-log10p)^{-0,55}M^{-0,5}\Phi_w^{0,67}$$

Dans le cas présent p représente la pression réduite =$p=\frac{p}{p_{crit}}$, sachant que la pression critique de l'eau est de 221,2 bars.

Le nombre d'ébullition s'écrit:

$$B_0=\frac{\Phi_w}{h_lg G}$$

$$E=1+24000B_0+1,37(\frac{1}{X_{tt}}^0,86)$$

$$ S=\frac{1}{1+1,15*10^{-6}E²{Re_l}^1,17}$$

Le flux total s'écrit donc:

$$\Phi_w=\Phi_{fc}+\Phi_{nb}=h_{fc}(T_w-T_l)+h_{nb}(T_w-T_{sat})S$$

 

Corrélation de Liu et Winterton [5]

La corrélation de Liu et Winterton utilise les mêmes corrélations pour les coefficients d'échange que la corrélation de Gungor et Winterton. Elle a été validée sur la meme banque de données. En effet, seule la pondération des flux et la formule du flux total changent et se déclinent comme suit.

$$F=[1+xPr_l(\frac{\rho_l}{\rho_g}-1)]^{0,35}$$

$$S=\frac{1}{1+0,055F^0,1{Re_l}^0,16}$$

$$\Phi_w=\sqrt{(Fh_{fc}(T_w-T_l))²+(Sh_{nb}(T_w-T_{sat}))²}$$

 

 

Simulations et interprétation

Simulations et interprétation

 

La corrélation de Gungor et Winterton n’est pas valable dans le cadre d’une ébullition sous-refroidie, nous avons néanmoins choisi de la représenter pour estimer l’éventualité de l’extension de son domaine de validité.

Ainsi, la corrélation de Liu et Winterton serait la plus fiable dans toutes les plages de pressions balayées par l’étude, en donnant des résultats pertinents et proche des corrélations de références validées aux points choisis. La corrélation de référence, cependant, à la pression atmosphérique reste la corrélation de Chen.

Modèles mécanistes

Modèles mécanistes

Les codes de calcul, notamment le code NEPTUNE CFD, utilisent en général de modèles de prédiction de type mécaniste. En effet, ils modélisent la chaleur nécessaire pour réchauffer le liquide et la chaleur nécessaire pour la génération de bulle. Ainsi, au cours de cette étude, les modèles de Kurul et Powsky et Yeoh et Al seront particulièrement analysés.

Modèle de Kurul et Podowski

Modèle de Kurul et Podowski

Selon le modèle de Kurul et Podowsky, le flux pariétal est la somme de trois contributions :

•Un flux de chaleur monophasique par convection forcée $\Phi_{fc}$
•Un flux de chaleur par conduction instationnaire $\Phi_{tc}$
•Un flux net de chaleur par évaporation$\Phi_e$

Deux cas de figure peuvent se présenter lors de l’utilisation de ce modèle :

•La température de paroi $T_w$ est connue et le flux pariétal inconnu$\Phi_w$
•Le flux pariétal $\Phi_w$ est connu et la température de paroi $T_w$ inconnue

Ce modèle nécessite aussi la modélisation de quatre paramètres que sont :

  • Le coefficient d’échange monophasique par convection forcée $h_{fc}$
  • La densité de sites actifs de nucléation $N_a$
  • Le diamètre de décollage des bulles $D_d$
  • La fréquence de nucléation f

 

Modélisation du flux de convection forcée

Le flux de chaleur monophasique par convection forcée est calculé à partir de la modélisation des paramètres cités précédemment. Il s’écrit donc :

$$\Phi_{fc}=A_{fc}h_{fc}(T_w-T_l)$$

Où $A_fc=1-A_tc=1- \frac{\pi N_a D_l²}{4}$ la fraction d’aire de la paroi non-influencée par les bulles.

Modélisation du flux de conduction instationnaire

Le flux de conduction instationnaire est définie par :

$$\Phi_{tc}=A_{tc}t_w f \frac{2k_l(T_w-T_l)}{\sqrt{\pi \eta_l t_w}}$$

Où f représente la fréquence moyenne de décollage des bulles, le temps d’attente entre le décollage d’une bulle et l’apparition d’un nouvel embryon, et la durée de croissance de la bulle.

Ainsi,$f=\frac {1}{t_w+t_l}$

Cependant le modèle de Kurul et Podowski suppose que le temps de croissance est négligeable devant le temps d’attente.

Modélisation du flux net d'évaporation

Le flux net d’évaporation s’écrit :

$$\Phi_e=f\nu V_b\rho_g h_{lg} N_a$$

Où $V_b$ désigne le volume des bulles.

Modélisation des paramètres principaux

  • Le coefficient d’échange monophasique par convection forcée

Un modèle de type Dittus-Boelter [7] est utilisé pour fermer le modèle :

$$h_{fc}=0,023\frac{k_l}{D_{hyd}}Re_l^{0,8} Pr_l^{0,4}$$

  • La densité de sites actifs de nucléation

La modélisation de la densité de site de nucléation est prédite par la corrélation de
Lemmert et Chwala [8]:

$$N_a=[210(T_w-T_{sat})]^{1,8}$$

  • Le diamètre de décollage des bulles

L’expression du diamètre de décollage utilisée par Kurul et Podowski est :

$$D_l=10^{-4}(T_l -T_{sat})+0,0014$$

  • La fréquence de nucléation

La théorie de Kurul et Podowski utilise l’expression de Ceumern et Lindenstjerna [9] pour la fréquence de nucléation!

$$f=\sqrt{\frac{4}{3}g\frac{\rho_l-\rho_g}{\rho_l D_l}}$$

Simulation et interprétation

Ainsi, selon le graphique ci-après, le modèle de Kurul et Podowski est adapté aux écoulements à forte vitesse et faiblement surchauffés. En effet, pour des surchauffes importantes le modèle diverge. Par ailleurs, la modélisation du flux convectif est problématique ; car pour peu que la densité de site soit surévaluée la fraction d’aire sur la paroi influencée par les bulles devient négative.

Modèle de Yeoh et al

Modèle de Yeoh et al

Par rapport au modèle précédent, le modèle de  Yeoh et Al introduit la notion de glissement des bulles à la paroi. Ainsi, le flux de conduction instationnaire  est subdivisé en deux flux liés au détachement /décollage des bulles et au glissement de celle-ci le long de la paroi. La coalescence des bulles lorsqu’elles rencontrent le liquide froid n’est pas pris en compte, ce qui  est constitue une des limites du modèle.

Modélisation du flux par convection forcée

$$\Phi_{fc}=h_{fc}(T_w-T_l)$$

Modélisation du flux de conduction instationnaire lors  du détachement ou du décollement

$$ \Phi_{tc}=2 \sqrt{ \frac{k_l \rho_l c_{pl}}{\pi t_w}}(T_w-T_l) R_f N_a (K_{inf} \frac{ \pi D_d²}{4})t_wf+2 \sqrt{ \frac{k_l \rho_l c_{pl}}{\pi t_w}}(T_w-T_l)R_f N_a( \frac{ \pi D_d²}{4})(1-t_w f)$$

Où, $K_inf$ est pris égal à 1,8 selon les travaux de Judd et Hwang [11], $R_f$ est le facteur de réduction, il est défini par $R_f=\frac{1}{l_{sl}\sqrt{N_a}}$ et  représente la longueur de glissement de la bulle.

 

Modélisation du flux de conduction instationnaire lors du glissement

$$\Phi_{tc,sl}=2 \sqrt{ \frac{k_l \rho_l c_{pl}}{\pi t_w}}(T_w-T_l)R_f N_a l_{sl} K_{inf}D t_{wf}+2 \sqrt{ \frac{k_l \rho_l c_{pl}}{\pi t_w}}(T_w-T_l)R_f N_a f t_{sl}( \frac{\pi D_d²}{4})(1-t_w f)$$

Où $t_{sl}$ désigne le temps de glissement de la bulle

Modélisation du flux net d’évaporation

$$\Phi_e=\frac{pi}{6}\rho_g h_{lg}D_l^3R_f N_a f$$

Modélisation des paramètres principaux

Par rapport au modèle précédent, deux nouveaux éléments sont à modéliser, il s’agit de la longueur de glissement et du temps de glissement. Ainsi, l’ensemble des paramètres à modéliser sont :

Diamètre de détachement  $D_d$et diamètre de décollage $D_l$

Densité de sites de nucléation $N_a$

Coefficient d’échange monophasique par convection forcée $h_{fc}$

Le temps d’attente $t_w$, le temps de croissance $t_d$

Le temps de glissement $t_{sl}$

La longueur de glissement $l_{sl}$

Etude paramétrique

Etude paramétrique des temps caractéristiques

Etude paramétrique des temps caractéristiques

Etude paramétrique du temps d'attente

Les modèles de Basu et Al et Yeoh et Al donnent des modélisations du temps d’attente différentes. Dés lors, il serait intéressant de les comparer dans des conditions identiques pour estimer les écarts potentiels entre les deux modèles. De plus, le modèle de Yeoh et al dépend des angles de mouillage, il serait donc intéressant de voir sa sensibilité à ceux-ci.

Le temps d’attente selon la corrélation de Yeoh et Al s’écrit :

$$t_w=\frac{1}{ \pi \eta} (\frac{(T_w-T_l) C_1 r_c}{(T_w-T_{sat})-\frac{2 \sigma T_{sat}}{C_2 \rho_g h_{lg} r_c}})²$$

avec

$$r_c=(\frac{1}{C_1 C_2})^{\frac{1}{2}}(\frac{2\sigma T_{sat}k_l}{\rho_g h_lg \Phi_w})^{0,5}$$

$$C_1=\frac{(1+cos \theta)}{sin \theta}$$

$$C_2=\frac{1}{sin\theta}$$

Les simulations ont été effectuées à pression atmosphérique car les corrélations de Basu et al ne sont valides qu’entre 1 et 3,2 bars.

Les écarts entre les temps d’attente proposés par les deux modèles sont très conséquents. En effet, à faible surchauffe, le rapport des temps est de l’ordre de 10². Des études, notamment la thèse [12] de monsieur Montout avec les données de Maity[13], montrent qu’aucune de ces modélisations ne prédit fiablement le temps d’attente réel. Le modèle de Basu surestime fortement le temps réel et le modèle de Yeoh sous-estime celui-ci. Ensemble, ils permettent d’encadrer la fourchette d’ordre de grandeur.

Par ailleurs, il convient de remarquer que selon le modèle de Yeoh et al, une erreur de
30° entraine une erreur de l’ordre de 10% à forte surchauffe allant jusqu’à 30% à faible
surchauffe.

Comparaison des temps d'attente et de décollage selon Basu[14]

Selon le modèle de Kurul et Podowski, le temps de croissance est négligeable devant le temps d’attente. Dés lors, il parait intéressant de comparer ces deux temps pour estimer les limites de validité de cette hypothèse. Etant donné que nous ne disposons pas de données expérimentales, nous nous proposons de faire cette étude avec le modèle de Basu.

Ainsi, l’hypothèse est valide à faible surchauffe. En effet, à partir 10°C à pression atmosphérique, l’hypothèse est erronée car le temps de croissance devient supérieur au temps d’attente.

Densité des sites de nucléation

Densité de site de nucléation

La densité de site de nucléation permet de quantifier le nombre de sites actifs par unités de surface sur la paroi, elle dépend de l’état de la paroi et donc du temps. On dispose de très peu d’informations concernant celle-ci. Néanmoins, quelques corrélations permettent de la modéliser.

La corrélation de Lemmert et Chwala [8]

Dans la corrélation de Lemmert et Chawla la densité de site de nucléation actifs dépend uniquement de la surchauffe. Ainsi,

$$N_a=[210(T_w-T_{sat})]^{1,8}$$

La corrélation de Mikic et Rosenhow [15]

Elle donne le nombre de cavités de rayon supérieur à r 0 activée. Le paramètre m est égal à 4 ou 5.

$$N_a=(\frac{r_0}{r})^m$$

$$N_a=(\frac{r_0 (T_l - T_{sat}(P_l))h_{lg}}{2 \sigma (T_0) T_sat(P_0) v_v})^m$$

La corrélation de Basu et al. [14]

Lors du développement de leur modèle, Basu et al ont fourni cette corrélation empirique :

$N_a=0,34.10^4[1-cos(\Phi)] \Delta T_{sat}²$ si $\Delta T_{sat}<15$

$N_a=0,34.10^4[1-cos(\Phi)] \Delta T_{sat}^{5,3}$ si $\Delta T_{sat}\geq 15$

Simulation et interprétation

Les paramètres de la simulation pour la corrélation de Rosenhow est r0=250µm et m=4 pour avoir le même ordre de grandeur que les autres corrélations. Ainsi, on remarque que selon l’angle φ choisi, la corrélation de Basu varie trés fortement. Elle tend à rejoindre la corrélation de Lemmert et Chwala pour des angles assez importants. On choisit donc la corrélation de Lemmert et Chwala comme corrélation de référence.

Détachement et décollage de bulles

Détachement et décollement des bulles

Inventaire des forces

L’écoulement est de type ascendant dans un tube vertical. L’axe x est orthogonal à la paroi et l’axe y est tangent à la paroi et est orienté vers le haut (Cf. configuration de la bulle). Vu que l’on s’intéresse au détachement et au détachement de bulle, on s’intéressera particulièrement à la projection du produit fondamental de la dynamique suivant l’axe x. L’inventaire des forces s’exerçant sur la bulle est le suivant :

  • La force de flottabilité :  $F_b=V_b(\rho_l-rho_g)g$
  • La force de tension de surface : $F_s=2\int{\pi}^{0} \sigma r_c t(s) ds$

Elle est simplifiée grâce aux travaux de Klausner et al.  Ainsi, la composante intéressante de la force capillaire est approximée par :  $F_{sx}≈2 r_c \sigma \frac {pi}{\alpha _\beta}(cos \beta - cos  \alpha)$

$r_c$: représente le rayon du pied de la bulle

t le vecteur unitaire tangent à la surface de la bulle et perpendiculaire à la ligne de contact triple

s l’angle polaire le long du pied de la bulle

α et β les angles de contact (cf. configuration de la bulle)

  • La force de pression de contact : $F_{cp}=(P_g -P_l)\pi r_c²=\frac{2 \sigma}{r_r}\pi r_c²e_x$

$P_g -P_l$ la différence de pression au pied de la bulle

$r_r$ le rayon de courbure au pied de la bulle entre 5R et 2R (modèles)

  • La force de portance : $F_l=\rho_l C_l V_b (U_l-U_b)\land rot U_l$

Le coefficient de portance est modélisé grâce à la corrélation de Legendre et Magnaudet [16] qui s’écrit : $C_l=\sqrt{(C_l^{low})^2+(C_l^{high})^2}$

Avec $C_l^{low}=\frac{6}{\pi}\frac{2,255}{sqrt{Re_b S_r}(1+0,2 \frac{Re_b}{S_r})^{\frac{3}{2}}}$  , $C_l^{high}=\frac {1}{2}\frac{1+16/Re_b}{1+29/Re_b}$ et $S_r=\frac{2 \omega R}{|U_l-U_b|}$

  • La force de trainée : $F_d=\frac{1}{2}\rho_l \pi C_d R²||U_l-U_b||(U_l-U_b)$ 

Le coefficient de portance utilisé est celui développé par Legendre et al [16]: 

$$C_d=\frac{48}{Re_b}(1+g(s)-(1+s^{-3})\frac{2,211}{Re_b^{0,5}}+O(Re_b^{-5/6})$$

  $s=\frac{distance entre le centre de deux bulles}{rayon des bulles}$s=2 dans le cas de deux bulles en contact

$g(s)=s^{-3}+\frac{3}{4}s^{-6}+O(s^{-7})$:la correction de Kok

  • La force d’inertie : $F_I=- \rho_l \frac {d}{dt}(V_b(t)(U_b-U_l))-\rho_l C_{am1xy}\frac{d}{dt}(V_b(t)(U_b-U_l))- \rho_l V_b(t)(C_{am2} \frac{(dR/dt)²}{R}+C_{am3} (\frac{d²R}{dt²}))e_y$

D’après les travaux de Legendre et al [16] : $C_{am1x}=0,636$, $C_{am1y}=0,784$ ,$ C_{am2}=-\frac{1}{2}$, $C_{am3}=-\frac{1}{2}$

Diamétre de détachement

 Diamètre de détachement

  Lorsqu’ un germe de vapeur est piégé sur un site de nucléation, une bulle de vapeur d'eau est crée et grandis sur la paroi sous l'effet de la surchauffe pariétale. Ces bulles grandissent jusqu'à atteindre un diamètre maximal où elles se détachent de la paroi. On parle alors d'un diamètre de détachement. La connaissance de ce diamètre permet de connaître les temps caractéristiques des différents phénomènes d'échange de chaleur entre la paroi et le liquide.

                                                

                                                          Configuration de bulle étudiée

  Des expériences ont été réalisées par plusieurs auteurs et des corrélations ont été développées afin d'établir une relation entre le diamètre de détachement et les propriétés physiques de l'écoulement ambiant.
Dans cette partie, on présentera la corrélation de Basu et al. ainsi que la corrélation analytique établie à partir du bilan de forces agissant sur la bulle. On comparera ensuite ce résultat analytique la base expérimentale de Maity (2000) en étudiant l'influence des angles α et β.

  Corrélation de Basu et al.

  La corrélation principale utilisée pour sa modélisation est celle développée par Basu que l’on rappelle ici :

            

            

Avec,
lc : longueur capillaire
Jasup : Nombre de Jacob basé sur la sous saturation
Jasat : Nombre de Jacob basé sur la surchauffe pariétale
Ø : l’angle de contact statique pour lequel aucune fermeture n'est proposée.
Cette corrélation a été développée dans le domaine suivant :

 

                                                    

 

Relation obtenu à partir du bilan des forces

  Une autre relation a été obtenue par Montout (2009) à partir d'un bilan de forces agissant sur la bulle. L'expression trouvée par Montout est donnée par la relation suivante :

                        

  La fonction h(α,β) à la forme suivante :

 

                               

Comparaison avec les données expérimentales de Maity (2000)

 

  La relation précédente sera comparée avec les données expérimentales eau/vapeur développée en pression atmosphérique par Maity. Il s'agit des seules données disponibles pour le moment qui donnent les valeurs des diamètres de détachement.

 

                              

 

 

Calcul du paramètre K

  Le paramètre K utilisé dans la formule précédente intervient dans la loi de croissance de la bulle. La loi choisie est la loi de croissance de bulle contrôlée par les effets diffusifs proposée par Zuber rappelée ici :
  Les résultats expérimentales de la croissance de la bulle de Maity trouvés pour une vitesse du liquide de 0,25m/s et ΔTsub=0,2°C et ΔTsat=5,3°C sont présentée sur la figure suivante.

 

                                

Évolution du diamètre de la bulle en fonction du temps

  Le nombre de Jacob est calculé à partir de ΔTsat dans le tableau précédent, on trouve alors k=0,83.

Choix des angles α et β :

  Au moment de détachement, il est difficile de connaître les valeurs exactes des angles α
et β. Donc une étude paramétrique sera réalisée afin de voir l'influence de ces angles sur
les résultats obtenus.

Résultats :

  Le tableau suivant montre les valeurs obtenues pour le diamètre de détachement pour
plusieurs valeurs d'angles α et β :

 

   

Tableau récapitulatif des résultats du diamètre de détachement

  On remarque que le choix des angles de contact entre la bulle et la paroi a une influence énorme sur les valeurs du diamètre de détachement obtenues.
  Le bon accord avec le résultat de Maity est obtenu avec les valeurs de α=80° β=40° où l'erreur est inférieure à 10 %.

 

 

 

 

 

 

 

 

 

 

 

 

 

Diamètre de décollage

 Diamètre de décollage

 Après que la bulle ait atteint son diamètre de détachement, deux scénarios sont possibles :

  -  Décollage immédiat de la paroi.
  - Glissement sur la paroi et croissance de la bulle jusqu'à atteindre un autre diamètre plus grand où le décollage à lieu.
 Dans les deux cas on parle d'un diamètre de décollage.

 

                                                             

Glissement et décollage de la bulle après détachement de la paroi

 

Plusieurs corrélations sont proposées pour le diamètre de décollage. Nous allons citer les différentes corrélations utilisées et ensuite nous les comparerons avec la base expérimentale de Situ et al.

 

Corrélation de Basu et al.

  Elle a la même forme que celle du diamètre de détachement. Elle a été aussi développée à partir des données expérimentales eau/vapeur de Maity (2000)

  Les paramètres qui intervient dans cette expression ont déjà été définis dans la partie (diamètre de détachement )
 

Corrélation de Ünal et Boree et al.

  C'est la version finale des formules de Kurul et Podowski et celle de Ünal, elle s'écrit sous
la forme suivante :

                                              

 

où a et b et φ sont des paramètres définis par :

                      

 

 : une constante dépendant des propriétés physiques de la paroi et du fluide :

  où :
  kw : représente la conductivité thermique de la paroi, ρw la masse volumique de la paroi et Cw la chaleur spécifique de la paroi.

  Le domaine d’établissement de cette corrélation dimensionnelle est proposé dans le tableau suivant. Il concerne exclusivement les écoulements eau/vapeur.

 

                   

 Domaine de validité de la corrélation de Borée et al.

 

  Borée et al. ont étendu le domaine d'utilisation de cette corrélation pour comprendre les écoulement faiblement sous-saturé. Pour cela ils ont modifié l'expression de la constante b.

 

 avec  Stlim =0,0065

Relation obtenu à partir du bilan des forces

Vitesse de bulle

A partir d'un bilan de forces agissant sur la bulle pendant son glissement sur la paroi, Montout aboutit à une autre expression du diamètre de décollage. La projection de ce bilan de forces sur les deux axes est donnée par :

$$x \rightarrow -\frac {1}{8} \rho_l \nu_l²( \frac{5}{3} \frac{K^4 J_a^4}{Pr_l^2}+52,9 \frac {K² J_a²}{Pr_l})+ \frac{1}{2}C_l R²V_r =0$$

$$y \rightarrow -\frac {19}{24} \rho_l R² \frac{dV_r}{dt}- \frac{1}{4} \mu_l V_r (52,9+ \frac {19 K^2 J_a^2}{4 Pr_l})+ \frac{4}{3}(\rho_l -\rho_g) R² g =0$$

Ainsi, grâce à ces travaux, lorsque le paramètre K est supposé constante dans la loi de croissance de Zuber qui s’écrit : $R(t)=KJ_a sqrt(\eta t)$, il est possible d’obtenir une expression analytique de la vitesse relative des bulles en connaissant le diamètre de détachement qui fera office de condition initiale.

Ainsi,

$$V_r(t)=\frac{\epsilon_2}{1-\epsilon_1}t+\epsilon_3 t^{\epsilon_1}$$

Avec,

$$\epsilon_1=-(\frac{6}{19} \frac{52,9 \nu_l}{K²J_a²\eta_l}+\frac{3}{2})$$

$$\epsilon_2=\frac{32}{19}g$$

$$t_0=\frac {D_d²}{4K² J_a² \eta_l}$$

$$\epsilon_3=-(U_l\frac{D_d}{2}+\frac{\epsilon_2}{1-\epsilon_1}t_0)/t_0^{\epsilon_1}$$

Où $U_l$est la vitesse du liquide au centre de la bulle

Des simulations ont été réalisées grâce aux données de la thèse de Maity pour les diamètres de détachement, avec l’hypothèse d’un paramètre K constant dans la loi de Zuber. En lissant l’évolution des diamètres en fonction du temps, le paramètre K est de l’ordre de 0,83 en moyenne, ce qui semble assez cohérent avec les données de l’expérience numéro 2. Aux temps longs cependant, il convient de noter que le modèle diverge par rapport aux expériences. La vitesse débitante est de 0,25m/s.

  A partir d'un bilan de forces agissant sur la bulle pendant son glissement sur la paroi, Montout aboutit à l'expression suivante du diamètre de décollage :

 

                                  

 

Avec
Cl : le coefficient de portance
Vr : la vitesse relative.

Comparaison des résultats avec les données expérimentales de Maity(2000)

 On compare les valeurs du diamètre de décollage obtenue avec chaque corrélation avec les valeurs expérimentales de Maity (cf. Tableau de Maity). Le tableau suivant montre les résultats obtenus :

 

 Pour le cas i : le modèle local basé sur le bilan des forces réduit bien l'écart entre les corrélations et l'expérience, pourtant cet écart reste encore important (17,5%).

 Pour le cas ii : le modèle de Basu et al présente l'écart le moins important (9%).

 Comme on peut le remarquer, aucune des corrélations marche bien dans les deux cas donc il n'est pas évident de conclure sur ces résultats.

 

Comparaison des résultats avec les données expérimentales de Situ

  On comparera les corrélations précédentes avec la base de données de Situ et al établie à pression atmosphérique.

 

                                          

Données expérimentales de Situ et al.

 

 Le tableau suivant montre les résultats obtenus :

 

 

Comparaison des corrélations avec les données expérimentales de Situ et al.

 

 Les corrélations de Basu et al. et de Borée et al. sur-estiment beaucoup les diamètres de décollage trouvés par Situ et al. On retrouve la même remarque avec le modèle basé sur le bilan des forces pour les deux derniers cas.

 Pourtant, les valeurs trouvées à partir du bilan des forces dans les 3 premiers cas sont plus proches des données expérimentales de Situ et al. Les écarts trouvés ne dépassent pas 35 %. Ce modèle peut alors constituer une bonne amélioration des modèles précédents.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Conclusions

Conclusions

Ainsi, la modélisation des flux de chaleur en paroi de réacteur pour prédire la crise d’ébullition soulève de nombreuses problématiques, et par le biais de ce travail, nous avons pu examiner certains aspects de ce sujet. La corrélation, la plus fiable pour prédire le flux de chaleur à hautes pressions est celle de Liu et Winterton. Par ailleurs, le modèle de Kurul et Podowski est fiable à faible surchauffe, mais très peu précis à forte surchauffe. Le modèle de Yeoh et al couplé au modèle de Basu et al donne des résultats prometteurs pour la prédiction du flux de chaleur. En outre, l’expression de la vitesse relative donnée par le bilan des forces donne des résultats pertinents mais il semblerait que l’évolution du rayon de la bulle suive la loi de Zuber avec un paramètre K variable en fonction du temps, lorsqu’on étudie les vitesses à K variable. De plus, les hypothèses faites sur l’ordre de grandeur des temps caractéristiques ne sont pas toujours justifiées. Les simulations ont aussi montré que les écarts entre les prédictions de deux corrélations pour le temps d’attente ne sont pas non négligeables. En ce qui concerne les diamètres de décollage, les principales corrélations utilisées les surestiment. Ainsi, les perspectives de ce travail s’orienteraient vers une meilleure estimation  des temps de décollage et de détachement et des diamètres de décollage et de détachement puisque les modèles disponibles sont très peu précis et fiables en ce qui concerne ces paramètres.

Références

Reférences

[1] JENS, W.H., AND LOTTES, P. A. Analysis of heat transfer burnout, pressure drop, and density data for high-pressure water. Tech. rep., Argonne National Laboratories, 1951

[2]Thom, J. R. S., Walker, W. M., Fallon, T. A., and Reising, G. F. S. Boiling in subcooled water during flow up heated tubes or annuli. In Symposium on Boiling Heat Transfer in Steam Generating Units and Heat Exchangers (Institute of Mechanical Engineers,London, 1965)

[3]Chen, J. C. A correlation for boiling heat transfer for satured fluids in convective flows. Industrial and Engineering Chemistry Process Design and Development 5 (1966), 322–329

[4] Gungor, K. E., and Winterton, R. H. S. A general correlation for flow boiling in tubes and annuli. International Journal of Heat and Mass Transfer 29, 3 (1986), 351–358

[5] Liu, Z., and Winterton, R. H. S. A general correlation for saturated and subcooledflow boiling in tubes and annuli, based on a nucleate pool boiling equation. International Journal of Heat and Mass Transfer 34, 11 (1991), 2759–2766.

[6] Kurul, N., and Podowski, M. Multidimensional effects in forced convection subcooled boiling. In Proceedings of the 9th Heat Transfer Conference (1990), pp. 21–26.

[7] Dittus, F. W., and Boelter, L. M. K. Heat transfer in automobile radiators of the tubular type. University of California Publications on Engineering, Berkeley, California 2,13 (1930), 443–461.

[8] Lemmert, M., and Chwala, J. M. Influence of flow velocity on surface boiling heat transfer coefficient. Heat Transfer in Boiling (1977), 237–247

[9] Ceumern-Lindenstjerna, W. C. Bubble departure diameter and release frequencies during nucleate pool boiling of water and aqueous nacl solutions. Heat Transfer in boiling (1977).

[10] Yeoh, G. H., and Tu, J. Y. A unified model considering force balances for departing vapour bubbles and population balance in subcooled boiling flow. Nuclear Engineering and Design 235 (2005), 1251–1265.

[11] Judd, R. L., and Hwang, K. S. A comprehensive model for nucleate pool boiling heat transfer including microlayer evaporation. Journal of Heat Transfer (1976), 623–629

[12] Montout, Contribution au developpement d’une Approche Prédictive Locale de la crise d’ébullition, Doctorat de l’Université de Toulouse, 2009

[13] Maity, S. Effect of velocity and gravity on bubble dynamics. Master’s thesis, University
of California, Los Angeles, 2000

[14] Basu, N.,Warrier, G. R., and Dhir, V. K. Wall heat flux partitioning during subcooled flow boiling :
Part 1 - model development. Journal of Heat Transfer 127 (2005), 131–140.
Basu, N.,Warrier, G. R., and Dhir, V. K. Wall heat flux partitioning during subcooled flow boiling
Part 2 - model validation. Journal of Heat Transfer 127 (2005), 141–148

[15] Rohsenow, W. M. Heat transfer with evaporation. University of Michigan Press, 1953

[16] Legendre, D., Colin, C., and Coquard, T. Lift, drag and added mass of a hemispherical bubble sliding and growing on a wall in a viscous linear shear flow. Philosophical Transactions of the Royal Society A doi :10.1098/rsta.2008.0009 (2008).

Legendre, D., and Magnaudet, J. The lift force on a spherical bubble in a viscous linear shear flow.Journal of Fluid Mechanics 368 (1998), 81–126.

Legendre, D., Magnaudet, J., and Mougin, G. Hydrodynamic interactions between two spherical rising side by side in a viscous liquid. Journal of Fluid Mechanics 497 (2003), 133–166.

[17] Zuber, N. Hydrodynamic aspects of boiling heat transfer. PhD thesis, University of California, Los Angeles, 1959.

Zuber, N. The dynamics of vapor bubbles in nonuniform temperature fields. International Journal of Heat and Mass Transfer 2 (1961), 83–98.

[18] Situ, R., Mi, Y., Ishii, M., and Mori, M. Photographic study of bubble behaviors inforced convection subcooled boiling. International Journal of Heat and Mass Transfer 47 (2004), 3659–3667

[19] Unal, H. C. Bubble-departure diameter, bubble-growth time and bubble-growth rate during the subcooled nucleate flow boiling of water up to 177 bar. Tech. rep., Central Technical Institute TNO, P.O. Box 342, Apeldoorn, The Netherlands, 1973. Ref. No. 73-02829.

Unal, H. C. Maximum bubble diameter, maximum bubble-growth rate during the subcooled nucleate flow boiling of water up to 17.7 mn/m2. International Journal of Heat and Mass Transfer 19
(1976), 643–649.

[20] Boree, J., Charnay, G., Fabre, J., Legendre, D., and Magnaudet, J. Ecoulements diphasiques eau-vapeur avec changement de phase. Tech. rep., Rapport intermediaire IMFTInterface, 1992.