Davy Nayigizente

Thomas Schmitt
email: thomas.schmitt@ecp.fr

Sébastien Ducruix

Development of an interface thickening method for the direct numerical simulation of compressible liquid-vapor flows in the framework of the Second Gradient theory

Keywords: Second Gradient, Two-phase flow, Real gas, Direct Numerical Simulation

Development of an interface thickening method for the direct numerical simulation of compressible liquid-vapor flows in the framework of the Second Gradient theory

Davy Nayigizente, Schmitt Thomas, Sébastien Ducruix

The study of so called real gas flows in combustion chambers is still an active field of research, see [START_REF] Bellan | Supercritical (and subcritical) fluid behavior and modeling: drops, streams, shear and mixing layers, jets and sprays[END_REF], [START_REF] Dahms | Liquid jet breakup regimes at supercritical pressures[END_REF], [START_REF] Oefelein | Modeling high-pressure mixing and combustion processes in liquid rocket engines[END_REF], [START_REF] Okong'o | Direct numerical simulation of a transitional supercritical binary mixing layer: heptane and nitrogen[END_REF], , [START_REF] Schmitt | Large-eddy simulation of supercritical-pressure round jets[END_REF], [START_REF] Yang | Modeling of supercritical vaporization, mixing, and combustion processes in liquid-fueled propulsion systems[END_REF].

However, during the early moments of the ignition, fuel is injected at subcritical pressures often leading to the presence of a liquid phase, which necessitates an additional treatment to deal with two-phase flow effects such as the presence of interfaces, capillary phenomena and phase change. Two-phase flow models are usually separated in two classes with a distinct handling of the interface geometry and thermodynamics.

. Sharp Interface Methods (SIM) consider the interface as a discontinuity that must be accurately located during the simulation. They are comprised of two main approaches. Interface tracking methods such as Arbitrary Lagrangian-Eulerian in [START_REF] Chan | A generalized arbitrary lagrangian-eulerian method for incompressible flows with sharp interfaces[END_REF], [START_REF] Hirt | An arbitrary lagrangian-eulerian computing method for all flow speeds[END_REF], Boundary Integrals from [START_REF] Birkhoff | Helmholtz and taylor instability[END_REF][START_REF] Rocard | The formation of vortices from a surface of discontinuity[END_REF], Front-Tracking by [START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF], [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF] or Marker-And-Cell of [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF], [START_REF] Welch | The MAC method -A computing technique for solving vis-630 cous, incopressible, transient fluid-flow problems involving free surfaces[END_REF] rely on a specific mesh deforming with the flow or Lagrangian vortices/markers advected with the flow to locate its different constituents. These methods are usually extremely accurate for small and regular movements of the interface but struggle to deal with large scale and complex topological changes, in particular interface break-up and reconnection. In modern literature, interface capturing methods are often favored over the former, the most commonly used being the Level-Set [START_REF] Osher | Level set methods: an overview and some recent results[END_REF], [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF] and the Volume-of-Fluid (VOF) [START_REF] Hirt | An arbitrary lagrangian-eulerian computing method for all flow speeds[END_REF], [START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF] methods. Level-Set methods represent the interface as an iso-surface of the so-called level-set function, usually a signed distance function. Early Level-Set methods offered a simple yet precise representation of the interface topology but lacked conservative properties.

Since then, several improvements have been made to achieve conservative Level-Set methods such as in [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF] or [START_REF] Owkes | A discontinuous galerkin conservative level set scheme for interface capturing in multiphase flows[END_REF]. VOF methods reconstruct the interface geometry using a phase volume fraction advected on the computational mesh. They tend to offer better conservative properties than Level-Set methods but the interface reconstruction procedure displays an important complexity that grows rapidly with the dimension of the problem, the topology changes of the interface and the use of non-cartesian meshes. The coupling of VOF and Levet-Set is a promising way to take full advantage of both methods strengths while accommodating their drawbacks, see [START_REF] Balcázar | A coupled volume-of-fluid/level-set method for simulation of two-phase flows on unstructured meshes[END_REF], [START_REF] Sun | A coupled volume-of-fluid and level set (voset) method for computing incompressible two-phase flows[END_REF], [START_REF] Sussman | A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows[END_REF], [START_REF] Wang | A coupled level set and volume-of-fluid method for sharp interface simulation of plunging breaking waves[END_REF]. For both Level-Set and VOF methods, adjusted jump conditions allow an accurate treatment of heat and mass fluxes through the interface to address phase change while surface tension is usually incorporated using the Convected-Distributed Force (CDF) model of [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF] or the Continuous Surface Force (CSF) model of [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]. Both CDF and CSF models require an accurate location of the interface thus the importance allocated to this task in the sharp methods. Despite the mentioned hurdles, both Level-Set and VOF methods have allowed groundbreaking numerical simulations such as in [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF] or [START_REF] Ménard | Coupling level set/vof/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet[END_REF] and are now established as a standard to deal with incompressible flows. Moreover, their application to compressible flows has sustained substantial improvements such as the ones proposed in [START_REF] Hu | A conservative interface method for compressible flows[END_REF], [START_REF] Terashima | A front-tracking/ghost-fluid method for fluid interfaces in compressible flows[END_REF] which foresees sufficient robustness and efficiency.

. In the other class, the Diffuse Interface Methods (DIM), the interface is considered as a volumetric region within which the thermodynamic variables vary rapidly yet continuously. By principle, DIM do not require to locate the interface. The smoothing of the interface implies the creation of an interfacial region where classical thermodynamic treatments are no longer valid. The two approaches present among the DIM offer specific ways to treat the thermodynamics in the interface. They are particularly adapted to the study of compressible flows, as treated in this work, hence our focus on this class of methods. The Multi-Fluid methods (MFM), see [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF], [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (ddt) in reactive granular materials[END_REF], [START_REF] Kapila | Two-phase modeling of deflagration-todetonation transition in granular materials: Reduced equations[END_REF], [START_REF] Powers | Theory of two-phase detonation-part ii: structure[END_REF], [START_REF] Saurel | A finite volume scheme for two-phase compressible flows[END_REF], are based on a definition of an artificial mixture between the two phases in the interfacial region. Depending on the problem, additional equilibrium hypothesis can be enforced upon the two phases in the artificial mixture, for instance velocity, temperature, pressure and/or chemical equilibrium. From this Eulerian description ensues a set of four to seven equations to solve at once, usually involving a volume or a mass fraction. MFM have been introduced for shock capturing purposes and have been extended to other types of flows by [START_REF] Gaillard | Interfaces diffuses et flammes transcritiques LOX/H2[END_REF], [START_REF] Murrone | Modèles bi-fluides à six et sept équations pour les écoulements diphasiques à faible nombre de Mach[END_REF], [START_REF] Perigaud | A compressible flow model with capillary effects[END_REF]. In particular, they articulate efficiently with the CSF model and allow the use of separate and more efficient EoS to describe the different phases. The second approach is formed by the Phase-Field methods (PFM) built after the work in [START_REF] Allen | A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening[END_REF], [START_REF] Cahn | Free energy of a nonuniform system. i. interfacial free energy[END_REF], van der Waals (1893), among others. In these methods, a phase parameter c is used to discriminate both phases, in a fashion similar to that of the Level-Set approach, the values of c varying continuously in the interface. However, the phase parameter definition is usually thermodynamically motivated: a component concentration for non-miscible fluids, the density for a two-phase compressible flows, etc... As such, the convection-diffusion equation of the phase parameter differs from the purely advective equation of the level-set function. Additionally, the phase parameter c is explicitly used to express the capillary forces, following the work of [START_REF] Korteweg | Sur la forme que prennent les équations du mouvement des fluides si l'on tient compte des forces capillaires causées par des variations de densité considérables mais continues et sur la théorie de la capillarité dans l'hypothese d'une variation continue de la densité[END_REF]. Interestingly, PFM combine quite easily with real gas thermodynamics as shown in [START_REF] Kou | An adaptive finite element method for simulating surface tension with the gradient theory of fluid interfaces[END_REF] where the Peng-Robinson EoS has been used.

Other applications can be found in [START_REF] Lee | An efficient numerical method for simulating multiphase flows using a diffuse interface model[END_REF], [START_REF] Li | A conservative numerical method for the cahn-hilliard equation with dirichlet boundary conditions in complex domains[END_REF] or [START_REF] Li | A phase-field fluid modeling and computation with interfacial profile correction term[END_REF]. A key element for DIM is controlling the interface width. Given that for most of these methods, the physical process driving the interface thickness are unknown or lack control, two points of view confront. Traditionally in MFM, the interface is maintained as sharp as possible using appropriate numerical treatments, as these methods were initially designed for shock capturing. Conversely, in PFM, the purpose is rather to have an interface wide enough to be resolved on the computational mesh because the capillary forces, expressed from the phase parameter gradients, need the latter to be properly evaluated. In that respect, the interface diffusion intrinsic to the DIM can act both as a help and an inconvenience.

. The Second Gradient (SG) theory, introduced by van der Waals (1893) qualifies as a Phase-Field method where the density ρ of the flow is chosen as the phase parameter, adding thermodynamic consistency to its derivation.

This consistency shows through the values of the surface tension and the interface width that are physically predicted by the model, a feature that distinguishes it from other DIM, in particular PFM. Besides, the model can accommodate both subcritical and supercritical regimes indiscriminately. These particularities, along with the other distinctive features of the PF methods motivated the use of the SG theory for our study. Nonetheless, a major difficulty remains as the typical width of an interface, as predicted by the model, varies between a few Å and a few µm as soon as the temperature moves away from its critical value T c . If one is to perform a straightforward numerical simulation of such systems, one must ensure that the mesh allows a decent mapping of the interface, of which width is not known a priori. This problem is shown prohibitive in terms of computational cost, even in simple 2D academic configurations. In the framework of a cryogenic rocket-engine, the system mostly works at high temperatures and high pressures once initiated. However, in order to simulate the complete ignition of the engine, one would encounter reduced temperatures with significantly lower values of interface width, hence the attention given to this issue in the present work. As already shown in [START_REF] Fouillet | Généralisation à des mélanges binaires de la méthode du Second Gradient et application à la simulation numérique directe de l'ébulition nucléée[END_REF], [START_REF] Gaillard | Interfaces diffuses et flammes transcritiques LOX/H2[END_REF] or [START_REF] Jamet | Etude des potentialités de la théorie du Second Gradient pour la simulation numérique directe des écoulements liquide-vapeur avec changement de phase[END_REF] the model provides a framework to achieve the thickening of an interface, with thermodynamic consistency, hopefully allowing the latter to be captured on typical DNS meshes from the turbulence point of view. In this work, the SG theory is further investigated using a new thickening strategy based on a consistent thermodynamic modification of the interface.

. The Second Gradient theory is recalled in Section 2 where its predictive capabilities regarding the saturation values (liquid and vapor densities, vapor pressure) as well as the unsteady response to mechanical perturbations applied to the interface are demonstrated. In Section 3, an in-depth analysis of the macroscopic variables of interest is proposed, namely the surface tension and the interface width, since these are the variables one primarily wants to control when applying a thickening method to the interface. Our new thickening method is presented in Section 4 after a brief overview of the existing strategies. It is then validated with steady and unsteady simulations in 1D and 2D configurations in Section 5. Finally we discuss the conclusions and the prospects of our study in Section 6.

1. Thermodynamics of a capillary fluid 1.1. Equation of state . At very high pressures or in the liquid phase, the assumptions used to derive the ideal gas law, mostly that the components of the fluid merely interact which each other except through inelastic collisions, do not hold anymore.

Since van der Waals (1894) several equations of state (EoS) have been proposed, relying on the molecular and kinetic theories of fluid, to achieve a better description of the real gas thermodynamics. With respect to these EoS, the state variables not only depend on the temperature but also on a secondary intensive variable such as the pressure or the density. The most commonly used EoS are the cubic ones as proposed in [START_REF] Peng | A new two-constant equation of state[END_REF], van der Waals (1894), or the Soave-Redlich-Kwong (SRK) EoS in [START_REF] Soave | Equilibrium constants from a modified redlich-kwong equation of state[END_REF] of which expression is given by: . However, for temperatures lower than T c , despite the ability of the cubic EoS to describe both liquid and vapor phases, it allows the existence of the spinodal region where the thermodynamics becomes unstable, see [START_REF] Giovangigli | Multicomponent flow modeling[END_REF]. In practice, a state in this region corresponds to a two-phase regime where both the liquid and the vapor phases coexist, a distinctive feature for an interface. DIM offer multiple specific ways to describe the thermodynamics in the unstable region.

P EoS (ρ, T ) = ρrT 1 -bρ - a (T ) ρ 2 1 + bρ (1

Second Gradient theory

Van der Waals demonstrated in van der Waals (1893) that the apparent thermodynamic instability in the spinodal region can be cleared up by improving the internal description of the fluid. The pressure curve in this region bears an actual physical meaning as a partial description of the internal constraints in the fluid, that must be completed. His strategy has been adapted afterwards to binary fluids in [START_REF] Cahn | Free energy of a nonuniform system. ii. thermodynamic basis[END_REF], [START_REF] Cahn | Free energy of a nonuniform system. i. interfacial free energy[END_REF] and [START_REF] Cahn | Free energy of a nonuniform system. iii. nucleation in a two-component incompressible fluid[END_REF]. It was later used for DNS in [START_REF] Fouillet | Généralisation à des mélanges binaires de la méthode du Second Gradient et application à la simulation numérique directe de l'ébulition nucléée[END_REF], [START_REF] Jamet | Etude des potentialités de la théorie du Second Gradient pour la simulation numérique directe des écoulements liquide-vapeur avec changement de phase[END_REF] or [START_REF] Seppecher | Etude d'une modelisation des zones capillaires fluides: interfaces et lignes de contact[END_REF] and is recalled in the next paragraphs.

1.2.1. Thermodynamic description . In the Second Gradient theory, the thermodynamics of the interface depends on the temperature T , the density ρ and additionally the density gradient ∇ρ. The volumetric free energy F reads:

F (T, ρ, ∇ρ) = F EoS (T, ρ) + λ 2 (∇ρ) 2 (2)
where F EoS is the volumetric free energy as directly given by the EoS and λ is called the capillary coefficient.

A more precise thermodynamic description, accounting for both density and entropy gradients can be found in 125 [START_REF] Casal | Invariance properties of inviscid fluids of grade n[END_REF] but will not be discussed in the present work.

From Eq. ( 2) one can write the modified Gibbs relation given by Eq. ( 3) where s is the specific entropy and derive the expression of the other variables like the thermodynamic pressure P th in Eq. ( 4), the specific internal energy e s in Eq. ( 5) and the specific free energy f in Eq. ( 6). Noticeably the expressions of the specific enthalpy h, the chemical potential µ and the specific isochoric heat capacity C v are not modified by this new thermodynamic description under the assumption, later justified, that λ is constant.

de s = T ds + P ρ 2 dρ + λ∇ρ • d∇ρ P th = P EoS -λ (∇ρ) 2 e s = e s EoS + λ 2ρ (∇ρ) 2 f = f EoS + λ 2ρ (∇ρ) 2 (3) (4) (5) (6)
1.2.2. Mechanical equations . To remain consistent, this new thermodynamic description of the interfacial region must be coupled with a higher order modeling of the internal mechanical constraints as shown by [START_REF] Korteweg | Sur la forme que prennent les équations du mouvement des fluides si l'on tient compte des forces capillaires causées par des variations de densité considérables mais continues et sur la théorie de la capillarité dans l'hypothese d'une variation continue de la densité[END_REF], [START_REF] Rocard | The formation of vortices from a surface of discontinuity[END_REF] or [START_REF] Seppecher | Etude d'une modelisation des zones capillaires fluides: interfaces et lignes de contact[END_REF]. Using an Hamiltonian approach as in [START_REF] Casal | Relation entre l'équation de l'énergie et l'équation du mouvement en théorie de korteweg de la capillarité[END_REF] or the Virtual Power Principle as in [START_REF] Jamet | A numerical description of a liquid-vapor interface based on the second gradient theory[END_REF] it is possible to derive the new equations of motion for the fluid satisfying the first and second principles of thermodynamics:

∂ρ ∂t = -∇ • ρv ∂ρv ∂t = -∇ • ρv ⊗ v + pI + λ∇ρ ⊗ ∇ρ -τ d ∂ρE ∂t = -∇ • (ρE + p) v + λ (∇ρ ⊗ ∇ρ) • v + λρ∇ρ (∇ • v) -τ d • v + q (7a) (7b) (7c)
where v, p, τ d , q are respectively the fluid velocity, mechanical pressure, viscous stress tensor and energy flux.

Once given the expressions of τ d , q and the choice of an EoS, the previous system is closed. The mechanical pressure p and the pressure P EoS from the EoS are linked by Eq. ( 8) and this denomination, inherited from Jamet (1998), is justified by Eq. ( 7b) where p appears as the isotropic part of the non-viscous stress tensor.

p = P EoS - λ 2 (∇ρ) 2 -ρ∇ • (λ∇ρ) (8)
To simplify this study, the thermal conduction and the viscous constraints are not taken into account as their proper integration into the method presented in this study is still an ingoing work. Here, we primary focus on the mechanical and thermodynamic behavior of the interface. Phase change will not be treated as it would require a discussion beyond the scope of the current paper.

Canonical isothermal interfaces

For an isothermal 1D planar interface at equilibrium (v = 0) with no viscosity and thermal conduction, the system Eq. (7a) -Eq. (7c) reduces to the momentum equation in the form:

∂P EoS ∂x = λρ ∂ 3 ρ ∂x 3 (9)
By solving this equation at a chosen temperature T 0 with a stationary solver, it is possible to retrieve the profiles of thermodynamic variables for a given interface.

. . Figure 2 shows density profiles acquired trough the resolution of Eq. ( 9) for different temperatures using a Newton solver. The interface width, calculated with Eq. ( 28) is also provided for a wide range of reduced temperature T r = T /T c . The resulting curves in Fig. 2 raise the following observations. The values of the saturation densities ρ l and ρ v are not given a priori to the solver but result from the computations. Besides, these values match those one would find with equilibrium calculations based on chemical potential balancing. These values solely depend on the temperature T 0 and not on the capillary coefficient λ. The width of the interface directly depends on the values of T 0 and λ. This result was expected and will be explained analytically in the next section. If T 0 is chosen above the critical temperature of the fluid T c , the system independently converges towards a constant value for the density profile. This result concurs with the fact that no capillary effects spontaneously arise in trans/supercritical conditions given that, in these conditions, the heterogeneous system is no longer stable.

The results from the stationary solver can then be used as initial solutions for non stationary simulations

where the system Eq. (7a) -Eq. ( 7c) is solved in 1D and also act as a reference for the expected final solutions.

Starting from an interface given by the stationary solver, a stretching or a compression is applied to introduce a mass unbalance. The interface can be convected at a constant speed. Eventually, both the deformation and the convection can be combined. In all those cases, we verified that the interface behaves as expected and effectively goes back to its initial width. These simulation served as a first validation of both the model and our numerical implementation (the numerical setup used is detailed in Sec. 4.1).

Surface tension and interface thickness

2.1. Surface tension for capillary fluid 2.1.1. Thermodynamic approach . Considering a reference temperature T 0 < T c, and the associated saturation densities ρ l and ρ v , it is wellestablished, see [START_REF] Carey | Liquid Vapor Phase Change Phenomena: An Introduction to the Thermophysics of Vaporization and Condensation Processes in Heat Transfer Equipment[END_REF], that when the capillary forces are not accounted for, an unstable thermodynamic state instantaneously evolves toward the closest stable one between the liquid and vapor phases. For a 1D planar interface, this results in a discontinuity on the density profile between the values ρ l and ρ v . This is the description retained by the sharp interface methods. The position of that discontinuity is noted here x int and is used to represent the position of the interface. The corresponding density profile, noted ρ and defined by

ρ (x) = ρ v if
x int and ρ l if x > x int , will be used as a reference. Following the idea proposed in [START_REF] Cahn | Free energy of a nonuniform system. i. interfacial free energy[END_REF] one can express the surface tension σ of the interface taking into account the capillary forces. It corresponds to the extra free energy acquired by the diffuse interface comparatively to the discontinuous profile ρ. This free energy gap writes:

σ = x l xv F * -F dx = x l xv F (ρ * ) -F EoS (ρ) dx (10)
where ρ * (x) is the equilibrium density profile of the interface when the capillary forces are considered, bounded between the positions x v and x l (meaning that ρ * (x) = ρ v for x x v and ρ l (x) = ρ l for x x l ), F * = F (ρ * ) the corresponding volumetric free energy profile and F = F EoS (ρ) the volumetric free energy profile of the reference sharp interface. This definition of the surface tension is based on the notion of excess values from Edwards et al. x l x i x v

Volumetric free energy F excess value . The equilibrium profile of the interface will be the one minimizing the total free energy F of the system. However, since the total free energy of the canonical interface F is a fixed quantity, the minimization of F is equivalent to the minimization of F -F leading to:

F EoS (ρ) (sharp interface) F(ρ * ) (diffuse interface)
min ρ F (ρ) = min ρ x l xv F (ρ) dx ⇐⇒ min ρ F (ρ) -F = min ρ x l xv F (ρ) -F dx (11)
Calling I the integrand in the right integral of Eq. ( 11) and using the expression of F in Eq. ( 2) one can write:

I = F (ρ) -F = F EoS (ρ) + λ 2 ∂ρ ∂x 2 -F (12)
Following [START_REF] Margenau | The mathematics of physics and chemistry[END_REF] the integral minimization problem in Eq. ( 11) can be reformulated in terms of a partial derivative equation:

I - ∂ρ ∂x     ∂I ∂ ∂ρ ∂x     ρ = a (13)
Here a is a constant to be determined. For that, one can notice that (∂ρ/∂x) 2 is the only term in I in Eq. ( 12)

with a formal dependency on the density gradient ∂ρ/∂x. The relation in Eq. ( 13) is satisfied by the equilibrium profile. In any of the bulk phases, the density gradient is equal to zero. Both the reference sharp profile and the equilibrium profile have the same values of boundary densities ρ l and ρ v . Consequently, F EoS and F are also equal in the bulk phases, resulting in a = 0 which allows to write:

F EoS (ρ * ) -F = λ 2 ∂ρ * ∂x 2 (14)
This new relation can then be injected in Eq. ( 10) to finally obtain the expression of the surface tension:

σ = x l xv λ ∂ρ ∂x 2 dx ( 15
)
This result can be physically interpreted as follows. The minimization of the integral of I, leading to the equilibrium, will result from a balance between two terms. The first one is the integral of ∆F EoS = F EoS (ρ) -F EoS (ρ) and translates as the energetic cost necessary to diffuse the interface rather than to maintain it as a discontinuity. The second one is the integral of ∆F capi = λ/2 (∂ρ/∂x) 2 and is the energetic cost associated to the capillary forces themselves. To lower ∆F capi , the interface will tend to diffuse itself to reduce its gradient but this will increase ∆F EoS , the energy gap with the discontinuous profile, and vice versa. The equilibrium is achieved when these two terms are equal.

Relation between the thermodynamical and mechanical points of view

. Eq. ( 15) has also been obtained in [START_REF] Jamet | A numerical description of a liquid-vapor interface based on the second gradient theory[END_REF] using Eq. ( 9). In particular, given the differential relation dµ EoS = -s dT + dP EoS /ρ which simplifies in dµ EoS = dP EoS /ρ in the isothermal case, Eq. ( 9) can be written as:

µ EoS -µ sat (x) = λ ∂ 2 ρ ∂x 2 (x) (16)
where µ EoS is the chemical potential of the fluid as given by the equation of state and µ sat is the value of this chemical potential at saturation. A classic result from thermodynamics is that at equilibrium, µ sat has the same value in the liquid and the vapor phases. Thanks to the differential relation dF EoS = -S dT + µ EoS dρ where S is the volumetric entropy, which simplifies in dF EoS = µ EoS dρ, it appears that Eq. ( 16) and Eq. ( 14) are equivalent.

. These two equivalent differential equations ensure the uniqueness of the non-constant equilibrium interface profile for a chosen temperature (and a chosen λ). This in turn defines unique values for the surface tension σ and for the interface width h allowing to retrieve that both are only functions of the temperature at equilibrium. This link is further developed by the following derivations.

. Firstly, one can notice that from Eq. ( 16):

∂ ∂x ∂ρ ∂x 2 = 2 µ EoS -µ sat ∂ρ ∂x (17)
and then introducing Υ by:

Υ (ρ) = 2 ρ ρv µ EoS -µ sat ( ) d (18)
one can write the following equation, substituting the spatial variable κ with the density variable :

x xv ∂ ∂κ λ ∂ ∂κ 2 (κ) dκ = 2 x xv µ EoS -µ sat (κ) ∂ ∂κ dκ = 2 ρ ρv µ EoS -µ sat ( ) d (19) λ ∂ ∂κ 2 (x) = Υ (ρ (x)) (20)
The surface tension can therefore be linked to the chemical potential gap with:

σ = x l xv Υ (ρ) dx (21)
Here, Υ has been expressed as an integral that contains formally no dependency on the geometry of the interface but only on its intrinsic thermodynamic description through the expression of µ EoS (ρ). However, the integration that is necessary to obtain σ (Eq. ( 21)) shows a clear spatial dependency. Eq. ( 21) also shows that the surface tension depends on the energetic behavior of the interface through its propensity to store capillary energy driven by µ EoSµ sat (or likewise P EoS -P sat ) and on its inherent geometry mostly controlled by λ and the interfacial constraints. This duality will play a preponderant role in the rest of the study.

Determination of the capillary coefficient

185

. Given its role in Eqs. ( 14) and ( 16), it is clear that λ impacts the equilibrium profile of the interface. Since experimentally it is well established that σ and h solely depend on the temperature for pure fluids, λ should only depend on the temperature as well. This is confirmed by the expression derived by [START_REF] Rocard | The formation of vortices from a surface of discontinuity[END_REF]:

λ = ad 2 η RM 2 (22)
where R is the gas constant, M the molar mass of the fluid and a, η, d coefficients depending on the microscopic distribution of the molecular central forces of interaction, which is indeed directly linked to the temperature of the fluid. However, this formula is not of practical use for us given the difficulty to estimate a, η and d.

. Another way to evaluate λ is to use Eq. ( 15) backwards exploiting experimental measurements or analytical developments that allow to directly estimate the surface tension for pure fluids at a given temperature T . Knowing this reference value σ (T ) it is then possible to iterate on the value of λ to compute 1D planar interfaces until a match is reached using Eq. ( 15) to calculate σ (λ). Such a comprehensive analysis has been proposed in [START_REF] Lin | Gradient theory modeling of surface tension for pure fluids and binary mixtures[END_REF] using the Volume Translated Peng-Robinson EoS and has even been extended to binary mixtures.

We performed the same calculation using a simplier framework with the classical Soave-Redlich-Kwong EoS. For instance Fig. ( 4) shows the evolution of λ with the temperature for pure oxygen. The surface tension has been evaluated using the formula from [START_REF] Curl | Volumetric and thermodynamic properties of fluids -enthalpy, free energy, and entropy[END_REF] based on the corresponding state principle.

. . Both Lin's and our correlations give the same trend. For a wide range of temperature, λ can be considered constant with a limited loss of accuracy. However, the capillary coefficient diverges as the temperature approaches its critical value. Physically, this is explained by the fact that the interface keeps on widening as the temperature increases and is eventually infinitely wide once the regime becomes supercritical. For these temperatures, the constant λ hypothesis does not hold anymore. However, since our primary concern is to ensure that the interface is wide enough to be captured on the simulation mesh, it can still be a useful simplification. Besides, this allows to transit from subcritical to supercritical temperatures without loss of generality. Therefore, for the rest of the study and for the sake of simplicity, λ will be considered as independent of the temperature.

Influence of the capillary coefficient

2.3.1. Auto-similarity of the plane interface profiles . Let us consider a 1D planar interface at equilibrium for a fixed temperature T 0 and ρ 0 the density profile of said interface calculated with the capillary coefficient λ 0 . The corresponding EoS pressure profile is P 0 = P EoS ρ 0 .

The interface is supposed to be centered at x int = 0. Let φ be a strictly positive constant and X = x/φ. Since ρ 0 satisfies Eq. ( 9), one can write:

∀X, ∂P 0 EoS ∂x (X) = λ 0 ρ 0 (X) ∂ 3 ρ 0 ∂x 3 (X) (23)
Let ρ φ be the density profile defined by:

∀x, ρ φ (x) = ρ 0 (X = x/φ) (24)
which is an expansion if φ > 1 (or a compression otherwise) of the initial profile ρ 0 by a factor φ. One can also define the associated pressure profile P φ = P EoS ρ φ (x) = P 0 (X) and obtain the following relations:

∂P φ ∂x (x) = 1 φ ∂P 0 ∂x (X) ∂ 3 ρ φ ∂x 3 (x) = 1 φ 3 ∂ 3 ρ 0 ∂x 3 (X) (25) (26)
which, injected into Eq. ( 23), lead to:

∂P φ ∂x (x) = φ 2 λ 0 ∂ 3 ρ φ ∂x 3 (x) (27)
It clearly appears that the profile ρ φ is the solution of Eq. ( 9) where the capillary coefficient has been multiplied by a factor φ 2 . φ will play a major role in our thickening strategy and its impact will be detailed in the rest of the study. 210

Impact on the surface tension and the interface thickness

The previous result can be inverted by saying that solving Eq. ( 9), at a fixed temperature T 0 , with two capillary coefficients λ and φ 2 λ will give two density profiles, the second being an expansion (or compression) of the first one by a factor φ.

The thickness of an interface can be defined in multiple ways. In this study we choose the definition given by Eq. ( 28), which bears a meaning regarding numerical simulations since increasing or decreasing this thickness directly means changing how well the stiff gradients are resolved on a given mesh:

h = (ρ l -ρ v ) max |∇ρ| (28)
With this definition, calling h 0 the thickness of the interface calculated with a capillary coefficient λ 0 , Eq.

. To evaluate the surface tension, one can start by transforming Eq. ( 15), substituting the variables x and ρ in order to get an integral, of which limits do not depend on the geometry of the interface (which is affected by a modification of the capillary coefficient). For instance, ρ l and ρ v only depend on the temperature, and:

σ = x l xv λ ∂ρ ∂x 2 dx = x l xv λ ∂ρ ∂x ∂ρ ∂x dx = ρ l ρv λ ∂ρ ∂x dρ (29)
Once again, calling σ 0 and σ φ the surface tensions of the interfaces calculated with the capillary coefficients λ 0 and λ φ , using Eq. ( 29) gives:

σ φ = ρ l ρv λ φ ∂ρ φ ∂x dρ = ρ l ρv φ 2 λ 0 1 φ ∂ρ 0 ∂x dρ = φ ρ l ρv λ 0 ∂ρ 0 ∂x dρ = φσ 0 (30)
Finally, these two results can be summarized by the following dependencies:

h ∝ √ λ σ ∝ √ λ (31)
These correlations have been retrieved analytically near the critical point in [START_REF] Jamet | Etude des potentialités de la théorie du Second Gradient pour la simulation numérique directe des écoulements liquide-vapeur avec changement de phase[END_REF] and numerically by [START_REF] Gaillard | Interfaces diffuses et flammes transcritiques LOX/H2[END_REF] and are confirmed in Fig. ( 5). In both series of curves displayed, the slope is 1/2 in logarithmic scale, as expected from Eq. ( 31). However, the results previously shown rely on very fine meshes and the width of the interfaces is still far under the size we can afford to resolve for industrial applications given the current available computational power. From . Given this observation, even for DNS configurations, the interface should be thickened to be resolved on more realistic meshes, even for academic applications. A first idea would be to increase the capillary coefficient λ.

However the relations given in Eq. ( 31) showed that while the interface would indeed be thickened, the surface tension would also be augmented modifying the macroscopic behavior of the interface. We present here a novel method in oder to thicken the interface while avoiding this major drawback, the objective being to reach typical DNS mesh sizes without modifying the surface tension.

. It should be noticed that the issue of thickening has already been addressed by Jamet. He proposed two different strategies of which principles are the foundations of the method we proposed in this study. In his first strategy, [START_REF] Jamet | Etude des potentialités de la théorie du Second Gradient pour la simulation numérique directe des écoulements liquide-vapeur avec changement de phase[END_REF] shows that Eq. ( 16) has an analytical solution. This allowed him to express both the interface width and the surface tension with:

       h = 4 ρ l -ρ v λ 2A (32a) σ = (ρ l -ρ v )
. The physical meaning of coefficient A can be found in [START_REF] Jamet | The second gradient method for the direct numerical simulation of liquid-vapor flows with phase change[END_REF]. It is related to the isothermal compressibility in the bulk phases. In Eqs. (32a) and (32b), the dependencies h ∝ √ λ and σ ∝ √ λ previously found are visible again. The key point here is the dependency of the interface width h and the surface tension σ with A, namely h ∝ 1/ √ A and σ ∝ √ A. This reversed dependency allows to play on both λ and A to obtain any desired values for σ and h simultaneously. This method, although simple to apply, presents three main drawbacks.

Firstly, in order to apply the relations in Eqs. (32a) and (32b) the values of the saturation densities ρ l and ρ v must be obtained beforehand. This is however achievable using for instance a thermodynamic equilibrium solver.

Secondly, the range of application for this method is extremely narrow around the critical conditions. In practice this method seems to fail for reduced temperatures below T r = 0.99 and therefore is not suited for industrial application in its current formulation. Lastly, modifying the coefficient A leads in fact to a modification of the compressibility coefficients in the bulk phases and may jeopardize the thermodynamic consistency since it may in turn modify the heat capacities and the sound speed. Through this study, Jamet showed that a second coefficient could be introduced to modulate the macroscopic behavior.

. Following this first idea, [START_REF] Jamet | The second gradient method for the direct numerical simulation of liquid-vapor flows with phase change[END_REF] proposed a more complex and complete procedure, by modifying directly the functional Ψ (ρ) =ρ = ∂ρ/∂x in the interface. Given this definition of Ψ, he derived the following six equations, the meaning of which is discussed hereunder:

Ψ (ρ v ) = 0 Ψ (ρ l ) = 0 dΨ dρ (ρ v ) 2 = 1 λ 1 ρ v dP dρ (ρ v ) dΨ dρ (ρ l ) 2 = 1 λ 1 ρ l dP dρ (ρ l ) max Ψ = ρ l -ρ v h λ ρ l ρv Ψ (ρ) dρ = σ (33a) (33b) (33c) (33d) (33e) (33f)
One must notice that the derivative dP / dρ is directly given by the chosen equation of state (T being fixed for the study) and its definition is not modified by the Second Gradient theory. First Eqs. (33a) and (33b) account for the fact that there is no density gradient in the bulk phases. Eqs. (33c) and (33d) impose the thermodynamics to be differentiable at the boundaries between the bulk phases and the binodal region. Finally Eqs. (33e) and (33f) link Ψ to the interface width as defined by Eq. ( 28) and the surface tension as expressed in Eq. ( 15)

. This method ensures the thermodynamic consistency of the system. In particular the heat capacities and the latent heat are continuous throughout the interface. However, one must find a profile function Ψ that satisfies all these constraints, which reveals to be intricate for Eqs. (33c), (33d) and (33e) where the use of a Newton-like solver is required, at each time step on every mesh points. Using this method, [START_REF] Jamet | The second gradient method for the direct numerical simulation of liquid-vapor flows with phase change[END_REF] studied several typical capillary phenomena such as the disappearance or the coalescence of bubbles and the dynamical hysteresis in the contact angle with a wall.

. This second thickening method is both powerful and consistent. However we spotted two major drawbacks when it comes to its application in our targeted simulations. Firstly, the method requires to geometrically locate and qualify the interface in order to be applied because it relies on a local geometrical modification of the interface shape, while the first method only needed the values of the saturation densities. Secondly, most of the thermodynamic variables have to be recalculated multiple times at every step for every point because of the calculations required by the Newton-like solving procedure. One could argue that some of the results could be tabulated to save computational time, however this tabulation would be mixture-dependent which also constitutes a drawback regarding future complex simulations.

The Thickened Interface Method (TIM)

Foundations of the method

Changing uniformly the value of λ has been shown to be insufficient to thicken the interface without augmenting the surface tension. The addition of a second parameter has proven to be a good idea in order to choose both σ and h as in Jamet's first thickening strategy. Using this paradigm, the purpose is therefore to derive a method that allows such flexibility. Moreover, in his second method, Jamet showed the necessity to consistently modify the thermodynamics jointly with the mechanics in the interfacial region. Essentially, we would like to introduce here a coefficient φ s so that h ∝ √ φ s and σ ∝ 1/ √ φ s (the same effect as when modifying A in Eqs. (32a) and (32b)). The impact of φ s should allow to mitigate the effect of modifying λ on both h and σ. From Eq.( 21), reducing σ by a certain factor can be done in a thermodynamically consistent manner by dividing Υ by a constant.

Moreover, from the definition given by Eq. ( 18), a straightforward way to do so is to divide a priori µ EoSµ sat by this constant. This modification has the advantage of not being performed on λ and therefore should avoid the associated drawbacks. As done to address the impact of a modification of λ in 2.3, the impact of the thermodynamic modification in Eq. ( 34) is studied in the next paragraph, where η is a constant of which value will be chosen depending on the desired outcome.

µ EoSµ sat → µ EoSµ sat /η (34)

Modification of the equation of state

The modification introduced by the substitution in Eq. ( 34) presumably allows to diminish the surface tension of the interface in a controlled manner. However, by itself, it suffers two shortcomings. Firstly, this modification inherently assumes the interface to be isothermal as µ sat is actually a function of the temperature. Secondly, this modification of the chemical potential must be carried on and applied to the other thermodynamic variables to achieve a consistent model. In essence, Eq. ( 34) is to be understood as a modification of the EoS describing the fluid. In the following, we present the expression of the thermodynamic variables induced by this modification of 280 the EoS. Here, for a variable χ, χ EoS still refers to the variable as provided directly by the initial EoS and χ η to the expression with the modified EoS induced by Eq. ( 34).

. The starting modification is applied upon the chemical potential. Its full expression showing density and temperature dependencies reads:

µ η (ρ, T ) = µ sat (T ) + µ EoS (ρ, T ) -µ sat (T ) η (35)
The chemical potential µ η and the volumetric free energy F η are linked by µ η = (∂F η /∂ρ) T . Enforcing the condition F η (ρ v (T ) , T ) = F EoS (ρ v (T ) , T ) and noticing F EoS (ρ v (T ) , T ) = ρ v (T )µ sat (T )P sat (T ), F η can be expressed by:

F η (ρ, T ) = F sat * (ρ, T ) + F EoS (ρ, T ) -F sat * (ρ, T ) η ( 36
)
where F sat * (ρ, T ) = ρµ sat (T ) -P sat (T ) (the same expression is obtained enforcing

F η = F EoS in the liquid).
The new EoS pressure is then defined by P η = ρµ η -F η . Given Eqs. ( 35) and (36), one gets:

P η (ρ, T ) = P sat (T ) + P EoS (ρ, T ) -P sat (T ) η (37)
The specific entropy is defined by s η = -(1/ρ)(∂F η /∂T ) ρ . Using Eq. ( 36) it can be written:

s η (ρ, T ) = s (T ) + s EoS (ρ, T ) -s (T ) η + 1 - 1 η 1 ρ - 1 ρ (T ) dP sat dT (T ) ( 38
)
where the notation χ (T ) designates the saturation value of a thermodynamical variable χ EoS calculated indifferently in the liquid or the vapor phase for the temperature T as long as the same criteria is used to chose this 285 value for all the thermodynamic variables.

Eventually, using e s η = µ η + T s η -P η /ρ, the specific internal energy can also be expressed by: e s η (ρ, T ) = ēs (T ) + e s EoS (ρ, T )ēs (T )

η + 1 - 1 η 1 ρ - 1 ρ (T ) T dP sat dT (T ) -P sat (T ) (39)
Practically, if for some computational node the value of χ (T ) is chosen in the liquid phase then s (T ), μ (T ), ēs (T ), etc... also have to be chosen in the liquid phase and vice versa. Noticeably though, the criterion to decide in which phase the value of χ is evaluated can differ from one point to the other in the computational domain and can change from one iteration to the other. The only condition is the local consistency of this choice throughout the calculation of the thermodynamic variables at one point and at one instant. Conveniently, this lead to the useful following results (where ρ k = ρ l or ρ v ) ensuring the continuity of the thermodynamics through the saturation bell independently of η:

lim ρ→ρ k e η (ρ, T ) = e EoS (ρ k (T ) , T ) lim ρ→ρ k s η (ρ, T ) = s EoS (ρ k (T ) , T ) (40)
(41)

The expression of other variables such as the specific isochoric heat capacity C V η , thermal expansion coefficient α η , the isothermal compressibility coefficient β η and the specific isobaric heat capacity C P η are also to be derived.

The aspect of P η , µ η , e s η and s η are given in Fig.

Impact on the isothermal interface

Now that all the thermodynamic variables have been expressed for the modified EoS, the momentum equation can be written for an isothermal interface at equilibrium, still using a reference capillary coefficient λ 0 . From Eq.

(9), it reads:

∂P η ∂x = λ 0 ρ η ∂ 3 ρ η ∂x 3 (42)
Here, we took care to reference the density profile by ρ η as it mays also be impacted by the EoS modification.

Since for an isothermal interface P sat (T ) is a constant, the previous equation simplifies in

∂P EoS (ρ η ) ∂x = ηλ 0 ρ η ∂ 3 ρ η ∂x 3 (43)
From the study presented in 2.3 and in particular Eq. ( 27), one can recognize from Eq. ( 43) that ρ η is the density profile solution of Eq. ( 9) with a capillary coefficient ηλ 0 , that is to say, the reference interface for a capillary coefficient λ 0 thickened by a factor √ η. In particular, this means that using the variable switching

X = x/ √ η, one has ρ η (X) = ρ (x)
. Moreover, Eq. ( 18) and Eq. ( 35) allow to write Υ (ρ η ) = Υ ρ 0 /η. With that, the new value of σ can also be obtained using Eq. ( 21), it gives:

σ η = Xv X l Υ (ρ η ) d X = xv x l Υ ρ 0 (X) √ η dx = σ 0 √ η (44)
The impact of the EoS modification induced by Eq. ( 34) and consistently carried over the other thermodynamic variables in the previous paragraph has been established for h and σ. The dependency on η observed for both values matches the one targeted for the parameter φ s . Therefore, this modification of the EoS will be adopted for our method by formally using φ s in place of η in the new equation of state (and the notation χ η by χ φs for a variable χ).

Complete modifications of the TIM

By itself, the EoS modification is not sufficient to achieve the desired action for the TIM. However, by introducing the coefficient φ l applied to the capillary coefficient such as λ = φ l λ 0 , one can obtain the dependencies h ∝ √ φ l and σ ∝ √ φ l . When combined with the modification associated with φ s , the final dependencies read:

h ∝ φ l φ s σ ∝ φ l /φ s (45)
As expected, the combined action of modifying the capillary coefficient and the EoS in the interface makes it possible to act independently on h and σ. An additional precaution as to be taken to complete the method.

Indeed, for non isothermal cases, the TIM must leave the behavior of the fluid unchanged in the bulk phases. In that respect, all the modifications of the thermodynamics and mechanics of the interface are to be restrained to the binodal region only. The final expression of the variables in the TIM method, refereed to with the exponent TIM , can be detailed. First, the variables that are not impacted by capillary contributions, such as P , s, µ, etc..., only endure the modification from the EoS. For instance, the new EoS pressure and specific entropy are given by:

P TIM =         
P sat + P EoS -P sat φ s = P φs in the binodal region

P EoS otherwise and s TIM =         
s φs in the binodal region

s EoS otherwise
As previously mentioned, the coefficient φ l is applied to the capillary coefficient leading to the expression:

λ TIM =          φ l λ 0 in the binodal region λ 0 otherwise
Therefore, thermodynamic variables impacted by both λ and the EoS, such as P th , p, e s , etc... endure both modifications. For instance, the new mechanical pressure and specific internal energy are given by

p TIM =          P φs - φ l λ 0 2 (∇ρ)
2 + 2ρ∆ρ in the bin. reg.

P EoS - λ 0 2 (∇ρ) 2 + 2ρ∆ρ otherwise and e s TIM =          e s φs + φ l λ 0 2ρ (∇ρ)
2 in the bin. reg.

e s EoS + λ 0 2ρ (∇ρ) 2 otherwise
All things considered, the new Euler equations can be formulated in the framework of the TIM. Practically, no formal modification is applied to the set Eqs.( 7a)-(7c). It is only the global thermodynamic modification of the TIM that has to be included in order to get the new system. This is achieved by substituting p, λ 0 and E with p TIM , λ TIM and E TIM = e s TIM + v 2 /2, to get the following equations:

∂ρ ∂t = -∇ • ρv ∂ρv ∂t = -∇ • ρv ⊗ v + p TIM I + λ TIM ∇ρ ⊗ ∇ρ ∂ρE TIM ∂t = -∇ • ρE TIM + p TIM v + λ TIM (∇ρ ⊗ ∇ρ) • v + λ TIM ρ∇ρ (∇ • v) (46a) (46b) (46c)
This result is of a high importance from a numerical point of view because it implies that only the thermodynamic routines of the code have to be modified to transit from the native Second Gradient to the TIM-modified 300 system, thus greatly reducing the implementation difficulties.

Additional considerations on the TIM

As both Jamet's thickening methods, the Thickened Interface Method (TIM) requires the calculation of the saturation values (pressure and densities) together with a manner to discriminate the stable/unstable nature of the current thermodynamic regime and detect whether the modification is to be applied or not. This is done using a thermodynamic equilibrium solver: the saturation values are retrieved by equalizing the chemical potential of the liquid and vapor states at a given temperature.

. Finally, since most of the time one wants to keep the surface tension unmodified while applying a thickening factor F to the interface, the dependencies Eq. ( 45) becomes φ l /φ s = 1 and √ φ l φ s = F. It leads to the following expressions for φ l and φ s which will remain unchanged in the rest of the study:

φ l = φ s = √ F (47)
4. Numerical validations of the TIM

Numerical setup

The Second Gradient theory model has been implemented in the semi-academic code AVBP jointly developed by CERFACS, IFPEN with specific modules from EM2C laboratory. It is designed to solve the Euler/Navier-Stokes equations on 2D/3D unstructured meshes for compressible flows with real gas thermodynamics. We used the SRK EoS given Eq. ( 1), a third order Galerkin spatial discretization and a second order Runge-Kutta time integration. The Navier-Stokes Characteristic Boundary Conditions from [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF] are used at inlets and outlets.

Validation in isothermal configurations

Isothermal equations

Isothermal configurations allow to work with controlled saturation values and comparisons with analytical results from the literature. For these reasons, first validation cases have been carried out in isothermal configurations for which the system (46a)-( 46c) can be simplified since the energy equation is not solved and P sat is constant in the domain, implying ∇P * = ∇ P EoS / √ F in the binodal region. This isothermal hypothesis does not bear any meaning outside the binodal region. Although points outside the binodal region may arguably be encountered in the simulations, their impact on the interface dynamics is assumed to be marginal. Therefore for isothermal cases, all points of the computational domain are treated as laying in the binodal region. This pushes the simplification further, implying that in all the domain λ * = φ l λ 0 and P * = P sat + P EoS -P sat / √ F. The system Eq. ( 46) reduces to:

∂ρ ∂t = ∇ • [ρv] ∂ρv ∂t = -∇ • ρv ⊗ v + P EoS √ F I + √ Fλ 0 - 1 2 (∇ρ) 2 + ρ∇ • (∇ρ) I + (∇ρ ⊗ ∇ρ) (48a) (48b)

Validation in 1D

As for the native Second Gradient model, the TIM has first been experimented on using 1D convected and/or deformed interfaces for several thickening factors. The result of one such case is shown in Fig. ( 8). As desired, the TIM does not alter the behavior of 1D isothermal interfaces. The method allows a convection of the interface with a return to equilibrium without any noticeable discrepancy. mechanically enlarged at T = 139.12 K with λ = 1.0 • 10 -10 m 7 .kg -1 .s -2 and convected at u 0 = 2.5 m.s -1 . The reference is given by ( )

Case design

. To validate the method in 2D configurations, the chosen case is an oscillating plane interface initially deformed by a harmonic longitudinal perturbation as illustrated in Fig. ( 9) for low thickening factors. The wavelength of the perturbation is reduced to one period in the domain so that the local curvature radius of the interface remains substantially higher than the interface thickness. Moreover, the amplitude A of the perturbation has been kept low in order to remain in the linear deformation regime. With these settings, the interface is expected to oscillate almost indefinitely since no dissipation, apart from the numerical diffusion, is integrated in the equations. From [START_REF] Fyfe | Surface tension and viscosity with lagrangian hydrodynamics on a triangular mesh[END_REF] one gets that the expected frequency of the oscillations ω obeys Eq. ( 49) where ρ u (resp. ρ d ) and h u (resp. h d ) are the density and width of the upper (res. lower) fluid, σ the surface tension of the interface, g the gravity acceleration and k the wavenumber of the deformation. Noticeably, it does not depend on the amplitude of the initial deformation.

ω 2 = (ρ d -ρ u ) gk + σk 3 ρ d coth kh d + ρ u coth kh u (49)
In our configuration, no gravity is considered, the upper and lower densities are the saturation densities ρ v and ρ l , only one period of a sine is used for the deformation therefore k = 2π/l where l is the domain width i.e. the length of the interface when no perturbation is applied and the upper and lower widths h u , h d have been chosen so that h u = h d = 1.5l. The simplified expression for the period is now given by Eq. ( 50).

τ = l 3/2 σ 1/2 coth (1.5) (ρ l + ρ v ) 2π (50)
. Using Eq. ( 50) we designed test cases, exposed in Fig. ( 9), in order to fulfill the previously stated hypothesis.

First, we established a reference case on a numerically affordable mesh, the amplitude A ≈ 0.05l allowed for substantial oscillations without stepping out of the linear domain. The other cases were then obtained by applying an expansion δ = F to the whole geometry while simultaneously thickening the interface by a factor F. The relative amplitude of the imposed deformation (5%) remained unchanged.

. In this figure, for the sake of clarity, the relative scales of the different configurations are respected and thus represented for small factors only. However, for the results presented hereafter, the methodology is used for factors up to F = 10 4

F = 1. 0 F = 1. 5 F = 2. 0 F = 2. 5 l 0 F • l 0
. Since we ensured that all the cases were homothetic, they are expected to mimic the reference case. Beside, by applying the same expansion to the mesh step size we ensured a constant computational time for the same number of iterations throughout all cases. Though at first sight this design is expected to alter the oscillation period, Eq. ( 50) allows to predict the evolution of said period with the modification of the geometry. This strategy was eventually adopted to validate the 2D implementation of the Second Gradient model in the solver but also the behavior of the TIM regarding the dynamics of the interface.

Results

. For these simulations, in addition to the settings listed in the previous paragraphs, a high order conservative selective filter, inspired by [START_REF] Mathew | An explicit filtering method for large eddy simulation of compressible flows[END_REF], has been applied to suppress the point to point numerical . From Eq. ( 31) one knows that thickening the interface by a factor F using only the coefficient φ l will result in the same factor being applied to the surface tension and therefore, from Eq. ( 50) it can be concluded that the period should behaves as shown by Eq. ( 51). This expected result is symbolized by the dashed line in Fig. ( 14).

The adjective "inconsistent" will be used in the following to qualify this method as it should not conserve the value of σ.

τ φ l τ 0 = (Fl 0 ) 3/2 (Fσ 0 ) 1/2 σ 0 1/2 l 0 3/2 = F (51)
. The same a priori estimation of the period evolution can be done with Eq. ( 50) for the TIM supposed to thicken the interface by a factor F without modifying its surface tension, obtaining Eq. ( 52). This expected behavior is symbolized with a solid line in Fig. (14).

τ TIM τ 0 = (Fl 0 ) 3/2 σ 0 1/2 σ 0 1/2 l 0 3/2 = F 3/2 (52)
. . The two behaviors obtained with Eqs. ( 51) and ( 52) are accurately retrieved in the simulations. The results for the "inconsistent" thickening method using only φ l further demonstrate the predictive capabilities of the native . Given the manner the TIM has been derived, its application is not restrained to pure fluid or to temperatures extremely close to T c . Besides, the shape of the pressure/chemical potential is not analytically prescribed in the binodal region. This ensures a straightforward applicability to mixtures for which the pressure/chemical potential cannot be obtained analytically due to the complex dependencies with the mixture fractions.

Validation in 2D non-isothermal configurations

Non isothermal oscillating interface

As a first step towards non-isothermal simulations, the same previous cases have been computed, using the same isothermal initial solutions but this time by solving the complete system of equations as expressed in Eqs. The first noticeable outcome is the curves superimpose satisfyingly, meaning that the correlation for the surface tension preservation from Eq. ( 52) is still observed by the model. This results is also visible in the left graph of The second noticeable result is the damping of the interface oscillations. The entropy balance in the domain shows that the second principle is respected, as expected from the model since capillary terms are non-dissipative by nature. Plus, the phenomenon is presumably not caused by the thickening procedure of the TIM since the non-thickened case displays the same behavior. The mechanism behind this damping of the interface oscillations is still under investigation and lack a proper physical justification as, to the best of our knowledge, no such cases without physical dissipation have been produced in the literature. Nonetheless, the overall behavior of the interfaces remains coherent and the surface tension is still conserved.

. To widen the range of application of the method, simulations with non-isothermal initial solutions have been carried out. The same framework as in Sec. 4.2.3 has been used however this time a non-constant temperature evolution is imposed outside of the interfacial zone so as to get an actual non-isothermal configuration as shown in Fig. ( 17). Essentially, the liquid is over-cooled while the gas is overheated to create bulk phases outside the edges of the binodal region. The initial solutions are designed so as to get a constant pressure outside the interface, equal to the saturation pressure associated to the interface temperature, to prevent acoustics driven movements.

. Notwithstanding the mentioned difficulties, the period of oscillation is found to follow the trend τ ∝ F 1.44 which is close to the canonical value F 1.5 of the isothermal cases. Overall, these simplified non-isothermal cases showcased a good agreement with the theory and further validated the TIM methodology. non-isothermal and the non-isothermal system is solved.

Periodic liquid jets

For the second case, we introduced more complexity by considering a symmetrical mixing layer that takes the aspect of a liquid jet in its own vapor. The interface is set at T i = 100 K, the liquid is slightly supercooled while the vapor is superheated and the initial solution is isobaric at P sat = 7.85 bar. The surrounding vapor is initially at rest whereas the liquid has a constant descending vertical velocity v l . The simulations are performed on a regular cartesian mesh. The domain, as depicted in Fig. For W e 1 = 50, the jet faces a strong macroscopic deformation with no clear breakup. From [START_REF] Baillot | Behaviour of an air-assisted jet submitted to a transverse high-frequency acoustic field[END_REF], this is expected for a regime between the Rayleigh non-axisymmetric and the shear breakup ones. Macroscopic structures form and eventually start separating and coalescing to form a small number of big droplets. The case W e 2 = 100 starts to display macroscopic ligaments. Droplets are also generated but they are quickly absorbed by coalescence in the larger structures. These observations are consistent with a regime between the shear breakup and the membrane regimes. With the case W e 3 = 200, a great modification occurs in the jet topology. Long macroscopic ligaments form and quickly break-up to form droplets among which the biggest also tend to experience a secondary break-up. This is reminiscent of the membrane regime as illustrated in [START_REF] Baillot | Behaviour of an air-assisted jet submitted to a transverse high-frequency acoustic field[END_REF].

In all cases, a small diffusion of the interface can be noted and is suspected to be caused by the filter used in the simulation. Overall, the three test cases qualitatively display the expected behavior in agreement with experimental studies. In particular the atomization regimes observed concur adequately with the results in [START_REF] Baillot | Behaviour of an air-assisted jet submitted to a transverse high-frequency acoustic field[END_REF].

Conclusions and future work

In this study, we confirmed the predictive capabilities of the Second Gradient theory regarding the simulation of real gas flows under subcritical conditions where capillary effects lead to the presence of interfaces. It is possible to retrieve the saturation for a given temperature, matching the values otherwise obtained through classical equilibrium considerations as well as the interface width and the surface tension.

. The TIM, a thermodynamically consistent thickening method, has been proposed and should allow to resolve interfaces on DNS (in the sense of turbulence) grids for a wide range of temperatures representative of cryogenic rocket engines at a relatively low computational cost. It preserves the interface dynamics through the conservation of its surface tension.

All the validation cases have been performed using the AVBP code solving the compressible Euler equations with real gas thermodynamics further cementing the applicability of both the Second Gradient theory and the TIM for more complex academic purposes and industrial configurations. The isothermal cases displayed the expected mechanical and thermodynamic behaviors, confirming the consistency of the model. The application to non isothermal cases also proved quite satisfactory for cases approaching that of an actual liquid jet injection. The results demonstrated the capability of the model to simulate such cases.

. Some relevant questions still have to be treated and stand as primary concerns regarding our future work. The results for non-isothermal cases showed a damping, even for non-thickened cases, that remains unexplained, in particular considering the fact that capillary terms are dissipation-free in nature. Additionally, the effect of the TIM on transitory phenomena should be further investigated, in particular on the characteristic relaxation time of the interface in 1D and 2D configurations. Moreover, the interaction of the method with thermal conduction and viscosity must be addressed to resolve the full set of Navier-Stokes equations and get closer to practical configurations. Finally, the behavior of the method during the transition between subcritical and supercritical regimes should be looked at since the applicability of the method to such a thermodynamic change has been a

driving concern for this study.

Figure 1 :

 1 Figure 1: Clapeyron's diagram of pressure versus specific volume for pure nitrogen N 2 with isothermal curves for different thermodynamic regimes calculated using the SRK equation of state

Figure 2 :

 2 Figure 2: Examples of density profiles for isothermal interfaces at different temperatures (left) and evolution of the interface width h with the reduced temperature Tr = T /Tc (right). Nitrogen N 2 with λ = 1.0 • 10 -10 m 7 .kg -1 .s -2

Figure 3 :

 3 Figure 3: Schematic representation of the excess volumetric free energy F

]Figure 4 :

 4 Figure 4: Evolution of the capillary coefficient λ with the temperature for oxygen O 2 using SRK equation of state

Figure 5 :

 5 Figure 5: Evolution of the interface thickness h and surface tension σ with the capillary coefficient λ for oxygen O 2 at different reduced temperature T /Tc in logarithmic scales

Fig. ( 6

 6 Fig. (6) the expected width using the Second Gradient theory would be comprised between a few nm and a few

Figure 6 :

 6 Figure 6: Evolution of the interface width h with the temperature T as predicted by the Second Gradient theory for pure oxygen O 2 , Tc = 158.54 K

Figure 7 :

 7 Figure 7: Pressure from EoS, chemical potential, specific internal energy and speficic entropy of nitrogen N 2 at T = 117 mathrmK as modified by Eq. (34) for different values of η

Figure 8 :

 8 Figure 8: Density (up) and pressure (down) profiles for the relaxation of a convected thickened (F = 1000) oxygen interface initially

Figure 9 :

 9 Figure 9: Evolution of the initial density profiles for interfaces with different thickening factors, a fixed relative interface perturbation and a proportionally expanded geometry.

  Fig. (13).

Figure 12 :Figure 13 :

 1213 Figure 12: Zoom on the normalized density gradient profiles at different instants for the oscillation of an initially harmonically perturbed plane nitrogen N 2 interface with T init = 119.88 K, λ 0 = 1.0 • 10 -16 m 7 .kg -1 .s -2 and F = 50

Figure 14 :

 14 Figure14: Evolution of the interface oscillation period normalized by T F=1 , with the thickening coefficient F (in log scale) for consistent and inconsistent thickening strategies, comparison with numerical results using the inconsistent and TIM methods

(Figure 15 :

 15 Figure15: Time evolution of the density at the center of the domain for the oscillation of an initially harmonically perturbed plane nitrogen N 2 interface with T i = 100 K, λ 0 = 1.0 • 10 -16 m 7 .kg -1 .s -2 . The intiial solution is isothermal but the non-isothermal system is solved for different thickening factors

Fig. ( 16

 16 Fig. (16) where the nondimensionalized periods have been extracted for different thickening factors.

Figure 16 :

 16 Figure 16: Evolution of the normalized interface oscillation period, with the thickening coefficient F (in log scale) for different types of initialization and equations solving

Figure 17 :

 17 Figure 17: 1D density (left) and temperature (right) profiles for isobaric nitrogen N 2 interface with P sat = 7.85 bar and λ 0 = 1.0 • 10 -16 m 7 .kg -1 .s -2 A central density oscillation curve for this cases is given Fig. (18) for F = 10. The same behavior of a damped

Figure 18 :

 18 Figure18: Time evolution of the density at the center of the domain for the oscillation of an initially harmonically perturbed plane nitrogen N 2 interface with T i = 100 K, λ 0 = 1.0 • 10 -16 m 7 .kg -1 .s -2 thickened by a factor F = 10. The initial solution is

Figure 19 :

 19 Figure 19: Schematic representation of the density field for the simulation configuration

  134 mN/m N y = N x /2 600 While the density, pressure and temperature profiles are the same for all cases, three different values have been used for the liquid velocity so as to get three different Weber numbers W e = ρ v v 2 l D/σ typical of three different atomization regimes (see Lasheras & Hopfinger (2000)): v 1 = 17.9 m/s for W e 1 = 50, v 2 = 25.4 m/s for W e 2 = 100 and v 3 = 35.9 m/s for W e 3 = 200. Results are given in Fig. (21).

  W e = 50, t = 0µs W e = 100, t = 0µs W e = 200, t = 0µs W e = 50, t = 4µs W e = 100, t = 4µs W e = 200, t = 4µs W e = 50, t = 8µs W e = 100, t = 8µs W e = 200, t = 8µs W e = 50, t = 12µs W e = 100, t = 12µs W e = 200, t = 12µs W e = 50, t = 16µs W e = 100, t = 16µs W e = 200, t = 16µs W e = 50, t = 24µs W e = 100, t = 24µs W e = 200, t = 24µs W e = 50, t = 32µs W e = 100, t = 32µs W e = 200, t = 32µs 50 100 150 200 250 300 350 400 450 500 550 600 650 700 Density [kg.m -3 ]

Figure 21 :

 21 Figure 21: Evolution in time of the density profiles for two-dimensional periodic Nitrogen liquid jets in their vapor. Three cases have been carried out for three different Weber numbers: W e = 50 (left), W e = 100 (center) and W e = 200 (right)

Table 1 :

 1 Simulation parameters

	T i	100.0 K	h	1.31 µm	∆ x = ∆ y	0.133 µm
	λ 0 19.6 mm 7 /g/s 2	D	20 µm	L y = L x /2	90.0 µm

(25) applied to the densities ρ 0 and ρ φ gives that max ∇ρ φ = max ∇ρ 0 /φ for the interface calculated with λ φ = φ 2 λ 0 . Therefore, with Eq. (28) it leads to h φ = φh 0 .

√ 2Aλ (32b)

Acknowledgements

The authors acknowledge funding from ANR through project ANR-14-CE22-0014 (SUBSUPERJET). This work was granted access to the HPC resources of CINES made available by GENCI (Grand Equipement National de Calcul Intensif) under the allocation x20172b6176. A part of this work was performed using HPC resources from the mésocentre computing center of Ecole CentraleSupélec and Ecole Normale Supérieur Paris-Saclay supported by CNRS and région Île-de-France. The authors would also like to thank D. Jamet, R. Saurel and M. Pelletier for fruitful discussions and insightful inputs.

Second Gradient model and its proper implementation in the code AVBP. Likewise, the results obtained with the TIM offer an additional validation for its mechanical consistency. Indeed, when no thickening is applied, the period expected from Eq. ( 50) is τ F=1th = 17.14µs while the period given by the simulation is τ F=1sim = 17.89µs leading to a satisfactory error below 5%. The error between the theoretical and numerical periods remain below 5% for all the thickening coefficients considered.