Ebullition dans les réservoirs cryogéniques de lanceurs Ariane V-A5ME

Guillaume Davy (ENSIACET) et Sébastien Renaudière de Vaux (ENSEEIHT)

 

 


Encadrants :


 

Le Bureau d'Etudes Industrielles (BEI) entre dans le cadre de la formation d'ingénieur en Mécanique des Fluides de l'ENSEEIHT. Le projet se déroule en partenariat avec un industriel, ici Air Liquide, et l'ENSEEIHT.

Il s'agit dans ce BEI de modéliser l'ébullition dans les réservoirs cryogéniques de lanceurs Ariane V-A5ME. L'étude est réalisée à l'aide du code NEPTUNE_CFD, dans lequel il faudra implémenter des modèles d'ébullition en micro-gravité.

Logo_partenaires

 

 

Une partie de cette étude (notamment les dimensions de la cellule) est confidentielle et ne sera pas présentée ici.

 

Si certaines figures ne sont pas lisibles, il suffit de cliquer dessus pour les agrandir.

Introduction

Contexte industriel

Le but de cette étude est de réaliser la simulation d’ébullition d’oxygène liquide. Ce carburant, est utilisé dans les réservoirs des lanceurs Ariane V-A5ME. L'oxygène en contact avec les parois chaudes du réservoir peut se vaporiser et entraîner une montée en pression du réservoir. Il est donc bien important de contrôler les phénomènes d'ébullition en paroi en condition de microgravité, lors de vols balistiques A partir de cas tests expérimentaux Air Liquide, nous essayerons de simuler des cas d’ébullition avec le code NEPTUNE_CFD. Cette étude entre dans le cadre  de l’expérience OLGA du CEA. Au final, les résultats OLGA permettraient de valider les simulations et de réutiliser le code pour simuler de l’ébullition dans des réservoirs cryogéniques de lanceurs spatiaux.

 

Moteur Vinci, image Snecma

Objet de l’étude

Dans un premier temps, nous réaliserons des cas tests de convection naturelle afin de valider le code sur des géométries et des cas tests connus. Après ça nous réaliserons de l’ébullition d’eau, basée sur les paramètres et conditions opératoires fournis par EDF/CEA directement intégrés dans NEPTUNE_CFD. On adaptera ensuite ces modèles à la géométrie des cas tests Air Liquide. Nous nous intéresserons notamment aux flux thermiques imposés de chauffe et à l’influence de la gravité sur la génération de l’ébullition et l’évolution de la température.

 

Banc OLGA, image Air Liquide

Présentation du code

La version de NEPTUNE_CFD utilisée est la version NEPTUNE_CFD V1.08@Tlse. Ce code résout les équations moyennées de conservation de la masse, quantité de mouvement et énergie avec un modèle à deux fluides (modélisation Euler-Euler). Il n'est donc pas possible de remonter à un diamètre de bulle, mais à un taux de présence de la phase liquide ou gazeuse. Il faut pour cela ajouter une équation de transport d'aire interfaciale.

NEPTUNE_CFD résout les équations suivantes :

Equations moyennes de conservation de la masse de la phase $k$ :

\begin{equation} \dfrac{\partial}{\partial t} (\alpha_k \rho_k) + \dfrac{\partial}{\partial x_j} (\alpha_k \rho_k U_{k,j} )= \Gamma_k \end{equation}

où :

Equations moyennes de quantité de mouvement de la phase $k$ :

\begin{equation} \dfrac{\partial}{\partial t} (\alpha_k \rho_k U_{k,i}) + \dfrac{\partial}{\partial x_j} (\alpha_k \rho_k U_{k,j}  U_{k,i})= \alpha_k \rho_k g_i + \dfrac{\partial}{\partial x_j} \left[\alpha_k (\tau_{k,ij}+\tau^t_{k,ij}) \right] + I_{kj,i} \end{equation}

où :

Equations moyennes de l'énergie de la phase $k$ :

\begin{equation} \dfrac{\partial}{\partial t} (\alpha_k \rho_k E_k) + \dfrac{\partial}{\partial x_j} (\alpha_k \rho_k U_{k,j}  E_{k})= \dfrac{\partial}{\partial x_j} \left[ \alpha_k (\tau_{k,ij} U_{k,j}-q_{k,j}-q^t_{k,j} ) \right] + \alpha_k \rho_k U_{k,j} g_j +  \Pi_k^E \end{equation}

où :

 

Validation du code en convection naturelle

Avant de lancer le calcul sur une situation réelle, il est nécessaire de valider le code sur des situations tests, soit en validant des cas théoriques simples, soit en retrouvant des résultats expérimentaux. Dans notre cas, il existe assez peu de résultats théoriques, mais de nombreux résultats expérimentaux ont été publiés. La convection naturelle entre deux plaques planes a longtemps été étudiée, dans diverses configurations. Nous étudierons deux cas particuliers de convection naturelle :

L'étude portera sur le comportement du nombre Nusselt à différents nombres de Rayleigh pour de l'air lors de la validation.

Nombre de Nusselt

Le nombre de Nusselt représente un flux de chaleur convectif adimensionné par le flux conductif correspondant au même écart de température et sera défini comme suit tout au long de notre étude :

\begin{equation} \boxed{Nu=\dfrac{1}{l^2} \int_0^l \int_{0}^l \left( u^* T^* - \dfrac{ \partial T^*}{\partial y^*} \right) dxdy}\end{equation}

 

où les grandeurs $.^*$ sont des grandeurs adimensionnées :

\begin{eqnarray} u^* &=& \dfrac{u l}{\kappa} \\ y^* &=& \dfrac{y}{l} \\ T^* &=& \dfrac{T-T_{froid}}{T_{chaud}-T_{froid}} \end{eqnarray}

Les grandeurs $l$, $\kappa$, $u$, $T_{froid}$ et $T_{chaud}$ réprésentent respectivement la hauteur de la boîte, la diffusivité thermique du fluide, la vitesse du fluides et les températures des parois.

Nombre de Rayleigh

Le nombre de Rayleigh compare les transferts thermiques par conduction aux transferts par convection. On le définit comme :

\begin{equation} \boxed{Ra=\dfrac{g \alpha l^3 }{\nu \kappa}}\end{equation}

avec $\alpha$ un coefficient de dilatation thermique et $g$ l'accélération de la gravité.

Loi de Gay-Lussac (air)

Pour un écoulement d'un gaz chauffé, les variations de densité doivent être prises en compte. Pour une transformation à pression constante, la loi des gaz parfaits donne :

\begin{equation}P=\rho r T \end{equation}

ce qui donne $\rho T= constante$. On obtient donc la loi de Gay-Lussac :

\begin{equation}\boxed{\rho_1 =\rho_0 \dfrac{T_0}{T_1}}\label{gaylussac}\end{equation}

Code source du fichier Fortran 77 (NCELET est le nombre de cellules) :


    DO IEL = 1, NCELET
        PROPHY(IEL,IROM(IPHAS)) = RO0(IPHAS) *T0(IPHAS) /W1(IEL)
    ENDDO

Modèle de Boussinesq (liquides)

L'écoulement est considéré comme incompressible mais les sources thermiques modifient la densité dans le terme de gravité des équations de Navier-Stokes uniquement.

\begin{equation} \boxed{\dfrac{\rho}{\rho_0} = 1-\alpha (T-T_0) }\label{boussinesq}\end{equation}

où $\rho_0$ est une densité de référence et $\alpha$ une constante qui dépend du fluide.

Code source du fichier Fortran 77 (NCELET est le nombre de cellules) :


    DO IEL = 1, NCELET
        PROPHY(IEL,IROM(IPHAS)) = RO0(IPHAS) *(1-
        & alpha*(W1(IEL)-T0(IPHAS)))
    ENDDO

Cavité carrée

Corrélation de Wright

Nous essayons de trouver une relation qui lie le nombre de Nusselt au nombre de Rayleigh. Dans notre cas de convection naturelle en cavité carré, la corrélation de Wright nous donne :

\begin{equation}   Nu=
\begin{cases}
    0,028154 Ra^{0,4134},& \text{si } Ra \le 5\cdot 10^4\\
    0.0673838Ra^{1/3}, & \text{si } 5\cdot 10^4 \le Ra  \le 10^6 \\
\end{cases} \label{wright_corr}
\end{equation}

Géométrie et maillage

Nous utilisons comme paramètres de calcul un maillage carré non structuré de 5 cm de côté, avec 20 cellules par centimètre.

 

Résultats qualitatifs

Nous observons le champ de température avec $Ra = 1371$ et $Ra = 1.6779\cdot 10^6$.

 

 

Champ température Champ température

 

La paroi chaude de gauche réchauffe le fluide, le dilate et le fait donc monter par gravité. Il est ensuite refroidi à la paroi de droite où plus dense, il descend. Les vitesses au milieu sont faibles, le transfert thermique se fait essentiellement par conduction.

Comparaison avec la corrélation

Quatre valeurs de Rayleigh ont été testées et comparées avec la corrélation de Wright. Les résultats de NEPTUNE semblent satisfaisants pour ce cas test.

Champ température

Plaques planes horizontales

Corrélation de Dropkin-Somerscales

On réalise une deuxième étude comparative qui donne encore une fois le nombre de Nusselt en fonction du Rayleigh. Le but de cette étude est de comparer les simulations NEPTUNE_CFD avec la corrélation expérimentale de Dropkin-Somerscales (1965) pour des plaques horizontales :

\begin{equation}Nu=0,069Ra^{1/3} Pr^{0,074} \end{equation}

Géométrie

Cette corrélation est valable pour deux plaques planes horizontales et parallèles de longueur $L$ et de hauteur $h$ à la limite $L/h \to \infty$. Dans notre étude, on impose $L/h=10$. Nous simulons en premier lieu avec de l'air. Les parois latérales sont adiabatiques. Le maillage utilisé comporte 100 mailles par côté.

 

Champ température

Résultats qualitatifs

Nous observons le champ de températures pour $Ra=1,1587\cdot 10^6$. Les rouleaux de convection caractéristiques de l'instabilité de Rayleigh-Bénard sont bien visibles.

Champ température

 

Comparaison avec la corrélation

Calcul du coefficent de compressibilité pour l'oxygène

Nous utilisons les tables thermodynamiques de l'oxygène liquide données par le National Institute of Standards and Technology, sur une plage allant de $55K$ à $154K$. Nous calculons le coefficient $\alpha$ du modèle de Boussinesq à $100K$ .

Boussinesq

L'hypothèse de Boussinesq est donc valable entre $55K$ et $125K$ et $\alpha=0.0049K^{-1}$.

Pour de l'air, l'état stationnaire est facilement atteint, grâce à la faible densité du gaz. La loi de Gay-Lussac est codée dans la subroutine usphyv.F, qui définit les propriétés physiques des phases.

Corrélation Ra-Nu

Les résultats de NEPTUNE semblent ici aussi tout à fait satisfaisants quelque soit le fluide, et sur une plage de nombre de Rayleigh relativement large. La convection naturelle étant donc validée, nous passons donc à l'étape d'ébullition.

Réalisation du cas d'ébullition

Après avoir validé le code en convection naturelle, il est maintenant nécessaire de passer à la géométrie réelle et au cas d'ébullition. Dans un premier temps, nous réalisons de l'ébullition d'eau (grace au module eau/vapeur de NEPTUNE).

Géométrie et maillage

Encore une fois nous sommes obligés de passer par des étapes de validation et de prise en main. Nous adoptons la géométrie du banc OLGA, avec une sortie modifiée. Le maillage 2D est relativement fin. Au centre, un élément chauffant sert à chauffer le fluide et à provoquer l’ébullition.

 

Maillage

Mise en place du cas sous NEPTUNE_CFD

Après avoir réalisé le maillage, il faut mettre en place le cas d'ébullition sous l'interface Edamox. Avant de commencer, l'utilisateur doit se placer en mode expert. L'interface de NEPTUNE_CFD est composée de plusieurs boutons et doit être parcourue comme suit :

interface

Modules spéciaux

NEPTUNE_CFD possède par défaut des modules spéciaux qui permettent de traiter des cas particuliers prédéfinis. Dans notre cas il faut sélectionner le module eau/vapeur qui est conçu pour traiter les écoulements bouillants en convection forcée. Néanmoins, nous verrons que le modèle reste pertinent pour de la convection naturelle.

 

modules

Propriétés du fluide et de l'écoulement

Ensuite, les propriétés propres aux deux phases eau et vapeur doivent être définies. Par défaut, avec le choix du module eau/vapeur, les valeurs sont entrées. Si ce n'est pas le cas, il faut entrer les valeurs à la main.

fluid

Attention, même si l'écoulement bouillant comporte des bulles dont le volume peut varier à cause de variation de pression et température, il ne faut en aucun cas cocher compressible, car l'écoulement n'est pas compressible.

Il faut ensuite ajouter des modèles de turbulence et de transferts entre phases :

fluid

Contrôle des entrées et sorties

Cet onglet permet le choix du nombre de pas de temps, de la fréquence des sorties pour le post-traitement etc. La modification de ces données ne change pas le calcul.

in_out

Généralités

C'est la partie la plus importante pour mettre en place un cas d'ébullition. Nous utilisons les tables CATHARES standard, pour de l'eau à 1bar.

generalities

Dans le cas d'un écoulement bouillant il faut toujours travailler en pression absolue et ne pas oublier de reporter la pression aux conditions limites de sortie.

 generalities

Pour les cas d'ébullition, il faut modéliser le flux de chaleur pariétal $\Phi_W$. Le modèle 4 flux utilisé est un modèle de Kurul et Podowski étendu (modèle de Seiler) :

\begin{equation*} \Phi_W = \Phi_C + \Phi_Q + \Phi_E + \Phi_{C2} \end{equation*}

  • $\Phi_C$ le flux de chaleur convectif entre la paroi et la phase liquide

  • $\Phi_{C2}$  le flux de chaleur convectif entre la paroi et la phase vapeur

  • $\Phi_Q$ le flux de Quenching : flux de condction instationaire

  • $\Phi_E$ le flux d'évaporation

Schémas numériques

La partie schémas numériques est laissée par défaut.

Scalaires

Lorsqu'on choisit le module eau/vapeur, deux scalaires sont automatiquement ajoutés. Il servent à transporter l'enthalpie des deux phases. Pour l'eau le scalaire 1 est défini comme nul à ($T_0$, $P_0$) et pour la vapeur le scalaire 2 est défini à enthalpie de saturation nulle à $P_0$, c'est-à-dire $H_{sat}(P_0)=0$.

scalars

Lorsque l'on impose $H_{sat}(P_0)=0$, cette condition impose la température de la phase vapeur à la température de saturation.

Conditions limites

C'est ici que les conditions limites doivent être définies, comme pour un cas classique.

boundary

Notre objectif est de chauffer une cellule en bas de notre cavité. Le liquide étant sous refroidi, nous voulons imposer des parois à une certaine température en dessous de celle d'ébullition. Imposer un flux de chaleur et des parois adiabatiques ne posent pas de problème, ni même qu'imposer une température de paroi et des parois adiabatiques. Il nous a cependant été impossible d'imposer des conditions mixtes, à savoir une température en paroi et un flux de chaleur sur une autre paroi adjacente. Pour contourner le problème, il faut définir la paroi à laquelle on veut imposer une température (sous forme de Dirichlet) comme étant une entrée, et imposer une vitesse nulle.

Contrôle des variables de sortie

C'est ici que l'utilisateur va définir les variables qu'il souhaite avoir pour son post-traitement.

Lancement du calcul

Nous pouvons à présent lancer nos calculs.

Résultats de simulation

Sans refroidissement

Dans un premier temps, nous nous sommes intéressés au changement de gravité. En prenant la géométrie précédente, à flux fixé, nous avons fait varier la gravité de 0 à 1g. Toutes les autres parois (hormis la cellule chauffante) sont adiabatiques. Nous chauffons à flux constant de 450 000 W/m² sans refroidir. L'eau est initialisé à $5K$ en-dessous de la température de saturation, soit à $367.14K$ à 1bar.

Il s’agit ici d’une phase de test. En effet il faut refroidir la paroi autour de la cellule chauffante. On voit cependant bien l’impact de la gravité sur les forces de flottabilité favorisant ainsi le détachement de petites bulles. Celles-ci quand elles existent se recondensent très rapidement dans l’eau qui est sous-refroidie. Nous pouvons aussi observer le champ de température :

De même, les effets de convection naturelle sont très présents pour les gravités fortes, alors qu'ils sont inexistants à gravité zéro.

Avec refroidissement

En faible gravité, le paquet de bulles a tendance à s'étaler. Pour éviter cela et rester dans le cadre de l'expérience OLGA, nous imposons une température de $367.14K$, soit un sous refroidissement de $5K$ sur les parois de la partie basse autour de la cellule chauffante.

Globalement le taux de vide reste similaire par rapport au cas non refroidi, mais la position des bulles est légèrement différente. Refroidir les parois atténue le fait que les bulles soient emportées par les gros tourbillons de convection. Seul en gravité 0 les bulles ont tendance à s'étaler plutôt qu'à être emportées par l'écoulement.

Comparaison des modèles de transfert de chaleur en paroi

Deux modèles sont disponibles pour modéliser le flux en paroi :

  • modèle à 3 flux de Kurul & Podowski,
  • modèle de Kurul & Podowski étendu à 4 flux.

Modèle à 3 flux

Dans ce modèle le flux de chaleur est composé de trois contributions :

\begin{equation} \Phi_W = \Phi_C + \Phi_Q + \Phi_E  \end{equation}

où :

  • $\Phi_C$ est le flux convectif pour le liquide,
  • $\Phi_Q$ est le flux de quenching (flux conductif instationnaire),
  • $\Phi_E$ est le flux d'évaporation.

Dans ce modèle, toute la chaleur sert à chauffer le liquide. Il n'y a pas de prise en compte des bulles en paroi.

Modèle à 4 flux

Ce modèle considère un quatrième flux $\Phi_{C2}$, qui est le flux convectif pour la vapeur. Ce modèle prend donc en compte les bulles à la paroi.

\begin{equation} \Phi_W = \Phi_C + \Phi_Q + \Phi_E  +\Phi_{C2} \end{equation}

Courbe d'ébullition

Nous simulons pour différents flux et pour chaque modèle, afin d'obtenir une courbe d'ébullition. La température de paroi est prise au centre de la cellule chauffante, afin de s'épargner des effets de bord. La figure de droite est un agrandissement de la zone avant flux critique.

On peut voir que le flux critique n'est pas atteint avec le modèle à 3 flux. Si on impose un flux plus élevé, il n'y a pas de convergence. Le modèle à 4 flux semble mieux se comporter. Nous pouvons aussi regarder l'évolution temporelle de $T_W$ pour un flux donné :

En régime établi, le modèle à 4 flux ne varie plus et $T_W$ est constante. Au contraire, le modèle à 3 flux a des variations périodiques. Celles-ci pourraient être dues à des détachements de bulles périodiques. Nous n'avons cependant pas d'explication pour les trois pics qui apparaissent à chaque fois en un point. Globalement le modèle à 4 flux impose une température en paroi plus chaude et la vapeur est plus surchauffée qu'avec le modèle à 3 flux.

Ebullition avec de l'oxygène liquide

Cette étape consiste en l'implémentation dans les sources du code de tables ou de corrélation pour les grandeurs physiques de l'oxygène liquide. Les tables CATHARE ne seront désormais plus utilisées (désactivées dans l'interface Edamox).

Les tranferts inter-phases d'énergie, de masse et de quantité de mouvement ne sont pas modifiés. Il s'agit de coder des lois d'évolution des grandeurs physiques dans la routine usphyv.F. C'est le tableau PROHY(NCELET,NPROP) qui doit être modifié.

Description du tableau PROHY(NCELET,NPROP)

Le tableau PROPHY contient les données physiques des deux phases. Il est structuré comme cela :

Description du tableau PROPHY
Entr​ée Signification
PROPHY(IEL,ITEMPK(IPHAS)) Température de la phase $i$
PROPHY(IEL,ICP(IPHAS)) Capacité calorifique de la phase $i$
PROPHY(IEL,IHSAT(IPHAS)) Enthalpie de saturation de la phase $i$
PROPHY(IEL,ITSAT) Température de saturation
PROPHY(IEL,IHLAT) Chaleur latente
PROPHY(IEL,ITHSIG)) Tension de surface à saturation
N/A viscosité, diffusivité de la chaleur

Ajout des propriétés dans le code

Nous avons implémenté un modèle de Boussinesq pour la phase liquide. Les $Cp$ sont supposés constants, et nous avons fixé l'enthalpie et la température à saturation.

Partie du code source du fichier Fortran 77 usphyv.F pour l'Oxygène liquide où PROPHY est une table de propriétés physiques, IEL est le numéro de la cellule :


C        Enthalpie de saturation
    hlsat = 178438.8
    hvsat = 102577.419
    Tsat = hlsat / Cp1
    PROPHY(IEL,IHSAT(1)) = hlsat
    PROPHY(IEL,IHSAT(2)) = hvsat
    PROPHY(IEL,ITSAT) = Tsat
C        Chaleur latente
    PROPHY(IEL,IHLAT) = abs(hlsat - hvsat)
C        Tension de surface
    PROPHY(IEL,ITHSIG) = 11.0d-3

Premiers résutats

Nous avons pour l’instant choisi de ne pas implémenter de loi d’évolution des grandeurs, mais de les laisser constantes pour les premières simulations. Seule la densité est variable avec un modèle de Boussinesq. Les conditions initiales sont les mêmes qu’au cas précédent, à savoir du LOx dans la cavité et de la vapeur à la sortie. Nous conservons le modèle à 4 flux car le modèle Kurul & Podowski de base échoue. En revanche, nous avons diminué le diamètre de détachement des bulles, en le passant de 1mm à 0,4mm. Ce diamètre correspond aux observations effectuées par le groupe chargé de la modélisation. Nous chauffons la cellule sans refroidir les parois adjacentes. Nous observons donc les évolutions des taux de vide pour différentes gravités. Au delà d'un certain temps néanmoins, les calculs divergent. 

 

Dans le dernier cas, en gravité zéro, les bulles ne montent plus mais coalescent pour former une grande bulle accrochée à la paroi :

De manière générale, les bulles sont peu refroidies après détachement et ne se recondensent pas.

Analyse numérique

Ces simulations nous ont donné accès à un certain nombre de paramètres : les différents flux du modèle 4-flux de NEPTUNE_CFD, la température de paroi Tw, les températures, les densités, les fractions volumiques ...

En ce qui concerne les flux, le flux total se stabilise rapidement avec le temps vers le flux imposé dans NEPTUNE. Ce qui est étonnant en revanche, c'est le faible flux de convection liquide qui représente 10% du flux total, et le faible flux d'évaporation (inférieur à 5%). La grande majorité (+80%) du flux total est contenu dans le flux de Quenching.

Évolution de la température de paroi en fonction de la gravité
Gravité (m/s²) Température de paroi Tw (K)
9,81 122.6
4,905 123.5
0,981 124.5
0 186

Comme dit précédemment, en 0g, toutes les bulles coalescent la cellule chauffe uniquement la vapeur, et la température de paroi croit fortement.

Ci-dessous, les profils verticaux de taux de vide sont représentés à différentes gravités. Plus la gravité diminue, plus la fraction de gaz en paroi est forte. En gravité 0,5g, des paquets de bulles se détachent et remontent périodiquement.

interface

Conclusion

Les résultats obtenus sont prometteurs pour la suite. En convection naturelle, NEPTUNE_CFD permet d'avoir de bon résultats qui correspondent aux différentes corrélations. Les modèles se comportent globalement bien, même s'ils sont sont conçus a priori pour des écoulements en convection forcée. Le modèle d'ébullition nucléée semble correct, même en très faible gravité. La première partie de cette étude a été faite en utilisant les fonctions CATHARE.

L'étape suivante a donc consisté en la modification de la lecture de ces tables, afin d'implémenter des tables ou corrélations pour l'oxygène liquide. La routine utilisateur à modifier est usphyv.F, qui contient les données physiques des fluides. Ainsi, il faut implémenter les comportements des densité, viscosité, diffusivité thermique, capacité calorifique et chaleur latente de chaque phase. Ce travail a déjà été réalisé pour un mélange LOx/Helium et doit être réajusté afin de concorder avec ce cas. Nous avons pour notre cas gardé des propriétés constantes pour le fluide, hormis la densité qui suit un modèle de Boussinesq. Le plus simple serait de reprendre les tables de l'oxygène liquide et de chercher une fonction polynomiale approchée, comme le modèle de Boussinesq qui est une approximation au premier ordre.

Remerciements

Nous tenons à remercier MM. Hervé Neau et Renaud Ansart pour leur aide précieuse tout au long du projet, les débuggages et les tutoriaux. M. Fabrice Mathey et la société Air Liquide pour leur confiance. Et Mme Catherine Colin pour nous avoir supervisé, son soutient et son expertise.

Références

[1] NIST chemistry webbook. Consulté : 13/02/2014.

[2] Air-Liquide : Modélisation des transferts par ébullition dans les réservoirs cryogéniques de lanceurs Ariane V-A5ME ; Marylou Garnier, Zacharie Vatimbella Consulté : 12/03/2014.

[3] Bassenghi, Federica. "Validation of the CFD code NEPTUNE for a full scale simulator for decay heat removal systems with in-pool heat exchangers." (2013).

[4] Dropkin, D., and E. Somerscales. "Heat transfer by natural convection in liquids confined by two parallel plates which are inclined at various angles with respect to the horizontal." Journal of Heat Transfer 87 (1965): 77.

[5] Incropera, Frank P., Adrienne S. Lavine, and David P. DeWitt. Fundamentals of heat and mass transfer. John Wiley & Sons Incorporated, 2011.

[6] Kurul, N., and Mm Z. Podowski. "Multidimensional effects in forced convection subcooled boiling." Proceedings of the Ninth International Heat Transfer Conference. Vol. 2. 1990.

[7] Montout, Michaël. Contribution au développement d'une Approche Prédictive Locale de la crise d'ébullition. Diss. INPT, 2009.

[8] Wright, John L. "A correlation to quantify convective heat transfer between vertical window glazings." TRANSACTIONS-AMERICAN SOCIETY OF HEATING REFRIGERATING AND AIR CONDITIONING ENGINEERS 102 (1996): 940-9

[9] Yeoh, G. H., et al. "Fundamental consideration of wall heat partition of vertical subcooled boiling flows." International Journal of Heat and Mass Transfer 51.15 (2008): 3840-3853.
 
[10] Seiler-Marie, Nathalie. "Modélisation et simulation des phénomènes d'ébullition et du transfert de chaleur dans la zone d'impact d'un jet sur une plaque chaude." (2004).