Stephane Mimouni
email: stephane.mimouni@edf.fr

Solene Fleau

Stéphane Vincent

Stephane Vincent
email: stephane.vincent@u-pem.fr

CFD CALCULATIONS OF FLOW PATTERN MAPS AND LES OF MULTIPHASE FLOWS

whether they are published or not. The documents may come

INTRODUCTION

Two-phase flows are featuring many industrial applications such as nuclear power plants, heat exchangers and chemical reactors. Based on the interface structure, topological classification commonly separates two-phase flow regimes into three main groups. The first group would be separated flows such as annular, film or stratified flows, the second the dispersed flows such as bubbly, droplet or particle flows, and the last group would be the flows such as bubbly annular, churn turbulent or slug flow.

Although these regimes have been experimentally confirmed since a few decades [START_REF] Vincent | Flow pattern transitions in vertical and upwardly inclined lines[END_REF], their numerical simulation is still challenging and a universal model remains to be established. One of the roles of computational modeling is to better understand these kinds of complex regimes. However, most of today's models are dedicated to a single flow regime which is generally either a dispersed flow or a separated-phase flow.

Several kinds of approaches for two-phase flow modeling can be found in the literature [START_REF] Bestion | Applicability of two-phase CFD to nuclear reactor thermalhydraulics and elaboration of best practice guidelines[END_REF]. Among these models, the Eulerian two-fluid model, which formalism has been extended to more than two fields, is yet capable of dealing with both modeling of dispersed bubbly flows [START_REF] Mimouni | Combined evaluation of second order turbulence model and polydispersion model for two-phase boiling flow and application to fuel assembly analysis[END_REF] and simulation of larger bubble interfaces (Denefle et al., 2015). The dispersed approach for bubbly flow simulation is perfectly adapted to the Eulerian two-fluid formalism. The resolution of each balance equation needs closure terms (Ishii, 1975). For bubbly two-phase flows, the momentum closure laws involve several correlation stem from empirical observation. Therefore, these correlations rely on assumption on the bubble shape, generally considered as spherical or slightly ellipsoidal. To treat accurately problems where these assumptions are not verified, an interface tracking model is preferred to follow the distortions of the bubbles.

Interface tracking methods

Several methods have been investigated over the last decades for the interface tracking. A first one is to treat the interface in a Lagrangian point of view [START_REF] Hyman | Numerical methods for tracking interfaces[END_REF], by adapting the grid points to the interface motion. The second category is front tracking methods [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF]: the velocity field is solved on a fixed grid, and the position of the interface is tracked by markers transported by the carrying fluid. The markers should remain evenly distributed and therefore a redistribution step may be necessary, especially when dealing with break up or coalescence. The third and last main category consists in capturing the interface through a volume function. Depending on the way this function is transported and the interface properties such as the local curvature or the interfacial vector are calculated, two main methods can be named, the Volume of Fluid (VOF) [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] and the level-set method [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF]. The VOF method provides a conservative transport of the interface, but its properties are difficult to calculate accurately. On the other side, the level-set method automatically deals with topological changes of the interfaces but is not conservative. [START_REF] Olsson | A conservative level set method for two phase flow[END_REF] introduced a level set method in which the advection of the level-set function is followed by an artificial compression step to ensure that the thickness of the interface layer is preserved, inducing a volume conservation.

Multi-field approach with interface tracking in NURESAFE

Rather intense development of the so-called two-phase CFD has started about a decade ago. In Europe a small EU project EUROFASTNET, which started in 2000, identified industrial needs for 3D two-phase nuclear thermal hydraulics. It was followed by three large European projects of the EU EURATOM research programmes: NURESIM (2005)(2006)(2007)(2008)(2009), NURISP (2009[START_REF] Bestion | Applicability of two-phase CFD to nuclear reactor thermalhydraulics and elaboration of best practice guidelines[END_REF]), and NURESAFE (2013-2015). Topics of these projects were actually multi-physics simulations; however, 3D thermal-hydraulics represented a crucial part of these efforts. Results of the thermal hydraulic research are available in a special issue of Multiphase Science and Technology [START_REF] Bestion | Two-phase CFD: The various approaches and their applicability to each flow regime[END_REF], which was dedicated to state of the art in two-phase flow CFD. Some of the deliverables of NURESIM and NURISP are available on the home page of the NURESAFE project, while all NURESAFE deliverables will be in open-access due to the new policy of the EU Commission.

A special work-package within the NURESAFE project was focused on coupling of the interface tracking models and multi-field models. The idea of interface tracking method coupled with the two-fluid model was suggested by Černe et. al. in 2001, who proposed the coupling technique and performed analyses of the VOF interface tracking and 2D two-fluid model. This specific coupling, which was based on single-phase equations solved on one side of the interface and two-fluid equations on the other side, was found to be too complex for implementation in 3D codes. However, his basic idea was later developed in several different ways. [START_REF] Štrubelj | Numerical simulations of stratified two-phase flows with two-fluid model and interface sharpening[END_REF][START_REF] Štrubelj | Two-fluid model with interface sharpening[END_REF]Štrubelj et. al.( ) and Denefle et al. (2012) ) proposed implementation of the conservative level set interface tracking within the framework of two-fluid model. They treated the vapor volume fraction equation as a level set function discussed. A similar type of approach was implemented in the Neptune_CFD code by [START_REF] Coste | A large interface model for two-phase CFD[END_REF] and is known as "Large interface tracking". The model is supposed to recognize sufficiently large interfaces, while the interfaces smaller than the filter-size are modelled with the two-fluid model. While the (semi)empirical models that predict disintegration of the large interface into the small interfaces can be relatively simple, the reverse modelling, that is, formation of a large interface from the dispersed flow (for example stratification of the dispersed flow), presents a bigger challenge. In the Algebraic Interfacial Area Density model [START_REF] Porombka | Drag and turbulence modelling for free surface flows within the twofluid Euler-Euler framework[END_REF] no explicit interface reconstructionis done, but the location of the interface is identidfied. Basing on the two-fluid approach an interface sharpening is activated and a blending function allows the application of different closures for the bubbly flow, interface and droplet flow regions.

Development of the interface tracking within the 3D two-fluid or multi-field models, which would enable two-phase simulations for all flow regimes, was an effort of the NURESAFE project, which remains to be completed in the near future. Success of this approach would be an equivalent to the introduction of the "flow regime maps" in 1D two-fluid model. Nevertheless, the 3D equivalent will hopefully have much simpler structure: interfaces are tracked where the model can do it, whereas dispersed flow closure relations are used elsewhere. This approach will eventually end up with a much smaller number of empirical correlations compared to the general 1D system-code models. Recently several concepts to combine disperse flow models with the ones for large interfaces were proposed by different groups. The objective of this paper is to give an overview of the recent advances in the Large Bubble Model developped in the NEPTUNE_CFD code.

Large Bubble Model

We propose here to model separately the small and spherical bubbles and the larger and distorted ones (see Figure 1). This three-field simulation consists in treating the small bubbles as a dispersed field where the closure laws are well known terms, and locating the interface of the larger bubbles which deformation cannot be taken correctly into account by non-dimensional numbers and empirical correlations. Although the achievement of the proposed approach is quite new, the concept of such three-field model has already been described by some papers. [START_REF] Bestion | Applicability of two-phase CFD to nuclear reactor thermalhydraulics and elaboration of best practice guidelines[END_REF] exposed the idea of a so called LES Hybrid approach that consists in a spatial scale that separates the located interfaces scales from the smaller ones that are filtered and treated with a subgrid model. This paper is organized as follows. First, the equations solved and the modelling is detailled. Validation cases for the three-field model are presented in section 5, 6 and 7. The important topic of LES of twophase flows is adressed in section 8. Finally, conclusions are drawn about our current capabilities to simulate two-phase flows and perspectives for future work are given.

Figure 1: Illustration of the multifield approach used to simulate a two-phase flow, the liquid phase is in blue, the gaseous phase in white, the three fields defined in the multifield approach correspond in blue to the continuous liquid field, in red to the continuous gas field and in green to the dispersed gas field.

LARGE INTERFACES TREATMENT

Introduction

The simulation of the dispersed fields with the two-fluid model has been widely studied with the NEPTUNE_CFD code by [START_REF] Mimouni | A second order turbulence model based on a Reynolds stress approach for two-phase boiling flow. Part 1: Adiabatic cases[END_REF][START_REF] Mimouni | A second order turbulence model based on a reynolds stress approach for two-phase boiling flow and application to fuel assembly analysis[END_REF][START_REF] Mimouni | Combined evaluation of second order turbulence model and polydispersion model for two-phase boiling flow and application to fuel assembly analysis[END_REF]. Thus, the modelling effort concerns large interfaces. A first requirement is the implementation of a surface tension model and a drag force law to couple the velocity of each continuous field at the interface since, with the two-fluid model, each field has its own velocity. Then, because the two-fluid equations are solved for each field, the large interfaces are spread. Thus, the final step necessary to simulate accurately such interfaces with this model is to implement an interface sharpening equation to limit the interface smearing. All these three requirements are referred as the "Large Bubble Model" (LBM) and presented in this section after detailing the twofluid model equations. LBM is only applied to the continuous fields.

Two-fluid model

The flow motion is followed using the two-fluid model of [START_REF] Ishii | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF] extended to n-phases. In this model, density, viscosity, volume fraction and local velocity are defined for each field in each cell. In the present work, we only consider isothermal and incompressible flows. Thus, density and viscosity are constant for a given field in the whole domain. Two conservation equations are solved for each continuous field k:

 The mass balance equation:

𝜕 𝑡 α k 𝜌 𝑘 + ∇. α k 𝜌 𝑘 𝒖 𝒌 = 𝛤 𝑘
(1) With α k the volume fraction of field k, 𝜌 𝑘 its density, 𝑢 𝑘 its velocity and 𝛤 𝑘 interfacial mass transfers.

 The momentum equation:

𝜕 𝑡 (α k 𝜌 𝑘 𝒖 𝒌 ) + ∇. α k 𝜌 𝑘 𝑢 𝑖,𝑘 𝑢 𝑗 ,𝑘 = ∇. α k 𝑆 𝑖𝑗 ,𝑘 -α k ∇𝑃 + α k 𝜌 𝑘 𝑔 𝑖 + 𝐼 𝑖,𝑘 + 𝐹 𝑖,𝑘
(2) With 𝜇 𝑘 the viscosity of field k, 𝑆 the viscous stress tensor, 𝑃 the pressure, 𝑔 the gravitational constant, 𝐼 the interfacial momentum transfers and 𝐹 extra source terms (coupling between the continuous liquid and the continuous gas fields: surface tension and drag force law described in section 1.3). 𝐼 represents the coupling between the continuous liquid phase and the dispersed bubbles, i.e. the forces exerted on dispersed bubbles : drag force, lift force, wall lubrification, virtual mass and turbulent dispersion force. The closure laws used here and more details of the dispersed modeling of bubbly flows can be found in [START_REF] Mimouni | Combined evaluation of bubble dynamics, polydispersion model and turbulence modeling for adiabatic two-phase flow[END_REF].

A common pressure is assumed for all fields. A Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) solver is used in the code NEPTUNE_CFD. An iterative coupling of these equations is applied to ensure mass conservation (and energy conservation for non-isothermal flows). The data structure is face-based to allow simulations on arbitrary-shaped cells including non-conforming meshes. NEPTUNE_CFD inherits the I/O and HPC capabilities of the EDF open-source CFD software Code_Saturne [START_REF] Archambeau | Code_Saturne: a finite volume code for the computation of turbulent incompressible flows-Industrial Applications[END_REF].

Large Bubble Model

The Large Bubble Model is composed by a surface tension model, a drag force model and an interface sharpening equation. It has been developed for the simulation of large interfaces within a two-fluid model and the multifield approach. Thus, the model integrates the possible existence of at least one dispersed field in the flow.

 Surface tension model: The Continuum Surface Force (CSF) model proposed by [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] is used to model surface tension. Since the CFD tool is based on the two-fluid model, a volumetric expression of this force is implemented:

𝑭 𝑪𝑺𝑭 = α k 𝜎 κ ∇α k Ω (3)
With 𝜎 the surface tension coefficient, Ω the cell volume and κ the local curvature:

κ = -∇. ∇α k ∇α k (4)
 Drag force model: The role of the drag force model is to couple the velocity of the two continuous fields at the interface. This coupling term is expressed as a drag force law to be consistent with the dispersed fields, for which a drag force is applied. Thus, the expression of the drag force law presented in Eqn (5) corresponds to a continuous approximation of the drag force expression applied to the dispersed fields. || . Far from them (regions containing a continuous liquid (resp. Gas) field and a dispersed gas (resp. Liquid) field), 𝑑 𝑝 is equal to the diameter of the dispersed bubbles/droplets.

 Interface sharpening equation:

The last component of the Large Bubble Model is the interface sharpening equation (Strubelj, 2009). This element is a requirement for the simulation of large interfaces with a two-fluid model since it is known to spread interfaces. The choice has been made to adapt the equation proposed by [START_REF] Olsson | A conservative level set method for two phase flow[END_REF]:

𝜕 𝜏 α k + ∇. α k 1 -α k 𝒏 = 𝜖∆α k (6)
With n the interface normal vector:

𝒏 = ∇α k ||∇α k || (7)
And ∆𝜏 and 𝜖 two parameters, which fixes the final thickness of the large interfaces:

∆𝜏 = ∆𝑥 𝑚𝑖𝑛 32 and 𝜖 = ∆𝑥 𝑚𝑖𝑛 2 (8)
With ∆𝑥 𝑚𝑖𝑛 = min 𝛺 1 3 . With this set of parameters, whatever the initial diffusion of the interface, the iterative resolution of the interface sharpening equation ensures that the final interface thickness will be always equal to 5 cells (Denefle, 2015).

The implementation of the interface sharpening equation in a two-fluid model is not straightforward. It is necessary to pay attention first to mass conservation and second to the limitation of spurious velocities. In fact, since the interface sharpening equation modifies the distribution of volume fractions around large interfaces, this motion has to be done without losing mass. Indeed, the mass fluxes corresponding to the term ∇. α k 𝜌 𝑘 𝒖 𝒌 and the compression mass fluxes corresponding to ∇. α k 1 -α k 𝒏 -𝜖∆α k must be solved jointly in the calculation of the final volume fraction 𝛼 𝑘 𝑛+1 at the next time step n+1.

Moreover, this spurious displacement of the interface must stay negligible when compared to the interface motion governing by physical phenomena. Thus, two criteria have been introduced: one to activate the resolution of the equation when large interfaces (larger than 5 cells) are detected and the other one to stop the iterative resolution when the final thickness of 5 cells is reached. The iterative resolution is also limited in zones in the computational domain where the cells are larger. All these criteria and details about the implementation of the interface sharpening equation are presented in [START_REF] Fleau | Conservative implementation of the interface sharpening equation within an incompressible isothermal multifield approach[END_REF].

THREE FIELD MODEL

As explained previously, the three field model consists in a multifield approach where the gas phase is split into two separated fields. To build this approach, mass transfers have to be implemented between these two gas fields. Among these transfers, we can identify the interaction between two present gas fields, and the creation terms. The interaction can be reduced in a first approximation as the coalescence of the dispersed field with the located one. As a first demonstrating case, this coalescence is here achieved during the recompression step, and consists in transferring to the located field the dispersed field present inside the located bubble, i.e. where  cg >0.5. The continuous gas field creation terms consist in first dealing with the initiation of a located bubble where the dispersed field presents characteristics such as the volume fraction of the dispersed field is over a critical value (

3 . 0   crit dg  
) and when (Hansch, 2012):

x crit dg      10 1   (9)
The transfer of located field into the dispersed field during breaking up is the second part of the creation terms. The small and resolved pieces of located interface, typically when >1/(8.x) where x is the spatial grid step, should be transfer in the dispersed field. A bubble with a diameter inferior to 8x will present indeed, because of the CSF scheme for the curvature calculation, an overlap of the information used to calculate the curvature at two opposite points of the interface. More precisely, the dispersed gas field creation is activated where [START_REF] Mimouni | Multifield approach and interface locating method for twophase flows in nuclear power plant[END_REF]:

20 x d dg      (10)
The latter term proved to be largely predominant in calculations and compensates discrepancies in the calculation of the continuous gas field creation. As already discussed previously, only small bubbles of spherical shape should remain in the dispersed gas field. This condition is fulfilled if the Eotvos number is inferior to 4. It allows to determine the maximum value of the bubble diameter in the dispersed phase and to construct a grid such that this quantity is inferior to the average cell size.

VALIDATION TEST CASES FOR TWO FIELDS

The modelling of dispersed bubbly flows has been extensively validated in previous works [START_REF] Mimouni | A second order turbulence model based on a Reynolds stress approach for two-phase boiling flow. Part 1: Adiabatic cases[END_REF][START_REF] Mimouni | A second order turbulence model based on a reynolds stress approach for two-phase boiling flow and application to fuel assembly analysis[END_REF][START_REF] Mimouni | Combined evaluation of second order turbulence model and polydispersion model for two-phase boiling flow and application to fuel assembly analysis[END_REF](Mimouni et al., , 2016)).

Moreover, the Large Bubble Model have been previously validated on a large range of isothermal test cases (Fleau, 2015[START_REF] Fleau | Conservative implementation of the interface sharpening equation within an incompressible isothermal multifield approach[END_REF]: single static, oscillating or rising bubbles and also free surface test cases with Kelvin-Helmholtz instabilities and Rayleigh-Taylor instabilities. These cases deal with the simulation of large interfaces between a continuous gas field and a continuous liquid field. The main models involved are the surface tension model, the drag force model and the interface sharpening equation. The simulation results have been compared to experimental and theoretical data and other simulations with other CFD tools. Reasonable agreements have been obtained for all the cases and the mesh convergence has been carefully examined. In the following, a dispersed gas field is added in order to assess the capabilities of a multifield approach: continuous liquid field and continuous and dispersed gas fields.

CASTILLEJOS EXPERIMENT WITH THE MULTIFIELD APPROACH

After presenting the Large Bubble Model and validating it on different large interfaces, we propose here to study a three field test case: the Castillejos' test case (Castillejos, 1986). In this test case, air is injected at the bottom of a cylindrical bubble column (diameter and height equal to 0.5 m) at a constant mass flow rate equal to 876 cm 3 .s -1 . The injector diameter is equal to 6.35 mm and the water height to 0.4 m. Before reaching the free surface, the injected bubbles break up and form smaller inclusions treated as a dispersed field, such that the flow contains bubbles with sizes from about 6 cm to smaller than 1 mm. Void fractions profiles obtained experimentally are available at different heights. In the simulation, these profiles are obtained by averaging the void fraction (sum of the volume fractions of the continuous gas and the dispersed gas fields) for 20 s after 5 s of physical time.

The air injected at the bottom of the column is defined as a continuous gas field. A non-uniform mesh with 260000 cells is used to describe better the bubble interfaces close to the inlet (see Figure 2). The time step is constant, equal to 5. 10 -5 𝑠. The simulation is performed with a Reynolds Stress Transport Model (RSTM) to model the turbulence effects in a two-phase flow [START_REF] Speziale | Modelling the pressure-strain correlation of turbulence: an invariant dynamical systems approach[END_REF])- [START_REF] Mimouni | A second order turbulence model based on a Reynolds stress approach for two-phase boiling flow. Part 1: Adiabatic cases[END_REF]. Figure 3 proposes a picture of the bubble plume at 1 s, which shows the existence of small spherical bubbles in red considered as a dispersed field. Figure 4 finally compares at two different heights the void fraction profiles (the same trends are observed on the other heights). The simulation predicts quite well the profile close to the injection but overestimates the void fraction far from it because the bubbles dispersion is slightly under-predicted. The dispersion rate depends on the turbulence effects discussed in section 8.

METERO EXPERIMENT WITH THE MULTIFIELD APPROACH

The test section, 5.40 m long, has an inner diameter D = 0.1 m (Figure 4). The air injection tubes have been set to ensure in the same time uniform bubble injection in the inlet section and low pressure drop for the required gas flowrates. The system also includes a serie of grids that ensures a low turbulence level at the inlet of the test section. Videos of various liquid and mass flowrates at inlet, allowed the building of the METERO [START_REF] Bottin | Experimental investigation of a developing two-phase bubbly flow in horizontal pipe[END_REF] flow pattern map at 40 diameters downstream of the injectors (Figure 5). Calculations have been performed with a three dimensional grid containing about 600000 cells. For all calculations the superficial gas velocity at inlet is J G = 0.1273 m/s.

In figures 7 to 15, the fluid is flowing from the right to the left. Figure 8 -Figure 10 -Figure 12 are the numerical simulations: first and second image show a side view of the dispersed and continuous gas field respectively and the bottom part of the picture gives a view from above of the continuous gas field.

For high values of liquid superficial velocities (J L = 4.42 m/s), a stratified bubbly flow is calculated (Figure 8) in accordance with the METERO flow pattern (Figure 5) and the experimental observations (Figure 7). In the calculations, the volume fraction of the continuous gas field is evaluated but is equal to 0 (Figure 8, bottom part).

When the water flowrate is decreased (J L = 2.12 m/s), the two phase flow enters an intermittent regime: the top bubbles coalesce to form plugs as can be seen in Figure 10 in accordance with the METERO flow pattern (Figure 5) and the experimental observations (Figure 9).

For smaller values of J L (J L = 1.06 m/s), the flow regime completely changes: a free surface is created but as a result of gas injection, Kelvin-Helmoltz instabilities lead the liquid to reach periodically the upper wall, generating a high velocity slug, as can be seen in Figure 12 in accordance with the METERO flow pattern (Figure 5) and the experimental observations (Figure 11).

Not only the multifield model proved to be of relevant interest to simulate flows containing a large range of bubble diameters but also to represent horizontal pipe two-phase flow patterns. This result is crucial in many industrial applications as for nuclear power plant where the two-phase flow occurring in steam generator are not determined up to now but are needed to solve tubes vibration induced by the flow. Moreover, the sensitivity to the turbulence model has been assessed for the plug flow. The Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF]) and the Wall-Adapting Local Eddy-Viscosity (WALE) model [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF] have been used in the simulations (Figure 13 -Figure 14). Even if the nature of the flow regime remains unchanged, the shape of the plug has been modified. The important topic of the turbulence modelling dedicated to two-phase flows is addressed in the section 8.

LIQUID / VAPOUR FLOW IN AN OBLIQUE TUBES BUNDLE GEOMETRY

The flow studied here is a 3D two-phase R114 Freon in oblique tubes bundle geometry (MAXI 2 experiment) [START_REF] Soussan | Local two-phase flow measurement in a oblique tube bundle geometrythe MAXI program[END_REF]. The test section is a tight channel holding tubes bundle with a squared flow section of 97.5 mm. This tank is equipped with 40 rows of 5 tubes (4 tubes plus 2 half-tubes on the wall sides). These adiabatic tubes have an external diameter of 13.5 mm and are inclined of 30° with the horizontal (Figure 15). For CPU time reasons, the 3D mesh used for these computations only include 11 rows of 5 tubes instead of 40 rows in the experiment (Figure 15). Indeed, previous computations [START_REF] Bouillet | Comparison between a simulation and local measurements of a two-phase flow across a 60° inclined tube bundle[END_REF] on this configuration have already shown that 11 rows of tubes are enough to correctly study the flow properties between the 7 th and 8 th row (instead of studying it between the 24 th and 25 th row in the experiment including 40 rows of tubes). The mesh has been generated with a particular attention paid to the mesh quality. The resulting grid is fully hexahedral, without any non-conformity (hanging nodes). The cells characteristic size continuity has also been taken into account to build an as much as possible uniform mesh. A sensitivity study has been done with three different refinements in (Merigoux et al., 2016). In accordance with this study, calculations are performed with the medium grid containing 6 500 000 cells.

RSTM model is used for the continuous fields. Two simulations have been performed: 1. Calculation using the three-field approach developed in the paper and where alp_tot is the sum of the void fraction for the continuous gas ( 2 ) and the void fraction for the dispersed gas ( 3 ). 2. Calculation with the standard two-field model (continuous liquid and dispersed gas fields) called "DISP". Figure 18 shows a comparison between experimental data and CFD results of void fraction and gas velocity profiles measured along the lines NS and WE. Both experimental and computation results have been averaged in time so as to get statistically converged profiles. By the way, it can be noticed that the Freon vapour distribution between the tubes is highly unsteady with the formation of different bubble sizes and shapes, moving along the rows of tubes (Figure 16 and Figure 17).

Globally, void fraction and gas velocity predictions are in good agreement with the experiment. In the NS direction, the maximum values of void fraction (and minimum values of velocity magnitude) behind the cylinders are well captured. Regarding the results obtained in the WE direction, the void fraction is a little bit over predicted in the East area (the vapour stays a little bit too much at the beginning of the tubes instead of sliding along them).

Regarding results obtained with the standard two-field approach (DISP), although this model give a good agreement with experimental data in many test cases [START_REF] Mimouni | A second order turbulence model based on a Reynolds stress approach for two-phase boiling flow. Part 1: Adiabatic cases[END_REF][START_REF] Mimouni | A second order turbulence model based on a reynolds stress approach for two-phase boiling flow and application to fuel assembly analysis[END_REF][START_REF] Mimouni | Combined evaluation of second order turbulence model and polydispersion model for two-phase boiling flow and application to fuel assembly analysis[END_REF], the agreement is quite bad in this case because of the formation of plugs and large bubbles.

LES OF TWO-PHASE FLOWS

The inversed phase benchmark is a phase separation test case, which have been previously used in different publications to study the role and to compare the order of magnitude of the two-phase subgrid terms [START_REF] Ishii | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF], [START_REF] Vincent | Numerical simulation of phase separation and a priori two-phase LES filtering[END_REF] and [START_REF] Larocque | Parametric study of LES subgrid terms in a turbulent phase separation flow[END_REF]. Contrary to the previous studies, the simulation is performed with a two-fluid model using the Large Bubble Model presented and validated above and in previous studies (Denefle,2015), (Fleau, 2015[START_REF] Fleau | Conservative implementation of the interface sharpening equation within an incompressible isothermal multifield approach[END_REF]. In this test case, an oil drop with a cubic shape (size of L/2) is initially placed in a cubic box (size L = 0.1 m) containing liquid water (see Figure 19). The evolution of the system is driven by the gravity forces. At the end of the simulation, the oil phase is supposed to be located in the top part of the box with the liquid water beneath. The fluid properties are: 𝜌 𝑜𝑖𝑙 = 900 kg. m -3 , 𝜇 𝑜𝑖𝑙 = 0.1 Pa. s, 𝜌 𝑤𝑎𝑡𝑒𝑟 = 1000 kg. m -3 and 𝜇 𝑤𝑎𝑡𝑒𝑟 = 1. 10 -3 Pa. s , 𝜎 = 0.045 N. m . The test case is simulated with three different mesh refinements: 128 3 cells, 256 3 cells and 512 3 cells. The time steps are kept constant and are respectively equal to 0.8 ms, 0.2 ms and 0.05 ms. The simulations have been performed with 144 cores for the first mesh and 1152 for the two others during respectively 7 hours, 47 hours and 2 months to reach 13 physical seconds.

Macroscopic behavior

First, we study different macroscopic quantities to validate the Large Bubble Model in this configuration. Thus, the evolution of the normalized kinetics energy (E 𝑐,𝑘 = 𝛼 𝑘 𝑟𝑜𝑡(𝑢 𝑘 )Ω 𝐼 , 𝑟𝑜𝑡(𝑢 𝑘 ) denoted the rotational of the velocity 𝑢 𝑘 ) is examined and reported in (Fleau, 2015). The results obtained with the Large Bubble Model reproduce the same trends observed with other CFD codes [START_REF] Vincent | Numerical simulation of phase separation and a priori two-phase LES filtering[END_REF]. Oscillations are observed for the potential and kinetics energies due to the sloshing motion of oil when it reaches the top of the box. Concerning the enstrophy, the location of the peak is also found to occur at around 3 in dimensionless time.

Filtered two-fluid model equations

As previously done in [START_REF] Ishii | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF] and [START_REF] Vincent | Numerical simulation of phase separation and a priori two-phase LES filtering[END_REF] with the single-fluid model and [START_REF] Lakehal | DNS and LES of turbulent multifluid flows[END_REF] with the two-fluid model including a dispersed field, we apply the LES filter to the twofluid model equations for two continuous fields within the Large Bubble Model.

 The mass balance equation:

𝜌 𝑘 𝜕 𝑡 α k + 𝜌 𝑘 ∇. α k 𝑢 𝑖,𝑘 + 𝜏 𝑖𝑛𝑡𝑒𝑟𝑓 = 0 (11)
With α k the filtered volume fraction of field k and 𝜏 𝑖𝑛𝑡𝑒𝑟𝑓 a subgrid terms related to the relationship between the filtered velocity 𝑢 𝑖,𝑘 and the interface topology (see Table 1).

With 𝜏 𝑐𝑜𝑛𝑣 , 𝜏 𝑑𝑖𝑓𝑓 the convective and diffusive subgrid terms and 𝜏 𝑝𝑟𝑒𝑠𝑠𝑢𝑟𝑒 , 𝜏 𝑠𝑢𝑝𝑒𝑟𝑓 and 𝜏 𝑑𝑟𝑎𝑔 three specific subgrid terms of the two-fluid model applied to two continuous fields (see Table 1 for the expressions).

We can notice that, with the Favre's averaging, three subgrid terms 𝜏 𝑖𝑛𝑡𝑒𝑟𝑓 , 𝜏 𝑡𝑒𝑚𝑝𝑠 and 𝜏 𝑑𝑖𝑓𝑓 are equal to zero. Therefore, the modeling effort concerns less subgrid terms.

-∇. α k 𝑆 𝑘 α k ∇𝑃 -α k ∇𝑃 𝜎 α k κ ∇α k -α k κ ∇α k Favre's averaging - - 𝜌 𝑘 ∇. α k 𝑢 𝑖,𝑘 𝑢 𝑗 ,𝑘 -∇. α k 𝑢 𝑖,𝑘 𝑢 𝑗 ,𝑘 - α k ∇𝑃 -α k ∇𝑃 𝜎 α k κ ∇α k -α k κ ∇α k
Moreover, a small contribution of the diffusion subgrid term 𝜏 𝑑𝑖𝑓𝑓 is observed for the two phases. The Favre's averaging does not affect the classification of the subgrid terms. Finally, as expected, we also observed a decrease of the magnitude of each term when the mesh is refined except for 𝜏 𝑠𝑢𝑝𝑒𝑟𝑓 in the oil phase (Fleau, 2015).

Turbulence model

Finally, we apply five different turbulence models to these subgrid terms: the Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF], the Wall-Adapting Local Eddy-Viscosity (WALE) model [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF], the Bardina's model [START_REF] Bardina | Improved subgrid scale models based on large eddy simulation[END_REF], the mixed Smagorinsky-Bardina's model [START_REF] Bardina | Improved turbulence models based on large eddy simulation of homogeneous, incompressible[END_REF]) and the Adaptative Deconvolution Model (ADM) [START_REF] Adams | A subgrid-scale deconvolution approach for schock capturing[END_REF] with an order of 6.

The relative error of each model with the DNS prediction of each subgrid term is calculated and highlights that the ADM is the most appropriate model for all the subgrid terms (Fleau, 2015) which is the key point of this section. Indeed, Figure 20 displays the dispersion and the slope obtained between each model and the DNS results for the convective subgrid term. ADM is the only model to present a slope close to 1 with a limited dispersion. The main outcome of the paper is the validation of the multifield approach against the Castillejos, METERO and MAXI2 experiments. A reasonable accuracy has been reached and it is worth noting that the flow pattern maps have been successfully predicted by the three-field approach. Regarding the simulation of the MAXI2 experiment, the standard two-field approach (dedicated to dispersed bubbly flows) fails and only the three-field approach can reach a reasonable agreement with the experimental data. Finally, the three-field approach can be seen as a generalization of the standard two-field approach.

Moreover, this article has been partly devoted to the LES study of the two-fluid equations with large interfaces configuration. The equations filtering highlighted new subgrid terms with non negligible orders of magnitude. Finally, different models have been applied and compared. The best correlation with the DNS results has been obtained with the ADM for all the subgrid terms.

Work is still in progress to extend now the model to interfaces with phase changes and to adapt the ADM in the Large Bubble Model.

Figure 2 :Figure 3 :

 23 Figure 2: Large bubbles occurring at different positions in the plume (at 110 mm from the inlet for (a) and (b), 200 mm from the inlet for (c) and (d))

Figure 4 :

 4 Figure 4: Averaged void fractions profiles after 20 s of averaging, the black crosses correspond to the experimental results and the red curve to the simulation. Computational data have been averaged over the last 10 s to get statistically converged profiles.

Figure 5 :

 5 Figure 5: Test section. Figure 6: METERO flow pattern for X/D = 40. Transition from slug to stratified flow (TSS): pink line; transition from plug to slug flow (TPS): orange line; transition from buoyant bubble flow to stratified bubble flow (TBBSB): green line; transition from stratified bubbles regime to plug (TSBP): purple line. Crosses correspond to video acquisitions for couples of liquid and gas superficial velocities JL and JG.

Figure 7 :

 7 Figure 7: Experimental observations: stratified bubbly flow .

Figure 8 :

 8 Figure 8: Simulation of a stratified bubbly flow .

Figure

  Figure 9: Experimental observations: plug flow Figure 10: Simulation of a plug flow

  Initially, the channel is filled with liquid Freon. At the inlet (bottom of the channel), the liquid and vapour (dispersed gas field) Freon are injected with respective velocities of 0.183 m/s and 0.319 m/s and 22% of void fraction. At the outlet (top of the channel), the pressure is set to 8.710 5 Pa. The fluid properties are taken at saturation and no mass transfer occurs. Void fraction and gas velocity are measured along the line NS defined by x = 48.75mm and z = 276.36mm, and the line WE defined by y = 48.75mm and z = 276.36mm (cf Figure15), i.e. between the 7 th and 8 th row tubes.

Figure

  Figure 15: Sketch of the MAXI 2 experiment

Figure 18 :

 18 Figure 18: Void fraction and gas velocity comparisons between experiment and NEPTUNE_CFD results

,

  𝐼 denoted the cell index and Ω the cell volume), potential energy (E 𝑝,𝑘 = 𝛼 𝑘 𝜌 𝑘 gzΩ 𝐼 ) and enstrophy (E 𝑠,𝑘 = 1 2

Figure 19 :

 19 Figure 19: Initial conditions of the inversed phase benchmark.

Figure 20 .

 20 Figure 20. Correlation between the turbulence models and the DNS results for the convective subgrid term for the oil phase

  𝑑 𝑝 is also defined as a function of this product. Close to large interfaces, 𝑑 𝑝 is evaluated by a

	local quantity	𝜇 𝑐𝑔 𝜇 𝑐𝑙	𝛼 𝑝 ||∇𝛼 𝑝
	𝛼 𝑔 < 0.3:		𝐅 𝐛𝐮𝐛𝐛𝐥𝐞 = 𝛼 𝑐𝑙 𝛼 𝑐𝑔	18𝜇 𝑙 𝛼 𝑐𝑙 𝑑 𝑝 2 (𝒖 𝒄𝒍 -𝒖 𝒄𝒈 )
	𝛼 𝑔 > 0.7:		𝐅 𝐝𝐫𝐨𝐩𝐥𝐞𝐭 = 𝛼 𝑐𝑙 𝛼 𝑐𝑔	18𝜇 𝑔 𝛼 𝑐𝑔 𝑑 𝑝 2 (𝒖 𝒄𝒍 -𝒖 𝒄𝒈 )	(5)
	0.3 ≤ 𝛼 𝑔 ≤ 0.7:		𝐅 𝐦𝐢𝐱 =	0.7 -𝛼 𝑐𝑔 0.7 -0.3	𝐅 𝐛𝐮𝐛𝐛𝐥𝐞 +	𝛼 𝑐𝑔 -0.3 0.7 -0.3	𝐅 𝐝𝐫𝐨𝐩𝐥𝐞𝐭
	With the subscript cl referring to a continuous liquid field and cg to the continuous gas field. The
	simulations can be also performed with two continuous liquid fields, replacing the subscript cg by another
	cl subscript.

This drag force is proportional to the product of the volume fractions 𝛼 𝑐𝑙 𝛼 𝑐𝑔 such that it has only non zero values at the interface. Let us note that mass, momentum and energy interfacial transfer term between the phases k and l are proportional to 𝛼 𝑙 𝛼 𝑘 in the models implemented in order to ensure a minimum-maximum principle for all resolved variables (numerical robustness). The characteristic length scale



  The momentum equation: 𝜌 𝑘 𝜕 𝑡 α k 𝑢 𝑖,𝑘 + 𝜏 𝑡𝑖𝑚𝑒 + 𝜌 𝑘 ∇. α k 𝑢 𝑖,𝑘 𝑢 𝑗 ,𝑘 + 𝜏 𝑐𝑜𝑛𝑣 = 𝜇 𝑘 ∇. α k 𝑆 𝑘 + 𝜏 𝑑𝑖𝑓𝑓 -α k ∇𝑃 -𝜏 𝑝𝑟𝑒𝑠𝑠𝑢𝑟𝑒 + α k 𝜌 𝑘 𝑔 𝑖 + 𝐹 𝐶𝑆𝐹 + 𝜏 𝑠𝑢𝑝𝑒𝑟𝑓 + 𝐹 𝐷𝑟𝑎𝑔 + 𝜏 𝑑𝑟𝑎𝑔

Table 1 Subgrid term appearing in the filtered two-fluid equations,

 1 𝐮 𝐢,𝐤 = These subgrid terms are compared in terms of order of magnitude to find the predominant and negligible ones. For this purpose, a top hat filter is applied to the simulation results extracted at the peak of enstrophy for the three grids. Only the first neighborhood of each cell (filter size of 2) is considered to obtain the value of the subgrid terms. Each subgrid term in the momentum equation is normalized by the convection resolved term. The interfacial subgrid term 𝜏 𝑖𝑛𝑡𝑒𝑟𝑓 appearing in the mass balance equation is normalized by its corresponding resolved part. The subgrid term 𝜏 𝑡𝑖𝑚𝑒 is not considered in this study since the analysis is proposed only for one time. Moreover, as documented in Labourasse et al. (2007),[START_REF] Vincent | Numerical simulation of phase separation and a priori two-phase LES filtering[END_REF],[START_REF] Larocque | Parametric study of LES subgrid terms in a turbulent phase separation flow[END_REF], the subgrid term 𝜏 𝑠𝑢𝑝𝑒𝑟𝑓 is predominant for the oil phase.

	𝛂 𝐤 𝐮 𝐢,𝐤 𝛂 𝐤	being the Favre's average of 𝐮 𝐢,𝐤

ACKNOWLEDGMENTS

The NEPTUNE project is funded by EDF (Electricité de France), CEA (Commissariat à l'Energie Atomique et aux Energies Alternatives), AREVA-NP and IRSN (Institut de Radioprotection et de Sûreté Nucléaire).