# A Diameterless Boiling-Flow Multiphase CFD Framework for Nuclear Reactor Conditions

Corentin Reiss, Antoine Gerschenfeld, Catherine Colin

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

## To cite this version:

Corentin Reiss, Antoine Gerschenfeld, Catherine Colin. A Diameterless Boiling-Flow Multiphase CFD Framework for Nuclear Reactor Conditions. 2024. hal-04527149

HAL Id: hal-04527149

https://hal.science/hal-04527149

Preprint submitted on 29 Mar 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Highlights

## A Diameterless Boiling-Flow Multiphase CFD Framework for Nuclear Reactor Conditions

Corentin Reiss, Antoine Gerschenfeld, Catherine Colin

* In a high-pressure developed boiling ascending pipe flow, changing the inlet temperature is equivalent to moving a probe up and down

* Boiling flows are simulated using a 3D map of the experimental bubble diameter

* Atmospheric-pressure adiabatic force closure models are not valid for high-pressure conditions

* A simplified set of closure models that doesn't require bubble diameter modeling is proposed
# A Diameterless Boiling-Flow Multiphase CFD Framework for Nuclear Reactor Conditions

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

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

<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

## 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 interfacial area closures, 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. Atmospheric-pressure closure terms are shown not to be able to reproduce measured void fraction profiles. We propose a new set of closures that is based on the hypothesis that bubbles are deformable in nuclear reactor conditions. This enables us to avoid bubble diameter modeling through 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 improved compared with the baseline set of closures.

Email addresses: corentin.reiss@cea.fr (Corentin Reiss),
antoine.gerschenfeld@cea.fr (Antoine Gerschenfeld),
catherine.colin@toulouse-inp.fr (Catherine Colin)
Keywords: Multiphase CFD, Interfacial forces, Interfacial heat flux

PACS: 64.70 fh, 44.05.+e, 44.35.+c, 47.27.T, 47.55.-t

2000 MSC: 80A20, 76T10

# 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)). Therefore, the models commonly used are mostly derived from simple situations very different from reactors, like the study of adiabatic single bubbles at atmospheric pressure (Tomiyama et al. (2002)).

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 ($<\alpha_v>= 0.7$). We therefore concentrate on modeling bubbly flow.

Our approach consists in exploiting existing experiments in such a way as to decouple the contributions of aforementioned terms. To avoid uncertainty related to interfacial area transport equations (Yao and Morel (2004)) or population balance (Krepper et al. (2008); Yuan et al. (2012)) methodologies, we use bubbly flows with measured mean Sauter diameters and enforce the experimental values as in Sugrue (2017). 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 figure 1 and section 3).
<table><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>

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

<figure><img src="image_4.png" /><figcaption>Figure 1: Diagrams of the two-phase experiments used for framework validation and heat transfer models evaluation.</figcaption></figure>

We then study the DEBORA experiment (Garnier et al. (2001)) (see figure 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, heat fluxes and inlet subcoolings. In order to use the measured diameters as inputs for the code, we show that in a high-pressure developed flow, changing the entrance temperature in a boiling pipe is equivalent to changing the distance from the entry where the flow characteristics are measured (section 4). We run simulations using the standard closure terms calibrated at atmospheric pressures, which yield unsatisfying results (section 5). To improve these models, we hypothesize that in nuclear reactor conditions, even small bubbles are deformable. This enables us to propose a set of closures that doesn't require bubble diameter modeling and improves simulation predictions (section 6).

## 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 can handle 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 \alpha_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} = \frac{\partial \alpha_k \rho_k e_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k e_k \vec{u}_k) = $$

Where $\alpha_k$ is the 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 and $T_k$ the fluid temperature.

In equation 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 term $\vec{F}_{ki}$ and the interfacial heat transfer $q_{ki}$. The wall heat transfer to phase k, $q_{kw}$, 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{align}
& \nu_t = \frac{k}{\omega} \qquad \overline{u_i' u_j'} = -\nu_t \nabla \vec{u}_l \qquad \overline{u_i' e_l'} = -\nu_t C_{pl} \nabla T_l \nonumber \\
& \partial_t(k) + \nabla \cdot (k \vec{u}_l) = \nu_t (\nabla \vec{u}_l + {}^t \vec{u}_l) \cdot \nabla \vec{u}_l - \beta_k k \omega + \nabla \cdot (\alpha_l (\nu_l + \sigma_k \nu_t) \nabla k) \tag{2} \\
& \partial_t \omega + \nabla \cdot (\omega \vec{u}_l) = \alpha_\omega (\nabla \vec{u}_l + {}^t \nabla \vec{u}_l) \cdot \nabla \vec{u}_l - \beta_\omega \omega^2 \nonumber \\
& \qquad + \nabla \cdot (\alpha_l (\nu_l + \sigma_\omega \nu_t) \nabla \omega) + \sigma_d \frac{1}{\omega} \max\{\nabla k \cdot \nabla \omega, 0\} \nonumber
\end{align}
$$

Where $k$ is the turbulent kinetic energy of the liquid, $\omega$ it's turbulent dissipation rate, $\nu_t$ it's turbulent viscosity and $\nu_l$ it's dynamic viscosity. The values of the constants are $\alpha_\omega = 0.5$, $\beta_k = 0.09$, $\beta_\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 u_\tau^2$.

The boundary condition on $k$ is $k = 0$ at the wall for $y_+ < 5$, where $y_+ = yu_\tau/\nu_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{align}
q_{wl} &= (T_w - T_l(y)) \frac{\rho_l C_{pl} u_\tau}{\Theta_+^w(y_+)} \nonumber \\
\beta_{\text{SP}} &= (3.85(Pr^{1/3}) - 1.3)^2 + 2.12 \log(Pr) \nonumber \\
\gamma &= \frac{0.01(Pr y_+)^4}{1 + 5Pr^3 y_+} \tag{3} \\
\Theta_+^w(y_+) &= \text{Pry}_+ \exp(-\gamma) + (2.12 \cdot \log(1 + y_+) + \beta_{\text{SP}}) \exp\left(-\frac{1}{\gamma}\right) \nonumber
\end{align}
$$

Where $T_w$ is the wall temperature, $C_{pl}$ is the heat capacity of the liquid, $Pr = \frac{\nu_l \cdot \rho C_{pl}}{\lambda_l}$ the Prandtl number, with $\lambda_l$ the liquid conductivity. This expression
is based on experimental measures for $y_+$ ranging from 0 to 300. Using the previously calculated expression of $u_τ$, we calculate $θ_+$ 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 their volumetric form. 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_i \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^{687}), \frac{8Eo}{3Eo+12} \right) $$

Where $C_D$ is the drag coefficient, $a_i = \frac{6\alpha_g}{d_b}$ is the interfacial area of the gas phase, $Eo = \frac{(\rho_l - \rho_g)gd_b^2}{\sigma}$ the Eotvos number and $Re_b = \frac{d_b ||\vec{u}_g - \vec{u}_l||}{\nu_l}$ the bubble Reynolds number.
Virtual mass force. The virtual mass writes:

$$ \vec{F}_{\text{VM}} = -C_{\text{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.

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. (2002)), the local void fraction can reach 0.7. In either above formulation, a liquid fraction of at least $0.7 \cdot 1/2 = 0.35$ would be entrained even though the total liquid fraction is 0.3. 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}_{\text{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) & \max(0, g(\alpha) &= 1.0155 - 0.0154\exp(8.0506\alpha)) \\ Wo &= Eo \frac{k}{\|\vec{u}_g - \vec{u}_l\|^2} & f(Wo) &= \min(0.03, 5.0404 - 5.0781Wo^{0.0108}) \end{aligned} \quad (9) $$

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

$$ \vec{F}_{\text{disp}} = -C_{\text{TD}} \rho_l k \nabla \alpha_g, \quad C_{\text{TD}} = \frac{3 C_D}{4 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:

$$C_L \rightarrow \begin{cases} 0 & \text{if } y/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{if } y/d_b \ge 1 \end{cases} \quad (11)$$

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

$$\vec{F}_{\text{wall}} = \begin{cases} C_{TD} \rho_l k \cdot \alpha_g \frac{1}{y} \frac{d_b - 2y}{d_b - y} \vec{n} & \text{if } y/d_b < 1/2 \\ 0 & \text{if } y/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$:

$$q_{ki} = a_i \frac{\lambda_l}{d_b} (T_g - T_l) Nu = \frac{6\alpha_v \lambda_l}{d_b} (T_g - 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 \text{Wm}^{-3}$) if $T_l > T_s$.

*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:

$$\begin{aligned} q_{KP} &= q_c + q_q + q_{w,l \rightarrow v} \\ q_{wl} &= q_c + q_q \end{aligned} \quad (15)$$
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) (equation 3):

$$ 
\begin{aligned}
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 g(\rho_l - \rho_g)}{3 \rho_l d_{b, det}}} \quad , \quad A_b = \min(1, \pi/4 \cdot N_s d_{b, det}^2)
\end{aligned} 
\quad (16) $$

This yields:

$$ 
\begin{aligned}
q_c &= (1 - A_b)q_{SP} \\
q_q &= 2A_b \lambda_l (T_w - T_l) \sqrt{\frac{f_{dep} \rho_l C_{pl}}{\pi \lambda_l}} \\
q_{w,l \to v} &= \frac{\pi}{6} f_{dep} d_{b, det}^3 \rho_g h_{lg} N_s
\end{aligned} 
\quad (17) $$

## 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 figure 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 figure 3-A,B). The simulated 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 figure 3-C). The liquid temperature from the experiment and the simulation have the same profile shape, though they are off by ~1°C.
<figure><img src="image_5.png" /><figcaption>Figure 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_\tau$ as a function of $y_+$. Black line: Reichardt (1951) solution. B. Non-dimensional turbulent kinetic energy $k_+ = k/u_\tau^2$ as a function of $y_+$. Black line: refined solution for $y_{+,1} = 1$. C. Non-dimensional dissipation rate $\omega_+ = \omega u_\tau / \nu$ as a function of $y_+$. Black line: Knopp et al. (2006) solution.</figcaption></figure>

This amounts to a 5% power loss on the 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. Each mesh can have 7, 14 or 28 radial elements, and 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, therefore to save computation time we only simulate slices in the rest of this paper.

### 3.4. *Two-phase adiabatic vertical tube*

Multiple experimental databases are available to study 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. We also select the Colin et al. (2012) experiments for downwards
<figure><img src="image_6.png" /><figcaption>Figure 3: Single-phase validation of the multiphase module of TrioCFD. A. Colin et al. (2012) single-phase axial velocity in a 4cm-diameter pipe. B. Colin et al. (2012) single-phase axial velocity fluctuations in a 4cm-diameter pipe. C. Temperature profiles from single-phase DEBORA runs for a 19.2mm diameter pipe, 2.61MPa pressure, 1996kg/(m²s) mass velocity and 74.4kW/m² heat flux with inlet temperatures of 19.84°C and 31.46°C (page 109 in Garnier et al. (2001)).</figcaption></figure>

and microgravity flow (see figure 1). We run different test cases using the interfacial force models described in section 2.3 (see figure 4). To avoid modeling the interfacial area, we enforce the radially-dependent steady-state experimental diameter in the simulations. The complete model is able to predict correctly void fraction profiles for low (figure 4-A) and high (figure 4-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 figure 4-B. Furthermore, in figure 4-D the gas velocities are well predicted by the model.

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 figure 1). It consists of a vertical Freon-12 pipe with a 1m-long inlet section, a 3.5m-long heated section and an instrumentation plane located 3.485m after the beginning of the heated section. During different campaigns, liquid temperatures, void fractions, gas velocities and Sauter mean diameters were measured at the outlet. 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)).
<figure><img src="image_7.png" /><figcaption>Figure 4: Two-phase adiabatic simulated and experimental void fraction distributions in a pipe. We use the interfacial force models described in section 2.3. 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>
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 total heat flux 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. In this section, we show that this hypothesis is correct 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 figure 5). Latent heat is $h_{lg}$. Then, as described in Pham et al. (2023):

$$X(z) = X_{in} + \frac{\mathcal{P}q_w}{G\mathcal{A}h_{lg}}z \qquad (18)$$

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

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

A diagram of this situation is presented in figure 5.

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 effect of hydrostatic pressure on saturation temperature, and the formation of the turbulent thermal boundary layer (section 4.2).

* a posteriori analysis: we run CFD simulations of different DEBORA test tubes and of a single tube with measures at equivalent altitude (section 4.3).
<figure><img src="image_8.png" /><figcaption>Figure 5: 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>

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

*Saturation temperature.* We perform a thought experiment. We consider a 5m long water test section, and two series of experiments at atmospheric and PWR pressures. For each, we assume two runs are performed: a hot run with an outlet subcoolings of 5°C and a cold run with 15°C subcooling. We suppose that the cold one run's equivalent altitude for the hot run configuration is at the inlet of the test section, and that the pressure is hydrostatic.

In the atmospheric pressure-experiment, the pressure at the bottom of the pipe is 1.5bar, resulting in a saturation temperature of 110°C, 10°C more than at 1bar (Bell et al. (2014)). Therefore, in the hot run at cold run-equivalent altitude the subcooling is 25°C, which is very different from the 15°C in the cold run. The wall boiling dynamics, and in turn the void fraction, liquid temperature and velocities. The equivalent altitude hypothesis is therefore not valid in this condition.

In the PWR pressure-run, the bottom of the pipe is at 155.5bar, resulting in a saturation temperature of 345.05°C, 0.27°C more than at the outlet (155bar and 344.78°C) (Bell et al. (2014)). This means that in the hot run at cold run-equivalent altitude, the subcooling is 15.27°C, very close to the 15°C in the cold run. Wall boiling dynamics will be very similar in a to the experimental condition.

Therefore, 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 before typical subcoolings.

The DEBORA experiments are conducted at a minimum pressure of
<figure><img src="image_9.png" /><figcaption>Figure 6: Void fraction as a function of equivalent altitude in the DEBORA 29 campaign. Each subplot represents a test tube. Vertical black lines represent $z/D_h=30$.</figcaption></figure>

14bar, with a heated length of 3.5m. The pressure variations between both conditions is smaller than 0.5bar. This changes the saturation temperature of R12-freon by ~1°C (Bell et al. (2014)), much smaller than the typical 20°C inlet subcoolings (Garnier et al. (2001)).

*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 to be similar at the outlet and at the equivalent altitude, the velocity and temperature fields must be the same 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 doesn't guarantee identical velocity and temperature profiles. For profiles to be the same, single-phase temperature and velocity profiles before the OSV must be developed. This condition is achieved after $z/D_h=30$, where the vertical black lines are placed in figure 6. To verify it, 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, the developed turbulence at OSV criterion is verified. However, this may not be the case in other databases, in particular when the heated section is much shorter.
*Use for experiment design.* To guarantee that the test tube hypothesis can be used in simulations, we recommend the following for future experiments:

* 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\text{bar}$).

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

4.3. *A posteriori analysis: numerical verification*

In order to verify that the equivalent altitude hypothesis is valid for the DEBORA experiment, we run two kinds of simulations for a single test tube. 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 figure 6, 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 14bar, 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.485m, the test section being 3.5m long.

Simulations are run using the models presented in section 2, with a 1mm-constant Sauter mean diameter. The mesh is a 1°-wide slice, with 20 radial and 500 axial cells. Simulation results are presented in figure 7.

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
<figure><img src="image_10.png" /><figcaption>Figure 7: 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>
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*<sub>j</sub> the outlet thermodynamic quality of run *j*. *j*<sub>max</sub> is the run of a given test tube at the highest inlet temperature *T*<sub>in,jmax</sub>. Using the test tube hypothesis, we launch calculations for each test tube with *T*<sub>in</sub> = *T*<sub>in,jmax</sub> and extract the physical parameters at the different *z*<sub>j</sub> so that:

$$X_{T_{\text{in}}=T_{\text{in},j\max}}(z_j) = X_j \qquad (20)$$

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*<sub>in,j</sub>. We then know the experimental value at any *z*<sub>j</sub>. 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 void fraction profiles are presented in figure 8, 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.
<table><caption>Table 1: 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)$.</caption><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,jmax</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>3.485</td><td>0.1091</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_11.png" /><figcaption>Figure 8: 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>

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 figure 9 for test tube II-G2P14W16. This test tube was selected 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 2-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)
<figure><img src="image_12.png" /><figcaption>Figure 9: 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 negative at high outlet qualities (see equation 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 lift coefficient
is near-zero, resulting in wall-peaked simulated profiles in all test tubes (see
figure 9).

*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 figure 9) 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 aren't core-peaked like the experiment because the vapor profile is wall-peaked, leading to a buoyancy effect at the wall (figure 9 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 (figure 9 bottom center), indicating that vapor is produced at the wall and migrates towards the center of the pipe.

*Force balance.* The force balances are at equilibrium because the flow is stationary. The axial one on the vapor phase (top right plot in figure 9) shows
that the pressure gradient is mainly compensated by the drag force, as ex-
pected. 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 as 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 figure 9) 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 α_v-dependent term in the Sugrue (2017) lift to be non-zero (see equation 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 pres-surized water reactor similarity conditions.

This isn't necessarily surprising, as most of these models are used outside
of their validity domain. Pressurized water reactor flows have high void
fractions, are extremely turbulent, have low surface tensions and density
ratios.

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*

Figure 10 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 equation 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: the velocity of bubbles in the near-wall cell is smaller
than the velocity of the liquid in the core. As bubbles migrate, their inertia prevents them from accelerating immediately, hence the negative relative velocity. This effect is exacerbated if vapor is produced with zero inertia, i.e. if an interfacial force $-\Gamma_{l \to v} \vec{u}_v$ is added in the near-wall cell (center column in figure 10). 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); Francois 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 figure 10). 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 figure 11). The first, $We_{drag}$, is defined using the relative velocity and accounts for the effect of the drag force on bubble deformation. The second, $We_\epsilon$, represents the impact of turbulent velocity fluctuations at the scale of the bubbles:

$$ 
\begin{aligned}
We_{drag} &= \frac{\rho_l d_b ||\vec{u}_g - \vec{u}_l||^2}{\sigma} \\
We_\epsilon &= \frac{\rho_l d_b (\epsilon d_b)^{2/3}}{\sigma}
\end{aligned}
\qquad (21) $$

In this computation $d_b$ comes from the experimental measures. The other physical quantities 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 drag Weber number $We_{drag}$ 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_\epsilon$ has an opposite behavior: it is large in the near-wall region and diminishes in the core. Bubble experience significant deformation when $We_\epsilon \ge 2$ (Hinze (1955); Risso and Fabre (1998); Masuk et al. (2021)). Therefore, on
<figure><img src="image_13.png" /><figcaption>Figure 10: 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_14.png" /><figcaption>Figure 11: Weber numbers from DEBORA experiment simulations using the experimental diameters as input but without the virtual mass force. See equation 21 for definitions, figure 8 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>
<figure><img src="image_15.png" /><figcaption>Figure 12: Distance between two spherical bubbles in a compact packing as a function of the void fraction.</figcaption></figure>

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 fraction, bubble diameter and distance between bubbles in then:

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

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 (23) $$

The dimensionless distance between bubbles is plotted as a function of $\alpha_v$ in figure 12. $d_{dist} \simeq d_b$ for $\alpha_v = 0.1$, and $d_{dist} \simeq 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 figure 8). In flow boiling, bubble agitation and turbulence mean that they necessarily collide, which will also lead to bubble deformation.

We call *deformed bubble hypothesis* the hypothesis that in flow boiling in PWR conditions, bubbles are deformed 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.* The drag force formulations most commonly used in the literature have a drag coefficient that depends on the bubble Reynold's number and diameter. As shown by Sugrue (2017), when the bubbles are sufficiently deformed the air-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)).

In addition to their often-used expressions for spherical bubbles, Ishii and Zuber (1979) also proposed a drag force formulation for deformed bubbles. We use this formulation for our set of closures:

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

*Lift force.* The lift coefficient is positive for shear flow around a spherical bubble, and become negative for deformable 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 2-fluid model. In figure 8, 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
<figure><img src="image_16.png" /><figcaption>Figure 13: The lift coefficient proposed in this work as a function of the void fraction.</figcaption></figure>

plotted in figure 13, 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} \quad (25)$$

*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 equation 10). Furthermore, we use the Ishii and Zuber (1979) deformable bubble drag force (see equation 24). 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 and anticipate.

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 (26)$$

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(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 (equation 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,\text{Thom et al.}} &= \left(\frac{T_w - T_s}{22.65} \exp(P/87)\right)^2 \\ q_{w,\text{Frost\&Dzakowic}} &= \frac{\lambda_{ls}(h_{gs} - h_{ls})\rho_v}{8\sigma T_s} \left(\frac{T_w - T_s}{Pr_s}\right)^2 \end{aligned} \quad (27) $$

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_{\text{Boil}} \\ q_{wl} &= 0 \\ q_{w,l \to v} &= q_{\text{Boil}} \end{aligned} \right. \quad (28) $$

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_{\text{SP}} \\ q_{wl} &= q_{\text{SP}} \\ q_{w,l \to v} &= 0 \end{aligned} \right. \quad (29) $$
5. Else we are between the onset of nucleate boiling and saturated boiling. The total heat flux is then $q_w = q_{\text{Boil}}$. To determine if we have passed the onset of significant void, we calculate $H_{l,OSV}$ (equation 26):

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

$$ \left\{ \begin{array}{l} q_w = q_{\text{Boil}} \\ q_{wl} = q_{\text{Boil}} \\ q_{w,l \to v} = 0 \end{array} \right. \qquad (30) $$

(b) Else we consider that the heat transfer to the liquid phase is the same as at OSV (equation 26), and that the rest of the energy is used for evaporation:

$$ \left\{ \begin{array}{l} q_w = q_{\text{Boil}} \\ q_{wl} = H_{l,OSV}(T_s - T_l(y_1)) \\ q_{w,l \to v} = q_{\text{Boil}} - H_{l,OSV}(T_s - T_l(y_1)) \end{array} \right. \qquad (31) $$

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_+) - 7} \right) \\ q_w &= \max(q_{SP}, q_{\text{Boil}}) \\ q_{wl} &= \min(q_w, q_{l,OSV}) \\ q_{w,l \to v} &= q_w - q_{wl} \end{aligned} \qquad (32) $$

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 anticipate and 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 of the simplified set of closures require bubble diameter modeling through and interfacial area transport equation (Yao and Morel (2004)) or population balance model (Krepper et al. (2008)). We therefore aim to build a correlation that doesn't 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(\rho_l-\rho_v)}}$, and all of them were between 0.25$L_c$ and 2$L_c$. For test tubes I, II and III, $L_c$ was 0.4, 0.65 and 0.45mm (compare with figure 9 bottom line). In all three test tubes, $d_b < L_c$ if $T_s - T_l > 1^\circ C$, i.e. in subcooled regions bubbles are smaller than $L_c$. We use the capillary length as the bubble diameter length scale used in the interfacial area ($a_i \sim 6\alpha_v/L_c$ in 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 (equation 23) 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_{\text{fit}} \frac{6\lambda_l}{L_c^2} \frac{1}{\min\left(1, \sqrt[3]{\frac{\pi}{3\sqrt{2}} \frac{1}{\min(\alpha_v, 0.6)}} - 1\right)} (T_v - T_l) \quad (33)$$

With $Nu_{\text{fit}}$ a fitted constant bubble Nusselt number. $Nu_{\text{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 equation 33 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 $\sim 10$ compared with low-volume fraction cases and in practice the liquid is at saturation temperature.

This correlation isn't 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); Kim and Park (2011). As discussed in section 2.4, all of them are based on equation 13. The specificity of each correlation then lies in the Nusselt number formulation, that can depend on the Prandtl, bubble Reynold's of
Jacob numbers:

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

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$ (Sugrue (2017)). 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_{ki} \propto d_b^{-2}$ that is visible in equation 13, particular in subcooled regions where $d_b < L_c$.

<table><caption>Table 2: Dependence of the interfacial heat transfer on the diameter for various correlations from the literature and for the current work.</caption><thead><tr><th>Reference</th><th>d<sub>b</sub> &lesssim; L<sub>c</sub></th><th>d<sub>b</sub> &gtrsim 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 figure 14. The inflection in the plots from the literature at $d_b \approx L_c$ come from the relative velocity becoming constant at 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 closure terms with 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. Figure B.19 in Appendix B 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 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: RELAP (Berry et al. (2018)), TRACE (NRC (2010)), CTF (Salko Jr 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 figure 15, and detailed results are shown for test tube II-G2P14W16 in figure 16. 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 us 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 (figure 16 lower left). The lift force is significantly larger in the radial force balance
<figure><img src="image_17.png" /><figcaption>Figure 14: 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 100bar.</figcaption></figure>
<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>8E<sub>o</sub></sup>/<sub>3E<sub>o</sub>+12</sub>)</td><td>Ishii and Zuber (1979)<br>deformable bubble<br>C<sub>D</sub> = <sup>2</sup>/<sub>3</sub>d<sub>b</sub>√<sup>g(ρ<sub>l</sub> - ρ<sub>v</sub>)</sup>/<sub>σ</sub></td></tr><tr><td>Lift</td><td>Sugrue (2017)<br>C<sub>L</sub> = f(W<sub>o</sub>) · 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> <sup>C<sub>D</sub></sup>/<sub>d<sub>b</sub></sub>|<u>u</u><sub>g</sub> - <u>u</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> <sup>C<sub>D</sub></sup>/<sub>d<sub>b</sub></sub>|<u>u</u><sub>g</sub> - <u>u</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>pl</sub>u<sub>τ</sub></sup>/<sub>Θ<sub>+</sub><sup>w</sup>(y<sub>+</sub>)</sub><br>N<sub>s</sub> = (210(T<sub>w</sub> - T<sub>s</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> = √<sup>4</sup>/<sub>3</sub> <sup>g(ρ<sub>l</sub> - ρ<sub>q</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>)√<sup>f<sub>dep</sub>ρ<sub>l</sub>C<sub>pl</sub></sup>/<sub>πλ<sub>l</sub></sub><br>q<sub>wl</sub> = q<sub>c</sub> + q<sub>q</sub><br>q<sub>w, 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, 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>pl</sub>u<sub>τ</sub></sup>/<sub>Θ<sub>+</sub><sup>w</sup>(y<sub>+</sub>)</sub><br>q<sub>Boil</sub> = (<sup>T<sub>w</sub> - T<sub>s</sub></sup>/<sub>22.65</sub> exp(P/87))<sup>2</sup><br>q<sub>l, OSV</sub> = max (<sup>0</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> = <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>Pr<sup>1/3</sup>)</td><td>Current work<br>q<sub>ki</sub> = <sup>6α<sub>v</sub>Nu<sub>fit</sub>λ<sub>l</sub>g(ρ<sub>l</sub> - ρ<sub>v</sub>)(T<sub>g</sub> - T<sub>l</sub>)</sup>/<sub>σ min(1, <sup>3</sup>/<sub>√3</sub> <sup>π</sup>/<sub>3√2</sub> min(α<sub>v</sub>, 0.6)</sub> - 1</sup></td></tr></tbody></table>
<figure><img src="image_18.png" /><figcaption>Figure 15: Void fraction simulation results for the set of closures proposed in this paper.</figcaption></figure>

than in the baseline closure (figure 16 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 figure 9 lower right).

Pressurized water reactor subchannels can be seen as the assembly of 3 different elementary shapes (see figure 17): a tube, in the center of a sub-channel; a 3mm-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 figure 18. 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.

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 isn't necessarily surprising. The Nus-
<figure><img src="image_19.png" /><figcaption>Figure 16: 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>

<figure><img src="image_20.png" /><figcaption>Figure 17: Diagram of a pressurized water reactor subchannel (fuel rods are in grey) and dimensions of the different experiments used for validation of CFD codes.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 18: 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.5mm inner diameter 27.2mm 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.8mm-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>
selt number of condensation formulations from the literature depend on fluid
characteristics through the Reynold's, Prandtl and Jacob numbers. More-
over, their characteristic Nusselt number are larger for water than for refrig-
erant fluids (see top row in figure 14). 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.

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 tem-
perature 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 isn't 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 doesn't 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 working on a critical heat flux criterion.

8. Acknowledgments

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

9. Copyright

For the purpose of Open Access, a CC-BY public copyright license has
been applied by the authors to the present document and will be applied to
all subsequent versions up to the Author Accepted Manuscript arising from
this submission.

<figure><img src="image_22.png" /></figure>
References

Alatrash, Y., Cho, Y.J., Yoon, H.Y., Song, C.H., Chu, I.C., 2022. Experimen-
tal and numerical investigation of local bubble parameters for subcooled
flow boiling in a pressurized annulus. International Journal of Heat and
Mass Transfer 194, 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.

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. URL: https://cea.hal.science/cea-02500815, doi:http://dx.doi.org/10.1615/978-1-56700-099-3.590.

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, 131–140. doi:https://doi.org/10.1115/1.1842784.

Béguin, C., Pelletier, E., Étienne, S., 2016. Void fraction influence on added mass in a bubbly flow. European Journal of Mechanics B/Fluids 56, 28–45.

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. Industrial & engineering chemistry research 53, 2498–2508.

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 and approach. *Computers and Chemical Engineering*.

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, 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, 877-890.

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 sub-cooled boiling flow in a pressurized vertical annulus channel. *Nuclear Engineering and Design* 312, 172-183.

Chuang, T.J., Hibiki, T., 2017. Interfacial forces used in two-phase flow numerical simulation. *International Journal of Heat and Mass Transfer* 113, 741-754.

Cole, R., 1960. A photographic study of pool boiling in the region of the critical heat flux. *AIChE Journal* 6, 533-538. doi:https://doi.org/10.1002/aic.690060405.

Colin, C., Fabre, J., Kamp, A., 2012. Turbulent bubbly flow in pipe under gravity and microgravity conditions. *Journal of Fluid Mechanics* 711, 469-515.

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. *International Journal of Heat and Mass Transfer* 28, 1907-1920. doi:https://doi.org/10.1016/0017-9310(85)90213-3.
Delhaye, J.M., 2008. Thermohydraulique des réacteurs. EDP Sciences.

Emonot, P., Souyri, A., Gandrille, J., Barré, F., 2011. Cathare-3: A new system code for thermal-hydraulics in the context of the neptune project. Nuclear Engineering and Design doi: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 debora-promoteur 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 mathis-h 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 Science and technology 23, 77-100. doi: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. Nuclear Engineering and Design 381, 111359. doi:https://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, C., Lance, M., Marie, J., 2002. Measurement of local flow characteristics in buoyancy-driven bubbly flow at high void fraction. Experimental Thermal and Fluid Science 26, 811-815.

Garnier, J., Manon, E., Cubizolles, G., 2001. Local measurements on flow boiling of refrigerant 12 in a vertical tube. Multiphase Science and Technology 13. doi: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. International Journal of Heat and Mass Transfer 44, 1869–1888.

Hinze, J.O., 1955. Fundamentals of the hydrodynamic mechanism of splitting in dispersion processes. AIChE journal 1, 289–295.

Hosler, E.R., 1967. Flow patterns in high pressure two-phase (steam-water) flow with heat addition. Technical Report WAPD-T-1824. Betis 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. AIChE Journal 25. doi:https://doi.org/10.1002/aic.690250513.

Jens, W.H., Lottes, P.A., 1951. Analysis of heat transfer, burnout, pressure drop and density change for high-pressure water. Technical Report No. ANL-4627. Argonne National Lab. (ANL), Argonne, IL (United States). doi:10.2172/4421630.

Kader, B., 1981. Temperature and concentration profiles in fully turbulent boundary layers. Int. J. Heat Mass Transfer 24, 1541–1544. doi:https://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. International Journal of Heat and Mass Transfer 54, 2962–2974.
Knopp, T., Alrutz, T., Schwamborn, D., 2006. A grid and flow adaptive wall-function method for rans turbulence modelling. Journal of Computational Physics 220, 19–40.

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. doi:https://doi.org/10.2514/2.1101.

Kommajosyula, R., 2020. Development and assessment of a physics-based model for subcooled 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 Science and Technology 23, 165. doi: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. Nuclear Engineering and Design 238, 1690–1702.

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. doi:http://dx.doi.org/10.1615/IHTC9.40.

Legendre, D., Magnaudet, J., 1998. The lift force on a spherical bubble in a viscous linear shear flow. J. Fluid Mech. 368, 81–126.

Leoni, N., Le Maitre, O., Rodio, M.G., Congedo, P.M., 2024. Bayesian calibration with adaptive model discrepancy. International Journal for Uncertainty Quantification 14. doi: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. Nuclear Engineering and Design 348, 1–13.

Lubchenko, N., Magolan, B., Sugrue, R., Baglietto, E., 2018. A more fundamental wall lubrication force from turbulent dispersion regularization for
multiphase CFD applications. International Journal of Multiphase Flow 98, 36–44.

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. Chemical Engineering Science 152, 579–585.

Martin, R., 1972. Measurement of the local void fraction at high pressure in a heating channel. Nuclear Science and Engineering 48, 125–138. doi:https://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. Journal of Fluid Mechanics 910, A21. doi:https://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, American Institute of Aeronautics and Astronautics. AlAA 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. Nuclear Engineering and Design 299, 28–36. doi: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. Nuclear Engineering and Design 239, 2365–2371. doi:https://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. Nuclear Engineering and Design 411, 112379. doi:10.1016/j.nucengdes.2023.112379.

Ranz, W.E., Marshall, W., 1952. Evaporation from droplets. Chem. Eng. Prog. 48, 141-146.

Reichardt, H., 1951. Vollständige darstellung der turbulenten geschwindigkeitsverteilung in glatten leitungen. Z. angew. Math. Mech. 31, 208-219. doi:https://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. Journal of Fluid Mechanics 372, 323-355. doi:https://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, 73-93. doi:https://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. Annals of Nuclear Energy 110, 79-87.

Thom, J.R.S., Walker, W., Fallon, T.A., Reising, G.F.S., 1965. Par-
per 6: boiling in sub-cooled water during flow up heated tubes or
annuli, in: Sage UK: London, E.S.P. (Ed.), Proceedings of the in-
stitution of mechanical engineers, conference proceedings, pp. 226-246.
URL: https://www.osti.gov/biblio/4263900, doi:https://doi.org/
10.1016/0301-9322(81)90040-9.

Thomas, R., 1981. Bubble coalescence in turbulent flows. International Journal of Multiphase Flow 7, 709-717.

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 International Journal Series B Fluids and Thermal Engineering.

Tomiyama, A., Tamai, H., Zun, I., Hosokawa, S., 2002. Transverse migration of single bubbles in simple and shear flows. Chemical Engineering Science 57, 1849-1858.

Vlček, D., Sato, Y., 2023. Sensitivity analysis for subcooled flow boiling using eulerian cfd approach. Nuclear Engineering and Design 405, 112194.
doi:https://doi.org/10.1016/j.nucengdes.2023.112194.

Wallis, G.B., 1974. The terminal speed of single drops or bubbles in an infinite medium. International Journal of Multiphase Flow 1, 491-511.

Yao, W., Morel, C., 2004. Volumetric interfacial area prediction in upward bubbly two-phase flow. International Journal of Heat and Mass Transfer 47, 307-328.

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 on Nuclear Reactor Thermal Hydraulics.
Yuan, C., Laurent, F., Fox, R., 2012. An extended quadrature method of moments for population balance equations. Journal of Aerosol Science 51, 1–23.

Zeitoun, O., Shoukri, M., Chatoorgoon, V., 1995. Interfacial heat transfer between steam bubbles and subcooled water in vertical upward flow. J. Heat Transfer 117, 402–407. doi:https://doi.org/10.1115/1.2822536.

Zuber, N., 1964. On the dispersed two-phase flow in the laminar flow regime. Chemical Engineering Science 19, 897–917.

Ü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. International Journal of Heat and Mass Transfer doi:https://doi.org/10.1016/0017-9310(76)90047-8.
$L$ length of the test section

$L_c = \sqrt{\frac{\sigma}{g(\rho_l - \rho_v)}}$ Capillary length

$N_s$ Nucleation site density (m⁻²)

$\mathcal{P}$ Perimeter of a test section

$q_{ki}$ Interfacial heat flux towards phase $k$ (Wm⁻³)

$q_{kw}$ Wall heat flux towards phase $k$ (Wm⁻²)

$q_{SP}$ Single-phase wall heat flux (Wm⁻²)

$q_w$ Total wall heat flux

$q_{w,l \to v}$ Evaporation wall heat

$r_+ = \frac{r}{R_{pipe}}$ Dimensionless radial position

$\mathcal{P}$ Pressure

$\mathcal{P}$ Heated perimeter of the test section

$T_k$ Temperature of phase $k$

$T_s$ Saturation temperature

$T_w$ Wall temperature

$T_{in}$ Inlet temperature

$X$ Thermodynamic quality of the flow

$X_{in}$ Inlet thermodynamic quality

$X_j$ Outlet thermodynamic quality of run $j$

$y$ Distance to wall

$y_+ = \frac{yu_\tau}{\nu_l}$ Dimensionless distance to wall

$\vec{u}_k$ Velocity of phase $k$

$u_{bulk}$ Bulk liquid velocity

$u'_{SP}$ Turbulent fluctuations of a single-phase flow with the same mass flux as the studied flow

$u_\tau$ Friction velocity at the wall

$z$ Axial position along the test tube

*Greek letters.*

$\alpha_k$ Fraction of phase $k$

$\Gamma_k$ Interfacial mass transfer towards phase $k$ (m⁻³s⁻¹)

$\lambda_k$ Thermal conductivity of phase $k$ (Wm⁻¹K⁻¹)

$\mu_k$ Dynamic viscosity of phase $k$ (Pas)

$\nu_k$ Kinetic viscosity of phase $k$ (m²s⁻¹)

$\nu_t$ Turbulent kinetic viscosity of liquid (m²s⁻¹)

$\rho_k$ Volume mass of phase $k$

$\sigma$ Surface tension (Nm⁻¹)
$\tau_{wf}$ Shear stress at the wall (Nm$^{-2}$)

$\omega$ Turbulent dissipation frequency of liquid

$\omega_+$ = $\omega u_\tau / \nu$ Non-dimensional dissipation rate

Dimensionless numbers.

$b = \frac{\alpha_g ||\vec{u}_g - \vec{u}_l||^2}{u_{SP}^2}$ Bubblance parameter

$E_o = \frac{(\rho_l - \rho_g)g d_b^2}{\sigma}$

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

$Nu$ Bubble Nusselt number

$Pe = \frac{d_b ||\vec{u}_g - \vec{u}_l|| \rho C_{pl}}{\lambda_l}$ Peclet number

$Pr = \frac{\nu_l \rho C_{pl}}{\lambda_l}$ liquid Prandtl number

$Re_b = \frac{d_b ||\vec{u}_g - \vec{u}_l||}{\nu_l}$ bubble Reynolds number

$Re = \frac{D_h u_{bulk}}{\nu_l}$ bulk Reynolds number

$We_{drag} = \frac{\rho_l d_b ||\vec{u}_g - \vec{u}_l||^2}{\sigma}$ Drag Weber number

$We_\epsilon = \frac{\rho_l d_b (\epsilon d_b)^{2/3}}{\sigma}$ Turbulent Weber number

$Wo = E_o \frac{\sigma}{||\vec{u}_g - \vec{u}_l||^2}$ Wobble number
# Appendix B. Links between major terms of a two-fluid model

Figure B.19 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_23.png" /></figure>

Figure B.19: 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.