# Benchmark DEBORA: Assessment of MCFD compared to high-pressure boiling pipe flow measurements

G. Bois, P. Fillion, F. François, A. Burlot, A. Ben Hadj Ali, A. Khaware, J.

Sanyal, M. Rehm, B. Farges, F. Vinauger, et al.

► To cite this version:

G. Bois, P. Fillion, F. François, A. Burlot, A. Ben Hadj Ali, et al., Benchmark DEBORA: Assessment of MCFD compared to high-pressure boiling pipe flow measurements. International Journal of Multiphase Flow, 2024, 179, pp.104920. 10.1016/j.ijmultiphaseflow.2024.104920. cea-04695258

HAL Id: cea-04695258

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

Submitted on 13 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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>

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

Check for
updates

# Benchmark DEBORA: Assessment of MCFD compared to high-pressure boiling pipe flow measurements

G. Bois<sup>a,*</sup>, P. Fillion<sup>a</sup>, F. François<sup>b</sup>, A. Burlot<sup>a</sup>, A. Ben Hadj Ali<sup>c</sup>, A. Khaware<sup>d</sup>, J. Sanyal<sup>e</sup>, M. Rehm<sup>f</sup>, B. Farges<sup>g</sup>, F. Vinauger<sup>g</sup>, W. Ding<sup>h</sup>, A. Gajšek<sup>i,j</sup>, M. Tekavčič<sup>i</sup>, B. Končar<sup>i</sup>, J.-M. Le Corre<sup>k</sup>, H. Li<sup>l</sup>, R. Härlin<sup>l,k</sup>, J. Jaseliūnaitė<sup>m</sup>, E. Baglietto<sup>n</sup>, R. Brewster<sup>n</sup>, A. Ding<sup>n</sup>, D. Vlček<sup>o</sup>, Y. Sato<sup>p</sup>, J. Xiong<sup>q</sup>, H. Wang<sup>q</sup>, H. Luo<sup>q</sup>, L. Vyskocil<sup>r</sup>, V. Hovi<sup>s</sup>

## ARTICLE INFO

Dataset link: https://doi.org/10.5281/zenodo.12551306

### Keywords:

Benchmark

MCFD

DEBORA experiment

High-pressure boiling flow

## ABSTRACT

A benchmark activity on two-fluid simulations of high-pressure boiling upward flows in a pipe is performed by 12 participants using different MCFD (Multiphase Computational Fluid Dynamics) codes and closure relationships. More than 30 conditions from DEBORA experiment conducted by CEA are considered. Each case is characterised by the flow rate, inlet temperature, wall heat flux and outlet pressure. High-pressure Freon (R12) at 14 bar and 26 bar is boiled in a 19.2 mm pipe heated over 3.5 m. Flow rates range from 2000 kg m⁻² s⁻¹ to 5000 kg m⁻² s⁻¹ and exit quality x ranges from single-phase conditions to x = 0.1 which leads to a peak void fraction of α = 70%. In these high pressure conditions, bubbles remain small and there is no departure from the bubbly flow regime (François et al., 2011; Hösler, 1968). However, different kind of bubbly flows are observed: wall-peak, intermediate peak or core-peak, depending on the case considered. Measurements along the pipe radius near the end of the heated section are compared to code predictions. They include void fraction, bubble mean diameter, vapour velocity and liquid temperature. The benchmark covered two phases. In the first phase of the benchmark activities, experimental data were given to the participants, allowing to compare the simulation results and to develop, to select or to adjust the models in the CMFD codes. The second phase included blind cases where the participants could not compare to the measurements. In between the two phases, possible additional model adjustments or calibrations were performed.

Overall, the benchmark involved very different closures and a wide range of models' complexity was covered. Yet, it is extremely difficult to have a robust closure for all conditions considered, even knowing experimental measurements. The wall-to-core peak transition is not captured consistently by the models. The degree of subcooling and the void fraction level are also difficult to assess. We were not capable of showing

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

<footnote><sup>b</sup> CEA, DES, IRESNE, Department of Nuclear Technology, CEA Centre de Cadarache, Saint-Paul-Lez-Durance, F-13108, France</footnote>

<footnote><sup>c</sup> Ansys Germany GmbH (ANSYS), Birkenweg 14a, Darmstadt, 64295, Germany</footnote>

<footnote><sup>d</sup> Ansys Software Pvt Ltd, Research and Development (ANSYS), Plot-34/1 Hinjewadi, Pune, 411057, India</footnote>

<footnote><sup>e</sup> Ansys Inc., Research and Development (ANSYS), 10 Cavendish Ct, Lebanon, 03766, USA</footnote>

<footnote><sup>f</sup> Framatome GmbH, Fuel Core Design, Paul-Gossen-Str. 100, Erlangen, 91052, Germany</footnote>

<footnote><sup>g</sup> Framatome SAS, Fuel Core Physics, 2 rue Professeur Jean Bernard, Lyon, 69007, France</footnote>

<footnote><sup>h</sup> Helmholtz Zentrum Dresden Rossendorf, Institute of Fluid Dynamics, Bautzner landstr. 400, Dresden, 01328, Germany</footnote>

<footnote><sup>i</sup> Jožef Stefan Institute, Reactor Engineering Division (R4), Jamova cesta 39, Ljubljana, 1000, Slovenia</footnote>

<footnote><sup>j</sup> University of Ljubljana, Kongresni trg 12, Ljubljana, 1000, Slovenia</footnote>

<footnote><sup>k</sup> Westinghouse Electric Sweden AB, Västerå, 72163, Sweden</footnote>

<footnote><sup>l</sup> Nuclear Engineering, Department of Physics, Royal Institute of Technology, 106 91, Stockholm, Sweden</footnote>

<footnote><sup>m</sup> Lithuanian Energy Institute, Laboratory of Nuclear Installation Safety, Breslaujos st. 3, Kaunas, 44403, Lithuania</footnote>

<footnote><sup>n</sup> Massachusetts Institute of Technology, Department of Nuclear Science and Engineering, 77 Massachusetts Avenue, Cambridge, MA, 02139, USA</footnote>

<footnote><sup>o</sup> Czech Technical University in Prague, Faculty of Nuclear Sciences and Physical Engineering, Department of Nuclear Reactors, V Holesovickach 2, Prague, 18000, Czech Republic</footnote>

<footnote><sup>p</sup> Paul Scherrer Institute, Division Scientific Computing, Theory and Data, Forschungsstrasse 111, Villigen PSI, 5232, Switzerland</footnote>

<footnote><sup>q</sup> Shanghai Jiao Tong University (SJTU), 800 Dongchuan Road Minhang, Shanghai, 200240, China</footnote>

<footnote><sup>r</sup> UJV Rez, a.s., Department of Safety Analyses, Hlavni 130, Husinec-Rez, 25068, Czech Republic</footnote>

<footnote><sup>s</sup> VTT Technical Research Centre of Finland Ltd., P.O. Box 1000, FI-02044 VTT, Finland</footnote>

<footnote>* Corresponding author.
E-mail address: guillaume.bois@cea.fr (G. Bois).</footnote>
superiority of some physical closures, even for part of the model. The interaction between mechanisms and their hierarchy are extremely difficult to understand.

Although departure from nucleate boiling (DNB) was not considered in this benchmarking exercise, it is expected that DNB predictions at high-pressure conditions depend strongly on the near-wall flow, temperature, and void fraction distributions. Therefore, the suitability of the closures also limits the accuracy of DNB predictions. The benchmark also demonstrated that in order to progress further in models development and validation, it is compulsory to have new measurements that include simultaneously as many variables as possible (including liquid temperature, velocity, cross-correlations and wall temperature); also, a better knowledge of the local bubble sizes distributions is the key to discriminate performances of interfacial area modelling (IATE, MUSIG or iMUSIG models, considering for instance the possibility of two classes of bubbles having totally different behaviour regarding the lift force).

Following this benchmark impulse, we hope that future activities will be engaged on high-pressure boiling water experiments with a continuation of models' comparisons and development.

# 1. Introduction: Benchmark's motivation and scope

Boiling flows can be of major interest to many industrial applications related to nuclear reactor safety, but also water-tubes boilers, heat exchangers for refrigerant systems (evaporators), boilers for chemical processes, etc. Dhir (1998). This kind of flow puts a strain on many aspects of the modelling, including dynamical (such as turbulence modelling and interfacial transfers), topological (laws and closures for the evolution of the bubble sizes distribution) and thermal aspects (interfacial transfers due to sub-cooling and wall heat flux partitioning). In the context of multiphase flow prediction, the Multiphase Computational Fluid Dynamics (MCFD) is a very complex field of research. For industrial applications, many computational tools have been developed and they undergo thorough validation processes. In particular, the two-fluid model relying on phase-averaged equations has attracted a lot of interest for more than a decade because of its wide range of applicability and because it is one of the few tools applicable in challenging industrial conditions.

However, the two-fluid model requires many closure relations and a great diversity of sub-models is being developed. So far, it seems that no set of closures has reached a universal description applicable to any flow condition. Each set of closures has its own specific validation database and their ranges of validity are difficult to compare. The empiricism of the models involved usually limits the applicability of MCFD to predominantly low pressure operating conditions (Muralidharan et al., 2016). This benchmark aims to lead the way towards a more unified method for testing and validating MCFD closures based on publicly available experimental data and not yet published data. We propose to start this massive joint effort by focusing on high-pressure upward flow boiling in a vertical cylindrical tube. For that, test conditions performed on a boiling flow of pressurised Freon experiment (DEBORA) are used to provide information on flow characteristics and form a relevant basis for codes and models comparisons. With DEBORA experiment, a particular focus is set upon the effect of high-pressure conditions on those mechanisms. Even though the geometry is very simple, the test is extremely challenging; it is a good configuration to assess code capabilities from a global perspective, in particular with respect to the void fraction level and spatial distribution. For instance, challenging conditions with high quality and local void fraction are available.

Parts of DEBORA database were used to validate two-phase CFD codes during the European project NURESIM, devoted to "improve the understanding and the predictive capabilities of the simulation tools for key two-phase flow thermal-hydraulic processes that can occur in nuclear reactors, focusing on two high priority issues, the Critical Heat Flux (CHF), and the pressurised thermal shock" (Bestion et al., 2009). This experimental data has also been widely used by the nuclear industry to develop and validate MCFD models (Mimouni et al., 2008; Vyskocil and Macek, 2008; Yun et al., 2010; Mimouni et al., 2010b; Krepper and Rzehak, 2011; Peltola and Pättikangas, 2012; Yun et al., 2012; Krepper et al., 2013; Krepper and Rzehak, 2014; Shaver and Podowski,

2015a; Tripathi and Buwa, 2014; Raj and Velusamy, 2017; Magolan et al., 2019; Peltola et al., 2019; Setoodeh et al., 2021a; Lubchenko, 2000; Kommajosyula, 2020; Mérigoux, 2022; Vlček and Sato, 2023; Gajšek et al., 2023; Wang et al., 2024), but also for application and assessment of sub-channels, drift-flux or system codes (Bestion et al., 2009; François et al., 2011). Thus, available DEBORA data were widely used in the literature to assess MCFD models. In this initiative, we collect a large number of those conditions and provide additional unpublished conditions to build an extensive basis for codes and models validations.

Many MCFD codes have reached a first level of modelling that satisfy some requirements. The benchmark activity provided a great opportunity to challenge these models in particular to attempt to reach finer predictions of the void-fraction level and cross-flow distribution. A large number of participants, codes and models were involved (see Table 1). The following comparison provides useful information on the maturity and capabilities of state-of-the-art MCFD models. One objective is to trigger a collective review of the pros and cons of the models and offer a mutual effort towards the improvement of their weaknesses. The entire community will therefore benefit from a large and critical review of the state-of-the-art models.

This paper starts with a description of the DEBORA experiment (Section 2). A clear compilation of the experimental data and conditions considered (including conditions not yet available in the literature) is provided in Section 3 along with numerical information on the computational domain, meshes and boundary conditions. A brief description of various models used by the participants to the benchmark is given in Section 4. Then, results are summarised, commented and compared in Section 5, allowing to conclude on the benchmark and to propose different ways to improve the simulations. Further specific analyses are presented in Section 6 before drawing the main conclusions of this benchmark (Section 7).

# 2. Presentation of DEBORA experiment

## 2.1. Scaling criteria

High pressure and high temperature steam/water experiments are expensive and difficult to perform. As an alternative, refrigerant fluids such as Freons or Novec are often used to simulate the operating conditions of a pressurised water reactor (PWR) (Yoo et al., 2016; Roy et al., 2001; Estrada-Perez and Hassan, 2010). Thanks to the physical properties of Freon,<fnref n="1" /> it is possible to simulate the prototypical conditions of a PWR but for much lower pressure, temperature and heat fluxes. This scaling is based on the following four non-dimensional parameters :

<footnote n="1">¹ In this paper, we will only consider Freon R12.</footnote>
<caption>**Table 1**</caption>

List of DEBORA benchmark participants and of the corresponding codes. The full names of the organisations and their affiliates are given on the first page.

<table><caption>**Table 2**</caption><thead><tr><th>Organisation</th><th>Acronym</th><th>Country</th><th>MCFD Code</th></tr></thead><tbody><tr><td>Ansys</td><td>ANSYS</td><td>Germany</td><td>ANSYS Fluent 22R2</td></tr><tr><td>Alternative Energies and Atomic Energy Commission</td><td>CEA</td><td>France</td><td>Neptune_CFD v7.2.0</td></tr><tr><td>Framatome</td><td>FRA</td><td>France/Germany</td><td>StarCCM+ v13.02</td></tr><tr><td>Helmholtz-Zentrum</td><td>HZDR</td><td>Germany</td><td>ANSYS CFX 18.0</td></tr><tr><td>Dresden-Rossendorf</td><td></td><td></td><td></td></tr><tr><td>Jožef Stefan Institute</td><td>JSI</td><td>Slovenia</td><td>ANSYS Fluent 21R2</td></tr><tr><td>Royal Institute of Technology/Westinghouse Electric Sweden</td><td>KTH-WSE</td><td>Sweden</td><td>OpenFOAM v9</td></tr><tr><td>Massachusetts Institute of Technology</td><td>MIT</td><td>USA</td><td>STAR-CCM+ 2021.2.1</td></tr><tr><td>Lithuanian Energy Institute</td><td>LEI</td><td>Lithuania</td><td>OpenFOAM v8</td></tr><tr><td>Paul Scherrer Institute/Czech Technical University in Prague</td><td>PSI-CTU</td><td>Switzerland/Czech Republic</td><td>ANSYS Fluent 21R2</td></tr><tr><td>Shanghai Jiao Tong University</td><td>SJTU</td><td>China</td><td>OpenFOAM v6</td></tr><tr><td>UJV Rez</td><td>UJV</td><td>Czech Republic</td><td>ANSYS CFX 18.0</td></tr><tr><td>VTT Technical Research Centre of Finland</td><td>VTT</td><td>Finland</td><td>OpenFOAM dev (April 2022)</td></tr></tbody></table>

Water and R12 corresponding thermal-hydraulic parameters range calculated using scaling criteria Eqs. (1) to (4).

<table><thead><tr><th>Control parameters</th><th>Water</th><th>R12</th></tr></thead><tbody><tr><td>Exit pressure (MPa)</td><td>10 to 18</td><td>1.4 to 3.0</td></tr><tr><td>Mass flux (kg m<sup>-2</sup> s<sup>-1</sup>)</td><td>500 to 5000</td><td>500 to 5000</td></tr><tr><td>Heat flux (MW m<sup>-2</sup>)</td><td>0.5 to 3</td><td>0.01 to 0.20</td></tr><tr><td>Inlet temperature (°C)</td><td>50 to 320</td><td>20 to 80</td></tr><tr><td>Equilibrium exit quality (-)</td><td>-0.3 to 0.5</td><td>-0.3 to 0.5</td></tr></tbody></table>

* Same vapour/liquid density ratio to determine the Freon R12 system pressure $P_s$,

$$ \left( \frac{\rho_l}{\rho_v} \right)_{\text{water}} = \left( \frac{\rho_l}{\rho_v} \right)_{\text{R12}} \qquad (1) $$

where $\rho_l$ and $\rho_v$ are the saturated liquid and vapour densities, respectively.

* Same Weber number (We) to determine the R12 mass flux $G$,

$$ \left( \frac{G^2 R}{\rho_l \sigma} \right)_{\text{water}} = \left( \frac{G^2 R}{\rho_l \sigma} \right)_{\text{R12}} \qquad (2) $$

where $\sigma$ and $R$ are the surface tension and the pipe radius, respectively.

* Same boiling number (Bo) to calculate the R12 corresponding heat flux $q_w$,

$$ \left( \frac{q_w}{G \mathcal{L}^{\text{vap}}} \right)_{\text{water}} = \left( \frac{q_w}{G \mathcal{L}^{\text{vap}}} \right)_{\text{R12}} \qquad (3) $$

where $\mathcal{L}^{\text{vap}}$ is the latent heat of vaporisation,

* Same inlet equilibrium quality $x_{eq,in}$ to determine the R12 inlet temperature

$$ \left( \frac{h_{l,in} - h_l^{\text{sat}}}{\mathcal{L}^{\text{vap}}} \right)_{\text{water}} = \left( \frac{h_{l,in} - h_l^{\text{sat}}}{\mathcal{L}^{\text{vap}}} \right)_{\text{R12}} \qquad (4) $$

where $h_{l,in}$ and $h_l^{\text{sat}}$ are the liquid enthalpy at the inlet of the test section and the saturated liquid enthalpy, respectively.

These criteria define the R12/water corresponding thermal-hydraulic parameters range that is listed in Table 2.

## 2.2. Test facility

DEBORA is a test facility designed and operated at CEA/Grenoble which aimed to provide a reliable database on boiling phenomena in a circular tube. The experimental rig which uses Freon (R12 when being designed then R134 A later) is made of two independent circuits: (i) the test section circuit and (ii) the coolant storage circuit which is not

represented on Fig. 1 but for which a complete description can be found in Garnier et al. (2001).

The test section circuit is composed of:

* A pump (2) with a flow rate of 4.5 kg s<sup>-1</sup> for a 100 m head, the flow rate in the test section being controlled by regulating valves,

* A 30 kW pre-heater (7) which imposes the inlet temperature of the test section,

* A test section (8) which is electrically heated by two power recitifiers each of them delivering 100 kW with a maximum intensity of 1250 A and a maximum voltage of 80 V,

* A condensing unit (9) which mixes the fluid exiting the test section with the incoming liquid from the parallel heat exchanger (4) to get liquid Freon at the inlet of the second heat exchanger (5),

* Two heat exchangers (4) (70 kW, 1.5 kg s<sup>-1</sup>) and (5) (200 kW, 6 kg s<sup>-1</sup>) to cool the Freon to a minimum temperature at the inlet of the test section equal to 20 °C,

* A thermal pressuriser (1) (0.14 m<sup>3</sup>, 6 kW) to impose the pressure at the outlet of the test section (maximum pressure is 4 MPa).

This experimental rig allows to control the four following parameters:

* outlet pressure $P_s$,

* inlet temperature $T_{in}$,

* mass velocity $G$,

* and the heat flux $q_w$.

A test condition is then fully characterised by those four parameters.

## 2.3. Control parameters uncertainties

Outlet pressure is measured using two pressure sensors based on capacity variations located on the top of the test section. The fluid temperature is measured at the inlet and the outlet of the test section using four Platinum sensors, two at the inlet and two at the outlet. Whereas the outlet sensor is located very close to the end of the heated length, the inlet sensor is located 0.5 m upstream of the entrance of the heated tube. According to the thermal insulation of the pipes (rockwool whose thickness is 6 cm) and the low temperature (lower than 100 °C), we assume that this temperature can be considered as the true inlet temperature. The mass flow rate is measured using a Coriolis flowmeter which ranges from 0.17 kg s<sup>-1</sup> to 1.7 kg s<sup>-1</sup>. This flowmeter has been calibrated for several temperatures ranging from 20 °C to 80 °C to cover the whole range of parameters encountered during operation. The electrical power injected in the test section is determined by measuring both the intensity and the voltage across the test section. The voltage
**Nomenclature**

**Acronyms**

<table><tr><td>BIT</td><td>Bubble Induced Turbulence</td></tr><tr><td>CHF</td><td>Critical Heat Flux</td></tr><tr><td>EOS</td><td>Equation Of State</td></tr><tr><td>FBM</td><td>Force Balance Model</td></tr><tr><td>IAC</td><td>Interfacial Area Concentration</td></tr><tr><td>IATE</td><td>Interfacial Area Transport Equation</td></tr><tr><td>iMUSIG</td><td>Inhomogeneous MUSIG</td></tr><tr><td>MCFD</td><td>Multiphase Computational Fluid Dynamics</td></tr><tr><td>MUSIG</td><td>MUltiple SIze Group</td></tr><tr><td>PBM</td><td>Population Balance Model</td></tr><tr><td>PWR</td><td>Pressurised Water Reactor</td></tr><tr><td>RANS</td><td>Reynolds-Averaged Navier Stokes</td></tr><tr><td>RPI</td><td>Rensselaer Polytechnic Institute</td></tr><tr><td>SGDH</td><td>Simple Gradient Diffusion Hypothesis</td></tr><tr><td>WHP</td><td>Wall Heat Partitioning</td></tr></table>

**Subscript**

<table><tr><td>k</td><td>k phase</td></tr><tr><td>l</td><td>Liquid phase</td></tr><tr><td>v</td><td>Vapour phase</td></tr></table>

**Greek symbols**

<table><tr><td>α</td><td>Void fraction</td></tr><tr><td>χ</td><td>Phase indicator function</td></tr><tr><td>η</td><td>Coalescence efficiency</td></tr><tr><td>η<sub>d</sub></td><td>Bubble departure frequency</td></tr><tr><td>Γ<sub>v</sub></td><td>Interfacial mass (evaporation) flux</td></tr><tr><td>λ</td><td>Conductivity</td></tr><tr><td>μ</td><td>Viscosity</td></tr><tr><td>ρ</td><td>Density</td></tr><tr><td>σ</td><td>Surface tension</td></tr><tr><td>τ<sub>v</sub></td><td>Time scale for vapour thermal equilibrium</td></tr></table>

**Latin**

<table><tr><td>c<sub>p</sub></td><td>Thermal capacity at constant pressure</td></tr><tr><td>L<sup>vap</sup></td><td>Latent heat of vaporisation</td></tr><tr><td>A</td><td>Bubble area of influence</td></tr><tr><td>T<sup>t</sup></td><td>Turbulent stress</td></tr><tr><td>T<sup>v</sup></td><td>Viscous stress</td></tr><tr><td>N<sup>e</sup><sub>b</sub></td><td>Nucleation site density</td></tr><tr><td>u</td><td>Velocity vector</td></tr><tr><td>I</td><td>Identity tensor</td></tr><tr><td>M<sub>l→k</sub></td><td>Interfacial transfer into phase k</td></tr><tr><td>q</td><td>Heat flux</td></tr><tr><td>C<sub>Co</sub></td><td>Coalescence factor</td></tr><tr><td>d<sub>d</sub></td><td>Bubble departure diameter</td></tr><tr><td>d<sub>32</sub></td><td>Sauter mean diameter</td></tr><tr><td>F<sub>b</sub></td><td>Break-up coefficient</td></tr><tr><td>F<sub>CB</sub></td><td>Buoyancy coalescence coefficient</td></tr><tr><td>F<sub>CT</sub></td><td>Turbulence coalescence coefficient</td></tr><tr><td>F<sub>LS</sub></td><td>Laminar shear coalescence coefficient</td></tr></table>

is directly measured by means of direct connection to the test section whereas the current is measured using a magnetic intensity converter located close to the power rectifiers. All those sensors are carefully

<table><tr><td>g</td><td>Gravitational acceleration</td></tr><tr><td>h</td><td>Enthalpy</td></tr><tr><td>k<sub>k</sub></td><td>Turbulent kinetic energy of phase k</td></tr><tr><td>Nu</td><td>Nusselt number</td></tr><tr><td>p</td><td>Pressure</td></tr><tr><td>Pr</td><td>Prandtl number</td></tr><tr><td>q<sub>e</sub></td><td>Evaporation heat flux</td></tr><tr><td>q<sub>q</sub></td><td>Quenching heat flux</td></tr><tr><td>T</td><td>Temperature</td></tr><tr><td>t</td><td>Time</td></tr></table>

**Mathematical symbols**

<table><tr><td>(∇φ)<sup>t</sup></td><td>Transposed gradient of φ</td></tr><tr><td>φ&#773;<sup>k</sup></td><td>Phase average of φ</td></tr><tr><td>φ&#773;</td><td>Statistical average of φ</td></tr><tr><td>φ'</td><td>Fluctuation of φ with respect to the phase average</td></tr><tr><td>D<sub>k</sub>/Dt</td><td>Material or substantial derivative</td></tr></table>

**Superscript**

<table><tr><td>sat</td><td>Saturation property</td></tr><tr><td>t</td><td>Turbulent</td></tr></table>

<figure><img src="image_4.png" /><figcaption>Fig. 1. General view of DEBORA loop. The flow direction in the test section is vertical upward.</figcaption></figure>

calibrated. Thanks to the insulation and the moderate temperature in operating conditions, heat losses along the heated pipe were assessed to remain below 60 W for all conditions. Thus, they were neglected.

In order to guarantee the quality of the measurements, the stability of the boundary conditions is checked. The time evolution of the aforementioned control parameters has been registered during seven hours for each given set of thermal-hydraulic conditions. The fluctuations of these parameters remain lower than ±0.5% for the pressure, the mass flow rate and the electrical power, whereas the inlet temperature fluctuations never exceed 0.1 °C. Uncertainties and stability criteria are summed up in Table 3:

## 2.4. Test section

The test section (Fig. 2) is a vertical pipe with a circular cross section. It is made of stainless steel (316Ti). The internal and external
<table><caption>**Table 3**
Measurement uncertainties.</caption><thead><tr><th></th><th>Sensor type</th><th>Accuracy</th><th>Stability</th></tr></thead><tbody><tr><td>Inlet temperature (°C)</td><td>Platinum sensor</td><td>±0.2 °C</td><td>±0.1 °C</td></tr><tr><td>Outlet temperature (°C)</td><td>Platinum sensor</td><td>±0.2 °C</td><td>±0.1 °C</td></tr><tr><td>Mass Flow rate (-)</td><td>Coriolis flow meter</td><td>±0.1%</td><td>±0.5%</td></tr><tr><td>Inlet pressure (bar)</td><td>Capacitive sensor</td><td>±0.01 bar</td><td>±0.5%</td></tr><tr><td>Outlet pressure (bar)</td><td>Capacitive sensor</td><td>±0.01 bar</td><td>±0.5%</td></tr><tr><td>Electrical power (-)</td><td>Voltmeter/ammeter</td><td>±3%</td><td>±0.5%</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** General view of the test section (Garnier et al., 2001).</figcaption></figure>

diameters are equal to 19.2 mm and 21.2 mm, respectively. The heated length is 3.5 m. Thanks to ultrasonic sensors, the thickness of the tube was measured along the vertical axis every 90° with a 100 mm axial step. Electrical resistance was also measured at the same locations to check the axial and azimuthal uniformity of the heat flux. The maximum variations of the local electrical resistance is 0.35% of the mean value. The wall temperature was measured using thermocouples (type K) welded on the external side of the heated tube. The wetted temperature is deduced from a conduction calculation through the heated wall. Four axial locations were monitored, located at 15 mm, 515 mm, 1015 mm and 2015 mm below the end of the heated length. The accuracy is estimated to ±0.2 °C. This data is not used in this benchmark.

Local measurements are performed near the end of the test section at $z = 3.485$ m using a specific device able to rotate in a central plane to cover a diameter, the flow being assumed to be axisymmetric. Thanks to an adapted specific calibration procedure, the accuracy of the probe positioning as well as the moving resolution were determined. They are equal to 50 µm and 10 µm, respectively.

## 2.5. Local measurements

Local measurements are performed by means of specific sensors moved along the diameter of the tube. Three sensors have been used, each of them leading to a specific databank. Space resolution of the measurements is close to 400 µm, which means that a full radial profile corresponds approximately to fifty radial locations.

The flow axisymmetry has been verified by checking the symmetry of the radial profiles measured along a whole diameter.

### 2.5.1. Liquid temperature measurement

The local liquid temperature has been measured by a thermocouple installed on the aforementioned moving device. This chromel/alumel thermocouple has an external diameter of 0.5 mm and an insulated junction. It has been shown by Garnier et al. (2001) that this sensor only measures the liquid temperature even for boiling conditions since it acts as a low-pass filter. This can be explained by the size of the sensor which leads to a large characteristic time compared to the characteristic scales of phase intermittency in the boiling flow.

To ensure good accuracy of these measurements, a great attention was paid to the calibration process. During this calibration process, adiabatic and isothermal liquid single-phase flow is imposed in the test section. Reference temperature is then obtained by considering inlet ant outlet temperatures that are measured using Platinum probes. The resulting accuracy of the liquid temperature measurement is ±0.2 °C. A specific campaign (C800) of liquid temperature measurements (Garnier et al., 2001) include single-phase flow conditions along with boiling ones.

### 2.5.2. Optical probes

Two campaigns have been performed, each of them using a specific sensor:

* Campaign C2900 using a single optical probe,

* Campaign C3000 using a double optical probe.

An optical probe is a sensor that measures the Phase Indicator Function (PIF). It acts as a light guide. A laser diode generates the light transmitted through the optical fibre. By means of a Y coupler, the reflected light is deviated to a photon converter. Then, an electronic amplifier transforms the intensity at the outlet of this converter into a high level voltage signal. Generally, the amplifier gain and offset are set in order to get no voltage for the liquid level and about 5 V for the vapour level. To obtain a signal close to a PIF, the amplifier bandwidth is higher than 1 MHz. Deviations from the theoretical binary signal are then only due to the interaction between the tip and the interfaces. Since recording the whole raw data was not possible, the amount of data to store being too large, the signal was processed in real time after approximately 150 s of acquisition (François et al., 2003).

For a single optical probe, such a treatment leads to the calculation of the following variables:

* the local void fraction $\alpha$ defined as:

$$ \alpha = \frac{1}{T_c} \sum_{i} t_{vi} \qquad (5) $$

where $T_c$ is the integration time and $t_{vi}$ is the vapour residence time corresponding to the $i$th bubble;

* the interference frequency $\nu$ defined as:

$$ \nu = \frac{N_p}{T_c} \qquad (6) $$

where $N_p$ is the number of bubbles that interfered with the sensor.

By assuming a model for the vapour velocity, and using stereological tools<fnref n="2" /> Garnier et al. (2001), it is then possible to access additional variables such as Sauter mean diameter $d_{32}$ or Interfacial Area Concentration (IAC, $a_i$). But one should notice that those quantities result from a modelling assumption. Their validity is then questionable and may strongly depend on the chosen model for the vapour velocity. Garnier et al. (2001) suggest to use a power law for the vapour velocity ($1/\gamma^{th}$) assuming a local homogeneous flow ($v_\nu(r) = v_i(r)$). This assumption on velocities also affects the diameter assessment; it can induce up to 30% variations locally on some cases (especially for low mass fluxes and low pressure tests).

<footnote n="2">² Some additional assumptions are necessary and they are fully described in Garnier et al. (2001).</footnote>
On the other hand, a double optical probe allows to determine,
under other assumptions, the vapour velocity and then to deduce addi-
tional parameters like Sauter mean diameter or IAC. Those assumptions
are widely discussed in Garnier et al. (2001). The vapour velocity
uncertainty is estimated by Garnier et al. (2001) to a value close
to ±10% whereas the uncertainty for Sauter mean diameter is ±20%.
Besides, a few liquid velocity measurements have been performed on
DEBORA later, between 2016 and 2018, before the loop was disman-
tled. Those measurements were performed using thermal-anemometry
and for more details, we suggest the reader to refer to the following
papers (François et al., 2020, 2021).

Consequently, it was decided for this benchmark to use only void
fraction, Sauter mean diameter and IAC as reference measurements for
campaign C2900 whereas axial vapour velocities were also included
from campaign C3000.

Moreover, for each kind of optical probe, in order to limit the time
duration of a test, experimentalists considered that the flow becomes
liquid when the void fraction gets lower than 0.1%; the probe motion
towards the pipe centre is then stopped without further measurements.

## 2.6. Available conditions

Roughly 50 cases are partly available for the benchmark from Garnier et al. (2001), many of which differ only by the inlet temperature. Some experimental conditions have been reproduced in different experimental campaigns to benefit from other measurement devices and provide additional information such as direct measurements of vapour velocity or liquid temperature profiles. These conditions are reported in Table 4. The work of Garnier et al. (2001) was widely used to validate MCFD codes, even though they only present selected variables that vary from one case to another. Here, the full available measurements are given for each selected case. Experimental and numerical data are reported in Bois et al. (2024).

## 2.7. Selected benchmarking conditions

The experimental conditions selected for the comparison are pre-
sented here. Their distributions are graphically represented on Fig. 3
in the plane pressure-mass flux (P, G). For the results analysis in
Section 5, cases are gathered in four groups G2P14W16, G2P26W23,
G3P26W16, G5P14 that can be differentiated on the plane of Fig. 3.
Note that the electrical power W is constant within each of the
three first groups, but evolves with some cases in the fourth group
G5P14. For the experimental conditions, each test-case is defined by
the set of parameters (G, P, qw, Tin) consisting of mass flux, outlet
pressure, wall heat flux and inlet temperature. We use a nomenclature
as 29G2P14W16Te25.3 where 29 refers to the campaign number
(29 with the single optical probe, 30 with the double-tipped optical
probe or 8 for thermal measurements), 2 is the truncated flow rate in
10³ kg m⁻² s⁻¹, 14 the truncated pressure in bar, 16 the total electrical
power supplied to the loop in kW and 25.3 the inlet temperature in °C.

The benchmark was performed in two steps whose conditions are
described separately in Tables 5 and 6. In the first phase, several
conditions are covered :

* From G2P14W16 (G, P, W) = (2000 kg m⁻² s⁻¹, 14.6 bar, 16 kW), four cases were studied in detail to describe the void distribution as the quality increases with the increase in inlet temperature: T<sub>e</sub> = 25.3 °C, 30 °C, 38.8 °C and 40 °C (29G2P14W16Te25.3, 30G2P14W16Te25.3, 29G2P14W16Te30, 30G2P14W16Te30, 30G2P14W16Te38.8, 29G2P14W16Te38.8, 30G2P14W16Te40).<fnref n="3" />

* 2 cases were provided at a higher pressure and flow rate to gather knowledge on the different evolution at a higher pressure P = 26 bar (G2P26W23).

* 4 cases with higher flow rate and heat-fluxes at a pressure P = 14.6 bar were also considered in G5P14 (29G5P14W29Te33.9, 29G5P14W29Te39.8, 29G5P14W29Te43.4, 29G5P14W29Te46). They illustrate the flow at similar thermodynamic equilibrium qualities but with different boiling numbers (Bo = qw/(GLvap), defined as in Eq. (3) by the heat flux intensity relatively to the flow rate).

In addition, the second phase of the benchmark included some
other open conditions from Garnier et al. (2001) with a particular
focus on liquid temperature measurements. Besides, new cases were se-
lected from the database and the experimental measurements were only
provided to the participants after collecting their results (conditions
prefixed B in Table 6 and referred to as *blind cases*).<fnref n="4" />

Overall, the benchmark covers challenging conditions with high quality at the exit (up to $x = 0.088$) and local void fraction above 50%. Due to the high pressure, bubbles remain small, and there is no transition from bubbly flow regime. This behaviour has been experimentally highlighted by Hösler (1968). He performed some visualisations of steam/water flows in a rectangular channel with one side electrically heated. He observed that for pressure above 70 bar, the flow tended to remain a bubbly flow whatever the value of the void fraction, whereas for lower pressure, transition from bubbly flow to slug then annular flow could be observed. Later, François et al. (2011) noticed a similar behaviour in Freon experiments where no flow transition from bubbly flow was observed even for void fractions as high as 90%.

The participants goal was to predict most accurately the devel-
opment of the boiling flow in the pipe for each set of experimental
conditions. In time, we hope to make this databank (Bois et al., 2024) a
part of a more complete one, useable for the validation and comparison
of CMFD codes.

## 3. Numerical setup and boundary conditions

This part describes the numerical setup and the conventions used to
describe the flow.

### 3.1. Simulation domain and boundary conditions

The simulation is performed on a portion of vertical pipe of internal radius $R = 9.6$ mm. The origin of the cylindrical coordinates system $(r, \theta, z)$ is taken at the pipe centre, at the beginning of the heated section, with the z-axis oriented upwards (see Fig. 4). Participants have discretised either a small wedge for 2D axisymmetric simulations or a 3D domain representing a quarter or a full pipe. Symmetry planes are used to cut the domain or at the central axis.

Test conditions are described by 4 parameters (G, P, W, T<sub>e</sub>) defined in Tables 5 and 6 for each benchmark phase respectively. Constant pressure P is set at the outlet boundary (z = 4.5 m). Constant heat flux qw (computed from the case total power W) is imposed at the wall r = R in the heated section for z ∈ [0; H], with the heated length H = 3.5 m. The axial domain is prolonged 1 m upstream and 0.5 m downstream of the heated length by adiabatic walls. Some mesh characteristics and resolution are given in Table 7. The inlet boundary

of prescribed loop conditions due to regulation (i.e., they result from different experimental campaigns). However, the tests are complementary as they use different radial resolutions and the tests in C30 provide explicit measurement of the vapour velocity.

<footnote n="3">³ The two cases 29G2P14W16Te25.3 and 30G2P14W16Te25.3 (or similarly for the pairs at Te30 or at Te38.8) are almost identical, up to the variations</footnote>

<footnote n="4">⁴ CEA had the experimental results before the organisation of the benchmark, but decided to use the reference Neptune_CFD models without any particular adjustment. In that sense, the CEA results can also be considered as blind (for both phases actually).</footnote>
<caption>**Table 4**</caption>

List of conditions and measurements presented in Garnier et al. (2001): x quality at the exit of the heated section, α: void fraction, d<sub>b</sub> = d<sub>10</sub>: bubble arithmetic mean diameter, as defined by Eq. (4.13) in Garnier et al. (2001), a<sub>i</sub>: interfacial area, d<sub>32</sub>: Sauter mean diameter, bcd: bubble centre density, u<sup>v</sup>: axial vapour velocity, T<sup>i</sup>: liquid temperature.

<table><thead><tr><th>#</th><th>Quality x</th><th>Case ref.</th><th>α</th><th>d<sub>b</sub></th><th>a<sub>i</sub></th><th>d<sub>32</sub></th><th>bcd</th><th>u<sup>v</sup></th><th>T<sup>i</sup></th></tr></thead><tbody><tr><td>1</td><td>-0.2298</td><td>29G1P30W14Te37.7</td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>2</td><td>-0.1612</td><td>29G1P30W14Te42.9</td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>3</td><td>-0.118</td><td>29G1P30W14Te46.3</td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>4</td><td>-0.0835</td><td>29G1P30W14Te48.8</td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>5</td><td>-0.0547</td><td>29G1P30W14Te50.8</td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>6</td><td>-0.032</td><td>29G1P30W14Te52.5</td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>7</td><td>-0.2165</td><td>29G1P30W12Te44</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>8</td><td>-0.0973</td><td>29G1P30W12Te52.7</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>9</td><td>0.0585</td><td>29G1P30W12Te63.6</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>10</td><td>0.1343</td><td>29G1P30W12Te68.9</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>11</td><td>0.2173</td><td>29G1P30W12Te74</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>12</td><td>-0.1127</td><td>8G2P14W16Te19.2</td><td></td><td></td><td></td><td></td><td></td><td></td><td>x</td></tr><tr><td>13</td><td>-0.0677</td><td>29G2P14W16Te22.6</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>14</td><td>-0.0043/-0.0431</td><td>30G2P14W16Te25.3</td><td>x</td><td></td><td>x</td><td></td><td></td><td>x</td><td></td></tr><tr><td>15</td><td>-0.0597</td><td>8G2P14W16Te25.6</td><td></td><td></td><td></td><td></td><td></td><td></td><td>x</td></tr><tr><td>16</td><td>-0.0297</td><td>29G2P14W16Te27.3</td><td>x</td><td>x</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>17</td><td>-0.0185</td><td>29G2P14W16Te28.7</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>18</td><td>-0.0071/-0.0043</td><td>29G2P14W16Te30</td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>19</td><td>-0.0043</td><td>30G2P14W16Te30</td><td>x</td><td></td><td>x</td><td></td><td></td><td>x</td><td></td></tr><tr><td>20</td><td>0.0014</td><td>29G2P14W16Te30.8</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>21</td><td>-0.0102</td><td>8G2P14W16Te31.4</td><td></td><td></td><td></td><td></td><td></td><td></td><td>x</td></tr><tr><td>22</td><td>0.0086</td><td>29G2P14W16Te31.8</td><td></td><td>x</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>23</td><td>0.0261</td><td>29G2P14W16Te33.8</td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>24</td><td>0.0319</td><td>29G2P14W16Te34.5</td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>25</td><td>0.0392</td><td>29G2P14W16Te35.3</td><td></td><td>x</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>26</td><td>0.046</td><td>29G2P14W16Te36</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>27</td><td>0.0596</td><td>29G2P14W16Te37.7</td><td>x</td><td>x</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>28</td><td>0.0687</td><td>29G2P14W16Te38.8</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>29</td><td>0.077</td><td>30G2P14W16Te38.8</td><td>x</td><td></td><td>x</td><td></td><td></td><td>x</td><td></td></tr><tr><td>30</td><td>0.0881</td><td>30G2P14W16Te40</td><td>x</td><td></td><td>x</td><td></td><td></td><td>x</td><td></td></tr><tr><td>31</td><td>0.1091</td><td>29G2P14W16Te43.5</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>32</td><td>-0.5308</td><td>8G2P26W16Te19.8</td><td></td><td></td><td></td><td></td><td></td><td></td><td>x</td></tr><tr><td>33</td><td>-0.4011</td><td>8G2P26W16Te31.5</td><td></td><td></td><td></td><td></td><td></td><td></td><td>x</td></tr><tr><td>34</td><td>-0.1851</td><td>8G2P26W16Te49.7</td><td></td><td></td><td></td><td></td><td></td><td></td><td>x</td></tr><tr><td>35</td><td>-0.0719</td><td>29G2P26W16Te58</td><td></td><td></td><td></td><td></td><td></td><td>x</td><td></td></tr><tr><td>36</td><td>-0.0463</td><td>29G2P26W16Te60</td><td></td><td></td><td></td><td></td><td></td><td>x</td><td></td></tr><tr><td>37</td><td>-0.0205</td><td>29G2P26W16Te62</td><td></td><td></td><td></td><td></td><td></td><td>x</td><td></td></tr><tr><td>38</td><td>0.006</td><td>29G2P26W16Te64</td><td></td><td></td><td></td><td></td><td></td><td>x</td><td></td></tr><tr><td>39</td><td>0.0324</td><td>29G2P26W16Te66</td><td></td><td></td><td></td><td></td><td></td><td>x</td><td></td></tr><tr><td>40</td><td>0.058</td><td>29G2P26W16Te68</td><td></td><td></td><td></td><td></td><td></td><td>x</td><td></td></tr><tr><td>41</td><td>0.0848</td><td>29G2P26W16Te70</td><td></td><td></td><td></td><td></td><td></td><td>x</td><td></td></tr><tr><td>42</td><td>0.0806</td><td>8G2P26W16Te70.3</td><td></td><td></td><td></td><td></td><td></td><td></td><td>x</td></tr><tr><td>43</td><td>-0.0519</td><td>29G3P26W23Te60.9</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>44</td><td>-0.0177</td><td>29G3P26W23Te63.2</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>45</td><td>0.0164</td><td>29G3P26W23Te65.7</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>46</td><td>0.0479</td><td>29G3P26W23Te68</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>47</td><td>0.077</td><td>29G3P26W23Te70</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>48</td><td>0.1005</td><td>29G3P26W23Te72</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>49</td><td>-0.0417</td><td>29G5P14W29Te33.9</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>50</td><td>0.0067</td><td>29G5P14W29Te39.8</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>51</td><td>0.0398</td><td>29G5P14W29Te43.4</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>52</td><td>0.0625</td><td>29G5P14W29Te46</td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr></tbody></table>

condition is defined by the flow rate G and inlet uniform temperature T<sub>e</sub> either based on a uniform flat liquid velocity profile at z = -1 m or using recirculations or predefined velocity profile at z = 0 m depending on the participant. Standard gravity g = -9.81 m s<sup>-2</sup> is applied along z.

Measurements were performed along r direction at z<sub>m</sub> = 3.485 m. Dimensionless coordinates are defined by r<sup>+</sup> = r/R and z<sup>+</sup> = z/H. Numerical predictions are assessed along these lines and near-wall streamwise evolution is also illustrated in Electronic Appendix for r<sup>+</sup> = 0.9375. Cross-section averaged void-fraction is also studied along the pipe streamwise coordinate z; it is defined by:

$$ \langle \alpha \rangle = \frac{2}{R^2} \int_{0}^{R} \alpha(r) r dr \qquad (7) $$

## 3.2. Fluid properties

The physical properties of the two phases of the R12 (Dichlorodifluoromethane) fluid used in the experiments are estimated independently by each participant (e.g., Bell et al., 2014; Lemmon et al., 2023) according to their tools capabilities. Some elements are given by each contributor in Appendix A and the effects of differences are discussed in Sections 6 and 5.3.

## 3.3. Figures of merit

There are many challenges in the prediction of high-pressure boiling flows. General interest concerns the local thermodynamical unbalance characterising the local subcooling and the intensity of void-fraction migration away from the wall to the core flow. They all govern the flow disequilibrium that is so difficult to capture in boiling flows. Local
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Physical conditions used for the benchmark compared to the conditions in Garnier et al. (2001). Coloured numbers indicate the number of conditions available either in Garnier et al. (2001) or in each phase from Tables 5 and 6.</figcaption></figure>

<caption>**Table 5**</caption>

List of conditions and measurements selected for the first phase of the benchmark. Profiles in Garnier et al. (2001) are represented by x, new profiles provided for the benchmark are noted o.

<table><caption>**Table 6**</caption><thead><tr><th>#</th><th>Case ref.</th><th>G<br>[kg m<sup>-2</sup> s<sup>-1</sup>]</th><th>P<br>[bar]</th><th>q<sub>w</sub><br>[kW m<sup>-2</sup>]</th><th>W<br>[kW]</th><th>T<sub>c</sub><br>[°C]</th><th>Quality<br>x [-]</th><th>α<br>[-]</th><th>a<sub>i</sub><br>[m<sup>-1</sup>]</th><th>u<sup>v</sup><br>[m s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>1</td><td>30G2P14W16Te25.3</td><td>2024.5</td><td>14.60</td><td>76.69</td><td>16.2</td><td>26.8</td><td>-0.0431</td><td>x</td><td>x</td><td>x</td></tr><tr><td>2</td><td>29G2P14W16Te30</td><td>2030.0</td><td>14.59</td><td>76.24</td><td>16.1</td><td>31.2</td><td>-0.0071</td><td>x</td><td>o</td><td></td></tr><tr><td>3</td><td>30G2P14W16Te30</td><td>2022.1</td><td>14.60</td><td>76.25</td><td>16.1</td><td>31.4</td><td>-0.00434</td><td>x</td><td>x</td><td>x</td></tr><tr><td>4</td><td>29G2P14W16Te38.8</td><td>2022.9</td><td>14.59</td><td>76.26</td><td>16.1</td><td>39.7</td><td>0.0687</td><td>x</td><td>x</td><td></td></tr><tr><td>5</td><td>30G2P14W16Te38.8</td><td>2007.7</td><td>14.59</td><td>77.68</td><td>16.4</td><td>40.0</td><td>0.0776</td><td>x</td><td>x</td><td>x</td></tr><tr><td>6</td><td>29G2P14W16Te40</td><td>2022.9</td><td>14.59</td><td>76.26</td><td>16.1</td><td>40.8</td><td>0.0788</td><td>o</td><td>o</td><td></td></tr><tr><td>7</td><td>30G2P14W16Te40</td><td>2007.8</td><td>14.59</td><td>77.69</td><td>16.4</td><td>41.2</td><td>0.0882</td><td>x</td><td>x</td><td>x</td></tr><tr><td>8</td><td>29G2P14W16Te43.5</td><td>2024.1</td><td>14.59</td><td>76.26</td><td>16.1</td><td>44.2</td><td>0.1091</td><td>x</td><td>x</td><td></td></tr><tr><td>9</td><td>29G3P26W23Te60.9</td><td>2992.6</td><td>26.16</td><td>107.52</td><td>22.7</td><td>60.9</td><td>-0.0519</td><td>x</td><td>x</td><td></td></tr><tr><td>10</td><td>29G3P26W23Te72</td><td>2984.3</td><td>26.17</td><td>107.52</td><td>22.7</td><td>72.5</td><td>0.1005</td><td>x</td><td>x</td><td></td></tr><tr><td>11</td><td>29G5P14W29Te33.9</td><td>5063.2</td><td>14.57</td><td>135</td><td>28.5</td><td>35.0</td><td>-0.0417</td><td>x</td><td>x</td><td></td></tr><tr><td>12</td><td>29G5P14W29Te39.8</td><td>5085.2</td><td>14.56</td><td>135</td><td>28.5</td><td>40.5</td><td>0.0067</td><td>x</td><td>x</td><td></td></tr><tr><td>13</td><td>29G5P14W29Te43.4</td><td>5063.0</td><td>14.56</td><td>135</td><td>28.5</td><td>44.1</td><td>0.0398</td><td>x</td><td>x</td><td></td></tr><tr><td>14</td><td>29G5P14W29Te46</td><td>5070.3</td><td>14.57</td><td>135</td><td>28.5</td><td>46.7</td><td>0.0625</td><td>x</td><td>x</td><td></td></tr></tbody></table>

List of operating conditions and measurements selected for the second phase of the benchmark. Profiles in Garnier et al. (2001) are represented by x, new profiles provided for the benchmark are noted o. Outlet quality is also assessed. For cases labelled B#, experimental measurements were provided after the participant's submissions (blind cases). Cases O22 and O25 are identical to cases 9 and 10 from phase 1 in Table 5.

<table><caption>**Table 7**</caption><thead><tr><th>#</th><th>Case ref.</th><th>G<br>[kg m<sup>-2</sup> s<sup>-1</sup>]</th><th>P<br>[bar]</th><th>q<sub>w</sub><br>[kW m<sup>-2</sup>]</th><th>W<br>[kW]</th><th>T<sub>c</sub><br>[°C]</th><th>Quality<br>x [-]</th><th>α<br>[-]</th><th>a<sub>i</sub><br>[m<sup>-1</sup>]</th><th>u<sup>v</sup><br>[m s<sup>-1</sup>]</th><th>T<sup>l</sup><br>[°C]</th></tr></thead><tbody><tr><td>O15</td><td>29G2P14W16Te30.8</td><td>2014.4</td><td>14.58</td><td>76.26</td><td>16.1</td><td>31.9</td><td>0.0014</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>O16</td><td>29G2P14W16Te31.8</td><td>2029.6</td><td>14.58</td><td>76.26</td><td>16.1</td><td>32.9</td><td>0.0086</td><td>o</td><td>o</td><td></td><td></td></tr><tr><td>B17</td><td>30G2P26W16Te62</td><td>2051.8</td><td>26.20</td><td>73.89</td><td>15.6</td><td>62.5</td><td>-0.0318</td><td>o</td><td>o</td><td>o</td><td></td></tr><tr><td>B18</td><td>30G2P26W16Te64</td><td>2050.4</td><td>26.20</td><td>73.89</td><td>15.6</td><td>64.6</td><td>-0.0050</td><td>o</td><td>o</td><td>o</td><td></td></tr><tr><td>B19</td><td>30G2P26W16Te66</td><td>2049.0</td><td>26.20</td><td>73.89</td><td>15.6</td><td>66.6</td><td>0.0213</td><td>o</td><td>o</td><td>o</td><td></td></tr><tr><td>B20</td><td>30G2P26W16Te68</td><td>2049.2</td><td>26.20</td><td>73.89</td><td>15.6</td><td>68.6</td><td>0.0473</td><td>o</td><td>o</td><td>o</td><td></td></tr><tr><td>B21</td><td>30G2P26W16Te70</td><td>2051.2</td><td>26.19</td><td>73.89</td><td>15.6</td><td>70.6</td><td>0.0743</td><td>o</td><td>o</td><td>o</td><td></td></tr><tr><td>O22</td><td>29G3P26W23Te60.9</td><td>2992.6</td><td>26.16</td><td>107.52</td><td>22.7</td><td>60.9</td><td>-0.0519</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>B23</td><td>29G3P26W23Te67.6</td><td>2998.1</td><td>26.19</td><td>109.11</td><td>23.0</td><td>68.0</td><td>0.0434</td><td>o</td><td>o</td><td>o</td><td></td></tr><tr><td>B24</td><td>30G3P26W23Te69</td><td>2986.9</td><td>26.18</td><td>109.00</td><td>23.0</td><td>69.5</td><td>0.0643</td><td>o</td><td>o</td><td>o</td><td></td></tr><tr><td>O25</td><td>29G3P26W23Te72</td><td>2984.3</td><td>26.17</td><td>107.52</td><td>22.7</td><td>72.49</td><td>0.1005</td><td>x</td><td>x</td><td></td><td></td></tr><tr><td>B26</td><td>29G5P14W42Te30.2</td><td>5066.8</td><td>14.58</td><td>197.12</td><td>41.6</td><td>31.18</td><td>0.0019</td><td>o</td><td>o</td><td></td><td></td></tr><tr><td>B27</td><td>29G5P14W33Te41</td><td>4994.9</td><td>14.58</td><td>152.83</td><td>32.3</td><td>41.64</td><td>0.0417</td><td>o</td><td>o</td><td></td><td></td></tr><tr><td>T28</td><td>8G2P14W16Te25.6</td><td>2009</td><td>14.6</td><td>73.9</td><td>15.6</td><td>25.65</td><td>-0.0597</td><td></td><td></td><td></td><td>x</td></tr><tr><td>T29</td><td>8G2P14W16Te31.4</td><td>2007</td><td>14.6</td><td>73.9</td><td>15.6</td><td>31.37</td><td>-0.0102</td><td></td><td></td><td></td><td>x</td></tr><tr><td>T30</td><td>8G2P26W16Te19.8</td><td>1996</td><td>26.1</td><td>74.4</td><td>15.7</td><td>19.84</td><td>-0.5308</td><td></td><td></td><td></td><td>x</td></tr><tr><td>T31</td><td>8G2P26W16Te31.5</td><td>1995</td><td>26.2</td><td>74.4</td><td>15.7</td><td>31.46</td><td>-0.4011</td><td></td><td></td><td></td><td>x</td></tr><tr><td>T32</td><td>8G2P26W16Te49.7</td><td>2000</td><td>26.2</td><td>73.9</td><td>15.6</td><td>49.72</td><td>-0.1851</td><td></td><td></td><td></td><td>x</td></tr><tr><td>T33</td><td>8G2P26W16Te70.3</td><td>1983</td><td>26.2</td><td>73.9</td><td>15.6</td><td>70.31</td><td>0.0806</td><td></td><td></td><td></td><td>x</td></tr></tbody></table>

Meshes characteristics. †: radial grading (step 0.5:1).

<table><thead><tr><th></th><th>ANSYS</th><th>CEA</th><th>FRA</th><th>HZDR</th><th>JSI</th><th>KTH-WSE</th><th>LEI</th><th>MIT</th><th>PSI-CTU</th><th>SJTU</th><th>UJV</th><th>VTT</th></tr></thead><tbody><tr><th>Dim.</th><td>2</td><td>2</td><td>3</td><td>2</td><td>2</td><td>2</td><td>2</td><td>3 (1/4)</td><td>2</td><td>2</td><td>2</td><td>2</td></tr><tr><th>N<sub>r</sub></th><td>40</td><td>40</td><td>≈ 22</td><td>40</td><td>15</td><td>40<sup>†</sup></td><td>40</td><td>≈16</td><td>40</td><td>24</td><td>26</td><td>40</td></tr><tr><th>N<sub>z</sub></th><td>1000</td><td>768</td><td>≈ 7000</td><td>400</td><td>500</td><td>400</td><td>500</td><td>1168</td><td>2000</td><td>200</td><td>2000</td><td>500</td></tr><tr><th>N<sub>tot</sub></th><td>40k</td><td>30.8k</td><td>19M</td><td>16k</td><td>7.5k</td><td>16k</td><td>20k</td><td>195k</td><td>80k</td><td>4.8k</td><td>56k</td><td>20k</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Fig. 4. Description of the computational domain and coordinate definition. Source: Adapted from Yao and Morel (2004).</figcaption></figure>

bubble sizes distribution is then a key factor in both momentum and heat interfacial transfers. A particular attention was set on several flow characteristics at the end of the heated section

* Maximum void fraction and its radial location,

* Near-wall void fraction value,

* Mean bubble diameter radial distribution,

* Liquid and vapour streamwise and radial velocities,

* Liquid temperature profiles,

* Turbulent kinetic energy,

or along the flow main axis:

* Area-averaged void fraction,

* Wall temperature,

* Evaporating heat flux.

All these quantities are provided in Electronic Appendix for each participant and compared to experimental measurements whenever available. Data is also provided as supplementary material. These profiles give a very thorough description of the flow modelled by each participant; they can support further analyses to comprehend the models complex interplay. Some of these flow characteristics are essential in nuclear safety as they are key elements related to the Critical Heat Flux (CHF) prediction.

## 4. Models

The aim of this section is to provide a thorough and accurate description of the physical closures used to perform the tests. Each participant used a single set of closures for the whole database as the goal of the benchmark is to progress towards a unified set of closures for MCFD models. First, Section 4.1 describes the global equation system and the main closure relations associated to it. Tables summarise the closures for each participant and more details can be found in dedicated Appendix A. Section 4.2 comes back on the most common choices, on some striking differences and on some important parameters to better emphasise the variety of approaches and the importance of calibrated parameters.

### 4.1. Brief summary of models

We describe here the general features of the models. All the participants rely on a Eulerian–Eulerian RANS approach with a single pressure. The equation system is established to govern phase-averaged quantities $\overline{\psi}_k^k = \chi_k \psi_k / \overline{\chi}_k$, where the statistical average $\overline{\psi}$ (with the classical properties of a Reynolds averaging) is weighted by the average $\overline{\chi}_k$ of the phase indicator function $\chi_k$ (where the subscript $k \in [l, v]$ stands for the liquid or vapour phase, respectively). The void fraction is then defined as $\alpha_v = \overline{\chi}_v$. The equation system relies on two sets of statistically-averaged conservation equations governing mass, momentum and enthalpy in each phase (Ishii and Hibiki, 2015; Morel, 2015) for $k \in [l, v]$:

$$ \frac{D_k (\alpha_k \rho_k)}{Dt} = \Gamma_k + \Gamma_k^w, \qquad (8a) $$

$$ \frac{D_k (\alpha_k \rho_k \bar{\mathbf{u}}_k^k)}{Dt} = -\alpha_k \nabla \bar{\rho}_l^l + \nabla \cdot (\mathcal{T}_k^v - \mathcal{T}_k^f) + \alpha_k \rho_k \mathbf{g} + \mathbf{M}_{l \to k}, \qquad (8b) $$

$$ \frac{D_k (\alpha_k \rho_k \bar{h}_k^k)}{Dt} = \frac{\partial \alpha_k \bar{\rho}_l^l}{\partial t} + \nabla \cdot (\alpha_k \lambda_k \nabla \bar{T}_k^k - \mathbf{q}_k^l) + q_k^l a_i + q_k^w, \qquad (8c) $$

where $D_k (\alpha_k \rho_k \bar{\phi}_k^k) / Dt = \partial (\alpha_k \rho_k \bar{\phi}_k^k) / \partial t + \nabla \cdot (\alpha_k \rho_k \bar{\mathbf{u}}_k^k \bar{\phi}_k^k)$ is the material or substantial derivative of any quantity $\bar{\phi}_k^k$ associated to the phase velocity $\bar{\mathbf{u}}_k^k$. $\mathcal{T}_k^v$ is the mean viscous stress approximated by $\mathcal{T}_k^v = \alpha_k \mu_k ( (\nabla \bar{\mathbf{u}}_k^k) + (\nabla \bar{\mathbf{u}}_k^k)^l )$. Here, $\mathbf{u}$, $\rho$, $h$, $T$, $\rho$, $\mu$ and $\lambda$ denote respectively the velocity, pressure, enthalpy, temperature, density, dynamic viscosity and conductivity. $\mathbf{g}$ is the gravitational acceleration. The sum of hydrostatic and dynamic pressures $p$ is considered in the liquid phase and local pressure equilibrium is assumed. Equations of state (EOS) are required for each phase to assess the latent heat of vaporisation $L^{\text{vap}}$, the surface tension coefficient $\sigma$ and phase properties (densities, viscosities, conductivities and heat capacities, of the form: $\rho_k = f(\bar{\rho}_k^k, \bar{\chi}_k^k), \dots$).

In addition to these sets of governing equations and EOS, several closure relations are required to fully describe the system. For the momentum transport equation, closure relations concern the turbulent stresses $\mathcal{T}_k^l = \alpha_k \rho_k \bar{\mathbf{u}}_k^l \bar{\mathbf{u}}_k^l$ in phase $k$ and the interfacial momentum transfer $\mathbf{M}_{l \to k}$. Transport equations are involved to predict turbulent quantities (turbulent stress $\mathcal{T}_k^l$, or turbulent kinetic energy $k_k = \text{tr}(\mathcal{T}_k^l)$ and eddy dissipation). Interfacial momentum transfer is modelled as a series of forces satisfying $\sum_k \mathbf{M}_{l \to k} = 0$. These forces classically include drag, added-mass, lift, turbulent dispersion and possibly wall repulsion and are given in Table 8. The lift force is expressed as follows

$$ \mathbf{F}_L = -C_L F_{\text{damp}} \rho_l \alpha_v (\bar{\mathbf{u}}_l^l - \bar{\mathbf{u}}_v^l) \times (\nabla \times \bar{\mathbf{u}}_l^l), \qquad (9) $$

with the possibility to consider a damping function $F_{\text{damp}}(y)$ near the wall.

For the phase mass conservation, the interfacial mass transfer $\Gamma_k$ satisfies the jump relation $\Gamma_l + \Gamma_v = 0$ ($\Gamma_v > 0$ is the interfacial mass flux due to evaporation). It is related to the interfacial heat transfer $q_k^l$ by $a_i \sum_k q_k^l = \Gamma_v L^{\text{vap}}$ where $a_i$ is the Interfacial Area Concentration (IAC). Interfacial heat fluxes $q_k^l = \lambda_l N u_k a_i (T^{\text{sat}} - T_k) / d_{32}$ are modelled by correlations for Nusselt numbers in each phase given in Table 10; they are then used to compute the interfacial mass flux $\Gamma_v$.

For the energy conservation, turbulent heat fluxes $q_k^l = \alpha_k \rho_k \bar{\mathbf{u}}_k^l \bar{\mathbf{h}}_k^l$ are also required; they are connected to turbulent stresses by the Standard or Generalised Gradient Diffusion Hypothesis (SGDH or GGDH) as given in Table 10. In addition, at the heated wall, a wall heat flux $q_k^w$ into each phase is accounted for. It results from the Wall Heat Partitioning (WHP) model given in Table 11 that splits the total wall flux $q_w$ known at a boundary condition into several parts

$$ q_w = q_l^w + q_v^w + q_e, \qquad (10) $$

where $q_e$ is the evaporation heat flux, $q_v^w$ the direct flux going to the vapour at the wall and $q_l^w = q_{SP} + q_q$ the wall heat flux heating the
<table><caption>**Table 9**
Interfacial area modelling.</caption><thead><tr><th colspan="13">INTERFACIAL AREA MODELLING</th></tr><tr><th>Model</th><th>ANSYS</th><th>CEA</th><th>FRA</th><th>HZDR</th><th>JSI</th><th>KTH-WSE</th><th>LEI</th><th>MIT</th><th>PSI-CTU</th><th>SJTU</th><th>UJV</th><th>VTT</th></tr></thead><tbody><tr><td>Type</td><td>MUSIG</td><td>IATE<br>RS07</td><td>IATE<br>YM04</td><td>iMUSIG<br>Kr13</td><td>MUSIG<br>sym.<br>24</td><td>MUSIG<br>23</td><td>MUSIG<br>27</td><td>S-γ<br>Va98<br>n.a.</td><td>MUSIG<br>20</td><td>iMUSIG<br>25</td><td>iMUSIG<br>10</td><td>iMUSIG<br>30</td></tr><tr><td>Diam.<br>grp.</td><td>15</td><td>n.a.</td><td>n.a.</td><td>20</td><td>24</td><td>23</td><td>27</td><td>n.a.</td><td>20</td><td>25</td><td>10</td><td>30</td></tr><tr><td>Velo.<br>grp.</td><td>1</td><td>n.a.</td><td>n.a.</td><td>2</td><td>1</td><td>1</td><td>1</td><td>n.a.</td><td>1</td><td>2</td><td>2</td><td>2</td></tr><tr><td>Coalescence</td><td>PB90</td><td>PB90<sup>†</sup></td><td>YM04</td><td>PB90</td><td>Li15</td><td>PB90</td><td>PB90</td><td>LZ09</td><td>PB90</td><td>Guo16</td><td>PB90</td><td>PB90</td></tr><tr><td>F<sub>CB</sub></td><td>0.02</td><td>0.0</td><td>0.0</td><td>0</td><td>n.a.</td><td>1.0</td><td>1</td><td>n.a.</td><td>0.2</td><td>n.a.</td><td>4.0</td><td>1.0</td></tr><tr><td>F<sub>CT</sub></td><td>0.02</td><td>1.0</td><td>C<sub>mk</sub> = 2.86</td><td>0.1</td><td>n.a.</td><td>2.014</td><td>4.53</td><td>n.a.</td><td>0.2</td><td>n.a.</td><td>4.0</td><td>4.028</td></tr><tr><td>F<sub>LS</sub></td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>n.a.</td><td>1.0</td><td>1</td><td>n.a.</td><td>0.0</td><td>n.a.</td><td>0.0</td><td>1.0</td></tr><tr><td>Break-up</td><td>LS96</td><td>Wu98</td><td>YM04</td><td>LS96</td><td>Li15</td><td>Le02</td><td>Le02</td><td>LZ09</td><td>LS96</td><td>LS96</td><td>LS96</td><td>Le02</td></tr><tr><td>F<sub>b</sub></td><td>0.5</td><td>n.a.</td><td>C<sub>mk</sub> = 1.6</td><td>0.</td><td>1.0</td><td>n.a.</td><td>n.a.</td><td>n.a.</td><td>0.5</td><td>0.5</td><td>1.0</td><td>n.a.</td></tr><tr><td colspan="13">Guo16 (Guo et al., 2016)</td></tr><tr><td colspan="13">LZ09 (Lo and Zhang, 2009)</td></tr><tr><td colspan="13">Va98 (Vaesson, 1998)</td></tr><tr><td colspan="13">Kr13 (Krepper et al., 2013)</td></tr><tr><td colspan="13">LS96 (Luo and Svendsen, 1996)</td></tr><tr><td colspan="13">Wu98 (Wu et al., 1998)</td></tr><tr><td colspan="13">Le02 (Lehr et al., 2002)</td></tr><tr><td colspan="13">PB90 (Prince and Blanch, 1990)</td></tr><tr><td colspan="13">YM04 (Yao and Morel, 2004)</td></tr></tbody></table>

n.a.: Not applicable; $F_{CB}$: Buoyancy coalescence coefficient; $F_{CT}$: Turbulence coalescence coefficient; $F_{LS}$: Laminar shear coalescence coefficient; Coalescence coefficients are introduced as fitting parameters in Eq. (21) of Prince and Blanch (1990) as

$$ \Gamma_{ij} = [F_{CT} \theta_{ij}^{T} + F_{CB} \theta_{ij}^{B} + F_{LS} \theta_{ij}^{LS}] \exp(-t_{ij}/\tau_{ij}) \quad (11) $$

$F_b$: Break-up coefficient in LS96 related to the coefficient $c_4$ in Eq. (27) of Luo and Svendsen (1996) by $c_4 = 0.923 F_b$.

$C_{Co}$: coalescence factor $\eta$: coalescence efficiency $\dagger$ Coefficients are modified in the IATE formulation due to integration over the diameter distribution function.

<table><caption>**Table 10**
Heat transfer models.</caption><thead><tr><th colspan="13">THERMAL MODELS</th></tr><tr><th>Model</th><th>ANSYS</th><th>CEA</th><th>FRA</th><th>HZDR</th><th>JSI</th><th>KTH-WSE</th><th>LEI</th><th>MIT</th><th>PSI-CTU</th><th>SJTU</th><th>UJV</th><th>VTT</th></tr></thead><tbody><tr><td>q'<sub>i</sub><br>T'<sup>sat</sup> &lt; T<sup>sat</sup></td><td>RM52</td><td>Ma00</td><td>RM52</td><td>RM52</td><td>RM52</td><td>RM52</td><td>RM52</td><td>RM52</td><td>RM52</td><td>RM52</td><td>RM52</td><td>RM52</td></tr><tr><td>q'<sub>i</sub><br>T'<sup>sat</sup> &gt; T<sup>sat</sup></td><td>RM52</td><td>Be83</td><td>RM52</td><td>KR11</td><td>RM52</td><td>RM52</td><td>RM52</td><td>RM52</td><td>RM52</td><td>RM52</td><td>RM52</td><td>RM52</td></tr><tr><td>q'<sub>v</sub></td><td>CTS</td><td>CTS</td><td>Nu<sub>v</sub></td><td>Nu<sub>v</sub></td><td>Nu<sub>v</sub></td><td>Nu<sub>v</sub></td><td>Nu<sub>v</sub></td><td>Nu<sub>v</sub></td><td>Nu<sub>v</sub></td><td>Nu<sub>v</sub></td><td>Nu<sub>v</sub></td><td>Nu<sub>v</sub></td></tr><tr><td>τ<sub>v</sub></td><td>0.05</td><td>0.05</td><td>n.a.</td><td>n.a.</td><td>n.a.</td><td>n.a.</td><td>n.a.</td><td>n.a.</td><td>n.a.</td><td>n.a.</td><td>n.a.</td><td>n.a.</td></tr><tr><td>Nu<sub>v</sub></td><td>n.a.</td><td>n.a.</td><td>26</td><td>∞</td><td>500</td><td>100</td><td>10</td><td>26</td><td>∞</td><td>1000</td><td>∞</td><td>10</td></tr><tr><td>q'<sub>i</sub></td><td>SGDH</td><td>SGDH</td><td>SGDH</td><td>SGDH</td><td>SGDH</td><td>SGDH</td><td>SGDH</td><td>SGDH</td><td>SGDH</td><td>SGDH</td><td>SGDH</td><td>SGDH</td></tr><tr><td>P<sub>r</sub><sub>i</sub></td><td>0.85</td><td>0.9</td><td>0.9</td><td>0.9</td><td>0.9</td><td>0.85</td><td>0.85</td><td>0.9</td><td>0.85</td><td>0.9</td><td>0.9</td><td>0.85</td></tr><tr><td colspan="13">Be83 (Berne, 1983)</td></tr><tr><td colspan="13">RM52 (Ranz and Marshall, 1952)</td></tr><tr><td colspan="13">Hu67 (Hughmark, 1967)</td></tr><tr><td colspan="13">Ma00 (Manon, 2000)</td></tr><tr><td colspan="13">KR11 (Krepper and Rzehak, 2011)</td></tr></tbody></table>

n.a.: Not applicable; CTS: Constant time scale ($\tau_v$); $\infty$: Zero resistance; $\tau_v$: Time scale for vapour thermal equilibrium; $Nu_v$: Nusselt number for vapour to interface heat-transfer; SGDH: Standard Gradient Diffusion Hypothesis; $P_{r_i}$: Turbulent Prandtl number.

<table><caption>**Table 11**
Wall heat flux partitioning models.</caption><thead><tr><th colspan="13">WALL HEAT FLUX PARTITIONING</th></tr><tr><th>Model</th><th>ANSYS</th><th>CEA</th><th>FRA</th><th>HZDR</th><th>JSI</th><th>KTH-WSE</th><th>LEI</th><th>MIT</th><th>PSI-CTU</th><th>SJTU</th><th>UJV</th><th>VTT</th></tr></thead><tbody><tr><td>Parti-tioning</td><td></td><td>Fa23</td><td></td><td>Di18</td><td>Se21b</td><td>La06</td><td>La06</td><td>MITB</td><td></td><td>La06</td><td>RPI</td><td>La06</td></tr><tr><td>q<sub>e</sub></td><td>ext. RPI</td><td>mod. KPI</td><td>KPI</td><td>Di18</td><td>Se21b</td><td>RPI</td><td>RPI</td><td>MITB</td><td>ext. RPI</td><td>Wa23</td><td>RPI</td><td>ext. RPI</td></tr><tr><td>q<sub>h</sub></td><td>DVK</td><td>DVK</td><td>DVK</td><td>Di18</td><td>Se21b</td><td>DVK</td><td>DVK</td><td>MITB</td><td>DVK</td><td>DVK</td><td>DVK</td><td>DVK</td></tr><tr><td>A</td><td></td><td>DVK</td><td>DVK</td><td>Di18</td><td>Se21b</td><td>DVK</td><td>Pe19</td><td>MITB</td><td>DVK</td><td>DVK</td><td>DVK</td><td>DVK</td></tr><tr><td>d<sub>d</sub></td><td>mod. TK70 &amp;Kr14</td><td>Un76</td><td>Un76</td><td>Di18</td><td>Se21a</td><td>TK70</td><td>TK70</td><td>MITB</td><td>mod. TK70</td><td>FBM</td><td>KI83</td><td>TK70</td></tr><tr><td>η<sub>v</sub></td><td>Co60</td><td>Co60</td><td>Co60</td><td>Co60</td><td>Co60</td><td>Co60</td><td>Co60</td><td>MITB</td><td>Co60</td><td>Ye14</td><td>Co60</td><td>KI95</td></tr><tr><td>N<sup>e</sup><sub>b</sub></td><td>mod. Co74</td><td>KPI</td><td>KPI</td><td>Di18</td><td>Se21b</td><td>LC77</td><td>LC77</td><td>MITB</td><td>Li18</td><td>LC77</td><td>KI83</td><td>LC77</td></tr><tr><td colspan="13">Co60 (Cole, 1960)</td></tr><tr><td colspan="13">DVK (Del Valle and Kenning, 1985)</td></tr><tr><td colspan="13">KI83 (Kocamustafaogullari and Ishii, 1983)</td></tr><tr><td colspan="13">Kr14 (Krepper and Rzehak, 2014)</td></tr><tr><td colspan="13">LC77 (Lemmert and Chawla, 1977)</td></tr><tr><td colspan="13">RPI (Kurul and Podowski, 1991)</td></tr><tr><td colspan="13">TK70 (Tolubinsky and Kostanchuk, 1970)</td></tr><tr><td colspan="13">Wa23 (Wang et al., 2023)</td></tr><tr><td colspan="13">Co74 (Cole, 1974)</td></tr><tr><td colspan="13">Fa23 (Favre, 2023)</td></tr><tr><td colspan="13">KI95 (Kocamustafaogullari and Ishii, 1995)</td></tr><tr><td colspan="13">La06 (Lavieville et al., 2006)</td></tr><tr><td colspan="13">MITB (Lahey et al., 2021)</td></tr><tr><td colspan="13">Se21a (Setoodeh et al., 2021b)</td></tr><tr><td colspan="13">Un76 (Unal, 1976)</td></tr><tr><td colspan="13">Di18 (Ding et al., 2018)</td></tr><tr><td colspan="13">HI03 (Hibiki and Ishii, 2003)</td></tr><tr><td colspan="13">KPI (Kurul and Podowski, 1990)</td></tr><tr><td colspan="13">Li18 (Li et al., 2018)</td></tr><tr><td colspan="13">Pe19 (Peltola et al., 2019)</td></tr><tr><td colspan="13">Se21b (Setoodeh et al., 2021a)</td></tr><tr><td colspan="13">Ye14 (Yeoh et al., 2014)</td></tr></tbody></table>

mod.: Modified; $d_d$: Bubble departure diameter; $\eta_d$: Bubble departure frequency; $q_e$: Evaporation heat flux; $q_h$: Quenching heat flux; $N_b^e$: Nucleation site density; $A$: Area of influence; FBM: Force Balance Model;

Regarding heat transfer (Table 10), interfacial model from Ranz and Marshall (1952) is used by all participant in subcooled conditions, mostly due to the fact that calibration efforts where oriented towards

void fraction predictions rather than liquid temperature. Only CEA and HZDR use a different model for the over-heated region. The influence of the vapour Nusselt number is negligible provided that it is sufficient
to ensure a quick return to saturation. Turbulent heat transfer was not specifically studied by the participants that all kept a classical SGDH with a turbulent Prandtl number in 0.85 to 0.9.

WHP has been the focus of longstanding research and it concentrated many calibration efforts also in this benchmark. Table 11 shows a small part of the complexity of this model. Some general trends can be observed, like (Del Valle and Kenning, 1985) for quenching heat flux and bubble area of influence, Tolubinsky and Kostanchuk (1970) for departure diameter, Cole (1960) for bubble departure frequency or Lemmert and Chawla (1977) to determine the nucleation site density. However, the details in Appendix A clearly demonstrates that this seeming unity is in fact broken by the different values chosen for adjustment parameters. More than all the previous models, the calibration of the WHP model is considered as a difficult topic, and its performance depends on the user's expertise. In particular, HZDR and MIT have a long research history in this field and they have both developed a very specific and global set of closure for the prediction of the partitioning.

WHP is also strongly linked to IAC predictions as the departure diameter and the vapour production source term provide boundary conditions to diameter evolution.

This brief analysis shows that the modelling trends differ between the participants, despite sharing certain references. This will enable a cross-analysis of the capabilities of the different models' sets. More detail specific to each participant's closures is given in Appendix A.

## 5. Results and analysis

The goal of the benchmark is to predict several quantities at the end of the heating section. Even though local measurements were performed at a unique axial location (z = 3.485 m), it is possible to consider that moving the sensor along the vertical axis is equivalent to a change in the mean equilibrium quality $x_{eq}$ resulting from a change of the inlet temperature, all the other boundary conditions being kept constant (up to the small variations of the loop control system). To be valid this assumption needs the resulting pressure drop variation to be weak and to have no effect on the flow pattern. Therefore, test series with the same triplet (G, P, W) are analysed together and can be considered as the flow evolution for the case with the highest inlet temperature.

A comparison is also performed further downstream in the adiabatic region to assess the variability between the participants, even though no experimental reference is available there.

Code comparisons are presented on the void fraction, Sauter mean diameter, vapour velocity and liquid temperature measurements when available. Additionally, other profiles such as liquid velocity, liquid turbulence, wall heat flux partitioning and wall temperature are also analysed in code comparison as their experimental measurements were impossible. Detailed illustrations are presented in Electronic Supplements A and B. The DEBORA cases have been organised in four series based on the operating pressure and coolant mass flux: G2P14W16, G2P26W16, G3P26W23 and G5P14. With the series defined, the comparison is more straightforward, and the influence of inlet subcooling for each series can be emphasised.

The current section is organised as follows. First, the effect of modelling choices is discussed (Section 5.1) with a focus on the consequences of model's variations on flow's predictions. Then, Section 5.2 presents the effect of flow parameters. The effect of progressive increase of inlet temperature is discussed for each of the four series. The effect of physical properties is analysed on a specific condition in Section 5.3. Finally, the specific role of two velocity groups for the vapour phase in iMUSIG approaches is presented in Section 5.4.

Additionally, Electronic Supplement A presents in greater detail individual analyses associated to each participant, followed by detailed figures on flow evolution. Then, Electronic Supplement B presents the

comparison of all the sets of closures considered for each physical condition, for several variables, both in radial and streamwise directions, on the cross-section average void fraction, and at the wall (temperature and heat partitioning). Numerical data are reported in Bois et al. (2024).

### 5.1. Effect of models

First, at low quality, liquid streamwise velocity and turbulent kinetic energy are similar in trend and amplitude for all participants. But very rapidly, differences in void fraction predictions distribute the buoyancy effect differently and the profiles gets different shapes. Specific production terms for BIT in CEA and HZDR models both induce higher turbulent kinetic energy than Sato and Sekoguchi (1975). Turbulence modelling was not found very influential and the relevance of addition of BIT is not certain.

Generally, the lift force inversion is one of the main focus of calibration efforts. It was often assessed mostly on G2P14W16 series. Both the intensity and the conditions of inversion of the lift coefficient were addressed and quite different closures were selected. The diverse behaviours are observable on radial void fraction profiles with higher void fraction in the core flow for models with the largest lift coefficient. Near-wall damping or additional wall force are also found important. In competition with this force detaching void from the wall, the differences of order of magnitude in the turbulent dispersion forces conducted to extremely scattered slopes for $a(r^+)$). The strongest turbulent dispersion calibrated by VTT lead to the flattest profiles, whereas the steepest ones are achieved by UJV and HZDR, for a combination of Tomiyama et al. (2002) model for lift without damping and a small turbulent dispersion from Burns et al. (2004) model.

Diameter prediction was definitely the most difficult variable to assess, with a lot of parameters to adjust, especially in the MUSIG or iMUSIG approaches. Different ranges where taken for the two pressures conditions, and the selection of the pivotal diameter separating the two velocity groups in iMUSIG approaches is very important. Near-wall values are mostly dictated by the closure for departure diameter. A quite good agreement with the experiment is observed, with a moderate scattering due to similar models for all participants. In contrast, the diameter radial evolutions $d_{32}(r^+)$ are very scattered both because (i) coalescence and break-up source terms reach large but opposite values and (ii) the models to predict IAC are strongly affected by variables difficult to determine for which a great dispersion among participant is observed (sub-cooling, turbulent kinetic energy, dissipation, ...). Large disparities in coalescence and break-up coefficients lead to significant spreading in diameter predictions. More in-depth analyses could be performed by investigating the influence on the shape of the bubble size distribution and not only its mean value.

Regarding heat transfer, the participants chose not to focus on the in-flow thermal aspects due to their higher expertise in near-wall effects (WHP models). Turbulent heat transfer and recondensation were not analysed deeply as most participants focused their previous research on the WHP modelling assumed dominant. Generally, models tend to over-predict the temperature, especially in the near-wall region where over-heated liquid is frequently predicted whereas it is not observed in the experimental measurements. Very similar slopes on $T^l(r^+)$ are observed as expected because participants use quasi-identical correlations for interfacial heat transfer (Ranz and Marshall, 1952) and for turbulent Prandtl number.

Lastly, the complexity and diversity of WHP models is striking on streamwise profiles of evaporation heat flux. It also affects the streamwise evolution of the mean void fraction $\langle\alpha\rangle(z^+)$ but the radial profile determined by the forces equilibrium also plays an important role in it. The difficulty of WHP models and the large number of degrees of freedom are clearly observed.
<figure><img src="image_8.png" /><figcaption>Fig. 5. Streamwise evolution of averaged void-fraction $\langle a \rangle$ and wall evaporation heat flux $q_e$ for case 07_30G2P14W16Te40, series G2P14W16.</figcaption></figure>

**IATE, MUSIG or iMUSIG.** Amongst the different models, one strong difference concerns the modelling of IAC that can be split into 3 categories: IATE, MUSIG or iMUSIG. In Electronic Supplement B, participants with a single IATE are represented by continuous lines whereas MUSIG and iMUSIG users are shown with dash-dotted and dashed lines respectively.

The results do not demonstrate a clear superiority of the iMUSIG model in predicting void fraction profiles even though it contains more complexity and embeds the potential for having opposite radial velocities for the two bubble groups. More detailed experimental post-processing is needed to inform on the actual local distribution of bubble sizes to know if two peaks are present in the distribution. Specific analyses of iMUSIG parameters is presented in Section 5.4.

**WHP and evaporation heat fluxes.** The WHP model and its contribution to the evaporation heat flux is an essential component of boiling models. It brings large complexity and strong differences between the approaches. As an illustration, Fig. 5(b) for case 07 represents the streamwise evolution of the evaporation heat flux that contributes to vapour generation at the wall. Differences by an order of magnitude are observed and trends are very different with some models predicting a decrease of evaporation heat flux at some point while others show a monotonous increase until the end. Even though case 07 has a very high outlet quality and measurements show that it is strongly core-peaked at the end of the heating section (see void fraction on Fig. 6(d)), the streamwise evolution is representative of all the flow evolution from the beginning of boiling. These large differences are expected to strongly affect the averaged void fraction profiles as illustrated on Fig. 5(a).

More generally on all cases, the averaged void fraction profiles show a very large spreading among the participants, including in particular very large differences in the prediction of the Onset of Significant Void (OSV). These differences are directly connected to differences in the evaporative fraction in the heat partitioning models.

**Models strength and weaknesses.** The models have to consider a large number of physical mechanisms and conduct to reasonable predictions. The most important one are the lift force and its near-wall modification, the turbulent dispersion, the parameters of PBM for IAC prediction, the departure diameter and WHP models.

Parameters for PBM are rather unknown and hard to set for boiling flows. For example, aggregation and turbulent break-up can have a large influence on bubble diameter but almost none on the vapour distribution (Gajšek et al., 2023). The potential gain with the iMUSIG expected from the supplementary degrees of freedom on velocity is not obvious. Unfortunately, experimental data on the bubble size distribution is lacking to determine whether a bimodal distribution exists and is the cause of complex lateral motion of vapour. Besides, lift force modification is not trivial and could not consistently lead to prediction

of experimental observations of centre peaks. Lift and dispersion forces dominate the radial motion of vapour. They are competing with the temperature distribution to determine the void fraction profile.

Calibration was mostly performed on G2P14W16 series and was sometimes found difficult to extend to different pressure and flowrate conditions. Detailed and complete experimental data is extremely difficult to obtain, especially in high-pressure conditions. This lack of information on mechanisms for IAC or temperature evolutions often conducts to the use of theoretical models from standard literature, without justifying their application in the condition studied; this default choice does not necessarily mean it is a good choice, nor that the actual mechanism has as weak an influence as it may seem.

WHP models need a great expertise and predict very different values.

## 5.2. Effect of flow parameters and physical conditions

Global comparisons of the participants on all test cases are presented in Electronic Supplement B. Significant spreading of the predictions are observed in particular on void fraction, diameter, liquid temperature or evaporation wall heat flux.

From an averaged point-of-view, the predictions of the streamwise flow evolution are very different too with some models sometimes generating a significant amount of void fraction since the beginning of the heated length. In general, it is very difficult to predict void fraction maximum and migration over the whole range of conditions explored. Few models are roughly capable of capturing the transition in G2P14W16 series: mostly MIT and PSI-CTU reach a correct shape of void fraction profile for case 07 (see Fig. 6(d)). For the case 08, it would be necessary to lower the lift force coefficient or turn it off completely for the PSI-CTU model, because the average void fraction exceeds 0.6 for which the concept of bubble lift force cannot be applied, thus the vapour void fraction is here strongly overpredicted. Yet, at higher pressure, for the highest quality case of G2P26W16 (see case B21 on Fig. 8(a)), these two models predict an intermediate peak that is not clearly present in the experiment. Besides, most models strongly overpredict void fraction in the core in this condition.

An other example of difficulties is given by case O25 from G3P26W 23 (see Fig. 7(d)) where all models over-predict void fraction and fail to detect the intermediate peak that starts to move away from the wall.

The comparison of cases 12 and 13 from G5P14W29 (see Figs. 9(a) and 9(b)) shows that the maximum void fraction suddenly moves towards the core. Several models where in close agreement with condition 12 but they all fail to capture the transition occurring in case 13.

Some analyses were performed to compare the participant's results and analyse the radial temperature profiles in view of the previously
<figure><img src="image_9.png" /><figcaption>Fig. 6. Radial void fraction profiles at $z = 3.485$m for some cases from G2P14W16.</figcaption></figure>

published experimental profiles. These new cases were compared to open cases where other experimental measurements were available (void fraction, diameter and vapour velocity). The results were compared according to the experimental resemblance and analyses were attempted to understand inter-connections between variables.

**G2P14W16.** In this lower pressure and low flow-rate conditions, the exit qualities range in $x_{eq}^{out} \in [-0.06; 0.11]$. Temperature measurements in T28 on Fig. 12(a) are performed in conditions close to case 01 (see Fig. 6(a) and Electronic Supplement Fig. C2) where void fraction, diameter and vapour velocity were measured. For this subcooled condition, a large dispersion of predictions is already observed on all variables. In particular, MIT and VTT models lead to the largest over-predictions of void fraction and vapour velocity whereas a reasonable agreement is reached on temperature.

For a slightly higher quality, cases T29, 02 and 03 are rather similar. The void fraction maximum starts to detach from the wall. It is detected by MIT, SJTU and VTT with various accuracies; and a flattening trend is also observed for LEI and PSI-CTU (see Figs. 6(b) and 6(c)). Temperature profiles are best predicted by ANSYS and JSI for T29 (see Fig. 12(b)).

Moving to saturated cases, the peak transition to the core is fully complete in case 07. Only MIT and PSI-CTU reach a correct shape of void fraction profile for this case (see Fig. 6(d)); they are the only models with negative global lift coefficients over the whole cross-section. In the case of iMUSIG, the coefficient is negative for the larger groups (i.e., the second velocity), but the contribution of the smaller bubbles remains too significant (and with a positive lift coefficient).

**G2P26W16.** For this higher pressure series, the exit qualities range in $x_{eq}^{out} \in [-0.53; 0.08]$. At $P = 26$ bar, the vapour phase remains finely dispersed and the maximum void fraction remains at the wall for all

conditions (Hösler, 1968; François et al., 2011). These conditions are better predicted by most models with still a large spreading in void fraction levels (see for instance case B18 on Fig. 7(a)) with predictions $\alpha \in [0; 0.2]$ at $r^+ = 0.8$. These differences are also visible on averaged void fraction increase on Fig. 11 exhibiting very different behaviour for each participant. For the highest quality, the experimental vapour velocity profile is flat for at least $r^+ < 0.8$ but this tendency is not captured by any model; they all predict a profile below experimental value and very close to liquid velocity which remains in a rather classical single-phase shape. This weakness could be attributed either to relative velocity modelling and drag force closure or to turbulence models in the liquid that are not accurate on the bubble induced effect.

Condition B21 is close to T33 with a similar exit quality. Most models over-predict liquid temperature (Fig. 12(c)), void fraction and bubble mean diameter (Fig. 8) but underestimate the vapour velocity.

**G3P26W23.** At a slightly higher flow-rate, with exit qualities ranging in $x_{eq}^{out} \in [-0.05; 0.10]$, we see once again that the models have difficulties to predict the rate at which the mean void fraction increase with increasing inlet temperature. The maximum void fraction starts to move away from the wall for case 10 (Fig. 7(c)) but numerical predictions are far from measurements; in particular, at the core, void fraction may be more than twice over-estimated. Profiles from MIT, SJTU, PSI-CTU and VTT have initiated a transition in void fraction regime. Again, participants' predictions of diameters are extremely dispersed.

**G5P14.** Tests in this series inform on the effect of the wall flux intensity on near-wall profiles. Four conditions (cases 11 to 14) where performed with a flux of $q_w = 135 \text{ kW m}^{-2}$ while cases B26 and B27 were realised at higher fluxes of 197.12 kW m$^2$ and 152.83 kW m$^2$ respectively.

Regarding outlet qualities, cases 12 and B26 are similar with $x_{eq}^{out} \approx 0$ and cases 13 and B27 with $x_{eq}^{out} \approx 0.05$ can also be compared.
<figure><img src="image_10.png" /><figcaption>Fig. 7. Radial void fraction profiles at $z = 3.485$ m for cases B18, B21, 10 and O25.</figcaption></figure>

Comparing case 12 to B26 is particularly interesting as it leads to very different conclusions depending on the participants. Whereas for most participants void fraction predictions are identical for cases 12 and B26, Fig. 9 shows that significant differences appear in the near-wall region for HZDR, MIT, PSI-CTU and SJTU models. Other variables are weakly affected (except the diameter and liquid temperature for HZDR, see Electronic Supplement Figs. C122 and C134). The most sensitive model is from HZDR with huge variations on the void fraction, largely overestimated (Fig. B63 in Electronic Supplement). MIT profiles show a strong increase of void fraction by 0.2 localised near the wall (Fig. B127). PSI-CTU has the best trends, very similar to experiment on these two conditions (Fig. B144). SJTU is slightly different and shows an intermediate peak (Fig. B160).

On cases 13 and B27, most participants have similar predictions. Again, MIT results show a strong increase of void fraction by 0.2 near the wall which is a similar trend as in the experiment. MIT results show the correct difference in void fraction between those four cases. HZDR results show a much smaller difference between those conditions than in the comparison of cases 12 and B26. Again, PSI-CTU has very good trends similar to experiment but do not detect the sharp near-wall increase in case B27. SJTU shows impact in an intermediate region $0.75 < r^+ < 0.95$ but the impact of flow condition is very limited closer to the wall where it is expected to be larger (as confirmed by the experimental results). SJTU also has slight differences on liquid temperature. These differences in SJTU predictions come from the damping function introduced in the lift model.

Cases 14 and B27 present two local maxima of void fraction. In those conditions, the bubbles are slowly migrating away from the wall but the second very sharp peak is appearing at the wall. This sharp increase is the signature of the proximity to CHF. The experiment reveals that it appears brutally between cases 13 and 14 and is also

present in B27 but it is not triggered similarly at higher wall flux in B26 where the void fraction remain simply wall peaked. These variations are extremely interesting. Unfortunately, no model is accurate enough to detect this kind of flow behaviour transition reliably; nonetheless, detecting these transitions is an essential requirement to the development of predictive capabilities for CHF detection using MCFD. This challenge still lies ahead.

**Temperature predictions.** Two kind of temperature predictions are generally observed, either a rather parabolic (almost laminar) profile or much sharper turbulent profile with a thin boundary layer. Besides, for many cases and participants, the saturation temperature is exceeded in the near-wall region; then, condensation is absent and the full heat flux goes into the liquid and contributes to large overheats (up to approximately 4K) and to a thick overheated layer (sometimes even in the whole cross-section). In contrast with these predictions, the few experimental measurements available demonstrate that the liquid remains subcooled almost everywhere (see for instance the $T^l$ profile of case T33 on Fig. 12(c) for an outlet quality $x = 0.0806$), except possibly in an extremely thin boundary layer. This inconsistency tends to indicate shortcomings in the force balance model as both heat-partitioning and condensation models are deactivated at the end of this case.

### 5.3. Effect of physical properties

In G2P26W16, cases T30 and T31 are in single phase conditions and case T32 is expected to be very close to OSV. Conditions T31 and T32 are selected to assess the effect of physical properties used by each participant depending on their own EOS tables. While the differences between the property treatments are largely negligible in
<figure><img src="image_11.png" /><figcaption>Fig. 8. Radial profiles at z = 3.485 m. Series G2P26W16. B21_30G2P26W16Te70.</figcaption></figure>

two-phase conditions, some differences were observed on those single-phase conditions. Indeed, Figs. 10(a) and 10(b) reveal differences of up to ±10% on streamwise liquid velocity and more than ±5 °C on liquid temperature (on Figs. 10(e) and 10(f)) even though most participants have very similar results. These variations are explained by small discrepancies of physical properties on Fig. 13. The largest differences are observed on viscosity.

Additional simulations were performed with the same constant properties given in Table 12 for all participant. These simulations referred to as U31 and U32 are presented on Figs. 10(c), 10(d), 10(g), and 10(h). As expected, liquid velocity and temperature predictions are identical for $r^+ < 0.9$ and differences arise near the wall depending on the meshes, wall laws and turbulent models used. Only SJTU prediction remains quite different. Temperature predictions for case U32 (with $T_e = 49.7$ °C) are in very good agreement with the experiment on Fig. 10(h) but differences by 5 °C to 10 °C are observed at the wall. In contrast, Fig. 10(f) for the same condition $T_e = 49.7$ °C performed with user-defined variable properties (case T32) shows significant differences in temperature profiles that can thus be mostly attributed to properties predictions. Unfortunately, these differences impact the predictions for the whole series and for instance at $T_e \approx 70$ °C in cases B21 and T33 (see void faction on Fig. 8(a) and liquid temperature on Fig. 12(c)). The large difference between HZDR and LEI temperature predictions at the incipience of boiling ($T_e = 49.7$ °C, case T32) reduces to less than 1 °C for the highest quality condition at $T_e \approx 70$ °C but the strong effect of physical properties remain is converted into a difference of 5% in void fraction prediction on Fig. 8(a).

<caption>Table 12</caption>

Constant properties at $T = 70$ °C and $P = 26$ bar from (Lemmon et al., 2023) used for cases U31 and U32.

<table><thead><tr><th></th><th>ρ, kg/m³</th><th>μ, Pa s</th><th>λ, W/(m K)</th><th>C<sub>p</sub>, J/(kg K)</th></tr></thead><tbody><tr><td>Liquid</td><td>1130.5</td><td>1.1729 × 10<sup>-4</sup></td><td>0.051 803</td><td>4.21 × 10<sup>3</sup></td></tr><tr><td>Vapour</td><td>170.69</td><td>1.6326 × 10<sup>-5</sup></td><td>0.014 762</td><td>2.077 × 10<sup>3</sup></td></tr></tbody><tfoot><tr><td colspan="5">σ = 0.0018 N/m, L<sup>vap</sup> = 86.48 kJ/kg</td></tr></tfoot></table>

## 5.4. iMUSIG analysis

Bubbly flow in DEBORA experiment is poly-disperse. Break-up, coalescence and condensation of bubbles lead to a wide bubble size distribution. The bubbles migrate from the heated wall to the pipe centre due to forces acting in this lateral direction, which is different from the direction of main drag force. These lateral bubble forces as well as the drag and the virtual mass force depend on the bubble size. The hydrodynamic lift force on bubbles even changes its sign as the bubble diameter changes (Tomiyama, 1998). As a result, the lateral migration of the bubbles leads to a radial separation of small and large bubbles and to a further coalescence of large bubbles that move towards the pipe centre. The iMUSIG model can consider different velocity fields for bubbles of different size. Within the iMUSIG model the dispersed gaseous phase is divided into several inhomogeneous velocity groups and each of these groups is subdivided into bubble size groups. Bubble break-up and coalescence between all bubble size groups are taken into account by appropriate models. Details of the model and its validation for vertical pipe flows are presented in paper by Kreppe et al. (2008). The differences between homogeneous MUSIG and iMUSIG models for
<figure><img src="image_12.png" /><figcaption>Fig. 9. Radial void fraction profiles $\alpha$ at $z = 3.485$ m. Series G5P14.</figcaption></figure>

subcooled boiling flow are shown in Fig. 14. In order to distinguish the movement of large and small bubbles, iMUSIG model was employed in the simulations by HZDR, SJTU and UJV.

The details of the implementation of iMUSIG model by SJTU are as follows. The bubbles are firstly divided into two velocity groups with the critical bubble diameters, i.e., 1 mm for cases at 26 bar and 1.5 mm for 14 bar according to the Tomiyama (1998) lift force model. The bubbles are further divided into 25 size groups. The groups 0 to 14 at 14 bar (groups 0 to 9 at 26 bar) are located in velocity group 0 while the groups 15 to 24 at 14 bar (groups 10 to 24 at 26 bar) are located in velocity group 1.

The void fraction in cases #1 and #5 in Table 5 are shown in Fig. 15. The void fraction of case #1 is relatively small, thus the collision process of the released bubbles from heated wall is weak. Few large bubbles form in velocity group 1, while most of bubbles

belong to velocity group 0 and accumulate in the near wall region. Therefore, void fraction of velocity group 1 approaches zero in case #1. When much more bubbles are nucleated on the heated wall, many large bubbles form and migrate into the pipe centre, their motion being induced by the lift force. However, the total void fraction profile still deviates from the experimental data shown in Electronic Supplement Fig. B148. This may be caused by the accuracy of the interphase force models since the effects of void fraction on the interphase force are not considered in the available models. Therefore, more efforts should be devoted to the development of interphase forces models.

UJV used the iMUSIG model with 10 bubble size groups and 2 velocity groups. The bubble diameter ranges in two velocity groups were set from 0 mm to 1.4 mm and from 1.4 mm to 2 mm. The effect of iMUSIG model on case #8 (Table 5) with the void fraction core peak is shown in Fig. 16. The figure shows the total void fraction, void fraction
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Radial profiles of liquid velocity $\mathbf{u}^l$ and temperature $T^l$ at $z = 3.485\,\text{m}$ for cases T31 and T32 with participant's physical properties and cases U31 and U32 where uniformised constant physical properties were imposed.</figcaption></figure>

of small bubbles in the velocity group 0 and the void fraction of large bubbles in the velocity group 1. Small bubbles remain close to the wall and the large bubbles are drawn to the core of the flow.

# 6. Sensitivity analyses

## 6.1. User, code and mesh effects

In view of the complexity of the physical models and associated numerical schemes, a working group from several benchmark participants

(JSI, LEI, KTH-WSE and UJV) was established to investigate potential effects of users, codes and computational meshes on the simulation results. The CFD codes used in the working group were as given in Table 1: ANSYS Fluent version 2021 R2 (JSI), OpenFOAM version 8 (LEI) and 9 (KTH-WSE), and ANSYS CFX version 18.0 (UJV). A baseline model was first selected by the working group comprised of physical models available in all codes. The baseline model is summarised in Table 13. For simplicity, a constant bubble diameter of $d_{32} = 0.45\,\text{mm}$ was selected. All other models and parameters found in Table 8 to 11 are either disabled or not applicable. The “zero-resistance” interfacial
<table><caption>**Table 13**
Baseline model. See Tables 8 to 11 for the references to the abbreviations.</caption><thead><tr><th></th><th>Model</th><th>Baseline</th></tr></thead><tbody><tr><td colspan="3">FLOW DYNAMICS</td></tr><tr><td rowspan="3">Turbulence</td><td>Continuous</td><td>k &minus; &omega; SST</td></tr><tr><td>BIT</td><td>SS75</td></tr><tr><td>Disperse</td><td>&mdash;</td></tr><tr><td rowspan="6">Interfacial forces</td><td>Drag</td><td>IZ79</td></tr><tr><td>C<sub>AM</sub></td><td>&mdash;</td></tr><tr><td>Lift</td><td>To02</td></tr><tr><td>Wall-damping</td><td>&mdash;</td></tr><tr><td>Wall force</td><td>&mdash;</td></tr><tr><td>TD</td><td>Be97</td></tr><tr><td colspan="3">INTERFACIAL AREA MODELLING</td></tr><tr><td colspan="3">Constant diameter</td></tr><tr><td colspan="3">THERMAL MODELS</td></tr><tr><td rowspan="6">q'<sub>i</sub></td><td>T<sup>l</sup> &lt; T<sup>sat</sup></td><td>RM52</td></tr><tr><td>T<sup>l</sup> &gt; T<sup>sat</sup></td><td>RM52</td></tr><tr><td>q'<sub>v</sub></td><td>Nu<sub>v</sub> = &infin;</td></tr><tr><td>q'<sub>f</sub></td><td>SGDH</td></tr><tr><td>Pr<sub>i</sub></td><td>0.9</td></tr><tr><td colspan="3">WALL HEAT FLUX PARTITIONING</td></tr><tr><td></td><td>d<sub>d</sub></td><td>TK70</td></tr><tr><td></td><td>q<sub>e</sub></td><td>RPI</td></tr><tr><td></td><td>Partitioning</td><td>RPI</td></tr><tr><td></td><td>&eta;<sub>d</sub></td><td>Co60</td></tr><tr><td></td><td>N<sup>e</sup><sub>b</sub></td><td>LC77</td></tr><tr><td></td><td>Area influence</td><td>DVK</td></tr><tr><td></td><td>q<sub>q</sub></td><td>DVK</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Streamwise evolution of area-averaged void-fraction ($\langle a \rangle$). Series G2P26W16. Case B21_30G2P26W16Te70.</figcaption></figure>

vapour heat transfer model corresponds to a sufficiently large $Nu_v$, which value does not affect the results.

Through several iterations within the working group, the empirical constants used in the models were carefully compared among participants and made equal in case of discrepancy. This includes $d_{Ref} = 0.6$ mm, $d_{Min} = 1$ µm, $d_{Max} = 1.4$ mm for TK70, $m = 210$, $n = 1.805$ for LC77 and $K = 4.8$, a waiting time ratio of 0.8, the liquid temperature $T^l$ at $y^+ = 250$ for DVK. However, the fluid properties considered by each code user came from different sources (NIST for JSI and UJV, CoolPack for LEI and CoolProp for KTH-WSE). It has been considered that only minor differences could be observed from this difference. In addition, some potentially relevant model parameters were not adjustable, for instance the maximum bubble area of influence in the wall heat flux partitioning model (set to 1 for all codes, except to 0.8 for CFX).

The base geometry and computation mesh was the same for all codes and consistent with the description provided in Section 3.1 (3.5 m heated with 1 m upstream and 0.5 m downstream adiabatic), except for KTH-WSE which did not consider the 1 m adiabatic upstream length

(however, inlet effects are considered negligible). The base computation mesh was set to 500 uniformly distributed axial cells and 40 uniformly distributed radial cells. The selected test case corresponds to case #13 29G5P14W29Te43.4 of Table 5. This is a low pressure (14.6 bar) high mass flux test case (5000 kg m⁻² s⁻¹) with slightly saturated output ($x_{out} = 0.0398$) but still near-wall peak void fraction distribution.

The results of the user and code effects are presented in Electronic Supplement C.1. When comparing the radial profiles of selected parameters at the end of heated length ($z = 3.485$ m on Electronic Supplement Fig. D2), non-negligible differences can be observed. While the profiles are generally similar, a clear difference in void fraction, phase axial velocities and turbulent kinetic energy can be observed due to code and user effects, larger than expected. The streamwise profiles at $r^+ = 0.9375$ (Fig. D3) and cross-section averaged (Fig. D4) confirm these observations. Note that one of the codes show significant instabilities for the selected reference models while others do not. Finally, the wall profiles (Fig. D5) confirm the overall non-negligible differences between submitted simulation results. Downstream the end of heated length, both the liquid and wall temperatures are predicted by all codes to be equal to the saturation temperature, as expected.

The results of the mesh effects are presented in Electronic Supplement C.2. For each considered code, various mesh configurations were tested. Overall, the base 500 × 40 computation mesh is sufficient to obtained converged results. Interestingly, it seems that ANSYS CFX (used by UJV) is more sensitive to mesh size while ANSYS Fluent (used by JSI) is less sensitive. Overall, however, the observed mesh effect is significantly lower than the observed user and code effect.

From the presented analysis, it becomes evident that the large number of required models in a two-phase flow convective boiling simulation, along with the large number of potential variations within the same model, can lead to significant differences that are challenging to fully track. While offering full flexibility, OpenFOAM was observed to be more error-prone without very careful review of the model input parameters. Other codes may suffer from lack of flexibility, making it challenging to establish a common baseline or to modify specific model parameters. In all cases, post-processing errors were also occurring mainly due to the lack of clarity regarding which parameters shall
<figure><img src="image_15.png" /><figcaption>Fig. 12. Radial temperature profiles at $z = 3.485$ m for cases T28 and T29 (similar to cases 01, 02 and 03) from G2P14W16 and T33 (similar to B21) from G2P26W16.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13. Physical properties at $P = 26.2$ bar. The terminology Ansys RGP gathers ANSYS, HZDR and PSI-CTU. Ansys Real Gas Property (RGP) tables contain liquid, vapour and saturation properties which are generated from NIST REFPROP Library (Lemmon et al., 2023).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 14. The framework of the MUSIG (a) and iMUSIG (b) model.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 15. Void fraction of velocity groups 0 and 1 at z = 3.485 m in cases #1 and #5 (Table 5) calculated by SJTU.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. Void fraction of velocity groups 0 and 1, and total void fraction at z = 3.485 m in case #8 (Table 5) calculated by UJV and comparison with total void fraction in experiment.</figcaption></figure>

be presented. Overall, this exercise should be consistently performed by teams aiming to adopt multiphase CFD, starting from a simpler adiabatic condition as presented for example by Sugrue et al. (2017), Marfaing et al. (2018), to reach a subcooled boiling inlet and finally a saturated boiling inlet leveraging the present work.

## 6.2. Sensitivity analysis of departure bubble parameters with the HZDR (CFX) and KTH-WSE (RPI in OpenFOAM) models

In the following section, we aim to present a sensitivity study of the parameters related to departure bubble dynamics on the heating wall in flow boiling simulations. This study is organised using two codes

with their respective wall boiling models. The HZDR CFX code has its own bubble dynamics model based on the force balance concept, while the KTH-WSE OpenFOAM code employs the conventional RPI wall heat partitioning model. The differences between these models have already been introduced in A.4 and A.6. The involved parameters are departure diameter $d_d$ (m11, m12), departure frequency $\eta_d$ (m13, m14), nucleate site density $N_b^e$ (m15, m16) and waiting time $t_w$ (m17, m18) (see Tab. D1 in Electronic Supplement). The total amount of departed bubble mass flux, corresponding momentum, and evaporating heat flux are calculated using these three parameters as follows:

$$ \dot{m}_v, q_e \propto d_d^3 \eta_d N_b^e \quad (12) $$

The radial distribution of the bubble Sauter diameter and void fraction distribution at the section of z = 3.485 m is first employed for comparison, followed by streamwise wall profiles at $r^+ = 1$ (see Electronic Supplement Figs. D2 and D5 respectively). Here again, case #13 29G5P14W29Te43.4 of Table 5 is considered. As stated in the results sections, both codes have good agreement with the experimental results using the default bubble dynamic parameters. Later, each parameter will be evaluated for sensitivity analysis by taking it as half and twice the original value.

**Departure diameter.** The Sauter diameter near the wall is 0.1 mm, 0.142 mm and 0.143 mm for departure diameters of $1/2 \times$, $1 \times$ and $2 \times$ of the original value in the HZDR model, while it is 0.525 mm, 0.615 mm and 0.7 mm in the KTH-WSE model. In the HZDR model, in the radial direction pointed to the centre, the Sauter diameter continues to increase when the departure diameter is 1/2 while it first increases, later decreases and finally increases again for the two other values of departure diameter. In the KTH-WSE model, it increases and decreases
for all. In the HZDR model, the wall void fraction of the 1/2× departure diameter is the largest while the one with the original departure is the smallest and the 2× is intermediate. The 1/2× case makes the local bubble smaller than the critical value of the lift force, and a large number of bubbles are pushed towards the wall. With an increase in the departure bubble to 2×, according to the stated equation (12), the local void fraction will increase as well. The void fraction peak appears on the wall in cases of 1/2× and 2× while showing a distance near the wall in the 1× case. The lift force effect on the wall void fraction is not shown in the KTH-WSE model which shows only an increase in the void fraction with an increase in the departure diameter. However, all void fraction peaks are at a distance near the wall.

For the streamwise wall evolution (see Electronic Supplement Fig. E3), in the HZDR model, a 2× departure bubble leads to the maximum Sauter diameter at $r^+ = 1$ appearing upstream, but the Sauter diameter decreases downstream later. Meanwhile, the 1/2× case has a peak of Sauter diameter much further downstream and the value is also the smallest among the three cases. The wall temperature along the pipe first increases and then decreases. The wall temperature peak appears earlier for the larger departure diameter (2×). However, compared to the other two cases (1/2×, 1×), it has the smallest value. The difference in the wall temperature peak between the 1/2× and 2× cases is 3 K. Similar to the wall temperature, the evaporative heat flux first increases and later decreases for all three cases. But with a larger departure diameter, the peak comes earlier and the value is larger. In the KTH-WSE model, the Sauter diameter always increases with an increase in the departure diameter along the pipe. The difference in Sauter diameter at the same position is nearly in the same proportion of 1/2×, 1×, and 2×. The wall temperature always increases and then decreases for all cases. The difference in peak temperature between cases 1/2× and 2× is 7 K. The evaporative heat flux in the KTH-WSE model consistently increases along the pipe. The peak value increases with an increase in departure diameter. The difference in Sauter diameter, wall temperature, and evaporative heat flux in the HZDR model compared to that of the KTH-WSE model is caused by considering bubble sliding.

**Frequency, nucleation site density and waiting time.** Compared to the departure diameter, the frequency has a much smaller impact. In the HZDR model, the Sauter diameter has nearly the same radial distribution (see Electronic Supplement Fig. E2). The largest difference in the Sauter diameter occurs at the centre line but variations from the baseline model remain weak. There is no impact of frequency shown in the KTH-WSE model. Again, due to the lift force, in the case of 2× of the HZDR model, smaller bubbles are pushed to the surface which leads to a void peak shown on the wall. Meanwhile, the 1/2× and 1× cases have a void fraction peak with a distance near the wall. Generally, the void fraction decreases with a decrease in frequency. In the KTH-WSE model, no impact of frequency on the void fraction distribution is found. Referring to Eq. (12), it is found that the impact of the nucleate site density is nearly the same as that of the frequency for both the HZDR and KTH-WSE models.

The axial distribution of the wall temperature and evaporative heat flux at the wall ($r^+ = 1$) is further evaluated (see Electronic Supplement Fig. E3). Compared to the impact of departure radius, the impacts of frequency and nucleate site density are much smaller, which agrees with the above-mentioned equation. In the HZDR model, the sharp increase in the Sauter diameter occurs later with a decrease in frequency. However, later, all cases have a similar Sauter diameter downstream until the outlet. The wall temperature follows the same rule as the departure diameter, but the difference is smaller at 0.8 K. The same applies to the evaporative heat flux.

In the KTH-WSE model, the axial distribution of the Sauter diameter is not impacted by frequency and nucleate site density. The wall temperature has a similar distribution as the case with different departure diameters, with only a small difference in the positions of peaks. The impact of frequency on the evaporative heat flux is much smaller than

that of the nucleate site density and departure diameter. The difference in the evaporative heat flux between 1/2× and 2× the frequency at the outlet is 6700 W m⁻² while it is 66798 and 31470 W m⁻² for that of departure diameters and nucleate site densities, respectively. Due to waiting time being a parameter related to frequency, 1/2×, 1×, and 2× waiting time is equivalent to 0.55:1:1.67 of frequency, so the analysis will not be repeated again.

**Conclusion.** In the wall boiling model, bubble dynamics parameters such as departure diameter, frequency, nucleate site density and waiting time have a significant impact not only on the near-wall region flow but also on the bulk flow. However, these parameters are empirically correlated or calibrated. In order to reduce the degree of freedom for the calibrations considering other tunable system parameters such as coalescence/break-up factors, group numbers, or others, it is necessary to consider a more physics-based mechanistic model.

## 7. Conclusion

In this paper, the simulation results of 12 participants are compared with each other and validated with experimental data. The participants used different two-fluid closure relations to simulate 30 different Freon (R12) boiling experiments in the DEBORA facility, conducted by CEA Grenoble in the 90's. The experiment consists of a pipe with an internal diameter of 19.2 mm and is heated over 3.5 m. Each case is characterised by the flow rate, inlet temperature, wall heat flux and outlet pressure. Operating conditions range with pressure from 14 bar to 26 bar and flow rates from 2000 kg m⁻² s⁻¹ to 5000 kg m⁻² s⁻¹. Exit quality ranges from single-phase conditions to $x = 0.1$ which leads to a peak void fraction of $\alpha = 70\%$. Due to the high pressure, bubbles remain small, and there is no transition from bubbly flow regime (Hösler, 1968; François et al., 2011). However, different kind of bubbly flows are observed: wall-peak, intermediate peak and core-peak, depending on the case considered. Measurements along the pipe radius near the end of the heated section are compared to code predictions. They include void fraction, bubble mean diameter, vapour velocity and liquid temperature.

This benchmark covered two phases. In the first open phase, participants could fine-tune their models to the experimental results, while in the second blind phase no experimental results were provided. The participants used the same models for open and blind phases.

Overall, the benchmark involved very different closures and a wide range of models' complexity was covered. Yet, it is extremely difficult to have a robust set of closures for a wide range of operating conditions, even when tuned to the experimental results. Predictions of near-wall or cross-section averaged void fractions are strongly scattered. The wall-to-core peak transition is not captured consistently by the models. The degree of subcooling, the bubble diameters, the void fraction distribution and its area-averaged mean value are also poorly predicted in many operating conditions. We were unable to show significant advantage of any particular model. The interactions between mechanisms and their hierarchy are difficult to understand. For instance, great care and detailed models can be put into the interfacial area modelling, with up to two-velocity groups, or in the wall heat-flux partitioning models. However, the more complicated models get difficult to calibrate as more and more information and a priori knowledge of the flow is required (e.g., the apparent contact angle in realistic flow conditions). As part of the sensitivity analyses, it was identified that user and code effects might still have some impact on the presented results, while mesh effects were found to be negligible. The complexity of the involved models makes it challenging to track all potential model implementation discrepancies, leading to code prediction differences based on a seemingly identical set of closure models. no model could clearly demonstrate superiority or that an element is a key to success. Instead, global calibration is essential and the model's overall capabilities are the results of the entirety of closure relationships.
In addition, it means that the experimental information is not sufficient to discriminate models and to prevent over-fitting or model's compensations. Therefore, in order to progress further in models development and validation, it is compulsory to have new measurements. Optical measurements of nucleation site density, bubble departure diameter and frequency could be used directly in the wall boiling model, e.g., MIT has performed such measurements for water at 10 bar (Kosso-lapov et al., 2024). The measurement of the temperature distribution on the heat transfer surface can improve the heat flux partitioning model, e.g., regarding the validity of the quenching heat flux, a new concept of heat flux due to sliding bubbles on the wall was introduced, etc. To improve the models of bubbly flow (bubble forces, coalescence and break-up) and the interfacial area models, measurements of the local bubble size distribution are required. Finer analysis of the void-fraction probes signal to extract probability density distribution is needed to better understand the flow topology and to really assess the relevance of iMUSIG approaches for these kind of flows. They are the key to discriminate performances of interfacial area modelling (IATE, MUSIG or iMUSIG models). Furthermore, in the case of a boiling flow measurement similar to DEBORA, it would be ideal to measure liquid temperature, phase velocities and void fraction simultaneously.

Besides, all liquid velocity predictions remain in a rather classical single-phase profile, even for the highest void-fractions whereas strong buoyancy effects were to be expected; liquid velocity measurements (including characterisation of the fluctuations) would help to better characterise and model bubble induced turbulence and its overall impact. Also, simultaneous measurements of liquid temperature with the other variables with close analysis on uncertainties would also be valuable to assess energy budget and consider the flow from a global point of view.

Another prospect is the prediction of Departure from Nucleate Boiling (DNB) in high-pressure conditions. Although DNB was not considered in this benchmarking exercise, it is expected that DNB predictions at high-pressure conditions depend strongly on the near-wall flow, temperature, and void fraction distributions. Therefore, the suitability of the closures also limits the accuracy of DNB predictions. These conditions remain very challenging for current models. Indeed, near-wall predictions of void-fraction are very dispersed; besides, none of the models was capable of predicting the sharp change in void-fraction profile in condition B27 that marks its close proximity to the critical heat flux. Further progresses are still required.

Following this benchmark impulse, it would be a success if future activities were engaged on high-pressure boiling water experiments with a continuation of models' comparisons and development. A lot of efforts is still required to develop a reliable two-fluid model with predictive capabilities.

### CRediT authorship contribution statement

**G. Bois:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Project administration, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.
**P. Fillion:** Writing – review & editing, Project administration, Funding acquisition.
**F. François:** Writing – review & editing, Writing – original draft, Validation, Investigation, Formal analysis, Conceptualization.
**A. Burlot:** Writing – review & editing, Validation, Investigation, Formal analysis, Conceptualization.
**A. Ben Hadj Ali:** Writing – review & editing, Supervision, Investigation.
**A. Khaware:** Writing – original draft, Investigation.
**J. Sanyal:** Supervision, Software.
**M. Rehm:** Investigation, Formal analysis, Data curation.
**B. Farges:** Writing – review & editing.
**F. Vinauger:** Writing – review & editing, Supervision.
**W. Ding:** Writing – review & editing, Supervision, Software, Investigation.
**A. Gajšek:** Writing – original draft, Software, Investigation.
**M. Tekavčič:** Writing – review & editing, Supervision, Software.
**B. Končar:** Writing – review & editing, Supervision.
**J.-M. Le Corre:** Writing – review & editing, Writing – original draft, Validation, Supervision, Software,

**Conceptualization.** **H. Li:** Validation, Supervision, Software.
**R. Härlin:** Validation, Software.
**J. Jaseliūnaitė:** Writing – review & editing, Writing – original draft, Validation, Software, Investigation.
**E. Baglietto:** Writing – review & editing, Supervision, Investigation.
**R. Brewster:** Writing – original draft, Writing – review & editing, Software, Investigation.
**A. Ding:** Writing – original draft, Software, Investigation.
**D. Vlček:** Writing – original draft, Software, Investigation.
**Y. Sato:** Writing – review & editing, Supervision, Software.
**J. Xiong:** Writing – review & editing, Supervision, Conceptualization.
**H. Wang:** Writing – original draft, Software, Investigation, Formal analysis.
**H. Luo:** Validation.
**L. Vyskocil:** Writing – review & editing, Writing – original draft, Investigation.
**V. Hovi:** Writing – review & editing, Writing – original draft, Validation, Software, Investigation.

### 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

Our data and processing code are shared at zenodo. https://doi.org/10.5281/zenodo.12551306.

### Acknowledgements

For the purpose of Open Access, a CC-BY public copyright licence 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_20.png" /></figure>

https://creativecommons.org/licenses/by/4.0/

The authors are grateful to the NEPTUNE project (gathering Electricité de France (EDF), Commissariat à l'Energie Atomique (CEA), Institut de Radioprotection et Sûreté Nucléaire (IRSN) and Framatome) for the financial support of the benchmark meeting and organisation as well as for providing new experimental data.

Authors from the Jožef Stefan Institute gratefully acknowledge the financial support provided by the Slovenian Research Agency through the grants P2-0026 and P2-0405.

VTT's work has been partly funded by the Finnish State Nuclear Waste Management Fund (VYR) as a part of the SAFIR2022 Programme (The Finnish Research Programme on Nuclear Power Plant Safety 2018–2022).

## Appendix A. Additional information on each set of models

This appendix provides additional information on each specific set of models.

### A.1. ANSYS model

The implicit Eulerian Multiphase Model in Ansys Fluent is used to investigate the DEBORA boiling cases at different operating conditions. The sub-grid wall mass transfer is modelled by an extended wall boiling model based on the RPI Wall boiling model (Kurul and Podowski, 1991) which accounts for the convective heat transfer to the vapour phase. This transfer is triggered by a volume fraction dependent blending function after Lavieville (Lavieille et al., 2006). To obtain grid independent results at the wall, Egorov's quenching correction method is applied. This method relies on the self-similarity of the temperature profile and leverage the turbulent thermal wall function to have an estimate for the liquid temperature at a fixed non-dimensional wall distance (Egorov and Menter, 2004).
The deployed RPI sub-model parameters, nucleation site density and bubble departure diameter, follow Krepper and Rzehak calibration efforts to model the DEBORA boiling cases (Krepper and Rzehak, 2014). The boiling model parameters were tuned to predict temperature profiles with higher accuracy. Accordingly, a modified bubble departure diameter by Tolubinsky and Kostanchuk (1970) was utilised. The reference diameters $d_{ref}$ were set to 0.35 mm and 0.24 mm at the operating pressure levels of $p = 1.46$ MPa and 2.6 MPa respectively. The Lemmert and Chawla model (Lemmert and Chawla, 1977) was utilised for modelling the nucleation site density. The model reference nucleation site densities $C_{ref}$ were set to $5 \times 10^6$ m$^{-2}$ and $3 \times 10^7$ m$^{-2}$ at the operating pressure levels of $p = 1.46$ MPa and 2.6 MPa respectively. The liquid-vapour mass transfer in the bulk flow is coupled to a heat balance at the interface. The heat transfer resistance on the liquid side of the interface is modelled after Ranz and Marshall (1952) whereas the constant time scale return to saturation method proposed by Laviéville et al. (2006) is utilised for the vapour phase resistance.

The interfacial momentum exchange included the two components of the instantaneous drag force, namely drag force and turbulent drift force, as well as turbulence interactions and the lift force. The Ishii and Zuber model is used for modelling the drag which considers the change in the gas bubble shape as a function of the bubble Reynolds and Eötvös numbers (Ishii and Zuber, 1979). The SST $k-\omega$ turbulence model was used to predict the liquid eddy viscosity. The turbulence in the gaseous phase has been neglected but its impact on the liquid phase has been considered. The continuous liquid phase fluctuations cause the gaseous particles to be transported from regions of higher concentration to regions of lower concentration. These effects are modelled with the Simonin turbulent dispersion force approach (Simonin and Viollet, 1990). The Sato enhanced viscosity approach is used to model turbulence modulation resulting from the interaction of dispersed bubbles with the continuous liquid phase (Sato and Sekoguchi, 1975). Finally, the momentum exchange due to lift force is modelled according to Hessenkemper et al. (2021), Ansys (2022).

The homogeneous discrete population balance method is used to accurately model interfacial area density $a_l$. This method assumes that the gas bubbles exist across a range of sizes and a transport equation is solved for each of these size fractions. The population balance model is coupled to the bulk flow through a Sauter mean diameter. In the current simulations all vapour size fractions are advected with the same vapour velocity. The vapour phase consists of 15 bins with sizes linearly distributed between 0.1 mm to 2.0 mm. The aggregation model of Prince and Blanch (Prince and Blanch, 1990) and the break-up model of Luo and Svendsen (1996) were chosen to model breakage and coalescence between the classes.

Tabulated R12 material properties, generated from NIST REFPROP library using Ansys Software, are used to provide the refrigerant thermochemical, transport and saturation properties. The governing equations are solved within the 2D axisymmetric double precision steady solver in Ansys Fluent 2022R2. To achieve a perfect mass balance, it has been verified that the sum of all bin fractions source terms are equal to the interfacial mass transfer.

## A.2. CEA model

CEA contribution is performed with Neptune_CFD V7.2.0 and the reference models from the V&V (Verification & Validation) (Guelfi et al., 2007; Douce et al., 2010; Mimouni et al., 2008, 2009, 2010a,b, 2011). In particular, a second order turbulence model based on SSG (Speziale et al., 1991) is used in the liquid phase with interfacial source terms $\Pi_{ij}^i$ for Reynolds stresses and $P_e^i$ for dissipation (Bois, 2017). Modified wall laws are used in the heated region to consider the impact of nucleating bubbles as an additional roughness (Mimouni et al., 2010a). Momentum transfers are described in Table 8 with a lift force adapted from Tomiyama et al. (Eq. 11, 2002) to use the maximum horizontal dimension of a deformed bubble given by Wellek

et al. (Eq. 34, 1966) to compute the modified Eötvös number $E_{o_d}$. Besides, the coefficient originally proposed only for $E_{o_d} < 4$ is used for any diameter:

$$C_L = \min[0.288 \tanh(0.121 Re_b(d_{32})) ; f(\min(E_{o_d}, 10.7))] \quad (A.1)$$

Turbulent dispersion is slightly different from Laviéville et al. (Eq. 13, 2015) and is expressed as:

$$F_{TD} = -\frac{2}{3} c_{TD} \rho_l k_l \nabla \alpha \quad \text{with:} \quad (A.2)$$

$$C_{TD} = (F_d \tau^l - 1) \left( \frac{b + \eta_r}{1 + \eta_r} + \frac{\alpha_g}{1 - \alpha_g} \right) + C_{AM} \frac{b^2 - b}{1 + \eta_r} + \frac{\rho_g}{\rho_l} \frac{b^2 + \eta_r}{1 + \eta_r}$$

$$\eta_r = \frac{\tau_{lg}^l}{\tau_{lg}^F}, \quad b = \frac{\rho_l + \rho_l C_{AM}}{\rho_g + \rho_l C_{AM}}, \quad \tau_{lg}^l = \frac{3}{2} C_\mu \frac{k_l}{\epsilon_l} \left( 1 + \beta \frac{V_r^2}{k_l} \right)^{-1/2},$$

$$\tau_{lg}^F = \frac{1}{F_d} \left( \frac{\rho_g}{\rho_l} + C_{AM} \right)$$

with $\beta = 2.7$ and $F_d = 3C_d \bar{u}_r/(4d_{32})$.

The closure of the IAC relies on a transport equation (Ruyer et al., 2007; Morel et al., 2010) that considers a parabolic bubble-size distribution and integrates (Prince and Blanch, 1990) model for coalescence and Wu et al. (1998) for break-up.

Interfacial heat transfers are described in Table 10 and detailed in Yao and Morel (2004), Favre et al. (2023). The non-dimensional liquid temperature $T^l$ in the first wall cell is modelled following a Van Driest formulation. The WHP model in Table 11 (Mimouni et al., 2016a,b) is combined with a wall function (Mimouni et al., 2010a) using Van Driest wall-law to assess liquid temperature $T^l$ and velocity $\bar{u}_l^l$ at a distance $d_d$ from the wall. The 3 classical fluxes introduced in Kurul and Podowski (1990, conduction, quenching and evaporation) acts over the fraction $f_\alpha$ of the surface whereas the remaining flux $q_{iv}^v$ is sent directly to the vapour for higher void fractions:

$$q_{iv}^v = h_{iv}^v (T_w - T^v); \quad f_\alpha = \begin{cases} 1 & \text{if } \alpha < 0.75; \\ 1 - 200(0.75 - \alpha)^2 & \text{if } 0.75 \le \alpha < 0.8; \\ 200(0.85 - \alpha)^2 & \text{if } 0.8 \le \alpha < 0.85; \\ 0 & \text{if } 0.85 \le \alpha. \end{cases} \quad (A.3)$$

Bubble sliding is not considered. The void fraction $\alpha$, physical properties and the vapour temperature $T_v$ are taken at the centre of the wall-adjacent cell. The heat transfer coefficient is computed by $h_v^v = 2\lambda_v/dr$ where $dr$ is the cell size in the wall normal direction. The departure diameter (or lift-off) $d_d$ is adapted from Unal (Eq. 5 and 13, 1976) to cope with saturated flow conditions by modifying $\Delta T_{sub}$ in their Eq. (5.2) as:

$$\Delta T_{sub} = \begin{cases} \max(T^{\text{sat}} - T^l, 1.0 \times 10^{-8}) & \text{if } 0 \le T \le T_{lim}; \\ q_{lw} / (\rho_l c_{pl} \bar{u}_l^l \text{St}_{lim}) & \text{otherwise.} \end{cases} \quad (A.4)$$

where the Stanton limiter $St_{lim} = 0.0065$ and $St = q_{lw} / (\rho_l c_{pl} \bar{u}_l^l d T)$.

A steady algorithm is used for time advancement of the solution with a velocity prediction algorithm coupling mean and relative velocities. Physical properties are computed from local thermodynamic parameters using in-house EOS for R12 (Cathare). Inlet turbulent profile (at $z = 0$ m) is established through a recycling boundary condition.

## A.3. FRA model

The computational domain for the DEBORA test was meticulously modelled using a sophisticated 3D trimmed mesh comprising approximately 19 million cells. A three-layer prism mesh was employed at the wall to enhance resolution and accuracy. The axial domain was extended by 1 m upstream and 0.5 m downstream from the heated section to ensure comprehensive flow development and capture of downstream effects. To accurately represent the adiabatic portion of the tube wall, a smooth wall assumption was utilised. Conversely, for the heated wall,
the Gilman (2014) wall roughness model was applied. This model was specifically chosen to account for the increased roughness resulting from bubbles adhering to the wall surface. The standard $k-\epsilon$ turbulence model was selected for its robustness and reliability in simulating turbulent flows. Flow dynamics were further refined by applying a power law shape wall damping for the lift force. Bubble dynamics, including size, break-up, and coalescence, were modelled using the IATE as proposed by Yao and Morel (2004). The turbulence coalescence coefficient $C_{bk} = 2.86$ and break-up coefficient $C_{bk} = 1.6$ were used as recommended by Yao and Morel (2004). For interfacial heat transfer, the Ranz-Marshall correlation was employed using a dispersed phase Nusselt number $Nu = 26$. Wall boiling phenomena were modelled using several well-established models: the Bubble Departure Diameter model by Unal (1976), the Bubble Departure Frequency model by Cole (1960), and the Nucleation Site Density model by Kurul and Podowski (1990). The Del Valle and Kenning (1985) model was employed for calculating the quenching heat flux calculation. The interfacial forces acting on bubbles were modelled using the following set of models: Ishii and Zuber (1979) model for drag forces, the Tomiyama et al. (2002) model for lift forces, the Antal et al. (1991) model for wall lubrication forces, and the Zuber (1964) model for virtual mass forces. Additionally, the turbulent dispersion force was accounted for using the inertial correction model from Software (2023). The latter model modifies the turbulent Prandtl number by incorporating particle relaxation time, eddy-turnover time, and the virtual mass coefficient. For material properties of R12, data was sourced from the comprehensive NIST database (Lemmon et al., 2023), ensuring the use of accurate and up-to-date thermophysical properties.

## A.4. HZDR model

Since many sub-models for the heat flux calculation applied in the conventional wall boiling model such as nucleation site density, bubble departure diameter and bubble release frequency (Setoodeh et al., 2021a) were derived from pool boiling experiments at ambient pressure (Kurul and Podowski, 1990; Tolubinsky and Kostanchuk, 1970), the universal application of the CFD approach is hindered which usually requires case dependent recalibration. Neglecting the bubble sliding on a heating wall is another debatable issue, while a bubble sliding on a heating wall disrupts the thermal layer in the wall region and enhances the single-phase heat transfer (Basu et al., 2005). To improve the applicability and the prediction accuracy of the CFD approaches in wall boiling flow simulation and to reduce the case dependency, mechanistic sub-models for bubble dynamics and wall partitioning models were developed in HZDR (Ding et al., 2018, 2019; Setoodeh et al., 2021b,a). This Eulerian-Eulerian MCFD approach includes three parts: a mechanistic model for bubble dynamics prediction, an updated heat partitioning model to consider the sliding and a nucleation site density model to consider the cavity size dependency on bubble activation superheat. The mechanistic model, which is based on the balance of forces applied on a single bubble, considers the evaporation of the microlayer underneath the bubble, thermal diffusion and condensation around the bubble cap as well as geometry changes and dynamic inclination and contact angles. The target of this model is to reduce the requirement of a recalibration of bubble dynamics-related parameters and consider the missing sliding in the conventional wall boiling models. However, because the bubble mechanistic model is dependent on the onset temperature, an adaptation of the current nucleation site activation model is required which discretises the conventional nucleation site density/superheat curve into several groups (Ding et al., 2019). Furthermore, we employ the population balance model and velocity groups (iMUSIG) in the approach. The critical diameter that separates velocity groups is based on the lift force from Tomiyama et al. to 1.5 mm for 14.6 bar and 1 mm for 26.5 bar (Krepper et al., 2013). Validation was made with experimental data from flow in a vertical pipe and an annulus. The experimental test cases cover a wide range of

operating parameters such as wall heat flux, fluid velocity, subcooling temperature and pressure (Setoodeh et al., 2021a). The approach was later extended and updated to consider the flow pattern development with the GENTOP (GENeralized-TwO Phase) concept (Setoodeh et al., 2022).

## A.5. JSI model

JSI used the Eulerian multiphase model in ANSYS Fluent version 2021 R2 with the standard RPI wall boiling model and population balance model for the bubble sizes. The $k-\omega$ SST (Shear Stress Transport) of Menter (1994) was used for the turbulence in liquid and the dispersed phase zero equation turbulence model is set for the vapour. Bubble induced turbulence was considered with the model of Sato and Sekoguchi (1975).

The following models for interfacial forces available in Fluent were used: Ishii and Zuber (1979) for the drag model, de Bertodano (1991) for the turbulent dispersion force, and a modified lift force model based on Tomiyama (1998). Wall lubrication force was not included. The symmetric model in Fluent was used for the IAC.

Initial simulations demonstrated that in certain cases there is not enough vapour transport into the core flow compared to the measurements. Based on these observations and performed parametric studies, the lift force coefficient was modified as

$$C_L = \begin{cases} 0.3 & \text{for } Eo' \le 3; \\ 0 & \text{for } 3 \le Eo' \le 7; \\ -0.3 & \text{for } Eo' \ge 7; \end{cases} \quad (A.5)$$

where $Eo'$ is the modified Eötvös number as defined by Tomiyama (1998). Compared to the original formulation, the key differences are the stepwise change of the lift coefficient and its sign reversal at slightly smaller bubble diameter.

The thermal phase change model with the two-resistance approach for calculation of interfacial heat transfer coefficient was used to model phase change in the flow (bulk flow evaporation or condensation). The model of Ranz and Marshall (1952) is used to calculate the heat transfer from the interface to the liquid. On the other side, the heat transfer from the interface to the vapour is set to a large value $Nu_v = 500$, assuming a saturated vapour phase.

The following parameters of the RPI wall boiling model were used by JSI: the model of Tolubinsky and Kostanchuk (1970) for the bubble departure diameter, Cole (1960) for the bubble departure frequency, Lemmert and Chawla (1977) for the nucleation site density, and Del Valle and Kenning (1985) model for area of bubble influence. Bubble waiting time coefficient was set to one. The quenching heat flux correction in Fluent was enabled with setting a fixed value for dimensionless wall distance $y^+ = 250$ for the liquid temperature used in the RPI model.

A population balance model was used to calculate the mean bubble diameter, using the homogeneous discrete method, with a single vapour velocity field for all bubble size classes. After a parametric study, 24 bubble size classes (bins) were selected, with exponential distribution of bubble sizes from 0.1–3.16 mm. Bubble coalescence and break-up were considered using the model of Liao et al. (2015), with all model coefficients set to one.

## A.6. KTH-WSE model

The two-phase CFD model used by KTH-WSE is based on the Eulerian-Eulerian framework implemented in OpenFOAM-9 multiphaseEulerFoam solver. The selection of closure models was performed following extensive sensitivity studies and optimisation documented in Härlin (2022).

A two-dimensional axisymmetric geometric model (1° angle) is used for the simulations with a 3.5 m uniformly heated section, followed by 0.5 m adiabatic. The computational mesh is made of 400 axial cells
and 40 radial cells with 0.5:1 grading to better capture near-wall flow variations. Pressure and temperature-dependent R12 properties were considered and obtained from Bell et al. (2014).

The final model selection and associated references can be found in Tables 8 to 11. The following additional information, relevant observations and necessary tunings during the model selection are provided below:

* A population balance model (set to 23 size group with 1 velocity group) was required to adequately capture the radial distribution of the bubble Sauter mean diameter. Attempts to use an interfacial area transport equation model were unsuccessful (Härlin, 2022).

* The wall partitioning model of Lavieville et al. (2006) is considered. For the investigated cases, however, the single-phase convective heat flux to the vapour phase is negligible and the heat is partitioned into single phase liquid convection, liquid evaporation and liquid quenching.

* The reference diameter $d^{Ref}$ of the bubble departure diameter model of Tolubinsky and Kostanchuk (1970) was set according to the equation from Kocamustafaogullari and Ishii (1995):

$$d^{Ref} = \sqrt{2.64 \times 10^{-5} \theta \left(\frac{\sigma}{g\Delta\rho}\right)^{0.5} \left(\frac{\Delta\rho}{\rho_v}\right)^{0.9}} \quad (A.6)$$

where $\Delta\rho = \rho_l - \rho_v$ and the contact angle $\theta$ is set to 13.5° through sensitivity studies (Härlin, 2022).

* The reference nucleation site density $N_b^{ref}$ of the nucleation site density model of Lemmert and Chawla (1977) was set according to the system pressure: $N_b^{ref} = 5 \times 10^6 \text{ m}^{-2}$ for pressure near 1.46 MPa and $N_b^{ref} = 3 \times 10^6 \text{ m}^{-2}$ for pressure near 2.62 MPa. This is one of the only ad-hoc parameters of the KTH-WSE model, which showed large sensitivity to the results and required careful tuning (Härlin, 2022).

* The quenching heat flux is based on the typical model of Del Valle and Kenning (1985) with a waiting time ratio of 0.8 and a liquid temperature evaluated at $y^+ = 250$.

* The thermal wall function of Jayatilleke (1966) was used (OpenFOAM default).

* A constant Nusselt number $Nu_v = 100$ was imposed for the interfacial heat transfer of the vapour phase, meant to be representative of no thermal resistance. Higher values led to code non-convergence, however a value of 100 is estimated to be sufficiently high (as compared to the thermal resistance on the liquid side set by Ranz and Marshall, 1952).

### A.7. LEI model

OpenFOAM v8 software and the models available in this code were used for the benchmark exercises. Turbulence in liquid was modelled by employing the $k-\omega$ SST model for dispersed bubbly flows with Sato bubble-induced turbulent viscosity model while the $k-\epsilon$ model for the gas phase in a two-phase system supporting phase inversion was used for gas turbulence. The (Ishii and Zuber, 1979) model was used for modelling drag force, Tomiyama et al. (2002) for lift force, Tomiyama (1998) for wall lubrication force with cosine wall damping, and de Bertodano (1991) for turbulent dispersion ($C_{TD} = 1$). No virtual mass model was used. The PBM with 27 size groups was used to calculate the non-uniform bubble size distribution. Size groups ranged from 0.13 mm to 11.4 mm. Only one velocity group was used. Prince and Blanch (1990) model was used for bubble coalescence, with a coefficient of 0.1, while a model of Lehr et al. (2002) was used for binary bubble break-up. The liquid side of heat transfer during condensation or evaporation was evaluated according to Ranz and Marshall (1952). $Nu_v = 10$ was applied for heat transfer from the interface to the vapour. The wall boiling correlations for the RPI model (Kurul and Podowski, 1991) was selected as follows: Tolubinsky and Kostanchuk (1970) for bubble departure diameter, Lavieville et al. (2006) for partitioning, Cole

(1960) for bubble departure frequency, Lemmert and Chawla (1977) for nucleation site density, Del Valle and Kenning (1985) for the area of bubble influence as well as for quenching heat flux.

### A.8. MIT model

The MIT benchmark contribution is based on models developed over the last ten years collectively called the MIT boiling (MITB) model. Much of this work was performed in collaboration with the CASL (Consortium for the Advanced Simulation of Light Water Reactors) program sponsored by the US Department of Energy. Both the momentum closures and the wall boiling representation in the MITB model leverage experimental and numerical databases recently produced by the multiphase heat transfer community to incorporate a more consistent representation of the physical phenomena across a wide range of conditions. Therefore, the DEBORA benchmark represents an excellent opportunity to assess the maturity of the closures without additional reformulation or dedicated calibration.

MIT's DEBORA computational model was constructed using the STAR-CCM+ hexahedral-dominant cell approach. One-quarter of the tube cross-section was modelled, with associated symmetry plane boundary conditions. The model consists of the full test section length plus regions upstream and downstream of the test section for proper flow development and outflow conditions, respectively. Six radial prismatic cell layers were used at the wall, with the core (inner) mesh consisting primarily of hexahedral cells. The overall model consists of approximately 200 000 control volumes, with about 195 000 of these comprising the test section domain and the remainder in the upstream and downstream regions. The thermophysical properties for the R12 refrigerant were obtained from the online NIST Chemistry WebBook, SRD 69 (Lemmon et al., 2023) and input to STAR-CCM+ as tables. Each simulation was run for a fixed number of iterations, with 30 000 iterations sufficient to obtain convergence in all cases.

The calculation of interphase momentum transport requires closures for bubble drag, lift, and other forces, which are generally specified using correlations. The strong effects of contamination on bubble drag were noted by Tomiyama et al. (1998), who developed correlations for “pure”, “slightly contaminated”, and “fully contaminated” systems. Because the fluid conditions for most experiments are unknown, the MITB momentum closures use the slightly contaminated Tomiyama drag correlation (Tomiyama et al., 1998, Eq. 32).

An accurate representation of the lift force is critical because it significantly affects the radial void fraction distribution. In particular, lift force inversion has been found to be strongly affected by bubble shape (Tomiyama et al., 2002; Feng and Bolotnov, 2018) and local vapour concentrations. Sugrue (2017) found that bubble shape effects in a turbulent flow field could be accounted for by a dimensionless “Wobble number” ($Wo$), defined as the product of the Eötvös number ($Eo$) and the liquid's non-dimensional kinetic energy. Sugrue (2017) also noted substantial evidence showing that the lift force behaviour changes significantly when considering mixed-sized bubble groups. As a result, earlier lift closures often considerably overestimated the lift coefficient by disregarding bubble wake interactions in a turbulent flow, limiting the lift force and corresponding lateral vapour migration.

Based on these observations, Sugrue (2017) developed a data-driven lift correlation as the product of two functions. The first function depends on the dimensionless “Wobble” number and accounts for bubble shape and turbulence effects, while the second function accounts for bubble interactions and depends on the void fraction. However, during the DEBORA experimental campaign, it was noted that the classic flow regime classifications do not necessarily apply to the density ratio represented in the experiment, in which the vapour phase remains dispersed up to very high void fractions. Because (Sugrue, 2017) introduced the void fraction correction based on the low-pressure transition from bubbly flow, this correction was eliminated for the DEBORA calculations, and the open phase of the benchmark confirmed this
assumption's applicability. It is also necessary to limit the predicted
strong lift force near walls to prevent nonphysical behaviour. Shaver
and Podowski (2015a) proposed a simple empirical lift coefficient
expression that effectively limits the near-wall lift force where the
velocity gradient of the continuous liquid phase velocity becomes quite
large. The (Shaver and Podowski, 2015a) limiting scheme and the (Sug-
rue, 2017) lift force correlation were used to model the lift force in the
MIT DEBORA simulations.

The balance between the lift force and turbulent dispersion effects is responsible for the lateral void fraction redistribution. The turbulent dispersion derived from the unsteady turbulent component of drag on bubbles by *Burns et al.* (2004) is employed by the MITB model. In addition to the lift and turbulent dispersion forces, a wall lubrication force is often adopted to recover the experimentally-observed wall-peaked void fraction profile for low void fraction conditions. However, *Lubchenko et al.* (2018) showed that this near-wall peaked profile results from the quadratic void fraction distribution caused by bubbles along the wall and proposed a geometric correction arising from the bubbles' spherical shape. After damping the lift force near the wall using the (Shaver and Podowski, 2015a) formulation, the turbulent dispersion is the only remaining interfacial force perpendicular to the wall. Therefore, a wall lubrication force based on the turbulent dispersion force must be added to maintain the desired void fraction distribution, resulting in equilibrium between the two forces. Using the turbulent dispersion force of *Burns et al.* (2004) yields the regularised expression for the wall lubrication force (*Lubchenko et al.*, 2018) used in the MITB model and in the MIT DEBORA calculations.

The MITB boiling closures encompass ten years of advancement in understanding subcooled flow wall boiling phenomena. Gilman and Baglietto (2017) identified all relevant wall phenomena and demonstrated the need to capture the interactions and evolution of the various mechanisms across flow and pressure conditions, where no a priori determination of the heat partition is possible. Gilman and Baglietto (2017) also pointed out the dominant role of single-phase heat transfer enhancement and proposed an initial representation of the various pertinent mechanisms. Eventually, Kommajosyula (2020) developed a complete set of closures, which has been validated across low- and high-pressure databases and is leveraged in this benchmark.

Because they can vary by an order of magnitude or more, modelling of bubble diameters as they depart from the wall is an important component of MCFD modelling. Kommajosyula (2020) developed the data-driven MITB bubble departure diameter correlation using an extensive data set applicable to a wide range of conditions, with significantly improved accuracy compared to previous models.

The bubble departure frequency is the inverse of the time for a bubble to nucleate at a cavity (wait time) plus the time spent by a bubble growing on the heated surface before departure (growth time). The MITB bubble wait time is modelled based on the energy storage limit of the thermal boundary layer (Kommajosyula, 2020). This limit is reached when the energy balance can no longer be sustained, so the excess energy causes bubble nucleation. The MITB bubble growth time uses an improved representation of the microlayer evaporation and the condensation on the bubble cap to compute the bubble departure frequency.

The MITB nucleation site density model is based on the spatial ran-
domness approach proposed by Gilman and Baglietto (2017) to predict
the number of active bubbles. Assuming that an active bubble site
suppresses all active sites within one bubble radius, the probability of
finding an active nucleation site within one bubble radius of a bubble-
generating site can be expressed as an implicit equation (Gilman and
Baglietto, 2017). Kommajosyula (2020) later proposed an analytical
solution to this equation based on Lambert's W-function and adopted
a piecewise approximation of the W-function used in the MITB model.

After a bubble leaves the boiling surface, it can condense in a sub-
cooled flow, so liquid- and vapour-phase heat transfer coefficient mod-
els are needed. The MITB boiling model uses the Ranz–Marshall cor-
relation (Ranz and Marshall, 1952) in the liquid phase and a constant

Nusselt number equal to $Nu_v = 26$ in the vapour phase. Turbulence in both the liquid and vapour phases is modelled using the standard $k-\epsilon$ model with logarithmic wall functions in the near-wall regions. The work of Lubchenko et al. (2018) and Magolan et al. (2019) have consistently demonstrated the applicability of this approach, where the turbulent flow is dominated by the shear contribution, so that bubble-induced turbulence models over-predict the turbulence levels.

The local bubble size (diameter) and interfacial area are used to calculate interfacial forces and heat transfer. The interfacial area is defined as the area per unit volume of the vapour-liquid interface. As bubbles depart the wall, their numbers and size can change due to condensation, coalescence, and breakup. Therefore, accurate bubble size and interfacial area modelling are important aspects of boiling predictions with MCFD codes, and the S-γ model, developed and implemented by Vaesson (1998) and Lo and Zhang (2009), provides a good compromise between accuracy and computational efficiency. The S-γ model adopts the methods of moments (MOM) approach and assumes log-normal distributions to minimise the number of transported moments. This transforms the problem into lower-order moments of the bubble size distribution and computes its parameters. The interfacial area can be readily computed once the bubble sizes and void fractions are calculated. The S-γ model requires wall boiling, bubble breakup, and bubble coalescence models as source terms for the various moments. The wall boiling source can be computed from the wall boiling closures discussed earlier, while details of the bubble breakup and coalescence models are provided by Lo and Zhang (2009).

A.9. *PSI-CTU model*

PSI-CTU propose a boiling flow model that follows the concept of the zero-closure model (Baglietto and Christon, 2013), in which simple and robust closure models are employed. The main feature of the PSI-CTU model is the emphasis on the lift force equation, which is adapted from the zero-closure model of Yoon et al. (2017), and together with the rest of the closure equations and boiling models, it is modified and tuned to the DEBORA experiments.

The Eulerian multiphase model, implemented in the Ansys Fluent release 21.2, is used for the simulations of convective boiling in an axisymmetric coordinate system. The k-ω SST model is used within the framework of the mixture turbulent model, and the bubble-induced turbulence interaction is covered by the Sato model (Sato and Sekoguchi, 1975).

The mass transfer rate in the bulk flow is calculated based on the heat fluxes computed with the two-resistance model. The heat transfer coefficients for the liquid and vapour phases are calculated using the Ranz–Marshall model and the zero resistance model, respectively.

A modified RPI subcooled non-equilibrium boiling model is used (Ansys, 2021), accounting for the convective heat flux between the vapour and the wall. The bubble area of the influence factor is calculated with the Del Valle and Kenning model (Del Valle and Kenning, 1985), and the Cole model (Cole, 1960) is employed for the bubble departure frequency. The nucleation site density is predicted with the unmodified Li model (Li et al., 2018), which is assumed to cover the system pressure from 0.1 MPa to 20 MPa and the wettability effects of the heated surface. To avoid numerical instability issues, a simple bubble departure diameter model, the Tolubinsky and Kostanchuk model (Tolubinsky and Kostanchuk, 1970), is employed with slight modifications: $d_d = \min(0.0014, d_1 \cdot \exp\left[\frac{-\Delta T_{ub}}{45.0}\right])$, where the constants $d_1$ and 45.0 were tuned for atmospheric pressure. The values for higher pressure recommended by Krepper and Rzehak (2011) are not used in the CTU-PSI model, as different closure models are employed. Through sensitivity analysis, we set the new value $d_1 = 5 \times 10^{-5}$ m.

The bubble size evolution is calculated by a discrete MUSIG Popula-
tion Balance Model (PBM) with 20 groups (bins): two bins for smaller
bubbles at diameters of 0.01 mm and 0.05 mm, and 18 equidistant bins
ranging from 0.1 mm to 3.5 mm at intervals of 0.2 mm. The bubble
coalescence model is taken from Prince and Blanch (1990), and the break-up model is from Luo and Svendsen (1996). The coalescence and breakage factors are set to 0.2 and 0.5, respectively, through sensitivity analysis. The interfacial area concentration is defined as $a_i = 6\alpha/d_{32}$, where $\alpha$ is the volume fraction of liquid, and $d_{32}$ the Sauter mean diameter calculated in PBM.

The bubble lift force model proposed by Yoon et al. (2017) is introduced, with the main advantage being its good ability to reproduce the physical effect of bubble migration. The lift coefficient, $C_L$, originally proposed to be 0.025, was modified with multiplication in high and low void regions by the coefficients $C_-$ and $C_+$, respectively, as can be seen in Eq. (A.7). The value of $\alpha_{crit}$ proposed by Yoon et al. (2017) is 0.25 and was not altered in this study.

$$C_L = \begin{cases} 0.025 \cdot C_+ & \text{for } \alpha \le \alpha_{crit} \\ -0.025 \cdot C_- & \text{for } \alpha_{crit} < \alpha \end{cases} \qquad (A.7)$$

The drag force coefficient is calculated with the Ishii-Zuber model (Ishii and Zuber, 1979), the Burns model (Burns et al., 2004) is used for the turbulent dispersion force, the wall lubrication force model implemented in Ansys (2021) is used with the constant coefficient $C_{WL} = 1$, and the virtual mass force model is not used to avoid numerical convergence issues.

In summary, the PSI-CTU model's strength lies in its capability to accurately predict void fractions in the centre of the tube, primarily due to our lift force closure model. However, it has a limitation in that it loses accuracy at very high void fractions, where the boiling regime transitions away from bubble-dominated behaviour, making it unsuitable for all experimental ranges. A more detailed description of the PSI-CTU model and its sensitivity studies are reported in Vlček and Sato (2023).

### A.10. SJTU model

According to the observations of available experiments, a five-component mechanistic heat flux partitioning model is proposed and implemented in OpenFOAM v6, the wall heat flux partitioning model can be given by,

$$q_{tot} = q_c + q_q + q_e + q_{tcsl} + q_{tclf} \qquad (A.8)$$

where $q_c$, $q_q$, $q_e$, $q_{tcsl}$, $q_{tclf}$ denote the heat flux components due to single-phase convection, solid quenching, evaporation, bubble sliding transient conduction and lift-off transient conduction. The bubble's life cycle can be divided into three stages. In the first stage, the bubble grows at the nucleation site and a dry spot appears under the growing bubble. Therefore, heat fluxes resulting the bubble growing and the superheat at the dry spot, i.e., solid quenching heat flux $q_q$, should be considered. In the second stage, the bubbles can slide along the heated wall, thermal boundary layer can be disturbed and the bubble grows during this period. The additional evaporation heat flux and transient conduction heat flux $q_{tcsl}$ should be considered. In the final stage, the bubble lifts off the heated wall and the subcooled water fills back the space lifted by the bubbles. The transient conduction flux $q_{tclf}$ induce by bubble lift-off should be considered. The evaporation heat flux $q_e$ can also be calculated by sum of the evaporation in the first and second stage. Some more details can be found in Wang et al. (2023). Force balance model exerting on the nucleated bubbles is employed to analyse the bubble departure diameter and lift-off diameter, according to the wall function in CFD simulations the sliding length can be further obtained.

The bubble diameter is simulated by iMUSIG model, Guo et al.'s Coalescence model (Guo et al., 2016) and Luo and Svendsen's break up model (Luo and Svendsen, 1996) is employed to simulate the bubble coalescence and break-up processes. The critical diameter separating velocity groups is set at 1.5 mm for cases at 14 bar and 1 mm for cases at 26 bar. Ranz-Marshall correlation is utilised to calculate the interfacial

heat transfer coefficient and mass transfer rate, e.g., condensation and evaporation processes. Tomiyama drag force model (Tomiyama et al., 1998), Tomiyama model (Tomiyama et al., 2002) for lift force with the cosine damping function and Burns model (Burns et al., 2004) for turbulent dispersion force are employed to model the interphase forces. The damping function in Eq. (9) is defined by:

$$F_{damp}(y) = \begin{cases} 0.0, & y/\Delta y \le 1. \\ \frac{1}{2} \left( 1 - \cos \left[ \pi \min \left( \frac{y-\Delta y}{d_{32}C_D}, 1 \right) \right] \right), & y/\Delta y > 1 \end{cases} \qquad (A.9)$$

where $\Delta y$ represents the height of the first cell near wall. With the damping function, the lift force does not act in the first layer of grid cells. In the direction away from the heated wall, the value of the damping function gradually tends to 1. The standard $k-\epsilon$ model is utilised to simulate the turbulence. To account for the bubble-induced turbulence, the Sato-Sekoguchi model (Sato and Sekoguchi, 1975) is employed to calculate the total turbulent viscosity. The two-phase wall function is applied to consider the effect of bubble on the wall turbulence based on the bubble lift-off diameter and evaporation heat flux.

### A.11. UJV model

UJV contribution is performed with ANSYS CFX 18.0 code and the models available in this code.

Liquid turbulence is modelled with SST k-ω turbulence model with buoyancy turbulence production and dissipation, and bubble induced turbulence by Sato and Sekoguchi (1975). Dispersed phase zero-equation model is used for vapour turbulence.

Momentum transfer is composed of drag (Ishii and Zuber, 1979), lift and turbulent dispersion forces. The lift force is modelled using Tomiyama (1998) for bubble diameters $d_b \le 1.4$ mm whereas a constant coefficient $C_L = -0.5$ is used for larger bubbles. Turbulent dispersion force is calculated from the Favre averaged drag force model by Burns et al. (2004) with the dispersion coefficient increased to 1.5.

Interfacial heat transfer on the liquid side is modelled with correlation by Ranz and Marshall (1952). Infinite heat transfer coefficient is assumed on the vapour side.

WHP model by Kurul and Podowski (1991) is used with the departure diameter $d_d$ and nucleation site density $N_b^e$ by Kocamustafaogullari and Ishii (1983). Bubble detachment frequency is calculated from correlation by Cole (1960). Bubble waiting time is equal to the reciprocal of the bubble detachment frequency. Quenching heat transfer coefficient is determined from Del Valle and Kenning (1985) correlation. Liquid temperature in the WHP model is evaluated at the fixed dimensionless wall distance $y^+ = 250$.

The closure of the IAC relies on the Inhomogeneous MUSIG model with 10 bubble size groups distributed in interval 0 mm to 2 mm and 2 bubble velocity groups with diameter ranges in 0 mm to 1.4 mm and 1.4 mm to 2 mm. Bubble coalescence is calculated from Prince and Blanch (1990) model with both buoyancy coalescence coefficient and turbulence coalescence coefficient increased to 4. Bubble break-up is determined from the model by Luo and Svendsen (1996).

Physical properties of R12 are computed from the tables by Lemmon et al. (2023). Liquid physical properties are defined as functions of temperature at constant pressure at the test section outlet. Vapour physical properties, saturation temperature and latent heat of vaporisation are all constant, they are determined for the pressure at the test section outlet.

### A.12. VTT model

A development version of the OpenFOAM foundation release, dated April 2022, was employed in the benchmark exercise. The mixture of liquid and gaseous Freon (R12) was described with a continuous liquid phase and two dispersed gas phases (or velocity groups).
Temperature-pressure tabulated material properties for R12 were obtained from NIST. Turbulence was modelled with the Lahey $k-\epsilon$ turbulence model (Lahey, 2005), which includes bubble induced turbulence sources for the liquid phase.

A PBM with 30 size groups in total ranging from 0.1 mm to 19.8 mm with a constant growth rate of 1.2 was selected for the transport of bubble size distributions. The size groups were assigned into the two velocity groups (small and large bubbles) based on the sign of the lift force, with separations at 1.28 mm and 0.892 mm for pressures of 14 bar and 26 bar respectively. The model proposed by Prince and Blanch (1990) was selected for bubble coalescence, with a coefficient of $C_{c0} = 0.356$, and the binary break-up model of Lehr et al. (2002) was selected for bubble break-up.

The interfacial forces were modelled with Ishii and Zuber (1979) for drag, de Bertodano (1991) for turbulent dispersion force ($C_{TD} = 1$), Hosokawa et al. (2002) for wall lubrication force, and Tomiyama et al. (2002) for lift (with cosine wall damping), where the bubble aspect ratio was taken according to Ziegenhein and Lucas (2017).

Liquid side interfacial heat transfer was calculated according to Ranz and Marshall (1952) and the gas side heat transfer was calculated based on a model which applies an analytical solution for heat transfer from the surface of a sphere to the matter within the sphere, accounting for the bubble aspect ratio (Ziegenhein and Lucas, 2017). The gas side heat transfer roughly corresponds to $Nu_v = 10$.

An extended version of the RPI wall boiling model (Kurul and Podowski, 1991) was selected to account for the wall heat transfer. The selected submodels include (Lavieville et al., 2006) for heat flux partitioning, Lemmert and Chawla (1977) for nucleation site density, Kocamustafaogullari and Ishii (1995) for departure frequency, and Tolubinsky and Kostanchuk (1970) for bubble departure diameter. The model parameters for nucleation site density, and bubble departure diameter were taken according to the model of Del Valle and Kenning (1985) with a bubble waiting time coefficient of 0.8. The bulk liquid temperature, used in the WHP submodels, is evaluated at a fixed distance of $y^+ = 250$.

## Appendix B. Participants one-by-one

See Electronic Supplement.

## Appendix C. Cases one-by-one

See Electronic Supplement.

## Appendix D. Code or user effects and mesh convergence

See Electronic Supplement.

## Appendix E. Sensitivity of wall parameters

See Electronic Supplement.

## Appendix F. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.ijmultiphaseflow.2024.104920.

## References

* Ansys, 2021. Ansys fluent® user's guide, release 21.2. ANSYS, Inc.

* Ansys, 2022. Ansys fluent® theory guide, release 2022. ANSYS, Inc.

* Antal, S.P., Lahey, R.T., Flaherty, J.E., 1991. Analysis of phase distribution in fully developed laminar bubbly two-phase flow. Int. J. Multiph. Flow 17, 635-652.

* Baglietto, E., Christon, M.A., 2013. Demonstration and assessment of advanced modeling capabilities for multiphase flow with sub-cooled boiling. CASL-U-2013-0181-001.

* Basu, N., Warrier, G.R., Dhir, V.K., 2005. Wall Heat Flux Partitioning during subcooled flow boiling: Part II - model validation. J. Heat Transfer 127 (2), 141-148. http://dx.doi.org/10.1115/1.1842785.

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

* Berne, 1983. Analyse critique des modèles d'autovaporisation utilisés dans le calcul des écoulements diphasiques en conduites (Ph.D. thesis). École Centrale de Paris.

* Bestion, D., Anglart, H., Caraghiaour, D., Péturaud, P., Smith, B., Andreani, M., Niceno, B., Krepper, E., Lucas, D., Moretti, F., Galassi, M.C., Macek, J., Vyskocil, L., Koncar, B., Hazi, G., 2009. Review of available data for validation of Nuresim two-phase CFD software applied to CHF investigations. Sci. Technol. Nucl. Installations 1-14. http://dx.doi.org/10.1155/2009/214512.

* Bois, G., 2017. Direct numerical simulation of a turbulent bubbly flow in a vertical channel: Towards an improved second-order Reynolds stress model. Nucl. Eng. Des. 321, 92-103. http://dx.doi.org/10.1016/j.nucengdes.2017.01.023.

* Bois, G., Burlot, A., Fillion, P., François, F., 2024. Benchmark DEBORA: experimental and MCFD databank for high-pressure boiling pipe flow. Zenodo, http://dx.doi.org/10.5281/zenodo.12551306.

* Burns, A.D., Frank, T., Hamill, I., Shi, J.-M.M., 2004. The Favre-averaged drag model for turbulent dispersion in Eulerian multiphase flows. In: Proceedings of the 5th International Conference on Multiphase Flow. ICMF'04, Yokohama, Japan, Japan, pp. May 30-June 4.

* Cole, R., 1960. A photographic study of pool boiling in the region of the critical heat flux. AIChE J. 6 (4), 533-538.

* Cole, R., 1974. Boiling nucleation. In: Advances in Heat Transfer, vol. 10, Elsevier, pp. 85-166.

* de Bertodano, M.L., 1991. Turbulent bubbly flow in a triangular duct (Ph.D. thesis). Rensselaer Polytechnic Institute, Troy, New York.

* de Bertodano, M.L., Saif, A.A., 1997. Modified k-epsilon model for two-phase turbulent jets. Nucl. Eng. Des. 172 (1-2), 187-196.

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

* Dhir, V.K., 1998. Boiling heat transfer. Annu. Rev. Fluid Mech. 30 (1), 365-401. http://dx.doi.org/10.1146/annurev.fluid.30.1.365.

* Ding, W., Krepper, E., Hampel, U., 2018. Quantitative prediction of critical heat flux initiation in pool and flow boiling. Int. J. Therm. Sci. 125, 121-131. http://dx.doi.org/10.1016/j.ijthermalsci.2017.11.022.

* Ding, W., Krepper, E., Hampel, U., 2019. The implementation of an activated temperature-dependent wall boiling model in an Eulerian-Eulerian Computational Fluid Dynamics approach for predicting the wall boiling process. Nucl. Technol. 205, 23-32. http://dx.doi.org/10.1080/00295450.2018.1496693.

* Douce, A., Baudry, C., Guingo, M., Laviéville, J., Mimouni, S., Boucker, M., 2010. Validation of NEPTUNE_CFD 1.0.8 for adiabatic bubbly flow and boiling flow. In: Proceedings of the CFD4NRS-3 Conference. Washington D.C., USA.

* Drew, D.A., 1978. The force on a small sphere in slow viscous flow. J. Fluid Mech. 88 (2), 393400. http://dx.doi.org/10.1017/S00221107208002177.

* Egorov, Y., Menter, F., 2004. Experimental implementation of the RPI wall boiling model CFX-5.6. In: Technical Report Ansys/TR-04-10. Technical Report, Ansys GmbH.

* Estrada-Perez, C.E., Hassan, Y.A., 2010. PTV experiments of subcooled boiling flow through a vertical rectangular channel. Int. J. Multiph. Flow 36 (9), 691-706. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2010.05.005.

* Favre, L., 2023. Modélisation et simulation de la crise d'ébullition dans les REP à l'échelle CFD (Ph.D. thesis). Université de Toulouse, URL https://oatao.univ-toulouse.fr/29686/.

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

* Feng, J., Bolotnov, I.A., 2018. Effect of the wall presence on the bubble interfacial forces in a shear flow field. Int. J. Multiph. Flow 99, 73-85.

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

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

* François, F., Garnier, J., Cubizolles, G., 2003. A new data acquisition system for binary random signal application in multiphase flow measurements. Meas. Sci. Technol. 14 (7), 929. http://dx.doi.org/10.1088/0957-0233/14/7/306.

* François, F., Kledy, M., Barre, S., Djeridi, H., 2020. Simultaneous measurements of liquid velocity and temperature for high pressure and high temperature liquid single phase flows. Nucl. Eng. Des. 362, 110587. http://dx.doi.org/10.1016/j.nucengdes.2020.110587.

* Frank, T., Zwart, P.J., Krepper, E., Prasser, H.-M., Lucas, D., 2008. Validation of CFD models for mono- and polydisperse air-water two-phase flows in pipes. Nucl. Eng. Des. 238 (3), 647-659. http://dx.doi.org/10.1016/j.nucengdes.2007.02.056.
Gajšek, A., Tekavčič, M., Končar, B., 2023. Parametric study of population balance model on the DEBORA flow boiling experiment. Nucl. Eng. Technol. http://dx.doi.org/10.1016/j.net.2023.10.039.

Garnier, J., Manon, E., Cubizolles, G., 2001. Local measurements on flow boiling of refrigerant 12 in a vertical tube. Multiph. Sci. Technol. 13 (1&2), URL http://www.dl.begellhouse.com/journals/5af8c23d50e0a883,18c85f7f5866b59a, 3dd586c1794b1dcf.html.

Gilman, L.A., 2014. Development of a general purpose subgrid wall boiling model from improved physical understanding for use in computational fluid dynamics (Ph.D. thesis). Massachusetts Institute of Technology.

Gilman, L., Baglietto, E., 2017. A self-consistent, physics-based boiling heat transfer modeling framework for use in computational fluid dynamics. Int. J. Multiph. Flow 95, 35–53. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2017.04.018.

Guelfi, A., Bestion, D., Boucker, M., Boudier, P., Fillion, P., Grandotto, M., Herard, J.M., Hervieu, E., Peturaud, P., 2007. NEPTUNE: A new software platform for advanced nuclear thermal hydraulics. Nucl. Sci. Eng. 156 (3), 281–324. http://dx.doi.org/10.13182/NSE05-98.

Guo, X., Zhou, Q., Li, J., Chen, C., 2016. Implementation of an improved bubble breakup model for TFM-pbm simulations of gas-liquid flows in bubble columns. Chem. Eng. Sci. 152, 255–266.

Härlin, R., 2022. Two-phase CFD modelling and validation: SH204X master thesis project report. (Master Thesis). In: TRITA-SCI-GRU, (2022:033), KTH, Nuclear Engineering, p. 101.

Hessenkemper, H., Ziegenhein, T., Rzehak, R., Lucas, D., Tomiyama, A., 2021. Lift force coefficient of ellipsoidal single bubbles in water. Int. J. Multiph. Flow 138, 103587. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2021.103587.

Hibiki, T., Ishii, M., 2003. Active nucleation site density in boiling systems. Int. J. Heat Mass Transfer 46 (14), 2587–2601.

Hösler, E.R., 1968. Flow patterns at high pressure two-phase (steam-water) flow with heat addition. AIChE Symp. Ser. 64 (82), 54–66.

Hosokawa, S., Tomiyama, A., Misaki, S., Hamada, T., 2002. Lateral migration of single bubbles due to the presence of wall. In: Fluids Engineering Division Summer Meeting, Volume 1: Fora, Parts A and B, pp. 855–860. http://dx.doi.org/10.1115/FEDSM2002-31148.

Hughmark, G.A., 1967. Mass and heat transfer from rigid spheres. AIChE J. 13 (6), 1219–1221.

Ishii, M., Hibiki, T., 2015. Thermo-fluid dynamics of two-phase flow. 2nd Edition, Springer Springer, pp. 1–462. http://dx.doi.org/10.1007/978-0-387-29187-1.

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

Jayatilleke, C.L.V., 1966. The influence of Prandtl number and surface roughness on the resistance of the laminar sub-layer to momentum and heat transfer. Department of Mechanical Engineering, Imperial College London, URL https://books.google.com/books?id=4zzZrQEACAAJ.

Kataoka, I., Serizawa, A., 1989. Basic equations of turbulence in gas-liquid two-phase flow. Int. J. Multiph. Flow 15 (5), 843–855. http://dx.doi.org/10.1016/0301-9322(89)90045-1.

Kocamustafaogullari, G., Ishii, M., 1983. Interfacial area and nucleation site density in boiling systems. Int. J. Heat Mass Transfer 26 (9), 1377–1387. http://dx.doi.org/10.1016/S0017-9310(83)00269-6.

Kocamustafaogullari, G., Ishii, M., 1995. Foundation of the interfacial area transport equation and its closure relations. Int. J. Heat Mass Transfer 38 (3), 481–493. http://dx.doi.org/10.1016/0017-9310(94)00183-V.

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.

Kossolapov, A., Hughes, M.T., Phillips, B., Bucci, M., 2024. Bubble departure and sliding in high-pressure flow boiling of water. J. Fluid Mech. 987, A35. http://dx.doi.org/10.1017/jfm.2024.405.

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

Krepper, E., Rzehak, R., 2011. CFD for subcooled flow boiling: Simulation of DEBORA experiments. Nucl. Eng. Des. 241 (9), 3851–3866. http://dx.doi.org/10.1016/j.nucengdes.2011.07.003, Seventh European Commission conference on Euratom research and training in reactor systems (Fission Safety 2009).

Krepper, E., Rzehak, R., 2014. CFD for subcooled flow boiling: Analysis of DEBORA tests. J. Comput. Multiph. Flows 6 (3), 329–359. http://dx.doi.org/10.1260/1757-482X.6.3.329.

Krepper, E., Rzehak, R., Lifante, C., Frank, T., 2013. CFD for subcooled flow boiling: Coupling wall boiling and population balance models. Nucl. Eng. Des. 255, 330–346. http://dx.doi.org/10.1016/j.nucengdes.2012.11.010.

Kurul, N., Podowski, M., 1990. Multidimensional effects in forced convection subcooled boiling. In: Proceedings of the 9th International Heat Transfer Conference. 1, Jerusalem, Israel, p. 21.

Kurul, N., Podowski, M., 1991. On the modeling of multidimensional effects in boiling channels. In: ANS Proceedings of the 27th National Heat Transfer Conference. Minneapolis, MN.

Lahey, R.T., 2005. The simulation of multidimensional multiphase flows. Nucl. Eng. Des. 235 (10), 1043–1060. http://dx.doi.org/10.1016/j.nucengdes.2005.02.020, Festschrift Edition Celebrating the 65th Birthday of Prof. Richard T. Lahey, Jr..

Lahey, R.T., Baglietto, E., Bolotnov, I.A., 2021. Progress in multiphase computational fluid dynamics. Nucl. Eng. Des. 374, 111018. http://dx.doi.org/10.1016/j.nucengdes.2020.111018.

Lamb, H., 1945. Hydrodynamics. University Press.

Laviéville, J., Mérigoux, N., Guingo, M., Baudry, C., Mimouni, S., 2015. A generalized turbulent dispersion model for bubbly flow numerical simulation in NEPTUNE_CFD. In: NURETH-16. pp. 4167–4181.

Laviéville, J., Quemerais, E., Mimouni, S., Boucker, M., Mechitoua, N., 2006. NEPTUNE_CFD V1.0, theory manual. (Technical Report). EDF H-181-2006-04377-EN, Nept_2004_L1.

Lehr, F., Millies, M., Mewes, D., 2002. Bubble-size distributions and flow fields in bubble columns. AIChE J. 48 (11), 2426–2443.

Lemmert, M., Chawla, J.M., 1977. Influence of flow velocity on surface boiling heat transfer coefficient. Heat Transf. Boil. 237 (247).

Lemmon, E.W., Bell, I.H., Huber, M.L., McLinden, M.O., 2023. Thermophysical properties of fluid systems. In: Linstrom, P., Mallard, W. (Eds.), National Institute of Standards and Technology, Gaithersburg MD, 20899, http://dx.doi.org/10.18434/T4D303.

Li, Q., Jiao, Y., Avramova, M., Chen, P., Yu, J., Chen, J., Hou, J., 2018. Development, verification and application of a new model for active nucleation site density in boiling systems. Nucl. Eng. Des. 328, 1–9. http://dx.doi.org/10.1016/j.nucengdes.2017.12.027.

Liao, Y., Oertel, R., Kriebitzsch, S., Schlegel, F., Lucas, D., 2018. A discrete population balance equation for binary breakage. Internat. J. Numer. Methods Fluids 87 (4), 202–215. http://dx.doi.org/10.1002/fld.4491.

Liao, Y., Rzehak, R., Lucas, D., Krepper, E., 2015. Baseline closure model for dispersed bubbly flow: Bubble coalescence and breakup. Chem. Eng. Sci. 122, 336–349.

Lo, S., Zhang, D., 2009. Modelling of break-up and coalescence in bubbly two-phase flows. J. Comput. Multiph. Flows 1 (1), 23–38.

Lubchenko, N., 2000. Near-wall modeling of bubbly flows (Ph.D. thesis). Massachusetts Institute of Technology, URL https://dspace.mit.edu/handle/1721.1/121709.

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

Luo, H., Svendsen, H.F., 1996. Theoretical model for drop and bubble breakup in turbulent dispersions. AIChE J. 42 (5), 1225–1233. http://dx.doi.org/10.1002/aic.690420505.

Magolan, B., Lubchenko, N., Baglietto, E., 2019. A quantitative and generalized assessment of bubble-induced turbulence models for gas-liquid systems. Chem. Eng. Sci. X 2, 100009. http://dx.doi.org/10.1016/j.cesx.2019.100009.

Manon, E., 2000. 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). École centrale de Paris.

Marfaing, O., Guingo, M., Laviéville, J., Mimouni, S., Baglietto, E., Lubchenko, N., Magolan, B., Sugrue, R., Nadiga, B.T., 2018. Comparison and uncertainty quantification of two-fluid models for bubbly flows with NEPTUNE_CFD and STAR-CCM+. Nucl. Eng. Des. 337, 1–16. http://dx.doi.org/10.1016/j.nucengdes.2018.05.028.

Menter, F.R., 1994. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 32 (8), 1598–1605. http://dx.doi.org/10.2514/3.12149.

Mérigoux, N., 2022. Multiphase Eulerian-Eulerian CFD supporting the nuclear safety demonstration. Nucl. Eng. Des. 397, 111914. http://dx.doi.org/10.1016/j.nucengdes.2022.111914.

Mimouni, S., Archambeau, F., Boucker, M., Laviéville, J., Morel, C., 2009. A second order turbulence model based on a Reynolds stress approach for two-phase flow - Part I: Adiabatic cases. Sci. Technol. Nucl. Installations 2009 (ID 792395).

Mimouni, S., Archambeau, F., Boucker, M., Laviéville, J., Morel, C., 2010a. A second order turbulence model based on a Reynolds stress approach for two-phase boiling flow. Part 1: Application to the ASU-annular channel case. Nucl. Eng. Des. 240 (9), 2233–2243.

Mimouni, S., Archambeau, F., Boucker, M., Lavieville, J., Morel, C., 2010b. A second order turbulence model based on a Reynolds stress approach for two-phase boiling flow and application to fuel assembly analysis. Nucl. Eng. Des. 240 (9), 2225–2232. http://dx.doi.org/10.1016/j.nucengdes.2009.11.020, Experiments and CFD Code Applications to Nuclear Reactor Safety (XCFD4NRS).

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

Mimouni, S., Benguigui, W., Lavieville, J., Merigoux, N., Guingo, M., Baudry, C., Marfaing, O., 2016b. New nucleation boiling model devoted to high pressure flows. In: Proceedings of the 9th International Conference on Multiphase Flow ICMF 2016. Firenze, Italy.

Mimouni, S., Boucker, M., Laviéville, J., Guelfi, A., Bestion, D., 2008. Modelling and computation of cavitation and boiling bubbly flows with the NEPTUNE_CFD code. Nucl. Eng. Des. 238 (3), 680–692. http://dx.doi.org/10.1016/j.nucengdes.2007.02.052.
Mimouni, S., Laviéville, J., Seiler, N., Ruyer, P., 2011. Combined evaluation of second order turbulence model and polydispersion model for two-phase boiling flow and application to fuel assembly analysis. Nucl. Eng. Des. 241 (11), 4523-4536.

Morel, C., 2015. Mathematical Modeling of Two-Phase Flow. Springer.

Morel, C., Ruyer, P., Seiler, N., Laviéville, J.M., 2010. Comparison of several models for multi-size bubbly flows on an adiabatic experiment. Int. J. Multiph. Flow 36 (1), 25-39.

Murallidharan, J.S., Prasad, B.V.S.S.S., Patnaik, B.S.V., Hewitt, G.F., Badalassi, V., 2016. CFD investigation and assessment of wall heat flux partitioning model for the prediction of high pressure subcooled flow boiling. Int. J. Heat Mass Transfer 103, 211-230. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.06.050.

Peltola, J., Pättikangas, T., 2012. Development and validation of a boiling model for OpenFOAM multiphase solver. In: CFD4NRS-4 Conference.

Peltola, J., Pättikangas, T., Bainbridge, W., Lehnigk, R., Schlegel, F., 2019. On development and validation of subcooled nucleate boiling models for OpenFOAM foundation release. In: 18th International Topical Meeting on Nuclear Reactor Thermal Hydraulics, NURETH 2019. Portland, Oregon, United States.

Prince, M.J., Blanch, H.W., 1990. Bubble coalescence and break-up in air-sparged bubble columns. AIChE journal 36 (10), 1485-1499. http://dx.doi.org/10.1002/aic.690481103.

Raj, M.N., Velusamy, K., 2017. CFD simulation of sodium boiling in heated pipe using RPI model. In: Proceedings of the 2nd World Congress on Momentum, Heat and Mass Transfer. MHMT'17, Barcelona, Spain.

Ranz, W.E., Marshall, W.R.J., 1952. Evaporation from drops. Chem. Eng. Prog. 48, 173-180.

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

Ruyer, P., Seiler, N., Weiss, M., Weiss, F.P., 2007. A bubble size distribution model for the simulation of bubbly flows.. In: Proceedings of the 6th International Conference on Multiphase Flows. Leipzig, Germany.

Sato, Y., Sekoguchi, K., 1975. Liquid velocity distribution in two-phase bubble flow. Int. J. Multiph. Flow 2 (1), 79-95. http://dx.doi.org/10.1016/0301-9322(75)90030-0.

Setoodeh, H., Ding, W., Lucas, D., Hampel, U., 2021a. Modelling and simulation of flow boiling with an Eulerian-Eulerian approach and integrated models for bubble dynamics and temperature-dependent heat partitioning. Int. J. Therm. Sci. 161, 106709. http://dx.doi.org/10.1016/j.ijthermalsci.2020.106709.

Setoodeh, H., Ding, W., Lucas, D., Hampel, U., 2021b. Prediction of bubble departure in forced convection boiling with a mechanistic model that considers dynamic contact angle and base expansion. Energies 12, 1950. http://dx.doi.org/10.3390/en12101950.

Setoodeh, H., Shabestary, A.M., Ding, W., Lucas, D., Hampel, U., 2022. CFD-modelling of boiling in a heated pipe including flow pattern transition. Appl. Therm. Eng. 204, 117962. http://dx.doi.org/10.1016/j.applthermaleng.2021.117962.

Shaver, D.R., Podowski, M.Z., 2015a. Modeling and validation of forced convection subcooled boiling. In: NURETH-16. Chicago, IL, USA.

Shaver, D.R., Podowski, M.Z., 2015b. Modeling of interfacial forces for bubbly flows in subcooled boiling conditions. Trans. Am. Nucl. Soc. 113 (10), 1368-1371.

Simonin, O., Viollet, P.L., 1990. Modelling of turbulent two-phase jets loaded with discrete particles. In: Phenomena in Multiphase Flows, vol. 1990, Hemisphere Publishing Corporation, pp. 259-269.

Software, S.D.I., 2023. Simcenter STAR-CCM+ user guide v. 2023.1.

Speziale, G.C., Sarkar, S., Gatski, T.B., 1991. Modelling the pressure strain correlation of turbulence: an invariant dynamical systems approach. J. Fluid Mech. 227, 245272. http://dx.doi.org/10.1017/S0022112091000101.

Sugrue, R.M., 2017. A robust momentum closure approach for multiphase computational fluid dynamics applications (Ph.D. thesis). Massachusetts Institute of Technology.

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

Tolubinsky, V.I., Kostanchuk, D.M., 1970. Vapour bubbles growth rate and heat transfer intensity at subcooled water boiling. In: Proceedings of 4th International Heat Transfer Conference. 23, Begel House Inc..

Tomiyama, A., 1998. Struggle with computational fluid dynamics. In: Third International Conference on Multiphase Flow. ICMF 98. Lyon, France.

Tomiyama, A., Kataoka, I., Zun, I., Sakaguchi, T., 1998. Drag coefficients of single bubbles under normal and micro gravity conditions. JSME Int. J. 41 (2), 472-479.

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

Tripathi, R.G., Buwa, V.V., 2014. Numerical simulations of gas-liquid boiling flows using OpenFOAM. In: IUTAM Symposium on Multiphase Flows with Phase Change: Challenges and Opportunities. 15, pp. 178-185. http://dx.doi.org/10.1016/j.piatam.2015.04.025.

Unal, H.C., 1976. Maximum bubble diameter, maximum bubble growth time and bubble growth rate during subcooled nucleate flow boiling of water up to 17.7 MN/m². Int. J. Heat Mass Transfer 19 (6), 643-649.

Vaesson, G.E.J., 1998. Coalescence and breakup modelling in a two-phase CFD code: Theoretical framework. (BRITE Report III-1).

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

Vyskocil, L., Macek, J., 2008. Boiling flow simulation in Neptune_CFD and Fluent codes. In: Workshop on Experiments and CFD Code Application To Nuclear Reactor Safety (XCFD4NRS). pp. 1-13, URL https://inis.iaea.org/search/search.aspx?orig_q=RN:44098671.

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

Wang, H., Xiong, J., Wang, J., 2023. Development and assessment of five-component wall boiling heat flux partitioning model. Int. J. Multiph. Flow 158, 104306.

Wellek, R.M., Agrawal, A.K., Skelland, A.H.P., 1966. Shape of liquid drops moving in liquid media. AIChE J. 12 (5), 854-862. http://dx.doi.org/10.1002/aic.690120506.

Wu, Q., Kim, S., Ishii, M., Beus, S.G., 1998. One-group interfacial area transport in vertical bubbly flow. Int. J. Heat Mass Transfer 41 (8-9), 1103-1112.

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

Yeoh, G.H., Vahaji, S., Cheung, S.C.P., Tu, J.Y., 2014. Modeling subcooled flow boiling in vertical channels at low pressures. part 2: Evaluation of mechanistic approach. Int. J. Heat Mass Transfer 75, 754-768. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.03.017.

Yoo, J., Estrada-Perez, C.E., Hassan, Y.A., 2016. Experimental study on bubble dynamics and wall heat transfer arising from a single nucleation site at subcooled flow boiling conditions part 1: Experimental methods and data quality verification. Int. J. Multiph. Flow 84, 315-324. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2016.04.018.

Yoon, S.-J., Agostinelli, G., Baglietto, E., 2017. Assessment of multiphase CFD with zero closure model for boiling water reactor fuel assemblies. In: Proc. of the 17th International Topical Meeting on Nuclear Reactor Thermal Hydraulics. (NURETH 17), Xi'an, China, p. 15.

Yun, B.J., Splawski, A., Lo, S., Song, C.-H., 2010. Advanced wall boiling model with wide range applicability for the subcooled boiling flow and its application into the CFD code. In: Transactions of the Korean Nuclear Society Spring Meeting. Pyeongchang, Korea.

Yun, B.-J., Splawski, A., Lo, S., Song, C.-H., 2012. Prediction of a subcooled boiling flow with advanced two-phase flow models. Nucl. Eng. Des. 253, 351-359. http://dx.doi.org/10.1016/j.nucengdes.2011.08.067, SI : CFD4NRS-3.

Ziegenhein, T., Lucas, D., 2017. Observations on bubble shapes in bubble columns under different flow conditions. Exp. Therm Fluid Sci. 85, 248-256. http://dx.doi.org/10.1016/j.expthermflusci.2017.03.009.

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