*Article* **Numerical Prediction of Two-Phase Flow through a Tube Bundle Based on Reduced-Order Model and a Void Fraction Correlation**

**Claire Dubot 1,2,\*, Cyrille Allery 1, Vincent Melot 2, Claudine Béghein 1, Mourad Oulghelou <sup>3</sup> and Clément Bonneau <sup>2</sup>**


**Abstract:** Predicting the void fraction of a two-phase flow outside of tubes is essential to evaluate the thermohydraulic behaviour in steam generators. Indeed, it determines two-phase mixture properties and affects two-phase mixture velocity, which enable evaluating the pressure drop of the system. The two-fluid model for the numerical simulation of two-phase flows requires interaction laws between phases which are not known and/or reliable for a flow within a tube bundle. Therefore, the mixture model, for which it is easier to implement suitable correlations for tube bundles, is used. Indeed, by expressing the relative velocity as a function of slip, the void fraction model of Feenstra et al. and Hibiki et al. developed for upward cross-flow through horizontal tube bundles is introduced and compared. With the method suggested in this paper, the physical phenomena that occur in tube bundles are taken into consideration. Moreover, the tube bundle is modelled using a porous media approach where the Darcy–Forchheimer term is usually defined by correlations found in the literature. However, for some tube bundle geometries, these correlations are not available. The second goal of the paper is to quickly compute, in quasi-real-time, this term by a non-intrusive parametric reduced model based on Proper Orthogonal Decomposition. This method, named Bi-CITSGM (Bi-Calibrated Interpolation on the Tangent Subspace of the Grassmann Manifold), consists in interpolating the spatial and temporal bases by ITSGM (Interpolation on the Tangent Subspace of the Grassmann Manifold) in order to define the solution for a new parameter. The two developed methods are validated based on the experimental results obtained by Dowlati et al. for a two-phase cross-flow through a horizontal tube bundle.

**Keywords:** steam generator; void fraction; mixture model; porous media approach; reduced-order model; Proper Orthogonal Decomposition (POD)

#### **1. Introduction**

Steam generators are heat exchangers used especially in nuclear propulsion. Water, heated by the reactor core, flows through a tube bundle, which is a closed circuit called the primary circuit. The heat of the primary fluid is diffused by conduction through metallic tube walls to the water, which flows outside the tubes. Water in the secondary circuit, also called the secondary fluid, enters in a liquid state and becomes a two-phase mixture of steam and water as heat transfer occurs along the heat exchanger. The steam is then used to generate electricity using rotating turbines.

A three-dimensional thermo-hydraulic analysis is essential to predict the performance of heat exchangers and their correct design, especially taking into account that the tube bundle, where there may be thousands of tubes, would require unacceptable computational

**Citation:** Dubot, C.; Allery, C.; Melot, V.; Béghein, C.; Oulghelou, M.; Bonneau, C. Numerical Prediction of Two-Phase Flow through a Tube Bundle Based on Reduced-Order Model and a Void Fraction Correlation. *Entropy* **2021**, *23*, 1355. https://doi.org/10.3390/e23101355

Academic Editor: Donald J. Jacobs

Received: 31 August 2021 Accepted: 13 October 2021 Published: 16 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

cost and time. Therefore, the whole code cited in this paper models the tube bundle as a porous medium. In the case of two-phase flows, the void fraction is the key parameter to characterize the flow. Indeed, it enables calculating the mixture density, the mixture viscosity, and the mixture velocity. Consequently, it plays an important role for computing pressure drops and heat and mass transfers.

Since the 1980s, thermal-hydraulic codes have been developed to understand the physical phenomena involved. One of the first codes, THIRST [1], was developed to compute three-dimensional, two-phase, and steady flow in steam generators. The twophase flow was solved by the homogeneous two-phase model, and the phase velocities were assumed to be equal, but this does not reflect what is really going on. To take into account the slip between phases, Navier–Stokes equations were solved for the mixture of the secondary fluid. The THYC code [2] gives the relative velocity thanks to a correlation using the drift flux model of Zuber and Findlay [3]. This model is based on the determination of drift-flux parameters for which there are many different empirical correlations [4–8]. Stevanovic et al. [9] used the two-fluid model to predict the thermal-hydraulic behavior in horizontal tube bundles. Navier–Stokes equations were solved for each phase. The accurate definition of the interfacial drag force, comprising a drag coefficient correlation, is important in order to predict the void fraction distribution. In their paper, the original drag correlation of Ishii and Zuber [10] was multiplied by 0.4. They validated this modification of the correlation with their experimental data. Nevertheless, most drag coefficient laws in the literature [11–15], including Ishii and Zuber's correlation, are made for different two-phase flow regimes (bubbly, slug, stratified, annular, or spray flow) inside a tube and not in tube bundles. In this study, this problem is tackled by using the mixture model, which is a simplified two-phase model where Navier–Stokes equations are solved for the mixture. The developed method involved formulating the relative velocity as a function of slip and then implementing a specific void fraction model for tube bundles derived from the literature.

The slip ratio is defined as the ratio between gas phase velocity and liquid phase velocity. Feenstra et al. [16] developed a slip ratio model based on their R-11 data for upward two-phase cross-flow through horizontal tube bundles. They identified the important variables that affected slip, and the application of the Buckingham pi theorem enabled them to reduce the slip ratio as a function of two dimensionless numbers, namely the Richardson number and the Capillary number. They demonstrated that it fitted well with experimental void fraction data in R-11 and air–water mixtures for a wide range of mass fluxes, qualities and pitch-diameter ratios. This model is not explicit for void fraction; indeed an iterative process is necessary to compute the void fraction. Likewise, Hibiki et al. [17] used the slip ratio correlation of Smith [18] to develop a correlation for the entrainment factor dependent on the tube's bundle arrangement. The entrainment factor is correlated with a dominant parameter such as nondimensional mass flux based on experimental data coming from various flow configurations and tube bundle arrangements. The developed correlation agrees both with parallel and crossflow in the tube bundle.

Moreover, the porous media approach implies adding a momentum sink to the governing momentum equation. This source term is defined by the Darcy–Forchheimer law [19,20]. A widely used and easy method is the use of correlations coming from the literature, such as Zukauskas et al.'s correlation [21] for a transverse flow to a tube bundle or Rhodes and Carlucci's correlation [22] for a parallel flow. These correlations result from experiments, it seems to be the most accurate method to determine the law but it is valid for a given design and it is expensive and time-consuming. However, for new, i.e., non-standard, tube geometry, these correlations are not available, and solving a reduced-order model to define them is suggested. From some CFD calculations of flow through a Representative Elementary Volume (REV) of the tube bundle, pressure and velocity fields were decomposed by POD (Proper Orthogonal Decomposition). Then, the non-intrusive reduced model method, Bi-CITSGM (Bi-Calibrated Interpolation on the Tangent Subspace of the Grassmann Manifold) [23,24], was applied in order to determine the solution for a new

parameter. In this method, spatial and temporal POD sampling bases are interpolated by ITSGM (Interpolation on the Tangent Subspace of the Grassmann Manifold) [25] and the temporal eigenvalues by usual methods such as Lagrange, IDW (Inverse Distance Weighting), or RBF (Radial Basis Function). Then, the reduced-order model was used to compute the source term of the porous media approach applied to the flow through a tube bundle.

This work was validated with the experimental results of Dowlati et al. [26] which is a two-phase cross-flow through a horizontal tube bundle. The summary of the developed methodology is presented in Figure 1. First, the mixture model was used, and we rewrote the slip velocity, *ugl*, as a function of the slip in order to implement a specific void fraction model. The implementation of Feenstra's correlation and Hibiki's correlation were compared to the usual formulation of Manninen et al. [27]. Moreover, the tube bundle was modelled by a porous medium. The ability to use a reduced-order model to determine the momentum sink term *Ft* was studied and compared to the usual correlation of Zukauskas et al. [21].

The remainder of the paper is organized as follows. In Section 2, the governing equations used to model the two-phase flow through a tube bundle by a porous media approach are presented. The rewriting of the relative velocity and the Darcy–Forchheimer term are also detailed. Section 3 deals with the methodology of the reduced-order model on the REV to compute the Forchheimer term. The ITSGM method and the non-intrusive approach Bi-CITSGM are reviewed. Section 4 validates the application of the proposed reformulation of the relative velocity, and Section 5 confirms the use of a ROM to determine the momentum sink term. Finally, conclusions are drawn.

**Figure 1.** Summary of the suggested approach. On the left side, the mixture model is modified by Feenstra's correlation to compute the void fraction. On the right side, the Bi-CITSGM method applied to a REV is introduced in the porous media approach to define the pressure drop of the tube bundle.

#### **2. Cross-Flow through a Horizontal Tube Bundle**

This study focuses on the simulation of the thermal-hydraulic behaviour of a twophase adiabatic crossflow through a horizontal tube bundle. The tube bundle is modelled by a porous medium that involves solving the conservation equations of mass and momentum and the turbulence model using the superficial velocity porous formulation. The flow is considered two-dimensional and only transverse to the tube bundle. All geometric notations pertaining to the tube bundle are illustrated in Figure 2.

**Figure 2.** Geometric definitions of the tube bundle.

#### *2.1. Governing Equations of the Two-Phase Flow*

In order to simulate the adiabatic two-phase flow, the mixture model is used. Only two phases are considered: the liquid phase *l* and the gas phase *g*. The mixture model allows the phases to be interpenetrating and to move at different velocities using the concept of slip velocities. The continuity equation for the mixture is

$$\frac{\partial}{\partial t}(\rho\_b) + \nabla \cdot (\rho\_b \vec{u}\_b) = 0 \tag{1}$$

where *ub* is the mass-averaged velocity defined by

$$
\vec{u}\_b = \frac{\alpha\_\mathcal{\mathcal{S}} \rho\_\mathcal{\mathcal{S}} \vec{u}\_\mathcal{\mathcal{S}} + \alpha\_I \rho\_I \vec{u}\_l}{\rho\_b}. \tag{2}
$$

*α<sup>k</sup>* is the volume fraction of phase *k*, and *ρ<sup>b</sup>* is the mixture density, defined by

$$
\rho\_b = \mathfrak{a}\_{\mathfrak{F}} \rho\_{\mathfrak{F}} + \mathfrak{a}\_l \rho\_l. \tag{3}
$$

The momentum equation for the mixture is obtained by summing the individual momentum equations of all phases. It takes the following form

$$\begin{split} \frac{\partial}{\partial t} (\rho\_b \vec{u}\_b) + \nabla \cdot (\rho\_b \vec{u}\_b \otimes \vec{u}\_b) &= -\nabla p + \nabla \cdot \left[ \mu\_b \left( \nabla \vec{u}\_b + \nabla \vec{u}\_b^T \right) \right] + \rho\_b \vec{g} \\ &- \nabla \cdot \left( \kappa\_\delta \rho\_\delta \vec{u}\_{d\tau,\emptyset} \vec{u}\_{d\tau,\emptyset} + \kappa\_l \rho\_l \vec{u}\_{d\tau,l} \vec{u}\_{d\tau,l} \right) + \vec{F}\_l. \end{split} \tag{4}$$

*μ<sup>b</sup>* is the mixture viscosity such as

$$
\mu\_b = \alpha\_\mathcal{S} \mu\_\mathcal{S} + \alpha\_l \mu\_l \tag{5}
$$

and *<sup>g</sup>* <sup>=</sup> <sup>−</sup>*g<sup>y</sup>* with *<sup>g</sup>* <sup>=</sup> 9.81 m/s<sup>2</sup> . The hydraulic resistance of tubes on the fluid is taken into account by the Darcy–Forchheimer term *Ft*. This is a source term due to the use of the porous media approach and it is defined in the following Section 2.3. *udr*,*<sup>k</sup>* is the drift velocity of the phase *k* defined by

$$
\vec{u}\_{dr,k} = \vec{u}\_k - \vec{u}\_b. \tag{6}
$$

Moreover, the volume fraction equation for phase *g* is built from the continuity equation for the gas phase, using the definition of the drift velocity (Equation (6)) to eliminate the phase velocity as:

$$\frac{\partial}{\partial t} \left( a\_{\mathcal{S}} \rho\_{\mathcal{S}} \right) + \nabla \cdot \left( a\_{\mathcal{S}} \rho\_{\mathcal{S}} \vec{u}\_b \right) = - \nabla \cdot \left( a\_{\mathcal{S}} \rho\_{\mathcal{S}} \vec{u}\_{dr,\mathfrak{E}} \right). \tag{7}$$

The drift velocity is related to the relative (or slip) velocity according to:

$$\vec{u}\_{dr\_{\vec{\mathcal{S}}}} = \vec{u}\_{\mathcal{J}} \left( 1 - \frac{\alpha\_{\mathcal{J}} \rho\_{\mathcal{S}}}{\rho\_{\mathcal{b}}} \right) \quad \text{and} \quad \vec{u}\_{dr\_{\vec{\mathcal{I}}}} = \vec{u}\_{\mathcal{J}} \left( \frac{\alpha\_{\mathcal{I}} \rho\_{\mathcal{I}}}{\rho\_{\mathcal{b}}} - 1 \right). \tag{8}$$

The most used algebraic slip formulation is the Manninen model [27]. With this formulation, the slip velocity is given by

$$\vec{u}\_{\mathcal{S}^l} = \frac{\tau\_{\mathcal{S}}}{f\_{drag}} \frac{\left(\rho\_{\mathcal{S}} - \rho\_b\right)}{\rho\_{\mathcal{S}}} \,\,\vec{a} \tag{9}$$

where *τ<sup>g</sup>* is the particle relaxation time,

$$
\pi\_{\mathcal{S}} = \frac{\rho\_{\mathcal{S}} d\_{\mathcal{S}}^2}{18\mu\_{\mathcal{S}}} \tag{10}
$$

*dg* is the gas particle diameter and*a* is the acceleration.

$$
\vec{a} = \vec{g} - (\vec{u}\_b \cdot \nabla)\vec{u}\_b - \frac{\partial \vec{u}\_b}{\partial t} \tag{11}
$$

The drag force *fdrag* is obtained with the model of Schiller and Naumann [11]. Commonly used, it is expressed as a function of the drag coefficient *CD* and the relative Reynolds number:

$$f\_{drag} = \frac{C\_D R \varepsilon}{24} \quad \text{and} \quad \text{Re} = \frac{\rho\_l ||\vec{u}\_{\mathcal{S}} - \vec{u}\_l||d\_{\mathcal{S}}}{\mu\_l} \tag{12}$$

This slip velocity formulation is not suitable for a two-phase flow in tube bundles because it is not designed for such a configuration and does not take into account the associated physical phenomena. In the next subsection, a formulation more adequate for the tube bundle configuration is suggested.

#### *2.2. Rewriting of the Slip Velocity*

The void fraction depends on the slip ratio *S* and quality *x* according to

$$\varepsilon = \frac{1}{1 + S \frac{\rho\_{\infty}}{\rho\_{\parallel}} \frac{1 - \chi}{\chi}} \tag{13}$$

Here, the slip velocity is reformulated as a function of slip in order to introduce a slip ratio model adapted to a two-phase through a tube bundle. The slip velocity is the velocity difference between the gas phase and the liquid phase.

$$
\vec{u}\_{\mathcal{S}^I} = \vec{u}\_{\mathcal{S}} - \vec{u}\_I \tag{14}
$$

The velocity components in Equation (14) can be written as

$$\begin{cases} u\_{\mathcal{S}l,x} = u\_{l,x} \left( \frac{u\_{\mathcal{S},x}}{u\_{l,x}} - 1 \right) = u\_{l,x} (\mathcal{S}\_x - 1) \\ u\_{\mathcal{S}l,y} = u\_{l,y} \left( \frac{u\_{\mathcal{S},y}}{u\_{l,y}} - 1 \right) = u\_{l,y} (\mathcal{S}\_y - 1) \end{cases} \tag{15}$$

*ul*,*<sup>x</sup>* (resp. *ul*,*y*) is the liquid velocity in the *x* (resp. *y*) direction. The same notation is used for the gas velocity and the slip velocity.

For an upward cross-flow to the tube bundle in the y direction, it can be assumed that the velocity ratio *Sx* is equal to 1 and that *Sy* is defined by a correlation coming from the literature. Finally, the drift velocity is written as

$$\begin{cases} u\_{dr,\emptyset,x} = u\_{l,x} (\mathcal{S}\_x - 1) \left( 1 - \frac{a\_{\mathcal{S}} \rho\_{\mathcal{S}}}{\rho\_b} \right) \\ u\_{dr,\emptyset,y} = u\_{l,y} (\mathcal{S}\_y - 1) \left( 1 - \frac{a\_{\mathcal{S}} \rho\_{\mathcal{S}}}{\rho\_b} \right) \end{cases} \tag{16}$$

The proposed approach needs to solve Equations (1), (4) and (7) with the modified definition of the drift velocity defined by Equation (16).

#### 2.2.1. Hibiki's Correlation (2017)

Smith [18] gives the slip ratio valid for all void fraction ranges and for vertical, inclined, and horizontal flows in a channel [28,29] as

$$S = \varepsilon + (1 - \varepsilon) \left( \frac{\frac{\rho\_l}{\rho\_\mathcal{g}} + \left( \frac{1 - \chi}{\chi} \right) e}{1 + \left( \frac{1 - \chi}{\chi} \right) e} \right)^{0.5} \tag{17}$$

where *e* is the ratio of the mass of liquid droplets entrained in the gas core to the total mass of liquid. In order to implement a void fraction correlation in a steam generator thermal-hydraulic code, Hibiki et al. defined this parameter for staggered tube bundles like Dowlati et al. as

$$
\varepsilon = \min(0.0637 N\_{\text{th,p}}^{0.871}, 1). \tag{18}
$$

*Nm*˙ ,*<sup>p</sup>* = *m*˙ *<sup>p</sup>*/*ρ<sup>g</sup> jg*,*crit* where the critical superficial gas velocity is

$$j\_{\rm g,crit} = \left(\frac{\Delta\rho\varrho\sigma}{\rho\_{\rm g}^2}\right)^{0.25} \left(\frac{\mu\_l}{\left(\rho\_l\sigma\sqrt{\frac{\sigma}{\rho\_{\rm g}\Delta\rho}}\right)^{0.5}}\right)^{-0.2}.\tag{19}$$

where *σ* is the surface tension. *m*˙ *<sup>p</sup>* is the pitch mass flux, which represents the mixture velocity between two tubes *ub*,*p* multiplied by the mixture density.

$$\dot{m}\_{\mathcal{P}} = a\_{\mathcal{S}} \rho\_{\mathcal{S}} \|\vec{u}\_{\mathcal{S},\mathcal{P}}\| + a\_{l} \rho\_{l} \|\vec{u}\_{l,p}\| \quad \text{where} \ \|\vec{u}\_{k,p}\| = \|\vec{u}\_{k}\| \frac{P}{P - D\_{ext}} \quad \text{for } k = \{l, \mathcal{g}\} \tag{20}$$

2.2.2. Feenstra's Correlation (2000)

Feenstra et al. [16] defined the slip ratio as

$$S = 1 + 25.7 \sqrt{Ri \ast Cap} \left(\frac{P}{D\_{ext}}\right)^{-1} \tag{21}$$

*Dext* is the outer diameter of the tubes, and *P* is the pitch, illustrated in Figure 2. The Richardson number is the ratio between the buoyancy force and the inertia force:

$$Ri = \frac{\left(\rho\_{\mathcal{S}} - \rho\_l\right)^2 \mathcal{g}(P - D\_{\rm crit})}{\left.m\_p\right|^2} \tag{22}$$

The Capillary number is the ratio between the viscous force and the surface tension force:

$$\text{Cap} = \frac{\mu\_l ||\vec{u}\_{\text{\textquotedblleft}p}||}{\sigma} \tag{23}$$

The gas phase velocity is based on the resulting void fraction:

$$|||\vec{u}\_{\mathcal{S},p}|| = \frac{\mathbf{x}\dot{m}\_p}{\varepsilon\rho\_{\mathcal{S}}}\tag{24}$$

#### *2.3. Definition of the Darcy–Forchheimer Term*

The source term *Ft* = *Ft*,*xx* + *Ft*,*yy* is added to the momentum equation because tube bundles are represented by a porous medium and is usually written as:

$$\begin{cases} \begin{array}{rcl} F\_{t, \mathbf{x}} &=& -\left( D\_{\mathbf{x} \mathbf{x}} \mu\_{b} \boldsymbol{u}\_{b, \mathbf{x}} + \mathcal{K}\_{\mathbf{x} \mathbf{x}} \frac{1}{2} \rho\_{b} ||\vec{\boldsymbol{u}}\_{b}|| \boldsymbol{u}\_{b, \mathbf{x}} \right) \\\ F\_{t, \mathbf{y}} &=& -\left( D\_{\mathbf{y} \mathbf{y}} \mu\_{b} \boldsymbol{u}\_{b, \mathbf{y}} + \mathcal{K}\_{\mathbf{y} \mathbf{y}} \frac{1}{2} \rho\_{b} ||\vec{\boldsymbol{u}}\_{b}|| \boldsymbol{u}\_{b, \mathbf{y}} \right) \end{array} \tag{25}$$

This term is composed of a viscous loss term and an inertial loss term resulting from Darcy–Forchheimer's law [19]. *ub* is the mixture velocity magnitude, *Dxx* and *Dyy* are the inverse of the permeability, and *Kxx* and *Kyy* are the correction terms of Forchheimer. In this study, the first term of Equation (25) is neglected because the flow is turbulent. As the term has the same dimensions as a pressure gradient, Equation (25) is rewritten by:

$$F\_{t, \mathbf{x}} = -\frac{\Delta P\_{f, \mathbf{x}}^{2\Phi}}{N\_{\mathbf{R}, \mathbf{x}} P\_{\mathbf{x}}} \quad \text{and} \quad F\_{t, \mathbf{y}} = -\frac{\Delta P\_{f, \mathbf{y}}^{2\Phi}}{N\_{\mathbf{R}, \mathbf{y}} P\_{\mathbf{y}}} \tag{26}$$

Δ*P*2<sup>Φ</sup> *<sup>f</sup>* ,*<sup>i</sup>* is the two-phase frictional pressure drop, *NR*,*<sup>i</sup>* is the number of tube rows, and *Pi* is the pitch in direction *i* (*x* or *y*) shown in Figure 2. From Equations (25) and 26, the unknown coefficients *Kxx* and *Kyy* are defined by:

$$\begin{cases} \begin{array}{rcl} K\_{xx} &=& -\frac{\Delta P\_{f,x}^{\Delta \Phi}}{\frac{1}{2}\rho\_b \|\|\vec{u}\_b\|\|u\_{b,x}N\_{R,x}} \frac{1}{P\_x} \\\ K\_{yy} &=& -\frac{\Delta P\_{f,y}^{\Delta \Phi}}{\frac{1}{2}\rho\_b \|\|\vec{u}\_b\|\|u\_{b,y}N\_{R,y}} P\_y \end{array} \end{cases} \tag{27}$$

From the method developed by Consolini et al. [30] to define the two-phase frictional pressure drop over horizontal tube bundles, Equation (27) is reduced to:

$$K\_{\rm xx} = -\frac{\lambda Eu}{P\_{\rm x}} \left( \frac{P\_{\rm y}}{P\_{\rm y} - D\_{\rm ext}} \right)^2 \quad \text{and} \quad K\_{\rm yy} = -\frac{\lambda Eu}{P\_{\rm y}} \left( \frac{P\_{\rm x}}{P\_{\rm x} - D\_{\rm ext}} \right)^2 \tag{28}$$

where the Euler number, *Eu*, can be given by Zukauskas et al. [21]. Zukauskas et al. defined a correlation for the Euler number resulting from their experiments and experimental results from the literature. This law enables to determine frictional pressure drops for in-line and staggered tube bundles with 1.25 <sup>≤</sup> *<sup>P</sup>*/*Dext* <sup>≤</sup> 2.5 and 10 <sup>≤</sup> *Re* <sup>≤</sup> <sup>10</sup>6. The Euler number is calculated as :

$$Eu = k\_1 \sum\_{i=0}^{4} \frac{c\_i}{Re^i} \tag{29}$$

where *ci* and *k*<sup>1</sup> are coefficients given in reference [21]. The two-phase multiplier coefficient is written by Consolini et al. [30] as:

$$
\lambda = \Lambda + (1 - \Lambda)(1 - 2x)^2 \quad \text{with} \quad \Lambda = \left(\frac{\dot{m}\_p}{400}\right)^{-1.5} \tag{30}
$$

It is important to note that the two-phase multiplier factor is equal to 1 when the quality tends towards 0 (only liquid phase) and 1 (only gas phase). This key argument is not available with the correlation of Ishihara et al. [31] based on the Lockhart–Martinelli approach [32], which was the method used by Dowlati et al. in [26].

For complex or non-standard geometries, such as helicoidal tube bundles or corrugated tubes [33–35], there are no suitable or reliable correlations in the literature to compute the Forchheimer force *Ft*. They can be obtained by experiments or numerical simulations by LES or RANS approaches but these methods imply a significant cost. In order to reduce computational time and cost, we suggest using a non-intrusive reduced model model (ROM). This parametric ROM simulates the flow in an REV (Representative Elementary Volume) of the tube bundle. The knowledge of this flow enables the Forchheimer term to be quickly obtained. In addition, the spatial distributions of the pressure and the velocity around each tube are precisely given by this approach. The ROM methodology is detailed in the next section.

#### **3. Reduced Ordel Model on the REV to Compute the Forchheimer Term**

#### *3.1. Representative Elementary Volume*

A Representative Elementary Volume (REV) of the tube bundle is considered in Figure 3, where the lower and upper boundary conditions are periodical. Governing equations for periodically fully developed flow are derived from the incompressible Navier–Stokes equations defined by:

$$\begin{cases} \nabla \cdot \vec{u} = 0\\ \frac{\partial \vec{u}}{\partial t} + \nabla \cdot (\vec{u} \otimes \vec{u}) = -\frac{\nabla p}{\rho} + \nu \nabla^2 \vec{u} + \frac{\vec{\beta}}{\rho} \end{cases} \quad \text{where} \quad \vec{\beta} = \begin{pmatrix} 0\\ \beta \end{pmatrix} \tag{31}$$

*u* is the periodic flow velocity:

$$\vec{u}(\mathbf{x}, y, t) = \vec{u}(\mathbf{x}, y + 2P\_{y\_\prime}t) \tag{32}$$

*P* represents the reduced pressure, which satisfies periodic boundary conditions,

$$P(\mathbf{x}, y, t) = P(\mathbf{x}, y + 2P\_{y'}t) \tag{33}$$

and the actual pressure is given by

$$p(\mathbf{x}, y, t) = -\beta y + P(\mathbf{x}, y, t) \tag{34}$$

according to Patankar [36]. *β* is the linear component of the pressure, which is to be calculated iteratively for a fixed mass flow rate [37,38].

**Figure 3.** REV of the tube bundle of Dowlati et al. [26].

#### *3.2. Definition of Reduced Bases*

Model reduction techniques make it possible to quickly and inexpensively obtain the temporal dynamics of a complex flow. The principle of the reduced-order model (ROM) is to approximate the solution in a small dimension sub-vector space, which enables capturing the dominant characteristics of the physical phenomenon studied. The solution *f*(*t*, *x*) is written as a linear combination of a finite number of spatial basis functions Φ*k*(*x*) as:

$$f(t, \mathbf{x}) \approx \sum\_{k=1}^{q} a^k(t) \Phi^k(\mathbf{x}) \tag{35}$$

*q* is relatively small compared to the problem size, and *a<sup>k</sup>* is the temporal coefficient. Here, *f* corresponds to the velocity *u* or the pressure *p*.

The most common method to compute the spatial basis Φ is the POD (Proper Orthogonal Decomposition) method [39,40]. However, the ability of the POD basis function to give the dynamic of the phenomenon studied is dependent of the information contained in the snapshots that form the basis. For instance, a POD basis built with snapshots for a Reynolds number *Re*<sup>1</sup> will not be able to predict the dynamics of the physical phenomenon for another Reynolds number *Re*2. To increase the validity domain of the POD basis, it is possible to interpolate a set of POD bases Φ1, ··· , Φ*<sup>N</sup>* built for different values of Reynolds numbers *Re*1, ··· , *ReN* in order to obtain the basis associated to the desired Reynolds number. Standard interpolation techniques (RBF, Lagrange, Spline, etc.) are not very efficient and are not generally representative of the phenomenon studied. To get around this difficulty, a basis interpolation approach based on the results of differential geometry and, more particularly, on the properties of the Grassmann manifold can be used [25,41–43]. In this work, we consider the approach offered by Amsallem et al. [25] that is subsequently noted as ITSGM (Interpolation on the Tangent Subspace of the Grassmann Manifold). The algorithm is given in Algorithm A1.

Once the spatial basis for the desired parameter is defined, the temporal coefficients are usually computed by solving a system of differential equations (ROM) resulting from the Galerkin projection of the full model on the basis functions Φ*k*(*x*) [44]. This method has the disadvantage of being costly and intrusive. Indeed, for each new value of Reynolds number, it is necessary to compute the coefficients of the ROM, which is costly. Moreover, derivative operators are difficult to assess using a commercial CFD code. Consequently, in this paper, we use the non-intrusive method Bi-CITSGM [23,45].

#### *3.3. Description of the Bi-CITSGM Method*

The methodology of the Bi-CITSGM is given in Algorithm A2 in the Appendix B. The first step of this method is to interpolate the spatial basis and the temporal basis, built by POD, using the ITSGM method. The singular value matrix is acquired by classical interpolation methods such as Lagrange, RBF, and Spline. Then, the ranking step aims to sort the interpolated spatial and time eigen modes with respect to the interpolated singular values. Indeed, as the spatial and temporal bases may not be in the same order as the singular values, orthogonal matrices need to be introduced. These calibration matrices are a solution to an optimization problem under constraints whose solution is analytically determined. For more details, see [23,24].

#### **4. Validation of the Using of the Modified Mixture Model on the Dowlati's Experiment** *4.1. Study Configuration*

The methodology, presented in Section 2, was validated with the experiment of Dowlati et al. They made void fraction and friction pressure drop measurements for vertical two-phase flow of air–water across staggered in-line tube bundles with different pitch-to-diameter ratios. Here, the tube bundle, illustrated in Figure 4, is made up of 20 tube rows in a staggered arrangement with five tubes in each row and the ratio *P*/*Dext* is 1.75. Geometry dimensions are given in Table 1. The estimated uncertainties in the data done by Dowlati et al. are detailed in Table 2. The present study is done with quality range between 1.3 <sup>×</sup> <sup>10</sup>−<sup>4</sup> and 3 <sup>×</sup> <sup>10</sup>−<sup>2</sup> and mass flux range between 164 kg/(m2· s) and 538 kg/(m2· s).

**Table 1.** Dimensions of the tube bundle.


**Table 2.** Estimated uncertainties in the data given by Dowlati et al. [26].


CFD calculations were performed with Ansys Fluent v2020R2. The tubes in the tube bundle were not represented. The tube bundle was modelled as a porous medium by adding a source term in the momentum equation. At the inlet, the homogeneous void fraction model is assumed. Thus, the homogeneous void fraction *ε<sup>H</sup>* and inlet phase velocities are written as follows,

$$\vec{u}\_{\vec{g},in} = \begin{pmatrix} 0\\ \frac{x\eta\_{in}}{\rho\_{\vec{g}}\varepsilon\_H} \end{pmatrix}\_{\left(\vec{x},\vec{g}\right)}; \quad \vec{u}\_{l,in} = \begin{pmatrix} 0\\ \frac{(1-x)\eta\_{in}}{\rho\_l(1-\varepsilon\_H)} \end{pmatrix}\_{\left(\vec{x},\vec{g}\right)}\tag{36}$$

$$\varepsilon\_{H} = \frac{1}{1 + \frac{\rho\_{\mathcal{X}}}{\rho\_{I}} \frac{1 - \mathbf{x}}{\mathbf{x}}} \tag{37}$$

where *x* is the quality. Calculations are initialized from the input boundary conditions. The mesh, illustrated in Figure 5, is defined in such a way as to ensure that the dimensionless wall distance *y*<sup>+</sup> is close to 1. The turbulent model k-*ω* SST [46] is used.

**Figure 4.** Staggered tube bundle from Dowlati's experiment (*P*/*Dext* = 1.75).

**Figure 5.** Mesh used for CFD calculations (size for the central cells = 3 mm/cell size along the border = 0.75 mm).

#### *4.2. Influence of the Slip Model on the Void Fraction Prediction*

In order to evaluate the influence of the slip model on the void fraction prediction, mean relative and absolute errors are introduced.

$$E\_{rel} = \frac{1}{N\_{pt}} \sum\_{i=1}^{N\_{pt}} \frac{\overline{\varepsilon}\_i^{calc} - \overline{\varepsilon}\_i^{exp}}{\overline{\varepsilon}\_i^{exp}} \tag{38}$$

$$E\_{\rm abs} = \frac{1}{N\_{pt}} \sum\_{i=1}^{N\_{pt}} \left( \overline{\varepsilon}\_i^{calc} - \overline{\varepsilon}\_i^{exp} \right) \tag{39}$$

*Npt*, *εcalc*, and *εexp* are, respectively, the number of data and the computed and experimental void fractions averaged over the tube bundle. To prove that we need to apply a void fraction model to the mixture model, calculations were performed with the slip velocity formulation of Manninen. Here, the value of the void fraction postprocessed with this formulation always matches with the homogeneous void fraction. Figures 6 and 7, for, respectively, *m*˙ *<sup>p</sup>* = 164 kg/(m2.s) and *m*˙ *<sup>p</sup>* = 401 kg/(m2.s), showed that the homogeneous void fraction, and thus the Manninen's formulation, do not take into account the slip between phases for flows in tube bundle. Indeed, the relative error is about *Erel* = 51.05% and the absolute error *Eabs* = 0.15. The results stress that it is important to implement a void fraction model appropriate for a two-phase cross-flow through a horizontal tube bundle. Moreover, for each mass flux, Figures 6 and 7 demonstrate the correct implementation of Feenstra or Hibiki's correlation in the CFD code. That is to say, the errors calculated subsequently in this paper derive only from the slip ratio model implemented and not from the numerical model. It can be seen in Figure 8, which depicts void fraction results as functions of quality and mass flux, that the void fraction increases with the quality and with the mass flux. Roser [47] justified this phenomenon by the upward movement of the gas phase against the liquid phase due to the buoyancy force making it all the more important that the mass flux is low. For higher mass flux, the gap with the homogeneous void fraction model is less significant. Indeed, the two phases are "well mixed" due to the increase in turbulence.

Figure 9 depicts CFD-computed void fractions versus experimental void fractions for all mass fluxes. Feenstra's correlation always underpredicts the void fraction; however, errors are acceptable with *Erel* = 17.49% and *Eabs* = 0.07. Errors seem to be more important for high mass fluxes and low qualities with some errors outside of range ±20%. Hibiki's correlation overpredicts the results for *ε* ≤ 0.5 and underpredicts it for *ε* > 0.5. The relative error is a little higher than Feenstra's correlation with *Erel* = 21.81%, but the absolute error

is better with *Eabs* = 0.05. Some void fraction points are outside the range ±20% for low void fractions; however, this correlation is better for higher void fractions.

**Figure 6.** Void fraction results for air-water flow through Dowlati's staggered tube bundle with *m*˙ *<sup>p</sup>* = 164 kg/(m2.s). (**a**) Feenstra's correlation; (**b**) Hibiki's correlation.

**Figure 7.** Void fraction results for air-water flow through Dowlati's tube bundle with *m*˙ *<sup>p</sup>* = 401 kg/(m2.s). (**a**) Feenstra's correlation; (**b**) Hibiki's correlation.

**Figure 8.** Void fraction as a function of quality and mass flux ( = results from CFD simulations and -= experimental measurements). (**a**) Feenstra's correlation; (**b**) Hibiki's correlation.

**Figure 9.** CFD computed void fraction versus experimental void fraction. (**a**) Feenstra's correlation; (**b**) Hibiki's correlation.

In order to improve the void fraction prediction, the Capillary number used in Feenstra's correlation is now expressed as a function of upstream mass flux. The upstream mass flux and the mass flux between two tubes are linked by:

$$
\dot{m}\_p = \dot{m}\_{upstream} \frac{P}{P - D\_{ext}} \tag{40}
$$

With this modified Feenstra's correlation, named "Upstream Feenstra's correlation", the results were always improved compared to the experimental results, the original Feenstra's correlation, and Hibiki's correlation. Indeed, Figure 10 shows that the results with the modified Feenstra's correlation were closer to experimental results than the original Feenstra's correlation. It is important to note that only 10% of the simulations had a relative error for the void fraction higher than 20% for the modified Feenstra's correlation. As can be seen in Figure 11, these points were located at low void fractions. For the higher void fraction, results were very close to the experiment. Compared to other correlations, this one has an absolute error of 0.03 and a relative error of 9.10%, which confirms the accuracy of the void fraction prediction.

Figure 12 plots the slip ratio obtained by the experiment and Feenstra, Hibiki and upstream Feenstra's correlations as a function of quality for each mass flux. As the void fraction increases with mass flux, it can be noticed that slip ratio decreases with the quality until reaching 1. For *<sup>m</sup>*˙ *<sup>p</sup>* <sup>≤</sup> 247 kg/(m2.s), the slip ratio turns out to be relatively constant regardless of the quality. However, no correlation studied here captures this phenomenon. For *<sup>m</sup>*˙ *<sup>p</sup>* <sup>≥</sup> 329 kg/(m2.s), the slip ratio increases with the quality, and correlations follow the same trend. Overall, the slip ratio obtained by the upstream Feenstra correlation is the closest to experimental results.

Now that we have verified the accurate prediction of the void fraction by rewriting the relative velocity and implementing a suitable void fraction model, we can focus on the pressure drop of the system. The porous media approach implies adding a source term to the momentum equation. We prove, in the next section, the possibility to use a non-intrusive parametric reduced model on an REV in order to compute this term.

**Figure 11.** CFD computed void fraction with the modified Feenstra's correlation versus experimental void fraction.

**Figure 12.** Comparison of slip ratio obtained by the different methods.

#### **5. Computation of the Pressure Drop in Tube Bundles by Using POD-ROM on an REV**

The porous media approach involves the implementation of the Euler number (Equation (28)) in the Forchheimer's correction term. Usually, this dimensionless number is computed with correlations coming from the literature. Here, we want to show that it is possible to compute this variable by a non-intrusive reduced-order model, Bi-CITSGM, applied to a REV of the tube bundle. The case of the tube bundle of Dowlati et al. is an case of application for which there is already a correlation resulting from the literature and for which the resolution of a reduced model is of little interest. However, we would like to extend this developed methodology to a case of a complex tube bundle for which there is no correlation from the literature. In addition, it is less expensive to define a correlation resulting from numerical calculations rather than from tests. First, the Bi-CITSGM method is validated on the REV, and then, the implementation on CFD simulations of the Dowlati's experiment is discussed.

#### *5.1. Validation of the Bi-CITSGM Method on the REV*

Figure 3 depicts the 2-D REV of the tube bundle of Dowlati et al. with symmetric and periodic boundary conditions. The dimensions of transverse pitches and the tube diameter are given in Table 1. CFD calculations were done with OpenFOAM. The governing equations were unsteady and dimensionless, and the k-*ω* SST turbulence model is used. The flow is computed until the fully developed flow is established. The density of the fluid was always equal to 1 kg/m3, the desired mass flow was 4.7625 <sup>×</sup> <sup>10</sup>−<sup>2</sup> kg/s, and the viscosity varied according to the desired Reynolds number.

All the considered training Reynolds numbers are {2000; 3000; ...; 29,000; 30,000}. This range of training Reynolds numbers fits with the Reynolds number in each cell of the Dowlati's simulation. For each Reynolds number, a thousand time steps are kept, and these snapshots are once and for all decomposed by POD. For each new Reynolds number, the Bi-CITSGM method gives the results almost instantly. In Figures 13–15, solving the ROM with five or ten POD modes are compared to the reference CFD calculation for each Reynolds number. The Bi-CITSGM method enables determining the pressure field with an accuracy less than 10%, except for Reynolds numbers less than 5000 (Figure 13). To increase the accuracy of pressure fields, the spacing between the training Reynolds numbers should be reduced in this range. However, the key parameter in this section is the pressure drop. Figures 14 and 15 show that the pressure drop prediction is much less than 10% even for Reynolds numbers less than 5000. For instance, at a Reynolds number of 24,500, the absolute error between the CFD reference pressure field and the interpolated pressure field at two given times is represented in Figure 16. The highest errors are rather local in the REV, and overall, the interpolated pressure field is very close to the CFD reference. Following the good results achieved by the Bi-CITSGM method on the REV, the results in the next subsection are plotted by keeping only five POD modes.

**Figure 13.** Relative error of the time- and area-weighted pressure average between Bi-CITSGM method and CFD calculation.

**Figure 14.** Time-weighted pressure drop average for Bi-CITSGM method and CFD calculations.

**Figure 15.** Relative error of the time weighted-pressure drop average between Bi-CITSGM method and CFD calculation.
