# Heat flux partition based on onset of significant void

Corentin Reiss, Antoine Gerschenfeld, Catherine Colin

► To cite this version:

Corentin Reiss, Antoine Gerschenfeld, Catherine Colin. Heat flux partition based on onset of significant void. International Journal of Multiphase Flow, 2024, 181, pp.104972. 10.1016/j.ijmultiphaseflow.2024.104972. hal-04837831v2

HAL Id: hal-04837831

https://ut3-toulouseinp.hal.science/hal-04837831v2

Submitted on 13 Dec 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
# Heat flux partition based on onset of significant void

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

## ARTICLE INFO

MSC:

80A20

76T10

Keywords:

Flow boiling

Heat flux partition

Onset of significant void

Multiphase CFD

## ABSTRACT

The thermal log-law $\Theta_+(y_+) = \beta + 2.12 \log(y_+)$ is valid in flow boiling with a value of $\beta$ that evolves as the flow develops. Using a multiphase flow cross-literature database, this constant is shown to be $\beta_{OSV} = -7$ at the point of onset of significant void (OSV). This means that at the OSV the liquid is at saturation temperature up to $y_+ \approx 30$. The OSV predictions using this model have a similar mean average error as the Saha and Zuber 1974 correlation for one less fitted constant for channel, pipe and annular flows for pressure from 1 to 147 bar and Peclet numbers from $3.5 \cdot 10^3$ to $4 \cdot 10^5$. This model is used to build a heat flux partitioning (HFP) inspired from system-scale codes (Lahey 1978). It predicts the distribution of the heat flux between the liquid phase and the evaporation term when the total heat flux is known. It does not give information on the total heat flux as a function of wall temperature and cannot be used to draw a boiling curve. In imposed flux conditions, this partition provides more coherent flux distribution between the evaporation and liquid terms than Kurul and Podowski (1990) based approaches and improves void fraction predictions in high-subcooling regions on the DEBORA database (Garnier et al. 2001) and on experiments by Bartolomei and Chanturiya (1967) and Bartolomei et al. (1982). When the wall temperature is imposed, it must be coupled with an empirical boiling total heat flux correlation to replace a traditional HFP. The prediction of the total heat flux is then as good as that of the correlation.

## 1. Introduction

Multiphase computational fluid dynamic models aim to simulate complex 3D flows without needing calibration for new geometries and flow conditions (Lahey et al., 2021). In that, they differ from 1D models that need to be adjusted with experimental data (Emonot et al., 2011; Berry et al., 2018). Most codes today use the two-fluid model (Ishii and Hibiki, 2006; Lahey et al., 2021). One such model is described in Appendix D. It is based on an averaging process on the phasic equations for mass, momentum and energy and their jump conditions. Once the equations are obtained, many interfacial transfer terms that require additional modeling appear.

For boiling-flow simulations, one of the most important closure laws is the wall heat transfer modeling in nucleate boiling. One must determine the total heat flux and its distribution between the phases thanks to heat flux partition (HFP) models. These use the wall, liquid, vapor and saturation temperatures and physical characteristics of the phases as inputs. They determine the heat flux that goes into each phase $q_l$ and $q_v$, and the evaporation heat flux $q_{l \to v}$ (Kurul and Podowski, 1990). The heat flux that enters the vapor phase $q_v$ is important only close to the critical heat flux (CHF) (Baglietto et al., 2019; Mimouni et al., 2016a). It will not be considered in this work. Furthermore, when

the liquid is at saturation temperature in the near-wall region, virtually all of the heat flux is evaporation: to the best of our knowledge, no liquid superheat has ever been measured in flow boiling (Garnier et al., 2001; Roy et al., 2002; Francois et al., 2021).

In practice, simulations are most often run with an imposed total heat flux $q_w$ at the wall (Favre, 2023). A Newton algorithm is then used to find the wall temperature so that $q_l(T_w) + q_{l \to v}(T_w) = q_v$, where $T_w$ is the wall temperature. The HFP is therefore mainly used to determine the heat flux *distribution*, and not the heat flux itself. Simulations can also be run with an imposed wall temperature, in which case the HFP predicts $q_w$.

**Heat flux partitioning in subcooled flow boiling.** The reference HFP for CFD codes was proposed by Kurul and Podowski (1990). It is a mechanistic model that includes three heat transfer mechanisms: one to the liquid phase that is based on single-phase flow, enabling a smooth transition as boiling picks off; one to the liquid phase that comes from the rewetting of the wall after bubble departure, called quenching; one evaporation term. The details of this model are given in Appendix C.2. In order to calculate these terms, many intermediate quantities are used like the nucleation site density, a bubble growth time scale and a bubble

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

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

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

Roman letters

<table><tr><td>A</td><td>area of the test section</td></tr><tr><td>C<sub>p</sub></td><td>heat capacity of the liquid</td></tr><tr><td>D<sub>h</sub> = 4A/P</td><td>hydraulic diameter of the test section</td></tr><tr><td>D<sub>he</sub> = 4A/P<sub>h</sub></td><td>heated diameter of the test section</td></tr><tr><td>G</td><td>mass flux</td></tr><tr><td>h<sub>k</sub></td><td>enthalpy of phase k</td></tr><tr><td>h<sub>ks</sub></td><td>enthalpy of phase k at saturation</td></tr><tr><td>h<sub>k,bulk</sub></td><td>bulk enthalpy of phase k</td></tr><tr><td>h<sub>OSV</sub></td><td>bulk liquid enthalpy at OSV</td></tr><tr><td>H<sub>l</sub></td><td>heat transfer coefficient toward the liquid phase</td></tr><tr><td>H<sub>l,SP</sub>(y<sub>+</sub>)</td><td>CFD-scale single-phase heat transfer coefficient</td></tr><tr><td>H<sub>l,OSV</sub>(y<sub>+</sub>)</td><td>CFD-scale heat transfer coefficient toward the liquid phase at OSV, enhanced by bubble agitation</td></tr><tr><td>P<sub>out</sub></td><td>outlet pressure of a test</td></tr><tr><td>P</td><td>perimeter of the test section</td></tr><tr><td>P<sub>h</sub></td><td>heated perimeter of the test section</td></tr><tr><td>q<sub>w</sub></td><td>total heat flux from the wall to the flow</td></tr><tr><td>q<sub>k</sub></td><td>heat flux transferred to phase k</td></tr><tr><td>q<sub>l&rarr;v</sub></td><td>evaporation heat flux</td></tr><tr><td>q<sub>l,OSV</sub></td><td>heat flux to liquid at OSV</td></tr><tr><td>q<sub>SP</sub></td><td>single-phase flow total heat transfer</td></tr><tr><td>q<sub>Boil</sub></td><td>boiling flow total heat transfer</td></tr><tr><td>r<sub>+</sub></td><td>dimensionless position in a tube: r<sub>+</sub> = 0 at the center and r<sub>+</sub> = 1 at the wall.</td></tr><tr><td>T<sub>l</sub></td><td>local liquid temperature</td></tr><tr><td>T<sub>s</sub></td><td>saturation temperature</td></tr><tr><td>T<sub>w</sub></td><td>wall temperature</td></tr><tr><td>T<sub>bulk</sub> = &lt;uT<sub>l</sub>&gt; / &lt;u&gt;</td><td>bulk liquid temperature</td></tr><tr><td>T<sub>*</sub> = q<sub>w</sub>/(ρ<sub>l</sub>C<sub>p</sub>u<sub>τ</sub>)</td><td>liquid temperature turbulent scale</td></tr><tr><td>u<sub>bulk</sub></td><td>bulk velocity of the fluid</td></tr><tr><td>u<sub>τ</sub></td><td>friction velocity</td></tr><tr><td>X = (h - h<sub>ls</sub>)/(h<sub>gs</sub> - h<sub>ls</sub>)</td><td>thermodynamic quality</td></tr><tr><td>X<sub>OSV</sub></td><td>thermodynamic quality at onset of significant void</td></tr><tr><td>y</td><td>distance to nearest wall</td></tr><tr><td>y<sub>+</sub> = yu<sub>τ</sub>/v<sub>l</sub></td><td>normalized distance to nearest wall</td></tr><tr><td>y<sub>+,s</sub></td><td>thickness of the liquid layer at T<sub>s</sub></td></tr><tr><td>y<sub>+,c</sub></td><td>critical thickness of the saturated liquid layer at which OSV is observed</td></tr><tr><td>y<sub>l</sub></td><td>size of the near-wall mesh element</td></tr><tr><td>y<sub>+,1</sub> = y<sub>l</sub>u<sub>τ</sub>/v<sub>l</sub></td><td>normalized size of the near-wall mesh element</td></tr></table>

Greek letters

<table><tr><td>α</td><td>void fraction</td></tr><tr><td>&lt;α&gt;</td><td>average void fraction on the test section</td></tr><tr><td>ΔT<sub>bulk</sub> = T<sub>s</sub> - T<sub>bulk</sub></td><td>liquid bulk subcooling</td></tr><tr><td>Θ<sup>w</sup><sub>+</sub>(y<sub>+</sub>) = (T<sub>w</sub> - T<sub>l</sub>)/T *</td><td>normalized liquid temperature with respect to wall temperature</td></tr><tr><td>Θ<sup>s</sup><sub>+</sub> = (T<sub>l</sub> - T<sub>s</sub>)/T<sub>*</sub></td><td>normalized liquid temperature with respect to saturation temperature</td></tr></table>

<table><tr><td>λ<sub>l</sub></td><td>liquid conductivity</td></tr><tr><td>ν<sub>l</sub></td><td>liquid kinematic viscosity</td></tr><tr><td>ρ<sub>k</sub></td><td>density of phase k</td></tr></table>

Subscripts

<table><tr><td>k</td><td>any phase</td></tr><tr><td>l</td><td>liquid phase</td></tr><tr><td>υ</td><td>vapor phase</td></tr><tr><td>s</td><td>saturation</td></tr><tr><td>w</td><td>wall</td></tr></table>

Dimensionless numbers

<table><tr><td>Nu = <sup>q<sub>w</sub>D<sub>h</sub></sup>&frasl;<sub>λ<sub>l</sub>(T<sub>s</sub> - T<sub>bulk</sub>)</sub></td><td>Nusselt number</td></tr><tr><td>Pe = <sup>λ<sub>l</sub></sup>&frasl;<sub>D<sub>h</sub>u<sub>bulk</sub>ρ<sub>l</sub>C<sub>p</sub></sub></td><td>Peclet number</td></tr><tr><td>Pr = <sup>ν<sub>l</sub>ρ<sub>l</sub>C<sub>p</sub></sup>&frasl;<sub>λ<sub>l</sub></sub></td><td>Prandtl number</td></tr><tr><td>Re = <sup>D<sub>h</sub>u<sub>bulk</sub></sup>&frasl;<sub>ν<sub>l</sub></sub></td><td>Reynolds number</td></tr><tr><td>St = <sup>q<sub>w</sub></sup>&frasl;<sub>GC<sub>p</sub>ΔT<sub>bulk</sub></sub></td><td>Stanton number</td></tr></table>

detachment diameter. Few measurements of these terms exist in the literature, particularly at high pressures and mass fluxes. Furthermore, the authors make the hypothesis that all bubbles depart at the same size, that all nucleation points behave in the same way and do not interact with each other and that they do not affect single-phase convection in locations where there is not bubble nucleation. Finally, no direct measures exist of the quenching heat flux, which is in this model a theoretical form based on conduction in the liquid phase, with no turbulent contribution. Again, the practical use of all of this modeling is mainly to determine which proportion of the heat flux enters the liquid phase through quenching and convection, and which proportion leads to evaporation.

More recent approaches, like those of Basu et al. (2005), Kommajosyula (2020) and Favre (2023), are also mechanistic. As time passes, the refinement of these models tends to increase by including more small-scale mechanisms, like bubble coalescence at the wall, interactions between nucleation sites, or the contribution of bubble sliding on the wall to the transfer to the liquid phase. This complexity makes the models very difficult to read and interpret without plotting the output heat fluxes as a function of the wall temperature or liquid temperature. It also makes code-to-code comparisons challenging and leads to long HFP calculation times. Though the number of closure terms increases significantly, some key mechanisms are missing, like the interaction between the bubble layer and detached bubbles which should influence the quenching flux.

To feed these models, experiments that allow a fine tracking of the bubble nucleation are being developed (Richenderfer et al., 2018; Tecchio, 2022). However, it is difficult to obtain high-pressure data. As far as we know, precise measurements were done at a maximum of 40 bar (Kossolapov, 2021), much lower than the 155 bar found in nuclear power plants (Todreas and Kazimi, 2021). An interesting takeaway from the work of Kossolapov (2021) is the huge variability of nucleation frequency between sites: for similar conditions, some are 1000 times greater than others. Bubble growth time and departure diameter also present big variations. This means that a rigorous mechanistic approach should also take into account distributions. To the best of our knowledge such a model has not been proposed for now.

Heat flux partition in system and component scale codes. These codes are used in the nuclear industry to simulate reactor cores, steam generators and primary and secondary circuits during steady-state operations and transients. In them, there is only one cell on the width of a channel or subchannel. They include, among others, TRACE (NRC,
2010), RELAP (Berry et al., 2018), Cathare (Emonot et al., 2011) and CTF (Salko Jr. et al., 2023). Many terms like the pressure drop or interfacial friction are calibrated using experiments specific to the geometry of the simulated section. In CTF and TRACE, the methodology from Lahey (1978) described hereafter can be selected by the user to determine the HFP in subcooled boiling. The first step consists in determining the total heat flux $q_w$ from the wall: either it is imposed, either it is calculated using a correlation that requires the wall and bulk liquid temperatures like that of Thom et al. (1965).

The bulk liquid enthalpy at the point of onset of significant void (OSV) $h_{OSV}$ is calculated in the second step. The OSV is defined as the point in a given flow where a noticeable increase in the void fraction of a flow takes place. $h_{OSV}$ is most often calculated with the Saha and Zuber (1974) correlation. More details on the definition of OSV and different approaches used in the literature can be found in Section 2.1.

In the Lahey (1978) methodology, the heat fluxes directed towards the liquid phase and the evaporation are then:

$$ \left\{ \begin{array}{ll} \text{If } h_{l,bulk} < h_{OSV} & \begin{cases} q_l = q_w \\ q_{l \to v} = 0 \end{cases} \\ \text{If } h_{l,bulk} > h_{OSV} & \begin{cases} q_l = q_w \frac{h_{ls} - h_{l,bulk}}{h_{ls} - h_{OSV}} \\ q_{l \to v} = q_w - q_l \end{cases} \end{array} \right. \quad (1) $$

Where $h_{ls}$ is the saturation enthalpy of the liquid phase.

Before OSV, all of the heat flux goes into the liquid. If the fluid is at saturation temperature, all of the energy is used for evaporation. There is a linear interpolation in between the two: as the fluid approaches saturation temperature, a larger part of the heat flux creates vapor.

This approaches has a few advantages compared with mechanistic approaches. It reduces the uncertainty on the total heat flux as a correlation adapted to the situation can be used. The known heat flux distributions before the OSV (single-phase only) and after the liquid is saturated (evaporation only) are respected, which cannot be guaranteed using a HFP. The behavior of the correlation is predictable: there are no hidden variables to calibrate and they are linked with a small number of simple equations. From a numerical point of view, no Newton algorithm is required to obtain the partition for imposed-heat flux conditions. This makes the calculation fast and the coding easy.

**Aims of this work.** We apply the 2-step system-scale HFP approach described previously to CFD codes to benefit from their advantages on physics and numerical standpoints. Numerous total heat flux correlations are already available in the literature (Jens and Lottes, 1951; Thom et al., 1965; Frost and Dzakowic, 1967), but they do not give the heat flux partition. Therefore, we strive to create a CFD-scale OSV criterion valid at high-Reynolds number to determine the heat flux distribution between the liquid phase and evaporation. We carry out the following steps:

* Review OSV prediction approaches from the literature (Section 2.1)

* Create a database of thermodynamic quality at OSV for different geometries and flow conditions (Section 2.2) and discuss preliminary results (Section 2.3)

* Using the structure of the temperature distribution in a boiling flow (Sections 3.1 &3.2), transform this system-scale data to local, CFD-scale data (Section 3.3)

* Create a local correlation for OSV (Section 3.4)

* Compare our results to the Saha and Zuber (1974) correlation, chosen as reference (Section 3.5)

* Use our correlation to build our own heat flux partition (Section 4)

* Compare the fraction of heat flux going to each phase for our HFP and Kurul and Podowski (1990)-based models (Section 4.4)

* Compare boiling 2-fluid simulations using our HFP and the original Kurul and Podowski (1990) model on the DEBORA (Cubizolles, 1996, Section 4.6) and Bartolomei (Bartolomei and Chanturiya, 1967; Bartolomei et al., 1982, Section 4.7) databases

## 2. Construction and analysis of an OSV database

### 2.1. Current approaches to OSV prediction

The onset of significant void (OSV), also called net vapor generation (NVG), is defined as the point where a noticeable increase in the void fraction in a boiling flow takes place. On a developing axial flow, this is a physical location. However, in general it is defined using the thermodynamic quality at this point, $X_{OSV}$. Cai et al. (2021) recently performed a literature review on the subject and compared different OSV correlations to data from the literature. Lee and Bankoff (1998) had previously done a similar exercise.

The methods most commonly used today to predict the point of OSV are global empirical approaches based on dimensionless numbers that do not go into the details of the physical mechanisms at play. Different mechanistic models have also been proposed, investigating bubble dynamics in the near-wall region.

**Empirical approaches.** The most commonly used OSV correlation is from Saha and Zuber (1974). It is a full-channel, empirical formula. Two regions are defined using the Peclet number $Pe = \frac{D_h u_{bulk} \rho_l C_p}{\lambda_l}$, where $D_h$ is the hydraulic diameter of the test section, $u_{bulk}$ the bulk velocity, $\rho_l$ the liquid density, $C_p$ the liquid heat capacity and $\lambda_l$ the liquid conductivity. For $Pe < 7 \cdot 10^4$, they propose that at the OSV vapor condensation and evaporation at the wall are at equilibrium. Evaporation is proportional to the heat flux $q_w$, and condensation to the local subcooling and conductivity $\lambda_l(T_s - T_{bulk})$, where $T_s$ is the saturation temperature and $T_{bulk}$ the bulk liquid temperature. This leads to a constant Nusselt number $Nu = \frac{q_w D_h}{\lambda_l(T_s - T_{bulk})}$. For high-$Pe$, they propose that bubble detachment is hydro-dynamically controlled. They propose that the Stanton number $St = \frac{q_w}{GC_p T_{bulk}}$ is constant at the OSV, where G is the mass flux. Fitting this law on data from 8 different experimental sources, they obtain that the flow is over the point of OSV if:

$$ \left\{ \begin{array}{ll} Nu = \frac{q_w D_h}{\lambda_l(T_s - T_{bulk})} \ge 455 & \text{if } Pe \le 7 \cdot 10^4 \\ St = \frac{q_w}{GC_p(T_s - T_{bulk})} \ge 0.0065 & \text{if } Pe > 7 \cdot 10^4 \end{array} \right. \quad (2) $$

In recent years, several refinements of this correlation have been proposed (Ha et al., 2020; Lee and Jeong, 2022). The formula and the results obtained are close to those of the original correlation. Both Cai et al. (2021) and Lee and Bankoff (1998) find that the Saha and Zuber expression and the ones that are similar outperform the others encountered in the literature. On their database, Cai et al. find a mean average error of ~20%, vs ~35% for the Levy (1967) correlation. We will consider the Saha and Zuber (1974) correlation, which will be noted SZ in the rest of this paper, as the reference formulation.

**Mechanistic approaches.** The Levy (1967) criterion is the first mechanistic approach on record. He determines a detachment radius for a bubble using an axial force balance between buoyancy, drag and the surface tension force holding back the bubble. He then assumes the point of OSV is reached when the thickness of the liquid layer that is at $T_s$ is greater than the detachment radius. The obtained correlation is local, and contains two fitted parameters on a database containing data from 5 experimental sources. The temperature and velocity of the liquid layer used in the calculation follow single-phase developed distributions, not accounting for bubble presence. Furthermore, in high-pressure flows bubbles slide along the wall as soon as they nucleate (Kossolapov, 2021), therefore the surface tension force holding back bubbles cannot be the key mechanism to determine the point of OSV.

Dix (1971) proposed a model based on radial, not axial, bubble movement. His experiments lead him to believe that the bubble ejection is hydro-dynamically controlled. He assumes the bubble radius at the wall is proportional to the heat transfer coefficient, and that at a critical
bubble diameter depending on the Reynolds number can be defined. A single fitting parameter is used, but only on one experimental data set. This limits the validity of the correlation to a 9.6 mm inner diameter, 18.6 mm outer diameter annular tube, with $1 \cdot 10^4 \le Pe \le 2 \cdot 10^4$.

Anne and Beattie (1996) argue that the OSV occurs when bubbles can survive in the center of the tube. They determine the amplitude of turbulent fluctuations as a function of the total heat flux. They consider that the OSV is reached when turbulent fluctuations enable bubbles to reach a large part of the tube in a liquid pocket that is at saturation temperature, i.e. when these fluctuations become larger than the bulk subcooling. This yields that at the OSV, $St = 0.088Pe^{-1/2}$, without the use of any fitting coefficient. This expression is close to the SZ correlation for $Pe \sim 3 \cdot 10^5$, but does not match the data for lower $Pe$.

Recently, Nguyen and Okawa (2024) observed OSV occurs when an important bubble coalescence takes place at the wall. They do not propose a quantitative model, but they believe that coalescence in the bubble layer is the key mechanism for OSV.

## 2.2. Building an OSV database

The OSV is dependent on the structure of the flow. In particular, it will be different if the single-phase thermal boundary layer before OSV is developed or not. In the experiments that we select to determine the OSV, sources are restricted to the cases where the heater is long enough for the single-phase thermal boundary layer to be developed.

One of the main difficulties in building an OSV database is the precise definition of the OSV. Most authors use plots of the void fraction $\alpha$ as a function of the thermodynamic quality $X$. The shape of the curve is often similar to a hyperbolic tangent (see Fig. 1). Depending on the authors, the OSV can be defined as the moment where the curve starts to take off, or after the void is already significant, leading to smaller subcooling at OSV. Many also extrapolated tanh-like functions from a small number of points and used the extrapolated plots to determine the point of OSV.

To harmonize the definition of the OSV, we select ourselves the point of OSV from $(X, < \alpha >)$ plots found in different sources. $X$ is the thermodynamic quality of the flow and $< \alpha >$ the average void fraction in the test section. We restrict ourselves to runs where no curve fitting and extrapolation is necessary to see the inflection, so there must be measures with $< \alpha >$ close to 0. We also only select points where the inflection in void fraction occurs for $< \alpha > < 5\%$, and where we estimate that the uncertainty on the OSV is smaller than 10%. As the OSV determination is done by hand and the OSV definition is not extremely precise, we believe it is difficult to increase this accuracy. An example of the point of OSV that we determine can be found in the red square point in Fig. 1. Some plots that we eliminate can be found there as well. The only exception is the data from Edelman and Elias (1981), that we choose to include even without access to the original $(X, < \alpha >)$ plots as we found no other database in the literature at so low $Pe$ numbers $(2 \cdot 10^3 \le Pe \le 4 \cdot 10^3)$.

In the end, we obtain 155 OSV data points in a wide range of geometries, pressures and Peclet numbers. This data can be written in

the form $X_{OSV}$ (Test section, $P, G, q_w$). Table 1 contains the list of all of the sources that we used. This database and the treatments we apply to it in the rest of this paper can be found at https://github.com/CoReiss/CFD_OSV.

We were not able to obtain the original data from other authors who compared different models, but in Appendix A we reproduce the plots presented by Saha and Zuber (1974) (Fig. A.16), Cai et al. (2021) (Fig. A.17) and Lee and Jeong (2022) (Fig. A.18). Our plots are close to those of these authors, which gives us a good level of confidence in our data collection.

## 2.3. System-scale analysis

Using the classical approach for OSV pioneered by Saha and Zuber (1974), we plot the Stanton number at OSV as a function of the Peclet number calculated using the hydraulic and heated diameter $D_h$ and $D_{he}$ for all of the runs of our database (Fig. 2).

$D_h$ and $D_{he}$ are different for annular geometries (Rouhani and Zeitoun) and channels heated on one side (Staub_Ch). Using $D_{he}$, the Saha and Zuber correlation remains valid for $Pe > 7 \cdot 10^4$, as it is a constant Stanton region and $St$ is independent of the diameter. This is no longer the case for $Pe < 7 \cdot 10^4$: data points are shifted to the right as $D_{he} \ge D_h$. While it seems logical that the characteristic length scale for high-$Pe$ region is $D_h$, we find it counter-intuitive in the low-$Pe$ region. According to SZ this region is dominated by thermal conduction and not turbulent effects, which is why the Nusselt number is constant. If this were the case, we expect $D_{he}$ to be the relevant length scale.

The transition from their so-called thermally controlled to detachment controlled regions is at $Pe = 7 \cdot 10^4$, though the temperature profile should be turbulent from $Pe \ge 1 \cdot 10^3$. Even if the transition is thermally controlled, the characteristic thermal diffusivity for $1 \cdot 10^3 \le Pe \le 7 \cdot 10^4$ is the turbulent viscosity, not the molecular diffusivity $\lambda_l/(\rho_l C_p)$, so a constant $Nu$ here is unexpected.

The characteristic temperature used in Eq. (2) is the bulk liquid temperature. One would expect the near-wall liquid temperature to play a great role in the OSV, as it is a near-wall phenomenon.

Finally, the density ratio, i.e. pressure, seems to have no impact on the outcome, though it is extremely important in the near-wall dynamics. At high pressures, bubbles no longer stick to the surface but start sliding as they nucleate (Kossolapov, 2021). Bubble diameters are very pressure-dependent (Kossolapov et al., 2024). The OSV is therefore not directly related to precise near-wall bubble dynamics.

# 3. CFD-scale approach

## 3.1. Literature review on the liquid temperature profile in boiling flow

Kader and Yaglom (1972) showed that in the turbulent boundary layer of a heated wall in single-phase flow, the dimensionless temperature profile writes:

$$ \Theta_{+}^{W}(y_{+}) = 2.12 \log(y_{+}) + \beta \quad (3) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Examples of $(X, < \alpha >)$ plots found in the literature that we use to determine the point of OSV. The red square marker in the first plot represents the identified point of OSV. Series without markers are disqualified for different reasons. Each run condition is given in the legend of the figure. Water is used in every one. Rouhani_A&B: Data from Rouhani (1966b). Rouhani_B is disqualified as we do not have enough points at low void fraction to see the inflection in $(X, < \alpha >)$. Sekoguchi_A: Data from Sekoguchi et al. (1980). This run is not considered as it does not have enough datapoints. Zeitoun_A&B: Data from Zeitoun (1994). Data in run Zeitoun_A is too irregular. In Zeitoun_B, $< \alpha > \sim 10\%$ at the inflection, proposed at $X = -0.022$ by the author. We find this too high for an OSV criterion.</figcaption></figure>
<caption>**Table 1**</caption>

Bibliographic sources used to calibrate our model. All runs are water, except the Staub_TuF runs which are R22 freon. Label column is used in subsequent figures. Geometry column gives the geometry of the test section (Pipe: circular cross-section pipe; Channel: rectangular cross-section, of which we give the short length; 1 W (2 W) means only 1 (2) wall(s) is (are) heated; Annulus: Annular cross-section, of which we give the inner diameter (ID) and the outer diameter (OD)). Press is outlet pressure. Nb Pts column contains the number of data points we took from the source. Total number of data points: 155. Used in column references other papers that used these data to calibrate or compare models (L: (Levy, 1967); SZ: (Saha and Zuber, 1974); LB: (Lee and Bankoff, 1998); C: (Cai et al., 2021); LJ: (Lee and Jeong, 2022)).

<table><thead><tr><th>Label</th><th>Reference</th><th>Geometry</th><th>Press (bar)</th><th>Peclet</th><th>Nb Pts</th><th>Used in</th></tr></thead><tbody><tr><td>Egen</td><td>Egen et al. (1957)</td><td>2.6 mm<br>2 W Channel</td><td>138</td><td>3 · 10<sup>4</sup> -<br>6 · 10<sup>4</sup></td><td>7</td><td>L, SZ,<br>LJ</td></tr><tr><td>Ferrell</td><td>Ferrell (1964)</td><td>11.6 mm Tube</td><td>4-16</td><td>4 · 10<sup>4</sup> -<br>8 · 10<sup>4</sup></td><td>11</td><td>L, LB,<br>C, LJ</td></tr><tr><td>Rouhani</td><td>Rouhani (1966a,b)</td><td>12 mm ID<br>25 mm OD<br>Annulus</td><td>10-40</td><td>1 · 10<sup>4</sup> -<br>10 · 10<sup>4</sup></td><td>6</td><td>L, SZ<br>LB, LJ</td></tr><tr><td>Bartolomei_1</td><td>Bartolomei and Chanturiya (1967)</td><td>15.4 & 24 mm Tube</td><td>15-45</td><td>9 · 10<sup>4</sup> -<br>15 · 10<sup>4</sup></td><td>10</td><td>SZ, C,<br>LJ</td></tr><tr><td>Bartolomei_2</td><td>Bartolomei et al. (1982)</td><td>12 mm Tube</td><td>30-147</td><td>9 · 10<sup>4</sup> -<br>30 · 10<sup>4</sup></td><td>16</td><td>C, LJ</td></tr><tr><td>Staub_Ch</td><td>Staub et al. (1969)</td><td>6.3 mm<br>1 W Channel</td><td>1.3-3</td><td>2 · 10<sup>4</sup> -<br>20 · 10<sup>4</sup></td><td>17</td><td>SZ, LB</td></tr><tr><td>Staub_TuF</td><td>Staub et al. (1969)</td><td>10.2 mm Tube</td><td>11-68</td><td>3 · 10<sup>4</sup> -<br>35 · 10<sup>4</sup></td><td>6</td><td>SZ, LB,<br>LJ</td></tr><tr><td>Staub_TuW</td><td>Staub et al. (1969)</td><td>10.2 mm Tube</td><td>11-68</td><td>7 · 10<sup>4</sup> -<br>11 · 10<sup>4</sup></td><td>2</td><td>SZ, LB,<br>LJ</td></tr><tr><td>Martin</td><td>Martin (1969, 1972)</td><td>2 & 2.8 mm<br>2 W Channel</td><td>78-138</td><td>2 · 10<sup>4</sup> -<br>16 · 10<sup>4</sup></td><td>22</td><td>SZ</td></tr><tr><td>Sabotinov</td><td>Sabotinov (1974)<sup>a</sup></td><td>11.7 mm Tube</td><td>68-108</td><td>8 · 10<sup>4</sup> -<br>20 · 10<sup>4</sup></td><td>8</td><td></td></tr><tr><td>Sekoguchi</td><td>Sekoguchi et al. (1980)</td><td>11, 13.6 &<br>15.8 mm Tube</td><td>1.3-16</td><td>3 · 10<sup>4</sup> -<br>16 · 10<sup>4</sup></td><td>16</td><td>LB, LJ</td></tr><tr><td>Edelman</td><td>Edelman and Elias (1981)</td><td>11.3 mm Tube</td><td>1</td><td>2 · 10<sup>3</sup> -<br>13 · 10<sup>3</sup></td><td>16</td><td>LB, C,<br>LJ</td></tr><tr><td>Labuntsov</td><td>Labuntsov et al. (1984)</td><td>12.1 & 34 mm Tube</td><td>5-70</td><td>8 · 10<sup>4</sup> -<br>40 · 10<sup>4</sup></td><td>4</td><td>C</td></tr><tr><td>Zeitoun</td><td>Zeitoun (1994)</td><td>12.7 mm ID<br>25.4 mm OD<br>Annulus</td><td>1.1-1.7</td><td>1 · 10<sup>4</sup> -<br>4 · 10<sup>4</sup></td><td>14</td><td>LJ</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** How our data fits the OSV the Saha and Zuber (1974) correlation. **Left:** Stanton number at OSV as a function of the Peclet number defined with the hydraulic diameter $D_h$. **Right:** Stanton number at OSV as a function of the Peclet number defined with the heated diameter $D_{he}$. The Saha and Zuber correlation is valid for the hydraulic diameter, not the heated diameter. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<footnote>a From Kolev (1985).</footnote>
Where $y_+ = yu_τ/v_l$, y is the distance to the nearest wall, $u_τ$ the friction velocity, $v_l$ the liquid kinematic viscosity, $\Theta_+^w(y_+) = (T_w - T_l(y))/T_*$ the dimensionless temperature, $T_l(y)$ the local liquid temperature, $T_* = q_w/(\rho_l C_p u_τ)$, and $\beta$ is a function of the Prandtl number. In particular, for low Prandtl numbers, $\beta$ can be negative.

In boiling flow, this log-law holds, as seen on various experiments. Roy et al. (2002) studied R113 in an annular channel, and found that on their data, where the local void fraction at the wall went up to 35% and bubbles had reached 40% of annulus width, $\Theta_+^w(y_+) = 1.95 \log(y_+) + 6.5$. Kledy (2018) performed measurements on R112 in a tube and for the 3 runs studied he finds $\Theta_+^w(y_+) = 2 \log(y_+) + 2.4$, $\Theta_+^w(y_+) = 1.9 \log(y_+) - 1.75$ and $\Theta_+^w(y_+) = 2.2 \log(y_+) - 8.05$. In these configurations, the void fraction was up to 40% at the wall and bubbles reached the center of the pipe. The ~2 prefactor of the log-law stayed valid for the temperature boundary layer, even though for the velocity field it changed dramatically, going from 2.5 to 7 between these 3 runs. This indicates that a log law of the form $\Theta_+^w(y_+) = 2.12 \log(y_+) + \beta$ remains valid for fully developed boiling flows, with $\beta$ depending on the flow regime, even though the velocity boundary layer changes form. We choose to keep the 2.12 pre-factor for consistency with the single-phase formulation.

In the OSV database that we have constructed, we do not have systematic measures of $T_w$ to define $\Theta_+^w(y_+)$. Furthermore, the saturation temperature is a clear reference temperature in flow boiling. Therefore, we chose to work in the remainder of this paper with $\Theta_+^s(y_+) = (T_s - T_l(y))/T_*$. This amounts to shifting the log-law by a constant, i.e. changing $\beta$, and the logarithmic profile should still be valid:

$$ \Theta_+^s(y_+) = 2.12 \log(y_+) + \beta \quad (4) $$

## 3.2. Analysis of the temperature profiles in the DEBORA experiment

To verify this, we use temperature measures from the DEBORA database (Garnier et al., 2001; Cubizolles, 1996). They consists in a series of experiments conducted on a 19.2 mm diameter 3.5m-long R12-filled tube with flow measurements at the outlet, where the flow is developed. Sets of (Mass flux, Pressure, Heat flux) conditions were selected. For each set, ~10 different entrance temperatures were used, thus changing the quality at the outlet. In this paper, all figures concerning data from the DEBORA experiment will use data from different entrance temperatures but same (Mass flux, Pressure, Heat flux) combinations. Each set of conditions is labeled G[ng]P[np]W[nw], where ng is the average mass flux (10³ kg/(m²s)), np the average pressure (bar) and nw the average total heat flux in the experiment (kW). In some test series the void fraction profiles were measured. These are used for comparison with our simulations in Section 4.6. In others the temperature profile at the end of the heated section were measured for single-phase and boiling flows. These are used in this section. To verify Eq. (4), the dimensionless temperature profiles are plotted as a function of the dimensionless distance to the wall in Fig. 3. $u_τ$ is calculated using the McAdams et al. (1942) correlation (see Eq. (9)).

In the five coolest runs of the G2P24W16 test (top row), the wall temperature is smaller than the saturation temperature. We have not yet reached the onset of nucleate boiling, and the corresponding plots of $\Theta_+^w(y_+)$ superimpose on the Kader and Yaglom (1972) single-phase profile. After the onset of nucleate boiling (ONB) is reached, when $T_{in} \ge 35.3C$, $\Theta_+^w$ and $\Theta_+^s$ shift towards lower values. $\beta$ decreases down to -15, becoming negative, until all of the liquid is at saturation temperature. The difference between the lower bound of $\Theta_+^w$ and $\Theta_+^s$ corresponds to the difference between wall and saturation temperatures. The G2P14W16 test (bottom row) behaves in a similar way, though there are no single-phase data points before ONB so no plateau of $\Theta_+^w$ is observed. For all of the plots in the figure, the log-law in Eq. (4) fitted for data points where $\Theta_+^s > 2$ remains valid, though $\Theta_+^s$ departs from the log-law when the flow is saturated.

Even for the high quality plots, we always have $\Theta_+^s > 0$. This is coherent with the fact that, to the best of our knowledge, no liquid temperature in boiling flow was ever measured above saturation temperature (Roy et al., 2002; Francois et al., 2021). We do not know the precise location of the OSV here, but Eq. (4) remains valid for $X_{out} \sim 0$ where we are in boiling flow. Therefore, the log-law is valid before ONB, between ONB and OSV and after OSV. The following equation is therefore a satisfying approximation for the dimensionless temperature on the complete data:

$$ \Theta_+^s(y_+) = \max(0, 2.12 \log(y_+) + \beta) \quad (5) $$

## 3.3. Transforming average-scale data to CFD-scale data

One way to interpret temperature profile in Eq. (5) is that $T_l(y) = T_s$ while $2.12 \log(y_+) + \beta < 0$, i.e. while:

$$ y_+ \le y_{+,c} = \exp \left( -\frac{\beta}{2.12} \right) \quad (6) $$

Our goal in this section is to identify a critical value $y_{+,c}$ so that the OSV occurs when the thickness of the saturated layer $y_{+,s}$ reaches $y_{+,c}$. An interpretation of such a criterion is presented in Fig. 4.

This approach is similar to the critical heat flux model developed by Nop et al. (2021), where the authors show that in highly subcooled atmospheric-pressure flows the boiling crisis is reached when a thick enough fluid layer reaches saturation temperature.

In order to work more easily with the temperature profiles, in practice we will look for:

$$ \beta_{OSV} = -2.12 \log(y_{+,c}) \quad (7) $$

So that at OSV the liquid temperature profile is:

$$ \Theta_+^s(y_+) = \max(0, 2.12 \log(y_+) + \beta_{OSV}) \quad (8) $$

We transform our database from the form $X_{OSV}$ (Test section, P, G, $q_w$) to the form $\beta_{OSV}$ (Test section, P, G, $q_w$). The first step is to calculate $u_τ$ for all of our data points. For tubes, the McAdams et al. (1942) correlation was used for high flow Reynolds number $Re = \frac{D_h u_{bulk}}{v_l}$. The Blasius (1913) formulation was used for low $Re$:

$$ \left\{ \begin{array}{ll} u_τ = u_{bulk}(0.184 Re^{-0.2}/8)^{1/2} & \text{if } Re > 3 \cdot 10^4 \\ u_τ = u_{bulk}(0.316 Re^{-0.25}/8)^{1/2} & \text{if } Re \le 3 \cdot 10^4 \end{array} \right. \quad (9) $$

For channels, we determine $u_τ$ so that the bulk velocity is correct using the log-law hypothesis along the width of the channel. For annular geometries, we run simulations with the $k-\omega$ turbulence model using the TrioCFD code (Angeli et al., 2015). The impact of an error on $u_τ$ on the prediction of the OSV is discussed in Section 3.4.

For given values of $\beta$, $u_τ$ and the physical properties of the fluids, we calculate the bulk liquid temperature. It is calculated by averaging the discharge over the flow's cross-section:

$$ T_{\text{bulk}} = \frac{\langle u_l(y) \rangle T_l(y) >}{\langle u_l(y) \rangle >} \quad (10) $$

For the local temperature field, we use Eq. (5). For tube and two-wall heated channels, the log-law is used for the liquid temperature across the whole width. For single-heated channels and annular flows, we take the temperature in the middle of the channel for the non-heated half.

For all flows, the local velocity used is the Reichardt (1951) single-phase adaptive wall law:

$$ \begin{aligned} u_l(y_+) &= u_\tau u_+(y_+) \\ u_+(y_+) &= \frac{1}{\kappa} \log(1 + 0.4y_+) + 7.8 \left( 1 - \exp\left(-\frac{y_+}{11}\right) - \frac{y_+}{11} \exp\left(-\frac{y_+}{3}\right) \right) \end{aligned} \quad (11) $$

To find $\beta_{OSV}$, the value of $\beta$ at OSV, for each point in our database, we perform a dichotomy on the value of $\beta$ so that the quality calculated for $\beta_{OSV}$ is $X_{OSV}$:

$$ \frac{h_l(T_{\text{bulk}}(\beta_{\text{OSV}})) - h_{ls}}{h_{gs} - h_{ls}} = X_{\text{OSV}} \text{ (Test section, P, G, } q_w). \quad (12) $$
<figure><img src="image_4.png" /><figcaption>Fig. 3. Dimensionless temperature profiles measured in the DEBORA experiment (Garnier et al., 2001; Cubizolles, 1996). A thermocouple is moved along the width of a 19.2 mm-diameter 3.5 m-long tube filled with refrigerant R12. Left column: Entrance temperature in the test section, and β fitted on Θ<sub>+</sub><sup>+</sup> on the corresponding run using Eq. (4). Central column: Dimensionless temperature difference between the liquid and the wall Θ<sub>+</sub><sup>-</sup>. There are two lines for each condition: one for each side of the pipe. The difference between both lines is an indicator of experimental error. Right column: Dimensionless temperature difference between the liquid and saturation temperature Θ<sub>+</sub><sup>+</sup> (full line) and fit on β using Eq. (4) (dashed line). Color scale: thermodynamic quality at the measuring point. Top: results for test number G2P26W16, conducted at G = 2 · 10<sup>3</sup> kg/(m<sup>2</sup>s), P<sub>out</sub> = 26.2 bar, q<sub>w</sub> = 73.9 kW/m<sup>2</sup> for various inlet temperatures. Bottom: results for test number G2P14W16, conducted at G = 2 · 10<sup>3</sup> kg/(m<sup>2</sup>s), P<sub>out</sub> = 14.6 bar, q<sub>w</sub> = 73.9 kW/m<sup>2</sup> for various inlet temperatures. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Physical interpretation of Eq. (6). After onset of nucleate boiling but before OSV, vapor is produced at the wall but in limited quantity and stays at small y<sub>+</sub>. At OSV, T<sub>i</sub>(y<sub>+</sub> < y<sub>+</sub><sub>c</sub>) ≈ T<sub>+</sub> and void can enter the flow. After OSV, vapor leaves the near-wall layer and the temperature can reach T<sub>+</sub> for higher y<sub>+</sub>.</figcaption></figure>

3.4. A simple CFD-scale correlation

As our applications are high-Pe flows in nuclear reactors, we con-
centrate on data where Pe > 5 · 10<sup>4</sup>. This limit was chosen as it is
slightly lower than the Pe = 7 · 10<sup>4</sup> transition in the Saha and Zuber
correlation. It should therefore encompass all turbulent-driven points
from their correlation. Furthermore, we select only geometries where
all walls are heated for the calibration step, as they have symmetrical
temperature fields. We keep tube and two-heated channel data. Annular
and single-heated channel data are used for the final model validation.

We strive to correlate β<sub>OSV</sub> with the local flow conditions and key
dimensionless numbers, i.e. u<sub>τ</sub>, Pr and ρ<sub>l</sub>/ρ<sub>v</sub>. As our goal is to determine
a CFD-scale correlation, we cannot use the Peclet or Reynolds numbers
for example. Fig. 5 shows β<sub>OSV</sub> as a function of u<sub>τ</sub>, with the color scale
representing the density ratio. The Prandtl number was left out as we
have only one source with another fluid than water: its variation is very
small across the database. No correlation is apparent in this plot.

<figure><img src="image_6.png" /><figcaption>Fig. 5. β at OSV vs friction velocity calculated for runs where Pe > 5 · 10⁴ and all wall of the test section are heated. The color scale represents the density ratio. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

In order to determine the optimal value of β<sub>OSV</sub>, rather than taking
the mean value we determine the value that minimizes the mean
average error (MAE). Results can be seen in Fig. 6. An interesting
takeaway from this figure is also that the MAE is not extremely sensitive
to the value of β<sub>OSV</sub>. Therefore, for simplicity’s sake, we will use:

$$
\beta_{OSV} = -7 \tag{13}
$$

To verify the robustness of this correlation in industrial codes,
where the friction velocity *u*<sub>*τ*</sub> is not necessarily well predicted up to a
few percent, we calculate the effect of a change on *u*<sub>*τ*</sub> on the predicted
<figure><img src="image_7.png" /><figcaption>Fig. 6. Mean average error on the high-Pe fully heated OSV database for different values of βOSV. The MAE is not extremely sensitive to βOSV. The optimal value is in the middle range of the ones seen in Fig. 5.</figcaption></figure>

$X_{OSV}$ (Fig. B.19 in Appendix B). Changing $u_t$ by 10% has a smaller effect on $X_{OSV}$, which means the correlation is robust.

More complex regressions performed using the Uranie platform (Blanchard et al., 2018) marginally improve the MAE on $X_{OSV}$ by a few %. However, we judge this improvement to small to justify the additional complexity. This also presents a risk of over-fitting the data, therefore we keep a constant $\beta_{OSV}$.

According to Eq. (6), the saturated liquid layer thickness at OSV is then:

$$y_{+,c} = 27 \approx 30 \qquad (14)$$

Before OSV, $T_l$ in this liquid sub-layer is smaller than $T_s$ and no void leaves it. At OSV, the boundary layer reaches $T_s$ and significant void can be produced. After OSV, the temperature can reach $T_s$ outside the sub-layer and vapor reaches regions that are far from the wall.

This $y_{+,c} = 30$ boundary layer size is coherent with the findings of Nop et al. (2021). They find that in a water channel with 25 K subcooling, the boiling crisis occurs when a boundary layer of thickness $\delta = \frac{80v_t}{u_s}$ reaches saturation temperature. This limit is three times higher than the one at which OSV is predicted in our model.

Recently, Kossolapov et al. (2024) measured bubble sizes at the wall in flow boiling for various pressures. The maximum dimensionless bubble size at the wall they observed was $\delta_+ = 30$. This is coherent with our finding of a saturated layer up to $y_{+,c} = 30$.

The single-phase velocity log-law in a near-wall region is valid for $y_+ > 30$ (Pope, 2000). That this is the same value as our $y_{+,c} = 30$ is remarkable but difficult to interpret.

### 3.5. Comparing our correlation to the literature

A comparison between the prediction of $X_{OSV}$ using different models for $Pe > 5 \cdot 10^4$ including annular and single-heated channel data is presented in Fig. 7. The MAE obtained on the partial database (15.77%) is comparable to that of the Saha and Zuber correlation (16.9%). Furthermore, as the error in the determination of the experimental $X_{OSV}$ is of around 10%, the observed MAE's are of an acceptable order of magnitude. We have therefore successfully built a local, CFD-scale OSV correlation valid at high $Pe$.

Fig. 8 plots the predicted results for $\beta_{OSV} = -7$ against those of Saha and Zuber (1974). The MAE (14.57%) is similar to that with the experimental data. This means that the $\beta_{OSV} = -7$ correlation is not simply a local version of the SZ correlation, or else the MAE would be significantly smaller than for the experimental comparison. This correlation is different and will predict different results in similar conditions.

To check the validity range of our correlation, we plot the results predicted by our model on the complete database and compare them with those of Saha and Zuber (Fig. 9). Some low-Pe runs with very small $X_{OSV}$ coming from Edelman and Elias (1981) are badly predicted using $\beta_{OSV} = -7$. For all of these points, grayed in Fig. 9, $Re < 2 \cdot 10^3$ and $Pe < 3.5 \cdot 10^3$. We eliminate these points for the calculation of the MAE. This gives us a lower boundary of the validity domain of the correlation. It is therefore valid for virtually all turbulent flows, which is more coherent than the $Pe = 7 \cdot 10^4$ limit found in the Saha and Zuber correlation.

In situations where $Re > 2 \cdot 10^3$, the MAE in our correlation is slightly larger than for the Saha and Zuber correlation: 21.85% vs 18.89%. However, there is only one fitted constant in our model, $\beta_{OSV}$, while there are two in SZ. The additional errors that we have come from low-Pe annular flows ($Pe \sim 1 \cdot 10^4$), i.e. the data of Rouhani (1966a,b) and Zeitoun (1994). However, high-Pe ($>10^5$) predictions for the Rouhani data are consistent with the experimental results.

## 4. Use for heat flux partitioning

### 4.1. Heat flux partitions used for comparisons with the literature

In this section, we will only study developed flows. The single-phase heat transfer law that we use for all heat flux partitions is from Kader (1981). A detailed description can be found in Appendix C.1.

<figure><img src="image_8.png" /><figcaption>Fig. 7. Comparison between models and experimental data for $Pe > 10^4$. Left: predicting $X_{OSV}$ using $\beta_{OSV} = -7$. Right: Saha and Zuber (1974) correlation. Color scale: $Pe$ number. Our correlation has a comparable MAE to the Saha and Zuber correlation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Comparison between $\beta_{OSV} = -7$ model and Saha and Zuber (1974) for $P_e > 5 \cdot 10^4$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

The first model that we choose as a reference is the original Kurul and Podowski (1990) formulation. A detailed description can be found in Appendix C.2.

The second model is a modified version of the Kurul and Podowski (1990) formulation, that is the default option in the Neptune_CFD code (Guelfi et al., 2007). The only difference with the original Kurul and Podowski formulation is the calculation of the departure diameter (Mimouni et al., 2016b). A detailed description can be found in Appendix C.3 and in Favre (2023).

## 4.2. Physical interpretation of a stationary developing boiling flow

Fig. 10 presents the different physical mechanisms occurring in a developing boiling flow. The flow enters as subcooled liquid (column ①, first row). The temperature profile is a single-phase profile (second row). The wall temperature is below the saturation temperature (third row), and increases progressively.

The ONB occurs when the wall becomes hotter than saturation temperature. The wall temperature quickly reaches a constant value and does not evolve as boiling picks off (Garnier et al., 2001). Between ONB and OSV (column ②), the temperature profiles follow Eq. (5), like those of Fig. 3, with $\beta$ becoming smaller and smaller. All of the heat flux keeps entering the liquid phase, and no vapor is produced.

After the OSV, the heat flux is split between heating the liquid and evaporation (column ③). The void fraction departs from 0. The temperature in the first element, at $y_+,1$, keeps increasing.

When it reaches saturation temperature, all of the heat flux is used for evaporation (column ④). At first, the flow remains subcooled boiling, as the core has not reached saturation temperature (pictured in the second row). The flow transitions to saturated boiling subsequently.

## 4.3. OSV-based heat flux partition

We use the physical mechanisms explained in the previous section to construct our HFP. When the wall temperature is imposed, the general philosophy is to separate the calculation of the total heat flux from that of the distribution of the heat flux between phases. Depending on the flow conditions, the former is determined using a classical single-phase heat transfer correlation or a total boiling heat flux correlation. Once the heat flux is known we use our OSV correlation to determine the heat flux distribution between phases. In imposed heat flux conditions we can directly determine the heat flux distribution without calculating the wall temperature.

In single-phase flow (region ① in Fig. 10), the heat flux at the wall is, for any distance to the wall y:

$$q_{\text{SP}} = \frac{(T_w - T_l(y))\rho_l C_p u_x}{\Theta_{+}^u(y_+)} = H_{l,\text{SP}}(y_+)(T_w - T_l(y)) \quad (15)$$

Where $H_{l,SP}(y_+)$ is the single-phase heat transfer coefficient calculated using the Kader (1981) correlation (see Appendix C.1).

Between the ONB and the OSV (region ②) the total heat flux $q_w$ will follow a boiling-flow correlation that depends of the wall temperature $T_w$. For water, we can use that of Jens and Lottes (1951) or Thom et al. (1965). For other fluids, the Frost and Dzakowic (1967) formulation can be used (from Delhaye, 2008). It was shown to give very consistent results with the DEBORA database (Gueguen, 2013). If another correlation is more appropriate for a given fluid or local conditions, as long as it only depends on near-wall quantities it can be used in our HFP without impacting the methodology. The aforementioned correlations read:

$$\begin{align}
q_{w,\text{Jens\&Lottes}} &= \left( \frac{T_w - T_s}{25} \exp(P/62) \right)^4 \nonumber \\
q_{w,\text{Thom et al.}} &= \left( \frac{T_w - T_s}{22.65} \exp(P/87) \right)^2 \tag{16}
\end{align}$$

$$q_{w,\text{Frost\&Dzakowic}} = \frac{\lambda_{ls}(h_{gs} - h_{ls})\rho_v}{8\sigma T_s} \left( \frac{T_w - T_s}{P_{rs}} \right)^2$$

Where $P_{rs}$ is the Prandtl number at saturation, $\lambda_{ls}$ the liquid conductivity at saturation and $\sigma$ the surface tension. The pressure $P$ must be in bar in the Jens and Lottes (1951) and Thom et al. (1965) formulations.

<figure><img src="image_10.png" /><figcaption>Fig. 9. Comparison between models and experimental data all experimental data. Left: predicting $X_{OSV}$ using $\beta_{OSV} = -7$. Right: Saha and Zuber (1974) correlation. Color scale: $P_e$ number. Grey points: $Re < 2000$. The MAE is calculated only using points where $Re > 2000$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 10. Physical interpretation of a stationary developing boiling flow. $T(y_+ = 30)$ is the average temperature in the boundary layer. $y_{+,1}$ is the dimensionless size of the first wall cell in a simulation. Onset of nucleate boiling (ONB) occurs when wall temperature passes saturation temperature. Onset of significant void (OSV) happens when the liquid temperature in the turbulent boundary layer reaches saturation temperature.</figcaption></figure>

We combine Eqs. (4) and (13): the heat transfer coefficient to the liquid at the OSV is:

$$H_{l,OSV}(y_+) = \frac{\rho_l C_p u_\tau}{\Theta_+^s (y_+, \beta = -7)} = \frac{\rho_l C_p u_\tau}{2.12 \log(y_+) - 7} \quad (17)$$

And the heat transfer towards the liquid phase at OSV is, at any distance y from the wall:

$$q_{l,OSV} = H_{l,OSV}(y_+)(T_s - T_l(y)) \quad (18)$$

Therefore, while the following inequality holds, the flow has not reached OSV and all of the heat flux remains in the liquid:

$$q_w < H_{l,OSV}(y_+)(T_s - T_l(y)) \quad (19)$$

The OSV occurs when inequality (19) no longer holds, i.e. when $T_l(y)$ has increased sufficiently. This is equivalent to being saturated for $y_+ < 30$. After this point (3), we use the same methodology as the system-scale models discussed previously (Section 1). We assume that after the OSV, the heat transfer coefficient towards the liquid stays the same. This is a strong hypothesis. We believe this is a lower bound on the heat transfer coefficient: as bubbles nucleate, grow and move away from the surface they are bound to increase the agitation and heat transfer efficiency compared with the situation at OSV. This enables us to calculate the heat flux towards the liquid phase. The heat flux

towards evaporation is then the difference with the total heat flux:

$$\begin{cases} q_l = H_{l,OSV}(T_s - T_l(y)) \\ q_{l \to v} = q_{lw} - H_{l,OSV}(T_s - T_l(y)) \end{cases} \quad (20)$$

Finally, once the first element has reached saturation temperature, all of the heat flux enters the evaporation term (region 4).

If $y_+ < 30$ in the first element, then vapor production will begin when $T_l = T_s$, before the OSV. We recommend that $y_+ \ge 100$ at in the near-wall cells to use the model. Similar limitations also exists for mechanistic HFP's. If the wall cell is smaller than the bubble departure diameter, the near-wall models described are no longer valid.

Combining all of these elements, we can build a heat flux partitioning algorithm. When the wall temperature is known, the inputs are $y_1$, $u_\tau$, $T_w$, $T_l(y_1)$ and the physical properties of the liquid. $y_1$ is the size of the first element. The steps of the algorithm are the following:

1. Construct a mesh so that $y_{+,1} \ge 100$

2. Calculate single-phase heat flux $q_{sp}$ using the Kader (1981) heat transfer coefficient (Eq. (15), see Appendix C.1 for details)

3. Calculate total boiling heat flux $q_{Boil}$ using a total heat flux correlation (Eq. (16))

4. If $T_l(y_1) \ge T_s: q_w = q_{Boil}$ goes into the evaporation term (region 4)

5. Else if $q_{sp} \ge q_{Boil}: q_{sp} = q_w$ goes into the liquid phase (region 1)
<figure><img src="image_12.png" /><figcaption>Fig. 11. Different steps of our heat flux partition algorithm for a fixed heat flux or a fixed wall temperature boundary condition. Inputs of the algorithm are distance to the wall in the near-wall cell $y_i$, friction velocity $u_\tau$, local liquid temperature $T_i(y_i)$ and either total heat flux $q_w$ either wall temperature $T_w$. Outputs are evaporation heat flux $q_{l→v}$ and heat flux entering liquid phase $q_l$.</figcaption></figure>

$$6. \text{ Else } q_w = q_{\text{Boil}}; \text{ calculate } H_{l,\text{OSV}} = \frac{\rho_l C_p u_\tau}{\Theta_+^s (y_{+1})}$$

(a) If $q_{l,OSV} = H_{l,OSV}(T_s - T_l(y_1)) > q_w$, $q_w$ goes into the liquid phase (region ②)

(b) Else the wall heat transfer towards the liquid phase is $H_{l,OSV}(T_s - T_l(y_1))$ and that towards the evaporation term is $q_w - H_{l,OSV}(T_s - T_l(y_1))$ (region ③)

The algorithm is shown in Fig. 11.

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

$$\begin{aligned} q_{l,OSV} &= \max \left( 0, \frac{\rho_c \mu_r (T_s - T_l(y_1))}{2.12 \log(y_{+1}) - 7} \right) \\ q_w &= \max(q_{SP}, q_{\text{Boil}}) \\ q_l &= \min(q_w, q_{l,OSV}) \\ q_{l→v} &= q_w - q_l \end{aligned} \qquad (21)$$

This presents a clear numerical advantage compared to classical HFP's for constant heat flux boundary conditions: we can skip steps 2, 3 and 5 of the algorithm. As $T_w$ is not needed for the other steps, this

means avoiding a Newton algorithm to determine the partition. This saves computation time and is easier to implement in a code.

From a physics standpoint, this guarantees a better calculation of the total heat flux than mechanistic models, as correlations directly fitted on experimental data are more precise. Furthermore, given the simplicity of the model it is easy to anticipate and interpret the outputs and the physical mechanisms at play.

The heat transfer coefficient is dependent on $u_\tau$. In this work, we only used models developed for single-phase flows. However, if reliable model for wall friction in bubbly flow are developed, comparable to that of Ramstorfer et al. (2005), $u_\tau$ will be affected. If $u_\tau$ changes between the ONB and the OSV, the situation could require either recalibrating the model, or keep using the single-phase $u_\tau$ for the thermal equation. If it changes after the OSV, this would affect the HFP after the OSV and the model would need additional verification. For now, we recommend to use this model solely with $u_\tau$ calculated from a single-phase correlation for the OSV and HFP calculations for the moment, though a multiphase formulation of $u_\tau$ could be used in the momentum equation.

### 4.4. Comparing our HFP with some from the literature

As the wall temperature, bulk liquid temperature and heat flux are the quantities most often measured in heat transfer experiments (Jens and Lottes, 1951; Thom et al., 1965; Garnier et al., 2001), the predicted wall temperature for a given heat flux is the main point of reference for a heat transfer model (Kommajosyula, 2020; Favre, 2023). Our model only predicts the partition and used external correlations for the heat transfer. It is therefore as good at predicting the wall temperature as the external correlation used, and we cannot validate it with this type of data.

Measurements of the destination of heat fluxes, to see if they go towards the liquid phase or evaporation, are few. To the best of our knowledge, such experiments were only conducted at MIT (Richenderfer et al., 2018; Kossolapov, 2021; Kossolapov et al., 2024) with a 10 mm×10 mm heater. However, our model is based on the assumption that the liquid thermal boundary layer is developed. These experiments therefore cannot be used to validate our partition.

The best we can do is to compare the predicted fluxes with models from the literature, i.e. the original Kurul and Podowski (1990) formulation and the Neptune_CFD Kurul and Podowski (Favre, 2023). To simplify notations, these will be called original KP and Neptune_CFD KP in the rest of this paper. We conduct these comparisons on data from the DEBORA experiment (Garnier et al., 2001). The liquid temperature used as an input is the one measured closest to the wall in the experiment. The y used to calculate $y_+$ and $\Theta_+(y_+)$ in our model and the single-phase heat transfer in all models is the one at which this liquid temperature is measured, i.e. $y = 0.55$ mm. The corresponding $y_{+1}$ are between 354 and 1140. The Frost and Dzakowic (1967) correlation is used for the total heat flux in our model.

Fig. 12 presents the effect of the wall temperature on the predicted heat flux for four different (Mass flux, Pressure, Liquid temperature) combinations. In the first row, one can see that the total heat flux predicted by our model (orange line) is much closer to the experimental measure (red point) than the other models. This is thanks to our use of the Frost and Dzakowic (1967) correlation. The difference between experimental heat fluxes and those predicted by the KP-based models is huge: the original KP is 3–10 times too high and the Neptune_CFD KP 1.5 to 5 times too low. Our model and the Neptune_CFD model follow the single-phase solution at low wall superheats on the section where they are superimposed. The inflection in the curve of our model, at $T_w - T_s \sim 1.5C$ in the first column for the G2P14W16 run, marks the activation of the boiling model, where $q_{SP} \approx q_{\text{Boil}}$. On the other hand, the original KP departs almost immediately from the single-phase solution.
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Impact of the wall temperature on the heat flux predicted by the original (Kurul and Podowski, 1990) and Neptune_CFD Kurul and Podowski model (Favre, 2023) and the current work vs data from the DEBORA experiment (Cubizolles, 1996). Each column represents a different (Mass flux, Pressure, Liquid temperature) combination given in the legend. **First row:** total heat flux as a function of wall superheat. **Second row:** fraction of heat flux entering the liquid phase as a function of the wall superheat. **Third row:** fraction of heat flux entering the liquid phase as a function of the total heat flux. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

The second row presents the effect of the wall temperature on the fraction of heat flux entering the liquid. The original KP model begins producing vapor as soon as the wall temperature exceeds saturation temperature. We believe this behavior to be non-physical. Our HFP needs a non-zero wall superheat before producing vapor, because of the $q_{l,OSV} = H_{l,OSV}(T_s - T_l(y)) > q_{lw}$ criterion. As expected, this occurs after the inflection for the total heat flux, at $T_w - T_s \sim 3$ °C in the first column (vs $T_w - T_s \sim 1.5$ °C). As the most common boundary condition for industrial applications is an imposed heat flux, we also plot the heat flux fraction entering the liquid phase as a function of the total heat flux (third row). This shifts the Neptune_CFD model to the left compared with the second row, as its total heat flux is smaller for a same wall superheat. For a same heat flux, our model predicts a much higher fraction of heat entering the liquid phase than the KP-based models.

**Fig. 13** focuses on the practical case of imposed-heat flux conditions. In the first row, we present the predicted wall temperatures as a function of the liquid subcooling. We can observe the transition from the single-phase regime at the right of each plot, where the subcooling is high and the wall temperature smaller than saturation temperature (black dashed line) to the boiling regime at the left of each plot. The

predicted single-phase wall temperatures are the same for all 3 models, as they all use the Kader (1981) law. In the boiling regime, the three curves separate. The original KP predicts a wall superheat two times too small. Our model is consistent with experimental data, but this is thanks to the use of the Frost and Dzakowic (1967) correlation. The Neptune_CFD KP over-predicts the wall superheat by a factor 2. The difference between experimental temperatures and temperatures predicted by the KP-based models is huge. On the experimental data, once the wall temperature has passed the saturation temperature there is a clear change of slope and the wall temperature becomes nearly constant (second and fourth columns, i.e. G2P26W16 and G5P26W24). The inflection in our model is slightly sharper than the experimental data, but much closer to the experiment than KP-based models. Our model is also the only one that predicts a constant wall temperature as the subcooling changes.

The second row presents the heat flux fraction entering the liquid as a function of the subcooling, also for a fixed heat flux. Our models transitions smoothly from sending all of the heat flux to the liquid phase at high subcoolings to evaporating the whole heat flux close to saturation temperature. The KP-based models keep heating the liquid
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Impact of the total heat flux on the wall temperature and heat flux partition predicted by the original (Kurul and Podowski, 1990) and Neptune_CFD Kurul and Podowski model (Favre, 2023) and the current work vs data from the DEBORA experiment (Cubizolles, 1996). Each column represents a different (Mass flux, Pressure, Heat flux) combination given in the legend. **First row:** Wall temperature a function of liquid subcooling at $y_+$ given in the legend. **Second row:** fraction of heat flux entering the liquid phase as a function of liquid subcooling. The dashed vertical lines represent the liquid subcooling at $y_+$ when the Saha and Zuber (1974) criterion is reached for these conditions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

after saturation temperature is reached. This should create a liquid overheat in the simulations, which was never measured in the DEBORA experiment. Both KP models begin evaporation as soon as the wall temperature goes above the saturation temperature: there is no zone between the ONB and the OSV where no vapor is produced. For our model, vapor production begins for much lower subcoolings, around 5 °C in the different plots and not 20 °C. We plot the local liquid temperature at the wall where the Saha and Zuber (1974) criterion is reached in the bulk. This temperature is calculated by interpolation of experimental data, knowing the thermodynamic quality and liquid temperature at the wall in the experiments. This is the purple vertical dashed line in Fig. 13. The beginning of vapor production in our model is consistent with the Saha and Zuber (1974) criterion, and is very far for the KP models. The latter are therefore expected to significantly over-predict the void fraction for high subcoolings compared with experimental results.

In this section, we have shown that our model predicts a more physical heat flux distribution between phases, especially at high and low subcoolings, and a coherence with the Saha and Zuber (1974) OSV model. It also enables the use of total boiling heat flux correlations that have good wall temperature predictions. More generally, we believe that to verify the physical coherence of a heat flux partitioning model, one must check that:

* Significant vapor production does not begin as soon as the wall temperature exceeds saturation temperature, but the transition is smooth.

* For a given wall temperature, the total heat flux in the boiling region is only marginally dependent on the subcooling, as is the case in experiments.

* When the liquid temperature tends towards saturation temperature, all the energy is transferred towards evaporation.

* When the liquid subcooling increases and the bulk enthalpy reaches the Saha and Zuber criterion, the fraction of heat evaporating should not be significant.

## 4.5. Bubble departure diameter calculation

As no mechanistic modeling is conducted in this model, no departure diameter is calculated. This is not an issue for our set of closures, as we do not have an interfacial area transport equation (Reiss et al., 2024). However, it is problematic for most two-phase sets of closures (Guelfi et al., 2007; Liao et al., 2018).

Both mechanistic and fitted departure diameter calculation methodologies need the wall temperature as an input (Ünal, 1976; Mazzocco et al., 2018; Kommajosyula, 2020; Favre et al., 2023). If, in the simulation, the wall temperature is given, then it is possible to calculate the departure diameter. If the heat flux is enforced, one can use a total heat flux correlation to calculate the wall temperature (Jens and Lottes, 1951; Thom et al., 1965; Frost and Dzakowic, 1967), and feed it in the departure diameter model.

## 4.6. DEBORA simulation result comparison

In order to compare the prediction of our HFP model with the original KP model, we use the two-fluid Euler–Euler module of the TrioCFD code (Angeli et al., 2015) to simulate the radial void fraction and liquid temperature in tubes from the DEBORA experiment. In these simulations, we use a $k-\omega$ single-phase turbulence model, a deformable Ishii and Zuber (1979) interfacial drag force, the Burns
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Comparison between the DEBORA database (Cubizolles, 1996), the original (Kurul and Podowski, 1990) HFP and the current work. Each column represents a different (Mass flux, Pressure, Heat flux) combination given in the legend. **Lines 1–3** present the experimental and simulated void fractions for the three coolest entrance temperatures for each combination. **Bottom row:** simulated temperatures for each combination. Black dashed line: saturation temperature. For the first two columns, experimental data is available and the results presented are for different entrance temperatures than the void fraction data. They respect the same color scale. For the two columns to the right, the simulated temperature is for the same entrance temperature as the simulated void fractions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

et al. (2004) turbulent dispersion force, and a dedicated formulation of the interfacial heat transfer in condensation and of the lift force (Reiss et al., 2024). The 2D axi-symmetric mesh used contains 20 radial and 400 axial elements. A detailed description of these closures can be found in Appendix D.

We select four (Mass velocity, Pressure, Heat Flux) conditions from the DEBORA database. We run simulations for different liquid entrance temperatures, for which the measured void fractions are small. Following the energy balance conducted by Favre (2023), we reduce by 5% the heat flux in the simulation compared with the experimental heat flux. For every simulation, the $y_+ > 100$ criterion in the near-wall cell to use our HFP is met. The predicted void fractions and liquid temperatures are presented in Fig. 14 as a function of the dimensionless radius $r_+(r_+ = 0$ at the center and $r_+ = 1$ at the wall).

The first three lines compare the experimental data (points) with the original KP model (full line) and our model (dashed line). For the lowest entrance temperatures (first row), the experimental void fraction at the wall is around 2%. The original KP model significantly over-predicts the void fraction for high subcoolings, as it sends a large

part of the heat flux in the vapor phase (see Fig. 13). Our HFP, on the other hand, does not yet predict vapor formation, apart for case G2P26W16 (second column). For higher entrance temperatures (second and third rows), our void fraction predictions are very satisfying, and much closer to the experimental data than the original KP model. As the entrance temperature increases, the difference between our model and the original KP decreases: the subcooling decreases and our HFP sends more and more energy in the evaporation term. As the $T_s - T_l$ keeps decreasing, the HFP will have less and less impact, and will no longer have any when $T_l = T_s$ in the first computation cell, where all of the energy will be used for evaporation.

The bottom row contains liquid temperature predictions for both models, as well as liquid temperature measures for runs G2P14W16 and G2P24W16 (first two columns). The liquid temperature is under predicted with the original KP model. With our HFP, more energy is injected in the liquid and less in the vapor phase. This leads to higher liquid temperatures in the bulk, by up to 1C. The temperature prediction is significantly improved.
<figure><img src="image_16.png" /><figcaption>Fig. 15. Comparison between the Bartolomei database (Bartolomei and Chanturiya, 1967; Bartolomei et al., 1982), the original (Kurul and Podowski, 1990) HFP and the current work. Each subplot represents a different (Mass flux, Pressure, Heat flux, Tube diameter) combination given in the legend.</figcaption></figure>

The heat flux distribution between liquid phase and evaporation as a function of the liquid subcooling is fairly similar between the original KP and the Neptune CFD KP models (see bottom row of Fig. 13). As this simulation has an imposed total heat flux, this distribution is the lone difference between HFP's and we expect the Neptune_CFD KP model to also significantly over predict the void fraction at high subcoolings, albeit less than the original KP.

## 4.7. Bartolomei simulation result comparison

R12, a refrigerant fluid, was used in the DEBORA setup so similar density ratio and boiling numbers as in pressurized water reactors (PWR) could be obtained, with operating pressures and heat fluxes respectively 5 and 10 times smaller (Cubizolles, 1996; Garnier et al., 2001). To verify that the results obtained in the previous section are still valid in PWR conditions, we would ideally use a database of local void fraction measures in pressurized flow boiling. However, such measures are not available in the literature. Therefore, we run simulations on the Bartolomei and Chanturiya (1967) and Bartolomei et al. (1982) experiments. These are part of our OSV database (see Section 2). The average void fraction was measured at different heights in 12, 15.4 and 24 mm diameter tubes at various mass fluxes, pressures and heat fluxes. We use the same set of closures as in Section 4.6, with 2D axi-symmetric meshes that have 20 radial and 400 axial elements. Results are presented in Fig. 15.

The OSV predicted by the Kurul and Podowski (1990) model occurs too early. The simulated void fraction at the experimental OSV is extremely high, around 10% in most cases. Our HFP predicts an OSV that is coherent with the experimental results, greatly improving the simulations compared with the Kurul and Podowski (1990) model. However, after the OSV the void fraction increases faster with our HFP than in the experiments. This could be an issue with our HFP: if the presence of a departed bubble layer increases the efficiency of the heat transfer towards the liquid, then the use of $H_{i,OSV}$ (Eq. (17)) as the liquid heat transfer coefficient underestimates the transfer towards the liquid and overestimates vapor production. Other terms could be the cause of this discrepancy: the condensation term could be insufficient in the bulk; a too low turbulent dispersion or lift force could prevent enough vapor from reaching the very subcooled regions near the core. For high thermodynamic qualities, once all of the liquid has reached saturation temperature, both HFP's and the experimental data converge to the same result.

## 5. Conclusion

We have developed an OSV model validated for tube, channel and annular geometries, for Peclet numbers ranging from $3.5 \cdot 10^3$ to $4 \cdot 10^5$ and pressures from 1 bar to 147 bar. This model is simple, has only one fitted coefficient, and performs as well as the Saha and Zuber (1974) correlation. When the total wall heat flux is known, it can be used to calculate the distribution of the heat flux between the transfer to the liquid and evaporation. This improves boiling-flow simulations compared with the Kurul and Podowski (1990) model for imposed-flux boundary conditions. When the wall temperature is known, our model can be used after an empirical boiling correlation has been employed to predict the heat flux.

Future work on this topic could include verifying if our model is valid for the prediction of OSV in transient power spikes, during reactivity insertion accidents (RIA) in nuclear reactors for example.

### CRediT authorship contribution statement

**Corentin Reiss:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Antoine Gerschenfeld:** Writing – review & editing, Supervision, Project administration, Methodology, Conceptualization. **Catherine Colin:** Writing – review & editing, Supervision, Methodology, Conceptualization.

### Declaration of competing interest

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

### Data availability

A link to a GitHub repository is given in the paper.

### Acknowledgments

We thank Anthony Quintin for assistance with the numerical analysis, and Guillaume Bois, Alan Burlot, Raksmy Nop and Elie Roumet for productive physics discussions.

### Appendix A. Comparing our data points to other references

Here we compare our datapoints with datapoints from the literature. Fig. A.16 is a comparison with Saha and Zuber (1974), Fig. A.17 with Cai et al. (2021) and Fig. A.18 with Lee and Jeong (2022).
<figure><img src="image_17.png" /><figcaption>**Fig. A.16.** Comparison between the data presented in the original (Saha and Zuber, 1974) article and our data collection. **Left:** Peclet vs Stanton at OSV on the Saha and Zuber database. Compare with figure 2 from Saha and Zuber (1974). **Right:** Peclet vs Stanton at OSV, using the datapoints that we collected and that are cited in Saha and Zuber (1974). We considered the void fraction measures in Maurer (1960) and Evangelisti and Lupoli (1969) not precise enough to determine *X*<sub>OSV</sub>. Furthermore, Dix (1971) defined the OSV as the point where bubbles are able to cross the channel, which is different from the one we selected, i.e. using the plot of the average void fraction as a function of *X*. This is why these 3 sources are absent from our database. Our points are consistent with those of Saha and Zuber, apart from a group of low-*St* low-*Pe* datapoints from Martin (1972) that come from high-pressure channels with a small hydraulic diameter.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. A.17.** Comparison between the data presented in the article from Cai et al. (2021) and our data collection. **Left:** *X*<sub>OSV</sub> prediction by the Saha and Zuber correlation vs experimental data on the Cai et al. database. Compare with figure 2.(a) from Cai et al. (2021). **Right:** *X*<sub>OSV</sub> prediction by the Saha and Zuber correlation vs experimental data, using the data points that we collected and that are used in Cai et al. (2021). Our points are consistent with those of Cai et al..</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. A.18.** Comparison between the (Lee and Jeong, 2022) article and our data collection. **Left:** Peclet vs Stanton at OSV on the Lee and Jeong database. Compare with figure 11 from Lee and Jeong (2022). 499 datapoints collected. **Right:** Peclet vs Stanton at OSV, using the datapoints that we collected and that are cited in Lee and Jeong (2022). We have much fewer datapoints (155 in total vs 972), as we chose not to use many references in their database. The points that are plotted seem consistent though we cannot identify which ones come from which source.</figcaption></figure>
## Appendix B. Sensitivity to friction velocity

Fig. B.19 presents the sensitivity of the $\beta_{OSV}$ calculation to the prediction of the friction velocity $u_r$.

<figure><img src="image_20.png" /><figcaption>Fig. B.19. Sensitivity of $X_{OSV}$ predicted using $\beta_{OSV} = -7$ when changing $u_r$.</figcaption></figure>

## Appendix C. Heat flux partitions used for comparisons with the literature

### C.1. Kader single-phase heat transfer

This model is proposed by Kader (1981). We place ourselves at a distance $y$ from the wall.

We calculate:

$$ \beta_{SP} = (3.85(Pr^{1/3}) - 1.3)^2 + 2.12 \log(Pr) \quad (C.1) $$

And a transition coefficient:

$$ \gamma = \frac{0.01(Pry_+)^4}{1 + 5Pr^3y_+} \quad (C.2) $$

We can then determine:

$$ \Theta_{+}^{w}(y_{+}) = Pry_{+} \exp(-\gamma) + (2.12 \cdot \log(1 + y_{+}) + \beta_{SP}) \exp\left(-\frac{1}{\gamma}\right) \quad (C.3) $$

And finally:

$$ q_l = (T_{lw} - T_l(y)) \frac{\rho_l C_\rho u_r}{\Theta_{+}^{w}(y_+)} \quad (C.4) $$

### C.2. Original Kurul and Podowski

This model is the one originally proposed by Kurul and Podowski (1990). We place ourselves at a distance $y$ from the wall.

The departure diameter is a linear interpolation between those of Ünal (1976) and Thomas (1981):

$$ d_{\text{departure}} = 10^{-4}(T_{lw} - T_s) + 0.0014 \quad (C.5) $$

The nucleation site density is from Del Valle and Kenning (1985):

$$ N_{\text{sites}} = (210(T_{lw} - T_s))^{1.8} \quad (C.6) $$

The area influenced by the presence of bubbles is:

$$ A_{\text{bubbles}} = \min \left( 1, \frac{\pi}{4} N_{\text{sites}} d_{\text{departure}}^2 \right) \quad (C.7) $$

The bubble departure frequency is given by Cole (1960):

$$ f_{\text{departure}} = \sqrt{\frac{4 \cdot 9.81(\rho_l - \rho_v)}{3\rho_l d_{\text{departure}}}} \quad (C.8) $$

The quenching heat flux is also from Del Valle and Kenning (1985), where the waiting period used is $t_{wait} = 1/f_{departure}$:

$$ q_{\text{quench}} = 2A_{\text{bubbles}}(T_{lw} - T_l(y))f_{\text{departure}} \sqrt{\frac{t_{\text{wait}}\lambda_l\rho_l C_\rho}{\pi}} \quad (C.9) $$

Using the Kader single-phase heat transfer (see Appendix C.1), this yields the following transfer toward the liquid phase:

$$ q_l = (1 - A_{\text{bubbles}})q_{\text{SP}}(y) + q_{\text{quench}}(y) \quad (C.10) $$

And the evaporation heat flux is:

$$ q_{l \to v} = \frac{\pi}{6} f_{\text{departure}} d_{\text{departure}}^3 \rho_v (h_{vs} - h_{ls}) N_{\text{sites}} \quad (C.11) $$

It can be noted that the evaporation heat flux in this formulation is independent of the local liquid temperature, while the quenching and single-phase heat flux are not.

The total heat flux is:

$$ q_{lw} = (1 - A_{\text{bubbles}})q_{\text{SP}}(y) + q_{\text{quench}}(y) + q_{l \to v} \quad (C.12) $$

### C.3. Neptune_CFD Kurul and Podowski

This model is a variation of the original Kurul and Podowski (1990) model found in the Neptune_CFD code. A detailed description and analysis can be found in Favre (2023) and Mimouni et al. (2016b). The difference with the original Kurul and Podowski model is the departure diameter. It requires the bulk liquid velocity $u_{bulk}$ as an input.

$$ d_{\text{departure}} = 2.4 \cdot 10^{-5} P^{709} \frac{a}{\sqrt{b\phi}} \quad (C.13) $$

$$ a = \frac{(T_{lw} - T_s)\lambda_l}{2\rho_v(h_{vs} - h_{ls})\sqrt{\frac{\pi\lambda_{lw}}{C_{pw}\rho_w}}} \quad (C.14) $$

Where $\lambda_w$ is the wall conductivity, $C_{pw}$ the wall heat capacity and $\rho_w$ the wall density.

$$ a = \frac{(T_{lw} - T_s)\lambda_l}{2\rho_v(h_{vs} - h_{ls})\sqrt{\frac{\pi\lambda_{lw}}{C_{pw}\rho_w}}} \quad (C.15) $$

$$ b = \begin{cases} \frac{T_s - T_l(y)}{2(1 - \frac{\rho_l}{\rho_v})} & \text{if } St \le 0.0065 \\ \frac{T_s - T_l(y)}{2(1 - \frac{\rho_l}{\rho_v})} \frac{(1 - A_{\text{bubbles}})q_{\text{SP}}(y) + q_{\text{quench}} + q_{l \to v}}{\rho_l C_\rho u_{\text{bulk}} 0.0065} & \text{if } St > 0.0065 \end{cases} \quad (C.16) $$

Where $St = \frac{(1-A_{\text{bubbles}})q_{\text{SP}}(y)+q_{\text{quench}}+q_{l \to v}}{\rho_l C_\rho u_{\text{bulk}}(T_s-T_l(y))}$ is defined using the bulk velocity.

$$ \phi = \max \left( 1, \left( \frac{u_{\text{bulk}}}{0.61} \right)^{0.47} \right) \quad (C.17) $$

## Appendix D. Two-fluid model used in our simulations

A detailed description of the two-fluid model that we used in our simulations is given in Reiss et al. (2024). A mass, momentum and energy equation are solved for the liquid and vapor phase (Ishii and Hibiki, 2006). The semi-conservative form of the momentum equation is used (Park et al., 2009). The equations used to govern a phase $k$ are:

$$ \frac{\partial \alpha_k \rho_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k \vec{u}_k) = \Gamma_k $$

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

$$ -P \left( \partial_t \alpha_k + \nabla \cdot (\alpha_k \rho_k \vec{u}_k) \right) + \nabla \cdot \left[ \alpha_k \lambda_k \nabla T_k - \alpha_k \rho_k \vec{u}_k' \vec{e}_k' \right] + q_{ki} \frac{\partial \alpha_k \rho_k e_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k e_k \vec{u}_k) = $$

$\vec{u}_k' \vec{u}_k'$ and $\vec{u}_k' \vec{e}_k'$ are the turbulent tensors, $\Gamma_k$ the mass transfer term, $\vec{F}_{ki}$ the interfacial force term and $q_{ki}$ the bulk interfacial heat transfer. The wall heat transfer to phase $k$, $q_{kiw}$, is a boundary condition of the energy equation.
A $k-\omega$ model is used for shear-induced turbulence (Kok, 1999), with a constant turbulent Prandtl number fixed at 0.9 for the energy equation:

$$v_t = \frac{k}{\omega} \qquad \overline{u'_i u'_j} = -v_i \nabla \vec{u}_i \qquad \overline{u'_i e'_i} = -v_i C_{p,i} \nabla T_i$$

$$\partial_t k + \nabla \cdot (k \vec{u}_i) = v_i (\nabla \vec{u}_i + \vec{t} \vec{u}_i) \cdot \nabla \vec{u}_i - \beta_k k \omega \\ + \nabla \cdot ((v_i + \sigma_k v_i) \nabla k) \qquad (D.2)$$

$$\partial_t \omega + \nabla \cdot (\alpha_i \omega \vec{u}_i) = \alpha_\omega (\nabla \vec{u}_i + \vec{t} \nabla \vec{u}_i) \cdot \nabla \vec{u}_i - \beta_\omega \omega^2 \\ + \nabla \cdot ((v_i + \sigma_\omega v_i) \nabla \omega) + \sigma_d \frac{1}{\omega} \max\{\nabla k \cdot \nabla \omega, 0\}$$

The main specificity of our set of closures is that it does not require an interfacial area transport equation. Therefore, no bubble diameter is used in the subsequent equations.

The interfacial force exerted by the liquid on the gas is $\vec{F}_{gl} = -\vec{F}_{il}$. All forces written here apply to the gas phase. We separate the interfacial force term in three different contributions:

$$\vec{F}_{gl} = \vec{F}_{drag} + \vec{F}_{lift} + \vec{F}_{TD} \qquad (D.3)$$

We implement the deformed bubble drag force of Ishii and Zuber (1979):

$$\vec{F}_{drag} = -\frac{3}{4} C_D \frac{\alpha_v \rho_l}{d_b} \|\vec{u}_v - \vec{u}_l\| (\vec{u}_v - \vec{u}_l) \\ C_D = \frac{2}{3} \frac{d_b}{L_c}, \quad L_c = \sqrt{\frac{\sigma}{g(\rho_l - \rho_v)}} \qquad (D.4)$$

The lift force used was calibrated in Reiss et al. (2024), it is:

$$\vec{F}_{lift} = -C_L \rho_l \alpha_v (\vec{u}_v - \vec{u}_l) \wedge (\underline{\nabla} \wedge \vec{u}_l) \\ C_L = \begin{cases} 0 & \text{if } \alpha_v < 0.25 \\ \max(-0.2, -0.7 \cdot (\alpha_v - 0.25)) & \text{if } 0.25 \le \alpha_v < 0.7 \\ \alpha_v - 0.9 & \text{if } 0.7 \le \alpha_v < 0.9 \\ 0 & \text{if } 0.9 \le \alpha_v \end{cases} \qquad (D.5)$$

We select the Burns et al. (2004) force, with a turbulent Prandtl number of 1:

$$\vec{F}_{disp} = -C_{TD} \rho_l k \nabla \alpha_v, \quad C_{TD} = \frac{3 C_D}{4 d_b} |\vec{u}_v - \vec{u}_l| \frac{1}{\omega} \left(1 + \frac{\alpha_v}{\alpha_l}\right) \qquad (D.6)$$

The interfacial heat transfer term used was also calibrated in Reiss et al. (2024). It reads:

$$L_c = \sqrt{\frac{\sigma}{g(\rho_l - \rho_v)}}, \quad Nu_B = 30 \\ q_{ki} = \alpha_v Nu_{fit} \frac{6\alpha_l}{L_c^2} \min\left(1, \sqrt{\frac{\pi}{3\sqrt{2} \min(\alpha_v, 0.6) - 1}}\right) (T_v - T_l) \qquad (D.7)$$

Where $L_c$ is the capillary length and $Nu_B$ the bubble Nusselt number.

## References

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

Anne, R.D., Beattie, D.R.H., 1996. Net vapour generation onset during high velocity subcooled boiling flow: A simple turbulence-based prediction method. In: Inc., B.H. (Ed.), Heat and Mass Transfer Australasia. pp. 503–509.

Baglietto, E., Demarly, E., Kommajosyula, R., 2019. Boiling crisis as the stability limit to wall heat partitioning. Appl. Phys. Lett. 114 (10), 103701. http://dx.doi.org/10.1063/1.5080724.

Bartolomei, G., Brantov, V.G., Molochnikov, Y.S., Kharitonov, Y., Solodkii, V.A., Batashova, G., Mikhailov, V.N., 1982. An experimental investigation of true volumetric vapor content with subcooled boiling in tubes. Therm. Eng. 29 (3), 132–135. http://dx.doi.org/10.2172/1262488.

Bartolomei, G.G., Chanturiya, V.M., 1967. Experimental study of true void fraction when boiling subcooled water in vertical tubes. Therm. Eng. 14 (2), 123–128.

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

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

Blanchard, J.B., Damblin, G., Martinez, J.M., Arnaud, G., Gaudier, F., 2018. The uranie platform: an open-source software for optimisation, meta-modelling and uncertainty analysis. arXiv preprint arXiv:1803.10656.

Blasius, H., 1913. Das aehnlichkeitsgesetz bei reibungsvorgängen in flüssigkeiten. In: Mitteilungen ÜBer Forschungsarbeiten Auf Dem Gebiete Des Ingenieurwesens: Insbesondere Aus Den Laboratorien Der Technischen Hochschulen. Springer, pp. 1–41.

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

Cai, C., Mudawar, I., Liu, H., Xi, X., 2021. Assessment of void fraction models and correlations for subcooled boiling in vertical upflow in a circular tube. Int. J. Heat Mass Transfer 171, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2021.121060.

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

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

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

Delhaye, J.-M., 2008. Thermohydraulique des réacteurs. EDP Sciences.

Dix, G.E., 1971. Vapor void fractions for forced convection with subcooled boiling at low flow rates. In: Technical Report. GE Report.

Edelman, Z., Elias, E., 1981. Void fraction distribution in low flow rate subcooled boiling. Nucl. Eng. Des. 66 (3), 375–382. http://dx.doi.org/10.1016/0029-5493(81)90167-9.

Egen, R.A., Dingee, D.A., Chastain, J.W., 1957. Vapor formation and behavior in boiling heat transfer. In: Technical Report. (BMI-1163), Battelle Memorial Inst., Columbus, Ohio.

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

Evangelisti, R., Lupoli, P., 1969. The void fraction in an annular channel at atmospheric pressure. Int. J. Heat Mass Transfer 12, 699–711. http://dx.doi.org/10.1016/0017-9310(69)90004-0.

Favre, L., 2023. Modeling and simulation of the boiling crisis within PWR at CFD scale (Ph.D. thesis). Institut National Polytechnique de Toulouse-INPT, URL https://theses.hal.sciences tel-04244931.

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 (ISSN: 0017-9310) 211, 124227. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2023.124227.

Ferrell, J.K., 1964. A study of convection boiling inside channels. In: Technical Report. North Carolina State University, Raleigh, North Carolina, http://dx.doi.org/10.2172/4598223.

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

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

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

Gueguen, J., 2013. Contribution à la modélisation multidimensionnelle des écoulements bouillants convectifs en conduite haute pression pour l'application au cas des réacteurs à eau pressurisée (Ph.D. thesis). Université de Grenoble, URL https://theses.hal.sciences.tel-01685226/.

Guelfi, A., Bestion, D., Boucker, M., Boudier, P., Fillion, P., Grandotto, M., Hérard, J.-M., Hervieu, E., Péturaud, 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.

Ha, T.-W., Yun, B.-J., Jeong, J.J., 2020. Improvement of the subcooled boiling model for thermal-hydraulic system codes. Nucl. Eng. Des. 364 (110641), http://dx.doi.org/10.1016/j.nucengdes.2020.110641.

Ishii, M., Hibiki, T., 2006. Thermo-fluid dynamics of two-phase flow. Springer Science and Business Media.

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

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

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

Kader, B.A., Yaglom, A.M., 1972. Heat and mass transfer laws for fully turbulent wall flows. Int. J. Heat Mass Transfer 15, 2329–2351. http://dx.doi.org/10.1016/0017-9310(72)90131-7.