The CFD code chosen to carry on the study was NEPTUNE_CFD, -version NEPTUNE_CFD V1.08@Tls- a highly precise code, in particular for heat transfer, based on state-of-the-art physical modelling and a fully unstructured finite volume solver.

 It was developed between EDF (France’s electricity), CEA (French Atomic Energy Commission), and AREVA-NP, in particular for nuclear reactor simulation tools, but the behaviour of different fluids made of several physical phases or components can be modelled using the general Eulerian multi-field balance equations. These several equations are obtained from the fundamentals conservation laws of physics, restricted to Newtonian mechanics: mass conservation, momentum conservation, and energy conservation, which are solved using the Reynolds Averaged Navier Stokes (RANS) approach. These relations are valid inside each fluid phase, whereas additional jump conditions are necessary for the interfaces between the phases.

The models implemented in this CFD tool have been validated for liquid forced convection and nucleate boiling. Subsequently, the Boussinesq model, presented in the previous section, was coded to be able to take into account the natural convection regime.

Mass balance equation

The multi-field mass balance equation for the phase k is written as follows:

$\frac{\partial}{\partial t} (\alpha_k \rho_k) + \frac{\partial}{\partial x_i} (\alpha_k \rho_k U_{k,i}) = \Gamma_k$   (6)

where $\alpha_k, \rho_k, U_{k,i}$ represent the volumetric fraction, the density and the mean velocity of phase k, respectively and $\Gamma_k$ is the interface mass transfer rate on phase k, sum of all the other phases contribution.

$\Gamma_k = \sum_{p≠k} {{{\Gamma} ^ {c}}_{p→k} + {\Gamma^ {nuc}}_{w→k}}$   (7)

With ${\Gamma^c}_{p→k}$ the interface mass transfer from phase p to phase k (bulk transfer) and ${\Gamma^{nuc}}_{w→k}$ the mass transfer contribution to phase k induced by wall nucleate boiling.

In the model used for the study, these terms are calculated near the heated wall for the two phases, liquid and vapor (k = 1, 2), thus

${\Gamma^{nuc}}_{w→1}+{\Gamma^{nuc}}_{w→2}=0$   (8)

where $\Gamma^{nuc}_{w→2} \ge 0$ since vapour is produced by nucleation.

Finally, conservation relations for mass and volume lead to

$\sum \alpha_k=1$   (9)

$\sum \Gamma_k=0$, since ${\Gamma^c}_{p→k}+{\Gamma^c}_{k→p}=0$   (10)

Momentum equation

The general multi-field momentum balance equation for the phase k is expressed as

$\frac{\partial}{\partial t} (\alpha_k \rho_k U_{k,i}) + \frac{\partial}{\partial x_j} (\alpha_k \rho_k U_{k,i} U_{k,j}) = \frac{\partial}{\partial x_j} (\alpha_k \tau_{k,ij}+\sum^{Re}_{k,ij})-\alpha_k\frac{\partial P}{\partial x_i}+\alpha_k \rho_k g_{i} +I_{p→k,i}+\alpha_k S_k$   (11)

Where $g_i$ represents acceleration due to gravity and P the mean pressure. The viscous stress tensor is defined by $\tau_{k,ij}=\mu_k(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}- \frac{2}{3}div(U)\delta_{ij})$ where $\mu_k$ is the dynamic viscosity and $\delta_{ij}$, the Kronecker delta.​

$\sum^{Re}_{k,ij}$ is the turbulent stress tensor, $S_k$ is the external source term and $I_{p→k,i}$ represents the average interface momentum transfer rate from phase p to phase k, that accounts for mass transfer, drag force, added mass force, lift and satisfies the local balance equations

$I_{p→k,i}+ I_{k→p,i}=0$   (12)

NEPTUNE_CFD code includes several assumptions:

  1. the surface tension force is neglected,
  2. the average pressure of the two phases on either side of the interface are equal.

Furthermore, in this study, the pressure drop and turbulent stress tensor are neglected, the last one because the flow is considered laminar.

Energy equations

Considering the total enthalpy of the phase k

$H_k=e_k+\frac{1}{2}u_k^2+\frac{P}{\rho_k}$   (13)

The energy equation is written

$\frac{\partial}{\partial t}(\alpha_k \rho_k H_{k})+\frac{\partial}{\partial x_j}(\alpha_k \rho_k H_{k}U_{k,j}) = \frac{\partial}{\partial x_j} (\alpha_k \tau_{k,ij} U_{k,i})-\frac{\partial}{\partial x_j}\alpha_kQ'_{k,j}-\alpha_k\frac{\partial P}{\partial x_i}+\alpha_k \rho_k U_{k,i} g_{i}+\Pi_k+Q_{wall→k}+I_{p→k,i}$   (14)

Where $Q'_k=\lambda_kT_k$ and $\lambda_k$ represents the thermal conductivity.

$Q{wall →k}$ denotes the heat exchange with boundaries and is described by the nucleate boiling model. It takes into account bubble creation and satisfies:

 $\sum_{k=1}^{n phase}Q{wall→k}=Q_{wall}$   (15)

where $Q_{wall}$ is the total imposed heat flux.

Furthermore, $\Pi_k$ represents the bulk interface heat transfer, sum of the interface transfer between phase p and phase k, which complies with conservation relation

$\Pi_k= \sum_{p≠k} \Pi_{p→k}$   (16)


$\Pi_{p→k}+\Pi_{k→p}=0$   (17)

It should be noted that the code does not consider the  terms and $\alpha_k \tau_{k,ij} U_{k,i}$ and $I_{p→k,i}$

Nucleate boiling model

A first model for nucleation mechanisms similar to the one developed by Frost and Dzakowik (1967) [5] is implemented in NEPTUNE_CFDThe single phase temperature profile in the viscous sub-layer allows to obtain the relation between the wall superheat at boiling incipience ($T_{wall}-T_{sat}$) and the wall heat flux ($Q_{wall}$), used as criteria for estimating the point at which nucleate boiling starts. For, $r_{Cmax}$ the maximum cavity radius available for activation on the wall, two criteria can be considered:

$T_{crit1}=T_{wall}-T_{sat}=(\frac{8 \sigma T_{sat } Q_{wall}}{\lambda_l \rho_{sat}H_{lat}})^{1/2}$   (18)

$T_{crit2}=\frac {Q_{wall}}{\lambda_l}r_{Cmax}+\frac{2\sigma T_{sat}}{h_{lg}\rho_{sat}}\frac{1}{r_{Cmax}}$   (19)

And the activated cavity radius:

$r_{cl}=\frac{\lambda_l T_{crit1}}{2Q_{wall}}$   (20)

Starting the boiling process follows the subsequent pattern:

​if $r_{cl}\le r_{Cmax}$

$T_{wall}-T_{sat}=T_{crit1}$   (21)

Else, for $r_{cl}\ge r_{Cmax}$

$T_{wall}-T_{sat}=T_{crit2}$   (22)

Three-flux model

In a first simplified approach, the boiling heat flux is calculated following the analysis of Kurul and Podowski (1990) [6], ​, where $Q$ it is split into three terms, a single phase flow convective heat flux $Q_c$ at the fraction of the wall area unaffected by the presence of bubbles, a quenching heat flux $Q_q$ where bubble departures periodically bring cold water in contact with the wall and finally, a vaporization heat flux $Q_e$ needed to form the vapor phase (figure 4).

$Q=Q_c+Q_q+Q_e$   (23)

Figure 4: Analysis of Kurul and Podowski (1990) for the boiling heat flux

The wall surface unit is split into two parts: an area influenced by bubble departure, neglecting the overlapping areas of influence between adjacent bubbles, defined as $A_q=min(1,\frac{\pi D_m^2n}{4})$ and a single phase area $A_c$ with the relation $A_c+A_q=1$$D_m$ represents the departure bubble diameter [m] and $n$, the active site density [number/m2].

The convective heat flux is calculated with the classical law of single phase flow heat transfer at the wall

$Q_c=A_c h_{log}(T_{wall}-T_{\delta})$   (24)

where $\delta$ refers to a point in the fluid at a dimensionless distance from the wall and $h_{log}$ represents the heat transfer coefficient within the thermal boundary layer.

The quenching heat flux is modelled as the mean value of a transient conductive heat flux supplied to a medium at external temperature $T_{\delta}$, during the period $t_Q$ between the departure of a bubble and the beginning of growth of the following one

$Q_q=A_Qt_Qf \frac{2\lambda_l(T_{wall}-T_{\delta})}{(\pi\alpha_lt_Q)^{1/2}}$   (25)

$\alpha_l$ and $f$ represent the liquid thermal diffusivity and the bubble departure frequency, respectively.

Finally, the vaporisation heat is expressed as follows

$Q_e=fV_b\rho h_{lg}n$   (26)

Proportional to the active site density $n$ and the bubble's volume $V_b=\frac{\pi D_m^3}{6}$

In order to run the simulation, closures laws as the mean bubble equivalent diameter, the mean bubble departure frequency and the active site spatial density, should be provided. The usual closures laws used in NEPTUNE_CFD are the following ones:

  • Mean bubble equivalent diameter (Unal -1976 [7])

$R_d=\frac{1,2.10^5P^{0.709}a}{\sqrt{b\Phi}}$   (27)


  • $a=\frac{T_{wall}- T_{sat}}{2\rho_gh_{lg}} \sqrt{\frac {\lambda_{wall}\rho_{wall}Cp_p}{\pi}}$,   (28)
  • $b=\frac{T_{sat}-T_l}{2(1-\frac{\rho_g}{\rho_l})}$   (29)
  • $\Phi=max(1,(\frac{U_l}{0.61})^{0.47})$   (30)


  • Frequency (Ceumern & Lindenstjerna -1977 [8]):

    $f=\sqrt{\frac{2g(\rho_l-\rho_v)}{3\rho_lR_d}}$   (31)

$t_Q = 1/f$   (32)

  • Active site spatial density (Lemmert et Chwala -1977 [9])

$n=(210(T_{wall}-T_{sat}))^{1.8}$   (33)

Four-Flux Model

The basic wall heat flux partitioning model assumes that the amount of liquid on the wall is sufficient to remove heat from the wall to be used for evaporation. Thus, superheating of the vapour that occurs at high void fractions is not modelled. Under these assumptions, the basic heat flux partitioning model cannot be used, and a generalised model is considered. In this model, the total heat flow is considered to be composed of four different contributions, however, the criteria for onset of nucleated boiling remains unchanged.

$Q=Q_c+Q_{c2}+Q_q+Q_e$   (34)

The main differences between the three-flux and the four-flux models are listed below:

  • The monophasic heat flux in the liquid phase $Q_c$ is not weighted by $A_c$ (the area of the wall unaffected by the bubbles) but by the volume fraction of the liquid $\alpha_l$.
  • A new term $Q_{c2}$ is considered, which represents the convective heat transfer between the wall and the vapour phase. It is weighted by the vapour volume fraction $\alpha_v$.

$Q_{c2}=h_{wfv}(T_{wall}-T_{v})$   (35)

Where $h_{​wfv}$ is the wall heat transfer coefficient calculated from the wall temperature function for the vapour phase and $T_v$ is the vapour temperature at the centre of the wall adjacent cell.

  • In the previous model, $Q_q$ was calculated by considering the liquid as a semi-infinite medium. In this new model, it is weighted by the volume fraction of liquid into the first mesh, since it would not exist if there is no new cold liquid in contact the wall.
  • The density of active cavities is modified to take into account the different boiling regimes. This equation has been obtained by modifying the equation given by Kurul & Podowski for nucleate boiling to: $n=min((210(T_{wall}-T_{sat})^{1.8},n_{max})$ where $n_{max}$ corresponds to a surface completely covered with bubbles.

This model seems to be more robust than the tree-flux model, and thus to provide better results when high heat flux are considered. The extension of the wall-heat-flux-partitioning model was used also to take into account the critical heat flux condition.