<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /></figure>

# Boiling-flow multiphase CFD simulations for nuclear reactor conditions without interfacial area transport equation

Corentin Reiss<sup>a,b,*</sup>, Antoine Gerschenfeld<sup>a</sup>, Catherine Colin<sup>b</sup>

## ARTICLE INFO

MSC:

80A20

76T10

Keywords:

Multiphase CFD

Interfacial forces

Interfacial heat flux

## ABSTRACT

We develop a two-fluid Euler-Euler CFD framework based on the PolyMAC numerical scheme (Gerschenfeld and Gorsse, 2022) in CEA's open-source TrioCFD code (Angeli et al., 2015). Interfacial momentum closure terms are selected and validated using bubbly adiabatic experiments on vertical flows (Colin et al., 2012; Hibiki et al., 2001). The local experimental bubble diameter is enforced to avoid the use of an interfacial area transport equation, as in Sugrue et al. (2017). Independently, it is shown that in a high-pressure developed boiling pipe flow, changing the entrance temperature while measuring flow characteristics at the outlet is equivalent to changing the distance from the inlet where the flow characteristics are measured. This enables us to simulate the DEBORA experiment (Garnier et al., 2001), an ascending boiling R12-freon flow in a tube, using a 3D map of the experimental diameter, avoiding distortions due to interfacial area modeling. We demonstrate that atmospheric-pressure closure terms are not able to reproduce measured void fraction profiles. We then consider that bubbles are deformed, i.e. non-spherical, in nuclear reactor conditions, characterized by high pressures, void fractions and flow velocities. Taking this deformation into account, we propose a new set of momentum and energy closures that is found to be independent of the bubble diameter. This enables us to run simulations without prior knowledge of the bubble diameter nor the need for an interfacial area transport equation or population balance model, as in system-scale codes extensively used in the nuclear industry (NRC, 2010; Berry et al., 2018). Void fraction predictions are more precise than with the baseline set of closures.

## 1. Introduction

Modeling multiphase flows is critical for nuclear applications (Delhaye, 2008; Todreas and Kazimi, 2021). However, these flows are extremely complex and a wide variety of flow patterns can exist (Ishii and Hibiki, 2006). Even restricting ourselves to bubbly flows, all configurations cannot be reliably simulated using existing computational fluid dynamics (CFD) codes (Lucas et al., 2015). Predicting the void fraction distribution in a boiling nuclear subchannel requires a knowledge of liquid turbulence, interfacial forces, wall boiling dynamics and transfers, condensation and bubble coalescence and fragmentation. Additionally, these terms interact with each other. For example, coalescence will affect bubble diameters, which will change their velocities and condensation rates, and the global mass, momentum and energy balances. This makes it very difficult to separate contributions and makes error compensation possible (Bestion et al., 2009).

In this paper, we contribute to building a CFD framework to study pressurized water reactor (PWR) conditions. Hosler (1967) and François et al. (2011) have shown that at PWR pressures and in PWR-similarity

conditions using a refrigerant fluid, the flow pattern remains bubbly, like an emulsion, at extremely high thermodynamic qualities ($X = 0.3$) and average vapor void fractions ($\langle\alpha_v\rangle = 0.7$). The two-fluid model is adapted to simulate such flows with a reasonable computational cost as it is not necessary to determine the precise position of the interfaces (Ishii and Hibiki, 2006). In the literature, due to the complexity of these flows and our inability to measure their precise dynamics, interfacial energy and momentum transfers closures are based on measures in simpler conditions (Chuang and Hibiki, 2017; Sugrue et al., 2017; Liao et al., 2019). We seek to evaluate the validity of these terms derived from single-bubble experiments (Tomiyama et al., 1998) or adiabatic flows (Sugrue, 2017) in high-pressure, high-void fraction and high mass flux conditions.

We implement standard single phase turbulence and bubbly flow closure terms (Chuang and Hibiki, 2017; Sugrue et al., 2017; Liao et al., 2019) (Section 2) and validate them on air-water bubbly flows in pipes at atmospheric pressure (Colin et al., 2012; Hibiki et al., 2001) (see Fig. 1 and Section 3). The prediction of the bubble diameters,

<footnote><sup>a</sup> Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique des Fluides, Gif-sur-Yvette, 91191, France</footnote>

<footnote><sup>b</sup> Institut de Mécanique des Fluides de Toulouse, Université de Toulouse, CNRS, INPT, UPS, Allée du Prof. Camille Soula, Toulouse, 31400, France</footnote>

<footnote>* Corresponding author at: Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique des Fluides, Gif-sur-Yvette, 91191, France.
E-mail addresses: corentin.reiss@cea.fr (C. Reiss), antoine.gerschenfeld@cea.fr (A. Gerschenfeld), catherine.colin@toulouse-inp.fr (C. Colin).</footnote>
Nomenclature

Subscripts

k Arbitrary phase

l Liquid phase

g Gas phase

v Vapor phase

s Saturation

b Bubble

w Wall

i Interfacial

j Number of a run in a test tube

r Radial

z Axial

Roman letters

A Area of a test section

A<sub>b</sub> Surface fraction of the wall occupied by bubbles

a<sub>i</sub> = 6α<sub>g</sub>/d<sub>b</sub> Interfacial area (m<sup>-1</sup>)

C<sub>p,k</sub> Heat capacity of phase k

C<sub>D</sub> Drag coefficient

C<sub>L</sub> Lift coefficient

C<sub>TD</sub> Turbulent dispersion coefficient

C<sub>V M</sub> Virtual mass coefficient

d<sub>b</sub> Bubble diameter, taken as Sauter mean diameter

d<sub>b, det</sub> Detachment bubble diameter

d<sub>dist</sub> Distance between two bubbles

D<sub>h</sub> Hydraulic diameter of the test section

e<sub>k</sub> Internal energy of phase k

F̃<sub>ki</sub> Interfacial momentum transfer to phase k (Nm<sup>-3</sup>)

f<sub>dep</sub> Bubble departure frequency

g̅ Gravity

G Flow mass flux (kgm<sup>-2</sup>s<sup>-1</sup>)

h<sub>lg</sub> Evaporation latent heat (Jkg<sup>-1</sup>)

H<sub>l,OSV</sub> Liquid wall heat transfer coefficient at onset of significant void (Wm<sup>-2</sup>K<sup>-1</sup>)

J<sub>k</sub> Superficial velocity of phase k (ms<sup>-1</sup>)

k Turbulent kinetic energy of liquid (m<sup>2</sup>s<sup>-2</sup>)

k<sub>+</sub> = k/u<sub>τ</sub><sup>2</sup> Non-dimensional turbulent kinetic energy

L Length of the test section

L<sub>c</sub> = √(σ/g(ρ<sub>l</sub> - ρ<sub>v</sub>)) Capillary length

N<sub>s</sub> Nucleation site density (m<sup>-2</sup>)

P Perimeter of a test section

q<sub>ki</sub> Interfacial heat flux towards phase k (Wm<sup>-3</sup>)

q<sub>kW</sub> Wall heat flux towards phase k (Wm<sup>-2</sup>)

q<sub>SP</sub> Single-phase wall heat flux (Wm<sup>-2</sup>)

q<sub>w</sub> Total wall heat flux

q<sub>w,l→v</sub> Evaporation wall heat

r<sub>+</sub> = r/R<sub>pipe</sub> Dimensionless radial position

P Pressure

P Heated perimeter of the test section

T<sub>k</sub> Temperature of phase k

T<sub>s</sub> Saturation temperature

T<sub>w</sub> Wall temperature

T<sub>in</sub> Inlet temperature

X Thermodynamic quality of the flow

X<sub>in</sub> Inlet thermodynamic quality

X<sub>j</sub> Outlet thermodynamic quality of run j

y Distance to wall

y<sub>+</sub> = yu<sub>τ</sub>/v<sub>l</sub> Dimensionless distance to wall

ũ<sub>k</sub> Velocity of phase k

u<sub>bulk</sub> Bulk liquid velocity

u'<sub>SP</sub> Turbulent fluctuations of a single-phase flow with the same mass flux as the studied flow

u<sub>τ</sub> Friction velocity at the wall

z Axial position along the test tube

Greek letters

α<sub>k</sub> Volume fraction of phase k

Γ<sub>k</sub> Bulk interfacial mass transfer towards phase k (m<sup>-3</sup>s<sup>-1</sup>)

Γ<sub>w,v</sub> Wall interfacial mass transfer towards phase k (m<sup>-2</sup>s<sup>-1</sup>)

λ<sub>k</sub> Thermal conductivity of phase k (Wm<sup>-1</sup>K<sup>-1</sup>)

μ<sub>k</sub> Dynamic viscosity of phase k (Pas)

ν<sub>k</sub> Kinematic viscosity of phase k (m<sup>2</sup>s<sup>-1</sup>)

ν<sub>t</sub> Turbulent kinetic viscosity of liquid (m<sup>2</sup>s<sup>-1</sup>)

ρ<sub>k</sub> Volume mass of phase k

σ Surface tension (Nm<sup>-1</sup>)

τ<sub>w,f</sub> Shear stress at the wall (Nm<sup>-2</sup>)

ω Turbulent dissipation frequency of liquid

ω<sub>+</sub> = ωu<sub>τ</sub>/v Non-dimensional dissipation rate

Dimensionless numbers

b = α<sub>g</sub> ||ũ<sub>g</sub> - ũ<sub>t</sub>||<sup>2</sup>/u<sub>SP</sub><sup>2</sup> Bubblance parameter

E<sub>o</sub> = (ρ<sub>l</sub> - ρ<sub>g</sub>)g<sub>b</sub>g<sub>b</sub>

Ja = ρ<sub>l</sub>C<sub>pl</sub>||T<sub>sat</sub> - T<sub>l</sub>||/ρ<sub>g</sub>h<sub>lg</sub>

Nu Bubble Nusselt number

Pe = d<sub>b</sub>||ũ<sub>g</sub> - ũ<sub>t</sub>||/ρC<sub>pl</sub>

Pr = v<sub>t</sub>·ρC<sub>pl</sub>/λ<sub>l</sub> liquid Prandtl number

Re<sub>b</sub> = d<sub>b</sub>||ũ<sub>g</sub> - ũ<sub>t</sub>||/v<sub>t</sub>

Re = D<sub>h</sub>u<sub>bulk</sub>/v<sub>t</sub> Bulk Reynolds number

We = ρ<sub>l</sub>d<sub>b</sub>||ũ<sub>g</sub> - ũ<sub>t</sub>||<sup>2</sup>/σ Weber number

We<sub>v</sub> = ρ<sub>l</sub>d<sub>b</sub>(cd<sub>b</sub>)<sup>2/3</sup>/σ Turbulent Weber number

Wo = Eo||ũ<sub>g</sub> - ũ<sub>t</sub>||<sup>2</sup>/k Wobble number

through interfacial area transport equations (Yao and Morel, 2004) or population balance (Krepper et al., 2008; Yuan et al., 2012) methodologies, is a major source of uncertainty and errors in two-fluid CFD

simulations (Leoni, 2022). We therefore use bubbly flow datasets with measured mean Sauter diameters and enforce the experimental values as in Sugrue (2017).

We then study the DEBORA experiment (Garnier et al., 2001) (see Fig. 1). This is a vertical heated tube filled with flowing Freon-12 designed to fill similarity criteria with PWR and boiling water reactor conditions. Void fraction, mean Sauter diameter and temperature measures were taken at the outlet for different pressures, flow rates,
<table><caption>Fig. 1. Diagrams of the two-phase experiments used for framework validation and heat transfer models evaluation.</caption><thead><tr><th>Experiment</th><th>Measured quantities</th><th>Models validated</th></tr></thead><tbody><tr><td>Hibiki et al. (2001)<br>Colin et al. (2012)</td><td>Gas & Liquid velocity<br>Gas & Liquid vel. fluct.<br>Void fraction<br>Sauter diam.</td><td>Lift force<br>Turbulent disp. force<br>Wall repulsion force</td></tr><tr><td>DEBORA<br>(Garnier et al. (2001))</td><td>Gas velocity & vel. fluct.<br>Liquid temperature<br>Void fraction<br>Sauter diam.</td><td>Heat flux partition<br>Interf. heat transfer</td></tr></tbody></table>

heat fluxes and inlet subcoolings. We show that in a high-pressure developed flow, changing the entrance temperature in a boiling pipe flow is equivalent to changing the distance from the entry where the flow characteristics are measured (Section 4).

This enables us to use a 3D map of the experimental diameters as simulation inputs to evaluate the standard closure terms validated at atmospheric pressures independently of the bubble diameter modeling (Section 5). These yield unsatisfying results, which shows that interfacial force closure laws validated in atmospheric-pressure adiabatic flows are no longer valid in pressurized water reactor similarity conditions even if bubble diameters are well predicted.

To improve the simulation results we show that in nuclear reactor conditions, even small bubbles are deformed, i.e. non-spherical. Accounting for this, we build a set of momentum and energy closures that improves simulation predictions compared with the baseline set of closures (Section 6). These new closure laws are found to be independent of the bubble diameter. We therefore do not require experimental diameter inputs or an interfacial area transport equation to simulate PWR conditions.

## 2. Physical model

### 2.1. Numerical framework and conservation equations

CEA is developing a multiphase Reynolds-averaged Navier-Stokes CFD module in its open-source CFD code, TrioCFD (Angeli et al., 2015). This module is based on the PolyMAC finite volume numerical scheme developed by Gerschenfeld and Gorsse (2022) for component-scale codes. PolyMAC is a staggered scheme that can handle arbitrary polyhedral conform or non-conform elements. Depending on the application, structured and unstructured meshes that can have prismatic near-wall layers are generated using the SALOME platform (Bergeaud and Lefebvre, 2010). The TrioCFD multiphase module solves mass, momentum and energy conservation equations for an arbitrary number of fluids in an Euler-Euler framework (Ishii and Hibiki, 2006). The semi-conservative form of the momentum equation is used (Park et al., 2009). The equations that govern a phase k are:

$$ \frac{\partial a_k \rho_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k \vec{u}_k) = \Gamma_k \\ \alpha_k \rho_k \frac{\partial \vec{u}_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k \vec{u}_k \otimes \vec{u}_k) - \vec{u}_k \nabla \cdot (\alpha_k \rho_k \vec{u}_k) = -\alpha_k \nabla P + \nabla \cdot [\alpha_k \mu_k \nabla \vec{u}_k - \alpha_k \rho_k \overline{u'_i u'_j}] + \vec{F}_{ki} + \alpha_k \rho_k \vec{g} \quad (1) $$

$$ -P \left( \partial_t \alpha_k + \nabla \cdot (\alpha_k \vec{u}_k) \right) + \nabla \cdot \left[ \alpha_k \lambda_k \nabla T_k - \alpha_k \rho_k \overline{u'_i e'_k} \right] + q_{ki} $$

Where $\alpha_k$ is the volume fraction of phase k, $\rho_k$ its density, $\vec{u}_k$ its velocity, P the local pressure, $\mu_k$ the dynamic viscosity, $\vec{g}$ the gravity vector, $e_k$ the internal energy, $\lambda_k$ the fluid conductivity and $T_k$ the fluid temperature.

In Eq. (1), the terms that need closure laws are the turbulent terms $u'_i u'_j$ and $u'_i e'_k$, the mass transfer term $\Gamma_k$, the interfacial forces accounting for momentum transfers between phases per unit of volume and time $\vec{F}_{ki}$ and the interfacial heat transfer $q_{ki}$. The wall heat transfer to phase k, $q_{ki}$, which is a boundary condition of the energy equation, also requires modeling.

In this paper, we work with two fluids: a continuous liquid phase l and a dispersed phase, written g for gas phase in adiabatic flow and v for vapor phase in boiling flow.

### 2.2. Turbulence modeling

**Shear-induced turbulence.** In rod bundles, two-equations turbulence models yield similar results as Reynolds stress models with faster calculation times (Franck et al., 2012). We therefore select the Kok (1999) $k-\omega$ turbulence model as it has similar properties to the more commonly used Menter (1993) $k-\omega$ model, i.e., a cross-diffusion term that is suppressed in the near-wall region, but is easier to implement. This yields:

$$ \begin{gathered} v_t = \frac{k}{\omega} \overline{u'_i u'_j} = -v_t \nabla \vec{u}_i \qquad \overline{u'_i e'_j} = -v_t C_{pl} \nabla T_i \\ \partial_t (k + \nabla \cdot (k \vec{u}_i)) = v_t (\nabla \vec{u}_i + \overline{u'_i}) \cdot \nabla \overline{u'_i} - \beta_k k \omega + \nabla \cdot (\alpha_l (v_l + \sigma_k v_t) \nabla k) \\ \partial_t \omega + \nabla \cdot (\omega \vec{u}_i) = \alpha_\omega (\nabla \vec{u}_i + \overline{\nabla \vec{u}_i}) \cdot \nabla \overline{u'_i} - \beta_\omega \omega^2 \\ + \nabla \cdot (\alpha_l (v_l + \sigma_\omega v_t) \nabla \omega) + \sigma_d \frac{1}{\omega} \max\{\nabla k \cdot \nabla \omega, 0\} \end{gathered} \quad (2) $$

Where k is the turbulent kinetic energy of the liquid, $\omega$ its turbulent dissipation rate, $v_t$ its turbulent viscosity, $v_l$ its kinematic viscosity and $C_{pl}$ is the heat capacity of the liquid. The values of the constants are $\alpha_\omega = 0.5$, $\beta_\omega = 0.09$, $\sigma_\omega = 0.075$, $\sigma_k = 2/3$, $\sigma_\omega = 0.5$ and $\sigma_d = 0.5$.

We implement an adaptive wall-law algorithm that begins by determining the friction velocity $u_\tau$ in the same way as in Carlson et al. (2015). The shear stress at the boundary is then computed and is used as a Navier boundary condition for the momentum equation: $\tau_{wf} = \alpha_l \rho_l \overline{u_\tau^2}$.

The boundary condition on k is $k = 0$ at the wall for $y_+ < 5$, where $y_+ = y_\tau / v_l$ and y is the distance between the wall and the first element center. For larger wall elements, it is a zero-flux condition. The transition is smoothed by a transition factor tanh $(y_+/10)^2)$.

For $\omega$, Knopp et al. (2006) give an analytical value in the near-wall region. A simple solution would be to enforce this value in the
first element. However, it was already used in TrioCFD and creates numerical issues for tetrahedron meshes. Instead, we calculate the analytical solution at a distance y/2 from the wall. We then enforce a fixed boundary condition at the wall: $\omega_{wall} = 2 \cdot \omega(y/2)$. This amounts to creating a virtual element between the first element and the wall in which we know the value of $\omega$.

The single-phase heat transfer coefficient that we have implemented is the one proposed by Kader (1981). He gives an expression of the wall heat flux towards the liquid $q_{wl}$:

$$ 
\begin{aligned}
q_{wl} &= (T_w - T_l(y)) \frac{\rho_l C_{pl} u_r}{\Theta_u^\omega(y_+)} \\
\beta_{SP} &= (3.85(Pr^{1/3}) - 1.3)^2 + 2.12 \log(Pr) \\
\gamma &= \frac{0.01(Pr_y)_+^4}{1 + 5Pr^3 y_+}
\end{aligned} 
\quad (3) $$

$$ \Theta_{+}^{\omega}(y_{+}) = Pr_{+} \exp(-\gamma) + (2.12 \cdot \log(1 + y_{+}) + \beta_{sp}) \exp(-\frac{1}{\gamma}) $$

Where $T_w$ is the wall temperature and $Pr = \frac{\nu_l \rho C_{pl}}{\alpha_l}$ the Prandtl number. This expression is based on experimental measures for $y_+$ ranging from 0 to 300. Using the previously calculated expression of $u_r$, we calculate $\theta_+$ and a wall heat transfer coefficient. We use the Kader expression as the convective heat transfer contribution in our wall heat flux partition model (see Section 2.4).

**Bubble-induced turbulence.** According to Almeras et al. (2017), liquid velocity fluctuations induced by bubble movements are small before those from shear if the bubblance parameter $b = \frac{\alpha_g \|\vec{u}_g - \vec{u}_l\|^2}{u_{sp}^2} < 0.5$, where $u'_{sp}$ are the turbulent fluctuations for a single-phase flow with the same mass flux. At PWR pressures, the bubble diameter $d_b < 1$ mm, which yields $\|\vec{u}_g - \vec{u}_l\| < 0.1$ m/s, and $\alpha_g < 0.5$. Using the Reichardt (1951) correlation for a bulk Reynolds number $Re = 10^5$, bubble-induced turbulence can be neglected in PWR's (i.e. $b < 0.5$) for bulk velocities $u_{bulk} > 1.5$ m/s. This is a low value as $u_{bulk} \sim 4$ m/s in operation (Delhaye, 2008), and $b \sim 0.1$. Therefore, we do not model bubble-induced turbulence.

## 2.3. Interfacial momentum exchanges modeling

The interfacial force exerted by the liquid on the gas is $\vec{F}_{gi} = -\vec{F}_{li}$. In this subsection, all forces written apply to the gas phase and are written in force per volume unit. We separate the interfacial momentum transfer term in five different contributions: the drag, virtual mass, lift and turbulent dispersion forces, and the wall correction.

$$ \vec{F}_{gi} = \vec{F}_{drag} + \vec{F}_{VM} + \vec{F}_{lift} + \vec{F}_{TD} + \vec{F}_{wall} \quad (4) $$

Apart from a modification to the virtual mass force, the momentum closures that we select come from Sugrue (2017).

**Drag force.** We implement the contaminated drag force of Tomiyama et al. (1998):

$$ 
\begin{aligned}
\vec{F}_{drag} &= -a_l \frac{1}{2} C_D \rho_l \|\vec{u}_g - \vec{u}_l\| (\vec{u}_g - \vec{u}_l) \\
&= -\frac{3}{4} C_D \frac{\alpha_g \rho_l}{d_b} \|\vec{u}_g - \vec{u}_l\| (\vec{u}_g - \vec{u}_l)
\end{aligned} 
\quad (5) $$

$$ C_D = \max \left( \frac{24}{Re_b} (1 + .15 Re_b^{0.87}), \frac{8E_0}{3E_0+12} \right) $$

Where $C_D$ is the drag coefficient, $a_l = \frac{6\alpha_g}{d_b}$ is the interfacial area of the gas phase, $E_0 = \frac{(\rho_l - \rho_g)g d_b^2}{\sigma}$ the Eotvos number and $Re_b = \frac{d_b \|\vec{u}_g - \vec{u}_l\|}{v_l}$ the bubble Reynolds number.

**Virtual mass force.** The virtual mass force writes:

$$ \vec{F}_{VM} = -C_{VM} \rho_l (\partial_t \vec{u}_g - \partial_t \vec{u}_l + \vec{u}_g \nabla \vec{u}_g - \vec{u}_l \nabla \vec{u}_l) \quad (6) $$

Where $C_{VM}$ is the virtual mass coefficient. It can be interpreted as the volume of liquid entrained by the accelerating or decelerating bubbles.

The most commonly used formulations are the constant coefficient $C_{VM} = \frac{1}{2} \alpha_g$ and the $C_{VM} = \frac{1}{2} \frac{1+2\alpha_g}{1-\alpha_g} \alpha_g$ (Zuber, 1964), that were both

derived theoretically. Recently, Béguin et al. (2016) performed potential flow simulations with random bubble positions and found $C_{VM} = \alpha_g (\frac{1}{2} + 0.34\alpha_g^2) \sim \frac{1}{2} \alpha_g$. Furthermore, in some DEBORA experimental runs (Garnier et al., 2001), the local void fraction can reach 0.7. If we used the standard constant coefficient $C_{VM} = \frac{1}{2} \alpha_g$, the liquid volume fraction entrained by the gas phase would be $0.7 \cdot 1/2 = 0.35$. The total liquid fraction, which is $\alpha_l = 1 - \alpha_g = 0.3$, would be less than the entrained liquid fraction. This is non-physical and leads to numerical stability issues. We therefore assume that at most 1/2 of the remaining liquid can be entrained by the gas, the value 1/2 being arbitrary. This leads to:

$$ C_{VM} = \min \left( \frac{1}{2} \alpha_g, \frac{1}{2} \alpha_l \right) \quad (7) $$

Our modification affects $C_{VM}$ for $\alpha_g > 0.5$.

**Lift force.** The general formulation for the lift force is:

$$ \vec{F}_{lift} = -C_L \rho_l \alpha_g (\vec{u}_g - \vec{u}_l) \wedge (\nabla \wedge \vec{u}_l) \quad (8) $$

The difference between lift force models is the lift coefficient $C_L$. A constant coefficient can be chosen by the user. The Sugrue (2017) formulation was also implemented, as it was designed to operate on high-void fraction ascending flows and not only single bubbles, contrarily to the Tomiyama et al. (2002) formulation. The Sugrue (2017) lift coefficient requires a so-called Wobble number $Wo$ and reads:

$$ 
\begin{aligned}
C_L &= f(Wo) \cdot g(\alpha) \quad , \quad g(\alpha) = \max(0, 1.0155 - 0.0154 \exp(8.0506 \alpha)) \\
Wo &= E_0 \frac{k}{\|\vec{u}_g - \vec{u}_l\|^2} \quad , \quad f(Wo) = \min(0.03, 5.0404 - 5.0781 Wo^{0.0108})
\end{aligned} 
\quad (9) $$

**Turbulent dispersion force.** We select the Burns et al. (2004) force :

$$ \vec{F}_{TD} = -C_{TD} \rho_l k \nabla \alpha_g \quad , \quad C_{TD} = \frac{3}{4} \frac{C_D}{d_b} |\vec{u}_g - \vec{u}_l| \frac{1}{\omega} \left( 1 + \frac{\alpha_g}{\alpha_l} \right) \quad (10) $$

Where $C_{TD}$ is the turbulent dispersion coefficient.

**Wall correction.** The main wall correction term implemented is the one proposed by Lubchenko et al. (2018). It is based on geometrical arguments. It suppresses lift and modifies turbulent dispersion close to the wall. The lift coefficient becomes:

$$ 
\begin{aligned}
C_L &\rightarrow \begin{cases} 0 & \text{ify}/d_b < 1/2 \\ C_L \left(3 \left(\frac{2y}{d_b} - 1\right)^2 - 2 \left(\frac{2y}{d_b} - 1\right)^3\right) & \text{if } 1/2 \le y/d_b < 1 \\ C_L & \text{ify}/d_b \ge 1 \end{cases} 
\end{aligned} 
\quad (11) $$

If $\vec{n}$ is the unit vector normal to the wall, the turbulent dispersion wall correction writes:

$$ \vec{F}_{wall} = \begin{cases} C_{TD} \rho_l k \cdot \alpha_g \frac{1}{y} \frac{d_b - 2y}{d_b - y} \vec{n} & \text{ify}/d_b < 1/2 \\ 0 & \text{ify}/d_b \ge 1/2 \end{cases} \quad (12) $$

## 2.4. Heat and mass transfers

**Interfacial heat transfer.** As we study bubbly flows, we select interfacial heat transfer formulations based on the calculation of the Nusselt number $Nu$ to describe the heat transfer from the liquid phase to the interface:

$$ q_{li} = a_l \frac{\lambda_l}{d_b} (T_s - T_l) Nu = \frac{6\alpha_v}{d_b} \frac{\lambda_l}{d_b} (T_s - T_l) Nu \quad (13) $$

the Ranz and Marshall (1952) model is implemented for condensation:

$$ Nu = 2 + 0.6 Re_b^{1/2} Pr^{1/3} \quad (14) $$

To the best of our knowledge, no liquid temperature was ever measured above the saturation temperature in flow boiling experiments (Roy et al., 2002; Garnier et al., 2001; Francois et al., 2021). Therefore, we do not allow the liquid to overcome the saturation
temperature by enforcing an extremely high heat transfer coefficient ($q_{ki} = 10^8$ Wm$^{-3}$) if $T_i > T_s$.

The heat transfer from the vapor phase to the interface also uses an extremely high heat transfer coefficient, as no vapor has been measured hotter than $T_s$ in bubbly flow:

$$q_{vi} = 10^8(T_s - T_v) \qquad (15)$$

**Wall heat transfer.** The original (Kurul and Podowski, 1990) model was selected. Though more recent and complex formulations have been proposed (Basu et al., 2005; Kommajosyula, 2020; Favre, 2023), it is used as the reference in the literature. It reads:

$$q_{KP} = q_c + q_q + q_{w,l \to v} \qquad (16)$$

$$q_{wl} = q_c + q_q$$

Where $q_c$ is the convective heat flux, $q_q$ the quenching heat flux, $q_{w,l \to v}$ the evaporation heat flux and $q_{wl}$ the total heat transfer towards the liquid. The nucleation site density $N_s$ is from Del Valle and Kenning (1985), the detachment bubble diameter $d_b, \det$ is a linear interpolation between those of Ünal (1976) and Thomas (1981), and the departure frequency $f_{dep}$ was proposed by Cole (1960). The surface fraction occupied by bubbles $A_b$ is calculated from these quantities, and the single-phase heat transfer $q_{SP}$ comes from Kader (1981) (Eq. (3)):

$$N_s = (210(T_w - T_s))^{1.8} \quad , \quad d_b, \det = 10^{-4} \cdot (T_w - T_s) + 0.0014$$

$$f_{dep} = \sqrt{\frac{4}{3} \frac{g(\rho_l - \rho_s)}{\rho_l d_b, \det}} \quad , \quad A_b = \min(1, \pi/4 \cdot N_s d_b^2, \det) \qquad (17)$$

This yields:

$$q_c = (1 - A_b)q_{SP} \qquad (18)$$

$$q_q = 2A_b \lambda_l (T_w - T_l) \sqrt{\frac{f_{dep} \rho_l C \rho_l}{\pi \lambda_l}}$$

$$q_{w,l \to v} = \frac{\pi}{6} f_{dep} d_b^3 \det \rho_g h_{lg} N_s$$

**Mass transfer.** The mass transfer between phases is entirely determined using the interfacial heat transfers. The mass transfer towards the vapor phase in the bulk of the flow $\Gamma_v$ and at the wall $\Gamma_{w,v}$ read:

$$\Gamma_v = -\Gamma_l = \frac{q_{vi} - q_{li}}{h_v - h_l} \quad , \quad \Gamma_{w,v} = \frac{q_{w,l \to v}}{h_{lg}} \qquad (19)$$

### 3. Framework validation

#### 3.1. Adiabatic single-phase flow

We run simulations of flow in a same 2D channel at $Re = 20,000$ with cartesian grids of varying refinements: $y_+$ in the first element ranges from 3 to 229 (see Fig. 2). The results are independent of $y_+$ in the first element and consistent with literature and refined solutions, except for $k_+$ in the near-wall region which is expected as we are transitioning from a wall-resolved to a wall-modeled solution. We also simulate pipe flow experiments with various tube diameters and fluids, from Colin et al. (2012) (see Fig. 3-A,B). The mesh used in these pipe simulations is shown at the bottom left of Fig. 4. It contains 80 vertical elements and 14 radial elements with mesh grading. The simulated mean velocity and velocity fluctuations match experimental results. This validates our implementation of single-phase turbulence and adaptive wall laws.

#### 3.2. Heated single-phase flow

One of the campaigns on the DEBORA experiment consisted in measuring the liquid R12 temperature for different flow rates, pressures and heating power (Garnier et al., 2001). We simulate single-phase heated flows from this campaign (see Fig. 3-C). The mesh used in these simulations is shown at the top left of Fig. 4. It contains 200 vertical elements and 20 radial elements of the same size. The liquid temperature from the experiment and the simulation have the same profile shape, though they are off by ~ 1 °C. This amounts to a 5% power loss on the

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** $k-\omega$ results for a turbulent 2D channel with $Re = 20,000$. The size of the first element at the wall, and therefore $y_+$, varies in each simulation. A. Non-dimensional velocity $u_+ = u/u_-$ as a function of $y_+$. Black line: Reichardt (1951) solution. B. Non-dimensional turbulent kinetic energy $k_+ = k/u_+^2$ as a function of $y_+$. Black line: refined solution for $y_{+,1} = 1$. C. Non-dimensional dissipation rate $\omega_+ = \omega u_+ / v$ as a function of $y_+$. Black line: Knopp et al. (2006) solution.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Single-phase validation of the multiphase module of TrioCFD. A. Colin et al. (2012) single-phase axial velocity in a 4 cm-diameter pipe. B. Colin et al. (2012) single-phase axial velocity fluctuations in a 4 cm-diameter pipe. C. Temperature profiles from single-phase DEBORA runs for a 19.2 mm diameter pipe, 2.61 MPa pressure, 1996 kg/(m$^2$ s) mass velocity and 74.4 kW/m$^2$ heat flux with inlet temperatures of 19.84 °C and 31.46 °C (page 109 in Garnier et al. (2001)).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Different meshes used in this paper. Top left: 2° slice with 20 radial elements of the same size. Bottom left: 45° slice with 14 radial elements and mesh grading. Right: 360° tube with 14 radial elements and mesh grading.</figcaption></figure>

experimental setup. In our boiling DEBORA simulations, we reduce the power boundary condition by 5% to account for this. A similar loss was found on the DEBORA experiment by Gueguen (2013) for single-phase flow and Favre (2023) for boiling flow.

### 3.3. Mesh refinement

To evaluate the mesh sensitivity of TrioCFD, we use test case U1 from Colin et al. (2012). Using SALOME software (Bergeaud and Lefebvre, 2010), we mesh a disk with quadrilaterals and extrude it to obtain a hexahedral pipe mesh. We run simulations on a full cylinder, a quarter of a cylinder and an eighth of a cylinder with symmetry boundary conditions on vertical planes, and on a 2° slice only one element wide. These meshes are shown in Fig. 4. The 2° slices have 20 or 40 radial elements of the same size. The others can have 7, 14 or 28 radial elements. They have 40, 80 or 160 vertical elements. There is a significant difference between 7 and 14 radial element results, but virtually none between 14 and 28. Calculations that run on a cylinder

or a slice give identical results. The details of the configuration used for the mesh refinement are given in Appendix A, along with Fig. A.21 that contains the mesh refinement results.

### 3.4. Two-phase adiabatic vertical tube

Multiple experimental databases are available to study adiabatic two-phase pipe flow. To validate the multiphase module, we select the Hibiki et al. (2001) database for upwards flow, as it covers a broad range of liquid and gas injection fluxes. The geometry used is a 3.06 m-long 5.08 cm-diameter pipe. We also select the Colin et al. (2012) experiments for downwards and microgravity flow in a 3.08 m-long 4 cm-diameter pipe. We run different test cases using the interfacial force models described in Section 2.3. To avoid modeling the interfacial area, we enforce the radially-dependent steady-state experimental diameter in the simulations. The diameter profiles of the presented simulations can be found in Fig. B.22, in Appendix B. The mesh used is a 45° slice with 14 radial elements, 80 axial elements and mesh grading (bottom left in Fig. 4). Air and water are injected at the bottom of the pipe at the mixture bulk velocity, with the void fraction set to respect experimental superficial velocities.

Simulation results are shown in Fig. 5. The complete model is able to predict correctly void fraction profiles for low (Fig. 5-A) and high (Fig. 5-C) liquid fluxes, in wall-peaked and core-peaked situations respectively. The prediction of the transition between both regimes can still be improved, as can be seen in Fig. 5-B. Furthermore, in Fig. 5-D the gas velocities are well predicted by the model. The downwards and micro-gravity profiles are also well predicted (Fig. 5-E and F).

## 4. The DEBORA database and the test tube hypothesis

The DEBORA loop was built at CEA/Grenoble to study boiling flow in reactor conditions by filling similarity criteria (Cubizolles, 1996; Garnier et al., 2001) (see Fig. 1). It consists of a vertical Freon-12

<figure><img src="image_6.png" /><figcaption>Fig. 5. Two-phase adiabatic simulated and experimental void fraction distributions in a pipe. We use the interfacial force models described in Section 2.3. $J_l$: superficial liquid velocity. $J_g$: superficial gas velocity. Lines: simulations. Circles: experimental results. A-C. Void fraction results in Hibiki et al. (2001) upwards flow experiments. D. Axial gas velocity results in Hibiki et al. (2001) upwards flow experiments. E. Void fraction results in Colin et al. (2012) downwards flow experiments. F. Void fraction results in Colin et al. (2012) microgravity experiments.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Diagram of the transition from thermodynamic quality to equivalent altitude. In run 1, fluid enters at quality $X_1$, and measures are performed at an altitude $L$. This is akin to measuring at an altitude $z_{1,\text{equivalent}} < L$ when the entrance quality is $X_2 > X_1$.</figcaption></figure>

pipe with a 1 m-long inlet section, a 3.5 m-long heated section and an instrumentation plane located 3.485 m after the beginning of the heated section. During different campaigns, liquid temperatures, void fractions, gas velocities and Sauter mean diameters were measured. Results from DEBORA have been extensively used to validate CFD boiling flows for nuclear applications (Končar et al., 2011; Favre et al., 2022; Vlček and Sato, 2023; Pham et al., 2023).

In the following, we will call a test tube a series of runs that have identical outlet pressures, mass flux and heating power. Each test tube consists of different runs in which the inlet temperature is varied. They are labeled G[ng]P[np]W[nw], where ng is the average mass flux (10³ kg/(m²s)), np the average pressure (bar) and nw the average heating power in the experiment (kW).

## 4.1. The test tube hypothesis

We call test tube hypothesis the assumption that reducing the inlet temperature is equivalent to measuring physical quantities at a lower point in the test section. This assumption has been used to simulate the DEBORA database in the past (Pham et al., 2023), but to the best of our knowledge its validity range has not been evaluated. In this section, we show that this hypothesis is accurate for the DEBORA database and for PWR conditions. This enables the use of Sauter mean diameter measures at the outlet for different inlet temperatures to enforce a 3D map of experimental diameters in simulations (see Section 5).

We assume a test section that has an area $A$, a mass flux $G$, a heated perimeter $P$, a length $L$ and a heat flux $q_w$. $z$ is the position along the test tube, and the instrumentation is located at the outlet at $z = L$. We assume two test runs, with inlet thermodynamic qualities $X_{1,in}$ and $X_{2,in}$ (see Fig. 6). Latent heat is $h_{lg}$. Then, through a simple heat balance:

$$X(z) = X_{in} + \frac{P q_w}{G A h_{lg}} z \quad (20)$$

And the equivalent altitude of run 1 for the conditions of run 2 is:

$$z_{1,\text{equivalent}} = L - \frac{G A h_{lg}}{P q_w} (X_{1,in} - X_{2,in}) \quad (21)$$

A diagram of this situation is presented in Fig. 6.

To assess the validity of the test tube hypothesis, we carry out two steps:

* a priori analysis: we evaluate the impact of two physical mechanisms that can have an impact on the validity of the hypothesis: the formation of the turbulent thermal boundary layer and the effect of hydrostatic pressure on saturation temperature (Section 4.2).

* a posteriori analysis: we run CFD simulations of different DEBORA test tubes and of a single tube with measures at equivalent altitudes (Section 4.3).

## 4.2. A priori analysis: single-phase turbulent boundary layers and saturation temperature

Single-phase turbulent boundary layers. In this paragraph, we analyze results from the DEBORA experiments (Cubizolles, 1996; Garnier et al., 2001). We calculate an equivalent altitude for each test of the database, using as a reference the runs with the highest inlet temperatures. For the boiling dynamics in the core of the flow to be similar at the outlet and at the equivalent altitude, the velocity and temperature fields must be the same at the point where vapor production substantially affects the flow, which is at the onset of significant void (OSV) (Saha and Zuber, 1974). The equivalent altitude methodology guarantees identical average enthalpies on the cross-section. However, it does not guarantee identical velocity and temperature profiles. For profiles to be the same, single-phase temperature and velocity profiles before the OSV must be developed.

Labunstov et al. (1974) conducted average void fraction measures in flow boiling after heated sections of varying lengths and found no difference for identical thermodynamic qualities when $z/D > 20$. In single-phase adiabatic flows, velocity profiles are developed after $z/D_h = 40$ (Doherty et al., 2007). As there is a non-heated section where the flow can develop before the heated section this criterion is always met. When the velocity profile is already developed, Abbrecht and Churchill (1960) show that the difference between the thermal boundary layer and the developed boundary layer is small when $z/D_h = 10$. Using the Al-Arabi (1982) correlation, the difference between the developed and developing Nusselt number is smaller than 3% at $z/D_h = 30$. $z/D_h = 30$ is therefore taken as an upper boundary of the reference length for the development of the thermal boundary layer in a turbulent flow. The vertical black lines are placed at that location in Fig. 7. We also plot the void fraction as a function of the equivalent altitude for each test tube. We can see that in the DEBORA database, OSV occurs for $z/D_h > 30$ and the developed turbulence at OSV criterion is verified. This may not be the case in other databases, in particular when the heated section is much shorter.

Saturation temperature. In Fig. 7, we can see that the equivalent altitude of some DEBORA cases is 3 m lower than the measuring plane. The hydrostatic pressure variation between the measuring plane in the experiment and the equivalent altitude plane in the simulation at these points is 0.4 bar. This changes the saturation temperature of R12-freon by ~1 °C (Bell et al., 2014). This is much smaller than the typical 20 °C inlet subcoolings (Garnier et al., 2001), and the difference in saturation temperature due to pressure variations should not affect the simulation.

However, if the experimental outlet pressure had been 1 bar, the difference between both pressures would have been 10 °C, which would significantly affect the flow. For the equivalent altitude hypothesis to be valid, the pressure difference between the different equivalent altitudes must be small enough for the saturation temperature variation to remain small with respect to typical subcoolings.

Use for experiment design. To guarantee that the test tube hypothesis can be used in simulations, we recommend the following for future experiments:

* Have a sufficiently long heated section to guarantee developed single-phase flows before the onset of significant void

* For a given test tube, reduce the hydrostatic pressure difference between the experimental outlet pressure and the one at the equivalent altitude. This can be done by increasing the outlet pressure by the hydrostatic pressure difference when the inlet temperature is reduced. This is key in low-pressure experiments ($P < 10$ bar).
<figure><img src="image_8.png" /><figcaption>Fig. 7. Void fraction as a function of equivalent altitude in the DEBORA 29 campaign. Each subplot represents a test tube. Vertical black lines represents $z/D_h = 30$.</figcaption></figure>

### 4.3. A posteriori analysis: numerical verification

In order to verify that the equivalent altitude hypothesis does not affect calculation results in the DEBORA experiment, we run two sets of simulations for a single test tube. The models presented in Section 2 are selected, with a 1 mm constant Sauter mean diameter chosen arbitrarily. The mesh is a 2°-wide slice, with 20 radial and 500 axial cells (top left in Fig. 4).

In the first set of simulations, the entrance temperature is the experimental temperature of a run. Physical quantities are extracted from the simulation at the outlet. This mimics the experimental conditions.

In the second set, the entrance temperature is the hottest experimental temperature in the test tube. Physical quantities are extracted from the equivalent altitudes for the experimental runs with lower entrance temperature.

We select the G2P14W16 series shown in Fig. 7, as it could be affected by both issues mentioned in the previous section. This series has among the earliest void production in the database. Furthermore, the outlet pressure is 14 bar, which should make it more sensitive to the saturation temperature variation evoked earlier than other series. We select 5 inlet temperatures for the simulation: 23.9, 29.8, 34.9, 39.7 and 44.2 °C. These amount to an equivalent altitude compared to the inlet of the test section for the hottest run of respectively .75, 1.52, 2.20, 2.86 and 3.485 m, the test section being 3.5 m long. Simulation results are presented in Fig. 8.

The two quantities simulated that are measured in the experiment are void fraction and liquid temperature. One can see that the difference between simulations the local void fraction is around 0.01, and the difference on liquid temperature under 1 °C. However, if we plot the difference between the liquid temperature and the saturation temperature to take into account the effect of pressure on saturation temperature, the difference is reduced to at most .3 °C. These are all well underneath the experimental uncertainty, which justifies the use of the equivalent altitude in the test tube of the DEBORA database where it had the most reasons not to be valid. Therefore, we validate the use of the equivalent altitude for the whole DEBORA experimental database.

## 5. Fixed-diameter boiling flow simulations and discussion

In this section, we take advantage of the validity of the test tube hypothesis on the DEBORA database to run simulations. We transform sets of experimental measures of the Sauter mean diameter for different inlet temperatures into measures at different altitudes. This enables us to run simulations with a 3D map of the experimental diameter. The

physical quantities at each equivalent altitude are compared with quantities measured at the outlet for different inlet temperatures, reducing the total number of calculations required. Finally, this enables us to picture the evolution of the flow along the whole boiling length, and not only its configuration at the outlet for a given inlet temperature.

**Selected runs from the DEBORA database.** We call $j$ a run number in a test tube and $X_j$ the outlet thermodynamic quality of run $j$. $J_{max}$ is the run of a given test tube at the highest inlet temperature $T_{in,j_{max}}$. Using the test tube hypothesis, we launch calculations for each test tube with $T_{in} = T_{in,j_{max}}$ and extract the physical parameters at the different $z_j$ so that:

$$X_{T_{in}=T_{in,j_{max}}} (z_j) = X_j \qquad (22)$$

The conditions of the runs that we select are presented in Table 1. These runs all come from Garnier et al. (2001). The Sauter mean diameter is measured along the radius of the channel for each $T_{in,j}$. We then know the experimental value at any $z_j$. We interpolate the experimental Sauter mean diameter at any $(r, z)$ point in the test section and enforce it in our simulations. We run simulations without having to predict the mean Sauter diameter in the flow.

We simulate the test tubes shown in Table 1 using our baseline closure laws. The mesh is a 2°-wide slice, with 20 radial and 500 axial cells (top left in Fig. 4). The void fraction profiles are presented in Fig. 9, along with the experimental Sauter mean diameters that were enforced in the simulation. The wall-peaked profiles are relatively well predicted for test tubes I and III. However, the experimental void fraction in test tube II is core-peaked and the simulation is far off.

More detailed simulation results, i.e. liquid temperatures, gas velocities and force balances, are presented in Fig. 10 for test tube II-G2P14W16. This test tube was selected for detailed analysis as it is the only one in which all different measures possible in the DEBORA setup were conducted.

**Lift force.** The experimental void fraction in test tube II is core-peaked. Furthermore, for $X_j = 0.1343$ in test tube I and $X_j \ge 0.0479$ in test tube III, the experimental void fraction peak moves away from the wall (this is very slight for test tube I). The only closure law in the two-fluid framework with which we work that can create a center-peaked void fraction profile is the lift force with a negative lift coefficient. This means that in order to simulate the three test tubes studied, a lift force with a negative lift coefficient must be used in the near-wall region.

In all three test tubes, the $f(W_o)$ contribution to the Sugrue (2017) lift coefficient is negative at high outlet qualities (see Eq. (9)). This strengthens the case of the use of a negative lift coefficient. However, the total lift coefficient is damped by the $g(\alpha)$ term and the resulting
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison between simulation performed in experimental configuration and using the equivalent altitude. No experimental data is pictured here. Full lines correspond to physical quantities at the equivalent altitude for the highest inlet temperature of our run, and dashed lines to the value at the outlet for the experimental conditions.</figcaption></figure>

<caption>**Table 1**</caption>

Flow conditions of the DEBORA cases we study. The DEBORA nomenclature for the test tube is given in the first column. $z_j$ is taken at the beginning of the heated length. Each $z_j$ matches a run in Garnier et al. (2001) in which $X_{\text{outlet}} = X(z_j)$.

<table><thead><tr><th>Test tube</th><th>Mass flux kg/(m²s)</th><th>Press. MPa</th><th>Heat flux W/m²</th><th>T<sub>in,max</sub> °C</th><th>z<sub>j</sub> m</th><th>X(z<sub>j</sub>)</th></tr></thead><tbody><tr><td rowspan="5">I G1P30W12</td><td rowspan="5">1007</td><td rowspan="5">3.01</td><td rowspan="5">58.2</td><td rowspan="5">73.7</td><td>0.72</td><td>-0.2165</td></tr><tr><td>1.48</td><td>-0.0973</td></tr><tr><td>2.48</td><td>0.0585</td></tr><tr><td>2.98</td><td>0.1343</td></tr><tr><td>3.485</td><td>0.2173</td></tr><tr><td rowspan="6">II G2P14W16</td><td rowspan="6">2016</td><td rowspan="6">1.458</td><td rowspan="6">76.26</td><td rowspan="6">44.21</td><td>0.75</td><td>-0.0677</td></tr><tr><td>1.52</td><td>-0.0185</td></tr><tr><td>1.81</td><td>0.0014</td></tr><tr><td>2.20</td><td>0.0261</td></tr><tr><td>2.49</td><td>0.046</td></tr><tr><td>2.86</td><td>0.0687</td></tr><tr><td rowspan="6">III G3P26W23</td><td rowspan="6">2994</td><td rowspan="6">2.618</td><td rowspan="6">109.3</td><td rowspan="6">72.49</td><td>1.64</td><td>-0.0519</td></tr><tr><td>2.02</td><td>-0.0177</td></tr><tr><td>2.41</td><td>0.0164</td></tr><tr><td>2.84</td><td>0.0479</td></tr><tr><td>3.14</td><td>0.077</td></tr><tr><td>3.485</td><td>0.1005</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Void fraction results of our simulations on the DEBORA setup using the baseline literature closure laws, and experimental Sauter mean diameters enforced in the simulation. The simulation configurations can be found in Table 1.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Detailed results of our simulations on test tube II-G2P14W16 of the DEBORA setup using the baseline set of closures. Simulated radial and axial vapor velocities were only plotted when $\alpha_v > .02$. The force balances shown are divided by $\alpha_v\rho_v$ so that regions with different void fractions can be compared with ease.</figcaption></figure>
lift coefficient is near-zero, resulting in wall-peaked simulated profiles in all test tubes (see Fig. 10).

**Turbulent dispersion force.** The turbulent dispersion force pilots the void fraction gradient in the simulations, i.e. the slope of the void fraction profiles. For $\alpha_v < .2$, the slopes of the simulated profiles are coherent with those of the experiments. However, this is no longer the case for $X_j = 0.1343$ in test tube I: the void peak at the wall is much larger and the overall slope is much steeper than in the experiment. This suggests that the turbulent dispersion could be higher than that of Burns at high void fractions. The $X_j = 0.2173$ profile in test tube I is difficult to interpret as the void fraction peak at the wall is reminiscent of the beginning of a boiling crisis. The simulated void fraction profiles in test tube II are too different from the experimental ones due to the issues with the lift force to infer anything on turbulent dispersion.

**Heat transfer and temperature profile.** The wall-peaked void fraction profiles are relatively well predicted for test tubes I and III, suggesting a globally correct condensation prediction in these flows. The temperature predictions for test tube II (see Fig. 10) are lower than the experiment for $X_j = -0.0102$. This is coherent with an overestimation of the void fraction for $X_j = -0.0185$.

**Gas velocities.** In test tube II, the axial vapor velocities are not core-peaked like the experiment because the vapor profile is wall-peaked, leading to a buoyancy effect at the wall (Fig. 10 top center). Axial gas velocities are significantly underestimated, especially at high void fraction for $X_j \ge 0.0775$. At these qualities, the liquid is at saturation temperature. As the total mass flow rate is a conserved quantity, the fluid as a whole cannot accelerate. Therefore the discrepancy can only be explained by a too small relative velocity, i.e. a too large drag coefficient in the simulation.

Vapor radial velocities are negative (Fig. 10 bottom center), indicating that vapor is produced at the wall and migrates towards the center of the pipe.

**Momentum balance.** The momentum balances are at equilibrium because the flow is stationary. The axial one on the vapor phase (top right plot in Fig. 10) shows that the pressure gradient is mainly compensated by the drag force, as expected. Two other terms play a role in the near-wall region. Wall friction at low void fraction is only applied to the liquid phase, and is calculated through the diffusion term in the code. However, when the near-wall void fraction in a cell exceeds 0.5, as is the case here, wall friction is also applied to the vapor for numerical stability which is why it is non-negligible. Convection plays a role since vapor is accelerated as it moves from the near-wall to the core region. This is discussed in detail in Section 6.1. In TrioCFD, the convection and virtual mass terms are coded together, so we cannot separate the numerical contributions.

The radial force balance (bottom right plot in Fig. 10) shows that the drag force is compensated by the turbulent dispersion. The lift force is non-zero only where the void fraction is small enough for the $\alpha_v$-dependent term in the Sugrue (2017) lift to be non-zero (see Eq. (9)). This force balance is very different from a developed adiabatic flow, where there is no radial gas velocity. The turbulent dispersion force is then at equilibrium with the lift force and the wall correction (Marfaing et al., 2016).

## 6. Model improvement

In the previous section, we have shown that interfacial force closure laws validated in atmospheric-pressure adiabatic flows are no longer valid in pressurized water reactor similarity conditions.

This is not necessarily surprising, as most of these models are used outside of their validity domain. Pressurized water reactor flows are extremely turbulent, have low surface tensions and density ratios. When

boiling in some accidental conditions they have high void fractions and can resemble emulsions (Hosler, 1967; François et al., 2011).

In this section, we propose a new set of closures for these flows, using a pragmatic physics-base approach that takes into account the particularities of pressurized water reactor flows.

### 6.1. Relative velocity inversion

Fig. 11 presents the impact of the virtual mass force on the gas-liquid relative velocity on test tube III. The first column contains the results for the baseline case, with $C_{VM}$ given by Eq. (7). The axial relative velocity is positive at the wall, becomes negative as the bubbles move away from the wall and becomes positive again at the core of the flow. This is due to bubble inertia in the simulation: the velocity of bubbles in the near-wall cell is smaller than the velocity of the liquid in the core. As bubbles migrate towards the pipe center, their inertia prevents them from accelerating immediately, hence the negative relative velocity. This effect is exacerbated if vapor is produced with zero velocity, i.e. if a momentum source $-\Gamma_{w,l \to v} \vec{u}_v$ is added to the vapor phase in the near-wall cell (center column in Fig. 11). This behavior has direct effects on many interfacial terms: it cancels out or inverts the sign of the lift force, and reduces the turbulent dispersion force and condensation. There are very few experiments where we have access to gas and liquid velocities in identical configurations (Roy et al., 2002; François et al., 2021). To the best of our knowledge, the gas velocities were never measured lower than the liquid ones underneath the experimental error. We therefore believe the relative velocity inversion to be non-physical and strive to prevent it in our simulations.

We therefore run test tube III without the virtual mass term to reduce the inertia of the vapor (right column in Fig. 11). The relative velocity remains positive, and all further simulations are carried out without virtual mass.

### 6.2. Deformed bubble hypothesis

In order to evaluate if bubbles are deformed or not in the flow, we plot two Weber numbers (see Fig. 12). The first, $We$, is defined using the relative velocity and accounts for the effect of the drag force on bubble deformation. The second, $We_c$, represents the impact of turbulent velocity fluctuations at the scale of the bubbles:

$$ 
\begin{align} 
We &= \frac{\rho_l d_b ||\vec{u}_e - \vec{u}_l||^2}{\sigma} \tag{23} \\ 
We_c &= \frac{\rho_l d_b (ed_b)^{2/3}}{\sigma} \nonumber 
\end{align} 
$$

In this computation $d_b$ comes from the experimental measures. The other physical quantities come from the simulations carried using the baseline set of closures, but without the virtual mass force, and enforcing the experimental diameter as in Section 5.

The Weber number $We$ is smaller than 1 in the near-wall region, which means that the bubbles should be spherical according to this criterion. In the core, it varies between 2 and 3. For such Weber numbers, bubbles should be deformed (Wallis, 1974). The turbulence-driven Weber number $We_c$ has an opposite behavior: it is large in the near-wall region and diminishes in the core. Bubble experience significant deformation when $We_c \ge 2$ (Hinze, 1955; Risso and Fabre, 1998; Masuk et al., 2021). Therefore, on all three test tube bubbles are in a region where they are deformed due to their rise or due to turbulence.

Furthermore, as the void fraction increases, for a same average bubble size, the distance between bubbles decreases. If we assume that bubbles are spherical and monodisperse, then the most dense possible configuration is face-centered cubic. The relationship between the void
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Void fraction and relative axial velocity results for test tube III-G3P26W23, for three different configurations. **Left:** with virtual mass and vapor formation at first-cell velocity. **Center:** with virtual mass and vapor formation at zero velocity. **Right:** without virtual mass and vapor formation at first-cell velocity.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Weber numbers from DEBORA experiment simulations using the experimental diameters as input but without the virtual mass force. See Eq. (23) for definitions, Fig. 9 for void fraction comparison with the experiment and Table 1 for simulation configurations. **Top row:** Weber numbers defined using the relative velocities. **Middle row:** Weber numbers defined using the turbulent velocity fluctuations at the scale of the bubbles.</figcaption></figure>
fraction, bubble diameter and distance between bubbles in then:

$$ \frac{\frac{\pi}{6} d_b^{3/6}}{\frac{\sqrt{2}}{12} (d_b + d_{dist})^3} = \alpha_v \qquad (24) $$

The distance between two bubbles $d_{dist}$ can then be expressed as a function of the void fraction:

$$ \frac{d_{dist}}{d_b} = \sqrt[3]{\frac{\pi}{3\sqrt{2}} \frac{1}{\alpha_v} - 1} \qquad (25) $$

The dimensionless distance between bubbles is plotted as a function of $\alpha_v$ in Fig. 13. $d_{dist} \approx d_b$ for $\alpha_v = 0.1$, and $d_{dist} \approx 1/2d_b$ for $\alpha_v = 0.2$. For such void fractions the bubbles are almost touching each other. $\alpha_v > 0.1$ in almost all of the DEBORA database (see Fig. 9). In flow boiling, bubble agitation and turbulence mean that they necessarily collide, which will also lead to bubble deformation.

<figure><img src="image_14.png" /><figcaption>Fig. 13. Distance between two spherical bubbles in a compact packing as a function of the void fraction.</figcaption></figure>

We call deformed bubble hypothesis the hypothesis that in flow boiling in PWR conditions, bubbles are deformed, i.e. non-spherical, for $\alpha_v > 0.1$ by the combined effects of turbulence, drag and crowding. As this void fraction is easily achieved, we will consider in the rest of this work that bubbles are always deformed in pressurized water reactor flows.

Finally, the crowding that occurs for $\alpha_v > 0.1$ will have impacts on interfacial force and heat transfer correlations compared with the most commonly used forms based on measures on single bubbles, though we lack the means to quantify them.

### 6.3. Interfacial forces

**Drag force.** Bubble rise velocities increases with bubble diameter until the later becomes similar to the capillary length. It then becomes independent of the diameter as the bubbles deform (Clift et al., 2013). The drag force formulations most commonly used in the literature take this effect into account by having a drag coefficient that depends on the bubble Reynolds number and diameter. As shown by Sugrue (2017), when the bubbles are sufficiently deformed the gas-water relative velocity becomes independent of the bubble diameter for all historical formulations (Ishii and Zuber, 1979; Tomiyama et al., 1998; Bozzano and Dente, 2001). The relative velocity for deformed bubbles is given by Ishii and Zuber (1979):

$$ u_{r,IZdef} = \|\bar{u}_g - \bar{u}_l\| = \sqrt{2} \left( \frac{g\sigma(\rho_l - \rho_g)}{\rho_l^2} \right)^{1/4} \qquad (26) $$

Fig. 14 shows the rise velocity for Tomiyama et al. (1998) drag force, the deformed (Ishii and Zuber, 1979) formulation, and a modification of the latter: $\min(1, d_b/L_c) \cdot u_{r,IZdef}$. All 3 formulations are very similar for large diameters and the contaminated (Tomiyama et al., 1998) relative velocity can be approached by:

$$ u_{r,Tcont} \approx \min(1, d_b/L_c) \cdot u_{r,IZdef} \qquad (27) $$

Ishii and Zuber (1979) also proposed a drag force formulation for deformed bubbles, that yields the bubble rise velocity given in Eq. (26). We use this formulation for our set of closures, as we consider that

<figure><img src="image_15.png" /><figcaption>Fig. 14. Bubble rise velocity for different drag formulations as a function of the dimensionless bubble diameter for 26 bar R12. All physical properties are taken at saturation temperature.</figcaption></figure>

bubbles are systematically deformed in nuclear reactor conditions:

$$ \begin{aligned} \vec{F}_{drag} &= -\frac{3}{4} C_D \frac{\alpha_v \rho_l}{d_b} \|\bar{u}_g - \bar{u}_l\| (\bar{u}_g - \bar{u}_l) \\ C_D &= \frac{2}{3} \frac{d_b}{L_c}, \quad L_c = \sqrt{\frac{\sigma}{g(\rho_l - \rho_g)}} \end{aligned} \qquad (28) $$

**Lift force.** The lift coefficient is positive for shear flow around a spherical bubble, and becomes negative for deformed bubbles (Legendre and Magnaudet, 1998; Tomiyama et al., 2002; Sugrue, 2017). Furthermore, as discussed in Section 5, the void fraction profiles in the DEBORA experiment can only be explained by the presence of a negative lift coefficient in the two-fluid model. In Fig. 9, the lift force seems to have an impact for $\alpha_v \gtrsim 0.3$. We therefore build a lift coefficient that is dependent on the void fraction. This approach is similar to the one of Yoon et al. (2017). It is equal to 0 up to $\alpha_v = 0.25$, as we have no proof of the impact of lift on the flow in this region. We then propose the following evolution for $C_L$. The coefficient drops linearly to a value of -0.2, in-between the minimal values proposed by Tomiyama et al. (2002) (-0.25) and Sugrue (2017) (-0.15). It returns to 0 for very high void fractions for numerical stability. The $C_L(\alpha_v)$ function is plotted in Fig. 15, and the coefficient is worth:

$$ C_L = \begin{cases} 0 & \text{if } \alpha_v < 0.25 \\ \max(-0.2, -0.7 \cdot (\alpha_v - 0.25)) & \text{if } 0.25 \le \alpha_v < 0.7 \\ \alpha_v - 0.9 & \text{if } 0.7 \le \alpha_v < 0.9 \\ 0 & \text{if } 0.9 \le \alpha_v \end{cases} \qquad (29) $$

<figure><img src="image_16.png" /><figcaption>Fig. 15. The lift coefficient proposed in this work as a function of the void fraction.</figcaption></figure>

**Turbulent dispersion force.** We have shown that the Burns et al. (2004) turbulent dispersion force in low-void fraction pressurized water reactor flows treats bubbles like a passive scalar diffused by turbulence (Reiss, 2024). For high void fractions, we believe that packing and bubble collisions will have an effect on the bubble dispersion. However, we lack experimental data to adjust such a model. We therefore continue to use the Burns et al. (2004) formulation, but our work enables us to better interpret the effect of the force on predicted void fraction distributions.

The Burns et al. (2004) turbulent dispersion coefficient contains a term in $C_D/d_b$ (see Eq. (10)). Furthermore, we use the Ishii and Zuber
(1979) deformed bubble drag force (see Eq. (28)). $d_b$ in these two terms cancel out and $C_{TD}$ is independent of the bubble diameter.

## 6.4. Heat transfer models

**Heat flux partition.** Heat flux partitions that improve predicted wall temperatures compared with the Kurul and Podowski (1990) model have been proposed in recent years (Basu et al., 2005; Kommajosyula, 2020; Favre, 2023). However, there are still issues with these models. They require many intermediate quantities to be closed, like the bubble wait time, departure frequency and diameter. Experimental data is still lacking, particularly at high pressure and mass flow rates. This myriad of intermediate quantities makes the models very difficult to read and the outputs hard to interpret.

We have therefore proposed a novel heat flux partition based on a CFD-scale onset of significant void (OSV) correlation (Reiss et al., 2024). It is shown that at the OSV, if the near-wall cell has a size $y_1 \ge 100$, the heat transfer to the liquid phase is:

$$q_{l,OSV} = H_{l,OSV}(y_{+,1})(T_s - T_l(y_1)) = \frac{\rho C_{pl} u_{\tau}}{2.12 \log(y_{+,1}) - 7} (T_s - T_l(y_1)) \quad (30)$$

This is used to create a heat flux partitioning algorithm. When the wall temperature is known, the inputs are $y_1$, $u_{\tau}$, $T_w$, $T_l$, $T_l(y_1)$ and the physical properties of the liquid. $y_1$ is the size of the first element. The steps of the algorithm are the following:

1. Calculate single-phase heat flux $q_{SP}$ using the Kader (1981) heat transfer coefficient (Eq. (3)):

2. Calculate total boiling heat flux $q_{Boil}$ using a total heat flux correlation of one's choosing (Jens and Lottes, 1951; Thom et al., 1965; Frost and Dzakowic, 1967 where the pressure $P$ must be in bar):

$$ 
\begin{aligned}
q_{w,Jens&Lottes} &= \left( \frac{T_w - T_s}{25} \exp(P/62) \right)^4 \\
q_{w,Thom&et al.} &= \left( \frac{T_w - T_s}{22.65} \exp(P/87) \right)^2 \\
q_{w,Frost&Dzakowic} &= \frac{\lambda_{Ts}(h_{gs} - h_{Ts})\rho_v}{8\sigma T_s} \left( \frac{T_w - T_s}{Pr_s} \right)^2
\end{aligned}
\quad (31) $$

3. If $T_l(y_1) \ge T_s$; we are in saturated boiling, therefore all of the energy is used for evaporation:

$$ 
\left\{ 
\begin{aligned}
q_w &= q_{Boil} \\
q_{wl} &= 0 \\
q_{w,l \to v} &= q_{Boil}
\end{aligned}
\right. 
\quad (32) $$

4. Else if $q_{SP} \ge q_{Boil}$; we are in a single-phase regime before the onset of nucleate boiling: $q_{SP} = q_w$ goes into the liquid phase:

$$ 
\left\{ 
\begin{aligned}
q_w &= q_{sp} \\
q_{wl} &= q_{sp} \\
q_{w,l \to v} &= 0
\end{aligned}
\right. 
\quad (33) $$

5. Else we are between the onset of nucleate boiling and saturated boiling. The total heat flux is then $q_w = q_{Boil}$. To determine if we have passed the onset of significant void, we calculate $H_{l,OSV}$ (Eq. (30)):

(a) If $q_{l,OSV} = H_{l,OSV}(T_s - T_l(y_1)) > q_w = q_{Boil}$, we have not yet reached the onset of significant void: all of the energy goes into the liquid phase:

$$ 
\left\{ 
\begin{aligned}
q_w &= q_{Boil} \\
q_{wl} &= q_{Boil} \\
q_{w,l \to v} &= 0
\end{aligned}
\right. 
\quad (34) $$

(b) Else we consider that the heat transfer to the liquid phase is the same than at OSV (Eq. (30)), and that the rest of the energy is used for evaporation:

$$ 
\left\{ 
\begin{aligned}
q_w &= q_{Boil} \\
q_{wl} &= H_{l,OSV}(T_s - T_l(y_1)) \\
q_{w,l \to v} &= q_{Boil} - H_{l,OSV}(T_s - T_l(y_1))
\end{aligned}
\right. 
\quad (35) $$

As we always have $q_{SP} < q_{l,OSV}$, a condensed way to write steps 3 to 5 of this algorithm is the following:

$$ 
\begin{aligned}
q_{l,OSV} &= \max \left( 0, \frac{\rho_l C_{pl} u_{\tau} (T_s - T_l(y_1))}{2.12 \log(y_{+,1}) - 7} \right) \\
q_w &= \max(q_{sp}, q_{Boil}) \\
q_{wl} &= \min(q_w, q_{l,OSV}) \\
q_{w,l \to v} &= q_w - q_{wl}
\end{aligned}
\quad (36) $$

Compared with (Kurul and Podowski, 1990)-base heat flux partitions, this methodology saves computation time and is easier to implement in a code. From a physics standpoint, it guarantees a better calculation of the total heat flux than mechanistic models, as correlations directly fitted on experimental data are more precise. Furthermore, given the simplicity of the model it is easy to interpret the outputs and the physical mechanisms at play.

**Condensation.** In high-void fraction flow boiling, bubbles are polydisperse, bump into each other, coalesce and break. These micro-scale mechanisms, which are not understood well enough to be modeled, impact the condensation and interfacial heat transfer terms.

None of the other closures that we have proposed in this section require bubble diameter modeling through an interfacial area transport equation (Yao and Morel, 2004) or population balance model (Kreper et al., 2008). We therefore aim to build a correlation that does not require bubble diameter modeling. In the DEBORA database, the measured bubble diameters were mostly close to the capillary length $L_c = \sqrt{\frac{\sigma}{g(p_l - p_v)}}$, and all of them were between 0.25$L_c$ and $2L_c$. For test tubes I, II and III, $L_c$ was 0.4, 0.65 and 0.45 mm (compare with Fig. 10 bottom line). In all three test tubes, $d_b < L_c$ if $T_s - T_l > 1$ °C, i.e. in subcooled regions where condensation is non-negligible bubbles are smaller than $L_c$. We choose the capillary length as the bubble diameter length scale used in the interfacial area ($a_l \sim 6a_v/L_c$ in Eq. (13)). This approach is used for bubbly-flow condensation in the system-scale codes RELAP-7 (Berry et al., 2018) and TRACE-5 (NRC, 2010), widely used in the nuclear industry. When the distance between bubbles (Eq. (25)) is larger than the capillary length, we use the former as the characteristic temperature diffusion length ($d_b$ in (13)). When it is smaller, we use the latter. This yields:

$$ q_{ki} = \alpha_v N u_{fit} \frac{6\lambda_l}{L_c^2} \min \left( 1, \sqrt{\frac{1}{3\sqrt{2} \min(\alpha_v, 0.6)}} - 1 \right) (T_v - T_l) \quad (37) $$

With $Nu_{fit}$ a fitted constant bubble Nusselt number. $Nu_{fit} = 30$ was found to be an optimal value on the DEBORA database. $\frac{\pi}{3\sqrt{2}} \approx 0.74$. As the cubic root term in Eq. (37) goes to 1 as $\alpha_v \to \frac{\pi}{3\sqrt{2}}$, we included $\min(\alpha_v, 0.6)$ for numerical stability. When $\alpha_v = 0.6$, the heat transfer coefficient is multiplied by ~10 compared with low-volume fraction cases and in practice the liquid is at saturation temperature.

This correlation is not equivalent to using the capillary length as bubble diameter and using a Nusselt correlation from the literature. We never fit the intermediate quantity that is the bubble diameter, but only the total interfacial heat flux. We can have a correct interfacial heat flux with bubble diameters significantly larger or smaller than $L_c$.

We compare this expression to classical correlations from the literature: Ranz and Marshall (1952), Chen and Mayinger (1992), Zeitoun et al. (1995) and Kim and Park (2011). As discussed in Section 2.4, all of them are based on Eq. (13). The specificity of each correlation then
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Comparison between different correlations used for interfacial condensation (Ranz and Marshall, 1952, Chen and Mayinger, 1992, Zeitoun et al., 1995, Kim and Park, 2011 and current work). The relative velocity between the phases was chosen so that the contaminated (Tomiyama et al., 1998) drag force compensates buoyancy. A 5 °C subcooling and 0.1 void fraction were used for the Zeitoun et al. (1995) and Kim and Park (2011) formulations. **Top row:** Nusselt number as a function of the bubble diameter nondimensionalized by the capillary length. **Bottom row:** Heat transfer coefficient over (Ranz and Marshall, 1952) heat transfer coefficient (chosen as reference) for different correlations. **Left:** Results in the conditions of test tube II-G2P14W16, i.e. boiling water reactor similarity. **Center:** Results in the conditions of test tube III-G3P26W23, i.e. pressurized water reactor similarity. **Right:** Results for water at 100 bar.</figcaption></figure>

lies in the Nusselt number formulation, that can depend on the Prandtl, bubble Reynolds or Jacob numbers:

$$ 
\begin{aligned}
Nu_{\text{Ranz\&Marshall}} &= 2 + 0.6 Re_b^{0.5} Pr^{0.33} \\
Nu_{\text{Chen\&Mayinger}} &= 0.185 Re_b^{0.7} Pr^{0.5} \\
Nu_{\text{Zeitoun et al.}} &= 2.04 Re_b^{0.61} \alpha^{0.328} Ja^{-0.308} \\
Nu_{\text{Kim\&Park}} &= 0.2575 Re_b^{0.7} Ja^{-0.2043} Pr^{-0.4564}
\end{aligned}
\qquad (38) $$

Where $Ja = \frac{\rho_l C_{pl} ||T_{sat}-T_l||}{\rho_g h_{lg}}$ is the Jacob number.

All of these correlations contain the bubble Reynolds number, i.e. the bubble diameter and relative velocity. If the Tomiyama et al. (1998) drag formulation is used in a simulation, $u_{gz} - u_{lz} \propto d_b$ when $d_b \lesssim L_c$ and $u_{gz} - u_{lz}$ is independent of $d_b$ when $d_b \gtrsim L_c$ (see Eq. (27)). Therefore, $Re_b \propto d_b^2$ when $d_b \lesssim L_c$ and $Re_b \propto d_b$ when $d_b \gtrsim L_c$. Table 2 compares the dependence of the interfacial heat flux on the bubble diameter for different correlations. The dependence on the bubble diameter remains, but is much less significant than $q_{kl} \propto d_b^{-2}$ that is visible in Eq. (13), particular in subcooled regions where $d_b < L_c$.

<caption>**Table 2**</caption>

Dependence of the interfacial heat transfer on the diameter for various correlations from the literature and for the current work.

<table><thead><tr><th>Reference</th><th>d<sub>b</sub> &le; L<sub>c</sub></th><th>d<sub>b</sub> &ge; L<sub>c</sub></th></tr></thead><tbody><tr><td>Ranz and Marshall (1952)</td><td>q<sub>ki</sub> &prop; d<sub>b</sub><sup>-1</sup></td><td>q<sub>ki</sub> &prop; d<sub>b</sub><sup>-1.5</sup></td></tr><tr><td>Chen and Mayinger (1992)</td><td>q<sub>ki</sub> &prop; d<sub>b</sub><sup>-0.6</sup></td><td>q<sub>ki</sub> &prop; d<sub>b</sub><sup>-1.3</sup></td></tr><tr><td>Zeitoun et al. (1995)</td><td>q<sub>ki</sub> &prop; d<sub>b</sub><sup>-0.78</sup></td><td>q<sub>ki</sub> &prop; d<sub>b</sub><sup>-1.39</sup></td></tr><tr><td>Kim and Park (2011)</td><td>q<sub>ki</sub> &prop; d<sub>b</sub><sup>-0.6</sup></td><td>q<sub>ki</sub> &prop; d<sub>b</sub><sup>-1.3</sup></td></tr><tr><td>Current work</td><td>q<sub>ki</sub> &prop; d<sub>b</sub><sup>0</sup></td><td>q<sub>ki</sub> &prop; d<sub>b</sub><sup>0</sup></td></tr></tbody></table>

The Nusselt as a function of the bubble diameter for different conditions is presented in the top row of Fig. 16. The inflection in the plots from the literature at $d_b \approx L_c$ come from the relative velocity becoming constant from this point. The bottom row contains

the heat transfer coefficient compared to that of the Ranz and Marshall (1952) correlation. Correlations from the literature yield very different Nusselt numbers and heat transfer coefficients: the Zeitoun et al. (1995) expression can be five times higher than that of Ranz and Marshall (1952) or (Kim and Park, 2011), with the Chen and Mayinger (1992) correlation in the middle of the ballpark. The relative heat transfer of our proposal compared with that of Ranz and Marshall (1952) has an inflection at $d_b \approx L_c$, again as the relative velocity plateaus. In all conditions, our proposed expression is in the inter-model uncertainty range.

## 6.5. Comparison of the proposed closure terms with the baseline model

Table 3 compares the different closure terms of the baseline model with our proposal. The expressions of all of the terms in the current proposal are more simple than those of the baseline set of closures. Furthermore, no bubble diameter modeling by an interfacial area transport equation or population balance model is required. Fig. C.23 in Appendix C is a diagram of the links between different unknowns and equations used in classical formulations and the ones that we no longer have in our simplified methodology.

This approach makes the numerical implementation of such a set of closures easier. More importantly, it significantly simplifies the interpretation of the results from a multiphase CFD simulation. Each individual term is easier to read and understand. Furthermore, there is less retroaction between terms than in the baseline set of closures. This eases the understanding the root cause of discrepancies between experimental and simulated results, reduces the risk of error compensation and will enable us to iteratively improve this set of models.

Having no explicit dependency on the bubble diameter may seem surprising for a CFD-scale code. To the best of our knowledge, no such approach has been proposed for a boiling flow. However, it is common in subchannel and system codes in the nuclear industry:
<table><caption>**Table 3**
Comparison of the closures from the baseline set of models with our proposal. TD: turbulent dispersion. VM: virtual mass. HFP: heat flux partition. Cond: condensation.</caption><thead><tr><th>Term</th><th>Baseline closure</th><th>Current proposal</th></tr></thead><tbody><tr><td>Drag</td><td>Tomiyama et al. (1998)<br>C<sub>D</sub> = max (<sup>24</sup>/<sub>Re<sub>b</sub></sub>(1 + .15Re<sub>b</sub><sup>687</sup>), <sup>8Eo</sup>/<sub>3Eo + 12</sub>)</td><td>Ishii and Zuber (1979)<br>deformed bubble<br>C<sub>D</sub> = <sup>2</sup>/<sub>3</sub>d<sub>b</sub>√<frac><sup>g(ρ<sub>l</sub> - ρ<sub>v</sub>)</sup><sub>σ</sub></frac></td></tr><tr><td>Lift</td><td>Sugrue (2017) C<sub>L</sub> = f(Wo)·g(α<sub>v</sub>)</td><td>Current work</td></tr><tr><td>TD</td><td>Burns et al. (2004)<br>C<sub>TD</sub> = <sup>3</sup>/<sub>4</sub>C<sub>D</sub>|ũ<sub>g</sub> - ũ<sub>l</sub>|<sup>1</sup>/<sub>ω</sub>(1 + <sup>α<sub>v</sub></sup>/<sub>α<sub>l</sub></sub>)</td><td>Burns et al. (2004)<br>C<sub>TD</sub> = <sup>3</sup>/<sub>4</sub>C<sub>D</sub>|ũ<sub>g</sub> - ũ<sub>l</sub>|<sup>1</sup>/<sub>ω</sub>(1 + <sup>α<sub>v</sub></sup>/<sub>α<sub>l</sub></sub>)</td></tr><tr><td>VM</td><td>Zuber (1964)<br>adapted in current work<br>C<sub>VM</sub> = min(<sup>1</sup>/<sub>2</sub>α<sub>g</sub>, <sup>1</sup>/<sub>2</sub>α<sub>l</sub>)</td><td>C<sub>VM</sub> = 0</td></tr><tr><td>HFP</td><td>Kurul and Podowski (1990)<br>q<sub>SP</sub> = (T<sub>w</sub> - T<sub>l</sub>(y))<sup>ρ<sub>l</sub>C<sub>pi</sub>u<sub>τ</sub></sup>/<sub>ρ<sub>v</sub>u<sub>v</sub>(y<sub>+</sub>)</sub><br>N<sub>s</sub> = (210(T<sub>w</sub> - T<sub>l</sub>))<sup>1.8</sup><br>d<sub>b, det</sub> = 10<sup>-4</sup>·(T<sub>w</sub> - T<sub>s</sub>) + 0.0014<br>f<sub>dep</sub> = √(<frac><sup>4</sup>/<sub>3</sub>g(ρ<sub>l</sub> - ρ<sub>g</sub>)</sup>/<sub>ρ<sub>l</sub>d<sub>b, det</sub></sub>)<br>A<sub>b</sub> = min(1, π/4·N<sub>s</sub>d<sub>b, det</sub><sup>2</sup>)<br>q<sub>c</sub> = (1 - A<sub>b</sub>)q<sub>SP</sub><br>q<sub>q</sub> = 2A<sub>b</sub>λ<sub>l</sub>(T<sub>w</sub> - T<sub>l</sub>)√<frac><sup>f<sub>dep</sub>ρ<sub>l</sub>Cρ<sub>l</sub></sup>/<sub>πλ<sub>l</sub></sub></frac><br>q<sub>w,c</sub> = q<sub>c</sub> + q<sub>q</sub><br>q<sub>w,c,l→v</sub> = <sup>π</sup>/<sub>6</sub>f<sub>dep</sub>d<sub>b, det</sub><sup>3</sup>ρ<sub>g</sub>h<sub>lg</sub>N<sub>s</sub><br>q<sub>w</sub> = q<sub>c</sub> + q<sub>q</sub> + q<sub>w,c,l→v</sub></td><td>Reiss et al. (2024)<br>q<sub>SP</sub> = (T<sub>w</sub> - T<sub>l</sub>(y))<sup>ρ<sub>l</sub>C<sub>pi</sub>u<sub>τ</sub></sup>/<sub>ρ<sub>v</sub>u<sub>v</sub>(y<sub>+</sub>)</sub><br>q<sub>Boil</sub> = (<frac><sup>T<sub>w</sub> - T<sub>l</sub></sup>/<sub>22.65</sub>exp(P/87))<sup>2</sup><br>q<sub>l,OSV</sub> = max(0, <frac><sup>ρ<sub>l</sub>C<sub>pi</sub>u<sub>τ</sub>(T<sub>s</sub> - T<sub>l</sub>(y))</sup>/<sub>2.12 log(y<sub>+</sub>) - 7</sub>)<br>q<sub>w</sub> = max(q<sub>SP</sub>, q<sub>Boil</sub>)<br>q<sub>l</sub> = min(q<sub>w</sub>, q<sub>l,OSV</sub>)<br>q<sub>w,l→v</sub> = q<sub>w</sub> - q<sub>l</sub></td></tr><tr><td>Cond</td><td>Ranz and Marshall (1952)<br>q<sub>ki</sub> = <frac><sup>6α<sub>v</sub>λ<sub>l</sub>(T<sub>g</sub> - T<sub>l</sub>)</sup>/<sub>d<sub>b</sub><sup>2</sup></sub>(2 + 0.6Re<sub>b</sub><sup>1/2</sup>P<sub>r</sub><sup>1/3</sup>)</td><td>Current work<br>q<sub>ki</sub> = <frac><sup>6α<sub>v</sub>Nu<sub>fl</sub>λ<sub>l</sub>g(ρ<sub>l</sub> - ρ<sub>v</sub>)(T<sub>g</sub> - T<sub>l</sub>)</sup>/<sub>σ min(1, √<frac><sup>π</sup>/<sub>3√2</sub><frac><sup>1</sup>/<sub>min(a<sub>v</sub>, 0.6)</sub></frac></frac>) - 1</sub></frac></td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Void fraction simulation results for the set of closures proposed in this paper.</figcaption></figure>

RELAP (Berry et al., 2018), TRACE (NRC, 2010), CTF (Salko et al., 2023) or CATHARE-3 (Emonot et al., 2011) can be used to simulate boiling flows and do not use an interfacial area transport equation in this regime.

## 6.6. Simulations with new set of closures

We simulate the test tubes presented in Table 1 with the set of closures proposed in this paper. Void fraction results are shown in Fig. 17, and detailed results are shown for test tube II-G2P14W16 in
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Detailed results of our simulations on test tube II-G2P14W16 of the DEBORA setup using our proposed set of closures. Simulated radial and axial vapor velocities were only plotted when $\alpha_v > .02$. The force balances shown are divided by $\alpha_v \rho_v$ so that regions with different void fractions can be compared with ease.</figcaption></figure>

Fig. 18. The simulation predictions are improved compared to the baseline model for test tubes II and III, and are of a similar quality for test tube I, despite not using the experimental diameter to close the system of equations. In particular, the void fraction predictions are improved in the near-wall region. This region is critical for the prediction of the critical heat flux with CFD codes (Mimouni et al., 2016). Furthermore, the liquid temperature and gas velocity predictions for test tube II are significantly improved (Fig. 18 lower left). The lift force is significantly larger in the radial force balance than in the baseline closure (Fig. 18 lower right). It pulls bubbles towards the center of the pipe, and is opposed to the drag and turbulent dispersion forces. The radial forces are much stronger than in the baseline set of closures (compare with Fig. 10 lower right).

Pressurized water reactor subchannels can be seen as the assembly of 3 different elementary shapes (see Fig. 19): a tube, in the center of a subchannel; a 3 mm-wide channel, between two rods; and an annulus, around a rod. We have partially validated our new set of closures in the tube geometry of the DEBORA database. However, validation is required for our set of closures for the two other elementary geometries. Chu et al. (2017) recently presented experimental results in a R134A annular channel, in which they were able to change the altitude of measuring plane. Void fractions, gas velocities and bubble diameters were measured using optical probes. Martin (1972) performed X-ray attenuation measures in a water-filled boiling 2.8-mm wide channel. He measured the average void fraction along the width of the channel for different inlet temperatures and conditions. We present void fraction simulation results for these experiments in Fig. 20. The test conditions are given in the legend of the figure.

The shape of the predicted void fraction profiles for all tests are coherent with experimental results, in geometries very different from those in which the closures that we proposed were fitted. This validates our approach of building a simple, high-pressure set of closure models.

<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Diagram of a pressurized water reactor subchannel (fuel rods are in gray) and dimensions of the different experiments used for validation of CFD codes.</figcaption></figure>

The void fraction is underestimated in both (Chu et al., 2017) tests, indicating a possible overestimation of condensation. In the Martin test, the void fraction is slightly overestimated. This is not necessarily surprising. The Nusselt number of condensation formulations from the literature depend on fluid characteristics through the Reynolds, Prandtl and Jacob numbers. Moreover, their characteristic Nusselt number are larger for water than for refrigerant fluids (see top row in Fig. 16). Nu = 30 could be too low for water. Work remains to be done on adapting our condensation correlation across different fluids by extending our flow boiling database to more high-pressure water conditions.
<figure><img src="image_21.png" /><figcaption>**Fig. 20.** Comparison of simulation results for the set of closures proposed in this work with experiments from the literature. Configurations are given in the legend of the figures. **Left and center:** experiments in a 9.5 mm inner diameter 27.2 mm outer diameter annulus from Chu et al. (2017) using R134A as a fluid. $r_+ = 0$ on the heated wall on the inside of the annulus, and $r_+ = 1$ on the opposing wall. **Right:** experiments in a 2.8 mm-channel from Martin (1972) using high-pressure water. $r_+ = 0$ at the center of the channel, and $r_+ = 1$ on the heated wall. There are two full lines for experimental results in each condition as measures were carried out on both sides of the channel. The difference between the two gives a measure of the experimental error.</figcaption></figure>

# 7. Conclusions

We developed a bubbly-flow CFD framework and selected a set of closure terms that we validated on adiabatic atmospheric-pressure data. We show that the test tube hypothesis, i.e. that changing the experimental inlet temperature is similar to moving up or down a measurement section in a flow, is valid in the conditions of the DEBORA experiment and nuclear reactors. This enables us to enforce the experimental diameter in our simulations. We then show that, independently of interfacial area modeling, this selection of models is not appropriate for the nuclear reactor-similarity conditions of the DEBORA experiment. We propose a modified set of closures, based on the hypothesis that bubbles are deformed in PWR conditions. It does not require bubble diameter modeling and improves simulation results on the DEBORA database.

To improve this framework, future plans include conducting a bayesian calibration on the Nusselt number and lift coefficient that we defined (Leoni et al., 2024), increasing the size of the high-pressure validation database and evaluating the model in other flow conditions, in particular atmospheric-pressure flows.

## CRediT authorship contribution statement

**Coretin Reiss:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Formal analysis, Conceptualization. **Antoine Gerschenfeld:** Writing – review & editing, Supervision, Software, Project administration, Methodology, Funding acquisition, Formal analysis, Conceptualization. **Catherine Colin:** Writing – review & editing, Supervision, Project administration, Methodology, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgments

The authors thank Elie Saikali and Yannick Gorsse for their help in developing the multiphase module of TrioCFD, and Alan Burlot and Guillaume Bois for discussions on the physics at play.

## Appendix A. Mesh refinement

Fig. A.21 presents the mesh refinements tests that we conducted.

## Appendix B. Experimental and interpolated bubble diameter profiles

Fig. B.22 presents the experimental and interpolated bubbles diameter profiles used in the adiabatic simulations in Section 3. These diameters were used in the simulation results presented in Fig. 5.

## Appendix C. Links between major terms of a two-fluid model

Fig. C.23 presents the links between the unknowns and the major terms for the proposal of the current work (full boxes and lines) and the standard approach (all boxes and lines) used by Favre et al. (2022), Alatrash et al. (2022), Pham et al. (2023) and Vlček and Sato (2023) for example. The dashed boxes and terms are links that we do not have in our set of closures. We greatly simplify the system of equations solved, but it remains complex.
<figure><img src="image_22.png" /><figcaption>**Fig. A.21.** Mesh refinement verification for the TrioCFD multiphase module. The case simulated is a 3.08 m-long 4 cm-diameter adiabatic pipe (Colin et al., 2012). The liquid and gas superficial velocities are respectively 0.27 m/s and 0.023 m/s. Top row: verification of the equivalence of meshing different angular openings. These meshes are presented in Fig. 4. The 45°, 90° and 360° meshes have 14 radial elements with mesh grading. The 2° mesh have 40 radial elements and the same near-wall cell size. Bottom row: Mesh refinement tests on the 90° slice. The difference between 14 and 28 radial elements is negligible.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. B.22.** Two-phase adiabatic experimental and interpolated bubble diameter distributions in a pipe. Each subfigure matches a subfigure in Fig. 5. *J*<sub>t</sub>: superficial liquid velocity. *J*<sub>g</sub>: superficial gas velocity. Circles: experimental results. Lines: interpolation. A-C. Void fraction results in Hibiki et al. (2001) upwards flow experiments. E. Void fraction results in Colin et al. (2012) downwards flow experiments. F. Void fraction results in Colin et al. (2012) microgravity experiments. In this last case, no experimental diameter profile was measured, only an average bubble diameter for each run.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. C.23.** Comparison of the links between major terms and unknowns for the current work (full boxes and lines) and more standard approaches (all boxes and lines; Favre et al., 2022; Alatrash et al., 2022; Pham et al., 2023; Vlček and Sato, 2023). Each box aligned on the left represents a conservation equation. Each equation has unknowns, specific source terms, and terms that are shared between multiple equations. The latter are given in the box to the right. Each arrow represents an input from an unknown or a generic source term to an equation or a source term. For readability, the physical properties of the fluids are not represented in the figure, but they play a key role in all terms and are functions of the temperature and pressure. grav the gravity, disp the turbulent dispersion force, prod and diss the turbulence production and dissipation terms and coal and frag the bubble coalescence and fragmentation terms.</figcaption></figure>

## References

Abbrecht, P.H., Churchill, S.W., 1960. The thermal entrance region in fully developed turbulent flow. AIChE J. 6 (2), 268-273. http://dx.doi.org/10.1002/aic.690060220.

Al-Arabi, M., 1982. Turbulent heat transfer in the entrance region of a tube. Heat Transf. Eng. 3 (3-4), 76-83. http://dx.doi.org/10.1080/01457638108939586.

Alatrash, Y., Cho, Y.J., Yoon, H.Y., Song, C.H., Chu, I.C., 2022. Experimental and numerical investigation of local bubble parameters for subcooled flow boiling in a pressurized annulus. Int. J. Heat Mass Transfer 194, 123040. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2022.123040.

Almeras, E., Mathai, V., Lohse, D., Sun, C., 2017. Experimental investigation of the turbulence induced by a bubble swarm rising within incident turbulence. J. Fluid Mech. 825, 1091-1112. http://dx.doi.org/10.1017/jfm.2017.410.

Angeli, P.E., Bieder, U., Fauchet, G., 2015. Overview of the TrioCFD code: Main features, v&v procedures and typical application to nuclear engineering. In: NURETH-16. http://dx.doi.org/10.1615/978-1-56700-099-3.590, URL https://cea.hal.science/cea-02500815.

Basu, N., Warrier, G.R., Dhir, V.K., 2005. Wall heat flux partitioning during subcooled flow boiling: Part 1—model development. J. Heat Transfer 127 (2), 131-140. http://dx.doi.org/10.1115/1.1842784.

Béguin, C., Pelletier, É., Étienne, S., 2016. Void fraction influence on added mass in a bubbly flow. Eur. J. Mech. B Fluids 56, 28-45. http://dx.doi.org/10.1016/j.euromechflu.2015.11.008.

Bell, I.H., Wronski, J., Quoilin, S., Lemort, V., 2014. Pure and pseudo-pure fluid thermophysical property evaluation and the open-source thermophysical property library CoolProp. Ind. Eng. Chem. Res. 53 (6), 2498-2508. http://dx.doi.org/10.1021/ie4033999.

Bergeaud, V., Lefebvre, V., 2010. SALOME. A software integration platform for multi-physics, pre-processing and visualisation. In: Proceedings of SNA + MC2010: Joint International Conference on Supercomputing in Nuclear Applications + Monte Carlo Tokyo.

Berry, R.A., Peterson, J.W., Zhang, H., Martineau, R.C., Zhao, H., Zou, L., Andrs, D., Hansel, J., 2018. Relap-7 Theory Manual. Technical Report, Idaho National Lab.(INL), Idaho Falls, ID (United States).

Bestion, D., Lucas, D., Boucker, M., Anglart, H., Tiselj, I., Bartosiewicz, V., 2009. Some lessons learned from the use of two-phase CFD for nuclear reactor thermalhydraulics. In: The 13th International Topical Meeting on Nuclear Reactor Thermal Hydraulics. NURETH-13.

Bozzano, G., Dente, M., 2001. Shape and terminal velocity of single bubble motion: a novel approach. Comput. Chem. Eng. http://dx.doi.org/10.1016/S0098-1354(01)00636-6.

Burns, A.D., Frank, T., Hamill, I., Shi, J.-M., 2004. The favre averaged drag model for turbulent dispersion in Eulerian multi-phase flows. In: 5th International Conference on Multiphase Flow, vol. 4, pp. 1-17, URL http://www.drthfrank.de/publications/2004/Burns_Frank_ICMF_2004_final.pdf.

Carlson, J.-R., Vatsay, V.N., Whitey, J., 2015. Node-centered wall function models for the unstructured flow code Fun3D. In: 22nd AIAA Computational Fluid Dynamics Conference. p. 2758.

Chen, Y.M., Mayinger, F., 1992. Measurement of heat transfer at the phase interface of condensing bubbles. Int. J. Muhiphase Flow 18 (6), 877-890. http://dx.doi.org/10.1016/0301-9322(92)90065-0.

Chu, I.-C., Lee, S.-J., Youn, Y.J., Park, J.K., Choi, H.S., Euh, D.-J., Song, C.-H., 2017. Experimental evaluation of local bubble parameters of subcooled boiling flow in a pressurized vertical annulus channel. Nucl. Eng. Des. 312, 172-183. http://dx.doi.org/10.1016/j.nucengdes.2016.06.027.

Chuang, T.-J., Hibiki, T., 2017. Interfacial forces used in two-phase flow numerical simulation. Int. J. Heat Mass Transfer 113, 741-754. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.05.062.

Clift, R., Grace, J.R., Weber, M.E., 2013. Bubbles, drops, and particles.

Cole, R., 1960. A photographic study of pool boiling in the region of the critical heat flux. AIChE J. 6 (4), 533-538. http://dx.doi.org/10.1002/aic.690060405.

Colin, C., Fabre, J., Kamp, A., 2012. Turbulent bubbly flow in pipe under gravity and microgravity conditions. J. Fluid Mech. 711, 469-515. http://dx.doi.org/10.1017/jfm.2012.401.

Cubizolles, G., 1996. Etude Stereologique de la Topologie des Ecoulements Diphasiques a Haute Pression (Ph.D. thesis). Ecole Centrale Lyon, URL https://inis.iaea.org/collection/NCLCollectionStore/_Public/48/007/48007225.pdf.

Del Valle, V.H., Kenning, D.B.R., 1985. Subcooled flow boiling at high heat flux. Int. J. Heat Mass Transfer 28 (10), 1907-1920. http://dx.doi.org/10.1016/0017-9310(85)90213-3.

Delhaye, J.M., 2008. Thermohydraulique Des Réacteurs. EDP Sciences.

Doherty, J., Ngan, P., Monty, J., Chong, M., 2007. The development of turbulent pipe flow. In: 16th Australasian Fluid Mechanics Conference. School of Engineering, The University of Queensland.

Emonot, P., Souyri, A., Gmandrille, J., Barré, F., 2011. CATHARE-3: A new system code for thermal-hydraulics in the context of the NEPTUNE project. Nucl. Eng. Des. http://dx.doi.org/10.1016/j.nucengdes.2011.04.049.

Favre, L., 2023. Modeling and Simulation of the Boiling Crisis Within PWR at CFD Scale (Ph.D. thesis). Institut National Polytechnique de Toulouse-INPT, URL https://theses.hal.science/tel-04244931.

Favre, L., Puget, S., Mimouni, S., Colin, C., 2022. Neptune_Cfd simulations of DEBORAPromoteur experiments: Boiling freon in a vertical pipe with mixing vanes. In: The 19th International Topical Meeting on Nuclear Reactor Thermal Hydraulics. NURETH-19.

Franck, T., Jain, S., Matyushenko, A., Garbaruk, A., 2012. The OECD/NEA MATHISH-benchmark - CFD analysis of water flow through a 5x5 rod bundle with spacer grids using ansys fluent and ansys CFX. In: CFD4NRS-4, Conference on Experimental Validation and Application of CFD and CMFD Codes in Nuclear Reactor Technology, OECD/NEA and IAEA Workshop.

François, F., Delhaye, J.M., Clément, P., 2011. The distribution parameter C0 in the drift modeling of forces convective boiling. Multiphase Sci. Technol. 23 (1), 77-100. http://dx.doi.org/10.1615/MultScienTechn.v23.i1.40.

Francois, F., Djeridi, H., Barre, S., Kledy, M., 2021. Measurements of void fraction, liquid temperature and velocity under boiling two-phase flows using thermal-anemometry. Nucl. Eng. Des. 381, 111359. http://dx.doi.org/10.1016/j.nucengdes.2021.111359.

Frost, W., Dzakowic, G., 1967. An extension of the method for predicting incipient boiling on commercially finished surfaces. ASME.

Garnier, J., Manon, E., Cubizolles, G., 2001. Local measurements on flow boiling of refrigerant 12 in a vertical tube. Multiphase Sci. Technol. 13 (1&2), http://dx.doi.org/10.1615/MultScienTechn.v13.i1-2.10.

Gerschenfeld, A., Gorsse, Y., 2022. Development of a robust multiphase flow solver on general meshes; application to sodium boiling at the subchannel scale. In: NURETH-19.
Gueguen, J., 2013. Contribution À La Modélisation Multidimensionnelle Des Écoulements Bouillants Convectifs En Conduite Haute Pression Pour L'application Au Cas Des Réacteurs À Eau Pressurisée (Ph.D. thesis). Université de Grenoble, URL https://theses.hal.science/tel-01685226/.

Hibiki, T., Ishii, M., Xiao, Z., 2001. Axial interfacial area transport of vertical bubbly flows. Int. J. Heat Mass Transfer 44, 1869–1888. http://dx.doi.org/10.1016/S0017-9310(00)00232-5.

Hinze, J.O., 1955. Fundamentals of the hydrodynamic mechanism of splitting in dispersion processes. AlChE J. 1 (3), 289–295. http://dx.doi.org/10.1002/aic.690010303.

Hosler, E.R., 1967. Flow Patterns in High Pressure Two-Phase (Steam-Water) Flow with Heat Addition. Technical Report WAPD-T-1824, Bettis Atomic Power Lab, Pittsburgh, Pa., URL https://www.osti.gov/biblio/4844120.

Ishii, M., Hibiki, T., 2006. Thermo-Fluid Dynamics of Two-Phase Flow. Springer Science and Business Media.

Ishii, M., Zuber, N., 1979. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AlChE J. 25 (5), http://dx.doi.org/10.1002/aic.690250513.

Jens, W.H., Lottes, P.A., 1951. Analysis of Heat Transfer, Burnout, Pressure Drop and Density Date for High-Pressure Water. Technical Report No. ANL-4627, Argonne National Lab.(ANL), Argonne, IL (United States), http://dx.doi.org/10.2172/4421630.

Kader, B., 1981. Temperature and concentration profiles in fully turbulent boundary layers. Int. J. Heat Mass Transfer 24 (9), 1541–1544. http://dx.doi.org/10.1016/0017-9310(81)90220-9.

Kim, S.J., Park, G.C., 2011. Interfacial heat transfer of condensing bubble in subcooled boiling flow at low pressure. Int. J. Heat Mass Transfer 54, 2962–2974. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2011.03.001.

Knopp, T., Alrutz, T., Schwamborn, D., 2006. A grid and flow adaptive wall-function method for RANS turbulence modelling. J. Comput. Phys. 220, 19–40. http://dx.doi.org/10.1016/j.jcp.2006.05.003.

Kok, J., 1999. Resolving the Dependence on Free-Stream Values for the K-Omega Turbulence Model. Technical Report NLR-TP-99295, National Aerospace Laboratory NLR, http://dx.doi.org/10.2514/2.1101.

Kommajosyula, R., 2020. Development and Assessment of a Physics-Based Model Forsubcooled Flow Boiling with Application to CFD (Ph.D. thesis). Massachusetts Institute of Technology, URL https://dspace.mit.edu/handle/1721.1/129051.

Končar, B., Morel, C., Mimouni, S., Vyskocil, L., Galassi, M., 2011. Computational fluid dynamics modeling of boiling bubbly flow for departure from nucleate boiling investigations. Multiphase Sci. Technol. 23 (2), 165. http://dx.doi.org/10.1615/MultScienTechn.v23.i2-4.40.

Krepper, E., Lucas, D., Frank, T., Prasser, H.M., Zwart, P.J., 2008. The inhomogeneous MUSIG model for the simulation of polydispersed flows. Nucl. Eng. Des. 238 (7), 1690–1702. http://dx.doi.org/10.1016/j.nucengdes.2008.01.004.

Kurul, N., Podowski, M., 1990. Multidimensional effects in forced convection subcooled boiling. In: Inc., B.H. (Ed.), International Heat Transfer Conference Digital Library. pp. 21–26. http://dx.doi.org/10.1615/IHTC9.40.

Labunstov, D.A., Lobachev, A.G., Zakharova, E.A., Kol'chugin, B.A., 1974. Influence of subcooling of liquid at the inlet of a heated channel on true volumetric steam content. Thermal Eng. (English Translation of Toploenergetica) 21 (9), 22–24.

Legendre, D., Magnaudet, J., 1998. The lift force on a spherical bubble in a viscous linear shear flow. J. Fluid Mech. 368, 81–126. http://dx.doi.org/10.1017/S0022112098001621.

Leoni, N., 2022. Bayesian Inference of Model Error for the Calibration of Two-Phase CFD Codes (Ph.D. thesis). Institut Polytechnique de Paris.

Leoni, N., Le Maitre, O., Rodio, M.-G., Congedo, P.M., 2024. Bayesian calibration with adaptive model discrepancy. Int. J. Uncertain. Quantif. 14 (1), http://dx.doi.org/10.1615/Int.J.UncertaintyQuantification.2023046331.

Liao, Y., Krepper, E., Lucas, D., 2019. A baseline closure concept for simulating bubbly flow with phase change: A mechanistic model for interphase heat transfer coefficient. Nucl. Eng. Des. 348, 1–13. http://dx.doi.org/10.1016/j.nucengdes.2019.04.007.

Lubchenko, N., Magolan, B., Sugrue, R., Baglietto, E., 2018. A more fundamental wall lubrication force from turbulent dispersion regularization for multiphase CFD applications. Int. J. Multiph. Flow 98, 36–44. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2017.09.003.

Lucas, D., Krepper, E., Rzehak, R., Liao, Y., Ma, T., Ziegenhein, T., 2015. Status and challenges of CFD-modelling for poly-disperde bubbly flow. In: NURETH-16.

Marfaing, O., Guingo, M., Laviéville, J., Bois, G., Méchitoua, N., Mérigoux, N., Mimouni, S., 2016. An analytical relation for the void fraction distribution in a fully developed bubbly flow in a vertical pipe. Chem. Eng. Sci. 152, 579–585. http://dx.doi.org/10.1016/j.ces.2016.06.041.

Martin, R., 1972. Measurement of the local void fraction at high pressure in a heating channel. Nucl. Sci. Eng. 48 (2), 125–138. http://dx.doi.org/10.13182/NSE72-A22466.

Masuk, A.U.M., Salibindla, A.K., Ni, R., 2021. Simultaneous measurements of deforming Hinze-scale bubbles with surrounding turbulence. J. Fluid Mech. 910, A21. http://dx.doi.org/10.1017/jfm.2020.933.

Menter, F.R., 1993. Zonal two equation k-cl, turbulence models for aerodynamic flows. In: American Institute of Aeronautics and Astronautics 24th Fluid Dynamics Conference. NASA-Ames Research Center, American Institute of Aeronautics and Astronautics, AIAA 93-2906.

Mimouni, S., Baudry, C., Guingo, M., Lavieville, J., Merigoux, N., Mechitoua, N., 2016. Computational multi-fluid dynamics predictions of critical heat flux in boiling flow. Nucl. Eng. Des. 299, 28–36. http://dx.doi.org/10.1016/j.nucengdes.2015.07.017.

NRC, U., 2010. TRACE v5.0 Theory Manual, Field Equations, Solution Methods, and Physical Models. Technical Report, United States Nucl. Regul. Comm. URL https://www.nrc.gov/docs/ML1200/ML120060218.pdf.

Park, I., Cho, H., Yoon, H., Jeong, J., 2009. Numerical effects of the semi-conservative form of momentum equations for multi-dimensional two-phase flows. Nucl. Eng. Des. 239, 2365–2371. http://dx.doi.org/10.1016/j.nucengdes.2009.06.011.

Pham, M., Bois, G., Francois, F., Baglietto, E., 2023. Assessment of state-of-the-art multiphase CFD modeling for subcooled flow boiling in reactor applications. Nucl. Eng. Des. 411, 112379. http://dx.doi.org/10.1016/j.nucengdes.2023.112379.

Ranz, W.E., Marshall, W., 1952. Evaporation from droplets. Chem. Eng. Prog. 48 (3), 141–146.

Reichardt, H., 1951. Vollständige darstellung der turbulenten geschwindigkeitsverteilung in glatten leitungen. Z. Angew. Math. Mech. 31 (7), 208–219. http://dx.doi.org/10.1002/zamm.19510310704.

Reiss, C., 2024. Burns turbulent dispersion considers the dispersed phase as a passive scalar. HAL archive. URL https://hal.science/hal-04494644.

Reiss, C., Gerschenfeld, A., Colin, C., 2024. Heat flux partition based on onset of significant void. HAL archive. URL https://hal.science/hal-04524455. Submitted to International Journal of Multiphase Flow.

Risso, F., Fabre, J., 1998. Oscillations and breakup of a bubble immersed in a turbulent field. J. Fluid Mech. 372, 323–355. http://dx.doi.org/10.1017/S0022112098002705.

Roy, R., Kang, S., Zarate, J., Laporta, A., 2002. Turbulent subcooled boiling flow—Experiments and simulations. J. Heat Transfer 124 (1), 73–93. http://dx.doi.org/10.1115/1.1418698.

Saha, P., Zuber, N., 1974. Point of net vapor generation and vapor void fraction in subcooled boiling. In: Inc., B.H. (Ed.), International Heat Transfer Conference Digital Library. pp. 175–179, URL https://www.nrc.gov/docs/ML1733/ML17338A800.pdf.

Salko, Jr., R., Avramova, M., Wysocki, A., Hizoum, B., Toptan, A., Hu, J., Porter, N., Blyth, T.S., Dances, C.A., Gomez, A., et al., 2023. CTF Theory Manual: Version 4.3. Technical Report, Oak Ridge National Laboratory (ORNL), Oak Ridge, TN (United States), URL https://www.osti.gov/biblio/1994732.

Sugrue, R., 2017. A Robust Momentum Closure Approach for Multiphase Computational Fluid Dynamics Applications (Ph.D. thesis). Massachusetts Institute of Technology.

Sugrue, R., Magolan, B., Lubchenko, N., Baglietto, E., 2017. Assessment of a simplified set of momentum closure relations for low volume fraction regimes in STAR-CCM+ and OpenFOAM. Ann. Nucl. Energy 110, 79–87. http://dx.doi.org/10.1016/j.anucene.2017.05.059.

Thom, J.R.S., Walker, W., Fallon, T.A., Reising, G.F.S., 1965. Paper 6: Boiling in subcooled water during flow up heated tubes or annuli. In: Sage UK: London, E.S.P. (Ed.), In: Proceedings of the Institution of Mechanical Engineers, Conference Proceedings, vol. 180–3, pp. 226–246. http://dx.doi.org/10.1016/0301-9322(81)90040-9, URL https://www.osti.gov/biblio/4263900.

Thomas, R., 1981. Bubble coalescence in turbulent flows. Int. J. Multiph. Flow 7 (6), 709–717. http://dx.doi.org/10.1016/0301-9322(81)90040-9.

Todreas, N.E., Kazimi, M.S., 2021. Nuclear Systems I: Thermal Hydraulic Fundamentals. CRC Press.

Tomiyama, A., Kataoka, I., Zun, I., Sakaguchi, T., 1998. Drag coefficients of single bubbles under normal and micro gravity conditions. JSME Int. J. Ser. B Fluids Therm. Eng. http://dx.doi.org/10.1299/jsmeb.41.472.

Tomiyama, A., Tamai, H., Zun, I., Hosokawa, S., 2002. Transverse migration of single bubbles in simple and shear flows. Chem. Eng. Sci. 57, 1849–1858. http://dx.doi.org/10.1016/S0009-2509(02)00085-4.

Ünal, H., 1976. Maximum bubble diameter, maximum bubble-growth time and bubble-growth rate during the subcooled nucleate flow boiling of water up to 17.7 MN/m2. Int. J. Heat Mass Transfer http://dx.doi.org/10.1016/0017-9310(76)90047-8.

Vlček, D., Sato, Y., 2023. Sensitivity analysis for subcooled flow boiling using Eulerian CFD approach. Nucl. Eng. Des. 405, 112194. http://dx.doi.org/10.1016/j.nucengdes.2023.112194.

Wallis, G.B., 1974. The terminal speed of single drops or bubbles in an infinite medium. Int. J. Multiph. Flow 1 (4), 491–511. http://dx.doi.org/10.1016/0301-9322(74)90003-2.

Yao, W., Morel, C., 2004. Volumetric interfacial area prediction in upward bubbly two-phase flow. Int. J. Heat Mass Transfer 47 (2), 307–328. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2003.06.004.

Yoon, S.J., Agostinelli, G., Baglietto, E., 2017. Assessment of multiphase CFD with zero closure model for boiling water reactor fuel assemblies. In: 17th International Topical Meeting OnNuclear Reactor Thermal Hydraulics.

Yuan, C., Laurent, F., Fox, R., 2012. An extended quadrature method of moments for population balance equations. J. Aerosol Sci. 51, 1–23. http://dx.doi.org/10.1016/j.jaerosci.2012.04.003.

Zeitoun, O., Shoukri, M., Chatoorgoon, V., 1995. Interfacial heat transfer between steam bubbles and subcooled water in vertical upward flow. In: Trans. ASME. J. Heat Transfer 117 (2), 402–407. http://dx.doi.org/10.1115/1.2822536.

Zuber, N., 1964. On the dispersed two-phase flow in the laminar flow regime. Chem. Eng. Sci. 19, 897–917. http://dx.doi.org/10.1016/0009-2509(64)85067-3.