In our work in 2008, we evaluated the aptitude of the code Neptune_CFD to reproduce the incidence of a structure topped by vanes on a boiling layer, within the framework of the Neptune project. The objective was to reproduce the main effects of the spacer grids. The turbulence of the liquid phase was modeled by a first-order - model. We show in this paper that this model is unable to describe the turbulence of rotating flows, in accordance with the theory. The objective of this paper is to improve the turbulence modeling of the liquid phase by a second turbulence model based on a - approach. Results obtained on typical single-phase cases highlight the improvement of the prediction for all computed values. We tested the turbulence model - implemented in the code versus typical adiabatic two-phase flow experiments. We check that the simulations with the Reynolds stress transport model (RSTM) give satisfactory results in a simple geometry as compared to a - model: this point is crucial before calculating rod bundle geometries where the - model may fail.

1. Introduction

In a pressurized water nuclear reactor (PWR), an optimum heat removal from the surface of the nuclear fuel elements (rod bundle with spacer grids) is very important for thermal margin and safety.

One important goal is to carry out sensitivity analyses on the angle of the vanes of the fuel assembly spacer grids. In [1], the critical heat flux (CHF) experiment on the effect of the angle and of the position of mixing vanes was performed in a rod bundle. The authors show that the mixing vanes increase the value of the CHF and the result is correlated to the magnitude of the swirl generated by the mixing vanes. If the angle of the mixing vanes is relatively small, the magnitude of the swirling flow is smaller because the rotating force created by the mixing vanes is weak. If the angle of the mixing vanes is relatively large, the mixing vanes play the role of flow obstacle under the departure from nucleate boiling (DNB) condition. Therefore, it is important that the turbulence modeling deals with rotation effects.

There have been several studies on flow mixing and heat transfer enhancement caused by a mixing-vane spacer grid in rod bundle geometry. Lee and Choi [2] simulate the flow field and heat transfer in a single-phase flow for a rod bundle with eight spans of mixing vanes. The FLUENT commercial code is employed and a Reynolds stress transport model (RSTM) is used for turbulence. According to the authors, RSTM is helpful. Ikeda et al. [3] study an assembly consisting of a heater rod bundle and eight specific mixing vane grids. For Ikeda et al., it might be insufficient to apply a standard model to swirl-mixing flow and narrow-channel flow conditions that include nonisotropic effects. Moreover, In et al. [4] have performed a series of CFD single-phase flow simulations to analyze the heat transfer enhancement in a fully heated rod bundle with mixing-vane spacers. For future work, In et al. recommend that a refined computational fluid dynamic (CFD) model be developed to include details of the grid structure and a higher-order turbulence model be employed to improve the accuracy of such simulations.

Additional two-phase effects like accumulation of bubbles in the center of subchannel or pockets of bubbles on the rods should be taken into account to improve the simulation of flow close to DNB. Indeed, single-phase simulations remain insufficient and boiling flows simulations are required. In [5], the authors describe CFD approaches to subcooled boiling and investigate their capability to contribute to fuel assembly design. A large part of their work is dedicated to the modeling of boiling flows and to forces acting on the bubbles. The authors note that the size of bubbles in the bulk is correlated to the local subcooling which is an important parameter (see [1]).

Considering flow in a PWR core in conditions close to nominal, when boiling occurs, a high-velocity steady flow takes place with very small times scales associated to the passage of bubbles (10-4 second −10-3 second) and with quite small bubble diameters (10-5 m to 10-3 m) compared to the hydraulic diameter (about 10-2 m). According to the synthesis of the work performed in WP2.2 of the NURESIM project [6], these are perfect conditions to use a time average or ensemble average of equations as usually done in the RANS approach. All turbulent fluctuations and two-phase intermittency scales can be filtered since they are significantly smaller than the scales of the mean flow.

The large-eddy simulation is also a possible approach. In the context of the NURESIM project, several studies have been carried out with a large-eddy simulation to study the axial development of air-water bubbly flows in a pipe. But in the synthesis of the work performed in WP2.2, Bestion [6] notes that several open modeling and numerical issues still remain. So, we will focus on the RANS approach in this paper.

According to all these recommendations, a better understanding of the detailed structure of a flow mixing and heat transfer downstream of a mixing-vane spacer in a nuclear fuel rod bundle has to be investigated with an RSTM.

Moreover, the overwhelming majority of industrial CFD applications today are still conducted with two-equation eddy viscosity model, especially the standard model, while RSTM remains exceptional.

As an example of RSTM development, an RSTM model adapted to bubbly flows is studied in [7] and used to perform simulations of three basic bubbly flows (grid, uniform shear, and bubbly wake). The authors decomposed the Reynolds stress tensor of the liquid into two independent parts: a turbulent part produced by the mean velocity gradient that also contains the turbulence of the bubble wakes and a pseudoturbulent part induced by bubble displacements; each part is predicted from a transport equation. This model is interesting but has not been selected here for the following reasons. Firstly, the computation effort is doubled (turbulent and pseudoturbulent parts). Secondly, considering the flow close to nominal PWR core conditions, when boiling occurs, a high-velocity steady flow takes place and the bubble diameter is quite small (10-5 m to 10-3 m), therefore the bubbles follow the liquid streamlines and so the modeling of the pseudoturbulent part induced by bubble displacements can be omitted. Thirdly, the two-phase flow modeling proposed in [7] does not tend to a single-phase flow formulation when the void fraction tends to zero. These three arguments have imposed the choice of the higher-order turbulence model described in the paper.

A second-order moment turbulence model for simulating a bubble column is also proposed in [8]. The authors defined a Reynolds tensor for each phase. Furthermore, following a similar method used in deriving and closing the Reynolds stress equations, a modeled transport equation of two-phase velocity correlation is also solved. For the same reasons as those mentioned above for the model proposed in [7], we have not adopted the model proposed in [8]. The turbulence modeling described in the present paper takes into account the Reynolds tensor for the liquid only, while a more basic modeling is used for the vapor phase.

However, to the authors’ knowledge, no industrial CFD approach for boiling flows with an RSTM approach is available in the context of the fuel assembly design. This may be due to the fact that numerical problems may occur when using an RSTM approach without caution. Furthermore, the turbulence modeling of boiling flows is not straightforward. The use of RSTM also requires finer meshes than eddy viscosity models (EVMs) and RSTM may therefore be more time and storage consuming. Developing an industrial RSTM approach is a quite challenging task but it is worth working at it: RSTM and EVM results will probably differ and it is important to determine what consequences it may have on thermal margin and safety of reactors.

In the framework of an R&D program carried out in the Neptune project (EDF, CEA, AREVA-NP, IRSN), the following strategy has been adopted:

(1)validation of the Neptune_CFD code with an RSTM approach on single-phase flow with mixing vanes and on more academic cases of air-water adiabatic bubbly flows in a pipe;(2)validation of the Neptune_CFD code with an RSTM approach on boiling flows in a pipe and sensitivity to the angle of the vanes for fuel assembly spacer grids performed in a rod bundle in a boiling flow configuration;(3)validation of the Neptune_CFD code with an RSTM approach for a rod bundle with mixing vanes currently used for commercial nuclear fuel. Step (1) is described in this paper. The second step should be finalized by the end of 2008. The step (3) is not yet started. Our main objective in this paper is to check that the simulation with the RSTM gives satisfactory results in a simple geometry as compared to an EVM: this point is crucial before calculating rod bundle geometries where the EVM model may fail.

This paper is organized as follows. In Section 2.2, the general model we use for adiabatic bubbly flow simulations is presented in details. In Section 2.3, we underline the weaknesses of the EVM models. In Sections 3.1 and 3.2, the second-moment closure model for high Reynolds number two-phase flows is presented. In Section 2.4, we give some examples of typical PWR problems that EVM models fail to represent. In Sections 4.1 and 4.2, respectively, the Liu and Bankoff case and the sudden expansion experiment are briefly described. The comparison of the results of Neptune_CFD calculations and the experimental data are presented. The sensitivity of the numerical results to the turbulence model for the fluid and to the most important models is studied.

Finally, conclusions are drawn about our current capabilities to simulate bubbly flows with an RSTM model and perspectives for future work are given.

2. The Neptune_CFD Solver and Physical Modeling

2.1. Introduction

Neptune_CFD is a three-dimensional two-fluid code developed more especially for nuclear reactor applications. This local three-dimensional module is based on the classical two-fluid one pressure approach, including mass, momentum, and energy balances for each phase.

The Neptune_CFD solver, based on a pressure correction approach, is able to simulate multicomponent multiphase flows by solving a set of three balance equations for each field (fluid component and/or phase) [912]. These fields can represent many kinds of multiphase flows: distinct physical components (e.g., gas, liquid, and solid particles), thermodynamic phases of the same component (e.g., liquid water and its vapor), distinct physical components, some of which split into different groups (e.g., water and several groups of different diameter bubbles), and different forms of the same physical components (e.g., a continuous liquid field, a dispersed liquid field, a continuous vapor field, and a dispersed vapor field). The solver is based on a finite volume discretization, together with a collocated arrangement for all variables. The data structure is totally face-based, which allows the use of arbitrary shaped cells (tetraedra, hexahedra, prisms, pyramids, etc.) including nonconforming meshes (meshes with hanging nodes).

2.2. Governing Equations and Physical Modeling

The CFD module of the Neptune software platform is based on the two-fluid approach [13, 14]. In this approach, a set of local balance equations for mass, momentum, and energy is written for each phase. These balance equations are obtained by ensemble averaging of the local instantaneous balance equations written for the two phases. When the averaging operation is performed, the major part of the information about the interfacial configuration and the microphysics governing the different types of exchanges are lost. As a consequence, a certain number of closure relations (also called constitutive relations) must be supplied for the total number of equations (the balance equations and the closure relations) to be equal to the number of unknown fields. We can distinguish three different types of closure relations: those which express the interphase exchanges (interfacial transfer terms), those which express the intraphase exchanges (molecular and turbulent transfer terms), and those which express the interactions between each phase and the walls (wall transfer terms). The balance equations of the two-fluid model we use for adiabatic bubbly flows and their closure relations are described in the following subsections.

2.2.1. Main Set of Balance Equations

The two-fluid model we use for our adiabatic bubbly flow calculations consists of the following balance equations.

Two mass balance equations where t is the time, denote the time fraction of phase k, its averaged density and velocity. The phase index k takes the values l for the liquid phase and g for gas bubbles.

Two momentum balance equations where p is the pressure, is the gravity acceleration, is the interfacial momentum transfer per unit volume and unit time, and and denote the molecular and turbulent stress tensors, the latter being also called the Reynolds stress tensor.

The interfacial transfer of momentum appearing in the RHS of (2) is assumed to be the sum of four forces:

The four terms are the averaged drag, added mass, lift, and turbulent dispersion forces per unit volume. Now we will give the expressions we use for these forces and for their coefficients.

(i)Drag force where CD is the drag coefficients for bubbles which can be determined experimentally.(ii)Added mass force where is the added mass coefficient which is equal to for a spherical bubble and the factor takes into account the effect of the bubbles concentration [15, 16].(iii)Lift force where CL is the lift coefficient. This coefficient is equal to 1/2 in the particular case of a weakly rotational flow around a spherical bubble in the limit of infinite Reynolds number [17].(iv)Turbulent dispersion force where Kl is the liquid turbulent kinetic energy and is a numerical constant of order 1. This expression was proposed by Lance and Lopez de Bertodano [18]. An alternative approach is proposed by [19, 20] to model the turbulence induced by bubbles: an algebraic model developed in the framework of Tchen’s theory, where the turbulent kinetic energy for the dispersed phase and the covariance are calculated from the turbulent kinetic energy of the continuous phase.

For the dispersed phase, the Reynolds stress tensor is closed using a Boussinesq-like hypothesis: where is the identity tensor, with , the turbulent viscosity for the dispersed phase, , the gas turbulent kinetic energy, , the covariance of the dispersed phase, , the ratio between the time scale of the continuous phase turbulence viewed by the dispersed phase (takes into account crossing trajectories effect) and the characteristic time scale of the momentum transfer rate between the liquid and dispersed phases: with is the turbulent Schmidt or Prandtl turbulent for the continuous phase, is the crossing trajectories coefficient taken equal to 1.8, and is the added mass coefficient; Turbulent dispersion force and Tchen’s model will be compared below for bubbly flows in a straight pipe and in a sudden expansion.

2.2.2. Turbulent Transfer Terms

The model describes energy processes in terms of production and dissipation, as well as transport through the mean flow or by turbulent diffusion. The Kolmogorov spectral equilibrium hypothesis also enables one to predict a large eddy length-scale. On the other hand, the anisotropy of the stresses is quite crudely modeled. First of all the EVM model assumes the Reynolds stress tensor is aligned with the strain rate tensor (Boussinesq approximation): where is the identity tensor, Kl is the liquid turbulent kinetic energy, and is the liquid turbulent eddy viscosity. The liquid turbulent eddy viscosity is expressed by the following relation: where . The turbulent kinetic energy Kl and its dissipation rate are calculated by using the two-equation approach.

2.3. EVM Weaknesses: Theoretical Approach

Flows encountered in vertical pipe are of great interest to validate the most important heat, mass, and momentum closure relations. However, some negligible effects in simple geometry sometimes become preponderant in complex geometries. For example, the modeling of two-phase flow in water-cooled nuclear reactors needs to take into account swirls and stagnation points. Applications in complex geometries also need to take into account the complex features of the secondary motions which are observed experimentally. These requirements highlight the need for meticulous turbulence modeling.

A reason for the persistent widespread use of low-level turbulence modeling in two-phase CFD is perhaps the fact that the use of two-phase CFD in complex industrial geometries is only starting. Moreover, many studies merely require “order of magnitude” or “good tendencies” answers.

However, extensive testing and application over the past three decades have revealed a number of shortcomings and deficiencies in EVM models, and among them the model, such as

(i)limitation to linear algebraic stress-strain relationship (poor performances wherever the stress transport is important, e.g., non equilibrium, fast evolving, separating, and buoyant flows),(ii)insensitivity to the orientation of turbulence structure and stress anisotropy (poor performances where normal stresses play an important role, e.g., stress-driven secondary flows in noncircular ducts),(iii)inability to account for extra strain (streamline curvature, skewing, rotation),(iv)poor prediction particularly of flows with strong adverse pressure gradients and in reattachment regions. In a plane strain situation, such as upstream of a stagnation point on a bluff body, the exact (as obtained by an RTSM) production and that obtained from an EVM are, respectively, [21] The difference between the normal stresses actually grows slowly on the short-time scale needed for the flow to travel around the stagnation point, so the production remains moderate, and in any case is bounded, whereas usually yields a severe over-prediction when the strain is high.

The simulation of swirling flow generated by the mixing vanes is our main goal since it plays an important role for the prediction of the CHF for the fuel assemblies. For this reason, the rotation effects are more specifically addressed hereafter.

It can be easily shown [22] that in the presence of an initially anisotropy turbulence, rotation will cause a redistribution of energy between normal components without affecting the value of this quantity. In fact, the angular velocity does not appear explicitly in the K-equation, obtained by adding the normal stresses: Thus, the model is totally blind to rotation effects. The swirling flows can be regarded as a special case of fluid rotation with the axis usually aligned with the mean flow direction so that the Coriolis force is zero. This aspect is crucial for the simulation of hot channel of a fuel assembly. In fact, mixing vanes at the spacer grids generate a swirl in the coolant water to enhance the heat transfer from the rods to the coolants in the hot channels and to limit boiling.

In the following section, we present some examples of large-scale industrial applications, performed using eddy-viscosity models, and subsequently discuss areas of weakness of the models, highlighting some improvements that can be obtained through the use of more advances stress transport closures.

2.4. EVM Weaknesses: Illustration on the AGATE-Mixing and DEBORA-Mixing Experiment

Keeping in mind the long-term objective (two-phase CFD calculations validated under typical pressurized water reactor (PWR) geometries and thermalhydraulic conditions), we started very recently to evaluate Neptune_CFD against spacer grid type experiments. An experimental device representing three mixing blades (Figure 2) was introduced in a heated tube (diameter = 19.2 mm) and used for two different programs.

(i) AGATE-mixing experiment [23]. Single-phase liquid water tests, with laser-doppler liquid velocity measurements upstream and downstream the mixing blades (for each of the 15 horizontal planes, the liquid velocity is measured along 12 different diameters and there are 12 points for each radius); the velocity at inlet is 3 m/s and the pressure is 2 bar.

(ii) DEBORA-mixing experiment [24]. Boiling R-12 freon tests on the same geometry but the total length of the calculation domain is 3.5 m; the tube is heated and the uniform wall heat flux is 109300 W/m2 which gives about 2% of vapor at outlet; the outlet pressure is 26.2 bar; the inlet liquid temperature is 63.3°C; the inlet liquid mass is 0.873 kg/s. The main physical phenomena to reproduce are wall boiling, entrainment of bubbles in the wakes, and recondensation. So, the prediction of the swirl is crucial.

For the mixing blades part (60 mm), 77000 cells are needed. This grid is considered as a reasonable compromise between the numerical accuracy and the computational effort. Figure 1 compares computed and experimental orthoradial (circular component in a horizontal plane) liquid velocity downstream the mixing vane (AGATE test). One can notice that the rotating flow is qualitatively well reproduced by Neptune_CFD although the velocity is underestimated. This is mainly due to the turbulence model (standard here) which is not optimum for this type of geometry.

The model underestimates orthoradial velocities downstream the blades, but the results remain qualitatively satisfactory. The model gives satisfactory results (Figure 1).

In the following section, we propose a second-moment closure model to take into account the liquid turbulence in order to validate in the long-term calculations in typical pressurized water reactor (PWR) geometries and thermalhydraulic conditions.

In the present paper, we suppose that RSTM is well known in single-phase flow [21]. Now, our objective is to test and if possible to improve our RSTM model adapted to bubbly flows as compared to experimental data and results. Indeed, we are interested by two-phase high Reynolds numbers flows, but beforehand, the mechanical models implemented in the Neptune_CFD code must be tested on the simpler cases of air-water adiabatic bubbly flows.

3. The Second-Moment Closure Model for High Reynolds Number Flows Dedicated to the Continuous Phase (Liquid)

In this section, we omit the subscript “l” for the liquid and “” is the void fraction for the sake of simplicity.

3.1. Equation on Rij

We have In this model, the Reynolds stress tensor of the continuous phase is split into two parts, a turbulent dissipative part produced by the gradient of mean velocity and by the wakes of the bubbles and a pseudoturbulent nondissipative part induced by the displacements of the bubbles. The displacements of the bubbles should be taken into account in experiments, where air is injected at the bottom of a water pool creating a large, axisymmetric bubble plume with a large-scale recirculation flow around the plume. But swirling flows and high Reynolds number characterize our industrial applications.

Hence, we neglected, in our approach and in first analysis, the nondissipative component called “pseudoturbulent.” We consider only the “turbulent” dissipative part. Within this framework, the term of production by the bubbles interfaces is written as [7] where n indicates the normal to the interface and a Dirac function on the interface. It was omitted in [7]. Indeed, according to [7], dissipation in the wakes is balanced by the interfacial production: the equation of transport of the Reynolds stress tensor has the same form as in the single-phase case and is given by (15). When the void fraction is vanishing, the two-phase flow modeling naturally degenerates to the single-phase flow modeling.

Some terms of the equation of transport of the Reynolds stress tensor cannot be computed directly and must be modeled. A modeling resulting from [21] is proposed below.

A common way to model the viscous destruction of stresses for high Reynolds number flows is The turbulent diffusion is of diffusive nature and the most popular model is the generalized gradient diffusion: Pressure fluctuations tend to disrupt the turbulent structures and to redistribute the energy to make turbulence more isotropic: with




(iv) where xn is the distance to the wall and the base vector normal to the wall.

(v) G is the production by body force.

3.2. Equation on

In the RTSM closures, the same basic form of model equation for ε is used as in the model, except that now is available, which has the following implications.

The production of kinetic energy (P and G) in the source term of is treated in exact form.

The generalized gradient hypothesis is used to model turbulent diffusion.

Hence, the model equation for has the form The coefficients of the model are shown in Table 2.

4. Validation on Adiabatic Bubbly Flow Cases

4.1. The LIU&BANKOFF Case [25]

In this section, we evaluate the modeling capabilities to simulate an upward bubbly flow in adiabatic conditions. We do not specifically optimize the coefficients and modeling of the momentum transfer terms to get results as good as possible. Our main objective is to validate the turbulence model for the fluid against the one.

The test section was a  mm long, vertical smooth acrylic tubing, with inside diameter of  mm.

The set of physical properties is the following: A uniform axial liquid profile is imposed at the inlet and is equal to 1.138 m/s. A uniform axial gas profile is imposed at the inlet and is equal to 1.333 m/s. The void fraction at the inlet is 0.045.

The interfacial momentum transfer term is assumed to be the sum of four different forces. The three first ones are simplified averaged expressions of the classical drag, added mass, and lift force. The fourth is the turbulent dispersion force.

The flow is assumed to be axisymmetric therefore a two-dimensional axisymmetric meshing is used. Computations have been performed on two kinds of meshing: a coarse grid (20 cells in the radial direction and 50 cells in the axial direction) and a fine grid (30 cells in the radial direction and 100 cells in the axial direction). Results are similar and computations are performed on the first grid.

At the measuring station (), we compare numerical results against experimental data for the axial liquid velocity and the void fraction.

As recommended in [26], the lift coefficient is taken to be equal to 0.1. We have also tested the turbulent dispersion force of Davidson model [26] written as This expression gives values negligible with respect to the Lopez de Bertodano expression [18] (the ratio is about 1000) and the void fraction profile is observed to be similar to the profile calculated without any turbulence model for the dispersed phase. If we take the bubble fluctuation into account with the Hinze-Tchen algebraic model of bubble turbulence, we get the same results with the turbulent dispersion force.

In [27], Grossetête considers that bubbles are deformed near the wall. To take into account this effect, the author proposes to put a negative lift coefficient near the wall and otherwise a positive one. Nevertheless, calculations (not presented in this paper) show a peak of void fraction at the wall. We have also tested the Tomiyama lift force [28, 29], but results are not improved. Moreover, a wall lubrification force [30] can push the bubbles away from the wall and improve the results.

Figures 3 and 4 show that the simulation results are in quite good agreement with the experimental data. In [31], computations performed with a turbulence model for the liquid produced comparable results. Our main objective in this paper is to check that the simulation with the turbulence model gives satisfactory results in a simple geometry, which is crucial before calculating industrial geometries, where the turbulence model is susceptible to fail.

But improvement of the modeling of the interfacial forces exerted on bubbles by the surrounding liquid is required. A strong sensitivity to the lift coefficient has been found in our calculations. Other forces, like the turbulent dispersion force, have also a crucial effect on the void fraction distribution. These forces depend on uncontrolled parameters like the bubble shape, the liquid turbulence, the bubbles collective effects, and so on [31].

4.2. Sudden Expansion Experiment

Bel Fdhila [32] investigated experimentally several upwards bubbly flows in a vertical pipe with a sudden expansion. The total length of the pipe was equal to 14 meters. The bottom part of the tube had a length equal to 9 meters and an internal diameter equal to 50 mm, the top part of the tube (5 m length) having an internal diameter equal to 100 mm. The fluids used were water and air under atmospheric pressure and ambient temperature. Six measuring sections were located upstream and downstream of the singularity. The first measuring section was located two centimetres before the singularity, the other five were located at 7, 13, 18, 25, and 32 cm above the singularity. In each measuring section, the radial profile of the void fraction was measured by means of a single optical probe, and two components of the liquid velocity were measured by means of a hot film anemometer. The time-averaged components of this velocity field and three components of the liquid Reynolds stress tensor were deduced (the flow being assumed axisymmetric). The bubble size was not measured in the experiment. According to the author, the observed bubble diameter was equal to a few millimetres, the bubbles remaining relatively small due to the strong turbulence existing in the liquid phase.

In our calculations, only a small part of the tube, containing the singularity and the six measuring sections, was reproduced. The radial profiles of the void fraction and the liquid mean and fluctuating velocities measured upstream of the singularity were used as inlet conditions. The length of the calculation domain is equal to 38 cm. The flow is assumed to be axisymmetric therefore a two-dimensional axisymmetric meshing is used. Several calculations have been done in order to test the sensitivity to the axial and radial grids. Four different grids have been tested in [33], the number of radial meshes in the largest section multiplied by the number of axial meshes being 10*38, 20*76, 20*152, and 40*152, respectively. The comparison of different calculations of the same two-phase flow, realized on these four different grids show that the calculations performed with the finest grid can be considered as converged. All the calculations presented here have been done on the finest 40*152 calculation grid. The flow studied here is characterized by the liquid and gas superficial velocities and the area-averaged void fraction in the two sections given in Table 1. It can be noted that the averaged void fraction has important values for this test (12%).

Cokljat [34] performed calculations of the sudden expansion experiment with the FLUENT code. Predictions were obtained using the standard K-ε model as well as an RSTM for the continuous phase, while the turbulence closure for the dispersed phase is achieved by the algebraic model of Tchen [19, 20]. With this approach, similar as ours, the authors indicate that both models produce similar results for the axial velocity but void fraction results are improved with the RSTM model.

We only consider the classical drag, added mass, and dispersion turbulent force. The dispersion force coefficient is equal to 2 in the computations. The bubble diameter is equal to 2 mm. Following [18, 33, 34], we underline the necessity to discard the lift force. In fact, the effect of the lift force is to produce sharp peaks near the wall because the classical modeling of the lift force seems not well adapted to this case.

The simulation results are in reasonable agreement with the experimental data for the void fraction profiles (Figure 5) and have been improved as compared to [33, 34]. But the profiles at  cm and  cm show that the void fraction is underestimated which mean that a better understanding of the physical mechanisms is still needed.

Especially for the axial and radial mean liquid velocity profiles(Figure 6,7), we have obtained a good agreement which means that the recirculation zone is well captured.

We have obtained only qualitatively good results for the RMS quantities (Figure 8) because the turbulence mechanisms in a bubbly flow are far from being fully understood [18]. But the turbulence modeling of the dispersed phase in a PWR core in conditions close to nominal is less crucial than the liquid turbulence modeling.

Finally, results with the turbulence model for the fluid are similar to the one, which is our main objective in this case, before calculating rod bundle geometries, where the model may fail.

5. Conclusion and Perspectives

An analysis of turbulence modeling for two-phase flows has been proposed. Indeed, the use of eddy viscosity models is widespread and may be sufficient for parallel flows in vertical pipes, but that type of model does not account for effects that are preponderant in complex geometries, especially when swirling flows are involved, for example, in pressurized water reactor cores downstream of mixing vanes and spacer grids. In accordance with the theory, it is demonstrated in the case of a flow downstream of a mixing vane that using a Reynolds stress model is an efficient way to improve the simulation of such complex flows. To demonstrate that the use of a Reynolds stress model is not bound to deteriorate the classical results obtained with an eddy viscosity model, a validation step on more analytical experiments is detailed (bubbly flows in a straight pipe and in a sudden expansion): the study shows that the Reynolds Stress model implemented in the multiphase 3D code Neptune_CFD satisfactorily reproduces the results obtained with the standard eddy viscosity model and both compare reasonably well with the experiments.

As concern the computational cost, we note that in the case of the DEBORA-mixing test which is under process, the time required by iteration is, respectively, 3.09 seconds and 2.81 seconds for RSTM and EVM. The time step is, respectively, 5 milliseconds and 5.4 milliseconds with a CFL equal to 1. In this particular case, the RSTM over-cost is about 18.8%.

Moreover, among the developments planned in the medium term, we have identified the need for a polydispersion model.

Besides, the Neptune project has set up a medium and long-term experimental program to acquire detailed measurements in simplified and real geometries, both in adiabatic and real conditions [9, 10].

Ai: Interfacial area concentration
Cd: Drag coefficient
dt: Numerical time step
g: Gravity acceleration
Kl: Liquid turbulent kinetic energy
: Interfacial momentum transfer per unit volume and unit time
p: Pressure
Prl: Liquid Prandtl number
: Reynolds stress tensor
Reb: Bubble Reynolds number
t: Time
: Fluctuation of the liquid velocity
: Averaged velocity of phase k
: Interfacial-averaged velocity
: Denotes the time fraction of phase k
: Dissipation rate
: Gas molecular viscosity
νl: Liquid kinematic viscosity
: Liquid turbulent eddy viscosity
ρk: Averaged density of phase k
: Surface tension
τw: Wall shear stress
: Molecular stress tensor.

l: Liquid state
g: Gas bubbles
k: Phase k = l or g.


This work has been achieved in the framework of the Neptune project, financially supported by CEA (Commissariat à l’Energie Atomique), EDF (Electricité de France), IRSN (Institut de Radioprotection et de Sûreté Nucléaire), and AREVA-NP.