Paul Chevalier

Patrick Bouchon

Michael Verdun

Emilie Steveler

Nathalie Bardou

Christophe Dupuis

Fabrice Pardo

Riad Haidar

Experimental demonstration of the optical Helmholtz resonance

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Experimental demonstration of the optical Helmholtz resonance

Paul Chevalier, 1, 2 Patrick Bouchon, 1, a) Michael Verdun, 1,2 Emilie Steveler, 2 Nathalie Bardou, 2 Christophe Dupuis, 2 Fabrice Pardo, 2 and Riad Haïdar 1, 3 1) DOTA, ONERA, Université Paris-Saclay, F-91123 Palaiseau -France 2) Centre de Nanosciences et de Nanotechnologies, CNRS, Université Paris-Sud, Université Paris-Saclay, avenue de la Vauve, 91120 Palaiseau, France Optical nanoantennas are studied to manipulate light and enhance light matter interactions. Here, we experimentally demonstrate the optical Hemholtz resonance in a metallic slit-box structure, which is predicted to be harmonic and to enhance the electric field intensity by several orders of magnitude. It is fabricated thanks to a two steps electron beam lithography, between which the box was filled with benzocyclobutene (BCB). Up to 80 % of the light is absorbed at a λ = 2.84 µm wavelength under a beam focused by a Cassegrain objective (NA=0.4), even if the dimensions of this resonator are deeply subwavelength for both the slit (width λ/55 and height λ/77) and the box (width λ/7 and height λ/37). As expected from the inductance nature of the box, the optical properties of the BCB filling the box has no influence on the resonance behavior.

Optical nanoantennas are used to manipulate the light at a subwavelength scale [1][2][3][4] . They are able to generate strong field intensity enhancement in hot spots or hot volumes 5 , which is an appealing feature for light matter interactions like non linear generation 6 , detection 7 or thermal emission 8 . Several nanoantennas designs have been inspired by their radiofrequency alter egos, e.g. the Yagi-Uda antenna 9 , the bow-tie antenna 10 , the patch antenna 11 . In these three examples, one of the dimension of the antenna is directly linked to its resonance wavelength, and these antennas exhibit harmonic resonances at lower wavelengths.

In contrast, the optical Helmholtz resonator was named by analogy with the 1D acoustic Helmholtz resonator [START_REF] Helmholtz | Theorie der Luftschwingungen in Röhren mit offenen Enden[END_REF] . It consists of a metallic structure made of a box which is coupled to the semi infinite medium through a slit aperture [START_REF] Chevalier | [END_REF] . Arrays of slits or boxes have been widely studied and are known to be the siege of Fabry Perot resonances 14 , and likewise resonances are predicted to occurred in either the cavity or the aperture of the optical Helmholtz resonator. However, the Helmholtz resonance is predicted to appear at a higher wavelength than Fabry Perot resonances, which depends on each geometric parameter, and to be harmonic since it relies on an LC resonance mechanism [START_REF] Chevalier | [END_REF] . It was suggested to be used for dual band absorption for infrared stealth applications 15 . Moreover, this structure has also been predicted theoretically to provide a strong field intensity enhancement 16 , which can be engineered from the near infrared to the THz wavelengths, and used for instance to enhance surface enhanced infrared absorption (SEIRA) 17 . Yet, there has been no experimental evidence of the Helmholtz resonance in optics, certainly due to the dimensions and geometry that are rather challenging to obtain.

In this letter, we report the experimental demonstration of the optical Hemholtz resonator in the mid-infrared a) Electronic mail: patrick.bouchon@onera.fr range. We developed a process for fabricating the resonator using two levels of electron beam lithography, and a planarizing dielectric for filling the box as an intermediate step. The reflectivity measurement in TM polariza- tion under light focused by a Cassegrain objective showed a strong resonant absorption at the theoretically predicted wavelength, and a fair agreement is obtained between the experimental and the computed spectra. The slit is not perfectly centered in the middle of the box, and we show that this alignment offset induces an additional resonance at a lower wavelength.

(a) (b) (d) (c)
The general shape of the optical Helmholtz resonators is recalled in Fig. 1(a). It is made of a box (height h b , width w b ) closed by a narrow slit (height h s , width w s ). In the following, the Helmholtz structure studied is long enough along the y-axis to be considered as infinite, and is repeated periodically along the x-axis with a period d. We choose to fill the cavity with benzocyclobutene (BCB), a planarizing dielectric, as shown in Fig. 1(a) to reduce the complexity of fabrication by avoiding freestanding parts, and to use gold as the metallic material for the rest of the structure. This structure was simulated using an exact Maxwell solver 18,19 , and the dielectric function of gold was computed using a Drude model which fits well the experimental data in the infrared 20 :

ε Au (λ) = 1 -1/[λ p /λ(λ p /λ + iγ)
] with γ = 0.055 and λ p = 159 nm. The BCB was described by a Lorentz model fitting experimental data 21 . The computed absorption spectrum, for a transverse magnetic (TM) polarized and normally incident light, is plotted in Fig. 1(b) for a given geometry (w b = 400 nm, h b = 77 nm, w s = 52 nm, h s = 37 nm, d = 450 nm ). At λ = 2.53 µm, 36 % of the light is absorbed due to the Helmholtz resonance. Other resonances appear at lower wavelengths between 500 nm and 700 nm. The small peak at λ = 3.4 µm is due to a vibrational absorption of BCB that is enhanced by the Helmholtz resonator 21 . For the sake of comparison, the absorption spectra are also plotted under the same illumination conditions for an array of slits (height 4h s , width w s ) and for an array of cavity filled with BCB (width w b , height h b ) that also exhibit resonant absorption, but at lower wavelengths.

The key properties of the Helmholtz resonance, that distinguishes it from others, are: (i) A giant field enhancement, typically two orders of magnitude higher than expected with gap plasmons resonances. (ii) A single resonance behavior with no harmonics. (iii) No influence of the permittivity of the material filling the cavity on the resonance wavelength. (iv) The resonance wavelength λ r depends on the square root of all four geometrical parameters according to the following formula:

λ r ∝ 2π ε s h b w b h s w s , (1)
where ε s is the permittivity of the material filling the slit. Noteworthily, the permittivity of the material filling the cavity has no influence on the resonance wavelength. We now consider an Helmholtz structure with a larger period d = 1.5 µm, in order to be closer to the critical coupling of the resonator, and the following goal dimensions: h b = 80 nm, w b = 400 nm, h s = 40 nm, w s = 50 nm.

The fabrication process starts with the deposition of 20 nm of chrome followed by 200 nm of gold on a doped silicon wafer. A first step of e-beam lithography is performed to define the cavity using ma-N 2405 resist spun at 3000 rpm on top of the gold layer. A grating pattern is written on a 0.5 × 0.5mm 2 area thanks to a Vistec EBPG 5000 electron beam pattern generator at 100 kV. The resist is developed using a tetramethyl ammonium hydroxide (TMAH) based developer (Clariant AZ 826 MIF) for 30 seconds. Then 80 nm gold is evaporated onto the sample and lifted-off with acetone (see Fig. 2(a)) to fabricate trenches of height 80 nm and width 400 nm. Benzocyclobutene (BCB) is spin coated on the sample and baked at 250 • C for 3 hours under nitrogen flow. A reactive ion etching (RIE) plasma (SF 6 + CHF 3 ) is done so that the surface of the BCB reaches the gold surface while keeping the cavity filled with dielectric (see Fig. 2(a)). A scanning electron microscope (SEM) image of the structure at this step is shown in Fig. 2(b), showing the gold trenches (light grey) filled with a dielectric (black). On top of theses trenches, 150 nm of polymethyl methacrylate (PMMA) is spin coated and baked at 160 • C degrees for 30 minutes and 60 nm of hydrogen silsesquioxane (HSQ) in spin-coated and baked on a hot plate at 110 • C for 20 minutes. A 20 nm pattern is written at 100 kV. The resist is developed for 10 seconds in a potassium hydroxide (KOH) based developer (Clariant AZ 400K diluted at 50%) and rinsed in distilled water. The PMMA is etched with a RIE O 2 plasma for 60 seconds to create a tall and thin wall of polymer centered on the cavity (see Fig. 2(c)). Finally 50 nm of gold is evaporated and lifted-off using trichloroethylen at 90 • C in ultrasound bath, which terminates the process. A SEM image showing the slits from the top of the structure at this final step is shown in Fig. 2(d), on this image a slight variation of the contrast allow to distinguish the position of the underlying cavity. The choice of a bi-layer resist stack for the final lift-off of the slit is motivated by the need for a fine resolution of the pattern; also the RIE etching step using the HSQ as a mask allows to create an undercut of the PMMA resist 22 which eases the following lift-off. The second level is aligned on the lower level with alignment markers positioned during the first step of lithography and detected by the e-beam pattern generator. The accuracy of this alignment on the sample size is about 50 nm, and a slight shift of the slit in regard to the cavity can be expected.

The reflectivity of the fabricated sample was measured using a Bruker Vertex 70 Fourier Transform Infrared (FTIR) spectrometer, with the Bruker Hyperion 2000 microscope using a Newport Cassegrain reflective objective (15x and NA 0.4). The measurements were done using the internal mid-infrared source of the spectrometer and an InSb photodetector. A square area of 100 µm×100 µm is imaged with the microscope using knives and the beam is linearly-polarized with the magnetic field along the slits using an IR polarizer. In the following, we define the magnetic parallel polarization, for which the magnetic field is in the plane parallel to the orientation of the grating (plane yz, i.e., H x = 0), to avoid any confusion with the TM polarization. The light is in the magnetic parralel polarization after going through the polarizer and the Cassegrain objective, and is exciting the TM modes of the metallic structure. A schematic of the beam on the sample surface recalling the polarization is shown in Fig. 3(a).

The magnetic-parallel polarized measured reflectivity of this sample, plotted as a function of the wavelength in Fig. 3(b) (blue continuous line), exhibits two resonance wavelengths at λ r = 1.56 µm and λ r = 2.82 µm. It should be underlined that light is incoming on a gold film patterned on less than 4% of its area, yet more than 70% of the incoming light is absorbed at the higher resonance wavelength. The simulated response of a Helmholtz resonator is also plotted (orange dashed line) with the following dimensions: h b = 77 nm, w b = 400 nm, h s = 37 nm, w s = 52 nm, and the slit was offset by 51 nm compared to the center of the box. The values taken for the box and slit widths, as well as the slit offset, are coherent with the SEM measurements (see Fig. 2). The simulated curve also corresponds to the integration of the reflectivity over various incidences: polar angles θ from 12 • to 23.6 • and azimuthal angles φ from 0 • to 360 • to take into account the microscope objective. A rather fair agreement is obtained between the computed and experimental spectra. As expected from Eq. 1, the resonance wavelength is close to the one predicted in Fig. 1, although it is red-shifted. However, the level of reflectivity is different. This is due to both the change of the period that modifies the coupling factor with the continuum (i.e., the level of absorption at resonance), and the change of the incidence configuration. Eventually, there is the apparition of a second resonance that can be attributed to the slit offset. In the following, we theoretically investigate the influence of both the focalization and of the slit offset on the response of the Helhmoltz resonator.

To understand how the integration of the response over various incidences influence the resonance, we plotted on the polar diagram in Fig. 3(c) the reflectivity response at the resonance wavelength (λ = 2.82 µm) as a function of the conical angle of incidence. The Helmholtz resonance is localized, and for θ < 25 • , the absorption is higher than 75% for any incidence. It should be noted that for the polar angle values used by the Cassegrain objective, the azimuthal angle has a slight influence on the actual reflectivity, but this impacts mainly the quality factor of the resonance. Noteworthily, it has been shown that the angular tolerance of the Helmholtz resonator can be improved by inserting a high refractive index material inside the slit 23 . In order to study the consequence of the misalignment between the slit and the box, we plotted in Fig. 4 the reflectivity spectrum of a Helmholtz structure with the slit offset. In the case of a perfectly centered slit, the reflectivity exhibits only a single resonance peak. As the slit moves away from the center of the structure, the main resonance is red-shifted, while a secondary resonance appear. The secondary resonance in the experimental curve of Fig. 3(b), and the red-shift of the Helmholtz resonance are indirect proofs that the slit was slightly misaligned during the fabrication process and it justifies the offset value used for the fit.

To conclude, a process to fabricate optical Helmholtz resonators has been presented, and used to experimentally demonstrate the optical Helmholtz resonance thanks to measurements in the near and mid-infrared were measured. A secondary resonance was measured at a lower wavelength, but it is explained by a defect on the symmetry of the resonator. The agreement between the theoretical and experimental far field optical properties of the resonators is a first step towards the experimental demonstration of the high field intensity enhancement 16 . An experimental proof could be given by measuring the enhancement of light matter interactions inside the slit, either by putting in it molecules with absorption features (SEIRA) or a nonlinear material.

3 )

 3 École Polytechnique, Département de Physique, 91128 Palaiseau, France (Dated: 11 April 2018)

FIG. 1 .

 1 FIG. 1. (a) General schematic of the optical Helmholtz resonator, made of a box and a slit in metal, both infinite along the y direction. (b) Computed reflectivity spectrum of the Helmholtz resonator array (period d = 450 nm) made of a slit (width ws = 52 nm, height hs = 37 nm) and a box (width w b = 400 nm, height h b = 77 nm), compared with the reflectivity spectra of the array of slits (height 4hs, width ws, orange curve) and boxes (height h b , width w b ,green curve). The impinging light is normally incident and TM-polarized.

FIG. 2 .

 2 FIG. 2. Fabrication steps for a periodic grating of optical Helmholtz resonators with the corresponding top view SEM image: (a) Gold trenches are created by the lift-off of 80 nm of gold after a first e-beam lithography (d = 3 µm).(b) A dielectric (BCB) is spin coated on the surface and etched by RIE to reach the gold surface (d = 3 µm). (c) A two layer lithography process is performed using PMMA and HSQ (d = 1.5 µm). (d) Gold is lifted-off resulting in the fabrication of the optical Helmholtz resonator (d = 1.5 µm). The narrow slits in the gold surface above the BCB-filled cavity, which is barely distinguishable.

FIG. 3 .

 3 FIG. 3. (a) Schematic representation of the experimental measurement of reflectivity: the incident beam is polarized and focused by a 15x Cassegrain objective and the area of interest limited to a square area of 100 µm × 100 µm. (b) Reflectivity measurement of the structure with h b = 77 nm, w b = 400 nm, hs = 37 nm, ws = 52 nm (blue line) under light magnetic-polarized focused by a Cassegrain objective with NA=0.4. The simulated response is obtained by integrating the polar angle θ from 12 • to 23.6 • and the azimuthal angle φ from 0 • to 360 • (orange dashed line). (c) Computed angular diagram at the resonance wavelength λ = 2.84 µm for a magnetic parallel polarized wave.

FIG. 4 .

 4 FIG. 4. Computed reflectivity diagram as a function of the wavelength and of the offset of alignment between the slit and the box. The reflectivity is computed for a TM-polarized normally incident plane wave.

ACKNOWLEDGMENTS

We acknowledge financial support from the ONERA through the MOLIERE project and from the Agence Nationale de la Recherche (ANR) (ANR-16-CE39-0007).


S Vincent

M Tavares

Solène Fleau

S Mimouni

M Ould-Rouiss

J.-L Estivalezes

A priori filtering and LES modeling of turbulent two-phase flows Application to phase separation

Keywords: phase inversion, two-phase flows, a priori filtering, LES modeling, one-fluid model, multi-field approach

, Wale [6], Bardina [7], Mixed [8] and ADM [9] are used to account for two-phase subgrid effects. These models are compared to filtered DNS results.

Introduction

Turbulent two-phase flows occur in many academic and applied fluid mechanics problems such as boiling crisis in nuclear plants, chemical reactors, material coating by plasma projection, bubbles in pipes, wave breaking, oil extraction in porous media, powder and fluidized beds processes, fuel injection in engines to cite a few. Exhaustive presentation of various turbulent multi-phase flows is given for example in [START_REF] Crowe | Multiphase flow handbook[END_REF]. As soon as turbulence and interface interact in a non linear way with macroscopic interfacial deformations inducing ligaments, coalescence or rupture, the experimental characterization of these flows is difficult due to the heterogenous character of the multi-phase medium. Modeling and numerical simulation thus represent an interesting way to study the physical processes that control these flows.

Figure 1: Example of multi-scale interfacial character of phase separation for a liquid-liquid two-phase system. The interface is plotted in blue and the black lines represent the wall of the cubic box. Following [START_REF] Larocque | Parametric study of LES subgrid terms in a turbulent phase separation flow[END_REF][START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF][START_REF] Vincent | Test-case n • 15: phase inversion in a closed box (pc)[END_REF], Re=7000 and We=11.

Various models are available in literature to tackle with two-phase flows. For dispersed flows, when a scale separation can be assumed between the car-rier fluid and the dispersed phase, i.e. droplets or bubbles, Eulerian-Eulerian or Eulerian-Lagrangian models exist [START_REF] Kleinstreuer | Two-phase flow: theory and applications[END_REF] in which the continuous phase, also called the carrier fluid, is treated with the classical fluid mechanics equations and the dispersed phase is taken into account through an Eulerian concentration variable or discrete Lagrangian positions and velocities. Specific source terms are added to the mass and momentum equations in order to represent the interaction effects between both phases. When the characteristic size of the interfacial structures is comparable to the macroscopic size of the problem or to the local mesh size, the previous scale separation models are no more valid.

The present work is concentrated on two-phase flow problems involving a wide range of interfacial structures that can be large and small at the same time, depending on the zone of the flow that is considered. Among others, examples are the primary atomization of a liquid jet or sheet [START_REF] Lebas | Numerical simulation of primary break-up and atomization: DNS and modelling study[END_REF][START_REF] Zuzio | Numerical simulation of primary and secondary atomization[END_REF][START_REF] Trontin | Direct numerical simulation of a freely decaying turbulent interfacial flow[END_REF] or the liquid-liquid phase separation [START_REF] Larocque | Parametric study of LES subgrid terms in a turbulent phase separation flow[END_REF][START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF][START_REF] Vincent | Test-case n • 15: phase inversion in a closed box (pc)[END_REF] presented in figure 1. As a starting point, we consider here multi-scale two-phase flows with separate phases in unsteady and turbulent regime. It is assumed that no scale separation exists between large interfacial structures and small turbulent scales. For large Reynolds numbers, performing a Direct Numerical Simulation (DNS) of the flow, i.e. a simulation with numerical time and space scales that are smaller than the physical characteristics of turbulence and interface, is almost impossible as it would require computer resources not available even on the most powerful parallel computers. As a consequence, for a kind of mesoscopic modeling issue, it has to be assumed that all the interfacial scales will be resolved whereas turbulence will be modeled by means of a Large Eddy Simulation (LES) approach [START_REF] Sagaut | Large Eddy Simulation for incompressible flows -An introduction[END_REF]. This way, relatively large grids and time steps could be considered that give access to deterministic and unsteady simulations on meshes coarser than with a DNS. However, two basic points have to be specified to lead our mesoscopic representation of the turbulenceinterface interaction: on one hand, the choice of models that allow the solving of deformable interfaces between non miscible fluids and on the other hand, the LES filtering of these models in the framework of two-phase flows. These two points are the main theme of our work. On a general point of view, one can adress the following remarks in order to build a LES two-phase flow modeling with resolved scale interfaces:

• concerning the modeling of two-phase flows with separated phases and non miscible fluids, the most widely used model in the literature is the single fluid or one-fluid model (OFM) [START_REF] Kataoka | Local instant formulation of two-phase flow[END_REF][START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF][START_REF] Prosperetti | Computational methods for multiphase flows[END_REF]. This Eulerian model is based on the use of a regular kernel of control volumes on which the in-terface is located by means of an auxiliary Eulerian variable. This variable can be the local Volume fraction of one phase in the other phase, i.e. the Volume Of Fluid (VOF) [START_REF] Hirt | Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries[END_REF], the signed distance function to the interface, i.e. the Level Set approach [START_REF] Osher | Level set methods: An overview and some recent results[END_REF], or the VOF function resulting from the projection of a Lagrangian tracking of the interface, i.e. the Front Tracking method [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]. In these three approaches, the standard incompressible Navier-Stokes equations are considered, with an additive specific capillary source term that accounts for the normal jump of constraints at the interface. In the VOF and Level Set, an additional advection equation of the interface is solved in an Eulerian way while a Lagrangian interface marker equation is considered with the Front Tracking. The one-fluid model allows the representation of twophase flows with separated fluids and deformable interfaces. Another class of models can be considered to have the same kind of two-phase flow representation. It is called the multi-field approach (MFA) [START_REF] Denèfle | Multifield hybrid approach for two-phase flow modeling -part 1: Adiabatic flows[END_REF][START_REF] Fleau | Simulations of twophase flows with a multifield approach[END_REF][START_REF] Fleau | Validation of a multifield approach for the simulations of two-phase flows[END_REF] or segregated phase modeling [START_REF] Prosperetti | Computational methods for multiphase flows[END_REF]. This type of two-phase flow modeling was initially designed for a dispersed phase that is represented through an Eulerian variable such as the local concentration of the dispersed phase in the surrounding carrier fluid. However, this Euler-Euler model can be extended to large scale and resolved interfaces by introducing a sharpening equation together with the transport of the dispersed phase concentration [START_REF] Denèfle | Multifield hybrid approach for two-phase flow modeling -part 1: Adiabatic flows[END_REF][START_REF] Fleau | Validation of a multifield approach for the simulations of two-phase flows[END_REF]. This sharpening equation is coupled with mass and momentum equations in order to ensure mass conservation to almost zero computer error. The obtained multi-field model allows simulating large scale interfaces in a similar approach as phase field models [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF], with an interface that is given a prescribed thickness of several grid cells. The OFM and MFA will be considered in the present article to discriminate the modeling strategy that is most suited for the simulation of turbulent two-phase flows with large interfaces.

• regarding the LES modeling of turbulence, this is an area of fluid mechanics that has been widely studied and developed for single-phase flows past twenty years [START_REF] Sagaut | Large Eddy Simulation for incompressible flows -An introduction[END_REF][START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF]. For cons, the LES of two-phase flows with large scale interfaces remains an understudied area. Among the few existing references include the work on LES formulation for the one fluid model and subsequent a priori LES filtering of DNS simulations for the estimate of new specific two-phase contributions in the mass, momentum and VOF advection equations [START_REF] Larocque | Parametric study of LES subgrid terms in a turbulent phase separation flow[END_REF][START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF][START_REF] Liovic | Multi-physics treatment in the vicinity of arbitrarily deformable gas-liquid interfaces[END_REF][START_REF] Liovic | Interface-turbulence interactions in large-scale bubbling processes[END_REF][START_REF] Chesnel | Subgrid modeling of liquid atomization[END_REF][START_REF] Labourasse | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF][START_REF] Toutant | Dns of the interaction between a deformable buoyant bubble and a spatially decaying turbulence: a priori tests for les two-phase flow modelling[END_REF]. Two main conclusions arise from these studies: new subgrid terms have to be modeled when the LES of two-phase flows is undertaken and the magnitude of these subgrid terms strongly varies depending on the considered physical configuration so that none of these subgrid terms can be neglected a priori, except perhaps the viscous subgrid term. Concerning the modeling of the two-phase LES subgrid terms, most advanced works have been devoted to the inertial and capillary terms with formulations based on Smagorinsky or deconvolution like approaches [START_REF] Chesnel | Subgrid modeling of liquid atomization[END_REF][START_REF] Aniszewski | A new approach to sub-grid surface tension for les of two-phase flows[END_REF]. None of these LES models have clearly demonstrated that it was possible to perform a LES of resolved scale interfaces with success, except perhaps [START_REF] Chesnel | Subgrid modeling of liquid atomization[END_REF] for the atomization of a fuel jet. In addition, no LES model have been proposed for the viscous, advection or mass conservation subgrid terms.

If one refers to the previous literature review, it can be noted that there is a real need to make a comprehensive and systematic a priori study of the two-phase LES subgrid terms appearing in the one-fluid and even more multifield formalisms. In addition, representative models have to be proposed for all subgrid terms and validated against DNS of turbulent two-phase flows. These two lines of research are the main goal of the present work. The article is structured as follows. In the next section, the OFM and MFA are shortly presented with associated numerical methods. A phase inversion benchmark configuration is also proposed for the LES modeling of two-phase flows with multi-scale interfaces. Section 3 is devoted to a priori LES filtering of OFM and MFA as well as comparison of subgrid term magnitude and classification. LES models for two-phase subgrid terms of OFM and MFA is proposed in section 4. These models are compared and discussed regarding DNS of phase inversion. Conclusions and perspectives are finally drawn in section 5.

Modeling and simulation of phase inversion 2.1 Phase inversion benchmark

Various configurations of phase inversion benchmark have previously been proposed and used for a priori LES filtering, comparison of interface tracking methods, characterization of multi-scale interfacial flow and comparison of DNS two-phase flow codes [2-4, 11, 33, 34]. A previous configurations was chosen here for its turbulent regime and large range of interfacial scales. The initial condition of the problem is described in Figure 2. An initial cubic blob of light liquid, referred to as fluid 1, is placed in the bottom part of a cubic box filled with a heavier liquid, referred to as fluid 2. The size of the box is (H, H, H), while the size of the blob of light fluid is (H/2, H/2, H/2), with H = 1m. All outer walls are considered as free-slip impermeable walls. The gravity is oriented along the vertical axis and its magnitude is -9.81m.s -2 . The interest of the present configuration lies in simplicity of initial and final flow configurations, i.e. all light fluid is inside a cubic blob initially and all light fluid is in the top part of the cavity at long times. Another advantage is that phase inversion naturally provides multiple coalescence, break-up and fragmentation events in a kind of atomization process, while not depending on complex initial or forcing conditions such as those encountered in real turbulent two-phase flow configurations [START_REF] Lebas | Numerical simulation of primary break-up and atomization: DNS and modelling study[END_REF] [START_REF] Tomar | Multiscale simulations of primary atomization[END_REF]. The characteristics of fluid 1 and 2 are ρ 1 = 900kg.m -3 , ρ 2 = 1000kg.m -3 , µ 1 = 0.1P a.s and µ 2 = 0.001P a.s. The constant surface tension coefficient is chosen equal to 0.45N.m -1 .

The fluids and flow properties that control the phase inversion dynamics can be expressed in terms of Atwood, Reynolds and Weber numbers. These numbers are defined as

At = ρ 2 -ρ 1 ρ 1 + ρ 2 , (1)
Re 1 = ρ 1 HU g µ 1 and Re 2 = ρ 2 HU g µ 2 , (2)
W e 1 = ρ 1 HU 2 g σ and W e 2 = ρ 2 HU 2 g σ , (3)
where the gravitational velocity U g is obtained as the balance between the pressure forces on the light fluid and the net gravity force on this fluid [START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF][START_REF] Vincent | A phase inversion benchmark for multiscale multiphase flows[END_REF]. It reads

U g = ρ 2 -ρ 1 ρ 1 gH 2 = At 1/2 ρ 1 + ρ 2 ρ 1 gH 2 (4)
Based on U g and a characteristic space scale of H/2, the characteristic time scale is

t c = H 2U g = At -1/2 ρ 1 (ρ 1 + ρ 2 ) H 2g . (5)
The characteristic parameters of phase inversion are At = 0.053, U g = 0.246m.s -1 , Re 1 = 2.214 • 10 3 , Re 2 = 2.460 • 10 5 , W e 1 = 1.210 • 10 2 , W e 2 = 1.345 • 10 2 , and t c = 2.032s. An example of simulation with these parameters is illustrated in figure 3 for a 512 3 grid with a one-fluid model. Due to large Reynolds and Weber numbers, the flow is first dominated by inertial effects induced initially by gravity forces with negligible capillary forces. For larger times, inertia decreases and surface tension forces become predominant. It is observed that initially the light fluid blob goes to the top of the cavity while being subjected to large interface deformation with rupture of ligaments into droplets. In a second step, a kind of atomization process occurs with the phase separation being active inducing coalescence of many droplets with the top layer of light fluid.

One-fluid model (OFM) 2.2.1 A single velocity formulation for two-phase flows

As is now well established, incompressible two-phase flows involving fluidfluid interfaces and Newtonian fluids can be modeled by a single set of incompressible Navier-Stokes equations with variable density and viscosity and possibly extra interfacial forces (e.g. the capillary force), together with the transport equation of the phase function C. The resulting model takes implicitly into account the mass and momentum jump relations at the interface [START_REF] Delhaye | Jump conditions and entropy sources in two-phase systems. Local instant formulation[END_REF] [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF], whereas the continuity of the fluid-fluid and fluid-solid interfaces are taken into account by the C equation. The entire set of equations reads:

∇ • u = 0 , (6)
ρ( ∂u ∂t + u • ∇u) = -∇p + ρg + ∇ • µ ∇u + ∇ t u + F st , (7)
∂C ∂t + u • ∇C = 0 , ( 8
)
where p is the pressure, F st is the interfacial force per unit volume and ρ and µ are the local density and viscosity of the two-phase medium, respectively. In the rest of the article, the viscous stress tensor will be referred to as S = 1 2 ∇u + ∇ t u . Motion of the interface in a given velocity field u is considered in equation ( 8) through the VOF function C. In the present work, it is chosen instead to use the signed distance function to the interface φ, that satisfies the following equation:

∂φ ∂t + u • ∇φ = 0 (9)
Capillary effects are inserted in the source term F st in the form F st = σκn i δ i , where σ denotes the surface tension, κ is the local mean curvature of the interface, i.e. κ = -∇ • n = -∇ • ∇C ∇C , n i is the unit vector normal to the interface and δ i is the interface Dirac function [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF], such that

n i δ i = ∇C.
The above OFM is almost identical to the classical incompressible Navier-Stokes equations, except that it involves an extra interfacial force and the local properties ρ and µ of the equivalent fluid depend on C, as does the capillary force. Localizing the interface requires solving the additional advection equation for the Level Set φ such that φ = 0 at the interface and φ > 0 (resp. < 0) in fluid 2 (resp. 1). In this case, C is defined as the Heaviside function C = H(φ) [START_REF] Osher | Level set methods: An overview and some recent results[END_REF]. It is assumed that all space and time scales present in the flow are resolved, so that no turbulence model is required.

Numerical methods

The in-house computational fluid dynamics library developed at ONERA called DyJeAT (Dynamic of Jet ATomization) is used. Velocity/pressure coupling is ensured with classical projection methods [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] [39] on a staggered grid. The Ghost Fluid method [START_REF] Osher | Level set methods: An overview and some recent results[END_REF][START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterials flows (the ghost fluid method)[END_REF] is implemented to deal with surface tension forces as well as density and viscosity jumps. Instead of inserting a surface tension volume force in the momentum equations as was proposed by [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF], with the Ghost fluid method, the jump relations are directly used in prediction and projection steps so as to handle density and viscosity contrasts at the interface as well as capillary effects. In the LES formalism, the delta formulation of [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] is kept only to provide an interpretation of the under resolved capillary forces. The Level Set approach is also used for tracking interfaces, as previously specified. To avoid singularities in the φ field, the fifth order conservative WENO [START_REF] Shu | Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws[END_REF] scheme is applied to discretize space advection terms together with a 3 rd order Runge-Kutta TVD scheme for time derivatives [START_REF] Shu | Efficient implementation of essentially nonoscillatory shock-capturing schemes[END_REF]. When the Level Set advection is solved, high velocity gradients can generate artificial spreading or stretching of φ which then no longer remains a distance function. A redistancing algorithm [START_REF] Sussman | Improved level set method for incompressible two-phase flows[END_REF] is thus considered to restore the distance property of φ, i.e. ∇φ = 1. Details of implementation and validation can be found in [START_REF] Trontin | Direct numerical simulation of a freely decaying turbulent interfacial flow[END_REF][START_REF] Trontin | Detailed comparisons of front-capturing methods for turbulent two-phase flow simulations[END_REF][START_REF] Zuzio | An efficient block parallel AMR method for two phase flow simulations[END_REF].

2.3 Multi-field approach (MFA)

A two-velocity model for two-phase flows

The two-fluid model of Ishii [START_REF] Ishii | Thermo-fluid dynamics, theory of two-phase flow[END_REF] is used for a two-phase flow. Initially, the Ishii Eulerian model for two-phase flows was considering the interaction between a carrier fluid (field 1) and a dispersed phase (field 2) under scale separation assumption. In our approach, the density, the viscosity and the local velocity are defined for each field as for Ishii model. However, field k = 1 corresponds to the first fluid and k = 2 to another fluid that is immiscible inside field 1. This study is restricted to incompressible and isothermal cases with a constant density in each field. The following governing equations are solved for each field k:

∂α k ρ k ∂t + ∇ • (α k u k ρ k ) = 0 , (10)
∂α k ρ k u k ∂t + ∇ • (α k ρ k u k ⊗ u k ) = -α k ∇p k +α k ρ k g + ∇ • [α k µ k S k ] + F st,k + F d,k , (11)
With α k the local volume fraction of field k, S k the viscous stress tensor of

field k, i.e. S ij,k = ∂u i,k ∂x j,k + ∂u j,k ∂x i - 2 3
∂u i,k ∂x i δ ij and F d,k a volume force ensuring
the coupling between the continuous fields through a drag force. This drag force is crucial since it allows coupling the velocity of the two continuous fields at the interface. Thus, a new drag force expression has been developed to deal with large interfaces within the MFA. Details can be found in [START_REF] Fleau | Validation of a multifield approach for the simulations of two-phase flows[END_REF][START_REF] Mimouni | Multifield approach and interface locating method for two-phase flows in nuclear power plant[END_REF]. The modeling of large interfaces and surface tension forces with Ishii model is based on a volumetric formulation [START_REF] Denèfle | Multifield hybrid approach for two-phase flow modeling -part 1: Adiabatic flows[END_REF][START_REF] Fleau | Validation of a multifield approach for the simulations of two-phase flows[END_REF][START_REF] Bartosiewicz | A first assessment of the NEPTUNE CFD code: Instabilities in a stratified flow comparison between the vof method and a two-field approach[END_REF] taken from [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]. It reads

F st,k = α k σκ k n i,k δ i with κ k = -∇ • n k = -∇ • ∇α k ∇α k and n i,k δ i = ∇α k .

Numerical methods

The code N EP T U N E CF D is based on an Eulerian approach with a finite volume discretization. In the code N EP T U N E CF D, the assumption of a common pressure for all fields is made. The solver SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) is used [START_REF] Patankar | A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows[END_REF] with a collocated arrangement for all variables. An iterative coupling of the equations is implemented to ensure mass conservation. Details are given in [START_REF] Fleau | Validation of a multifield approach for the simulations of two-phase flows[END_REF][START_REF] Fleau | Conservative implementation of the interface sharpening equation within a multifield approach[END_REF].

By nature, the Ishii model is diffuse as it represents the interaction between two phases (or fields) through a mixed Eulerian representation of the different fluids in a given volume. Here, the idea is to tackle with sharp interfaces, even if the MFA is a priori not adapted to a discontinuous representation of a two-phase flow. In order to limit the interface smearing induced by the resolution of the two-fluid equations, an interface sharpening equation, initially proposed by [START_REF] Olsson | A conservative level set method for two phase flow[END_REF] and adapted to the two-fluid formulation, is solved for each continuous field [START_REF] Fleau | Conservative implementation of the interface sharpening equation within a multifield approach[END_REF]:

∂α k ∂t + ∇ • (α k [1 -α k ] n i,k ) = ǫ∇ 2 α k (12)
With n i,k the interface normal vector pointing outside phase k. The value of ǫ is chosen equal to half the space step ∆x while the time step of the interface sharpening equation is taken equal to ∆x/32. With these parameters, the obtained final interface thickness is always equal to 5 cells whatever the initial interface diffusion [START_REF] Denèfle | Multifield hybrid approach for two-phase flow modeling -part 1: Adiabatic flows[END_REF][START_REF] Strubelj | Numerical simulations of stratified two-phase flows with twofluid model and interface sharpening[END_REF].

Model comparison on macroscopic quantities

Several macroscopic quantities characterizing the evolution of the flow field are of primary interest to check whether the two models describe the same two-phase motion. These quantities can be expected to be a priori essentially independent of the choice of the model and interface tracking techniques used in the computations. In all of what follows, we use the characteristic or color function C that is defined with respect to fluid 2, so that C = 1 in fluid 2 and C = 0 in fluid 1. Among these physically relevant quantities and in agreement with previous work on the phase inversion benchmark [START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF][START_REF] Vincent | A phase inversion benchmark for multiscale multiphase flows[END_REF], we select:

• the domain integrals of potential energies E p,2 = Ω Cρ 2 gydV and E p,1 = Ω (1-C)ρ 1 gydV in fluid 1 and 2 (with g = g ) and the domain integrals of kinetic energies E k,2 = 1 2 Ω Cρ 2 u 2 dV and E k,1 = 1 2 Ω (1 -C)ρ 1 u 2 dV will also be stored at each time step. They will help monitoring the conversion of potential energy into kinetic energy. they provide a sort of characterization of the stratification inside the fluids. Due to the simple topology of the interface in the initial and final stages, it can be shown that E p,1 = ρ 1 g H 4 32 and E p,2 = ρ 2 g 31H 4 32 for t = 0 and E p,1 = ρ 1 g 15H 4 128 and E p,2 = ρ 2 g 49H 4 128 for t → ∞. The typical magnitude of the kinetic energy in both fluids may be estimated by using the gravitational velocity U g , yielding

E k,1 = O 1 16 ρ 1 H 3 U 2 g and E k,2 = O 1 16 ρ 2 H 3 U 2 g .
• the time evolution of the volume integral of the enstrophy in both fluids will also be recorded. This quantity is defined as E r,1 = 1

2 Ω (1 -C)w 2 dV and E r,2 = 1 2 Ω Cw 2 dV , respectively, w = ∇ × u denoting the vorticity.

The simulated values are made dimensionless by using the previous reference estimates of time as well as potential and kinetic energies. For enstrophy, the time maximum of the value obtained with the one fluid model is used, as was proposed in [START_REF] Vincent | A phase inversion benchmark for multiscale multiphase flows[END_REF]. A synthesis of the dimensionless parameters is given in table 1. The results are presented in figure 4. For potential and kinetic energies in fluid 1 and 2, it is observed that the two models are in rather good agreement. They capture the same time evolution with peak of energies being located at the same instants. For intermediate times, when droplet generation is maximum, some discrepancies are noticed. They are due to the different interface tracking techniques and also to single velocity field representation of the two-phase flow by the OFM, compared to twovelocity description brought by the MFA. Concerning enstrophy in fluid 1 and 2, the OFM brings two to three times larger magnitudes than MFA. The peak of enstrophy is located at the same dimensionless time for both models, i.e. t * ≈ 2.5 in fluid 1 and t * ≈ 4 in fluid 2. The major differences observed concerning enstrophy were also reported in [START_REF] Vincent | A phase inversion benchmark for multiscale multiphase flows[END_REF] concerning the comparison of different codes on different grids (up to 2048 3 mesh). It was demonstrated that on a vorticity point of view, a 512 3 grid is not enough to perform a real DNS for the phase inversion benchmark whereas potential and kinetic energies can be assumed converged. In the rest of the work, it will be considered that a 512 3 grid simulation is enough for performing a DNS of the phase separation benchmark and provide a representative simulation of this complex two-phase flow, illustrated in figure 3. If different information is not specified, all the results presented in the following sections are investigated with the OFM and MFA on a 512 3 regular Cartesian mesh.

A priori filtering of turbulent two-phase flows

The Large Eddy Simulation filtering of motion equation for two-phase flows is now explained. Thanks to the phase function C indicating phase k, a lowpass frequency filtering operator G defining the spatial filtering of a given variable Ψ is written as the convolution product of G with Ψ as Ψ = GoΨ. The phase-weighted filtering of the velocity u, i.e. called Favre average, is Parameter Value Unit Table 1: Parameters used to define the dimensionless variables for phase inversion benchmark. also defined as ũ = ρu ρ for the OFM whereas it reads ũk = α k u k α k for the MFA.

t * = t t c = tU g H = t(ρ 2 -ρ 1 ) ρ 1 g 2H t 2.
E * cin,1 = E cin,1 1/16ρ 1 U 2 g H 3 E cin,1 3.4063 - E * cin,2 = E cin,2 1/16ρ 2 U 2 g H 3 E cin,2 3
It is assumed here that the commutation between the spatial filter and the derivative is satisfied, as in single-phase LES modeling [START_REF] Sagaut | Large Eddy Simulation for incompressible flows -An introduction[END_REF][START_REF] Vasilyev | A general class of commmutative filters for les in complex geometries[END_REF]. This assumption is reasonable as the commutation error depends more on the topology of the mesh than on the characteristic of the flow [START_REF] Dakhoul | Improved averaging method for turbulent flow simulation. part 1: theoretical development and application to burger's transort equation[END_REF][START_REF] Dakhoul | Improved averaging method for turbulent flow simulation. part 2: calculation and verification[END_REF]. In our simulation, the mesh is regular and isotropic so that commutation errors should be negligible.

By applying the filtering operators to the one fluid and multi-field models, specific subgrid terms appear in the LES formalism of two-phase flow with resolved scale interfaces. These subgrid terms have previously been discussed in [START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF][START_REF] Liovic | Multi-physics treatment in the vicinity of arbitrarily deformable gas-liquid interfaces[END_REF][START_REF] Labourasse | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF][START_REF] Toutant | Dns of the interaction between a deformable buoyant bubble and a spatially decaying turbulence: a priori tests for les two-phase flow modelling[END_REF][START_REF] Tavares | A priori study for the modeling of LES subgrid scale terms[END_REF] for the OFM and in [START_REF] Fleau | LES modeling with a multifield approach[END_REF] for the MFA. These two-phase LES terms are associated respectively to acceleration, inertia, viscous effects and interface presence. Their presence is due to the filtering and averaging operations of the mass and momentum equations. The common terms for the OFM and MFA are presented in table 2. The formulations are written as they should be estimated by a priori filtering. In fact, for the sake of simplicity, for all ∇Ψ like terms, the commutation between the filtering and derivation operations has been assumed, such that it has been considered that for any variable Ψ, ∇Ψ = ∇Ψ. For OFM, the filtered curvature of the interface κ is equal to -∇ • ∇C ∇C according to the filtering of Brackbill CSF force [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]. An equivalent expression is used for the filtered curvature . For the MFA, as was demonstrated by [START_REF] Fleau | LES modeling with a multifield approach[END_REF], two others subgrid terms exist associated to pressure gradient and drag force. They are not considered in the present work as they do not exist in the OFM. For more details, the reader is referred to [START_REF] Fleau | LES modeling with a multifield approach[END_REF].

Filter • Favre average • One-fluid model Time derivative τ t ρ ∂u ∂t -ρ ∂u ∂t (ρ 2 -ρ 1 )ũτ ia Convection τ c ρu • ∇u -ρ u • ∇u ρ u • ∇u -ρ ũ • ∇ũ Diffusion τ d ∇ • [2µS] -∇ • 2µS ∇ • [2µS] -∇ • 2µ S Surface tension τ st σ κ∇C -κ∇C σ κ∇C -κ∇C Interface advection τ ia u • ∇C -u • ∇C u • ∇C -ũ • ∇C Multi-field model Time derivative τ t ρ k ∂α k u k ∂t - ∂α k u k ∂t - Convection τ c ρ k ∇ • (α k u k ⊗ u k ) ρ k ∇ • (α k u k ⊗ u k ) -∇ • (α k u k ⊗ u k )) -∇ • (α k ũk ⊗ ũk )) Diffusion τ d µ k ∇ • [α k S k ] - -∇ • α k S k Surface tension τ st σ α k κ k ∇α k σ α k κ k ∇α k -α k κk ∇α k -α k κk ∇α k Interface advection τ ia ρ k ∇ • (α k u k ) -∇ • (α k u k ) -
Table 2: Two-phase LES subgrid terms obtained with one-fluid and multifield models by A priori filtering.

It can be observed that 5 subgrid LES terms exist with the standard filter whereas less terms appear with a Favre average, i.e. 4 for the OFM and 2 for the MFA. It has to be kept in mind that the subgrid terms τ c , τ d , τ st and τ ia appear in the mass, momentum and interface tracking equations as extra terms, as soon as a real LES simulation is investigated. For example, if τ ia is not zero, the flow is no more divergence free in the cells cut by the interface, even if no phase change is considered. In the present work, the time derivative subgrid term τ t has not been considered. Indeed, its a priori filtering would require to store at least two successive fields of the 512 3 DNS simulations. This was not achieved due to its numerical cost.

The first interesting issue is the estimate of the relative order of magnitude of the specific two-phase subgrid LES terms which appear in mass, momentum and interface advection equations. An a priori filtering of the DNS results have been provided for τ c , τ d , τ st and τ ia at the enstrophy peak, Filter i.e. for t * ≈ 3. The averaged magnitude of the subgrid terms is normalized by the maximum, at each time step, of the filtered resolved contributions of the flow (inertia for the momentum equation and advection for the phase function equation). The width of the low-pass filtering operator G is 2 cells in each direction. From a global point of view, the magnitude of the vertical component of the subgrid terms is larger than the contributions in the x and y directions. In the z-direction or gravity direction, the flow is strongly anisotropic. This is induced by the key phenomenon of the phase separation, i.e. the buoyancy force. Moreover, it was observed that the magnitude of the subgrid terms is almost the same in the x and y directions. In the rest of this section, only the vertical component of the subgrid terms will be considered. According to table 3, a classification of LES subgrid terms can be deduced:

• contrary to what could be intuitively expected, the inertial or convective subgrid term τ c is not the predominant term. With the OFM, its averaged magnitude is 7.5% the filtered resolved convective term with the top hat filtering and 3% with the Favre average. Concerning the MFA, the averaged estimated values of the subgrid convective term are 2.2% with filter • and 0.9% with the Favre average. This term cannot be neglected in LES modeling of resolved scale two-phase flows. For both OFM and MFA, using a Favre average reduces the order of magnitude of the subgrid contribution by a factor of 2.

• concerning the subgrid viscous contribution τ d , it is clearly the lower one, compared to convective, surface tension and interface advection, with both the OFM and MFA modeling approaches. This conclusion was previously given by [START_REF] Larocque | Parametric study of LES subgrid terms in a turbulent phase separation flow[END_REF][START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF][START_REF] Liovic | Interface-turbulence interactions in large-scale bubbling processes[END_REF][START_REF] Labourasse | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF] for various physical configurations. Its order of magnitude is two-times less than the filtered resolved convective terms. This term could be discarded in future LES modeling of two-phase flows.

• the capillary subgrid contribution τ st is larger than τ c , τ d and τ ia . This conclusion is different from what was obtained for example by [START_REF] Toutant | Dns of the interaction between a deformable buoyant bubble and a spatially decaying turbulence: a priori tests for les two-phase flow modelling[END_REF] for an isolated bubble interacting with a Homogeneous Isotropic Turbulence. This is mainly due to the fact that in our phase separation problem, a kind of atomization process is observed for t * = 3. In our configuration, many droplets are generated. They are poorly resolved on the 512 3 grid with only a few grid cells (3 to 6) over a droplet diameter. This resolution effect is measured by the magnitude of τ st . Indeed, as soon as the poorly resolved interfacial structures are filtered, the corresponding subgrid terms associated to interface, i.e. the surface tension contribution here, becomes significant compared to the filtered convective DNS terms. The same behavior arises for OFM and MFA. As a conclusion, for dispersed flows without scale separation, the subgrid surface tension LES term cannot be neglected. This is the predominant filtered effect. Contrary to what was noticed for the subgrid convective contribution, using a Favre average does not allow to decrease the magnitude of τ st . This term is highly depend on the interfacial structure and the associated fluid properties. The Favre average uses the density in the filtering Subgrid term 128 3 DNS mesh 256 3 DNS mesh 512 3 operations. In the present work, no particular attention has been paid to the local anisotropy of the interface. This way, when a given quantity is filtered in a cell cut by the interface, the density used can be ρ 1 , ρ 2 or an isotropic average of the two. It has been demonstrated by [START_REF] Sagaut | On the filtering paradigm for les of flows with discontinuities[END_REF] that for resolved scale two-phase flows, anisotropic filters should be used in order to build correct LES models in an heterogeneous medium such as the phase separation flow.

• the interface advection subgrid term is more important with OFM than with MFA. The diffuse character of the interface in MFA is certainly responsible for the lower order of magnitude of τ ia with MFA compared to OFM. As for τ st , using a Favre average instead of filter • does not decrease the magnitude of the interface advection subgrid term. In fact, the Favre average reduces the importance of the subgrid terms as soon as the filtering is applied to the velocity field. This effect does not hold when filtering interface depending terms. The subgrid term τ ia is clearly not negligible with OFM while it could be discarded with MFA.

On a general point of view, the magnitude of the two-phase LES subgrid terms is lower with MFA than with OFM. This is a direct consequence of the diffuse character of the two-fluid model. In the OFM, every change of interfacial position is directly correlated to velocity and volume fraction gradients so that the subgrid terms τ c , τ st and τ ia are always significative in the phase inversion problem when the fragmentation mechanisms are observed. In table 4, the width of the low-pass filtering operator G is changed from 2 cells witdh to 4 cells. The effect of this change on the relative magnitude (in %) of two-phase LES subgrid terms is considered with OFM. As was previously observed in a priori filtering works [START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF][START_REF] Labourasse | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF], the larger the compact support of the filter, the higher the magnitude of the subgrid terms is. This is particularly noticeable for τ c and τ ia .

To finish with a priori filtering of DNS results, three simulations are investigated with OFM on 128 3 , 256 3 and 512 3 grids. Even if with the two coarser grids, we are not really able to perform true "DNS", the main goal here is to estimate the effect of implicit LES, i.e. under resolved DNS, on two-phase subgrid contributions. The results are reported in table 5. With top hat filter •, refining the mesh size reduces the magnitude of the convective and interface advection subgrid terms, whereas the diffusive viscous and surface tension contributions increase. These observations are intuitive for τ c and τ ia , not for τ d and τ st . In fact, the filtering effect on velocity decreases with the local mesh size. This is a classical tendency that is observed in single phase flows. This effect is recovered for two-phase flows and resolved scale interfaces as in the present phase separation problem, the flow is isothermal without phase change and so the velocity field is continuous across interface. On a kinematic point of view, the problem is continuous and an isotropic filtering is well adapted for LES filtering. On the contrary, when filtered volume fraction terms are put into play, such as in τ st , the behavior of the subgrid term magnitude according to grid refinement is totally different. Indeed, the more you refine, the more small structures you resolve. As a consequence, the magnitude of the subgrid term is increased as more interfacial structures of larger curvature are obtained. In this case, using a phase-conditioned filter [START_REF] Sagaut | On the filtering paradigm for les of flows with discontinuities[END_REF] should perhaps be more suited to decrease the filtering effect when refining the simulation mesh. Concerning the viscous term, it can be assumed almost constant over all grids. A clear explanation is not straightforward. The same tendencies are observed all subgrid terms with the Favre average. The same conclusions hold. 6: Literature review of LES models that have been investigated for two-phase subgrid LES terms.

LES modeling of turbulent two-phase flows

The understanding and characterization of the coupling between turbulence and fluid interfaces is a difficult task. As we reported in the introduction section, few experimental [START_REF] Lacanette | Macroscopic analysis of gas-jet wiping: Numerical simulation and experimental approach[END_REF][START_REF] Larocque | Macroscopic analysis of a turbulent round liquid jet impinging on an air/water interface in a confined medium[END_REF] and numerical works [START_REF] Larocque | Parametric study of LES subgrid terms in a turbulent phase separation flow[END_REF][START_REF] Liovic | Interface-turbulence interactions in large-scale bubbling processes[END_REF][START_REF] Chesnel | Subgrid modeling of liquid atomization[END_REF][START_REF] Aniszewski | A new approach to sub-grid surface tension for les of two-phase flows[END_REF][START_REF] Calmet | High-schmidt number mass transfer through turbulent gas-liquid interfaces[END_REF][START_REF] Calmet | Statistical structure of high-reynoldsnumber turbulence close to the free surface of an open-channel flow[END_REF][START_REF] Calmet | Large-eddy simulation of high-Schmidtnumber mass transfer in a turbulent channel flow[END_REF] exists to investigate the formulation of LES models for two-phase flows involving resolved scale interfaces. A synthesis of literature works is proposed in table 6. Among the few existing publications, it can be observed that none of the authors have studied all the subgrid terms for a same two-phase problem. In addition, none of the authors have considered τ d , that is always assumed negligible. The goal of the present section is to contribute, for a given phase inversion benchmark problem, to the evaluation of LES models for all two-phase subgrid terms that have been formulated in section 4, except the time subgrid term τ t .

Based on single-phase LES [START_REF] Sagaut | Large Eddy Simulation for incompressible flows -An introduction[END_REF], functional models are first considered for two-phase LES subgrid terms. By nature, these functional models have been designed for inertial (or convective) terms. They rely on an energy cascade assumption and a dissipation modeling of turbulence at small scale. In this way, they can only be applied to subgrid terms τ c and τ ia as follows:

• the Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations. i : The basics experiment[END_REF] is first considered for its simplicity and robust behavior. For the convective term, we classically obtain

τ c = -2∇ • (C s ∆) 2 S f S (13)
By analogy, the interfacial subgrid term is written as

τ ia = -2∇ • (C s ∆) 2 S ∇C (14)
With this approach, a subgrid scale LES viscosity can be defined as ν s = (C s ∆) 2 S . The Smagorinsky constant C s is generally chosen between 0.1 and 0.24.

• 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] which is a modification of the Smagorinsky model in order to better predict turbulence behavior near solid boundaries. In particular, this model recovers the expected asymptotic values of the subgrid viscosity near the wall. It reads

τ c = -2∇ • ν w f S (15)
τ ia = -2∇ • ν w ∇C (16)
with

ν w = (C w ∆) 2 Z : Z 3/2
S : S 5/2 + Z : Z 5/4 and C w in the range 0.55 -0.6. In addition,

Z = SS + ΩΩ - 1 3 S : S -Ω : Ω Id (17)
with Ω the rotation tensor defined as

Ω = 1 2 ∇u -∇ t u (18)
In all LES models, f is 1 for the OFM whereas f = α k for MFA. In addition, for this last model, a k index have to be added in all LES expressions for velocity depending terms. To finish with, terms based on C have to be switched into terms depending on α k .

A second class of LES approaches can be used to design models for twophase subgrid terms. It is called structural analysis [START_REF] Sagaut | Large Eddy Simulation for incompressible flows -An introduction[END_REF]. The main interest of the structural LES is that it does not rely on an explicit understanding and formulation of the coupling between resolved and filtered scales of turbulence (and interfaces in the case of two-phase flows). The objectives and strategy of these approaches are to approximate the inverse of the LES filter and to build models with same mathematical structure as the subgrid term. A very interesting aspect with these formulations is that they use only the resolved filtered variables. Once G -1 and the filtered variables are known, the real unknowns can be reconstructed, as well as the subgrid contributions. Two main structural approaches have been selected:

• The Bardina model [START_REF] Bardina | Improved subgrid scale models for large eddy simulation[END_REF] considers a scale similarity which is to assume that the statistical structure of subgrid terms is similar to the terms evaluated from the smaller resolved scales. With this assumption, the filtering of a product is simply the product of filtered variables, i.e. ab = ab. After some manipulations, the two-phase subgrid LES terms read for the OFM:

τ c = C b ρu • ∇u -ρ u • ∇u (19)
τ d = C b ∇ • 2µS -∇ • 2µS (20)
τ st = C b σ κ∇C -κ∇C (21)
τ ia = C b u • ∇C -u • ∇C (22)
with κ = -∇ • ∇C ∇C . For the MFA, it reads:

τ c = C b ρ k ∇ • (α k u k ⊗ u k ) -∇ • α k u k ⊗ u k (23)
τ d = C b µ k ∇ • α k S k -∇ • α k S k (24)
τ st = C b σ α k κk ∇α k -α k κ k ∇α k (25)
τ ia = C b ρ k ∇ • (α k u k ) -∇ • α k u k (26)
with

κ k = -∇ • ∇α k ∇α k . The Bardina constant C b is adjustable. A value
of 1 has been used here.

• The Approximate Deconvolution Model (ADM) [START_REF] Adams | A subgrid-scale deconvolution approach for shock capturing[END_REF] is a general framework for building models that respect the mathematical structure of the subgrid terms. First, an approximation of the inverse G -1 is build with a given order of accuracy, in a kind of Taylor expansion. Based on G -1 , reconstructed unfiltered variables are build and used directly in the native structure of the subgrid LES terms as follows for the OFM:

τ c = C a ρu * • ∇u * -ρ u * • ∇u * (27)
τ d = C a ∇ • [2µS * ] -∇ • 2µS * (28)
τ st = C a σ κ * ∇C * -κ * ∇C * (29)
τ ia = C a u * • ∇C * -u * • ∇C * (30)
And for the MFA:

τ c = C a ρ k ∇ • (α k * u k * ⊗ u k * ) -∇ • (α k * u k * ⊗ u k * ) (31)
τ d = C a µ k ∇ • [α k * S * k ] -∇ • α k * S k * (32)
τ st = C a σ α k * κ k * ∇α k * -α * k κk * ∇α k * (33)
τ ia = C a ρ k ∇ • (α k * u k * ) -∇ • (α k * u k * ) (34)
with

u * ≈ G -1 7 ou, C * ≈ G -1 7 oC and α k * ≈ G -1 7 oα k .
The approximate inverse filter of order k is such that

G -1 k = l=0,k (1 -G) l (35)
γ * = G -1 oγ ≈ l=0,k (1 -G) l oγ ≈ γ + γ -γ + γ -γ + γ + ... (36)
We have used a 7 th order of approximation for G -1 , instead of a classical 5 th order found in literature [START_REF] Adams | A subgrid-scale deconvolution approach for shock capturing[END_REF][START_REF] Aniszewski | A new approach to sub-grid surface tension for les of two-phase flows[END_REF], because we have observed better accuracy of modeled subgrid terms compared to filtered DNS terms. For higher orders, no significative improvements have been found. These results are presented in the next section. It has to be stressed that the Bardina model is no more than the ADM model of 2 nd order.

The last type of model that has been considered is the mixed Smagorinsky-Bardina model [START_REF] Bardina | Improved turbulence models based on large eddy simulation of homogeneous, incompressible, turbulent flows[END_REF]. It combines the positivity property of the functional model with the structural representativeness of the Bardina model. Due to the intrinsic construction of the functional part of the model, it can only be applied to convective and interfacial subgrid terms. It reads for the OFM:

τ c = 1 2 f ∇ • -2(C s ∆) 2 S S + C b ρu • ∇u -ρ u • ∇u (37)
τ ia = 1 2 f ∇ • -2(C s ∆) 2 S ∇C + C b u • ∇C -u • ∇C (38)
And for the MFA: Thanks to DNS simulations, the two-phase LES subgrid terms can be directly calculated by applying filter G to the various variables of the phase separation problem. Moreover, the various LES models applied to τ c , τ d , τ st and τ ia can also be calculated with the DNS results. The main objective of the present section is to evaluate the representativeness of LES models presented in the previous section compared to filtered DNS subgrid terms. Results are presented in figure 5 in terms of correlation between a given model (Smagorinsky, Wale, Bardina, Mixed, ADM), for a given subgrid LES term, i.e. τ c , τ d , τ st and τ c , and the filtered DNS of a corresponding subgrid term. The evaluation of LES models is provided for both OFM and MFA. On a general point of view, it is observed that compared to OFM, the models Figure 5: Spatial correlation of two-phase LES subgrid terms -The modeled terms are plotted according to the filtered DNS subgrid terms. A perfect modeling would be on the y=x line. Five LES models are considered: Smagorinsky [START_REF] Smagorinsky | General circulation experiments with the primitive equations. i : The basics experiment[END_REF], WALE [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF], Bardina [START_REF] Bardina | Improved subgrid scale models for large eddy simulation[END_REF], Mixed model [START_REF] Bardina | Improved turbulence models based on large eddy simulation of homogeneous, incompressible, turbulent flows[END_REF] and ADM [START_REF] Adams | A subgrid-scale deconvolution approach for shock capturing[END_REF]. Left column is for the one-fluid model whereas right column is for the multi-field approach.

τ c = 1 2 ∇ • -2(C s ∆) 2 S k α k S k + C b ρ k ∇ • (α k u k ⊗ u k ) -∇ • α k u k ⊗ u k (39)
τ ia = 1 2 ∇ • -2(C s ∆) 2 S k ∇α k + C b ρ k ∇ • (α k u k ) -∇ • α k u k (40
Figure 6: Comparison of equivalent turbulent viscosities obtained with DNS and LES models, i.e. Smagorinsky [START_REF] Smagorinsky | General circulation experiments with the primitive equations. i : The basics experiment[END_REF], WALE [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF], Bardina [START_REF] Bardina | Improved subgrid scale models for large eddy simulation[END_REF], Mixed model [START_REF] Bardina | Improved turbulence models based on large eddy simulation of homogeneous, incompressible, turbulent flows[END_REF] and ADM [START_REF] Adams | A subgrid-scale deconvolution approach for shock capturing[END_REF]. Left column is for OFM and right column for MFA. The bottom line is a zoom on the middle line values. The results are plotted in a vertical line centered on [x = 0.6, y = 0.6] for OFM and [x = 0.8, y = 0.8] for MFA. The volume fraction is also plotted to locate interface position as C = 0.5. applied to MFA involve less dispersion of the correlation between LES models and filtered DNS. Indeed, the phase average α k is present in all terms for MFA. It tends to eliminate the LES modeling far from the interfacial zone and it damps the gradients or filtering operations in the interfacial zone. On the contrary, in the OFM, all values of the models are a priori sharper and present everywhere, including zones far from the interface. In this case, a larger dispersion of correlation between LES models and filtered DNS subgrid terms is induced. If attention is now paid to the quality of models for a given subgrid term, several observations arise:

• only structural models are able to model all subgrid terms. Indeed, functional based models such as Smagorinsky, Wale or Mixed do rely on a physical interpretation of the behavior and interaction of under-resolved scales with large scales of turbulence (energy cascade of Kolmogorov, dissipative effect of small scale turbulent structures, ...). They are mainly based on interpretation of kinetic energy and inertial effects in turbulent flows. However, they do not consider the multi-scale character of interfaces when two-phase flows are dealt with. They so can only be applied to advection or inertial terms in conservation equations. This is the reason why functional LES approaches are not possible to be investigated for τ d and τ st . If we pretend to perform real LES simulations accounting for all two-phase subgrid terms, using functional models seems to be inappropriate as they do not apply for all subgrid terms. This discussion is valid for both OFM and MFA.

• for the inertial and interfacial terms, all LES models can be compared. It is observed that for both one-fluid and multi-field approaches, the Smagorinsky and Wale models are not representative of filtered DNS terms. The Mixed model is better as it incorporates a contribution from the Bardina model that is working pretty well. The ADM approach is clearly the best for all terms and all modeling approaches (OFM and MFA). The functional based models involve a lot of dispersion and they do not respect the y = x slope that should be obtained if a perfect correlation between LES models and filtered DNS will be reached. The interesting point here is that ADM based models (and also Bardina but with poorer results) are working for all subgrid terms and all two-phase flow modeling approaches (OFM and MFA). The average percentage error that is obtained for all models, all subgrid terms and all two-phase flow models is given in table 7. It can be clearly inferred that Smagorinsky and Wale model are not suitable for two-phase subgrid terms, even concerning inertia. The representativeness of Mixed model is better than Smagorinsky, as it incorporates a part of Bardina modeling. The Bardina model is performing quite nicely for all subgrid terms. It behaves less well for τ d and τ ia when OFM is considered whereas its quality is worse for τ d and τ st with regard to MFA. It could be imagined to play on the Bardina constant C b so as to improve the performance of this model and to get closer to y = x curve on the correlation plots. Finally, the ADM model is the best for all subgrid terms and all two-phase flow modeling approaches. The maximum error with MFA is at least 20% whereas this error is in the range 17 to 32% for OFM, with a peak of error of 58% for τ ia .

If true LES were simulated, the use of LES models for convective subgrid terms will be equivalent to add a turbulent viscosity in the viscous term of the momentum equations. It is interesting to evaluate the equivalent turbulent viscosity resulting from the use of a given LES model, compared to the equivalent viscosity that the filtering of the DNS provides when the two-phase convective LES subgrid term is considered. Classically [START_REF] Sagaut | Large Eddy Simulation for incompressible flows -An introduction[END_REF][START_REF] Labourasse | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF], we can write the turbulent viscosity µ t for the OFM as follows:

µ t = τ c : ∇u S : ∇u (41)
And for the MFA:

µ t = α 1 τ c,1 : ∇u 1 S 1 : ∇u 1 + α 2 τ c,2 : ∇u 2 S 2 : ∇u 2 (42)
Equivalent LES turbulent viscosities compared to equivalent filtered DNS one are plotted in figure 6. For OFM, the instantaneous turbulent viscosities are extracted in a vertical slice starting from point (0.6, 0.6, 0) whereas for MFA, the strating of the slice is (0.8, 0.8, 0). A time t * ≈ 3 is considered for the peak of enstrophy. The instantaneous simulations of OFM and MFA at a given time are not exhibiting the same interfacial and flow structures. The slices for turbulent viscosities have been chosen in order to cut different times the interface. If we focus on the equivalent DNS viscosity that we should obtained with LES models, we observe that negative values are obtained. This have previously been been reported by [START_REF] Labourasse | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF][START_REF] Toutant | Dns of the interaction between a deformable buoyant bubble and a spatially decaying turbulence: a priori tests for les two-phase flow modelling[END_REF]. This important feature of the effect of inertial subgrid term cannot be represented by functional LES approaches as they are intrinsically built on positive terms. This is illustrated in the line of figure 6. This way, Smagorinsky and Wale models are again discriminated against structural ones. If we now compare Bardina, ADM and Mixed models, we notice that the first two are much better as expected. Concerning equivalent turbulent viscosities, we observe that ADM is globally nearer to DNS than Bardina, even if locally, Bardina can exhibit closer values to DNS reference.

To conclude on LES modeling of two-phase subgrid LES terms, it has been demonstrated that the structural approaches are the correct way of building models. The ADM formalism allows to recover the negative values of equivalent turbulent viscosities in the good zones of the flow and they respect the behavior of filtered DNS subgrid terms with less than 30% discrepancy (50% for τ ia ). We have to keep in mind that ADM modeling allows to reconstruct the fields between the explicit filter that we have and the implicit filter brought by the simulation mesh. In real LES, an additive relaxation term will have to be incorporated in the momentum equations [START_REF] Adams | A subgrid-scale deconvolution approach for shock capturing[END_REF] in order to obtain the effect of under-resolved scales smaller than the implicit filter size.

Summary and future work

Based on direct numerical simulation of two-phase flow with resolved scale interfaces, i.e. a phase inversion problem, a priori filtering of two-phase subgrid LES terms has been investigated. The one-fluid model and multi-field approaches were considered and filtered. In addition, different LES models taken from the turbulent single-phase flow literature were utilized to propose functional and structural LES models for two-phase LES subgrid terms. The order of magnitude of all a priori filtered subgrid terms and LES models for these terms were compared to filtered DNS values for the first time on a same two-phase problem. The formulation of LES models for all subgrid terms was also proposed and discussed for the first time, in particular for a multi-field approach.

Concerning a priori filtering, the main conclusion is that the inertia term is not predominant in two-phase flows with fragmentation and rupture of interface that exhibit a kind of atomization process. This conclusion is different from that of the studies of [START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF][START_REF] Liovic | Interface-turbulence interactions in large-scale bubbling processes[END_REF][START_REF] Chesnel | Subgrid modeling of liquid atomization[END_REF][START_REF] Labourasse | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF][START_REF] Toutant | Dns of the interaction between a deformable buoyant bubble and a spatially decaying turbulence: a priori tests for les two-phase flow modelling[END_REF]. As soon as small-scale interfacial structures are generated, correlation terms induced by the curvature in the subgrid surface tension terms or between the local interface presence and the velocity are predominant against inertial terms. The diffusive viscous term is always small in our phase inversion problem, whatever the two-phase flow model or the kind of filtering.

Five different LES models applied to two-phase subgrid LES terms have been compared against filtered DNS results. Two classes of models have been considered, the functional (Smagorinsky, Wale) and the structural (Bardina, ADM) approaches. A combined functional-structural model, called Mixed Bardina-Smagorinsky, have also been investigated. Conclusion are clear for both OFM and MFA: functional LES modeling do not correlate to filtered DNS results whereas structural approaches do. Bardina and ADM are clearly the good LES framework to consider for two-phase flows with resolved scale interfaces. ADM is clearly better than Bardina in our study. It has to be noticed that we do not have trigered the Bardina constant to try to improve the Bardina results. This could be usefull to try to get closer to ADM. The improvement of Bardina, if possible, would be a nice point as ADM of 7 th is clearly more CPU time consuming. The Mixed approach, as expected, is bet-ter than functional models and worse than structural LES techniques. For the first time, we have been able to demonstrate that a common LES framework could be used for modeling all two-phase subgrid terms with resolved scale interfaces, i.e. the approximate deconvolution method. The main interest of structural approaches is their capability of representing negative equivalent turbulent viscosities whereas functional models are generating intrinsically positive viscosities. We have shown that filtered equivalent DNS viscosities are clearly negative in specific zones of the flow, not only near the interface. This negative feature is expected from LES models. Future works will be oriented in two ways:

• performing a priori LES filtering and model comparisons on finer DNS meshes. A 2048 3 grid is under consideration for phase inversion. The same analysis will be used on other two-phase flows with separated phases, such as the atomization of a liquid sheet in a HIT flow [START_REF] Trontin | Detailed comparisons of front-capturing methods for turbulent two-phase flow simulations[END_REF]. The idea is to control the turbulence forcing and to vary the properties of the fluid and interface properties to try to extract general conclusions on LES modeling frameworks for two-phase flows.

• investigating real ADM LES with both OFM and MFA. A coarse LES grid will be used in this case. The idea is to see if compared to DNS, the LES are able to capture the main features of the flow. The other important point is to compare implicit LES, i.e. LES without explicit modeling of the subgrid terms, and ADM models for all two-phase subgrid terms. In particular, the role of the introduction of a relaxation term [START_REF] Adams | A subgrid-scale deconvolution approach for shock capturing[END_REF] in the momentum equations will be discussed and evaluated compared to subgrid LES contributions. The starting test problem could be the phase inversion benchmark, for which different Re and W e can be chosen in order to adjust inertial and capillary effects compared to viscous and gravity contributions. The interest of this case is that DNS can be generated and used for LES comparisons.

Figure 2 :

 2 Figure 2: Definition sketch for the phase inversion benchmark in a closed box [2, 3, 11].

Figure 3 :

 3 Figure 3: Simulation of a phase inversion in a closed box with a one-fluid model -The interface between the two liquids is plotted in orange -from left to right and top to bottom, the dimensionless time is 0, 0.246, 0.492, 0.738, 0.984, 1.329, 1.969, 3.396 and 15.305.

  .7847 -Maximum of enstrophy in fluid 1 (one-fluid model) 56.72 m 3 .s -2 Maximum of enstrophy in fluid 2 (one-fluid model) 224.03 m 3 .s -2

Figure 4 :

 4 Figure 4: Time evolution of macroscopic quantities for one-fluid and multifield models.

Table 3 :

 3 Relative magnitude (in %) of the vertical component of two-phase LES subgrid terms obtained with one-fluid and multi-field models after a priori filtering. The terms are made dimensionless by the magnitude of the vertical component of the largest resolved DNS term, i.e. the convective contribution. The width of the low-pass filtering operator G is 2 cells in each direction.

			• Favre average	•
	One-fluid model
	τ c	7.54	3.01
	τ d	1.19	0.14
	τ st	37.81	37.81
	τ ia	10.29	9.61
	Multi-field model
	τ c	2.24	0.86
	τ d	0.14	-
	τ st	5.29	5.29
	τ ia	1.07	-

Table 4 :

 4 Effect of width of the low-pass filtering operator G on the relative magnitude (in %) of the vertical component of two-phase LES subgrid terms obtained with one-fluid model after a priori filtering. The terms are made dimensionless by the magnitude of the vertical component of the largest resolved DNS term, i.e. the convective contribution.

	Subgrid term Filter • with 2 cells width Filter • with 4 cells width
	τ c	7.54	13.24
	τ d	1.19	1.96
	τ st	37.81	44.18
	τ ia	10.29	16.17

Table 5 :

 5 Effect of the DNS mesh size on the relative magnitude (in %) of the vertical component of two-phase LES subgrid terms obtained with onefluid model after a priori filtering. The terms are made dimensionless by the magnitude of the vertical component of the largest resolved DNS term, i.e.

	DNS mesh

the convective contribution. The width of the low-pass filtering operator G is 2 cells in each direction.

Table

	5.1 Large Eddy Simulation (LES) models for two-phase
	flows involving resolved scale interfaces
	Physical case	τ t	τ c	τ d	τ st	τ ia
	THI interacting	-	Smagorinsky [5]	-	-	-
	with a 2D droplet		Wale [6], Bardina [7]
	[31]		Mixed [8]
	THI interacting Bardina	Bardina	-	-	Bardina
	with a 3D droplet	[7]	[7]			[7]
	[32]
	Liquid jet	-	Smagorinsky	-	-	Smagorinsky
	atomization		Lilly			Lilly
	[30]		[5]			[5]
	Phase inversion	-	-	-ADM	-
	in a 3D cavity				[9]
	[32]

Table 7 :

 7 Average relative error (in %) of the modeled two-phase LES subgrid terms compared to the filtered DNS subgrid terms. Both one-fluid and multifield models are considered. Filter • is applied to DNS, with a width of the low-pass filtering operator G equal to 2 cells in each direction.

	)

Acknowledgements

The authors are grateful for access to the computational facilities of the French CINES (National computing center for higher education) and CCRT (National computing center of CEA) under project number x20152b6115 and x20162b6115. The Agence Nationale de la Recherche (ANR) through the MODEMI project is associated to this work. Thank EDF for the financial support it has given to this work.