# Momentum and energy budgets of high-pressure subcooled boiling flows using two-fluid RANS simulations

Etienne Chavasse-Fretaz, Alan Burlot, Guillaume Bois, Catherine Colin

## ► To cite this version:

Etienne Chavasse-Fretaz, Alan Burlot, Guillaume Bois, Catherine Colin. Momentum and energy budgets of high-pressure subcooled boiling flows using two-fluid RANS simulations. International Journal of Heat and Mass Transfer, 2025, 244, pp.126839. 10.1016/j.ijheatmasstransfer.2025.126839. cea-05273366

HAL Id: cea-05273366

https://cea.hal.science/cea-05273366v1

Submitted on 22 Sep 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Momentum and energy budgets of high-pressure subcooled boiling flows using two-fluid RANS simulations

Etienne Chavasse-Frétaz<sup>a,b</sup>, Alan Burlot<sup>a,*</sup>, Guillaume Bois<sup>a</sup>, Catherine Colin<sup>b</sup>

## ARTICLE INFO

**Keywords:**

Boiling flows

Multiphase CFD

Interfacial forces

Interfacial heat flux

## ABSTRACT

This study presents the simulation of a high-pressure upward boiling flow in a vertical pipe, using the two-fluid model with Neptune_CFD standard closure laws. Radially-dependent experimental bubble diameters are enforced in the simulations, reducing the number of closure laws and helping to decouple some phenomena in the analysis. The goal of this paper is to highlight the predominant physical phenomena as a function of flow conditions and as a function of the distance to the wall, in order to help determine which model should be improved as a priority. Analysis of the momentum balance for the vapor phase shows that in the radial direction, the predominant effects are turbulent dispersion moving the vapor away from the wall, balanced by the drag and lift forces. An increase in mass flow rate or thermodynamic quality increases liquid shear rate, turbulent viscosity and turbulent dispersion. In addition, liquid enthalpy budget analysis reveals that the liquid temperature is determined by an equilibrium between the radial turbulent diffusion of wall heat, axial inertial thermal effects and an interphase exchange by condensation. Simulations predict significant evaporation near the wall, although this phenomenon has not been observed experimentally, which is a source of error in the mass and enthalpy balances. We demonstrate the improvement in predictions in the near-wall region induced by the use of a limiter to the liquid heat flux in wall Heat Flux Partitioning (HFP) model. Condensation reaches its maximum in an intermediate radial position between the center and the wall, where interfacial area and subcooling are large. This position depends on the thermodynamic quality.

## 1. Introduction

High pressure boiling flows are encountered in many industrial applications such as nuclear reactors or heat exchangers. They are characterized by a strong coupling between thermal and dynamic phenomena, making them difficult to model on a local scale.

The two-fluid Euler–Euler RANS framework has aroused great interest in recent decades to model fluid interactions due to its potential for industrial cases. It requires several closure relations, which is why a wide variety of models has been developed for calculation codes.

Generally, the method adopted for selecting these models is to test them by comparison with experimental data, or even to fit the model parameters (as in [1–3]) without relying on physical arguments. Furthermore, the validation of these models is based solely on void fraction, temperature and velocity profiles, which makes it impossible to analyze the risk of error compensation. However, the DEBORA benchmark [4], which aimed to progress toward a unified set of closures for high-pressure boiling flows, showed that this method is ineffective for selecting a general set of models to accurately reproduce the physical

quantities of interest such as void fraction or liquid temperature under different flow conditions. Moreover, experimental data are not always available for parameter calibration. Other research groups are taking a different approach by building a baseline of coherent mechanistic models, such as the Massachusetts Institute of Technology [5,6] or the Helmholtz-Zentrum Dresden-Rossendorf [7]. Their aim is to develop individual closure laws with better understanding of the underlying physics and to get rid of tuning parameters. However, this approach multiplies the number of submodels, complicates the understanding of the behavior of the global set of models and increases the risk of error compensation. A better understanding of the behavior of solved equations and the interactions between models is therefore needed to guide the modeling work.

In this paper, we propose a different approach by establishing guidelines for the analysis of computational results. To this end, we focus on the modeling of a high-pressure upward flow in a vertical pipe, as in [4]. The DEBORA experimental database [8] covers a wide range of high-pressure flow conditions. This experiment, carried out with a

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

<footnote><sup>b</sup> Institut de Mécanique des Fluides de Toulouse, Université de Toulouse, CNRS, INPT, UPS, Toulouse, 31000, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: alan.burlot@cea.fr (A. Burlot).</footnote>
refrigerant under scale similarity conditions, is representative of the thermal-hydraulic conditions of pressurized water reactors and enables calculations to be compared with an experimental reference. Here, standard models are used to reproduce the experimental trends, preferring to focus on physical phenomena rather than adjusting coefficients in closure laws to better match the experiments.

Because of the high interaction between thermal and dynamic phenomena, we need to reduce the complexity of the calculations to make such an analysis possible. In particular, the bubble size is governed by dynamic phenomena (coalescence, breakup) and thermal effects (evaporation and condensation). However, the DEBORA benchmark has shown that this physical quantity, usually modeled by an interfacial area transport or population balance equations, is among those whose prediction by industrial codes has the greatest dispersion. As in [9,10], radially-dependent experimental bubble diameters are enforced in our simulations (see Section 3.2). This not only reduces the complexity of the calculation by limiting the number of models used, but also eliminates the error in calculating interfacial terms due to poor prediction of bubble diameter.

The goal of this paper is to highlight the predominant physical phenomena as a function of flow conditions and as a function of the distance to the wall, in order to help determining which model should be improved as a priority. The governing equations of the two-fluid model are presented in Section 2, along with the closure laws used. The specifics of the DEBORA experimental database are presented in Section 3. Section 4 describes the numerical setup used for the calculations. A mesh convergence study is also conducted. In Section 5, we analyze the simulations for different flow conditions. We focus specifically on the momentum and enthalpy balances to better understand the interactions between the involved phenomena.

## 2. Governing equations and closure laws

In the one-pressure two-fluid Euler–Euler RANS approach [11,12], a set of conservation equations for mass, momentum and energy is solved for each phase $k \in [l, v]$, where $l$ designates the liquid phase and $v$ the vapor phase:

$$ \frac{\partial}{\partial t} (\alpha_k \bar{\rho}_k) + \nabla \cdot (\alpha_k \bar{\rho}_k \mathbf{u}_k) = \Gamma_k \quad (1a) $$

$$ \frac{\partial}{\partial t} (\alpha_k \bar{\rho}_k \hat{\mathbf{u}}_k) + \nabla \cdot (\alpha_k \bar{\rho}_k \hat{\mathbf{u}}_k \otimes \hat{\mathbf{u}}_k) = -\alpha_k \nabla \bar{\rho} + \nabla \cdot \left[ \alpha_k (\bar{\mathcal{T}} + \mathcal{T}^T) \right] + \alpha_k \bar{\rho}_k \bar{\mathbf{g}} + \Gamma_k \hat{\mathbf{u}}_\sigma + \mathbf{M}_k \quad (1b) $$

$$ \frac{\partial}{\partial t} (\alpha_k \bar{\rho}_k \hat{H}_k) + \nabla \cdot (\alpha_k \bar{\rho}_k \hat{H}_k \hat{\mathbf{u}}_k) = \alpha_k \frac{\partial \bar{\rho}}{\partial t} + \alpha_k \bar{\rho}_k \bar{\mathbf{g}} \cdot \hat{\mathbf{u}}_k + \nabla \cdot \left[ \alpha_k \left( \frac{\lambda_k}{C_{p,l}} + \frac{\mu_k^T}{P r^T} \right) \nabla \hat{H}_k \right] + \Pi'_k + \Gamma_k \hat{H}_\sigma \quad (1c) $$

In these equations, $u, T, \rho, \mu, \lambda, C_p$ denote respectively the velocity, temperature, density, dynamic viscosity, conductivity and thermal capacity at constant pressure. $H_k$ is the total enthalpy of phase $k$, related to enthalpy $h_k$ by $H_k = h_k + u_k^2/2$. $\alpha_v$ (or $\alpha$) is the void fraction, $g$ is the gravitational acceleration, $\hat{u}_\sigma$ and $\hat{H}_\sigma$ are respectively the velocity and total enthalpy of the interface, $p$ corresponds to the pressure assumed locally equal for both phases. The Favre average of a quantity $\phi$ is denoted $\hat{\phi}$, the phase average is denoted $\hat{\phi}$.

In Eqs. (1a), (1b) and (1c), the source term $\Gamma_k$ represents the mass transfer from phase $k$ to the interface. In Eq. (1b), $\bar{\mathcal{T}}$ is the viscous stress tensor and $\mathcal{T}^T$ is the Reynolds stress tensor defined by $\mathcal{T}^T = -\rho_k \mathbf{u}_k' \otimes \mathbf{u}_k'$. The velocity of the interface $\hat{u}_\sigma$ is arbitrarily taken equal to the vapor velocity.

$\mathbf{M}_k$ represents interfacial momentum transfer. It is modeled as follows:

$$ \mathbf{M}_k = \mathbf{M}_k^{\text{drag}} + \mathbf{M}_k^{\text{lift}} + \mathbf{M}_k^{\text{AM}} + \mathbf{M}_k^{\text{TD}} + \mathbf{M}_k^{\text{WL}} \quad (2) $$

where the right-hand side terms represent the drag, lift, added mass, turbulent dispersion and wall lubrication effects respectively.

In particular, the drag contribution writes:

$$ \mathbf{M}_v^{\text{drag}} = -\mathbf{M}_1^{\text{drag}} = -\frac{1}{8} a_l \bar{\rho}_l C_D |\hat{\mathbf{u}}_r| |\hat{\mathbf{u}}_r| \quad (3) $$

with $\hat{\mathbf{u}}_r = \hat{\mathbf{u}}_v - \hat{\mathbf{u}}_l$. The drag coefficient $C_D$ is computed according to the Ishii and Zuber model [13].

$a_l$ is the interfacial area concentration:

$$ a_l = \frac{6\alpha_v}{D_s} \quad (4) $$

with $D_s$ the Sauter mean diameter.

The lift contribution writes:

$$ \mathbf{M}_v^{\text{lift}} = -\mathbf{M}_1^{\text{lift}} = -C_L \alpha_v \bar{\rho}_l \hat{\mathbf{u}}_r \wedge (\nabla \wedge \hat{\mathbf{u}}_l) \quad (5) $$

The lift coefficient $C_L$ is computed with an adaptation of the Tomiyama et al. model [14] as described in [4].

The turbulent dispersion contribution writes:

$$ \mathbf{M}_v^{\text{TD}} = -\frac{2}{3} C_{TD} \bar{\rho}_l k_l \nabla \alpha_v \quad (6) $$

where $k_l$ is the turbulent kinetic energy of the liquid phase, defined as $k_l = \text{tr}(\mathcal{T}^T)$. $C_{TD}$ is modeled by the Generalized Turbulent Dispersion model (GTD) [15,16]. Other closure models for the interfacial momentum transfer are described in Table 1.

A second order turbulence model based on the $R_{ij} - \epsilon$ Elliptic-Blending Reynolds Stress Model (EBRSM) [17] is used in the liquid phase with interfacial source terms for Reynolds stresses and dissipation. Turbulence is neglected in the vapor phase.

In Eq. (1c), the pressure contribution $\alpha_k \frac{\partial \bar{\rho}}{\partial t}$ is zero in the steady-state conditions of interest. The term $\nabla \cdot (\alpha_k \hat{\mathbf{u}}_k \cdot \bar{\mathcal{T}})$ representing the work of viscous forces can be neglected in comparison with heat transfers for low-viscosity fluids, as in [11]. In the same way, the gravity contribution $\alpha_k \bar{\rho}_k \bar{\mathbf{g}} \cdot \hat{\mathbf{u}}_k$ can be neglected for the flow conditions of the DEBORA experiment. In addition, the mechanical contribution to total enthalpy can be neglected, so $H_k \approx h_k$.

As described in [18] and neglecting the term $\Gamma_k (\hat{H}_\sigma - \hat{H}_l)$, the enthalpy budget (1c) for the liquid phase can be rewritten as follows:

$$ \alpha_l \frac{\partial \hat{H}_l}{\partial t} = \underbrace{-\alpha_l \bar{\rho}_l \hat{\mathbf{u}}_l \cdot \nabla \hat{H}_l}_{\text{inertia}} + \underbrace{\nabla \cdot \left[ \alpha_l \left( \frac{\lambda_l}{C_{p,l}} + \frac{\mu_l^T}{P r^T} \right) \nabla \hat{H}_l \right]}_{\text{diffusion}} + \underbrace{\Pi'_l}_{\text{interfacial heat transfer}} \quad (7) $$

$\mu_l^T$ denotes the turbulent dynamic viscosity and $Pr^T$ is the turbulent Prandtl number. The interfacial heat transfer $\Pi'_l$ is modeled by a correlation for the Nusselt number $Nu$ that depends on the bubble Reynolds number $Re_b$, the liquid Prandtl number $Pr_l$ and the Jakob number $Ja$:

$$ \Pi'_l = a_l Nu \frac{\lambda_l}{D_s} (T_{sat} - T_l) \quad (8) $$

with $T_{sat}$ the saturation temperature.

The wall Heat Flux Partitioning (HFP) is computed by the Mimouni et al. model described in [19,20]. It is based on the classical decomposition of the wall heat flux into three contributions (convection heat flux in the liquid $\phi_c$, quenching flux $\phi_q$ and evaporation flux $\phi_e$) introduced by Kurul and Podowski [21]. The HFP model calculates both these fluxes and the wall temperature $T_w$ using a fixed-point algorithm. The departure diameter $D_d$ is adapted from Ünal correlation [22]:

$$ D_d = 2.42 \times 10^{-5} p^{0.709} \frac{a}{\sqrt{b \left[ \max(u_l/0.61, 1) \right]^{0.47}}} \quad (9) $$

$$ a = \frac{\max(T_w - T_{sat}, 0)}{2 \rho_v L^{\text{vap}} \sqrt{\pi \alpha_w}}; \quad (10) $$
$$b = \frac{\rho_l \Delta T_{sub}}{2 (\rho_l - \rho_v)}; \qquad (11)$$

$$\Delta T_{sub} = \begin{cases} \max(T_{sat} - T_l, 1.0 \times 10^{-8}) & \text{if } 0 \le St \le St_{lim}; \\ q_{w}/(\rho_l c_{p,l} u_l S t_{lim}) & \text{otherwise.} \end{cases} \qquad (12)$$

where the wall diffusivity is $\alpha_w = \lambda_w / (\rho_w c_{p,w})$, $St = q_w / (\rho_l c_{p,l} u_l \Delta T_{sub})$, $St_{lim} = 0.0065$. Correlations from Cole [23] and Kurul and Podowski [21] are respectively used for the bubble departure frequency $\eta$ and nucleation site density $N_b''$:

$$\eta = \sqrt{\frac{4g}{3D_d} \left(1 - \frac{\rho_v}{\rho_l}\right)} \qquad (13)$$

$$N_b'' = (210 (T_w - T_{sat}))^{1.8} \qquad (14)$$

Eqs. (9), (10), (11), (12), (13) and (14) show that these submodels are strongly coupled to each other. Consequently, an erroneous prediction of a quantity can have a significant impact on the prediction of the others, and ultimately on the HFP. Error compensation is also possible. However, an in-depth study of these submodels is beyond the scope of this paper. More details on these closures are available in [24–28].

Additionally, as in [29], a Maximum Oversaturation Control (MOC) parameter is used to control the liquid overheating in the wall cell. It acts as a correction factor that diminishes the part of the wall heat flux that goes to the liquid phase. Indeed, at each HFP computation iteration, if $T_l > T_{sat}$, $\phi_c$ and $\phi_q$ are multiplied by a coefficient

$$C = \max \left( \frac{1 + \cos \pi r}{2}, 0.05 \right) \text{ with } r = \frac{T_l - T_{sat}}{\text{MOC}} \qquad (15)$$

In the following, MOC is taken equal to 0.1 °C or 1 °C; its value is discussed in Section 5.1.3. The liquid temperature $T_l$ in the first wall cell is calculated following a Van Driest law [30]. A wall-law accounting for bubble-roughness described in [31] is used to compute the liquid velocity tangential to the wall $u_{lz}$ in the first wall cell as follows:

$$\frac{u_{lz}}{u_{\tau}} = u^{+} = \frac{1}{\kappa} \ln (y^{+}) + B - \Delta u^{+} \qquad (16)$$

The last term corrects the single-phase liquid velocity by taking into account the bubbles size and density at the wall:

$$\Delta u^{+} = \begin{cases} 0 & \text{if } k_{r}^{+} \le 11.3 \\ \frac{1}{\kappa} \ln (1 + C_{kr} k_{r}^{+}) & \text{if } k_{r}^{+} > 11.3 \end{cases} \qquad (17)$$

where $\kappa = 0.41$, $B = 5.3$, $C_{kr} = 0.5$, $k_{r}^{+} = \rho_{l} k_{r} \sqrt{u_{\tau} u^{*}} / \mu_{l}$, $k_{r} = \alpha D_{S}$, $u_{\tau} = \sqrt{\tau_{w}/v_{l}}$, $u^{*} = C_{\mu}^{1/4} \sqrt{\bar{k}_{l}}$, $C_{\mu} = 0.09$.

The wall HFP model uses fluid properties at a distance of one bubble diameter from the wall.

The thermodynamic properties of the fluids are taken from the CATHARE tables.

The Neptune_CFD code [32] is a numerical simulation code for multiphase flows developed by EDF, CEA, Framatome and the Institut de Radioprotection et de Sûreté Nucléaire (IRSN), mainly for nuclear applications. It uses finite volume spatial discretization for unstructured meshes. In this study, version 7.2.0 is used with the models recommended in the Neptune_CFD Verification and Validation (VnV) note [29], widely used for this type of flow (e.g. [33,34]), with a few specifics:

* The $R_{ij} - \epsilon$ EBRSM turbulence model [17] is used because it provides better mesh convergence than the $R_{ij} - \epsilon$ SSG model. Indeed, the cell size of the finer mesh in the direction normal to the wall is too small for the SSG model ($y_{w}^{+} = 30$, see Section 4.1).

* No Interfacial Area Transport Equation (IATE) is used (see Section 3.2).

* No wall-lubrication force is used as no significant difference is observed by deactivating Tomiyama wall-lubrication force [35].

Table 1 presents a summary of the models used for the calculations.

<figure><img src="image_3.png" /><figcaption>Fig. 1. General view of the test section.</figcaption></figure>

Source: Adapted from [8].

# 3. Experimental database

## 3.1. Presentation of DEBORA experiment

DEBORA is a two-phase thermal-hydraulic experimental facility designed by CEA and EDF in the 1990s. A complete description of the experiment can be found in [8]. It consists in a high-pressure Freon (R12) upward flow in a vertical pipe and covers a wide range of inlet temperatures, pressures (14 to 30 bar), flow rates and wall heat fluxes. It is designed to represent pressurized water reactor conditions (100 to 180 bar) via a scale similarity. The specificity of these experimental measurements is that they provide radial profiles of void fraction, liquid temperature, bubble diameter and vapor velocity rather than only a cross-sectional average.

A diagram of the test section is shown on Fig. 1. The pipe internal diameter is 19.2 mm, and its wall is electrically heated over a length of 3.5 m. All measurements are carried out at $Z = 3.485$ m, below the end of the heated length.

Each test in the DEBORA experiment is identified according to the physical conditions of the flow, and noted as follows $C_{xx}G_{yy}P_{zz}W_{aa}Te_{bb}$ with:

* $C_{xx}$ the experimental campaign;

* $G_{yy}$ the fluid mass flux in $10^3 \times \text{kg m}^{-2} \text{s}^{-1}$;

* $P_{zz}$ the pressure in bar;

* $W_{aa}$ the wall heat power in kW;

* $Te_{bb}$ the inlet temperature in °C.

For each set of experimental flow conditions ($G, P, W$), a series of measurements is carried out at different inlet temperatures $T_l$, with a degree of subcooling between 15 and 30 °C. To give an order of magnitude, at 26 bar $T_{sat} = 87$ °C.

Campaign C29 (see [4]) provides information on void fraction, Sauter diameter and interfacial area profiles. A mono-optical probe was used for the measurements and a model was assumed for the vapor velocity. Campaign C8 provides liquid temperature measurements.
<table><caption>**Table 1**
Models used in the simulations.</caption><thead><tr><th>Phenomenon or quantity</th><th>Model</th></tr></thead><tbody><tr><td>Liquid phase turbulence</td><td>R<sub>ij</sub> - ε EBRSM [17] with turbulent reverse coupling</td></tr><tr><td>Drag force</td><td>Ishii and Zuber [13]</td></tr><tr><td>Added mass force</td><td>Zuber [36]</td></tr><tr><td>Lift force</td><td>Tomiyama et al. [14]</td></tr><tr><td>Wall lubrication force</td><td>-</td></tr><tr><td>Turbulent dispersion force</td><td>GTD [15,16]</td></tr><tr><td>Wall HFP</td><td>Modified Mimouni et al. [19], Kurul and Podowski [21] with MOC = 0.1 °C</td></tr><tr><td>D<sub>S</sub></td><td>Experimental measurements interpolation (see Section 3.2)</td></tr><tr><td>Condensation</td><td>Manon [37] Nu = 2 + 0.6Re<sub>b</sub><sup>1/2</sup>Pr<sub>l</sub><sup>1/3</sup></td></tr><tr><td>Evaporation</td><td>Berne [38] Nu = max(2, <sup>12</sup>/<sub>π</sub>Ja, <sup>2</sup>/<sub>√ε</sub>(Re<sub>b</sub>Pr<sub>l</sub>)<sup>1/2</sup>)</td></tr><tr><td>Liquid turbulent heat transfers</td><td>Standard Gradient Diffusion Hypothesis (SGDH)</td></tr><tr><td>Pr<sup>T</sup></td><td>1.03</td></tr></tbody></table>

## 3.2. Test tube hypothesis

The models presented in Section 2 strongly depend on the bubble diameter $D_S$. Consequently, an erroneous diameter prediction will have a strong impact on the prediction of the void fraction, velocities and liquid temperature. This is why the experimental measurements are used in this study to avoid dependence on a diameter prediction model.

Similarly as in [6,10], we make the assumption that reducing the inlet temperature is equivalent to measuring physical quantities at a lower height in the test section, all other quantities being equal. This allows us to interpolate a two-dimensional field of diameters from experimental profiles and the corresponding axial abscissa. The experimental measurements have a radial resolution of approximately 0.2 mm. The data is smoothed using hyperbolic tangent functions (see Fig. 17(b) and Fig. 24(b)) that enable extrapolation of the measurements near the wall. In the axial direction, a linear interpolation is done between the corresponding abscissa. A linear extrapolation to $D_S = 0$ at $z = 0$ m is performed upstream of the first profile. This field is then enforced in the simulations.

This method not only reduces the complexity of the calculation by limiting the number of models used, but also eliminates the error in calculating interfacial terms due to poor prediction of bubble diameter. It is worth noting that the condensation model described in Eq. (8) uses an interfacial heat transfer proportional to the interfacial area defined by Eq. (4). Hence, if an incorrect void fraction is predicted, the interfacial heat transfer would also be wrong even if the experimental diameter is enforced. The same observation is true for mass and momentum transfers. Enforcing the experimental interfacial area instead of the diameter would have had the same benefit in terms of reducing the complexity of the model. However, a majority of the submodels used have been developed for single bubbles, and therefore depend directly on diameter.

In the following, according to [10], we will call a *test tube* a series of experimental runs that have the same outlet pressure, mass flow rate and heating power, and in which the inlet temperature $T_l$ is varied. We selected the test tubes G3P26W23 and G2P26W16 of C29 campaign because they cover a wide range of physical conditions (see Fig. 2). As explained in [10,39], at each inlet temperature corresponds an outlet thermodynamic quality $x_{eq,out}$. A model taken from [39] of the flow quality $x$ is plotted for test tubes G3P26W23 and G2P26W16 as a function of cross-section average void fraction $\langle\alpha\rangle = 1/S \iint_S \alpha ds$ alongside the experimental values of the equilibrium thermodynamic quality $x_{eq}$. This model depends on the phase velocity slip ratio $\langle u_v \rangle_v / \langle u_l \rangle_l$, with $\langle u_k \rangle_k = \langle u_k \rangle / \langle \alpha_k \rangle$ and $\langle u_k \rangle = 1/S \iint_S \alpha_k u_k ds$ (see e.g. [40]). Flow conditions range from quasi-single-phase to an average void fraction of 0.3, which corresponds to a large proportion of the industrial conditions addressed by the two-fluid Euler–Euler RANS model. The equilibrium thermodynamic quality should be slightly lower than the flow quality, as the liquid remains subcooled at the core of the flow. For $\langle\alpha\rangle \le 5\%$, the flow is globally subcooled. For high void fractions, it is necessary to assume a significant slip ratio $u_v/u_l$ to verify the inequality $x_{eq} \le x$. This

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Physical conditions of the considered test tubes.</figcaption></figure>

ratio results from the combination of local relative velocity and void fraction distribution. Furthermore, for the same mean void fraction $\langle\alpha\rangle$, the equilibrium thermodynamic quality of the high flow rate case G3P26W23 tend to be higher than the equilibrium thermodynamic quality corresponding to test tube G2P26W16, meaning that the liquid is closer to saturation. This shows that a higher mass flow rate brings the flow closer to thermal equilibrium due to a turbulent mixing effect. This effect will be discussed at the local scale in Section 5.2.

When available, we plot the experimental measurements along with the calculations to ensure that the simulation gives the right orders of magnitude for the variables of interest $(\alpha, T_l, u_v)$.

## 3.3. Experimental results

In this section, we discuss some of the experimental results of the DEBORA experiment for the considered conditions. The experimental measurements available for the test tube G2P26W16 are displayed on Fig. 3. The radial profiles are plotted versus $r^+ = r/R$, between the pipe axis ($r^+ = 0$) and the heated wall ($r^+ = 1$). Each color corresponds to a different inlet temperature $T_l$.

By definition, the void fraction $\alpha$ increases with thermodynamic quality, and the zone where $\alpha > 0$ expands (Fig. 3(a)). For the highest inlet temperature, the void fraction peak is further away from the wall.

The experimental diameters used in the calculations are shown on Fig. 3(b). Near the wall, coalescence is the predominant phenomenon. The decrease of the diameter away from the wall in the subcooled zone is due to condensation. As inlet temperatures $T_l$ increase, condensation tends to be lower and diameters globally increase.
<figure><img src="image_5.png" /><figcaption>Fig. 3. Experimental data for test tube 29G2P26W16 (a, b) and 8G2P26W16 (c).</figcaption></figure>

The liquid temperature $T_l$ increases with the inlet temperature, until reaching $T_{sat}$ close to the wall for the highest inlet temperature (Fig. 3(c)). However the liquid temperature does not exceed $T_{sat}$. The core of the flow remains subcooled due to the high mass flow rate of liquid.

The measurement uncertainty is ±1% for void fraction and ±0.2 °C for liquid temperature [8]. The high radial resolution and accuracy of these measurements provide an accurate database for comparing simulation results.

# 4. Numerical setup

In this section, we describe the numerical configuration used for the calculations.

## 4.1. Meshing and boundary conditions

The cylindrical geometry of the DEBORA experimental device (see Section 3) allows us to model only an angular fraction of this cylinder with an angle of 0.1 radians (such as in [4]) and an internal radius $R = 9.6$ mm (see Fig. 4). We use symmetry boundary conditions on the vertical faces that delimit this angle and at the central axis. The right boundary at $r = R$ is a solid wall where a constant heat flux $q''$ is imposed in the heated section (see Section 2 for details on the HFP model). The inlet boundary condition at $z = 0$ is defined by the mass flow rate $G$ and inlet uniform temperature $T_i$. The inlet turbulent velocity profile is established through a recycling boundary condition. A constant pressure $P$ is set at the outlet ($z = 4$ m). Standard gravity $g = -9.81 \text{ m/s}^2$ is applied along $z$.

<figure><img src="image_6.png" /><figcaption>Fig. 4. Description of the computational domain.
Source: Adapted from [41]</figcaption></figure>

Table 2 shows the characteristics of the different meshes used. Cell size is regular along both axes. The calculations are initially conducted on the coarsest mesh, denoted as m2, until a steady state is established, which initializes a new calculation performed on the finer mesh m4.
<table><caption>**Table 2**
Mesh details.</caption><thead><tr><th>Mesh name</th><th>Radial cell number</th><th>Axial cell number</th><th>Total cell number</th></tr></thead><tbody><tr><td>m2</td><td>20</td><td>200</td><td>4000</td></tr><tr><td>m4</td><td>40</td><td>400</td><td>16 000</td></tr><tr><td>m8</td><td>80</td><td>800</td><td>64 000</td></tr></tbody></table>

In the same way, the steady state reached with m4 initializes the calculation on m8.

The Neptune_CFD steady algorithm is used for time advancement of the solution with a velocity prediction algorithm coupling mean and relative velocities.

## 4.2. Mesh convergence

The results obtained with the meshes m2, m4 and m8 are presented in this section for the case G3P26W23Te69. Here, the Maximum Over-saturation Control (MOC) is taken equal to 1°C. Physical quantities of interest are shown on Fig. 5. The mesh convergence is satisfactory. The Reynolds stress $R_{zz}$ shows some discrepancies near the wall, but it does not show a significant impact on the convergence of the velocities. However, the lower turbulent diffusivity $\mu_i^T / Pr^T$ for m8 (see Fig. 5(f)) slightly reduces the radial enthalpy diffusion and increases the liquid temperature gradient, which leads to imperfect convergence of the liquid temperature at the core (see Fig. 5(b)).

Fig. 6 shows radial profiles of the different terms of momentum balance (1b) for the vapor phase with the meshes m2, m4 and m8. These budgets will be analyzed in detail in Section 5. The axial components are plotted on Fig. 6(a) and the radial components on Fig. 6(b). The mesh convergence is excellent in the core of the pipe (for $r^+ < 0.9$). Near the wall, the lift and turbulent dispersion terms $M_v^{\text{lift}}$ and $M_v^{\text{TD}}$ do not converge well. Here, the radial size of a m4 cell is 0.24 mm, which is the order of a bubble diameter (see Fig. 3(b)); for m8 the radial cell size is the order of a bubble radius. However, the statistical approach used by these models is flawed for small cells relative to bubble size. A wall-damping function (e.g. [42]) or a Lubchenko regularization function [43] could be used to address the divergence of these forces near the wall. Such functions are not yet available in the Neptune_CFD code. Models better adapted to near-wall averaging are also developed in [44,45].

In addition, the phase change term $\Gamma_v \hat{u}_\sigma$ diverges and gets predominant in the vicinity of the wall due to the high evaporation phenomenon. Indeed, for the set of models used, the simulated cases systematically present superheated liquid in the region close to the wall. This discrepancy between simulations and experimental data in relation to the MOC parameter will be discussed in Section 5.1.3.

To conclude, a mesh convergence study has been carried out, proving the equivalence of the results for the m4 and m8 meshes for the model set used (see Section 2). The results presented below are obtained with the m4 mesh, due to its lower computational cost.

## 5. Physical analysis of flows using momentum and enthalpy balances

In this section, we highlight the predominant physical phenomena as a function of flow conditions (thermodynamic quality, mass flow rate), and as a function of the distance to the heated wall (near-wall region versus core of the flow). This study enables us to compare the relative weight of different contributions in order to determine the cause of discrepancies between simulation and experiment. We focus on the impact of the thermodynamic quality and mass flow rate because the flows that we aim to accurately reproduce are characterized by wide ranges of variation in these quantities. In addition, this is how the

DEBORA experiment was designed so it enables a comparison between the simulation and the experimental data.

The selected cases are presented in Table 3. They are chosen for comparison with similar flow conditions except for one parameter. The reference case is G2P26W16Te68. The effect of a higher mass flow rate is studied by comparison with G3P26W23Te69. The comparison of the reference case with G2P26W16Te64 and G2P26W16Te58 will show the effect of the outlet thermodynamic quality.

We only selected cases for which the void fraction distribution is wall-peaked. Indeed, as it has been previously demonstrated in the DEBORA benchmark [4], the Tomiyama et al. model [14] fails to predict the lift inversion required to obtain core-peaked void fraction distributions.

In this section, MOC is taken equal to 0.1°C, as justified in Section 5.1.3.

### 5.1. Reference case: G2P26W16Te68

Here, we analyze the reference case G2P26W16Te68 by comparison to the experimental data. The quantities relative to the vapor phase are plotted only if $\alpha \ge 0.01$. Fig. 7 shows the radial profiles of void fraction and liquid temperature for the entire test tube G2P26W16, i.e. for each inlet temperature $T_i$. The general trend of the void fraction is fairly well predicted, however its slope is too steep for most profiles. The average void fraction on the pipe cross-section is highly overestimated, as is its value at the wall, by a factor of around 30% for the highest inlet temperatures (see Table 4). For the reference case G2P26W16Te68, the void fraction is underestimated in the center ($r^+ < 0.2$). On the following figures, the saturated zone is delimited by diamond markers on the simulated profiles. For this case, the liquid is superheated in the region $r^+ > 0.9$, although the experimental temperature does not exceed $T_{sat}$. The temperature at the core is underestimated by about 1°C. The liquid temperature profiles for the test tube flatten out as the inlet temperature increases, for both measurements and simulation. These axial trends can also be observed on the 2D plots Fig. 8.

#### 5.1.1. Dynamic effects

The vapor axial velocity $\hat{u}_{v,z}$ is plotted on Fig. 9(a) alongside the experimental model of [8] based on several hypotheses, including an hypothesis of homogeneous equilibrium. These assumptions are discussed in [46]; in particular, this profile should be flat at $r^+ = 0$. The code predicts a higher axial vapor velocity compared to the experimental model, confirming the necessity of the slip ratio assumption as demonstrated in Section 3.3. The vapor velocity is slightly higher than the liquid velocity (Fig. 9(b)) due to the buoyancy effect.

Fig. 10 shows the radial and axial relative velocities for the reference case. At $r^+ \approx 0.1$, the relative axial velocity $\hat{u}_{r,z} = \hat{u}_{v,z} - \hat{u}_{l,z}$ is an order of magnitude higher than the relative radial velocity $\hat{u}_{r,r} = \hat{u}_{v,r} - \hat{u}_{l,r}$. $\hat{u}_{r,z}$ is around 3 times higher than $\hat{u}_{l,r}$ in the rest of the flow. The negative value of $\hat{u}_{r,r}$ shows that the vapor produced at the wall or in the overheated region globally goes toward the core, while the liquid radial velocity $\hat{u}_{l,r}$ absolute value is generally low over the entire cross-section with respect to the vapor radial velocity $\hat{u}_{v,r}$ due to the density difference between phases.

Axial and radial momentum balances for the vapor phase (Eq. (1b) with $k = v$) are presented on Fig. 11. In the axial direction, we find the classic equilibrium between the drag and pressure gradient contributions. The convection and gravity terms are not negligible. The gravity contribution would be negligible for a flow with a higher phase density ratio $\bar{\rho}_l / \bar{\rho}_v$ which occurs for lower pressure flows, such as the 14 bar DEBORA tests. The lift force also plays a role in the axial balance, as the radial and axial components of the relative velocity are close in magnitude. As mentioned in Section 4.2, the phase change term $\Gamma_v \hat{u}_\sigma$ is predominant in the vicinity of the wall due to the high evaporation phenomenon. However, the physical relevance of this term is questionable because there is no experimental overheating.
<figure><img src="image_7.png" /><figcaption>Fig. 5. Mesh convergence of physical quantities – G3P26W23Te69.</figcaption></figure>

<table><caption>Table 3
Flow conditions of selected cases.</caption><thead><tr><th>Case</th><th>Outlet quality x<sub>eq,out</sub> [ ]</th><th>Outlet pressure P [bar]</th><th>Mass flow rate G [kg/m<sup>2</sup>/s]</th><th>Wall heat flux q'' [kW/m<sup>2</sup>]</th><th>Inlet temperature T<sub>i</sub> [°C]</th></tr></thead><tbody><tr><td>G3P26W23Te69</td><td>0.06</td><td>26.17</td><td>2981.8</td><td>107.52</td><td>69.72</td></tr><tr><td>G2P26W16Te68</td><td>0.06</td><td>26.15</td><td>1986.0</td><td>73.89</td><td>68.52</td></tr><tr><td>G2P26W16Te64</td><td>0.01</td><td>26.15</td><td>1991.8</td><td>73.89</td><td>64.64</td></tr><tr><td>G2P26W16Te58</td><td>-0.07</td><td>26.15</td><td>2012.1</td><td>73.89</td><td>58.80</td></tr></tbody></table>

In the radial direction, with our set of models, the predominant terms are drag, lift and turbulent dispersion. The first two push the vapor phase toward the wall ($M_v^{drag} \ge 0$, $M_v^{lift} \ge 0$), while turbulent dispersion tends to reduce the gradient of the void fraction by pushing the vapor away from the wall ($M_v^{TD} \le 0$). Indeed, the latter force is proportional to the opposite of the void fraction gradient and to

the turbulent kinetic energy (see Eq. (6)). The sign of the lift force is due to the positive lift coefficient $C_L$ computed by the Tomiyama et al. model [14]. This force is especially important in the near-wall region due to the high shear rate of the liquid. It is worth noting that the lift model is developed for single bubbles so its prediction in these flow conditions may be rather unreliable. The GTD turbulent dispersion
<figure><img src="image_8.png" /><figcaption>Fig. 6. Vapor momentum balances – G3P26W23Te69.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Physical quantities for test tube G2P26W16.</figcaption></figure>

model [15,16] is derived through theoretical analysis and the order of magnitude of this effect is difficult to estimate a priori. A poor prediction of the lift force would lead to a compensation phenomenon by drag and turbulent dispersion that could induce a deformation of the void fraction profile. The role of the pressure gradient is of secondary importance in the radial direction.

### 5.1.2. Thermal effects

The interfacial heat transfer $\Pi'_i$ is shown on Fig. 12(a). As described in Eq. (8), this source term of the enthalpy budget is proportional to the interfacial area $a_i$ and to the degree of subcooling $\Delta T_i = T_{sat} - T_i$. Hence, $\Pi'_i > 0$ when the liquid subcooling constrains the vapor to condense, which means that the liquid is heated by the condensation. $\Pi'_i < 0$ reflects the evaporation phenomenon, i.e. when the superheated liquid evaporates. The condensation reaches a maximum for $r^+ = 0.4$. For $r^+ < 0.4$, the interfacial area tends to zero because the void fraction does; for $r^+ > 0.4$, $\Delta T_i$ decreases. The liquid radial velocity $\hat{u}_{l,r}$ reverses at $r^+ = 0.2$ because the vapor disappears and the liquid takes its place (see Fig. 12(b)). Thus, the liquid tends to go toward an area of high condensation due to mass conservation.

As discussed in Section 4.2, the evaporation phenomenon is intense near the wall ($r^+ < 0.9$), as the liquid temperature reaches 0.1 °C of superheating. Actually, the radial transport of the heat in the liquid does not compensate for the energy input at the wall, so the heat accumulates in the liquid, which reaches $T_i > T_{sat}$.

The enthalpy budget (7) is shown on Fig. 13(a). The inertial heat transfer contribution is proportional to the product $\hat{u}_l \cdot \nabla H_l$, which is determined by the axial gradients in the core, because the radial enthalpy gradient tends to zero in this region:

$$ \hat{u}_l \cdot \nabla H_l = \hat{u}_{l,z} \times \frac{\partial H_l}{\partial z} + \underbrace{\hat{u}_{l,r} \times \frac{\partial H_l}{\partial r}}_{\text{negligible at } r^+ = 0} \quad (18) $$

This effect illustrates the coupling between the dynamic and thermal effects and the coupling between axial and radial evolutions. In the superheated zone, the liquid velocity tends to zero at the wall, which leads to a negligible inertial term. On the other hand, the high enthalpy gradient due to the evaporation effect causes the diffusion term to be significant and it balances the evaporation term.

As mentioned in Section 4.2, each term of the enthalpy budget has been scaled by a factor $V_{pipe}/(S_{wall} \times q'') = R/2q''$. Consequently, the
<figure><img src="image_10.png" /><figcaption>Fig. 8. 2D fields of physical quantities for test tube G2P26W16.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Axial velocities – G2P26W16Te68. The experimental vapor velocity is reconstructed with a model [8]. The quantities relative to the vapor phase are plotted only if $\alpha \ge 0.01$.</figcaption></figure>

value of the diffusion term for $r^+ < 0.1$ is close to 1, which means that in this single-phase part of the flow, all the power of the wall is transmitted to the liquid in the core by heat diffusion. This result is coherent with the direct numerical simulations conducted in [47]. This heat input is compensated for by the convection effect, which removes heat from this zone by supplying subcooled liquid along the axial direction.

Fig. 13(b) shows the enthalpy balance integrated in the radial direction. Integral terms $\Phi$ are computed for each term $\phi$ of the budget as:

$$ \Phi(r^+) = \frac{1}{Rq''} \int_0^{r^+} r'\phi(r')dr' \qquad (19) $$

Hence $\Phi(r^+ = 1)$ is the cross-section average of $\phi$. In the whole cross-section, the evaporation phenomenon is small compared with the condensation. However, the balance between these two effects in the near-wall region is largely due to the HFP model. This is why the

sensitivity of one of its main parameters, the MOC, is studied in the next section.

### 5.1.3. MOC parameter sensitivity study

Here, we discuss the value chosen for the MOC parameter defined Eq. (15), in order to address the liquid overheating phenomenon. A comparison between simulations of case G2P26W16Te68 carried out with three values of this parameter is shown on Fig. 14. For stability reasons, this parameter cannot be zero, the reference case uses MOC = 0.1 °C. We compare it to the standard value MOC = 1 °C recommended in [29] and to an infinite MOC, which corresponds to the Mimouni et al. HFP model [19] without any correction.

As expected, the liquid overheating reduces from 0.5 °C to 0.1 °C, and the region where $T_l > T_{sat}$ shrinks. The rest of the flow remains unchanged. The impact on the void fraction and on the other quantities is negligible.
<figure><img src="image_12.png" /><figcaption>Fig. 10. Relative velocities – G2P26W16Te68.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Vapor momentum balances – G2P26W16Te68.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Liquid quantities showing the inversion of the radial liquid velocity in relation to the interfacial heat transfer – G2P26W16Te68.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 13. Radial enthalpy balance – G2P26W16Te68.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Physical quantities – G2P26W16Te68.</figcaption></figure>

Fig. 15(a) shows the axial profile of the distribution of total wall heat flux $\phi_l$ split into the three contributions $\phi_c$, $\phi_e$ and $\phi_q$ (see Section 2) for the test tube G2P26W16, where $z^+ = z/Z$ and $Z = 3.485$ m denotes the abscissa of the measurement plane. At approximately $z^+ = 0.6$, the liquid temperature near the wall exceeds the saturation temperature (vertical line on Fig. 15(a)) so the MOC model is activated. A low superheating tolerance MOC = 0.1 °C leads to greater flux correction. Indeed, the liquid fluxes $\phi_c$ and $\phi_q$ are reduced while, symmetrically, $\phi_e$ is increased. This correction results in a total evaporation flux that is 20% higher at the end of the heated section (see Fig. 15(b)).

There is almost no change in the momentum balance (not shown).

As shown on Fig. 16, in line with expectations, the evaporation contribution to the enthalpy balance strongly diminishes in the near-wall region when the MOC is reduced. The change in the evaporation term is balanced by a change in the diffusive term. The inertial term is determined by the velocity and enthalpy fields, which differ little overall, and is therefore insensitive to such small local temperature variations. This balance for MOC = 0.1 °C is in better agreement with the experimental data for the liquid temperature profile and leads to a restriction of saturated liquid much closer to the wall. Thus,

this modification of the wall HFP model is used to reduce the liquid superheating. The thermal effects near the wall have little impact on the rest of the flow which means that the radial transport effects are relatively low.

### 5.1.4. Conclusion

To summarize the analyses of the reference case, the flow can be divided in three zones with different predominant effects.

In the near-wall region, the gradients of liquid temperature and velocity are high, which induces a high lift and thermal turbulent diffusion effects, along with the high evaporation due to the wall HFP.

In an intermediate zone, the condensation phenomenon is predominant in the thermal balance. In the meantime, the vapor phase spreads out to the pipe axis due to the turbulent dispersion force.

In the center, the flow is still single-phase. The liquid temperature is the result of an equilibrium between the wall heat flux and the axial inertial thermal effects.
<figure><img src="image_17.png" /><figcaption>Fig. 15. Axial HFP with different MOC values for test tube G2P26W16. The vertical line represents the axial coordinate where $T_i = T_{sat}$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Radial profile of the enthalpy balance – G2P26W16Te68.</figcaption></figure>

## 5.2. Mass flow rate effect

In this section, we compare two cases where the mass flow rate and the wall heat flux change, all other quantities being as similar as possible. We select cases G2P26W16Te68 and G3P26W23Te69, which have the same quality $x_{eq,out} = 0.06$ and pressure $P = 26$ bar, but the second case has approximately 50% higher mass flow rate and therefore the bulk Reynolds number increases from $Re = 4.1 \times 10^5$ to $Re = 6.2 \times 10^5$. Their flow conditions are described in Table 3. In the following, MOC is taken equal to 0.1 °C.

Experimentally, a higher mass flow rate results in a flatter void fraction profile (see Fig. 17(a)). The mean void fraction in the pipe cross-section is roughly the same in both cases (as shown in Table 4), but is lower at the wall and higher in the center for the G3P26W23Te69 case. The simulations give the right order of magnitude with respect to the experimental data. The mean void fraction is overestimated for

both cases, but the void fraction at the core is well predicted for the G3P26W23Te69 case. For both cases, the simulations predict a void fraction gradient that is too high.

At the core, the experimental diameter $D_S$ of the two cases are equal (see Fig. 17(b)). The higher degree of subcooling of the G2P26W16Te68 case induces a slightly lower Sauter mean diameter at the center of the pipe than at $r^+ = 0.6$ partly because condensation shrinks the bubbles. For both cases, the liquid is subcooled for $r^+ \le 0.9$.

### 5.2.1. Dynamic effects

Relative velocities are presented on Fig. 18 and momentum balances for the vapor phase on Fig. 19. For both cases, the prediction of the relative axial velocity $\bar{u}_{r,z}$ is almost the same. Indeed, this quantity is classically considered to be the result of an equilibrium between the
<figure><img src="image_19.png" /><figcaption>Fig. 17. Physical quantities – G2P26W16Te68 and G3P26W23Te69.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 18. Relative velocities (plotted only if $\alpha \ge 0.01$) – G2P26W16Te68 and G3P26W23Te69.</figcaption></figure>

drag and buoyancy forces, as described in [48]:

$$ \hat{u}_{r,z} = \sqrt{\frac{4 \Delta \bar{\rho}}{3} \frac{g D_S}{\bar{\rho}_l} C_D (\alpha, D_S)} \quad (20) $$

Depending on the coefficient $C_D$, $\hat{u}_{r,z}$ only depends on $D_S$ and $\alpha$ (which are similar for both cases) and $\Delta\bar{\rho}$. $\Delta\bar{\rho}$ is mainly determined by the pressure, which is the same for both cases. Similarly, $\bar{\rho}_c\delta$ is identical, so the gravity source term in the axial momentum balance (see Fig. 19(a)) is the same for both cases. Even if the relative velocity is the same for both cases, the drag contribution increases significantly for the G3P26W23Te69 case in the axial balance. This is due to the fact that this term is proportional to the interfacial area (see Eq. (3) and Fig. 23(a)). This higher drag is mainly compensated by a higher convection term because the vapor velocities are higher. At the pipe center, the convection and drag contributions for case G3P26W23Te69 are not zero because some vapor remains at this location.

In the radial direction, the drag is greater for case G3P26W23Te69 for the same reason (see Fig. 19(b)). This higher drag is balanced by an increased turbulent dispersion. Indeed, the higher shear rate of the liquid (see Fig. 20(a)) increases the Reynolds stress $R_{zz}$ shown on Fig. 20(b). This in turn increases the turbulent dispersion contribution. The consequence of this modification of the radial momentum balance is a relative radial velocity $\hat{u}_{r,r}$ 60% larger for the G3P26W23Te69 case (see Fig. 18(b)), which flattens the void fraction profile.

### 5.2.2. Thermal effects

The turbulent diffusivity and liquid subcooling are presented on Fig. 21. In addition to the effect of mass flow rate on turbulent dispersion force, turbulent diffusivity increases with mass flow rate (see Fig. 21(a)). The result is an increase of the radial transport of the heat via the diffusion term of the liquid enthalpy balance $\nabla \cdot [a_l \left( \frac{\dot{\lambda}_l}{C_{p,l}} + \frac{\mu_l^T}{P_{r,l}} \right) \nabla \hat{H}_l]$, which flattens the liquid temperature profile shown on Fig. 21(b). Indeed, for the same enthalpy gradient, the diffusion effect is amplified by a higher turbulent diffusivity.

The liquid enthalpy balance terms are shown on Fig. 22. The effect of a higher mass flow rate on condensation is complex due to the coupling of the different factors in Eq. (8). The diameter $D_S$ of the two cases are similar, and the mass flow rate has almost no effect on the relative velocity magnitude (see Fig. 18). Hence the bubble Reynolds number $Re_b$ and Nusselt number $Nu$ are similar too, as shown on Fig. 23(b). As mentioned before, the interfacial area $a_l$ is slightly higher for the G3P26W23Te69 case (Fig. 23(a)). But the degree of local subcooling $\Delta T_l = T_{sat} - T_l$ is 2.5 times lower due to the enthalpy diffusion, so all put together, the condensation term $\Pi'_l$ is lower when the mass flow rate is higher.

Here, condensation is not zero at the center of the flow for the G3P26W23Te69 case due to the spreading of the vapor phase (see Fig. 22(a)). In the enthalpy budget, the value of the diffusion term at $r^+ = 0$ is lower than 1 because part of the wall power goes to the vapor.
<figure><img src="image_21.png" /><figcaption>Fig. 19. Vapor momentum balances – G2P26W16Te68 and G3P26W23Te69.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 20. Dynamic quantities – G2P26W16Te68 and G3P26W23Te69.</figcaption></figure>

Compensation by the inertial term is less important in absolute value, as the factors $\alpha_l$ and $\nabla \hat{H}_l$ are smaller than for the G2P26W16Te68 case.

### 5.2.3. Conclusion

The higher turbulence level of the G3P26W23Te69 case induces higher turbulent dispersion and turbulent diffusivity, which respectively flatten the liquid temperature and the void fraction profiles. In addition, the condensation in the core is higher, which flattens the liquid temperature profile even more.

For both cases, the simulated gradients of the liquid temperature and of the void fraction are too steep with respect to the experimental ones. The cause of this seems to be an underestimation of turbulent quantities $R_{zz}$ and $\mu_l^T$ and/or an overestimation of the turbulent Prandtl number $Pr^T$ which play a key role in the radial transport of liquid heat and vapor mass. Furthermore, the lift model may be partly responsible for the accumulation of the vapor near the wall. With improved radial transport, the condensation would increase in the core region and the

cross-section average of the void fraction would be more in line with the measurements.

### 5.3. Thermodynamic quality effect

In this section, we compare cases where the only parameter that varies is the thermodynamic quality. We select three cases of the G2P26W16 test tube: G2P26W16Te58 is subcooled ($x_{eq} < 0$), the quality of G2P26W16Te64 is higher ($x_{eq} \approx 0$) and G2P26W16Te68 is near saturation ($x_{eq} > 0$). Their flow conditions are described in Table 3.

As expected, the experimental and simulated void fraction increases with the thermodynamic quality, and the area where $\alpha > 0$ expands, as shown on Fig. 24(a). The experimental diameter (see Fig. 24(b)) also increases as condensation tends to be weaker. For these cases, the cross-section average void fraction is largely overestimated, as shown in Table 4. Similarly as in Section 5.2, the void fraction is especially overestimated near the wall, and underestimated in the core.
<figure><img src="image_23.png" /><figcaption>Fig. 21. Thermal diffusion – G2P26W16Te68 and G3P26W23Te69.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 22. Radial profile of the enthalpy balance – G2P26W16Te68 and G3P26W23Te69.</figcaption></figure>

Qualitatively, the experimental void fraction profile of case G2P26W16Te58 is approximately reproduced by a much more sub-cooled case, G2P26W16Te54 which corresponds to $x_{eq} = -0.14$. This shows that the code predicts the Onset of Nucleate Boiling (ONB) at a too low thermodynamic quality.

### 5.3.1. Dynamic effects

Fig. 25 shows the axial velocities for the G2P26W16Te58, G2P26W16Te64 and G2P26W16Te68 cases. The higher void fraction reduces the liquid cross-sectional area, hence increases the axial liquid velocity by conservation of mass flow rate (see Fig. 25(a)). Also, it induces a higher axial vapor velocity thanks to the density difference, as shown on Fig. 25(b). Hence, the Reynolds stress $R_{zz}$ increases with the thermodynamic quality (see Fig. 26(a)) because the liquid velocity gradient is higher.

Axial and radial momentum balances for the vapor phase are presented on Fig. 27 for cases G2P26W16Te64 and G2P26W16Te68. In the radial direction, the turbulent dispersion force is increased by the higher turbulence (see Fig. 26(a) and Fig. 27(b)). It increases the gas

radial velocity toward the core for $r^+ < 0.6$ (see Fig. 28(b)). Near the wall, the increase in liquid shear rate (Fig. 25(a)) induces a more predominant lift force in the radial balance. It decreases the gas radial velocity toward the core for $r^+ > 0.6$. The radial drag force is of the same magnitude due to the similar simulated interfacial area and radial relative velocity, as shown on Figs. 26(b) and 28(b).

### 5.3.2. Thermal effects

Fig. 29 shows the interfacial heat transfer and the liquid subcooling for the G2P26W16Te58, G2P26W16Te64 and G2P26W16Te68 cases. The superheated region slightly expands when the thermodynamic quality increases (see Fig. 29(a)). In addition, we observe a shift in the maximum of the condensation term toward the core, as part of the vapor phase has migrated in this direction (see Fig. 29(b)). In the subcooled cross-section region of the flow, the average condensation is almost twice as high for the G2P26W16Te64 case as for the G2P26W16Te68 case, as shown on the integrated balance Fig. 30(b). For the G2P26W16Te64 case, total evaporation in the superheated zone is lower because this zone is narrower and the interfacial area is smaller in it.
<figure><img src="image_25.png" /><figcaption>Fig. 23. Condensation terms – G2P26W16Te68 and G3P26W23Te69.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 24. Physical quantities – G2P26W16Te58, G2P26W16Te64 and G2P26W16Te68.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 25. Axial velocities – G2P26W16Te58, G2P26W16Te64 and G2P26W16Te68.</figcaption></figure>

In the region $r^+ \le 0.4$ of Fig. 30(a), the flow is single-phase for the case G2P26W16Te64, hence there is no condensation. The

single-phase zone expands to $r^+ \le 0.6$ for case G2P26W16Te58. Here, the inertial and diffusive terms are constant, as described in Section 5.1.
<figure><img src="image_28.png" /><figcaption>Fig. 26. Physical quantities – G2P26W16Te58, G2P26W16Te64 and G2P26W16Te68.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 27. Vapor momentum balances – G2P26W16Te64 and G2P26W16Te68.</figcaption></figure>

For $0.6 < r^+ < 0.95$ of the case G2P26W16Te64, the liquid is strongly heated by condensation while inertial cooling is low due to the low velocity of the liquid. As a result, the diffusion effects are reversed and contribute to the cooling of the flow by diffusing the heat into the subcooled zone. This effect is amplified for case G2P26W16Te58 in the region $0.9 < r^+ < 1$ due to the more pronounced condensation peak.

### 5.3.3. Conclusion

Thermodynamic quality mainly changes the void fraction and consequently the axial velocities of the two phases. An increase in quality increases turbulence and hence enhances radial transport of liquid enthalpy and vapor phase. Evaporation is intense for all three cases close to the wall, especially for the high quality case.

## 6. Conclusions and prospects

This paper presents the modeling and simulations using a two-fluid Euler–Euler RANS model of the DEBORA experiment, which is

<table><caption>Table 4
Experimental and predicted cross-section average void fractions.</caption><thead><tr><td>Case</td><td>x<sub>eq,out</sub></td><td>⟨α<sub>exp</sub>⟩</td><td>⟨α<sub>NCFD</sub>⟩</td></tr></thead><tbody><tr><td>G3P26W23Te69</td><td>0.06</td><td>0.22</td><td>0.28</td></tr><tr><td>G2P26W16Te68</td><td>0.06</td><td>0.22</td><td>0.29</td></tr><tr><td>G2P26W16Te64</td><td>0.01</td><td>0.09</td><td>0.17</td></tr><tr><td>G2P26W16Te58</td><td>-0.07</td><td>0.01</td><td>0.03</td></tr></tbody></table>

representative of the thermal-hydraulic conditions of pressurized water reactors. The approach adopted here is to constrain the calculation by enforcing experimental diameters in the Neptune_CFD simulations. This eliminates some of the coupling mechanisms between thermal and dynamic phenomena, thus reducing the complexity of the calculations and improving confidence in the modeling. Standard models have been used, with the exception of the diameters and the wall Heat Flux Partitioning (HFP), which is constrained to limit liquid superheating in line with experimental observations. Here, we have chosen to reproduce
<figure><img src="image_30.png" /><figcaption>Fig. 28. Relative velocities – G2P26W16Te58, G2P26W16Te64 and G2P26W16Te68.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Fig. 29. Thermal quantities – G2P26W16Te58, G2P26W16Te64 and G2P26W16Te68.</figcaption></figure>

the experimental trends, preferring to focus on physical phenomena rather than adjusting coefficients in closure laws to better match the experiments.

Analysis of the momentum balance for the vapor phase shows that in the axial direction, we find the classic equilibrium between the drag and pressure gradient contributions. With our set of models, in the radial direction, the predominant effects are turbulent dispersion moving the vapor away from the wall, balanced by the drag and lift forces. The role of the pressure gradient is of secondary importance in this direction. An increase in mass flow rate or thermodynamic quality increases liquid shear rate and the turbulent quantities, in turn modifying the equilibrium between these effects. In addition, the enthalpy balance terms have been physically analyzed. In general terms, the liquid temperature is determined by an equilibrium between the radial turbulent diffusion of wall heat, the axial inertial thermal effects and an interphase exchange by condensation or evaporation. The radial diffusion is important at any distance from the wall, whereas the inertial effect is only non-negligible in the core. The evaporation effect is predicted to be significant near the wall, although this phenomenon has not been observed experimentally, which is a source of simulation error in the mass and enthalpy balances. Nevertheless, the impact of this discrepancy on the prediction of void fraction distribution is small. Condensation reaches its maximum in an intermediate radial position between the center and the wall, which depends on the thermodynamic quality. The predominant effects by zone are summarized in Fig. 31.

The simulated quantities are compared to the experimental data. The main point is that the simulated profiles of temperature and void fraction gradients are steeper than the experimental ones. The turbulent quantities driving radial transport of the liquid heat (via turbulent diffusion) and vapor mass (via turbulent dispersion) toward the center of the flow appear to be underestimated. The lift force contributes too strongly to the accumulation of the vapor phase near the wall. Under the conditions considered, liquid overheating is partially treated by modifying the MOC parameter without affecting the rest of the flow.

This paper shows that it is necessary to identify precisely the source discrepancies between calculations and experience to know how to address them. To understand what leads to a prediction of a quantity of interest, it is necessary to present and analyze other quantities, in particular relative velocities, turbulent stress, turbulent viscosity, interfacial area and Nusselt number. Analyzing the momentum and enthalpy balances give valuable additional information on phenomena interaction or compensation. This methodology should be applied to the selection of models for calculation codes.

This study will be continued by investigating the effect of the pressure on the momentum and enthalpy balances. Indeed, the larger and more deformable bubbles that appear in lower-pressure cases such as the 14 bar DEBORA tests tend to migrate toward the axis until the void fraction is core-peaked. A prerequisite for this study is the availability of a lift model capable of predicting this migration. Furthermore, this study suggests that the turbulent diffusivity model using a fixed value
<figure><img src="image_32.png" /><figcaption>Fig. 30. Radial profile of the enthalpy balance – G2P26W16Te58, G2P26W16Te64 and G2P26W16Te68.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Fig. 31. Summary of the main effects by zone.</figcaption></figure>

$P_r^T \approx 1$ is insufficient to reproduce the increase in mixing due to the presence of bubbles. Consequently, the development of a more universal model for turbulent heat transfer would considerably improve the prediction capability of the model.

### CRediT authorship contribution statement

**Etienne Chavasse-Frétaz:** Writing – review & editing, Writing – original draft, Visualization, Software, Methodology, Investigation, Conceptualization. **Alan Burlot:** Writing – review & editing, Supervision, Methodology, Conceptualization, Data curation. **Guillaume Bois:** Writing – review & editing, Supervision, Methodology, Conceptualization. **Catherine Colin:** Writing – review & editing, Supervision, Methodology, Conceptualization.

### Declaration of competing interest

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

### Data availability

Data will be made available on request.

### References

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

[2] H. Wang, H. Luo, J. Xiong, Iterative screening methodology for optimal modeling of bubble coalescence/breakup and interphase force in CMFD simulation of flow boiling, Int. J. Multiph. Flow 170 (2024) 104644, https://doi.org/10.1016/j.ijmultiphaseflow.2023.104644.

[3] N. Shaparia, U. Pelay, D. Bougeard, A. Levasseur, N. François, S. Russeil, Investigation of wall boiling closure, momentum closure and population balance models for refrigerant gas-liquid subcooled boiling flow in a vertical pipe using a two-fluid Eulerian CFD model, Energies 17 (17) (2024) 4225, https://doi.org/10.3390/en17174225.

[4] G. Bois, P. Fillion, F. François, A. Burlot, A.B.H. Ali, A. Khaware, J. Sanyal, M. Rehm, B. Farges, F. Vinauger, W. Ding, A. Gajšek, M. Tekavčič, B. Končar, J.M.L. Corre, H. Li, R. Härlin, J. Jaseliūnaitė, E. Baglietto, R. Brewster, A. Ding, D. Vlček, Y. Sato, J. Xiong, H. Wang, H. Luo, L. Vyskocil, V. Hovi, Benchmark DEBORA: Assessment of MCFD compared to high-pressure boiling pipe flow
measurements, Int. J. Multiph. Flow 179 (2024) 104920, https://doi.org/10.1016/j.ijmultiphaseflow.2024.104920.

[5] E. Baglietto, E. Demarly, R. Kommajosyula, N. Lubchenko, B. Magolan, R. Sugrue, A second generation multiphase-CFD framework toward predictive modeling of DNB, Nucl. Technol. 205 (1-2) (2019) 1-22, https://doi.org/10.1080/00295450.2018.1517528.

[6] M. Pham, G. Bois, F. Francois, E. Baglietto, Assessment of State-of-the-art multiphase CFD modeling for subcooled flow boiling in reactor applications, Nucl. Eng. Des. 411 (2023) 112379, https://doi.org/10.1016/j.nucengdes.2023.112379.

[7] D. Lucas, R. Rzehak, E. Krepper, T. Ziegenhein, Y. Liao, S. Kriebitzsch, P. Apanasevich, A strategy for the qualification of multi-fluid approaches for nuclear reactor safety, CFD4NRS-5, Nucl. Eng. Des. 299 (2016) 2-11, https://doi.org/10.1016/j.nucengdes.2015.07.007.

[8] J. Garnier, E. Manon, G. Cubizolles, Local measurements on flow boiling of refrigerant 12 in a vertical tube, Multiph. Sci. Technol. 13 (1-2) (2001) 111, https://doi.org/10.1615/MultScienTechn.v13.i1-2.10.

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

[10] C. Reiss, A. Gerschenfeld, C. Colin, Boiling-flow multiphase CFD simulations for nuclear reactor conditions without interfacial area transport equation, Nucl. Eng. Des. 428 (2024) 113453, https://doi.org/10.1016/j.nucengdes.2024.113453.

[11] M. Ishii, T. Hibiki, Thermo-Fluid Dynamics of Two-Phase Flow, Springer, New York, NY, 2011, https://doi.org/10.1007/978-1-4419-7985-8.

[12] C. Morel, Mathematical modeling of disperse two-phase flows, in: Fluid Mechanics and Its Applications, vol. 114, Springer International Publishing, Cham, 2015, https://doi.org/10.1007/978-3-319-20104-7.

[13] M. Ishii, N. Zuber, Drag coefficient and relative velocity in bubbly, droplet or particulate flows, AIChE J. 25 (5) (1979) 843-855, https://doi.org/10.1002/aic.690250513.

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

[15] P.-A. Haynes, Contribution à La Modélisation de La Turbulence Pour Les Écoulements à Bulles : Proposition d'un Modèle (k-Epsilon) Multi-Échelles Diphasique (Ph.D. thesis), Toulouse, INPT, 2004.

[16] J. Laviéville, N. Mérigoux, M. Guingo, C. Baudry, S. Mimouni, A Generalized turbulent dispersion model for bubbly flow numerical simulation in NEPTUNE_CFD, 16th International Topical Meeting on Nuclear Reactor Thermal Hydraulics, Nucl. Eng. Des. 312 (2017) 284-293, https://doi.org/10.1016/j.nucengdes.2016.11.003.

[17] R. Manceau, K. Hanjalić, Elliptic blending model: A new near-wall Reynolds-stress turbulence closure, Phys. Fluids 14 (2) (2002) 744-754, https://doi.org/10.1063/1.1432693.

[18] EDF-R&D, NEPTUNE CFD Version 7.2.0 Theory Guide, Technical Report, EDF R&D, Chatou, 2022, p. 65.

[19] S. Mimouni, W. Benguigui, J. Laviéville, N. Merigoux, M. Guingo, C. Baudry, O. Marfaing, New nucleation boiling model devoted to high pressure flows, in: 9th International Conference on Multiphase Flow, ICMF, Firenze, Italy, 2016, pp. 1-6.

[20] S. Mimouni, C. Baudry, M. Guingo, J. Lavieville, N. Merigoux, N. Mechitoua, Computational multi-fluid dynamics predictions of critical heat flux in boiling flow, CFD4NRS-5, Nucl. Eng. Des. 299 (2016) 28-36, https://doi.org/10.1016/j.nucengdes.2015.07.017.

[21] N. Kurul, M.Z. Podowski, Multidimensional effects in forced convection subcooled boiling, in: Proceeding of International Heat Transfer Conference 9, Begellhouse, Jerusalem, Israel, 1990, pp. 21-26, https://doi.org/10.1615/IHTC9.40.

[22] H.C. Ünal, Maximum bubble diameter, maximum bubble-growth time and bubble-growth rate during the subcooled nucleate flow boiling of water up to 17.7 MN/M2, Int. J. Heat Mass Transfer 19 (6) (1976) 643-649, https://doi.org/10.1016/0017-9310(76)90047-8.

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

[24] C.S. Brooks, T. Hibiki, Wall nucleation modeling in subcooled boiling flow, Int. J. Heat Mass Transfer 86 (2015) 183-196, https://doi.org/10.1016/j.ijheatmasstransfer.2015.03.005.

[25] Y.-C. Lin, Y. Zhao, M. Ishii, J.P. Schlegel, K.J. Hogan, J.R. Buchanan, Assessment of nucleation boiling models and improvement by the chen correlation for Two-fluid Model CFD, Int. J. Heat Mass Transfer 175 (2021) 121363, https://doi.org/10.1016/j.ijheatmasstransfer.2021.121363.

[26] L. Favre, C. Colin, S. Pujet, S. Mimouni, An updated force balance approach to investigate bubble sliding in vertical flow boiling at low and high pressures, Int. J. Heat Mass Transfer 211 (2023) 124227, https://doi.org/10.1016/j.ijheatmasstransfer.2023.124227.

[27] Y. Zhao, Y.-C. Lin, K. Tang, M. Ishii, J.R. Buchanan-Jr, Investigation of nucleate boiling and bubble dynamics at multiple nucleation sites in subcooled boiling flow, Int. J. Heat Mass Transfer 219 (2024) 124886, https://doi.org/10.1016/j.ijheatmasstransfer.2023.124886.

[28] Y. Zhao, Z. Dang, K. Tang, M. Ishii, J.R. Buchanan, Characterization of departure diameter, departure frequency and active nucleation site density in subcooled boiling flow, Appl. Therm. Eng. (2025) 125497, https://doi.org/10.1016/j.applthermaleng.2025.125497.

[29] EDF-R&D, Neptune_cfd 7.0 Verification and Validation : Four-Party NEPTUNE Report, Technical Report, EDF R&D, Chatou, 2021, p. 65.

[30] L. Favre, Modélisation et Simulation de La Crise d'ébullition Dans Les REP à l'échelle CFD (Ph.D. thesis), Toulouse, INPT, 2023.

[31] S. Mimouni, F. Archambeau, M. Boucker, J. Laviéville, C. Morel, A second-order turbulence model based on a Reynolds stress approach for two-phase flow—Part I: adiabatic cases, Sci. Technol. Nucl. Install. 2009 (2008) e792395, https://doi.org/10.1155/2009/792395.

[32] A. Guelfi, D. Bestion, M. Boucker, P. Boudier, P. Fillion, M. Grandotto, J.-M. Hérard, E. Hervieu, P. Pé PURAUD, NEPTUNE: A new software platform for advanced nuclear thermal hydraulics, Nucl. Sci. Eng. 156 (3) (2007) 281-324, https://doi.org/10.13182/NSE05-98.

[33] S. Mimouni, M. Guingo, J. Lavieville, N. Mérigoux, Combined evaluation of bubble dynamics, polydispersion model and turbulence modeling for adiabatic two-phase flow, Multi-Scale Multi-Physics Analysis of Reactor Transients in the NURESAFE Project, Nucl. Eng. Des. 321 (2017) 57-68, https://doi.org/10.1016/j.nucengdes.2017.03.041.

[34] T. Herry, B. Raverdy, S. Mimouni, S. Vincent, Modeling and simulation of bubble condensation using polydispersion approach with bubble collapse model, Int. J. Multiph. Flow 176 (2024) 104846, https://doi.org/10.1016/j.ijmultiphaseflow.2024.104846.

[35] A. Tomiyama, Struggle with computational bubble dynamics, Multiph. Sci. Technol. 10 (4) (1998) 369-405, https://doi.org/10.1615/MultScienTechn.v10.i4.40.

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

[37] E. Manon, Contribution à l'analyse et à La Modélisation Locale Des Écoulements Bouillants Sous-Saturés Dans Les Conditions Des Réacteurs à Eau Sous Pression (Ph.D. thesis), Ecole centrale de Paris, Châtenay-Malabry, 2000.

[38] P. Berne, Contribution à La Modélisation Du Taux De Production De Vapeur Par Autovaporisation Dans Les Écoulements Diphasiques En Conduite (Ph.D. thesis), Ecole Centrale des Arts et Manufactures, Châtenay-Malabry, 1983.

[39] J.-M. Delhaye, Thermohydraulique Des Réacteurs: Edition Révisée, EDP Sciences, 2021, https://doi.org/10.1051/978-2-7598-1147-2.

[40] N.E. Todreas, M.S. Kazimi, Nuclear Systems: Thermal Hydraulic Fundamentals, CRC Press, Taylor & Francis Group, 2021.

[41] W. Yao, C. Morel, Prediction of parameters distribution of upward boiling two-phase flow with two-fluid models, in: 10th International Conference on Nuclear Engineering, 3, Arlington, Virginia, USA, 2002, pp. 1-8, https://doi.org/10.1115/ICONE10-22463.

[42] D. Shaver, M. Podowski, Modeling of interfacial forces for bubbly flows in subcooled boiling conditions, in: ANS Winter Meeting 2015, Washington, DC, 2015, p. 4.

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

[44] H. Lyu, D. Lucas, R. Rzehak, F. Schlegel, A particle-center-averaged Euler-Euler model for monodisperse bubbly flows, Chem. Eng. Sci. 260 (2022) 117943, https://doi.org/10.1016/j.ces.2022.117943.

[45] H. Lyu, F. Schlegel, R. Rzehak, D. Lucas, Euler-Euler model of bubbly flow using particle-center-averaging method, Nucl. Sci. Eng. 197 (10) (2023) 2602-2619, https://doi.org/10.1080/00295639.2022.2131344.

[46] P. Rival, F. Francois, V. Faucher, An improved model for the axial vapor velocity in the DEBORA databank, Int. J. Multiph. Flow (2025) (in revision).

[47] M. Tanaka, Numerical study on flow structures and heat transfer characteristics of turbulent bubbly upflow in a vertical channel, in: J. Zhu (Ed.), Computational Simulations and Applications, InTech, 2011, pp. 119-142, https://doi.org/10.5772/24405.

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