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.
The multi-ﬁeld 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)
The general multi-ﬁeld 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 deﬁned 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:
Furthermore, in this study, the pressure drop and turbulent stress tensor are neglected, the last one because the ﬂow is considered laminar.
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)
where
$\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}$
A ﬁrst model for nucleation mechanisms similar to the one developed by Frost and Dzakowik (1967) [5] is implemented in NEPTUNE_CFD. The single phase temperature proﬁle 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 ﬂux ($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 ﬁrst simpliﬁed approach, the boiling heat ﬂux is calculated following the analysis of Kurul and Podowski (1990) [6], , where $Q$ it is split into three terms, a single phase ﬂow convective heat ﬂux $Q_c$ at the fraction of the wall area unaffected by the presence of bubbles, a quenching heat ﬂux $Q_q$ where bubble departures periodically bring cold water in contact with the wall and finally, a vaporization heat ﬂux $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 inﬂuenced by bubble departure, neglecting the overlapping areas of inﬂuence 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/m^{2}].
The convective heat ﬂux is calculated with the classical law of single phase ﬂow heat transfer at the wall
$Q_c=A_c h_{log}(T_{wall}-T_{\delta})$ (24)
where $\delta$ refers to a point in the ﬂuid at a dimensionless distance from the wall and $h_{log}$ represents the heat transfer coefficient within the thermal boundary layer.
The quenching heat ﬂux is modelled as the mean value of a transient conductive heat ﬂux 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:
$R_d=\frac{1,2.10^5P^{0.709}a}{\sqrt{b\Phi}}$ (27)
Where,
Frequency (Ceumern & Lindenstjerna -1977 [8]):
$f=\sqrt{\frac{2g(\rho_l-\rho_v)}{3\rho_lR_d}}$ (31)
$t_Q = 1/f$ (32)
$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:
$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.
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.