# Void drift in narrow rectangular channels for highly advective bubbly flows: a macroscopic drift-flux model derived from two-fluid local simulations

Guillaume Bois, Robin Barbera, Matthieu Amate-Vignon

## ► To cite this version:

Guillaume Bois, Robin Barbera, Matthieu Amate-Vignon. Void drift in narrow rectangular channels for highly advective bubbly flows: a macroscopic drift-flux model derived from two-fluid local simulations. International Journal of Multiphase Flow, 2024, 182, pp.105038. 10.1016/j.ijmultiphaseflow.2024.105038. cea-04834989

HAL Id: cea-04834989

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

Submitted on 12 Dec 2024

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

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

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

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

# Void drift in narrow rectangular channels for highly advective bubbly flows: A macroscopic drift-flux model derived from two-fluid local simulations

G. Bois*, R. Barbera, M. Amate-Vignon

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

## ARTICLE INFO

**Keywords:**

Void drift

Two-phase flow

Volume averaging

Homogenisation

Up-scaling

Sub-channel modelling

Drift-flux model

Two-fluid model

Turbulent dispersion

Effective diffusion

Void fraction transport

Constitutive relationships

## ABSTRACT

This article addresses the issue of reduced models to describe turbulent two-phase flows in industrial applications. This work is an application based on Bois (2021) theoretical derivations to connect the spatially-averaged drift-flux model (DFM) to the local Reynolds-averaged two-fluid model (TFM). It presents new opportunities to calibrate closure laws or propose new models based on finer-scale descriptions. Highly convective steady bubbly flows with wall-peaking void fraction profiles are considered. Dispersed bubbles are weakly deformable. A new model dedicated to void fraction dispersion in flat rectangular channels is developed for application into the DFM. The model is calibrated based on reference local Euler–Euler two-fluid simulations of pressurised water and steam mono-dispersed bubbles, in adiabatic conditions. 64 conditions are considered to cover regularly the 3D parameter space in void fraction, Reynolds and Eötvös numbers. The new formulation proposed shows a very good fit with the post-processed predictions of local CFD simulations along with a substantial reduction in grid requirement. The model contains sub-filter correlations between vapour concentration and mixture or relative velocities. This model derived from local simulations is validated a posteriori in an industrial component scale code where the full space-filtered DFM is implemented and resolved. This work proves the benefits and feasibility of one-by-one model development and calibration based on two-fluid simulations.

## 1. Introduction

Accurate predictions of multiphase flows is very complex due to the wide range and time- and space-scales naturally involved. Yet, for industrial applications such as transport, chemical reactors or heat exchangers, it is essential to provide an acceptable description of the macroscopic system at a reasonable computational cost. To reach this goal, macroscopic models have been developed for engineering tools. They rely on time- and space-averaged equations mostly resorting to empirical correlations to describe interfacial or wall transfers, statistical cross-correlations (such as Reynolds Stresses in turbulent flows) and dispersion due to space-averaging or wall-laws.

Local two-fluid models (TFM, Drew and Passman, 1999; Ishii and Hibiki, 2010; Morel, 2015) have been developed based on Reynolds-Averaged Navier-Stokes (RANS) equations where the evolution of each phase is described separately by a set of conservation equations. Averaging conservation principles leads to an under-determined set of equations and additional local relations are necessary to achieve a determined set; they are usually inferred from experimental measurements (see e.g. Ishii and Zuber, 1979, for the determination of the drag force). Most often, experiments can only provide partial insights

into flow characteristics and the closure relations had to be developed and calibrated based on assessment of their global effect. For instance, models for drag and lift forces are evaluated based on the overall capability of a complete solver to provide reasonable predictions of void-fraction distribution and level for a range of physical conditions covered by the experiments. Consequently, experiments struggle to clearly isolate the individual role of each phenomena and of each model, in realistic flow conditions. In the past, due to strong convection in industrial applications, efforts were mostly directed towards modelling of dominant effects in the streamwise direction, neglecting the importance of cross-flow lateral interactions (lift, wall effect and turbulent dispersion).

Yet, local two-fluid models are still very challenging and hard to use in full-scale industrial applications; their validation remains difficult and incomplete. As a consequence, further reduction of computational costs are achieved by **averaging techniques and homogenisation** (Whitaker, 1967; Anderson and Jackson, 1967; Whitaker, 1969; Quintard and Whitaker, 1994; Whitaker, 1999) to incorporate the effect of (time-averaged) local gradients and replace the need to accurately

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: guillaume.bois@cea.fr (G. Bois).</footnote>
Nomenclature

Acronyms

<table><tr><td>CFD</td><td>Computational Fluid Dynamics</td></tr><tr><td>CHF</td><td>Critical Heat Flux</td></tr><tr><td>DFM</td><td>Drift-Flux Model</td></tr><tr><td>HEM</td><td>Homogeneous Equilibrium Model</td></tr><tr><td>IAC</td><td>Interfacial Area Concentration</td></tr><tr><td>MCFD</td><td>Multiphase Computational Fluid Dynamics</td></tr><tr><td>PDE</td><td>Partial Differential Equation</td></tr><tr><td>PWR</td><td>Pressurised Water Reactor</td></tr><tr><td>RANS</td><td>Reynolds-Averaged Navier Stokes</td></tr><tr><td>REV</td><td>Representative Elementary Volume</td></tr><tr><td>TFM</td><td>Two-Fluid Model</td></tr></table>

Subscript

<table><tr><td>F</td><td>Favre</td></tr><tr><td>f</td><td>Fluid</td></tr><tr><td>k</td><td>k phase</td></tr><tr><td>l</td><td>Liquid phase</td></tr><tr><td>M</td><td>Macroscopic (mixture and filtered)</td></tr><tr><td>m</td><td>Mixture</td></tr><tr><td>v</td><td>Vapour phase</td></tr><tr><td>w</td><td>Wall</td></tr></table>

Greek symbols

<table><tr><td>α</td><td>Void fraction</td></tr><tr><td>χ</td><td>Phase indicator function</td></tr><tr><td>χ<sub>f</sub></td><td>Indicator function of the fluid phase</td></tr><tr><td>δ<sub>w</sub></td><td>Dirac delta function associated to the wall</td></tr><tr><td>γ<sub>0</sub></td><td>Isothermal friction multiplier for two-phase flows</td></tr><tr><td>μ</td><td>Viscosity</td></tr><tr><td>φ</td><td>Porosity</td></tr><tr><td>ρ</td><td>Density</td></tr><tr><td>σ</td><td>Surface tension</td></tr></table>

Latin

<table><tr><td>K<sub>M</sub></td><td>Tensor set of coefficients to close momentum in Eq. (13)</td></tr><tr><td>T<sup>t</sup></td><td>Turbulent stress</td></tr><tr><td>T<sup>v</sup></td><td>Viscous stress</td></tr><tr><td>T<sup>dr</sup></td><td>Diffusion stress</td></tr><tr><td>u&#773;<sub>D</sub></td><td>Macroscopic drift velocity, u&#773;<sub>D1</sub> + u&#773;<sub>D2</sub>, Eq. (9)</td></tr><tr><td>u&#773;<sub>m</sub></td><td>Barycentric velocity</td></tr><tr><td>u&#773;<sub>r</sub></td><td>Relative velocity</td></tr><tr><td>u&#773;</td><td>Velocity</td></tr><tr><td>D<sup>r</sup></td><td>Macroscopic diffusion tensor (⟨T<sup>dr</sup>⟩<sub>f</sub>)</td></tr><tr><td>D<sub>M</sub></td><td>Macroscopic dispersion</td></tr><tr><td>F<sub>M</sub></td><td>Macroscopic wall friction</td></tr><tr><td>I</td><td>Identity tensor</td></tr><tr><td>M<sub>l→k</sub></td><td>Interfacial transfer into phase k</td></tr><tr><td>S<sub>M</sub></td><td>Deviatoric part of the macroscopic stress tensor, Eq. (13)</td></tr><tr><td>T</td><td>Macroscopic turbulent tensor (⟨T<sup>t</sup>⟩<sub>f</sub>)</td></tr><tr><td>V</td><td>Macroscopic viscous tensor (⟨T<sup>v</sup>⟩<sub>f</sub>)</td></tr><tr><td>A<sub>d</sub></td><td>Parameter for the transverse diffusion, Eq. (17)</td></tr></table>

<table><tr><td>a<sub>i</sub></td><td>Interfacial area concentration</td></tr><tr><td>c</td><td>Vapour mass fraction</td></tr><tr><td>C<sub>f</sub></td><td>Single-phase friction coefficient</td></tr><tr><td>d<sub>b</sub></td><td>Sauter mean diameter</td></tr><tr><td>D<sub>h</sub></td><td>Hydraulic diameter</td></tr><tr><td>D<sup>α</sup><sub>d</sub></td><td>Drift coefficient</td></tr><tr><td>D<sub>sc</sub></td><td>Drift diffusion coefficient</td></tr><tr><td>e</td><td>Gap dimension or channel width</td></tr><tr><td>E<sub>o</sub></td><td>Eötvös number</td></tr><tr><td>g</td><td>Gravitational acceleration</td></tr><tr><td>K<sub>0</sub></td><td>Parameter for the transverse diffusion, Eq. (17)</td></tr><tr><td>K<sub>c</sub></td><td>Transverse diffusion coefficient, Eq. (17)</td></tr><tr><td>l<sub>ξ</sub></td><td>Microscopic scale of variation of the variable ξ</td></tr><tr><td>L<sub>⟨ξ⟩</sub></td><td>Macroscopic scale of variation of the filtered variable ⟨ξ⟩</td></tr><tr><td>p</td><td>Pressure</td></tr><tr><td>q<sub>q</sub></td><td>Quenching heat flux</td></tr><tr><td>r<sub>0</sub></td><td>Characteristic length-scale of the filter kernel</td></tr><tr><td>Re</td><td>Reynolds number</td></tr><tr><td>t</td><td>Time</td></tr><tr><td>V</td><td>Averaging volume</td></tr><tr><td>y</td><td>Wall-normal coordinate</td></tr></table>

Mathematical symbols

<table><tr><td>δ<sub>F</sub>φ</td><td>Deviation of φ from the Favre filtering ⟨φ⟩<sub>F</sub>, Eq. (A.8b)</td></tr><tr><td>δ<sub>f</sub>φ</td><td>Deviation of φ from the fluid average ⟨φ⟩<sub>f</sub></td></tr><tr><td>δ<sub>s</sub>φ</td><td>Deviation of φ from the volume average ⟨φ⟩</td></tr><tr><td>⟨⟩<sub>F</sub></td><td>Favre averaging (fluid- and density-weighted filtering), Eq. (A.8a)</td></tr><tr><td>φ&#773;<sup>k</sup></td><td>Phase average of φ</td></tr><tr><td>φ&#773;</td><td>Statistical average of φ</td></tr><tr><td>⟨⟩<sub>f</sub></td><td>Weighted volume average - intrinsic filtering</td></tr><tr><td>⟨⟩</td><td>Volume average or space filter</td></tr><tr><td>∇<sup>T</sup></td><td>Transposed gradient</td></tr><tr><td>φ'</td><td>Fluctuation of φ with respect to the phase average</td></tr><tr><td>D<sub>k</sub>/Dt</td><td>Material or substantial derivative with respect to phase k</td></tr></table>

Superscript

<table><tr><td>t</td><td>Turbulent</td></tr><tr><td>v</td><td>Viscous</td></tr></table>

capture boundary layers and the effect of complex geometries by additional models. These homogenisation techniques are at the basis of macro-scale or macroscopic descriptions where filtered (space-averaged) variables describe the system. Homogenisation permits drastic reduction of the requirement on the mesh degrees of freedom. When it is used to average over the whole cross-sectional area, it leads to one-dimensional models popular in engineering tools. This methodology leads to models referred to as sub-channel or component models and codes. Several works (Drew and Passman, 1999; Ishii and Hibiki, 2010; Wallis, 1969; Delhaye and Achard, 1976; Delhaye et al., 1981a,b; Yeoh and Tu, 2016) present a specific derivation for targeted applications with appropriate simplifications. Here, we consider averaging per unit length over a single cross-flow direction.
In conjunction to homogenisation, if the coupling between the two-phases is strong enough to be modelled by algebraic relations, the two-fluid approach can be replaced by a drift-flux model (DFM, Hibiki and Ishii, 2003) or even an homogeneous equilibrium model (HEM) if no relative velocity is considered between the phases. Then, only mixture balance equations are considered and non-equilibrium between the phases for the velocity can be captured by additional transport equations or by algebraic closures (see e.g. Zuber, 1964, for one of the first application of the homogenised DFM to gas-liquid flows in vertical or horizontal pipes). The DFM is obtained as a simpler formulation; most applications and modelling were performed in one dimension and generalisation of the DFM's closures to multi-dimensional cases is not straightforward.

An innovative strategy was proposed in Bois (2021). It relies on the strong connections between these various approaches. Showing explicitly the relations between local and filtered approaches on one hand, and between two-fluid and mixture quantities on the other, Bois (2021) derived a set of equations for a macroscopic mixture model for two-phase turbulent flows based on local conservation laws and principles. Even though this equation system is under-determined (it is not closed in the sense that the number of unknown is larger than the number of equations) and as such, cannot be solved directly, it provides a sound basis for the derivation of independent closures for each model; indeed, each physical phenomena that requires modelling in the macroscopic DFM is expressed rigorously in terms of local two-fluid variables. These revealed definitions open the way to new techniques for model development that can be based on local two-fluid model predictions.

The main advantages of the approach are that it enables separate model calibration, even in complex flow configurations where concurrent mechanisms occurs, and that it ensures that the model proposed is consistent with local principles granted by local conservation equations. However, the inherent limit is that it requires and assumes a valid local two-fluid model. Following one of the prospects presented in Bois (2021), this paper presents an application of this up-scaling methodology. As an example of particular interest, this paper is focused on vapour cross-flow distribution and spreading under dynamical mechanisms. As a first step, adiabatic conditions are considered and physical properties are assumed constant and taken representative of high-pressure boiling conditions as encountered in industrial Pressurised Water Reactor (PWR) applications. This void drift mechanism is of particular concern to the nuclear safety for which vapour accumulation at the wall is determinant to the boiling regime (and possibly to determine conditions of Critical Heat Flux, CHF, occurrence). One advantage of numerical simulations is that the problem can be partially simplified and that the necessary complexity introduced in the problem is controlled. Here for instance, the case considered is simpler than true boiling applications where validation of the local two-fluid model is very challenging. Adiabatic conditions eliminate the need for the energy equation and complex wall-heat partitioning models. Besides, in high-pressure conditions, the flow remains made of dispersed almost-spherical bubbles; thus we consider in this first step a mono-dispersed flow, without needing complex interfacial area transport equations nor coalescence or break-up models. Extension to phase-change can be considered afterwards.

In this article, the two-step up-scaling strategy from the local TFM to the macroscopic (sub-channel) drift-flux description is applied to the void-fraction dispersion mechanism in large aspect-ratio rectangular channels with upwards mono-dispersed bubbly flows. To cope with limited validation of the TFM, the study is restricted to adiabatic highly convective flows with wall peaking void fraction profiles. The dispersed phase is made of weakly deformable bubbles. Our simulations reveal unsteady flow when natural convection takes the lead on advection. This unsteady regime is excluded from the up-scaling analysis while experimental investigations are conducted to confirm its physical existence.

Section 2 presents a review on void drift modelling, mostly in the context of a different sub-channel geometry. Then, Section 3 gives a brief reminder of the major elements established in Bois (2021) on which the methodology and modelling are based. It is followed by the description of the local TFM model and the volume-averaged DFM with theirs closures implemented in Neptune_CFD and Flica5 respectively. Neptune_CFD is used to perform the reference simulations in the conditions described in Section 5. Physical analysis of these results is presented in Section 6. Then, post-processing local TFM results provides reference data for a priori validation and calibration of void-drift closures for the filtered DFM (Section 7). Section 8 is dedicated to a posteriori assessment of the newly derived models and their comparison to other void drift models. Finally, conclusions are drawn in Section 9.

## 2. Review on void drift modelling

The concept of space-averaging applied to void fraction dispersion models was mostly derived for the nuclear industry in conjunction with a specific filter corresponding to subchannels in rod-bundle geometries. Cross-flow redistributions occur along the channel axis due to 3 mechanisms (Collier, 1979; Lahey and Moody, 1977): (i) lateral pressure difference between subchannels (Bowring, 1968; Rowe and Angle, 1969), (ii) void drift (due to a cross-flow occurring in the non-equilibrium flow irrespective of the lateral pressure difference) and (iii) turbulent mixing (without net mass or volume exchanges between adjacent subchannels for each phase). This subdivision introduced by Lahey and Moody (1977) was meant to help develop subchannel models but still very few models for cross-flow contributions exist.

Void drift can be represented by the void-settling model (Lahey et al., 1972; Gonzalez-Santalo and Griffith, 1972) to predict the volume flow rate of the gas phase moving from subchannel i to j

$$ \frac{dQ_{gij}}{dz} = \bar{D}\Delta_{ij}(\alpha - \alpha_{eq}) \qquad (1) $$

where $dQ_{gij}$ is the volume flow rate of the gas phase moving from subchannel i to j over the interval $dz$ and $\Delta_{ij}\alpha = \alpha_i - \alpha_j$. The coefficient of void diffusion $\bar{D}$ and the equilibrium void fraction difference $\Delta_{ij}\alpha_{eq}$ have to be determined experimentally.

Sadatomi et al. (1994) studied flow redistribution due to void drift in two-phase flows between two connected subchannels. A tentative subchannel analysis was presented based on experimental data on redistribution along the channel axis due to void drift under a condition without lateral pressure difference between the subchannels. In their analysis, the effects of void drift and turbulent mixing between subchannels were taken into account simultaneously by incorporating both the void-settling model (Lahey et al., 1972; Gonzalez-Santalo and Griffith, 1972) from Eq. (1) and a turbulent mixing rate. Their experimental data supports the fact that the void settling model in Eq. (1) is debatable as it requires both the subchannel void fractions equilibrium states $\delta\alpha_{ijEQ}$ and calibration of $\bar{D}$. The only model to assess void fractions equilibrium proposed by Levy (1963) is totally discredited by Sadatomi et al. (1994) data; besides, the calibration of $\bar{D}$ to a single dimensional value of $\bar{D} = 1.9 \times 10^{-3} \text{ m}^2 \text{ s}^{-1}$ is specifically dependent on the geometrical setup and fluids used. Even though they attempted to acquire experimental data on the flow redistribution process approaching an equilibrium flow owing to void drift, Sadatomi et al. (1994) concluded that further studies are needed to understand the mechanism of void drift. In their following works (Sadatomi et al., 2004; Kawahara et al., 2008), void drift phenomena is still modelled by a turbulent mixing, calibrated on their specific experiments, and by the void settling model for lateral mass flux of gas phase due to the void drift with an experimentally determined void diffusion coefficient.

Similar conclusions are drawn later by Otero et al. (2011) that performed rod bundle mixing experiments. They concluded that sub-channel codes are developed to predict cross-flow between neighbouring subchannels, and yet they require a number of assumptions and
empirical adjustments. Their results corroborate the great uncertainty associated with the turbulent mixing determination and their study insists on the need for further knowledge for subchannel systems.

For some industrial codes, the void drift model is presented. Carver et al. (1987) detail the DFM in the ASSERT subchannel code. The constitutive relationships for the relative velocity are presented as the heart of the successful model. It comprises the effects of (i) cross-section averaging, (ii) gravity separation and (iii) turbulent void diffusion. While (i) and (ii) are classically represented by Zuber and Findlay's model (Zuber and Findlay, 1965) with a distribution coefficient $C_0$ and a weighted mean velocity $V_{gj}$, turbulent diffusion (also called void drift) (iii) is expressed by a diffusion term. Thus, relative velocity $\bar{u}_R$ is modelled as

$$ \bar{u}_R = \left[ (C_0 - 1) j + V_{gj} - \frac{\epsilon}{l\alpha} \Delta_{ij} (\alpha - \alpha_{eq}) \right] / (1 - \alpha) \quad (2) $$

where j is the volumetric flux of the mixture and l is a characteristic (undefined) length between channels. Notably, Carver et al. (1987) explain that according to the usual practice in subchannel analysis, void diffusion is considered dominant in the lateral momentum equation, and only the last two terms in the numerator of Eq. (2) are used, as $C_0$ is taken as unity in this direction without further justification.

The void diffusion term (iii) includes two effects: the turbulent diffusion of void fraction and the diffusion toward a nonuniform void. The equilibrium void fraction $\alpha_{eq}$ is fitted on experimental observations to account for void migration towards larger subchannels (Lahey, 1974). The void eddy diffusivity coefficient $\epsilon$ is modelled as a Peclet number

$$ P e_{\alpha} = \frac{\epsilon}{VD} = 0.075 (\alpha/0.6)^6 \quad (3) $$

without specific precision on the definition of characteristic velocity $V$ and length $D$ scales. This model seems to rely on Rudzinski's work (Rudzinski, 1970) who produced two-phase turbulent mixing data in an experiment of air-water flows in adjacent triangular array subchannels.

Rudzinski (1970) found the water turbulent mass flux between channels proportional to the mean driving force related to the disequilibrium of void fraction between the subchannels. Fick's first law of diffusion was used to define the mass transfer coefficient, introducing the concept of eddy diffusivity of mass or heat $\epsilon$ and relating it to the ratio between bulk velocity and the hydraulic diameter. Prandtl mixing length theory is advocated to support the fact that the mixing coefficient is inversely proportional to the gap width.

In Rudzinski (1970), transverse mixing flow rate due to turbulent mixing $W^t$ is assumed to be given by $W^t = \rho_l \epsilon / d$ where d is an estimation of a characteristic distance between two sub-channels similar to l in Eq. (2). Water flows from the low to the high quality subchannel whilst air goes in the opposite direction. Water mass flux is proportional to the disequilibrium in liquid fraction for all conditions whereas the air turbulent mass flux is only proportional to the void fraction disequilibrium between channels when it is below 8% and become independent from it afterwards.

Furthermore, Carlucci et al. (2004) provide a review of general observations for two-phase turbulent mixing and experimental datasets; they propose an original formulation for the intersubchannel mixing rate per unit distance which involves complex expressions for incremental mixing parameters and flow obstruction factors. The phase diffusion coefficients $\epsilon_k$ are calculated from the corresponding mixing rates towards an equilibrium void fraction. Contrary to Tincq et al. (1991), mixing rate models depends on the void fraction instead of quality.

In the french thermal-hydraulic code for steam generators THYC, Tincq et al. (1991) model the relative velocity by Zuber and Findlay model (distribution coefficient $C_0$ and velocity $V_{gj}$) for (i) and (ii) while void drift (iii) is modelled by a diffusion coefficient $D_{sc}$ and the gradient of vapour mass concentration $c$:

$$ \bar{u}_R = f(C_0, V_{gj}) - \frac{D_{sc}}{c(1-c)} \nabla c \quad (4) $$

No further detail on the calibration of $D_{sc}$ could be found, except that it is generally considered proportional to the single-phase turbulent viscosity (Caremoli et al., 1993). It seems to be dependent on the channels' geometry and mixing elements (baffles, ...). In Jouhanique and Rascle (1995), $D_{sc}$ is taken equal to the single-phase turbulent viscosity by a similarity with mass diffusion and a Schmidt number.

Carlucci et al. (2004) stress that both the level of understanding and the supporting data on two-phase mixing are significantly less than for single-phase flow (Weisman and Bowring, 1975), where constitutive relationships already showed a wide spreading of predictions for the same condition. In their recent review, Moorthi et al. (2018) provide a comprehensive list of codes and experimental datasets for subchannels analyses, but void drift modelling is very briefly discussed.

It is also interesting to point out that in the thorough derivation of the area-averaged DFM by Ishii and Hibiki (2010), the drift mechanism is totally overlooked by the one-dimensional assumption in Chap. 14. The effects of local drift velocity (noted $V_{dj}$) and its area-averaged contribution (noted $\langle\langle V_{dj} \rangle\rangle$) are proved to be negligible compared to the effect of the distribution parameter $C_0$ in view of experimental evidences; however, the reasoning only considers the streamwise velocity component and is not valid for the analysis of the lateral migration. In the lateral direction, both the distribution parameter $C_0$ and the averaged drift velocity $\langle\langle V_{dj} \rangle\rangle$ are important in Eq. (14.22) of Ishii and Hibiki (2010) but this topic is not discussed further. Still, for a part the local drift velocity, Ishii and Hibiki (2010) postulates in Eq. (13.32) the existence of a drift coefficient $D_d^a$ in a linear constitutive law to account for the concentration gradient of the form: $-D_d^a \nabla \alpha / \alpha$. Here again, it is said that the effect of gravity and forces is usually the dominant part of the drift velocity and the contribution of this drift coefficient can be neglected but this consideration is applied in the streamwise and gravity direction. In a brief mention of horizontal flows, the argument of lack of knowledge leads to the assumption $\langle\langle V_{dj} \rangle\rangle \approx 0$ without real justification.

To conclude this brief overview, we recall that constitutive equations are very specific to the geometry and space-filtering considered. Most considerations seem to have been related to sub-channels in rod-bundles geometries, where disequilibrium is encountered in some parts of a nuclear fuel assembly. Diffusion of void fraction is scarcely analysed from an averaged point-of-view in the literature and we did not find applications to different geometries. In this work, we consider a very different averaging, over the gap of a narrow rectangular channel. To the best of our knowledge, void drift modelling has not been applied to this kind of geometry in the literature.

The filtering and homogenisation proposed here can be applied in free regions of the flow and not necessarily to porous regions where the goal is to eliminate the need for the description of the actual solid structure topology and to replace it by modelled sources. Indeed, our objective is to describe the flow at the macroscopic scale (i.e., to capture only the gradients of spatially averaged fields) and to account for high-shear regions indirectly by means of models.

### 3. Up-scaling concept

This section presents the up-scaling strategy as illustrated on Fig. 1. Two important aspects define each system: (i) the consideration of the fluid phases and (ii) the temporal and spatial scale examined. To go from the local TFM to the volume-averaged DFM, two paths can be taken (see Fig. 1). Here, following Bois (2021), we consider the local DFM derivation by variable reduction, and then, the application of volume-averaging to obtain the final macroscopic description (volume-averaged DFM). The sets of governing PDEs for these three descriptions are described in Sections 3.1, 3.2, and 3.3, respectively. Statistical averaging, volume-averaging and variable reduction are presented in Appendix A with the necessary variable definitions and specific relations. The other path on Fig. 1 leads to the same relations with the natural introduction
<figure><img src="image_3.png" /><figcaption>Fig. 1. Up-scaling description with the scales considered, the fluid phases and the operations connecting them. Definitions of scales: local instantaneous, statistical and macroscopic (both statistical and volume-averaged). Consideration of the two phases: two-fluid or drift-flux. Connecting mechanisms: statistical averaging (dotted horizontal arrow), volume averaging (continuous horizontal arrow), model reduction (dashed vertical arrow).</figcaption></figure>

of the macroscopic relative velocity $\bar{u}_R$ defined from phase-velocities instead of the drift velocity $\bar{u}_D$.

Then, computer codes solve the local TFM (Neptune_CFD) or the volume-averaged DFM (Flica5) and post-processing produces reference data for both local and volume-averaged DFMs. The closures used for each code are given in Section 4 and the description of post-processing is detailed in Section 7 and on Fig. 11.

In the context of this work, physical properties are considered constant and the flow is adiabatic. Thus, the energy equation is not considered.

### 3.1. Local TFM

Based on the local description of conservation equations for each phase and jump conditions at the interface (Kataoka, 1986), two-fluid models can be established by application of a statistical averaging operator (Ishii and Hibiki, 2010; Morel, 2015). At the local scale, for adiabatic conditions without phase change, they rely on two sets of statistically-averaged conservation equations governing mass and momentum for $k \in [l, v]$:

$$ \frac{D_k (\alpha_k \rho_k)}{Dt} = 0, \qquad (5a) $$

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

where $\alpha_v = \bar{\chi}_v$ is the void fraction, $D_k(\alpha_k \rho_k \bar{\phi}_k^k)/Dt = \partial(\alpha_k \rho_k \bar{\phi}_k^k)/\partial t + \nabla \cdot (\alpha_k \rho_k \bar{u}_k^k \bar{\phi}_k^k)$ is the material or substantial derivative of any quantity $\bar{\phi}_k^k$ associated to the phase velocity $\bar{u}_k^k$. Phase-averaged quantities $\bar{\psi}_k^k = \bar{\chi}_k \psi_k / \bar{\chi}_k$ result from statistical (Reynolds) average $\bar{\psi}$ weighted by the average $\bar{\chi}_k$ of the phase-indicator function $\chi_k \cdot \mathbf{u}, p, \rho$ and $\mu$ denote respectively the velocity, pressure, and constant density and dynamic viscosity of each phase. $g$ is the gravitational acceleration. Local pressure equilibrium is assumed. $\bar{\mathcal{T}}_k^v$ is the mean viscous stress approximated by $\bar{\mathcal{T}}_k^v = \alpha_k \mu_k (\nabla \bar{u}_k^k + \nabla^T \bar{u}_k^k)$. Closure relationships are

required to model the turbulent Reynolds stresses $\bar{\mathcal{T}}_k^l = \rho_k \bar{u}_k^l \bar{u}_k^k$ and the interfacial momentum transfer $\mathbf{M}_{l \to k}$.

The equation system (5) is associated to closure relations in Section 4.1 to produce reference results on void drift in Section 7. An alternative (but equivalent) description of the two-fluid model based on relative and mixture velocities can be derived (Bois, 2021). When the PDE for the relative velocity is replaced by a physics-based algebraic closure relation, the local drift-flux approach (local DFM) is obtained.

### 3.2. Local DFM

The focus is now set on the local mixture system which results of the manipulation of the local two-fluid equation system (5); the void fraction $\alpha$ is replaced by the vapour mass fraction $c$, the phase densities $\rho_l$ and $\rho_v$ by the mixture density $\rho_m$, and the phase velocities $\bar{u}_l^l$ and $\bar{u}_v^v$ by the mixture and relative velocity from vapour to liquid, respectively $\bar{u}_m$ and $\bar{u}_r = \bar{u}_v^v - \bar{u}_l^l$.

The mixture mass $\rho_m$ and momentum $\rho_m \bar{u}_m$ governing equations are obtained from the equation system (5). Mass conservation for the vapour phase (Eq. (5a) for $k = v$) is kept to characterise the void fraction transport while the relative velocity $\bar{u}_r$ is provided by an appropriate algebraic closure relation. Differential equations driving this system consists of the mixture mass conservation, the void fraction transport and the mixture momentum conservation (Ishii and Hibiki, 2010)

$$ \frac{\partial \rho_m}{\partial t} + \nabla \cdot (\rho_m \bar{u}_m) = 0, \qquad (6a) $$

$$ \frac{\partial \rho_m c}{\partial t} + \nabla \cdot (\rho_m c \bar{u}_m) = -\nabla \cdot (\rho_m c \bar{u}_{v \to m}), \qquad (6b) $$

$$ \frac{\partial \rho_m \bar{u}_m}{\partial t} + \nabla \cdot (\rho_m \bar{u}_m \bar{u}_m) = -\nabla \rho_m + \rho_m g \qquad (6c) $$

$$ -\nabla \cdot \left( c \rho_m \bar{u}_v^k \bar{u}_v^l + (1-c) \rho_m \bar{u}_l^k \bar{u}_l^l \right) - \bar{\mathcal{T}}_l^v - \bar{\mathcal{T}}_v^v - c(1-c) \rho_m \bar{u}_r^k \bar{u}_r^l \qquad \begin{array}{l} \tau^l: \text{Turbulent} \\ -\tau^v: \text{Viscous} \\ \tau^{dr}: \text{Diffusion} \end{array} $$
where the diffusion velocity $\bar{\mathbf{u}}_{v \to m} = \bar{\mathbf{u}}_v^v - \bar{\mathbf{u}}_m = (1-c)\bar{\mathbf{u}}_r$ is defined by the difference of the vapour velocity to the mixture velocity. A closure for this velocity $\bar{\mathbf{u}}_{v \to m}$ provides both closures to the void dispersion (last term in Eq. (6b)) and to consistently express the diffusion stress due to relative velocity $\mathcal{T}^{dr}$, as defined in Eq. (6c). The system must also be completed with expressions for the viscous and turbulent Reynolds stresses, $\mathcal{T}^v$ and $\mathcal{T}^t$ respectively (see Bois, 2021, for more details). The surface tension effect and the pressure difference are neglected.

Here, we consider a flow of dispersed bubbles where the two phases are strongly coupled dynamically in the streamwise direction; this DFM is then particularly relevant and it is also known as diffusion model or mixture model (Ishii and Hibiki, 2010, p. 103). Lastly, a macroscopic model can be derived by volume averaging.

### 3.3. Volume-averaged DFM

In Bois (2021), we rigorously derived the macroscopic mixture model by application of a general filter. Under current approximations and focusing on a filter restricted to the fluid domain (no porosity), the mixture model described by the equation system (6) can be filtered as

$$ \frac{\partial \rho_M}{\partial t} + \nabla \cdot (\rho_M \bar{\mathbf{u}}_M) = 0, \quad (7a) $$

$$ \frac{\partial \rho_M c_M}{\partial t} + \nabla \cdot (\rho_M c_M \bar{\mathbf{u}}_M + \rho_M \bar{\mathbf{u}}_D) = 0, \quad (7b) $$

$$ \frac{\partial \rho_M \bar{\mathbf{u}}_M}{\partial t} + \nabla \cdot (\rho_M \bar{\mathbf{u}}_M \bar{\mathbf{u}}_M) = -\nabla p_M + \rho_M g - \nabla \cdot (\mathbf{D}_M + \mathbf{V} + \mathbf{T} + \mathbf{D}^r) - \mathbf{F}_M. \quad (7c) $$

where the main variables are $\rho_M = \langle \rho_m \rangle_f$ the filtered mixture density, $c_M = \langle \rho_m c \rangle_f / \rho_M$ the vapour concentration, $\bar{\mathbf{u}}_M = \langle \rho_M \bar{\mathbf{u}}_m \rangle_f / \rho_M$ the centre of mass velocity and $\rho_M = \langle \rho_m \rangle_f$ the filtered pressure. The simple filter $\langle \xi \rangle$, the fluid filter $\langle \xi \rangle_f$, the Favre-averaged filter $\langle \xi \rangle_f$ and their associated deviations $\delta_\xi^\xi$, $\delta_f^\xi$ and $\delta_\xi^\xi$ are defined in Appendix A by Eqs. (A.2), (A.5), (A.8a) and Eq. (A.8b).

Besides, $\mathbf{D}_M = -\rho_M \langle \delta \bar{\mathbf{u}}_m \delta \bar{\mathbf{u}}_m \rangle_F = \rho_M \bar{\mathbf{u}}_M \bar{\mathbf{u}}_M - \langle \rho_M \bar{\mathbf{u}}_m \bar{\mathbf{u}}_m \rangle_f$ is the macroscopic dispersion due to local mean velocity profiles (correlation between mean velocity deviations $\delta \bar{\mathbf{u}}_m$ arising in the convective term), $\mathbf{D}^r = \langle \mathcal{T}^{dr} \rangle_f$ is the diffusion due to the relative velocity, $\mathbf{V} = \langle \mathcal{T}^v \rangle_f$ and $\mathbf{T} = \langle \mathcal{T}^t \rangle_f$ correspond respectively to the macroscopic viscous and turbulent diffusion tensors.

At the wall, viscous stresses $\mathcal{T}^v$ and a pressure term are the only remaining contributions as velocities ($\bar{\mathbf{u}}_m$, $\bar{\mathbf{u}}_r$ and $\bar{\mathbf{u}}_{v \to m}$) and turbulent stresses $\mathcal{T}^t$ go to zero due to the no-slip condition. They result in wall friction $\mathbf{F}_M$ given by

$$ \mathbf{F}_M = \langle \delta_f \rho_M \delta_w \mathbf{n} \rangle_f + \langle \mathcal{T}^v \delta_w \mathbf{n} \rangle_f \quad (8) $$

where $\delta_f \rho_M = \rho_M - \langle \rho_m \rangle_f$ is the pressure deviation at the wall.

In Eq. (7b), the sum in the divergence is a construction of the vapour mass flux $\rho_M c_M \bar{\mathbf{u}}_v$. The macroscopic drift velocity $\bar{\mathbf{u}}_D$ defined by

$$ \bar{\mathbf{u}}_D = c_M (1 - c_M) \bar{\mathbf{u}}_R = \frac{\langle \rho_M c \delta \bar{\mathbf{u}}_m \rangle_f}{\rho_M} + \frac{\langle \rho_M c (1-c) \bar{\mathbf{u}}_r \rangle_f}{\rho_M} \quad (9a) $$

$$ \bar{\mathbf{u}}_D = \underbrace{\langle \delta c \delta \bar{\mathbf{u}}_m \rangle_F}_{\bar{\mathbf{u}}_{D1}} + \underbrace{\langle c(1-c)\bar{\mathbf{u}}_r \rangle_F}_{\bar{\mathbf{u}}_{D2}} = \langle \delta c \delta \bar{\mathbf{u}}_m \rangle_F + \langle c \bar{\mathbf{u}}_{v \to m} \rangle_F \quad (9b) $$

should be modelled to represent both correlations between deviations of vapour concentration and mean velocity, and correlations of vapour concentration to the relative velocity.

This paper focuses on the modelling of the macroscopic drift velocity $\bar{\mathbf{u}}_D$. Additional work is required to fully close the equation system (7) by providing closures to $\mathbf{D}_M$, $\mathbf{V}$, $\mathbf{T}$, $\mathbf{D}^r$ and wall friction $\mathbf{F}_M$. As demonstrated in Bois (2021), the viscous term $\mathbf{V}$ can be closed solely resorting to the separate closures of $\bar{\mathbf{u}}_{D1}$ and $\bar{\mathbf{u}}_{D2}$.

<table><caption>**Table 1**
Flow dynamics models in the TFM of Neptune_CFD.</caption><thead><tr><th></th><th>Model</th><th>Neptune_CFD</th></tr></thead><tbody><tr><td rowspan="2">Turbulence</td><td>Continuous</td><td>R<sub>ij</sub> - ε SSG</td></tr><tr><td>Bubble-induced turbulence</td><td>S<sub>R<sub>ij</sub></sub>, S<sub>e</sub> (Bois, 2017)</td></tr><tr><td rowspan="7">Interfacial forces</td><td>Disperse</td><td>None</td></tr><tr><td>Drag</td><td>Ishii and Zuber (1979)</td></tr><tr><td>Added-mass</td><td>Zuber (1964)</td></tr><tr><td>Lift</td><td>Tomiyama et al. (2002)</td></tr><tr><td>Wall-damping</td><td>None</td></tr><tr><td>Wall force</td><td>Tomiyama (1998)</td></tr><tr><td>Turbulent dispersion</td><td>Laviéville et al. (2015)</td></tr></tbody></table>

## 4. Codes and models

The equation system of the different scales have been introduced in their open exact form in the previous section. Now, we present the sets of closure relations used in the RANS Euler-Euler TFM of Neptune_CFD and in the volume-averaged DFM of Flica5.

### 4.1. Local statistical description: the RANS Euler-Euler TFM (Neptune_CFD)

The local TFM description is based on equation system (5). This work relies on the local two-phase flow solver for a Euler-Euler RANS approach with a single pressure implemented in Neptune_CFD V6.3.0 (developed since Guelfi et al., 2007). Reference models recommended in the code VnV (Verification & Validation) document are used (see Table 1).

Transport equations predict liquid turbulent stresses $\mathcal{T}^t$ and eddy dissipation. A second order turbulence model based on SSG (Speziale et al., 1991) is used in the liquid phase (Mimouni et al., 2009, 2010a,b, 2011; Bois, 2017) with two interfacial source terms $\Pi_{ij}^t$ and $P_e^t$ for Reynolds stresses and dissipation respectively (see e.g. Bois, 2017, for their formulations). The wall law of Reichardt (1951) and the Mac Adams correlation for the friction velocity (McAdams and McAdams, 1954) are used. Turbulence is neglected in the vapour phase.

Besides, interfacial momentum transfer is modelled as a series of forces satisfying $\sum_k M_{i \to k} = 0$, thus neglecting the possible mixture momentum due to surface tension $\mathbf{M}_m$. These forces include drag (Ishii and Zuber, 1979), added-mass and Basset (Zuber, 1964), lift (Tomiyama et al., 2002) and wall lubrication (Tomiyama, 1998) forces. The lift force model of Tomiyama et al. (Eq. 11, 2002) is based on the maximum horizontal dimension of a deformed bubble given by the empirical correlation of Wellek et al. (Eq. 34, 1966) for the calculation of the modified Eötvös number $E_{o_d}$. The lift coefficient originally proposed only for $E_{o_d} < 4$ is used for any diameter:

$$ C_L = \min \left[ 0.288 \tanh \left( 0.121 \mathrm{Re}_b(d_b) \right); f(E_{o_d}) \right] \quad (10) $$

Turbulent dispersion is accounted for by an additional force inspired by the derivation in Laviéville et al. (Eq. 13, 2015) with some differences:

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

$$ C_{TD} = (F_d \tau^t - 1) \left( \frac{b + \eta_r}{1 + \eta_r} + \frac{\alpha_g}{1 - \alpha_g} \right) + C_{AM} \frac{b^2 - b}{1 + \eta_r} + \frac{\rho_g}{\rho_l} \frac{b^2 + \eta_r}{1 + \eta_r} \\ \eta_r = \frac{\tau_{lg}^t}{\tau_{lF}^t}, \quad b = \frac{\rho_l + \rho_l C_{AM}}{\rho_g + \rho_l C_{AM}}. $$

$$ \tau_{lg}^t = \frac{3}{2} C_\mu \frac{k_l}{\epsilon_l} \left( 1 + \beta \frac{V_r^2}{k_l} \right)^{-1/2}, \quad \tau_{lF}^t = \frac{1}{F_d} \left( \frac{\rho_g}{\rho_l} + C_{AM} \right) \quad (12) $$

with $\beta = 2.7$, $C_{AM} = 0.5$, $V_r \approx \bar{\mathbf{u}}_r$ and $F_d = 3C_d \bar{\mathbf{u}}_r/(4d_b)$ ($C_d$ is the drag coefficient given by Ishii and Zuber (1979)).

In this work, interfacial mass transfer is not considered and adiabatic flows are studied. The Interfacial Area Concentration (IAC) $a_i$ is simply given by the relation $a_i = 6\alpha/d_b$. To simplify the models
required in adiabatic conditions, we consider a constant diameter $d_b$ set in agreement with the Eötvös number of each case (see Section 5.4).

An adaptive time-stepping algorithm is used (constrained by a CFL below unity and Fourier numbers below 10) with a velocity prediction algorithm coupling mean and relative velocities. Inlet turbulent profile (at $z = 0$ m) is established through a recycling boundary condition prescribing the average flow-rate.

## 4.2. Local DFM

To the best of our knowledge, there is no tool available to resolve the local DFM, in particular due to the lack of appropriate closures for all the components of the local relative velocity. However, data provided by the local TFM is processed to compute mixture variables and provide references for this description. Then, space-averaging is used to create references for the macro-scale solution and models.

## 4.3. Volume-averaged DFM (Flica5)

At the component scale (space-averaged resolution), variables $c$, $\rho_m$, $\delta c$, $\delta \bar{\mathbf{u}}_m$ and $\bar{\mathbf{u}}_r$ are no longer accessible because they are local variables combined into subgrid terms before the averaging operator. Similarly, the difference between macroscopic vapour and liquid velocities $\bar{\mathbf{u}}_R = \bar{\mathbf{u}}_V - \bar{\mathbf{u}}_L$ (with $\bar{\mathbf{u}}_K = \langle \alpha_k \bar{\mathbf{u}}_k^k \rangle_f / \langle \alpha_k \rangle_f$ for $K \in [L, V]$) cannot be computed either (only the macroscopic mixture velocity $\bar{\mathbf{u}}_M$ is computed).

The filtered DFM description is resolved in the code Flica5 V1.0.11; it is a component code for nuclear applications developed by CEA. It is mostly used for nuclear core modelling but can handle various configurations and flow behaviours. Considering the strong influence of gravity on two-phase flows and for convection dominated flows, streamwise and transverse directions are very different. For the discussions here, we consider an upward flow in the direction $z$ aligned with gravity, and a filtering operator applied to the thin direction of the cross-section $y$. Direction $x$ is referred to as lateral or transverse direction. The equation system solved in Flica5 corresponds to the $y$-averaged mixture approach given by the equation system (7) in Section 3.3. Besides, due to the strongly convective conditions, streamwise velocity gradients are neglected ($\partial \bar{\mathbf{u}}_{Mz}/\partial z \approx 0$, $\partial(\rho_M \bar{\mathbf{u}}_{Dz})/\partial z \approx 0$, etc.).

In the adiabatic case, models for turbulent dispersion of mixture momentum, for wall friction and for the drift velocity are considered. Concerning dispersion and turbulence in the momentum equation (7c), several mechanisms are modelled together based on a single-phase mixing length model

$$ \mathbf{D}_M + \mathbf{V} + \mathbf{T} + \mathbf{D}^r = \mathbf{K}_M \mathbf{S}_M = \frac{\mathbf{K}_M}{2} (\nabla \bar{\mathbf{u}}_M + \nabla^T \bar{\mathbf{u}}_M) \quad (13) $$

where $\mathbf{K}_M$ is a tensor set of coefficients and $\mathbf{S}_M$ is the deviatoric part of the macroscopic stress tensor.

Besides, wall friction is modelled by the product of a single-phase coefficient $C_f$ by an isothermal friction multiplier for two-phase flows $\gamma_0$

$$ \mathbf{F}_M = -\frac{C_f \gamma_0 \rho_M \bar{\mathbf{u}}_M \|\bar{\mathbf{u}}_M\|}{2D_h} \quad (14) $$

where $D_h \approx 2e$ is the hydraulic diameter and $e$ is the gap dimension in direction $y$.

Finally, the void fraction transport equation needs a closure for the drift velocity. Streamwise and transverse drift velocities are considered separately. Their exact definition is given by Eq. (9) which shows that it is composed of two subgrid terms. In Flica5, the drift velocity $\bar{\mathbf{u}}_D$ in the void propagation equation (7b) is approximated as a diffusive gradient of the vapour mass concentration for its transverse component while the axial component is handled via experimental correlations for relative velocity. For historical reasons (component codes were originally considering only 1D sub-channels), in Flica5, the model for relative velocity does not include transverse terms but only streamwise components. To the best of our knowledge, there is no correlation

<table><caption>**Table 2**
Models proposed for the transverse drift velocity. fit means that the coefficient is a function of the case, i.e., $(\alpha_{in}, Re_l, Eo)$.</caption><thead><tr><th>Model</th><th>Label</th><th>Comment</th><th>C<sub>0</sub></th><th>K<sub>0</sub></th><th>B<sub>M</sub></th><th>C<sub>M</sub></th><th>A<sub>d</sub></th></tr></thead><tbody><tr><td>0</td><td>std.</td><td>Historical</td><td>1</td><td>K<sub>0</sub><sup>F5</sup></td><td>B<sub>M</sub><sup>F5</sup></td><td>C<sub>M</sub><sup>F5</sup></td><td>0</td></tr><tr><td>1</td><td>Uniform</td><td>Constant coefficient for each case, optimised at z = 1.5 m</td><td>0</td><td>fit</td><td>0</td><td>0</td><td>0</td></tr><tr><td>2</td><td>Local K<sub>0</sub></td><td>Variation of K<sub>0</sub> with z</td><td>0</td><td>K<sub>0</sub>(z, fit)</td><td>0</td><td>0</td><td>0</td></tr><tr><td>3</td><td>Local A<sub>d</sub></td><td></td><td>0</td><td>0</td><td>0</td><td>0</td><td>fit</td></tr></tbody></table>

in the literature able to predict relative velocity in the cross-flow direction (see review in Section 2). Indeed, a void-drift model specific to the geometry and the averaging is used in the lateral direction. In practice, this void-drift model is also computed in the streamwise direction without sound justifications but this contribution is negligible ($\partial(\rho_M \bar{\mathbf{u}}_{Dz})/\partial z \approx -\partial(K_c \partial c_M/\partial z)/\partial z \approx 0$).

Therefore, the void propagation equation (7b) is written

$$ \frac{\partial \rho_M c_M}{\partial t} + \nabla \cdot (\rho_M c_M \bar{\mathbf{u}}_M) + \frac{\partial (\rho_M c_M (1 - c_M) \bar{\mathbf{u}}_{Rz})}{\partial z} - \nabla \cdot (K_c \nabla c_M) = 0 \quad (15) $$

where the term $-K_c \nabla c_M$ accounts for the whole contribution $\rho_M \bar{\mathbf{u}}_D$ in the transverse direction ($x$) based on an approximation by a diffusion coefficient analogous to Fourier's or Fick's laws, while a separate model is used to predict $\bar{\mathbf{u}}_{Rz}$ in the direction of gravity. In addition, note that due to the strong convective velocity in this direction, the streamwise gradients ($\partial z$) are essentially negligible.

The streamwise component $\bar{\mathbf{u}}_{Rz}$ is modelled via the correlation (Ishii, 1977) for an averaged prediction of the drift velocity due to buoyancy

$$ \bar{\mathbf{u}}_{Rz} = \frac{(C_0 - 1) \bar{\mathbf{u}}_{Mz} + \bar{\mathbf{u}}_{z,lim}}{1 - c_M + (c_M - \alpha_M) C_0} \quad (16) $$

where expressions for the distribution coefficient $C_0$ and the bubble terminal velocity $\bar{\mathbf{u}}_{z,lim}$ are recalled in Appendix B, with respectively Eqs. (B.5a) and (B.5b) for the bubbly flow regime. The transverse closure of the drift velocity is the main focus of this work and several models are presented at the end of this section (Table 2). Model 0 refers to the historical version used by default in Flica5. New proposals are based on proportionality either to the Reynolds number (1 & 2 using $K_0$ parameter) or to the void fraction (3, using $A_d$ parameter). They are assessed a priori in Section 7 and validated a posteriori in Section 8.

The transverse component $\bar{\mathbf{u}}_{Dx}$ is modelled by means of the coefficient $K_c$ which is put in the general form:

$$ K_c = C_0 + K_0 (Re_{local} - Re_M)^{B_M} \gamma_0^{C_M} + A_d \alpha_M (1 - \alpha_M)^2 \quad (17) $$

with $Re_{local} = \rho_M \|\bar{\mathbf{u}}_M\| D_h / \mu_l$ the local Reynolds number and $C_0$, $K_0$, $Re_M$, $B_M$, $C_M$ and $A_d$ user defined coefficients. The contribution related to $K_0$ comes from a crude analogy to the scalar diffusion in single-phase flows while the contribution proportional to $A_d$ is inspired by empirical studies such as Sharma and Nayak (2016), Carlucci et al. (2004) and Sadatomi et al. (2004) for totally different geometries. Linear combination of the contributions in Eq. (17) is possible; however, our analysis is limited to separate assessment of the different parts as shown in Table 2 to avoid simultaneous fit of $K_0$ and $A_d$ with an additional degree of freedom. Note that by construction, this transverse closure for $\bar{\mathbf{u}}_{Dx}$ could also be used in any term involving the transverse relative velocity by Eq. (9a).

## 5. Application

The up-scaling strategy is applied to adiabatic two-phase upward flows in a thin rectangular channel with the aim to provide a macroscopic model (space-averaged) for void-fraction dispersion (void drift). The numerical set-up is described in Section 5.1 followed by fluid properties and mesh characteristics. Finally, the physical parameter space is given in Section 5.4.
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Illustration of the thin rectangular channel with the boundary conditions for wall, outlet and mixture injection *I* (blue), or liquid injection in the area *II* (red). Along *y*-axis, the injection is uniform and vapour rapidly transition either to a wall-peaked (red, left curve) or to a core-peaked (blue, right curve) void fraction profile. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 5.1. Geometry and simulation domain

By construction, macroscopic models depend on the geometry considered; they are applicable only for a given geometry; here, we consider the geometry of a thin rectangular channel in a vertical configuration, homogenised over the direction of the small gap $e = 3$ mm. The gravity applies downwards, in the streamwise counter-flow direction. This channel is elongated and thin with a large aspect ratio in the cross-flow plane $l_x/e = 30$ and a large elongation $l_z/l_x \approx 20$ ($l_z = 2$ m). Liquid is injected at a uniform velocity $u_{l,in}$ at the inlet boundary ($z = 0$ m) and vapour is injected with a constant void-fraction $\alpha_{in}$ and constant diameter $d_b$ on a 2 mm wide spanwise centred step of the inlet condition. For simplicity and easier consistency between the two approaches, vapour inlet velocity is prescribed equal to the liquid velocity: $u_{v,in} = u_{l,in}$. The physical mechanism considered is triggered by this inhomogeneous void-fraction injection. Constant pressure is prescribed at the outlet. Other boundaries are adiabatic walls (no-slip in local CFD simulation and slip with a wall-friction model in the volume-averaged DFM). The geometry is given in Fig. 2.

Boundary conditions have been chosen for their simplicity and their ability to trigger various flow regimes but the discontinuity implied by this choice may cause nonphysical behaviour in the lower part of the domain. As the goal of this work is to find general models for real flow configurations (where void fraction gradients are smooth), the lower part of the domain where $z < 0.5$ m ($z/D_h < 80$) is generally not considered in the analysis.

Simulations are performed until steady-state is reached.

## 5.2. Fluid properties

Fluid properties are taken as a mixture of high-pressure water (liquid phase) and steam (gas phase). In the present work, in order to simplify the study and decorrelate phenomenon, the work is conducted in adiabatic conditions with incompressible phases and constant physical properties given in Table 3. They correspond to saturated water at $P = 100$ bar.

<table><caption>**Table 3**
Physical properties of the two phases at $P = 100$ bar.</caption><thead><tr><th>Phase</th><th>Density [kg m<sup>-3</sup>]</th><th>Viscosity [Pa s]</th><th>Surface tension [N m<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Water</td><td>688.424</td><td>8.1718 &times; 10<sup>-5</sup></td><td></td></tr><tr><td>Steam</td><td>55.46</td><td>2.0194 &times; 10<sup>-5</sup></td><td>1.19 &times; 10<sup>-2</sup></td></tr></tbody></table>

## 5.3. Mesh characteristics

As the domain and boundary conditions present symmetries, quarter and half channel simulations with symmetric boundary conditions were investigated in Neptune_CFD. It appears that quarter channel computations are not able to properly catch the flow behaviour due to some unsteady fluctuations in specific conditions (as discussed in Section 6 and seen in Fig. 4). Thus, Neptune_CFD simulations are performed using only (xz) symmetry plane at the middle of the channel ($y = 0$). The CFD mesh aspect is highly elongated in the streamwise direction where the flow is dominated by a strong convection but the elongation remains significantly lower than the maximal code recommendation to prevent numerical instabilities ($\Delta z/\Delta y = 12.5 \ll 100$). In the cross-flow direction, an aspect ratio of 2:1 is used without significant effect on preliminary tests. Mesh convergence is assessed using four levels of refinement ($m_4$ for $f \in \{1, 2, 4, 8\}$) with consecutive doubling of elements in each direction. The mesh convergence analysis shows that the mesh $m_4$ provides a converged solution; it was selected for the analysis and it contains a total of 5 529 600 cells. The non-dimensional wall distance $y^+$ approximately ranges between 20 and 150 depending on the liquid Reynolds number.

For Flica5 simulations, the computational cost is very small and the full domain is simulated. Meshes are given in Table 4 and mesh convergence is achieved with mesh $m_4$ as shown in Section 8. Table 4 also illustrates the tremendous gain in computational resources achieved by the models developed. The number of mesh elements is decreased by a factor of at least 1000.
<figure><img src="image_5.png" /><figcaption>Fig. 3. Numerical cases in the parameter space.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Comparison of the instantaneous void fraction in a quarter channel mesh (left, the simulated domain has been duplicated in the illustration using the symmetry condition at x = 0 m) and a half channel mesh (right) for $\alpha_{in} = 0.10$, $Re_l = 15 000$, $Eo = 0.05$. Direction x is enlarged to improve the figure readability.</figcaption></figure>

<table><caption>**Table 4**

Meshes characteristics for simulations in Neptune_CFD (local TFM) and in Flica5 (averaged DFM).</caption><thead><tr><th>Model</th><th>Dim.</th><th>N<sub>x</sub></th><th>N<sub>y</sub></th><th>N<sub>z</sub></th><th>f</th><th>N<sub>tot</sub><sup>f=8</sup></th></tr></thead><tbody><tr><td>TFM</td><td>3</td><td>180f</td><td>3f</td><td>160f</td><td>1, 2, 4, 8</td><td>44.2Mons</td></tr><tr><td>ave. DFM</td><td>2</td><td>30f</td><td>1</td><td>25f</td><td>1, 2, 4, 8</td><td>48 000</td></tr></tbody></table>

## 5.4. Physical parameter space: $\alpha_{in}$, $Eo$, $Re_l$

In order to explore a wide range of flow conditions, 3 global parameters are selected to account for specific characteristics:

* the Eötvös number $Eo$ defined by Eq. (18a) measures the deformability of the bubbles in the flow by comparing the buoyancy of the bubbles to the surface tension force. It is a key parameter governing the void-fraction distribution in the thin gap direction $y$ by its effect on migration forces (lift) in the CFD simulations. More deformable bubbles (with a larger Eötvös number) tend to migrate more towards the core-flow.

$$Eo = \frac{(\rho_l - \rho_v) g d_b^2}{\sigma} \qquad (18a)$$

* the liquid phase Reynolds number $Re_l$ based on the injection velocity $u_{l,in}$ characterises the turbulence level

$$Re_l = \frac{\rho_l u_{l,in} D_h}{\mu_l} \qquad (18b)$$

* the injection void fraction $\alpha_{in}$ in the central region determines the buoyancy strength and the bubble-induced turbulent intensity.

The combined effect of Eötvös and Reynolds numbers on the transition from wall-peaked to core-peaked flows in adiabatic conditions is discussed in Kommajosyula et al. (2019). It will be determinant on the cross-correlations relevant for the drift-velocity modelling. Here, the inlet boundary condition on diameter $d_b$ and velocity $u_{l,in}$ controls the Eötvös and Reynolds numbers, respectively.

The 3D parameter space ($\alpha_{in}$, $Eo$, $Re_l$) is regularly paved by four values for each dimension, leading to a total of 64 conditions of simulations:

$$Eo \in \{0.05, 0.1, 0.2, 0.4\}, \qquad (19a)$$

$$Re_l \in \{15 000, 30 000, 70 000, 150 000\}, \qquad (19b)$$

$$\text{and } \alpha_{in} \in \{0.05, 0.10, 0.15, 0.20\}. \qquad (19c)$$

A representation of the parameter space is given in Fig. 3. The limit defined by the curved surface corresponds to the flow transition criterion observed by Kommajosyula et al. (2019), mostly for atmospheric air/water conditions. All simulations are defined below the limit to transition to core-peaked conditions because we have demonstrated that TFM and their closures have difficulties to reliably predict core-peaked conditions (Bois et al., 2024). In particular, the model currently implemented in Neptune_CFD always predict wall-peaked flows. In these conditions, dispersed bubbles are weakly deformed.

Numerical experiments are efficient to cover effectively the parameter space (including experimentally-challenging industrial conditions such as for instance high-pressure flows). They provide a way to control
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Mesh convergence of TFM Neptune_CFD simulation for case $α_{in} = 0.10$, $Re_l = 150 000$, $Eo = 0.05$ at $z = 1.5$ m. (a) Void fraction at $y/e ≈ 0.14$. (b) Local drift velocity $ρ_m \bar{u}_D$ at $x = 0.01$ m.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Mesh convergence of TFM Neptune_CFD simulation for case $α_{in} = 0.10$, $Re_l = 150 000$, $Eo = 0.40$ at $z = 1.5$ m. (a) Void fraction at $y/e ≈ 0.14$. (b) Local drift velocity $ρ_m \bar{u}_D$ at $x = 0.01$ m.</figcaption></figure>

separately the physical parameters relevant to the modelled phenomena and they offer the possibility to limit and control the mechanisms involved.

## 6. Physical analysis of local TFM simulations

This section briefly discusses the local TFM results.

All simulations were performed on a half domain in the gap direction (only for $y \in [0 : e/2]$). It is not possible to use the symmetry with the (yz) plane present in the boundary conditions because the solution for some cases demonstrated unsteadiness and oscillatory behaviour (see Fig. 4). If (yz) symmetry plane were used, these oscillations could not develop. Results on mesh m4 defined in Table 4 are converged as illustrated for instance on the void fraction on Fig. 5(a) or on the local drift velocity $ρ_m \bar{u}_D$ (computed from Eq. (9)) by volume-averaging local TFM quantities on Fig. 5(b) for case $α_{in} = 0.10$, $Re_l = 150 000$, $Eo = 0.05$. Results for the largest Eötvös number $Eo = 0.40$ are also well converged on mesh m4 (Figs. 6(a) and 6(b)).

Comparing these two cases reveals that when the Eötvös number increases, the void fraction gets more concentrated in the core of the channel where lateral (x-direction) drift velocities are stronger; hence, on average, the void drift is increased and the void fraction gets more spread in the lateral direction when the Eötvös number is increased (see Fig. 7).

Even though most conditions converged to a steady-state, some conditions below the dashed lines in Fig. 8 lead to unsteady behaviour and the development of oscillations in the void-fraction distribution. These oscillations are visible only after a distance from the inlet, which varies with the flow conditions. They are the strongest at the lowest Reynolds number, probably because they would require a much longer domain to develop significantly at higher Reynolds numbers. Once time-averaged, the mean void-fraction distribution is illustrated

for two unsteady cases on Fig. 10. It can be compared to stable cases on Fig. 9 at either larger Reynolds or Eötvös numbers. When the void-fraction is increased, the instability is triggered more easily and more strongly, thus leading to a wider region affected in the (Re_l, Eo) plane. This destabilising effect is probably related to the stronger buoyancy force leading to stronger velocity gradient similar to chimney effects in natural convection problems. It could be expected as similar oscillations were observed experimentally in mixing zones with void fraction gradients (Auletta and Manca, 2002; Muilwijk and Van den Akker, 2021; De Tournemine and Roig, 2010), in half-sparged configurations. Vortex roll-up are observed in the transition region. The nature of the primary instability of a bubbly mixing layer is associated to two primary mechanisms by De Tournemine and Roig (2010): the buoyancy and a transverse mass transfer due to liquid miscibility in the mixing layer. The induced buoyancy effect generates longitudinal accelerations, which are at the origin of self-excited large-scale oscillations under certain conditions. However, Muilwijk and Van den Akker observed large buoyancy driven vortex roll-up structures in some conditions where De Tournemine and Roig reported steady boundaries. They suspect that the higher liquid velocity and larger bubbles in their case could trigger flow instabilities. In our simulations, the confinement (ratio of bubble diameter to gap size $D_b/e$) and the liquid Reynolds number are larger; they can introduce other mechanisms. Besides, we found that the Eötvös number increase stabilises the flow but the reason remains unclear as it corresponds to more deformable bubbles which tend to be less confined into wall boundary layers and more mobile. This effect is not discussed in the literature to the best of our knowledge; further analyses are necessary to understand the mechanisms involved.

From the averaged prospective, the instability conducts to much stronger void-drift velocity. Besides, this transition is seen to be brutal and non-linear. The triggering might depend on the distance from
<figure><img src="image_9.png" /><figcaption>Fig. 7. Mesh convergence of TFM Neptune_CFD simulation for cases $\alpha_{in} = 0.10$, $Re_l = 150000$, $Eo = 0.05$ and $Eo = 0.40$. y-averaged void fraction $\langle \alpha \rangle$ profile at $z = 1.5$ m. (a) $Eo = 0.05$. (b) $Eo = 0.40$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Flow oscillations in the parameter space. Dashed-lines estimate the unstable region that is expending as the void-fraction increases.</figcaption></figure>

the injection with a characteristic time to develop the oscillations. Experimental confirmation is required to push the analysis further and to comfort the simulation results.<fnref n="1" /> For these reasons, we decided to limit the up-scaling analysis to the Reynolds numbers $Re_l = 70000$ and $150000$ where the flow is always steady. Therefore, we assume that the validity domain of the underlying TFM is restricted to advection-dominated wall-peaked flows, in agreement with the recent work in Bois et al. (2024) for high-pressure boiling flows.

## 7. Volume-averaging and a priori analysis: void drift models' assessment

We consider local TFM simulations as a reference in the restricted sub-space of physical parameters considered, that is we assume the local models to be capable of accurate predictions of mechanisms involved such as cross-flow void fraction distribution. Unfortunately, we are lacking adequate relevant experimental data to directly support the models derived and we fully rely on the global validation of Neptune_CFD code (Mimouni et al., 2010a,b; Laviéville et al., 2015; Marfaing et al., 2017; Mimouni et al., 2017; Mérigoux et al., 2017; Mérigoux, 2022).

To inform averaged models for void drift in volume-averaged DFM, the theoretical derivation has proved that macroscopic drift velocity

$\bar{u}_D$ has two parts: one from the closure of the two-fluid model itself and the other provided by variations of local two-fluid variables. These two parts are accessed by post-processing as illustrated in Step 1 on Fig. 11. By clarifying the exact void drift definition, we provide a better understanding of the complex evolution of the relative velocity and void-fraction distribution in the three-dimensional case. Then, several formulations to model void drift can be assessed a priori and calibrated (Step 2 on Fig. 11). This step can be simply a minimisation of the difference between the reference void-drift and a model by recalibration of coefficients or a more general expression can be sought for. Section 7.1.2 demonstrates the benefits of a coefficient based on the void-fraction itself. Finally, recalibrated models or new models have to be verified in actual volume-averaged DFM resolution (Step 3 on Fig. 11). This is presented in Section 8.

### 7.1. Up-scaling models for filtered DFM: void drift closure

Models to close the subgrid term $\bar{u}_D$ (or equivalently $\bar{u}_R$ as $\bar{u}_D = c_M(1-c_M)\bar{u}_R$) are assessed a priori, i.e., purely based on local TFM simulations. Neptune_CFD data is used at the volume-averaged scale to fit the reference drift velocity (post-processed from Neptune_CFD) and its divergence using only macroscopic variables.

Because of the apparition of oscillations and unsteady behaviour on certain cases, the complete simulation space made of 64 simulations is restricted to 32 simulations (limited to Reynolds numbers of 70000 and 150000). Models based on Eq. (17) and defined in Table 2 are considered. In the next sections, we evaluate separately the abilities of models based on $K_0$ or on a third order polynomial in $\alpha_M$ with a single coefficient $A_d$.

Note that some random noise is present in the reference averaged results. All analyses are performed on regularised data using Savitzky-Golay filters (see Appendix C).

#### 7.1.1. Models based on coefficient $K_0$ (model 1 and 2 in Table 2)

Firstly, we search for the best value of $K_0$ to approximate $\rho_M \bar{u}_{Dx}$ by $-K_0 \partial_x c_M$. Here, $K_0$ is the only non-zero parameter kept in Eq. (17). For each simulation defined by a triplet $(\alpha_{in}, Eo, Re_l)$, the value of $K_0$ is determined to minimise

$$f: K_0 (z_{fit}) \rightarrow \int_{-l_x/2}^{l_x/2} [\rho_M \bar{u}_{Dx}(x, z_{fit}) + K_0 \partial_x c_M] dx. \quad (20)$$

at a given elevation $z = z_{fit}$. This leads to a correlation based on 4 parameters $K_0$ ($z_{fit}$, $\alpha_{in}$, $Eo$, $Re_l$). Fig. 12 illustrates the evolution of this calibration as a function of the optimisation height $z_{fit}$ depending on the case considered. Black markers indicate the reference case used as an example for calibration (a priori) and verification (a posteriori, Section 8), also presented in Fig. 13.

The inlet void-fraction injection is uniform in the gap direction (y) and the vapour has no relative velocity at the inlet ($u_{v,in} = u_{f,in}$).

<footnote n="1">¹ An experimental program is undertaken at CEA to support our simulations and to increase our understanding of these unsteady conditions.</footnote>
<figure><img src="image_11.png" /><figcaption>Fig. 9. Flow development in stable conditions. (a) $\alpha_{in} = 0.10$, $Re_l = 15000$, $Eo = 0.20$. (b) $\alpha_{in} = 0.10$, $Re_l = 15000$, $Eo = 0.10$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Cases with unsteady flow behaviour. (a) $\alpha_{in} = 0.10$, $Re_l = 15000$, $Eo = 0.05$. (b) $\alpha_{in} = 0.10$, $Re_l = 15000$, $Eo = 0.05$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Description of data analysis *a priori* (Step 2) and *a posteriori* (Step 3).</figcaption></figure>

Consequently, the flow needs some time to reach an equilibrium void distribution in this direction. At $z = 0.5 \text{ m} \approx 80D_h$, this distribution is considered established and we assume independence to the inlet boundary condition to perform the analysis. After that, the variability of $K_0$ along the streamwise direction is still significant and the

value of $K_0$ keeps decreasing along the channel height. Predictions are significantly affected by the height of the optimisation for the model, whether it is optimised at the height $z_{fit} = 0.5 \text{ m}$ (noted $K_{05} = K_0 (z_{fit} = 0.5 \text{ m}, \alpha_{in} = 0.1, Eo = 0.1, Re_l = 1.5 \times 10^5)$), or at the height $z_{fit} = 1.5 \text{ m}$ (noted $K_{15}$), or using the averaged value over
<figure><img src="image_14.png" /><figcaption>Fig. 12. Evolution of $K_0$ coefficients along the channel height for 16 cases. Intermediate cases ($Eo = 0.1$ and 0.2) follow the same trends (not shown) Line-styles from dotted, dash-dotted, dashed to line for increasing $\alpha_{in} \in [0.05, 0.1, 0.15, 0.2]$. Symbols + mark the example case used to illustrate calibration and validation.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13. Profiles of $\rho_M \bar{u}_{Dx}$ and its lateral derivative at two elevations $z = 0.5$ and 1.5 m. Models based on $K$ (×) or on $A_d$ (o) calibrated at $z_{fit} = 0.5$ or 1.5 m are compared to the reference (continuous black line). Evolution of $\rho_M \bar{u}_{Dx}$ is anti-symmetric in $x$ while its derivative $\partial_x(\rho_M \bar{u}_{Dx})$ is symmetric. Case $\alpha_{in} = 0.10$, $Re_l = 150000$, $Eo = 0.10$. (a) $\rho_M \bar{u}_{Dx}$ at $z = 0.5$ m. (b) $\rho_M \bar{u}_{Dx}$ at $z = 1.5$ m. (c) $\partial_x(\rho_M \bar{u}_{Dx})$ at $z = 0.5$ m. (d) $\partial_x(\rho_M \bar{u}_{Dx})$ at $z = 1.5$ m.</figcaption></figure>

the sub-domain with $z \in [0.5; 1.5]$ m (noted $K_{05-15}$). These effects are notable both on models for $\rho_M \bar{u}_{Dx}$ or on their lateral derivative $\partial_x(\rho_M \bar{u}_{Dx})$ in Fig. 13. A uniform intermediate value leads to over-diffusion in the lower part and under-diffusion in the higher part of the channel. This sensitivity is illustrated on Fig. 13 with the curved marked with crosses for this particular case. This figure also reveals that the location of maximal drift velocity prediction is slightly shifted. The model also diffuses too strongly at the edges of the injected void fraction step.

There are two main issues related to the non-locality of the model to predict $K_0$ that reduces the applicability of this calibration. Indeed, the calibration in Eq. (20) depends on global parameters associated to the inlet boundary and on the distance to it. This makes the model very case-dependent and inapplicable to other configurations. A more practical model should be calibrated on local flow properties. To achieve that, we consider $z_{fit}$ set to $z_{fit} = 1.5$ m and obtain a uniform model (Model 1); its value does not depend on the position but the calibration is still based on the boundary condition $\alpha_{in}$ and $Re_l$. It is assessed a

posteriori in Section 8 with the value $K_{15}(\alpha_{in} = 0.10, Re_l = 150000, Eo = 0.10) = 1.95 \times 10^{-2} \text{ kg m}^{-2} \text{ s}^{-1}$.

### 7.1.2. Models based on coefficient $A_d$ (model 3 in Table 2)

To alleviate the limitations mentioned above, an other formulation based on a void fraction third order polynomial with a single constant coefficient $A_d$ is evaluated. Here, $A_d$ is the only non-zero parameter kept in Eq. (17). The minimisation is such that:

$$ 
\begin{aligned}
& \int_{-l_x/2}^{l_x/2} \rho_M \bar{u}_D(x, z_{fit}) \bar{e}_x \, dx \\
& = - \int_{-l_x/2}^{l_x/2} A_d \alpha_M(x, z_{fit}) (1 - \alpha_M(x, z_{fit}))^2 \nabla c_M \bar{e}_x \, dx
\end{aligned} 
\quad (21) 
$$

As previously, the coefficient $A_d$ depends a priori on the elevation considered for the fit $z_{fit}$ and on the case boundary conditions ($\alpha_{in}$, $Eo$, $Re_l$). But this coefficient is less sensitive to the height of the calibration in the upper part of the channel where the void-fraction has reached an equilibrium distribution in the gap direction (see Fig. 14 $z > 0.5$ m). More importantly, it is not sensitive to the void fraction at the inlet and
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Evolution of $A_d$ coefficients along the channel height for 16 cases. Intermediate cases ($Eo = 0.1$ and 0.2) follow the same trends (not shown) Line-styles from dotted, dash-dotted, dashed to line for increasing $\alpha_{in} \in [0.05, 0.1, 0.15, 0.2]$. Symbols + mark the example case used to illustrate calibration and validation.</figcaption></figure>

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

Mean value and standard deviation of $A_d$ coefficient. Averages and standard deviations are computed over the different values of $\alpha_{in}$ and for $z \ge 0.5$ m.

<table><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="4">Eo</th></tr><tr><th>0.05</th><th>0.1</th><th>0.2</th><th>0.4</th></tr></thead><tbody><tr><td rowspan="2">Re</td><td>70k</td><td>0.081 ± 0.019</td><td>0.173 ± 0.019</td><td>0.309 ± 0.018</td><td>0.512 ± 0.013</td></tr><tr><td>150k</td><td>0.291 ± 0.024</td><td>0.488 ± 0.033</td><td>0.719 ± 0.082</td><td>1.003 ± 0.052</td></tr></tbody></table>

the whole dependency to local parameters is covered by the local void fraction polynomial in Eq. (17).

Besides, compared to the previous model based on $K_0$ coefficient (Section 7.1.1), the new proposal leads to very satisfactory local predictions both in terms of direct prediction of $\rho_M \bar{u}_D$ or in terms of its divergence; the accuracy remains very good throughout the channel for different elevations (Fig. 13). Thus, averaging the optimised values of $A_d$ for $z \in [0.5 : 1.5]$ m does not introduce significant additional approximation; the $A_d$ parameter which depends on height, injected void fraction, liquid Reynolds number and Eötvös number can be averaged along the vertical axis in order to assess an optimal value for each numerical cases. It is also clear from Fig. 14 that the optimal $A_d$ is very similar for the four conditions at different $\alpha_{in}$ but with the same Reynolds and Eötvös numbers. This property is important for the model to be local, based on local values of void fraction only. Thus, $A_d(\alpha_{in}, Re_l, Eo)$ can be further reduced to $\bar{A}_d^{\alpha_{in}}(Re_l, Eo)$ by averaging it over several values of $\alpha_{in}$. For the case ($\alpha_{in} = 0.10, Re_l = 150000, Eo = 0.10$), the coefficient $A_{15} = 0.445 \text{ kg m}^{-2} \text{ s}^{-1}$ is obtained.

Averaging the optimised values for the four $\alpha_{in}$ and for $z \ge 0.5$ m gives the reference optimised coefficients in Table 5 with the small corresponding standard deviations. On Fig. 15, a clear behaviour can be observed and an analytical function based on $Re_l$ and $Eo$ with second order polynomials is fitted:

$$ \bar{A}_d^{\alpha_{in}} = \sum_{i \le 2} \sum_{j \le 2} C_{ij} (Re/10^5)^i Eo^j \quad (22) $$

For this work, due to the low discretisation of parameter space, the polynomial is in excellent agreement with the data (as seen in Fig. 15). Thus, the coefficients (summarised in Table 6) are limited to the current ranges of Reynolds $Re_l$ and Eötvös $Eo$ numbers; we assume the validity range to be restricted to advection-dominated wall-peaked flows to cope with limitations of current TFM closures. However, with more data and experimental support, this methodology could easily be applied to determine a model for $\bar{A}_d^{\alpha_{in}}$ with a broader validity range. The agreement of the fit with the reference values does not introduce additional error (see Fig. 15). This model calibration is referred to as Model 3 in the following a posteriori assessment.

<caption>**Table 6**</caption>

Fitted polynomials' coefficients $C_{ij}$ for $\bar{A}_d^{\alpha_{in}}$.

<table><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="3">j</th></tr><tr><th>2</th><th>1</th><th>0</th></tr></thead><tbody><tr><td rowspan="3">i</td><td>2</td><td>0</td><td>0.7798</td><td>0.05862</td></tr><tr><td>1</td><td>-3.45</td><td>0.8042</td><td>0.02664</td></tr><tr><td>0</td><td>0.9355</td><td>0.9423</td><td>-0.0536</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** $\bar{A}_d^{\alpha_{in}}$ versus Eötvös number for $Re_l = 70k$ (blue) and $Re_l = 150k$ (red) with the proposed fit (dashed line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 7.2. A priori verification

The distribution of the error made with the prediction of the third model is illustrated a priori on the mass flow term $\rho_M \bar{u}_{Dx}$ with the calibrated model for $\bar{A}_d^{\alpha_{in}}$ given by Eq. (22).

The Weighted Error (WE) defined in Appendix D is shown on Fig. D.22.

Maximum errors are restricted to the lower part of the channel due to the boundary conditions. This justifies to consider the error above 0.5 m to evaluate the model and not the one introduced by the discontinuous boundary condition. For the case presented, the WE remains under 20% and is localised on both sides of the void fraction zone (Fig. D.22(b)).

In summary, a general a priori analysis of the models has been conducted. It showed that models based on $K_0$ coefficient are not practical because the effect of the injected void fraction cannot be eliminated; they are unusable in practical application because of this
<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Void fraction distribution in the cross flow direction (x) for $\alpha_{in} = 0.10$, $Re_l = 150000$, $Eo = 0.10$ with $C_\infty = C_\infty^{opt}$ and the 3rd model for the transverse diffusion ($A_d = 0.466 \text{ kg m}^{-2} \text{ s}^{-1}$). (a) $z = 0.5$ m. (b) $z = 1.5$ m.</figcaption></figure>

non-locality; Besides, these models are less accurate than those based on $A_d$ coefficient with a local void fraction polynomial. More importantly, $A_d$-models can be reduced to depend only on local quantities for the prediction of the drift velocity. A very accurate expression of $A_d$ as a function of $Re_l$ and $Eo$ was found (the dispersion of $A_d$ reference values to the polynomial fit is below 4%). Very good predictive capabilities for void drift are expected from a priori analysis.

## 8. DFM a posteriori verification

This section assesses the effect of different closures for void-drift modelling in the volume-averaged DFM when the equation system (7) is resolved with the closure presented in Section 4.3. Several closures given in Table 2 are evaluated.

### 8.1. Mesh resolution

The void-drift model is intended to represent the physical effects of the distribution in the gap direction. Therefore, it has to be evaluated on converged solutions of the global system, relying on sufficiently fine meshes. Indeed, the calibration proposed is not meant to account for numerical diffusion due to coarse resolution of the discretised equations. Therefore, a mesh convergence analysis is carried out to determine the optimal mesh size. The 4 meshes used are presented in Table 4. Note that a single cell is always used in the gap-wise direction (y), in agreement with the modelling hypotheses. Discretisation of DFM and spatial averaging affect different directions. As a consequence, the mesh should be determined to capture all scales in (x, z) directions as the spatial averaging only reduces the mesh dimension in y direction (to a uniform field in that direction).

The results of mesh convergence are illustrated on Fig. 16 on the macroscopic void fraction $\alpha_M$ distribution at two elevations ($z = 0.5$ m and $z = 1.5$ m) for case $\alpha_{in} = 0.1$, $Re_l = 150000$ and $Eo = 0.1$. They demonstrate that mesh $m_4$ is sufficient to provide converged results so that mesh $m_4$ is used in the following analysis. On the example setup used in this work, the discontinuous inlet leads to a sharp requirement for mesh convergence which would be gentler on continuous physical cases.

### 8.2. Models comparison

The comparison is performed between averaged TFM results and Flica5 predictions with different models, on several sets of input parameters. Comparisons are done on different variables of interest (notably the void fraction distribution).

First, before assessing the void-drift models, we demonstrate the predominant role of the relative velocity prediction in the streamwise direction. Its effect on void-fraction profiles is such that it is compulsory

to adjust this model before performing the comparison of void-drift models. Once this model is adequately calibrated, the capabilities of void-drift models are assessed.

#### 8.2.1. Adaptation of Ishii's correlation for relative velocity prediction

The macroscale simulation is very sensitive to the injection boundary condition combined with the model for relative velocity, as shown on Fig. 17. Indeed, the value of $C_\infty$ affects the void fraction to reach the same vapour mass influx. Besides, some strong lateral drift appears close to the entrance, near ($x = 0.01$ m), because of this discontinuity at the inlet. It is not an important effect in physical continuous conditions.

Predictions with the standard and widely used model from (Ishii, 1977) differ significantly from the averaged TFM results in the conditions investigated. The streamwise relative velocity is not accurate compared to Neptune_CFD reference data, very close to the inlet (i.e., on the first row of cells of the domain in Fig. 17(d)) but also in the whole domain (see Fig. 17(f) at $z = 1.5$ m). Because of mass conservation, this misprediction causes a systematic deviation (bias) on the void fraction. In order to solve this issue, the distribution coefficient $C_\infty$ is calibrated to minimise the deviation between Neptune_CFD relative velocity and an a priori assessment of Ishii's correlation based on Neptune_CFD data. This leads to a recalibration of the $C_\infty$ coefficient from the initial value 1.35 for rectangular channels to a value $C_\infty^{opt}$ slightly below unity. With the optimised coefficient, Fig. 17 shows that a correct injection can be achieved with a good control over the void fraction level. Note that the value of $C_\infty^{opt}$ is calibrated for a single case here, with $\alpha_{in} = 0.10$, $Re_l = 150000$ and $Eo = 0.10$, where the void fraction is maximum near the wall. It cannot be considered valid in a general scope. As an indication, our flows are described as concave void-fraction profiles observed in the developing two-phase flow region in Ishii (1977) and Ishii and Hibiki (2010) for boiling flows. For dynamical aspects, an adiabatic wall-peaked flow can be compared to a subcooled boiling to some extent. For boiling flows, Ishii and Hibiki observe that the maximal distribution value is reached in the developed region when the void fraction increases and the profile is convex. On Fig. 14.6 of Ishii and Hibiki (2010), the data scattering in the developing flow due to boiling is very large in the region where $\alpha < 0.2$. The calibration used here for $C_\infty^{opt}$ falls within the points scattering ($\alpha, C_0$) = (0.1, 0.635). A more appropriate comparison can be observed on Fig. 7 of Hibiki and Ishii (2002) for adiabatic flows in a vertical pipe or in Rassame and Hibiki (2024) for rectangular channels. Larger values of $C_\infty$ are observed experimentally (above 0.9 for Hibiki and Ishii (2002) and close to 1.35 in Rassame and Hibiki (2024)). The difference can be partly attributed to the fixed bubble diameter used in our simulations. The current calibration of $C_\infty^{opt}$ is not meant for general use.
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Void fraction distribution and streamwise relative velocity for Neptune_CFD and Flica5 with default ($C_\infty = 1.35$) or adjusted ($C_\infty = C_\infty^{opt}$) coefficient for Ishii's correlation and the 3rd model for the transverse diffusion ($A_d = 0.466 \text{ kg m}^{-2} \text{ s}^{-1}$). (a) $\alpha_M$ at $z = 0$ m. (b) $\bar{u}_{Rz}$ at $z = 0.1$ m. (c) $\alpha_M$ at $z = 0.1$ m. (d) $\bar{u}_{Rz}$ at $z = 0.1$ m. (e) $\alpha_M$ at $z = 1.5$ m. (f) $\bar{u}_{Rz}$ at $z = 1.5$ m.</figcaption></figure>

### 8.2.2. Void-drift assessment

Now, the transverse drift velocity and its derivative are analysed. This analysis is performed on DFM Flica5 results where the streamwise relative velocity is modelled by Eq. (B.4) for which $C_0$ comes from the combination of Eqs. (B.5a) and (B.6a) with $C_\infty = C_\infty^{opt}$. With the proper coefficient for the streamwise relative velocity, the transverse component of the drift velocity $\rho_M \bar{u}_{Dx} = -K_c \partial c_M / \partial x$ is a relevant indicator to determine the best model. Fig. 18 represents the comparison between Neptune_CFD reference and Flica5 simulations. The reference transverse drift velocity and the equivalent modelled contribution in Flica5 processed from the three different models are compared on Figs. 18(c) and 18(d) for elevations $z = 0.5$ m and $z = 1.5$ m respectively. It is clear that the historical model (based on the local Reynolds number) is inaccurate compared to Neptune_CFD's reference as the transverse drift velocity is over-estimated. It does not match the reference because void-drift in two-phase flow cannot be characterised only by turbulent effects, and thus the Reynolds number. Previous calibration of $K_0$ and $C_M$ parameters did not rely on sufficient configurations to be adapted to the current conditions.

The first (global $K_c = K_0 = 1.95 \times 10^{-2} \text{ kg m}^{-2} \text{ s}^{-1}$) and third (local $K_c$ with the local void fraction dependency given by the polynomial with

$A_d = 0.466 \text{ kg m}^{-2} \text{ s}^{-1}$) models are much more accurate for both elevations presented. But as observed during the a priori analysis (Fig. 13), the 1st model's velocity prediction is best suited for the elevation $z = 1.5$ m and it is under-estimated lower in the domain.

As expected, a good prediction of the drift velocity term leads to a good agreement on the divergence term too. The standard (historical) model is overestimating the divergence contribution in the whole domain, thus leading to over-diffusion. The lateral derivative in the third model presents large discrepancies in the thin region where the void fraction $\alpha_M$ goes to zero, but they are possibly attributed to numerical problems and coupling (in particular, smoothness of the solution when the vapour mass fraction goes to zero); in fact, these differences occur in a region with little impact as the void fraction is already very small there. The first model also has reasonable predictive capabilities but it is overestimating the diffusion on the sides of the peak and underestimating it in the middle.

Regarding the void fraction distribution, the standard model is affected by the over-estimation of diffusion. First and third model both show a good agreement with the reference. The void fraction distribution using the first model has slight differences in shape. In the middle of the domain, the void fraction is less diffused than for
<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Void fraction $\alpha_M$, drift velocity $-K_c \partial_x c_M$ and drift velocity derivative $-\partial_x (K_c \partial_x c_M)$ obtained with different Flica5 models compared with the reference void fraction or transverse component of the drift velocity $\rho_M \bar{u}_{Dx}$ from Neptune_CFD. $\alpha_{in} = 0.10$, $Re_l = 150000$, $Eo = 0.10$. (a) $\alpha_M$ at $z = 0.5$ m. (b) $\alpha_M$ at $z = 1.5$ m. (c) $\bar{u}_{Dx}$ at $z = 0.5$ m. (d) $\bar{u}_{Dx}$ at $z = 1.5$ m. (e) $\partial_x \bar{u}_{Dx}$ at $z = 0.5$ m. (f) $\partial_x \bar{u}_{Dx}$ at $z = 1.5$ m.</figcaption></figure>

Neptune_CFD whereas it is slightly higher on the sides as seen in Fig. 18(b).

With the third model, the void fraction diffusion is in very good agreement with Neptune_CFD's reference void fraction profiles over the whole computational domain. A deviation is observed only on the sides of the injected Heaviside, as expected from the drift velocity and divergence in Section 8.2.2. The modelled divergence with the third model is much sharper than the reference one. This leads to a void fraction profile which does not perfectly diffuse as expected on the sides. This phenomenon is only a minor issue as most of the void fraction profile is perfectly predicted. Improvement in this region might be obtained by a calibration combining $K_0$ and $A_d$ in a formulation depending on both Reynolds number and void fraction.

Finally, the third model is an excellent alternative to the usual model in Flica5; it has demonstrated a good precision in both a priori and a posteriori tests.

Absolute errors between Neptune_CFD and Flica5 on $\alpha_M$ are computed for the 3 different Flica5 models. Fig. 19(c) confirms the third model superiority. As it constrains the error in a tiny band on the edge of the void fraction profile, the first model also shows a good behaviour, but as discussed earlier, it is unusable in real applications because it is

specific to the configuration studied (in particular with the knowledge of $\alpha_{in}$).

## 8.3. Sensitivity analysis on $A_d$

As stated previously, Flica5 does not have any Eötvös number implemented. Thus, the fitted model for $A_d$ coefficient showed in Fig. 15 and in Eq. (22) has a strong dependency to the Eötvös number and therefore to the bubble sizes. This information is scarcely assessed in classical volume-averaged DFM codes. Yet, a strong uncertainty on $A_d$ is introduced in the absence of Eötvös number knowledge as seen in Fig. 15. To evaluate its impact, the sensitivity of the void fraction dispersion is shown on Fig. 20 where 3 values of $A_d$ are compared. A ratio of 2 between two coefficients is considered based on the observed variations of $A_d$ on the considered Eötvös range. Fig. 20 shows that the uncertainty on the Eötvös number can lead to huge over or under estimation of the diffusion coefficient, which then leads to higher or lesser void fraction diffusion throughout the entire domain. The knowledge of the Eötvös number is thus an essential information to control the void fraction distribution in the gap-wise direction. This
<figure><img src="image_21.png" /><figcaption>Fig. 19. Absolute error on the macroscopic void fraction $\alpha_M$ between Flica5 and Neptune_CFD for several transverse drift velocity model. (a) std. model (historical). (b) 1st model. (c) 3rd model.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 20. Void fraction distribution for Neptune_CFD and Flica5 for the 3rd model with $A_d = 0.5 A_d^{ref}$, $A_d = A_d^{ref}$ and $A_d = 2A_d^{ref}$, ($A_d^{ref} = 0.466 \text{ kg m}^{-2} \text{ s}^{-1}$). (a) $z = 0.5$ m. (b) $z = 1.5$ m.</figcaption></figure>

information lost in volume-averaged DFM should be the focus of further modelling efforts.

In summary, the models proposed by calibration from TFM results shows a significant improvement in the predictions of the void fraction spreading in the lateral channel direction. Both approaches (based on either $K_0$ or $A_d$ coefficients) lead to similar results with a good accuracy but the value of $K_0$ in the first model strongly depends on the simulation parameters which probably makes this model less generic and inapplicable in other geometrical configurations. Besides, its variations along the channel's height are less accurate.

On the contrary, the third model with $A_d$ coefficient is very accurate and more general (independent of the height in the channel and of the injected void fraction $\alpha_{in}$). Furthermore, the dependency to the injected Reynolds number can also be reduced in the future by scaling this coefficient on a local Reynolds number. The variation of this coefficient with the Eötvös number is very important and it needs an estimation of the mean bubble diameter. This new model is a major improvement compared to the model currently used in Flica5 simulations.

## 9. Conclusions and prospects

In this paper, we applied an averaging strategy to local two-fluid simulations of upward bubbly flows to derive a closure for the void-drift velocity in volume-averaged drift flux models. The work focuses on highly convective steady bubbly flows, for adiabatic conditions (without phase-change), with wall peaking void fraction profiles. It was conducted as follows:

(i) We demonstrate that the void-drift mechanism quantified by a lateral velocity $\rho_M \bar{u}_D$ is theoretically connected to the local flow characteristics by gap-averaging over the thin direction of the cross-section of rectangular channels.

(ii) We considered 3 dominant parameters ($\alpha$, $E_o$, $Re_l$) and we regularly paved this parameter space whilst keeping configurations where the void fraction is expected to remain mostly trapped in the near-wall regions based on experimental observations of wall-peaked and core-peaked two-phase flows.
(iii) For each condition covering the dimensionless space, local two-fluid model (TFM) is solved and the convection-diffusion of the central void injection is observed.

(iv) These TFM results are analysed. In particular, unsteady flow oscillations occur in some low-Reynolds conditions. These oscillations are developing along the flow and have direct consequences on the void drift which is strongly enhanced. These unsteady conditions are excluded from the drift modelling while experimental support is prepared.

(v) The reference TFM results for stable configurations are averaged and post-processed to determine reference values for the drift velocity $\rho_M \bar{u}_D^{\text{ref}}(\alpha_{in}, E_o, Re_l)$ and for other macroscopic fields.

(vi) Several closures are assessed a priori providing either calibration to existing formulations or suggesting alternate formulations. We actually demonstrate that a diffusion model $-K_c \nabla c_M$ is appropriate provided that the coefficient depends on the void fraction itself under the form $K_c = A_d \alpha_M (1 - \alpha_M)^2$. Then, the coefficient $A_d$ evolved monotonically with the Eötvös and local Reynolds numbers. The dependency of the coefficient to the bubble size (via the Eötvös number) is an important result with critical implications on DFM application to industrial flows. It reveals that a reduced model for the drift must be coupled to a diameter assessment.

(vii) Lastly, the models were also assessed a posteriori with a complete resolution of the volume-averaged drift-flux model and it demonstrated very good capabilities.

The proper definition of the drift velocity based on local fields variations reveals that it has two contributions: one originates from cross-correlations between deviations of mean velocity and vapour mass concentration while the other is based on correlations between vapour mass concentration and the vapour centre-of-volume-velocity with respect to the mixture. We expressed the drift velocity as a function of three dimensionless parameters (void fraction, Reynolds number and Eötvös number) characterising the local conditions at the macroscopic scale. We showed that this new model is compatible with a volume-averaged DFM resolution and that it brings large improvement on void-fraction prediction. This illustrates how to apply the up-scaling methodology on an important mechanism for industrial applications, namely the void fraction dispersion.

The a posteriori validation of the void-drift model also revealed the capital influence of the model for streamwise relative velocity. Predictions with the very common model of (Ishii, 1977) differs significantly from the averaged TFM results in the adiabatic mono-dispersed conditions investigated.

This new methodology to derive industrial models for multiphase flows is very beneficial because it is possible to isolate each mechanism in the modelling and thus prevent model compensation during calibration. This is a new way to access information that would be difficult to obtain experimentally, as long as CFD simulations are trusted. Therefore, we plan to use this approach to assess and eventually revise other closures for the volume-averaged Drift-Flux Model (streamwise relative velocity, friction, mixture momentum diffusion, etc.). This paper also demonstrates that the transverse closure for the drift velocity could be used consistently in any other term involving the transverse relative velocity and in the correction of the macroscopic viscous term.

In addition, an experimental program (named CARIBE, operational in 2025) is launched at CEA to support simulations and to improve the models' validations and capabilities, in particular with respect to advanced variables such as secondary flows or unsteady void-fraction fluctuations that are central elements to the new procedure developed. Future work including TFM validation will support and extend the present analysis, with respect to instabilities, relative velocity and void drift in adiabatic conditions. The experimental program will also provide information on the differences between adiabatic and boiling flows in wall-peaked configurations. This work can also be extended

to other flow conditions such as core-peaked configurations, poly-dispersed and boiling flows, and to different orientations of gravity, channel size, etc. This way, the robustness of the model can be assessed. Also, bubble sizes or deformability is found very influential and should be predicted in volume-averaged DFM solvers but models are very difficult to establish.

### CRediT authorship contribution statement

G. Bois: Writing – review & editing, Writing – original draft, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. R. Barbera: Writing – original draft, Visualization, Validation, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. M. Amate-Vignon: Visualization, Methodology, Investigation, Formal analysis, Data curation, 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.

### Acknowledgements

The authors are grateful to Commissariat à l'Energie Atomique (CEA) for founding support and to some colleagues and in particular Dr. du Cluzeau, Burlot, François and Royer for fruitful discussions and for their careful proofreading of a previous version of this article.

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

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

## Appendix A. Averaging and variable reduction

### A.1. Local statistical averaging

Considering the quantity $\xi$, one defines the phase $k$ statistical average $\bar{\xi}_k^k = \bar{\xi}_k^k / \bar{\chi}_k$ where $\bar{\xi}$ is the statistical average of the quantity $\xi$ (with the classical properties of a Reynolds averaging). The average of the vapour indicator function is called the void fraction $\alpha = \alpha_g = \bar{\chi}_g$. The fields can now be split into a phase averaged field and a fluctuation field : $\xi_k^k = \bar{\xi}_k^k + \xi'_k$. This is an extension of the Reynolds Averaging to multiphase flows, and it has the same properties

$$ \bar{\xi}_k^k = \bar{\xi}_k^k \text{ and } \bar{\xi}_k^f = 0. \qquad (A.1) $$

where the subscript $k \in [l, v]$ stands for the liquid or vapour phase respectively.

### A.2. Volume-averaging — Fluid filtering

Homogenisation is used to describe the two-phase flow at a larger scale. It relies on space-averaging the local system of equations to describe it at the macroscopic scale (Whitaker, 1967). We consider the space filter $\langle \rangle$ of any variable $\xi$ given by

$$ \langle \xi \rangle (x, t) = \frac{1}{V} \int_V \xi(x', t) dV, \qquad (A.2) $$

where the averaging volume $V$ is independent of the position. We also introduce the indicator function of the fluid phase $\chi_f$ equal to unity in the fluid domain and zero otherwise. Besides, walls are associated
to a Dirac delta-distribution $\delta_w$ and their unit normal $\mathbf{n}_w$ (oriented outward from the fluid into the wall) is defined by $\nabla\chi_f = -\delta_w\mathbf{n}_w$. The porosity $\phi = \langle\chi_f\rangle = V_f/V$ gives the ratio of the fluid volume $V_f$ to the averaging volume $V$. Similarly to Whitaker (1999), the intrinsic filtering is defined from the following weighted-average:

$$ \langle\xi\rangle_f(\mathbf{x}, t) = \frac{\langle\xi\rangle_{\chi_f}}{\langle\chi_f\rangle} = \frac{1}{V_f} \int_{V_f} \xi(\mathbf{x}', t) dV. \quad (A.3) $$

Finally, any quantity $\xi$ is decomposed into its filtered quantity $\langle\xi\rangle_f$ and a spatial deviation $\delta_f\xi$: $\xi = \langle\xi\rangle_f + \delta_f\xi$. Contrary to the statistical average, the filter does not commute with space derivatives and it is not idempotent (Whitaker, 1999). For motionless walls, the following rules apply (Howes and Whitaker, 1985):

$$ \phi \langle\nabla\xi\rangle_f = \phi \nabla \langle\xi\rangle_f + \phi \langle\delta_f\xi\delta_w\mathbf{n}\rangle_f, \quad (A.4a) $$

$$ \langle\chi_f \nabla\xi\rangle = \phi \langle\nabla\xi\rangle_f = \nabla(\phi \langle\xi\rangle_f) + \phi \langle\xi\rangle \delta_w \mathbf{n}_f, \quad (A.4b) $$

$$ \phi \left\langle \frac{\partial \xi}{\partial t} \right\rangle_f = \phi \frac{\partial \langle \xi \rangle_f}{\partial t} = \frac{\partial \phi \langle \xi \rangle_f}{\partial t}. \quad (A.4c) $$

The filter is idempotent and the deviation is centred with respect to it if there exists a scale $r_0$ such that

$$ l_{\xi} \lesssim r_0 \lesssim L_{\langle\xi\rangle} $$

where the characteristic sizes $r_0$, $l_\xi$ and $L_{\langle\xi\rangle}$ correspond to respectively the filter kernel, the microscopic scale of variation of the variable $\xi$ and the macroscopic scale of variation of its filtered counterpart $\langle\xi\rangle$.

For space-averaging, a Representative Elementary Volume, REV, or averaging control volume $V$ must be defined. Here, we are interested in converged solutions of the equations given in Section 4. The space-averaging is performed only in the y-direction such that the macroscopic problem is bidimensional in the (x, z) plane.

As a particular example of space-averaging operator $\langle\cdot\rangle$, we consider in this paper a flow in a narrow rectangular channel, with parallel walls at $y = -e/2$ and $y = e/2$. The fluid indicator function is then $\chi_f = 1$ in the fluid domain, i.e., $-e/2 \le y \le e/2$, and $\chi_f = 0$ in the solid structures. In that case, the space-averaging operator is defined over the small channel gap $e$:

$$ \langle\xi\rangle_f(\mathbf{x}, z, t) = \frac{1}{e} \int_{-e/2}^{e/2} \chi_f \xi(x, y, z, t) dy = \frac{1}{e} \int_{-e/2}^{e/2} \xi(x, y, z, t) dy \quad (A.5) \\ = \langle\xi\rangle(\mathbf{x}, y = 0, z, t). $$

After the application of this space-averaging operator, the problem description reduces to two dimensions. Profiles along the channel width are unresolved and sub-grid correlations appear due to non-linearities, products of variables or gradient operators. Due to the dimensionality reduction of this particular filter, the filter is naturally idempotent as one can easily show that

$$ \langle\langle\xi\rangle\rangle_f = \langle\xi\rangle_f \quad (A.6) $$

as a result of the independence of $\langle\xi\rangle_f$ to $y$, which also means that the deviation $\delta_f\xi(x, y, z)$ is centred with respect to that averaging:

$$ \langle\delta_f\xi\rangle_f = 0 \quad (A.7) $$

In addition, this filtering on y-direction naturally commutes with x- and z- derivatives.

### A.3. Favre averaging: density- and fluid-weighted average

As the mixture density strongly varies as a result of variations in vapour concentration, it is more convenient to define a Favre averaging operator weighted by both the mixture density $\rho_m$ and the fluid phase indicator function $\chi_f$, along with the definition of the corresponding deviation:

$$ \langle\xi\rangle_F = \frac{\langle\chi_f \rho_m \xi\rangle}{\langle\chi_f \rho_m\rangle} = \frac{\phi \langle\rho_m \xi\rangle_f}{\phi \langle\rho_m\rangle_f} = \frac{\langle\rho_m \xi\rangle_f}{\rho_M} \quad (A.8a) $$

$$ \delta_F \xi = \xi - \langle\xi\rangle_F \quad (A.8b) $$

$\rho_M = \langle\rho_m\rangle_f$ is the macroscopic mixture density (identified by the capital subscript $M$). From these definitions, we can demonstrate that for two fields $\xi$ and $\zeta$, we have under the assumption of scale separation or for an idempotent filter:

$$ \langle\chi_f \rho_m \xi^\zeta\rangle = \phi \rho_M \langle\xi\rangle_F \langle\zeta\rangle_F + \phi \langle\rho_m \delta\xi \delta\zeta\rangle_F \\ = \phi \rho_M (\langle\xi\rangle_F \langle\zeta\rangle_F + \langle\delta\xi \delta\zeta\rangle_F) \quad (A.9) $$

### A.4. From phase velocities to mixture and relative velocities

Locally, phase and mixture variables are linked by

$$ \rho_m c = \rho_v \alpha, \quad (A.10a) $$

$$ \rho_m = \alpha \rho_v + (1 - \alpha) \rho_l, \quad (A.10b) $$

$$ \bar{\mathbf{u}}_m = c \bar{\mathbf{u}}_v^\nu + (1-c) \bar{\mathbf{u}}_l^\nu, \quad (A.10c) $$

$$ \rho_m \bar{\mathbf{u}}_m = \alpha \rho_v \bar{\mathbf{u}}_v^\nu + (1 - \alpha) \rho_l \bar{\mathbf{u}}_l^\nu, \quad (A.10d) $$

$$ \bar{\mathbf{u}}_r = \bar{\mathbf{u}}_v^\nu - \bar{\mathbf{u}}_l^\nu. \quad (A.10e) $$

where $\rho_m$ is the mixture density and $c$ is the vapour mass fraction. $\bar{\mathbf{u}}_m$ is the centre of mass velocity called barycentric velocity. $c$ is related to the void fraction (noted $\alpha$ instead of $\alpha_v$ to lighten the notations) by the relations

$$ c \rho_m = \alpha \rho_v \quad \text{or} \quad c = \frac{\alpha \rho_v}{\alpha \rho_v + (1 - \alpha) \rho_l} \quad (A.11a) $$

$$ (1-c)\rho_m = (1-\alpha)\rho_l \quad \text{or} \quad 1-c = \frac{(1-\alpha)\rho_l}{\alpha\rho_v + (1-\alpha)\rho_l} \quad (A.11b) $$

and the mixture density can also be expressed as:

$$ \rho_m = \frac{\rho_l \rho_v}{(1 - c) \rho_v + c \rho_l} \quad (A.11c) $$

Conversely, phase velocities can be derived from the mixture and relative velocities by:

$$ \bar{\mathbf{u}}_v^\nu = \bar{\mathbf{u}}_m + (1-c)\bar{\mathbf{u}}_r, \quad (A.11d) $$

$$ \bar{\mathbf{u}}_l^\nu = \bar{\mathbf{u}}_m - c \bar{\mathbf{u}}_r. \quad (A.11e) $$

### A.5. Definition of macroscopic variables

Similarly to the local mixture variables (A.10), filtered quantities weighted by the centre of mass are introduced. Filtered mixture density $\rho_M$, vapour concentration $c_M$, centre of mass velocity $\bar{\mathbf{u}}_M$ and filtered pressure $p_M$ are defined as

$$ \rho_M = \langle\rho_m\rangle_f \quad (A.12a) $$

$$ \rho_M c_M = \langle\rho_m c\rangle_f \quad (A.12b) $$

$$ \rho_M \bar{\mathbf{u}}_M = \langle\rho_m \bar{\mathbf{u}}_m\rangle_f \quad (A.12c) $$

$$ \rho_M = \langle\rho_m\rangle_f \quad (A.12d) $$

by application of Favre filtering (Eq. (A.8a)). Macroscopic quantities are sub-scripted $M$. One can finally note the following property for two variables $\xi$ and $\zeta$:

$$ \langle\rho_m \xi \zeta\rangle_f = \phi \rho_M \langle\xi\rangle_F \langle\zeta\rangle_F + \phi \langle\rho_m \delta\xi \delta\zeta\rangle_f = \phi \rho_M (\langle\xi\rangle_F \langle\zeta\rangle_F + \langle\delta\xi \delta\zeta\rangle_F). \quad (A.13) $$

Besides, macroscopic variables for the phase's velocities are defined. These variables are not used directly in the resolution of the volume-averaged DFM but are needed in the relative velocity models in Appendix B.

From the statistical-averaged local phase velocities $\bar{\mathbf{u}}_l^\nu$ and $\bar{\mathbf{u}}_v^\nu$ and the local vapour mass concentration $c$ and mixture density $\rho_m$, volume
<figure><img src="image_24.png" /><figcaption>Fig. C.21. Noise reduction with Savitzky-Golay filters. Profiles of $\alpha_M$, $\rho_M \bar{u}_{Dx}$, $\partial_x \rho_M \bar{u}_{Dx}$, and their assessment with $A_d$ model at the elevation $z = 1.5$ m. Evolution of $\rho_M \bar{u}_{Dx}$ is anti-symmetric in $x$ while its derivative $\partial_x (\rho_M \bar{u}_{Dx})$ is symmetric. Case $\alpha_{im} = 0.10$, $Re_l = 70000$, $Eo = 0.40$. (a) $\alpha_M$. (b) $\rho_M \bar{u}_{Dx}$. (c) $\partial_x (\rho_M \bar{u}_{Dx})$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. D.22. Weighted error (WE) comparing local $A_d$ variations to model 3 $\bar{A}_d^{a_0} = 0.466 \text{ kg m}^{-2} \text{ s}^{-1}$ from Eq. (22) for case $\alpha_{im} = 0.10$, $Re_l = 150000$, $Eo = 0.10$. The mass flow term is on the left and the divergence on the right of each subfigure. (a) Local $A_d(z)$ (left $\rho_M \bar{u}_{Dx}$, right: $\partial_x(\rho_M \bar{u}_{Dx})$). (b) Model 3: $A_d = 0.466 \text{ kg m}^{-2} \text{ s}^{-1}$ (left $\rho_M \bar{u}_{Dx}$, right: $\partial_x(\rho_M \bar{u}_{Dx})$).</figcaption></figure>

and phase-averaged macroscopic phase velocities $\bar{u}_L$ and $\bar{u}_V$ are such that

$$ \bar{u}_K = \frac{\langle \rho_m c_k \bar{u}_k^k \rangle_f}{\langle \rho_m c_k \rangle_f} = \frac{\langle \alpha_k \bar{u}_k^k \rangle_f}{\langle \alpha_k \rangle_f} \quad \text{for } K \in [L, V] \quad (A.14) $$

because phase densities are considered constant and using Eqs. (A.11a) and (A.11b).

Thus, $\bar{u}_V$ and $\bar{u}_L$ are defined from $\bar{u}_v^v$ and $\bar{u}_l^l$ spatially-pondered by the void fraction while $\bar{u}_M$ is defined from $\bar{u}_m$ spatially-pondered by the mixture density as given by Eq. (A.12c)

$$ \bar{u}_M = c_M \bar{u}_V + (1 - c_M) \bar{u}_L \quad \text{while} \quad \bar{u}_m = c \bar{u}_v^v + (1 - c) \bar{u}_l^l \quad (A.15) $$

and we also get the macroscopic relative velocity $\bar{u}_R = \bar{u}_V - \bar{u}_L$ while the local one is $\bar{u}_R = \bar{u}_v^l - \bar{u}_l^v$.

This formulation for the macroscopic phase velocities also preserve the reversed relations between $\bar{u}_L$, $\bar{u}_V$ and $\bar{u}_M$, $\bar{u}_R$

$$ \bar{u}_L = \bar{u}_M - c_M \bar{u}_R, \quad (A.16a) $$

$$ \text{and } \bar{u}_V = \bar{u}_M + (1 - c_M) \bar{u}_R \quad (A.16b) $$

while the lower scale relations are given by Eqs. (A.11d) and (A.11e).

## Appendix B. Flica5 model for the filtered relative velocity (Ishii's correlation)

This section describes the closure of the streamwise relative velocity $\bar{u}_{Rz}$ in Flica5.

The correlation for the relative velocity is proposed in Ishii (1977). Although only low void fraction is considered in this work (below 40%), Ishii's correlation is able to handle higher void fractions in the annular regime. For cases considered here, the flow regime is always bubbly but in Flica5, Ishii's correlation is built so that there is a continuity between regimes which slightly affects the predictions even in bubbly configurations.

By definition, the mean drift velocity is equal to

$$ \bar{u}_{V \to J} = \frac{\langle \alpha \bar{u}_v^v \rangle_f - \langle j \rangle_f}{\langle \alpha \rangle_f} = \bar{u}_V - J = (1 - \alpha_M) \bar{u}_R \quad (B.1) $$

with $j = j_l + j_v = \alpha_l \bar{u}_l^l + \alpha_v \bar{u}_v^v$. This definition can be reformulated introducing the distribution coefficient $C_0$ the expression

$$ \bar{u}_{V \to J} = \frac{\langle \alpha \bar{u}_{V \to J} \rangle_f}{\langle \alpha \rangle_f} + (C_0 - 1) J \quad \text{with} \quad C_0 = \frac{\langle \alpha j \rangle_f}{\langle \alpha \rangle_f \langle j \rangle_f} \quad (B.2) $$
to finally get

$$ \bar{\mathbf{u}}_V = \frac{\langle j_v \rangle_f}{\langle \alpha \rangle_f} = \frac{\langle \alpha \bar{\mathbf{u}}_{v \to j} \rangle_f}{\langle \alpha \rangle_f} + C_0 \mathbf{J} \quad (B.3) $$

$$ \bar{\mathbf{u}}_R = \frac{(C_0 - 1) \bar{\mathbf{u}}_M + \bar{\mathbf{u}}_{Vlim}}{1 - c_M + (c_M - \alpha_M) C_0} \quad (B.4) $$

where the vapour velocity limit $\bar{\mathbf{u}}_{Vlim} = \frac{\langle \alpha \bar{\mathbf{u}}_{v \to j} \rangle_f}{\langle \alpha \rangle_f}$ is the terminal

vapour velocity when $\mathbf{J} \to 0$. Ishii and Zuber (1979) tested various configuration of multiphase upward flows in pipes to derive correlations for $C_0$ and $\bar{\mathbf{u}}_{Vlim}$. In Flica5, the Ishii's correlation is designed to be a continuous blending between bubbly and annular flows. Even though we are only interested in the bubbly flow regime here, this blending slightly affects the models predictions when $\alpha_M > 0.1$.

### Bubbly flow

$$ C_{0,b} = \left( C_{\infty} + (1 - C_{\infty}) \sqrt{\frac{\rho_v}{\rho_l}} \right) (1 - e^{-\zeta \alpha_M}) \quad (B.5a) $$

$$ \bar{\mathbf{u}}_{Vlim,b} = C_l \left( \frac{(\rho_l - \rho_v) g \sigma}{\rho_l^2} \right)^{1/4} (1 - \alpha_M)^{7/4} \quad (B.5b) $$

with $C_\infty = 1.35$ for rectangular channel and 1.2 for cylindrical pipes and $C_l = 4.5$ for pipe section and $\sqrt{2}$ for other cases.

### Annular flow

$$ C_{0,a} = \frac{1 + \mathbf{B} \sqrt{\frac{\rho_v}{\rho_l}}}{\alpha_M + \mathbf{B} \sqrt{\frac{\rho_v}{\rho_l}}} \quad \text{with} \quad \mathbf{B} = \frac{\sqrt{1 + 75(1 - \alpha_M)}}{\sqrt{\alpha_M}} \quad (B.6a) $$

$$ \bar{\mathbf{u}}_{Vlim,a} = \frac{(1 - \alpha_M)^{3/2}}{\alpha_M + \mathbf{B} \sqrt{\frac{\rho_v}{\rho_l}}} \sqrt{\frac{g D_{h xi} (\rho_l - \rho_v)}{0.015 \rho_l}} \quad (B.6b) $$

These relations provide a prediction for $\bar{u}_{Rz}$ from local values of $\bar{u}_M$, $\alpha_M$, $c_M$ and fluid properties.

## Appendix C. Noise reduction on reference data using Savitzky-Golay filters

Noise reduction has been a major issue in this work and is succinctly presented here. On some Neptune_CFD simulations, local and averaged void-fractions present numerical noise, as illustrated on Fig. C.21(a). As a consequence, direct derivation of $c_M$ and void fraction multiplication to construct the *a priori* assessment of the third model conduct to very noisy data (see Figs. C.21(b) and C.21(c)). To reduce this numerical noise, the averaged void-fraction $\langle \alpha \rangle$ was smoothed using Savitzky-Golay filters (SG) from Savitzky and Golay (1964) implemented in method savgol_filter from Scipy module. A third order polynomial composed of 41 elements is used. The outcome is also illustrated on Fig. C.21(a), with the curve SG($\alpha_M$). After this, a priori assessment is greatly improved, as shown on the drift velocity and its approximation on Fig. C.21(b). Calibration of $A_d$ coefficient can now be fitted properly.

## Appendix D. Error criteria

A weighted error WE is defined to represent the discrepancy between model and reference terms:

$$ WE = \frac{\| NCFD - Model \|}{\max_{z \ge 0.5} (\| NCFD \|)} \times 100 \quad (D.1) $$

The weighted error normalisation is calculated only for a height over 0.5 m.

A global error criteria MWE is built from WE:

$$ MWE = mean_{\alpha_M \ge 0.0001 \alpha_{in}} (WE) \quad (D.2) $$

This criterion shows the mean weighted error for each case, calculated only when the local void fraction $\alpha_M$ is greater than a minimum value $\alpha_M \ge 0.0001 \alpha_{in}$, to eliminate regions where terms vanish.

## Data availability

The authors do not have permission to share data.

## References

Anderson, T.B., Jackson, R., 1967. Fluid mechanical description of fluidized beds. Equations of motion. Ind. Eng. Chem. Fundam. 6 (4), 527–539. http://dx.doi.org/10.1021/i160024a007.

Auletta, A., Manca, O., 2002. Heat and fluid flow resulting from the chimney effect in a symmetrically heated vertical channel with adiabatic extensions. Int. J. Therm. Sci. 41 (12), 1101–1111. http://dx.doi.org/10.1016/S1290-0729(02)01396-0.

Bois, G., 2017. Direct numerical simulation of a turbulent bubbly flow in a vertical channel: Towards an improved second-order Reynolds stress model. Nucl. Eng. Des. 321, 92–103. http://dx.doi.org/10.1016/j.nucengdes.2017.01.023, Multi-scale multi-physics analysis of reactor transients in the NURESAFE project.

Bois, G., 2021. Derivation of a macroscopic mixture model for two-phase turbulent flows. Int. J. Heat Mass Transfer 178, 121500. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2021.121500.

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

Bowring, R.W., 1968. A Computer Programme for the Subchannel Analysis of the Hydraulic and Burnout Characteristics of Rod Clusters. Part 2. The Equations. Technical Report AEEW-R582, Atomic Energy Establishment, Winfrith, URL: https://www.osti.gov/biblio/4486717.

Caremoli, C., Rascle, P., Aubry, S., Olive, S., 1993. THYC, a 3D Thermal-Hydraulic Code for Rod Bundles: Recent Developments and Validation Tests. Technical Report, Électricité de France, Direction des Etudes et Recherche, URL: https://inis.iaea.org/search/search.aspx?orig_q=RN:26001958.

Carlucci, L., Hammouda, N., Rowe, D., 2004. Two-phase turbulent mixing and buoyancy drift in rod bundles. Nucl. Eng. Des. 227 (1), 65–84. http://dx.doi.org/10.1016/j.nucengdes.2003.08.003.

Carver, M., Judd, R., Kiteley, J., Tahir, A., 1987. The drift flux model in the ASSERT subchannel code. Nucl. J. Can. J. Nucl. Can. 1 (2), 153–165, URL: https://canteach.candu.org/Content%20Library/NJC-1-2-08.pdf.

Collier, J.G., 1979. Two-phase gas-liquid flows within rod bundles. In: Institution of Mechanical Engineers, Conference Publications. Vol. 2, pp. 1041–1055, URL: https://www.scopus.com/inward/record.uri?eid=2-s2.0-0018713775&partnerID=40&md5=58d8f4222db2f925591ab5434996df8c, in S. Kaka and D.B. Spalding (eds.), Turbulent Forced Convection in Channels and Bundles, Hemisphere, London.

De Tournemine, A.L., Roig, V., 2010. Self-excited oscillations in buoyant confined bubbly mixing layers. Phys. Fluids 22 (2), 023301. http://dx.doi.org/10.1063/1.3327290, arXiv:https://pubs.aip.org/aip/pof/article-pdf/doi/10.1063/1.3327290/15792425/023301_1_online.pdf.

Delhaye, J.-M., Achard, J.-L., 1976. On the Averaging Operators Introduced in Two-Phase Flow Modeling. Centre d'études nucléaires de Grenoble.

Delhaye, J.-M., Collier, J.-G., Hewitt, G., Bergles, A.E., 1981b. Two-Phase Flow and Heat Transfer in the Process and Power Industries. Hemisphere Publishing Corporation, McGraw-Hill, New York.

Delhaye, J.-M., Giot, M., Riethmuller, M.L., 1981a. Thermohydraulics of Two-Phase Systems for Industrial Design and Nuclear Engineering. Hemisphere Publishing Corporation.

Drew, D.A., Passman, S.L., 1999. Theory of Multicomponent Fluids, springer ed. Springer, http://dx.doi.org/10.1007/b97678.

Gonzalez-Santalo, J., Griffith, P., 1972. Two Phase Flow Mixing in Rod Bundle Subchannels. ASME Paper, 72-WA/NE-19, p. 1, URL: https://www.scopus.com/inward/record.uri?eid=2-s2.0-85024352521&partnerID=40&md5=6f2c1f51905310264b5e449f2a2fe2f.

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

Hibiki, T., Ishii, M., 2002. Distribution parameter and drift velocity of drift-flux model in bubbly flow. Int. J. Heat Mass Transfer 45 (4), 707–721. http://dx.doi.org/10.1016/S0017-9310(01)00195-8.

Hibiki, T., Ishii, M., 2003. One-dimensional drift-flux model and constitutive equations for relative motion between phases in various two-phase flow regimes. Int. J. Heat Mass Transfer 46 (25), 4935–4948. http://dx.doi.org/10.1016/S0017-9310(03)00322-3.

Howes, F.A., Whitaker, S., 1985. The spatial averaging theorem revisited. Chem. Eng. Sci. 40 (10), 1387–1392.
Ishii, M., 1977. One-Dimensional Drift-Flux Model and Constitutive Equations for Relative Motion Between Phases in Various Two-Phase Flow Regimes. Technical Report ANL-77-47, Argonne National Lab., Ill. (USA), http://dx.doi.org/10.2172/6871478.

Ishii, M., Hibiki, T., 2010. Thermo-Fluid Dynamics of Two-Phase Flow, second ed. Springer, New York, NY, pp. XVII-518. http://dx.doi.org/10.1007/978-1-4419-7985-8.

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

Jouhanique, T., Rascle, P., 1995. A Fifth Equation to Model the Relative Velocity in the 3-D Thermal-Hydraulic Code THYC. Technical Report, Direction des études et recherches. Production d'énergie (hydraulique, thermique et nucleaire).

Kataoka, I., 1986. Local instant formulation of two-phase flow. Int. J. Multiph. Flow 12 (5), 745-758. http://dx.doi.org/10.1016/0301-9322(86)90049-2.

Kawahara, A., Sadatomi, M., Higuchi, T., 2008. Effects of surface tension on two-phase void drift between triangle tight lattice subchannels. J. Eng. Gas Turbines Power 131 (1), http://dx.doi.org/10.1115/1.2983083, arXiv:https://asmedigitalcollection.asme.org/gasturbinespower/article-pdf/131/1/012903/5498637/012903_1.pdf, 012903.

Kommajosyula, R., Bois, G., Burlot, A., Rodio, M.-G., Cariteau, B., Baglietto, E., 2019. A criterion to classify void fraction profiles in adiabatic bubbly flow based on averaged flow quantities for use in subchannel codes. In: 18th International Topical Meeting on Nuclear Reactor Thermal Hydraulics, NURETH 2019. Vol. 592-28889, pp. 4720-4731.

Lahey, R.T., Moody, F.J., 1977. The Thermal Hydraulics of a Boiling Water Nuclear Reactor. p. 631.

Lahey, R., Shiralkar, B., Rarcliffe, D., Polomik, E., 1972. Out-of-pile subchannel measurements in a nine-rod bundle for water at 1000 PSIA. In: Hetroni, G., Sideman, S., Hartinett, J. (Eds.), Proceedings of the International Symposium on Two-Phase Systems. Pergamon, pp. 345-363. http://dx.doi.org/10.1016/B978-0-08-017035-0.50023-9.

Lahey, Jr., R.T., 1974. Two-Phase Flow in Boiling Water Nuclear Reactors. Vol. NEOO-13388, General Electric Co., URL: https://www.osti.gov/biblio/4158019.

Laviéville, J., Mérigoux, N., Guingo, M., Baudry, C., Mimouni, S., 2015. A Generalized Turbulent Dispersion Model for bubbly flow numerical simulation in NEPTUNE_CFD. In: NURETH-16. pp. 4167-4181.

Levy, S., 1963. Prediction of two-phase pressure drop and density distribution from mixing length theory. J. Heat Transfer 85 (2), 137-150. http://dx.doi.org/10.1115/1.3686033, arXiv:https://asmedigitalcollection.asme.org/heattransfer/article-pdf/85/2/137/5753561/137_1.pdf.

Marfaing, O., Guingo, M., Laviéville, J.M., Mimouni, S., 2017. Analytical void fraction profile near the walls in low Reynolds number bubbly flows in pipes: Experimental comparison and estimate of the dispersion coefficient. Oil Gas Sci. Technol. 72 (1), http://dx.doi.org/10.2516/ogst/2016027.

McAdams, W.H., McAdams, W.H., 1954. Heat Transmission, vol. 3, McGraw-hill, New York.

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

Mérigoux, N., Laviéville, J., Mimouni, S., Guingo, M., Baudry, C., Bellet, S., 2017. Verification, validation and application of NEPTUNE_CFD to two-phase Pressurized Thermal Shocks. Nucl. Eng. Des. 312, 74-85. http://dx.doi.org/10.1016/j.nucengdes.2016.06.041, 16th International Topical Meeting on Nuclear Reactor Thermal Hydraulics.

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

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

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

Mimouni, S., Guingo, M., Laviéville, J., Mérigoux, N., 2017. Combined evaluation of bubble dynamics, polydispersion model and turbulence modeling for adiabatic two-phase flow. Nucl. Eng. Des. 321, 57-68. http://dx.doi.org/10.1016/j.nucengdes.2017.03.041.

Mimouni, S., Laviéville, J., Seiler, N., Ruyer, P., 2011. Combined evaluation of second order turbulence model and polydispersion model for two-phase boiling flow and application to fuel assembly analysis. Nucl. Eng. Des. 241 (11), 4523-4536.

Moorthi, A., Kumar Sharma, A., Velusamy, K., 2018. A review of sub-channel thermal hydraulic codes for nuclear reactor core and future directions. Nucl. Eng. Des. 332, 329-344. http://dx.doi.org/10.1016/j.nucengdes.2018.03.012.

Morel, C., 2015. Mathematical Modeling of Two-Phase Flow. Springer, http://dx.doi.org/10.1007/978-3-319-20104-7.

Muilwijk, C., Van den Akker, H.E., 2021. The effect of liquid co-flow on gas fractions, bubble velocities and chord lengths in bubbly flows. Part II: Asymmetric flow configurations. Int. J. Multiph. Flow 138, 103562. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2021.103562.

Otero, F., Veloso, M.A., Pereira, C., Fortini, A., Lombardi, A., 2011. Assessment of 4x4 rod bundle subchannel mixing experiments. In: Proceedings of the International Nuclear Atlantic Conference - INAC. Belo Horizonte, MG, Brazil, pp. 1-10, URL: https://inis.iaea.org/search/search.aspx?orig.q=RN:43046460.

Quintard, M., Whitaker, S., 1994. Transport in ordered and disordered porous media II: Generalized volume averaging. Transp. Porous Media 14, 179-206. http://dx.doi.org/10.1007/BF00615200.

Rassame, S., Hibiki, T., 2024. Drift-flux model for dispersed adiabatic and boiling two-phase flows in rectangular channels. Int. J. Heat Mass Transfer 224, 125270. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2024.125270.

Reichardt, H., 1951. Vollständige darstellung der turbulenten geschwindigkeitsverteilung in glatten leitungen. ZAMM - J. Appl. Math. Mech. / Z. Angew. Math. Mech. 31 (7), 208-219. http://dx.doi.org/10.1002/zamm.19510310704, URL: arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/zamm.19510310704.

Rowe, D.S., Angle, C.W., 1969. Crossflow Mixing Between Parallel Flow Channels During Boiling. Part III. Effect of Spacers on Mixing Between Two Channels. Technical Report BNWL-371 PT 3, Battelle Pacific Northwest Labs., Richland, WA (United States), http://dx.doi.org/10.2172/4823782.

Rudzinski, K.F., 1970. Two-Phase Turbulent Mixing for Air-Water Flows in Adjacent Triangular Array Subchannels (Ph.D. thesis). University of Windsor, URL: https://scholar.uwindsor.ca/etd/6626, Electronic Theses and Dissertations.

Sadatomi, M., Kawahara, A., Kano, K., Sumi, Y., 2004. Single- and two-phase turbulent mixing rate between adjacent subchannels in a vertical 2x3 rod array channel. Int. J. Multiph. Flow 30 (5), 481-498. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2004.03.001.

Sadatomi, M., Kawahara, A., Sato, Y., 1994. Flow redistribution due to void drift in two-phase flow in a multiple channel consisting of two subchannels. Nucl. Eng. Des. 148 (2), 463-474. http://dx.doi.org/10.1016/0029-5493(94)90126-0.

Savitzky, A., Golay, M.J.E., 1964. Smoothing and differentiation of data by simplified least squares procedures. Anal. Chem. 36 (8), 1627-1639. http://dx.doi.org/10.1021/ac60214a047, arXiv:https://doi.org/10.1021/ac60214a047.

Sharma, M., Nayak, A., 2016. Experimental investigation of two phase turbulent mixing rate under bubbly flow regime in simulated subchannels of a natural circulation pressure tube type BWR. Exp. Therm Fluid Sci. 76, 228-237. http://dx.doi.org/10.1016/j.experthrmflusci.2016.03.016.

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

Tincq, D., David, F., Rascle, P., 1991. THYC: A 3D Thermal-Hydraulic Code for Steam Generators, Heat Exchangers and Condensers: Recent Developments and Validation Tests. Technical Report, Direction des études et recherches. Production d'énergie (hydraulique, thermique et nucleaire), URL: https://inis.iaea.org/search/search.aspx?orig.q=RN:23063454.

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

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

Wallis, G., 1969. One-Dimensional Two-Phase Flow, McGraw Hill ed. Springer-Verlag, p. 408.

Weisman, J., Bowring, R.W., 1975. Methods for detailed thermal and hydraulic analysis of water-cooled reactors. Nucl. Sci. Eng. 57 (4), 255-276. http://dx.doi.org/10.13182/NSE75-A15419, arXiv:https://doi.org/10.13182/NSE75-A15419.

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

Whitaker, S., 1967. Diffusion and dispersion in porous media. AIChE 3 (13), 420-427.

Whitaker, S., 1969. Advances in theory of fluid motion in porous media. Ind. Eng. Chem. 61 (12), 14-28. http://dx.doi.org/10.1021/ie50720a004.

Whitaker, S., 1999. Theory and Applications of Transport in Porous Media: The Method of Volume Averaging. Kluwer Academic Publishers.

Yeoh, G., Tu, J., 2016. Basic theory and conceptual framework of multiphase flows. pp. 1-47. http://dx.doi.org/10.1007/978-981-4585-86-6 1-1.

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

Zuber, N., Findlay, J., 1965. Average volumetric concentration in two-phase flow systems. J. Heat Transfer 87, 453-468.