NUMERICAL PREDICTION OF TWO-PHASE FLOW IN TUBE BUNDLES WITH A CFD POROUS MEDIA APPROACH BASED ON MIXTURE MODEL AND A VOID FRACTION CORRELATION

Being able to predict the void fraction is essential for a numerical prediction of the thermohydraulic behaviour in steam generators. Indeed, it determines two-phase mixture density and affects two-phase mixture velocity which enable to evaluate the pressure drop of heat exchanger, the mass transfer and heat transfer coefficients. In this study, the flow is modelled by coupling Ansys Fluent with an in-house code library where a CFD porous media approach is implemented. In this code, the two-phase flow has been modelled so far using the Eulerian model. However, this two-phase model requires interaction laws between phases which are not known and/or reliable for a flow within a tube bundle. The aim of this paper is to use the mixture model, for which it is easier to implement suitable correlations for tube bundles. By expressing the relative velocity, as a function of slip, the void fraction model of Feenstra et al. developed for upward cross-flow through horizontal tube bundles is introduced. With this method, physical phenomena that occur in tube bundles are taken into consideration in the mixture model. The developed approach is validated based on the experimental results obtained by Dowlati et al. NOMENCLATURE ?⃗? dr Drift velocity m/s ?⃗? gl Relative velocity m/s ?⃗? b,p Mixture velocity between two tubes m/s ?̇?p Mass flux between two tubes kg/(m2.s) F t Source term N/m 3 ΔPf 2Φ Two-phase frictional pressure drop Pa S Slip (Velocity ratio) τg Particle relaxation time s fdrag Drag function CD Drag coefficient ε Void fraction α Volume fraction εH Homogeneous void fraction λ Two-phase multiplier coefficient NR Number of tube rows Pi Pitch in the direction i m Dext Outer diameter of a tube m dg Gas particle diameter m INTRODUCTION Steam generators are heat exchangers used in particular 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 of 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 the heat exchanger and its correct dimensioning. However, modelling and simulating these heat exchangers by taking into account the tube bundle in detail, where there may be thousands of tubes, would require unacceptable computational cost and time. To decrease the computing time, the tube bundle can be modelled using the porous media theory. E3S Web of Conferences 321, 01002 (2021) ICCHMT 2021 https://doi.org/10.1051/e3sconf/202132101002 © The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/). Being able to predict precisely the void fraction is important to characterize the two-phase flow in heat exchangers. It is the key parameter for calculating the mixture density, the mixture viscosity and the mixture velocity. Therefore, it plays an important role for computing pressure drops and heat and mass transfers. Due to the relevance of the void fraction, many prediction methods have been proposed [1-2]. Since the 1980s, thermal-hydraulic codes have been developed to understand the physical phenomena involved. One of the first codes, THIRST [3], was developed to compute three-dimensional, two-phase and steady flow in steam generators. The two-phase flow was solved by the homogeneous two-phase model, the phase velocities were assumed to be equal but this does not even reflect what is really going on. To take into account the slip between phases, NavierStokes equations were solved for the mixture of the secondary fluid. The THYC code [4] gives the relative velocity thanks to a correlation using the drift flux model of Zuber and Findlay [5]. This model is based on the determination of drift-flux parameters for which many different empirical correlations exist [6-10]. Moreover, a porous medium approach was also used in this development and the primary fluid was solved with the one-dimensional heat equation. The porous media concept and the two-fluid model were used by Stosic and Stevanovic [11] 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. Stosic and Stevanovic [11] have modified the original correlation of Ishii and Zuber [12] by multiplying it by 0.4 and they validated it with their experimental data. Nevertheless, most of drag coefficient laws in the literature [13-15], including the Ishii and Zuber’s correlation, are made for different two-phase regime flows (bubbly, slug, stratified, annular or spray flow) in a tube and not in tube bundles. This problem can be tackled by using the mixture model, which is a simplified two-phase model where Navier-Stokes equations are solved for the mixture. In this study, the developed method involved formulating the relative velocity as a function of slip and consequently it was possible to implement a specific void fraction model for tube bundles derived from the literature. Two types of flows can be found in a U-tube steam generator: cross flow and parallel flow. This work was first validated for the cross-flow configuration with the experimental results of Dowlati et al. [16]. Feenstra et al. [17] proposed a slip ratio model developed for two-phase cross-flow in horizontal tube bundles. Thus, we rewrote the slip velocity used in the mixture model as a function of the slip in order to implement this void fraction model. In addition, the bundle was also modelled with a porous media approach. ANALYSIS AND MODELLING In order to simulate the thermal-hydraulic behaviour in the tube bundle, CFD calculations were performed with a porous medium approach. This involves solving the conservation equations of mass and momentum and the turbulence model using the superficial velocity porous formulation. Moreover, the two-phase flow was modelled using the mixture model of Ansys Fluent v2019R2. Governing equations of the two-phase flow: In this section, two phases are considered, the liquid phase l and the gas phase g. The flow is twodimensional and only transverse to the tube bundle. The continuity equation for the mixture is: ∂ ∂t (ρb) + ∇ ⋅ (ρb?⃗? b) = 0 (1) where ?⃗? b is the mass-averaged velocity: ?⃗? b = αgρg?⃗? g + αlρl?⃗? l ρb (2) αk is the volume fraction of phase k and ρb is the mixture density: ρb = αgρg + αlρl (3) The momentum equation for the mixture is obtained by summing the individual momentum equations of all phases. It takes the following form: ∂ ∂t (ρb?⃗? b) + ∇ ⋅ (ρb?⃗? b?⃗? b) = −∇p + ∇ ⋅ [μb(∇?⃗? b + ∇?⃗? b )] + ρbg (4) −∇ ⋅ (αgρg?⃗? dr,g?⃗? dr,g + αlρl?⃗? dr,l?⃗? dr,l) + F t F t is the source term due to the use of porous medium approach. μb is the mixture viscosity: μb = αgμg + αlμl (5) ?⃗? dr,k is the drift velocity of the phase k: ?⃗? dr,k = ?⃗? k − ?⃗? b (6) E3S Web of Conferences 321, 01002 (2021) ICCHMT 2021 https://doi.org/10.1051/e3sconf/202132101002


INTRODUCTION
Steam generators are heat exchangers used in particular 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 of 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 the heat exchanger and its correct dimensioning. However, modelling and simulating these heat exchangers by taking into account the tube bundle in detail, where there may be thousands of tubes, would require unacceptable computational cost and time. To decrease the computing time, the tube bundle can be modelled using the porous media theory.
Being able to predict precisely the void fraction is important to characterize the two-phase flow in heat exchangers. It is the key parameter for calculating the mixture density, the mixture viscosity and the mixture velocity. Therefore, it plays an important role for computing pressure drops and heat and mass transfers. Due to the relevance of the void fraction, many prediction methods have been proposed [1][2].
Since the 1980s, thermal-hydraulic codes have been developed to understand the physical phenomena involved. One of the first codes, THIRST [3], was developed to compute three-dimensional, two-phase and steady flow in steam generators. The two-phase flow was solved by the homogeneous two-phase model, the phase velocities were assumed to be equal but this does not even 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 [4] gives the relative velocity thanks to a correlation using the drift flux model of Zuber and Findlay [5]. This model is based on the determination of drift-flux parameters for which many different empirical correlations exist [6][7][8][9][10]. Moreover, a porous medium approach was also used in this development and the primary fluid was solved with the one-dimensional heat equation. The porous media concept and the two-fluid model were used by Stosic and Stevanovic [11] 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. Stosic and Stevanovic [11] have modified the original correlation of Ishii and Zuber [12] by multiplying it by 0.4 and they validated it with their experimental data. Nevertheless, most of drag coefficient laws in the literature [13][14][15], including the Ishii and Zuber's correlation, are made for different two-phase regime flows (bubbly, slug, stratified, annular or spray flow) in a tube and not in tube bundles. This problem can be tackled by using the mixture model, which is a simplified two-phase model where Navier-Stokes equations are solved for the mixture. In this study, the developed method involved formulating the relative velocity as a function of slip and consequently it was possible to implement a specific void fraction model for tube bundles derived from the literature. Two types of flows can be found in a U-tube steam generator: cross flow and parallel flow. This work was first validated for the cross-flow configuration with the experimental results of Dowlati et al. [16]. Feenstra et al. [17] proposed a slip ratio model developed for two-phase cross-flow in horizontal tube bundles. Thus, we rewrote the slip velocity used in the mixture model as a function of the slip in order to implement this void fraction model. In addition, the bundle was also modelled with a porous media approach.

ANALYSIS AND MODELLING
In order to simulate the thermal-hydraulic behaviour in the tube bundle, CFD calculations were performed with a porous medium approach. This involves solving the conservation equations of mass and momentum and the turbulence model using the superficial velocity porous formulation. Moreover, the two-phase flow was modelled using the mixture model of Ansys Fluent v2019R2.
Governing equations of the two-phase flow: In this section, two phases are considered, the liquid phase and the gas phase . The flow is twodimensional and only transverse to the tube bundle.
The continuity equation for the mixture is: where ⃗ is the mass-averaged velocity: is the volume fraction of phase and is the mixture density: The momentum equation for the mixture is obtained by summing the individual momentum equations of all phases. It takes the following form: is the source term due to the use of porous medium approach.
is the mixture viscosity: ⃗ , is the drift velocity of the phase : The drift velocity is related to the relative (or slip) velocity according to: In the CFD code used, the algebraic slip formulation from Manninen [18] was implemented. With this formulation, the slip velocity is given by: where is the particle relaxation time, is the gas particle diameter and is the acceleration.
For the model of Schiller and Naumann [13], commonly used, is expressed as a function of the drag coefficient and the relative Reynolds number: This slip velocity formulation is not suitable for a twophase flow in tube bundles because it was not designed for such a configuration and it does not take into account the associated physical phenomena.
Here, we proposed to reformulate the slip velocity in order to introduce a model adapted to a two-phase flow within a tube bundle. The slip velocity is the velocity difference between the gas phase and the liquid phase.
The velocity components in the previous equation can be written as: where (resp. ) is the slip according to direction (resp. direction). For an upward cross-flow to the tube bundle in the direction (resp. direction), it can be assumed that , = 0 (resp. , = 0) and therefore = (resp. = ).
For cross-flow, the velocity ratio of Feenstra et al. [17] defined by: is introduced in equation (15). The Richardson number is the ratio between the buoyancy force and the inertia force.
̇ is the pitch mass flux which represents the mixture velocity between two tubes ‖⃗ , ‖ multiplied by the mixture density: The Capillary number is the ratio between the viscous force and the surface tension force.
The gas phase velocity is based on the resulting void fraction: ‖⃗ , ‖ =̇ From the continuity equation for the gas phase, using the definition of the drift velocity (6) to eliminate the phase velocity, the volume fraction equation for the phase is solution of: The source term is added to the momentum equation because tube bundles are represented by a porous medium. It is expressed in / 3 and is usually written as: This term is composed of a viscous loss term and an inertial loss term resulting from Darcy-Forchheimer's law [19]. ‖⃗ ‖ is the mixture velocity magnitude and , , and diagonal coefficients of the matrices and . In this study, the first term of equation (23) is neglected because the flow is turbulent.

Moreover,
and are obtained from empirical correlations of two-phase frictional pressure drop in tube bundles coming from the literature. Equation (23) is rewritten as: 2Φ is the two-phase frictional pressure drop, , is the number of tube rows and is the pitch in direction ( ou ). From equations (23) and (24), unknown coefficients and are defined by: The two-phase friction factor 2Φ and the two-phase frictional pressure drop Δ 2Φ are linked through the relation: Consolini et al. defined a two-phase multiplier coefficient as the ratio between the two-phase friction factor and the single-phase friction factor [20]. 2Φ =

(27)
Their correlation is written as: 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. [21] based on the Lockhart-Martinelli approach [22] which was the method used by Dowlati et al. in [16].
To compute the two-phase frictional pressure drop, the single-phase friction factor must be defined.

= 4 (29)
Zukauskas et al. [23] proposed 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 inline and staggered tube bundles with 1.25 ≤ / ≤ 2.5 and 10 ≤ ≤ 10 6 . The Euler number is calculated as : where and 1 are coefficients given in reference [9].   Table 1.

RESULTS AND DISCUSSION
At the inlet, the homogeneous void fraction model is assumed, the slip is thus 1 and the phase velocities are the same. The homogeneous void fraction is defined from equation (32) and phase velocities are written as follows:  Calculations are initialized from the input boundary condition. The mesh is defined in such a way as to ensure that + is close to 1 (Figure 2). The turbulent model k-SST was used. The black line represents the homogeneous model defined by Equation (32) with = 1. This model always overpredicts the void fraction. For an upward two-phase cross-flow through a horizontal tube bundle, the gas phase and the liquid phase do not have the same velocity and this can be noticed in Figure 3. It can also be seen that the void fraction increases with the quality and also with the mass flux. Roser [24] justified this phenomenon by the upward movement of the gas phase against the liquid phase due to the buoyancy force is all the more important that the mass velocity is low. For higher mass velocity, the gap with the homogeneous void fraction model is less significant. Indeed, the two phases are "well mixted" due to the increase in turbulence.
The developed method (plotted in Figure 3 with Δ) is in agreement with the experimental results (plotted in Figure 3 with ▲). Feenstra's correlation is implemented in the CFD code and it can be seen in Figure 4 that it always underpredicts the void fraction compared to the experimental results. However, errors are acceptable and Figure 4 shows that the relative error is more important for high mass flux and low quality.   Void fraction results for ̇= 247 kg/(m².s) with the original Feenstra's correlation Feenstra used the mass flux between two tubes ̇ in order to compute the Capillary number in Equation (20) and Equation (21). In Figure 7 and Figure 8, we showed that it was better to express the Capillary number from the upstream mass velocity ̇. The upstream mass velocity and the mass velocity between two tubes are linked by: With this modified Feenstra's correlation, the results were always improved compared to the experimental results and the original Feenstra's correlation. Indeed, Figure 7 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 8, these points were located at low void fractions. For higher void fraction, results were closer to the experiment.   Computed total pressure drop with the modified Feenstra's correlation Figure 9 represents the total pressure drop obtained by CFD calculations with the modified Feenstra's correlation. Unfortunately, these values cannot be compared with the experimental ones because the experiments did not explicitly give the value of the frictional pressure drop and the post processing method causes too many uncertainty of measurement. In Figure 10, the contributions of two-phase frictional and gravitional pressure drops into the total pressure drop are plotted. The gravitational pressure drop decreases when the void fraction increases and is barely dependent on the mass flux. On the contrary, the two-phase frictional pressure drop is highly dependent on the mass flux and less on the void fraction. The frictional pressure drop depends on the void fraction when it is larger than 0.6.