N Rutard

L.-H Dorey

C Le Touze

S Ducruix

Large-eddy simulation of an air-assisted liquid jet under a high-frequency transverse acoustic forcing

The present contribution lies on the high-frequency combustion instability issue occurring in liquid rocket engines, and focuses on subcritical operating conditions. Thanks to the development of high-performance computing and advanced numerical simulation codes, large-eddy simulation has become an essential analysis tool in the domain. Different unsteady numerical strategies exist to handle the wide range of time and length scales involved in two-phase flows related to subcritical operating conditions, but only a few are adapted to deal with reactive flows. What is more, most of them are developed without any consideration for possible acoustic effects on atomization processes while these may have a great impact on the stability of the engine. Therefore, the present contribution aims at evaluating the ability of one of these numerical strategies in particular to render acoustic effects on atomized liquid jets typical of what happens under unstable operating conditions. To do so, the numerical simulation of a non-reactive two-phase flow submitted to a transverse acoustic modulation is performed. Two cases are simulated, the first one with acoustic modulation, the jet lying in the vicinity of a transverse intensity anti-node, and the second one without any. In order to highlight the effect of acoustics on the liquid phase, the two cases are compared to one another as well as to experiments presented in the literature. The numerical strategy used in these simulations is based on the coupling between a diffuse interface method for the simulation of large liquid structures, and a kinetic-based Eulerian model for the description of droplets. It is found that the flattening of the liquid core under acoustic constraints is retrieved in the simulations. This flattening process thus induces a decrease of the liquid core length and an intensification of its stripping. The additional mass of liquid thus ripped from the liquid jet is transformed into droplets and the spray undergoes a drastic modification of its shape thanks to an appropriate primary atomization source term. Finally, periodic oscillations of droplets under the acoustic velocity are rendered thanks to the coupling between the gaseous phase and the spray through a drag source term. The numerical strategy thus appears adapted to deal with numerical simulations of coaxial two-phase flows under transverse acoustic modulation, in which the response of the gas, the liquid core and the spray are all linked to one another. This strategy can then be used for future numerical studies of high-frequency combustion instabilities.

Introduction

The occurrence of High-Frequency (HF) combustion instabilities is a major issue in the design of Liquid Rocket Engines (LRE). In such high power density engines, the coupling between the injection system, the propellant jets, their combustion and the acoustic modes of the thrust chamber may induce a rapid growth of heat release rate fluctuations and cause dramatic damages to the propulsion system. During the second half of the twentieth century, most advances in the domain have been made by means of full-scale engine tests, as for example during the Apollo program for which more than 2,000 trial and error tests have been made to solve combustion stability problems encountered in the F-1 engine [START_REF] Oefelein | Comprehensive Review of Liquid-Propellant Combustion Instabilities in F-1 Engines[END_REF]. Despite decades of substantial research [START_REF] Harrje | Liquid propellant rocket combustion instability[END_REF][START_REF] Anderson | Liquid Rocket Engine Combustion Instability[END_REF][START_REF] Culick | Unsteady Motions in Combustion Chambers for Propulsion Systems[END_REF], origins of this unstable phenomenon are still not fully understood, essentially due to 1) the complexity of coupling mechanisms involved, 2) the limited instrumentation access to the thrust chamber inherent in the severe operating conditions and 3) the significant cost of full-scale engine tests. As a result, to improve knowledge on this critical issue whilst meeting with present-day industrial's need to reduce engine design costs, a three-pronged approach emerged, combining experiments on lab-scale set-ups, high-performance Computational Fluid Dynamics (CFD) and low-order modelling [START_REF] Poinsot | Prediction and control of combustion instabilities in real engines[END_REF].

Last block of this triptych is widely used in industry for engine development stage and aims at predicting the stability of the propulsion system at low cost by modelling the unsteady heat release rate as a function of acoustic perturbations. This has been initiated by [START_REF] Crocco | Theory of combustion instability in liquid propellant rocket motors[END_REF] with the sensitive time lag formalism, also designated as the n-τ model. This linear approach is essentially based on the idea that the unsteady heat release rate depends on two parameters, namely the interaction index n and the time lag τ . The former represents the sensitivity of propellant combustion to pressure fluctuations, while the latter corresponds to the time lag between the injection time of a propellant element and the time at which it burns. Despite a good estimation of the initial growth of the unsteady heat release rate, this approach is not adapted to predict the non-linear saturation phenomenon. Therefore, in order to predict limit cycle amplitudes and frequencies, which are key parameters in the design of a propulsion system, nonlinear approaches have since been proposed and can be divided into two groups. The first one relies on the idea that combustion is the main non-linear element governing the saturation phenomenon. Approaches of this group are then based on the describing function framework, in which the relative unsteady heat release rate is expressed as a function of the relative upstream velocity fluctuations, the frequency and the amplitude of oscillations [START_REF] Gelb | Multiple-input describing functions and nonlinear system design[END_REF][START_REF] Noiray | A unified framework for nonlinear combustion instability analysis based on the flame describing function[END_REF]. On the opposite, the keystone of the second group is the idea that limit cycle parameters are set by acoustic non-linearities and can be predicted by using the unified approach of [START_REF] Zinn | Application of the Galerkin Method in the Solution of Nonlinear Axial Combustion Instability Problems in Liquid Rockets[END_REF] (see also [START_REF] Culick | Nonlinear behavior of acoustic waves in combustion chambers-I[END_REF]) based on the expansion of the acoustic field on normal modes of the system together with the use of the linear n-τ model to describe the flame response. Despite their fast, low cost and easy-to-use aspects, these low-order approaches suffer from a low degree of generality and require previous experimental or CFD analysis to determine either the flame describing function or n and τ parameters, usually casedependent.

To address this lack of generality, it is necessary to get a better insight on driving mechanisms of combustion instabilities for the entire scope of engine operating conditions. In such engines, two regimes are identified, namely sub-and transcritical operating conditions. The former correspond to low thrust and transient operating phases and are characterized by low pressure levels in the combustion chamber, resulting in a liquidstate injection of oxygen. On the contrary, transcritical conditions are related to nominal regimes during which pressure levels are higher than the critical pressure of propellants, which induces an injection of oxygen at a quasi gaseous-state. It is also known that the damping of instabilities in LRE is usually higher for longitudinal acoustic modes than for transverse ones due to the nozzle damping effect [START_REF] Crocco | Theory of combustion instability in liquid propellant rocket motors[END_REF], yielding to more detrimental transverse instabilities. This is why most studies only focus on the effect of transverse acoustic perturbations on cryogenic-rocket-like jets. In the case of transcritical conditions, experiments have been performed on lab-scale configurations for non-reactive jets [START_REF] Davis | Measurements in an Acoustically Driven Coaxial Jet under Sub-, Near-, and Supercritical Conditions[END_REF]) in order to observe how coaxial two-phase flows respond to such transverse acoustic perturbations. The same kind of experiments have been carried out later for reactive jets (see [START_REF] Hardi | LOx Jet Atomization Under Transverse Acoustic Oscillations[END_REF] for example) to understand how these acoustics/two-phase flow coupling mechanisms may impact the flame behaviour. In addition, complementary investigations have been successfully carried out by means of CFD, either in gas/gas [START_REF] Selle | Prediction and Analysis of Combustion Instabilities in a Model Rocket Engine[END_REF] or transcritical conditions (Hakim et al. 2015a,b;[START_REF] Gonzalez-Flesca | Large eddy simulations of a transcritical round jet submitted to transverse acoustic modulation[END_REF][START_REF] Urbano | Exploration of combustion instability triggering using Large Eddy Simulation of a multiple injector liquid rocket engine[END_REF]. All these contributions have led to a better understanding of coupling mechanisms between acoustic perturbations and transcritical flames, which is essential to express general low-order models.

Regarding subcritical operating conditions, many differences exist in the coupling mechanisms due to the presence of a liquid phase. Indeed, physical mechanisms inherent to two-phase flows, such as atomization, fragmentation, droplet dispersion and vaporization, are sensitive to pressure and velocity fluctuations and may participate to the thermoacoustic coupling [START_REF] Oefelein | Comprehensive Review of Liquid-Propellant Combustion Instabilities in F-1 Engines[END_REF][START_REF] Culick | Unsteady Motions in Combustion Chambers for Propulsion Systems[END_REF]. This issue has been investigated experimentally in the past by submitting round or coaxial liquid jets to acoustic modulation and different response mechanisms have been observed. First, it has been shown by [START_REF] Hoover | Pressure Oscillation Effects on Jet Breakup[END_REF], and confirmed later by [START_REF] Chehroudi | Interaction of acoustic waves with a cryogenic nitrogen jet at sub-and supercritical pressures[END_REF], that acoustics may flatten the central liquid core along the acoustic propagation axis. When dealing with coaxial jets, [START_REF] Baillot | Behaviour of an airassisted jet submitted to a transverse high-frequency acoustic field[END_REF] explained this phenomenon by the suction effect induced by the heterogeneous distribution of acoustic radiation pressure around the liquid core. For acoustic modulation of sufficiently high amplitudes, the authors even demonstrated that this may lead to a drastic change in primary atomization processes. The second effect of acoustics that has been noticed on the liquid core is the reduction of its length, first in non-reactive conditions [START_REF] Davis | Measurements in an Acoustically Driven Coaxial Jet under Sub-, Near-, and Supercritical Conditions[END_REF][START_REF] Baillot | Behaviour of an airassisted jet submitted to a transverse high-frequency acoustic field[END_REF]). In the case of reactive two-phase flows, [START_REF] Méry | Experimental investigation of cryogenic flame dynamics under transverse acoustic modulations[END_REF] observed a direct impact of this phenomenon on the flame length, which may affect the stability of the system. The third effect of acoustics being noticed is the transverse oscillation of the liquid core at the acoustic modulation frequency. First observed for free liquid jets by [START_REF] Miesse | The Effect of Ambient Pressure Oscillations on the Disintegration and Dispersion of a Liquid Jet[END_REF] and [START_REF] Buffum | Response of Turbulent Jets to Transverse Acoustic Fields[END_REF], this observation has been confirmed later for coaxial two-phase flows by [START_REF] Davis | Measurements in an Acoustically Driven Coaxial Jet under Sub-, Near-, and Supercritical Conditions[END_REF]. These oscillations have been interpreted as being caused by the periodic fluctuations of the gas phase velocity around the liquid jet. In the case of reactive two-phase flows, the same observations have been made for the spray [START_REF] Méry | Experimental investigation of cryogenic flame dynamics under transverse acoustic modulations[END_REF], resulting in a periodic transverse movement of the flames that may drive HF combustion instabilities. Regarding the spray, another noticeable impact of acoustics have been observed by [START_REF] Baillot | Behaviour of an airassisted jet submitted to a transverse high-frequency acoustic field[END_REF] and consists in a transverse deviation of droplets due to acoustic radiation forces. Finally, [START_REF] Baillot | Behaviour of an airassisted jet submitted to a transverse high-frequency acoustic field[END_REF] and [START_REF] Ficuciello | Investigation of airassisted sprays submitted to high frequency transverse acoustic fields: Droplet clustering[END_REF] observed a change in the shape of the atomized jet, under specific operating conditions, due to a drastic expansion of the spray angle in the direction orthogonal to the acoustic propagation axis. Given all these observations, it appears that the combustion stability may thus be strongly affected by the impact of acoustics on the atomized liquid phase. However, experiments alone are not sufficient to reach a complete understanding of all coupling mechanisms. As proved in transcritical conditions, CFD analysis of reactive or non-reactive jets under acoustic modulation are adapted to bring complementary observations, but a significant lack of contribution is noticed in this domain. This might be due to numerical difficulties arising from the disparity in time and length scales of various phenomena related to the presence of the liquid phase, such as atomization, fragmentation, droplet dispersion and vaporization.

Numerical simulations of cryogenic flames under subcritical conditions has been performed in the past by [START_REF] Jay | Combined surface density concepts for dense spray combustion[END_REF] in a Reynolds average Navier-Stokes approach. Good agreements between the results and experimental data were obtained, but such approach is not suitable for the simulation of unsteady phenomena such as combustion instabilities. As regards unsteady numerical analysis, direct numerical simulations have provided promising results thanks to the coupling of interface tracking or capturing methods for the primary atomization process to Lagrangian methods for the description of the spray [START_REF] Herrmann | A parallel Eulerian interface tracking/Lagrangian point particle multiscale coupling procedure[END_REF][START_REF] Herrmann | On simulating primary atomization using the refined level set grid method[END_REF][START_REF] Kim | Numerical simulation of the breakup of a round liquid jet by a coaxial flow of gas with a subgrid Lagrangian breakup model[END_REF][START_REF] Zuzio | An improved multiscale Eulerian-Lagrangian method for simulation of atomization process[END_REF]. These approaches are nevertheless not yet adapted to industrial reactive configurations. To this end, Large-Eddy Simulations (LES) seem more adapted, but only few contributions can be found in the literature. Among them, LES based on the coupling of diffuse interface methods with either a Lagrangian [START_REF] Zuzio | An improved multiscale Eulerian-Lagrangian method for simulation of atomization process[END_REF] or an Eulerian (Le Touze 2015; [START_REF] Gaillard | Numerical simulation of cryogenic injection in rocket engine combustion chambers[END_REF] spray description have proved their capabilities. Indeed, the latter, in which the fully Eulerian formalism offers a natural and robust coupling to address primary atomization, has already been applied to simulate an isolated cryogenic LOX/GH 2 flame and reproduce a subcritical operating point of the Mascotte test bench [START_REF] Murrone | Atomization of cryogenic rocket engines coaxial injectors[END_REF][START_REF] Gaillard | Numerical simulation of cryogenic injection in rocket engine combustion chambers[END_REF]).

Therefore, this numerical strategy has been retained in the present contribution to advance the understanding of acoustic impact on two-phase flows. This coupling strategy is presented in section 2. Since it has initially been formulated without any acoustic modulation consideration, the present contribution aims not only to advance the understanding of the coupling mechanisms between transverse acoustic perturbations and atomized two-phase flows, but also to assess the behaviour of the previously mentioned numerical strategy once submitted to high-frequency acoustic modulation. To do so, the numerical investigation has to be performed on atomization regimes representative of LRE. The experimental configuration of [START_REF] Ficuciello | Investigation of airassisted sprays submitted to high frequency transverse acoustic fields: Droplet clustering[END_REF] is then simulated and the numerical specifications of this simulation are presented in section 3. To the authors' knowledge, this represents the first numerical study of the unsteady response of a nonreactive shear coaxial flow to transverse acoustic modulation under subcritical operating conditions. Some results without acoustic modulation are then presented in section 4 as an illustration of the numerical strategy aforementioned. Finally, qualitative comparisons of numerical results with experimental data are discussed in sections 5 to 7, as well as the behaviour of models under acoustic modulation. Note that quantitative comparisons with experiments do not fall within the scope of the present contribution and are planned for future work.

Coupling strategy for atomized two-phase flow simulations

General description

The numerical strategy retained in this work to simulate atomized two-phase flows is represented on figure 1. First, gas and liquid propellants are injected into the computational domain thanks to a Separated Phases Solver (SPS) dealing with the multispecies compressible Navier-Stokes equations within an LES formalism. In the context of two-phase flows, the latter is nothing else than a diffuse interface model based on the homogeneous relaxation model, also referred to as the 4-equation model [START_REF] Downar-Zapolski | The non-equilibrium relaxation model for one-dimensional flashing liquid flow[END_REF][START_REF] Goncalvès | Modeling for non isothermal cavitation using 4-equation models[END_REF], with the difference that no mass transfer between the liquid and gas phases is considered. At the coaxial injector outlet, the dense liquid phase is submitted to strong shear stresses due to the high surrounding gas injection velocity. Hence, this induces the atomization of the dense phase and produces multiple liquid structures of different sizes and shapes. In theory, it would be possible to pursue the description of these structures with the SPS up to the smallest droplets, but this would require extremely refined meshes unsuitable for practical applications. Therefore, the choice is made in this work to fix the mesh to a reasonable refinement capable of describing large liquid structures and to use a statistical approach to deal with smaller ones that constitute the spray. A mass transfer from the SPS toward a kinetic-based Eulerian solver is then performed once the mesh becomes too coarse. Such Dispersed Phase Solver (DPS) is far less restrictive in terms of mesh refinement since it provides a subgrid modelling of the spray. In addition, a reverse transfer from the DPS to the liquid phase of the SPS accounts for a possible impingement of droplets onto the dense liquid structures. This coupling strategy has been developed by [START_REF] Gaillard | Numerical simulation of cryogenic injection in rocket engine combustion chambers[END_REF] and is implemented in Cedre, Onera's multi-physics simulation code [START_REF] Refloch | CEDRE Software[END_REF], with Cedre solvers Charme and Spiree as the SPS and the DPS respectively.

Governing equations

For the SPS, the system of governing equations can be written in vector form as:

∂ t Q + ∇ • (F -ϕ D -ϕ C ) = S G + S L + S F , (2.1)
where ∂ t is the time partial derivative and ∇ • () is the divergence operator. Q = (ρY g ρY l ρu ρe t ) t is the vector of conservative variables, where Y g and Y l stand for the mass fraction of the gas and the liquid phases, while ρ, u and e t correspond to the mixture density, velocity and specific total energy respectively. Note that the density of the mixture reads ρ = ρY g +ρY l = α g ρ g +α l ρ l , with ρ g and ρ l the densities of pure phases expressed through dedicated equations of state, while α g = ρY g /ρ g and α l = ρY l /ρ l are the phase volume fractions. Since both fluids are injected at ambient temperature in this contribution, no thermal expansion is considered for the liquid phase. Considering only compressibility, the equation of state of the liquid thus reads

ρ l (p) = ρ 0 [1 + β 0 (p -p 0 )],
where β 0 is the isothermal compressibility of the liquid at the reference state (ρ 0 , p 0 ) and p is the local pressure of the mixture. Regarding the gaseous phase, the ideal gas law is retained. The thermodynamic closure of the system is then performed thanks to a thermal and mechanical equilibrium assumption [START_REF] Goncalvès | Modeling for non isothermal cavitation using 4-equation models[END_REF][START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF].

In system (2.1), F = Q ⊗ u + p(0 0 I u) t is the vector of convective fluxes, with I the identity tensor and ⊗ the tensor product. The vector of laminar and turbulent diffusive fluxes reads ϕ D = (0 0 τ τ • u + λ∇T ) t , where λ and ∇T are the thermal conductivity and the temperature gradient of the mixture respectively. Note that classical Fourier's law is used to describe the laminar diffusion of heat even if it is not significant in this study, and that no mass diffusion is considered as the gas and liquid species belong to two different non-miscible phases. For the turbulence modelling, the Boussinesq's concept of eddy viscosity is adopted. The viscous stress tensor thus reads τ = 2(µ + µ t )D, where µ is the laminar dynamic viscosity of the mixture, µ t is the turbulent dynamic viscosity and D is the deviator of the resolved symmetric strain rate tensor expressed as:

D = 1 2 ∇ ⊗ u + (∇ ⊗ u) t - 1 3 (∇ • u)I, (2.2)
with (∇⊗u) t the velocity gradient tensor. The modelling of the turbulent viscosity inside the gaseous phase is performed thanks to the Wall-Damping version of the [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. the basic experiment*[END_REF] model expressed as µ t = ρ min(C s ∆, κd w ) 2 f t , where C s = 0.18, as initially proposed by [START_REF] Lilly | The representation of small-scale turbulence in numerical simulation experiments[END_REF], κ = 0.41 is the von Kármán constant, ∆ is the size of the local control volume and d w is the distance to the nearest wall. The turbulent frequency is expressed as f t = (2D : D) 1/2 , where ():() is the double dot product of two tensors.

In addition, the capillary fluxes read

ϕ C = (0 0 τ C τ C • u) t
, where τ C is the capillary stress tensor of the Continuum Surface Stress method [START_REF] Lafaurie | Modelling Merging and Fragmentation in Multiphase Flows with SURFER[END_REF]. The latter is a conservative alternative to the Continuum Surface Force approach of [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF], which provides a volume reformulation of surface tension forces. The capillary stress tensor is written as:

τ C = σ ∇α l ∇α l ∇α l ⊗ ∇α l ∇α l -I , (2.3)
with σ the surface tension of the liquid-gas interface. Finally, S G is the gravity source term written as S G = (0 0 ρg ρg • u) t , with g the gravitational acceleration , while S L and S F are source terms referring to the coupling between the SPS and the DPS. These coupling source terms will be discussed later.

Regarding the dispersed phase, its description at the highest level of accuracy is provided by the Williams-Boltzmann kinetic equation (WBKE) [START_REF] Williams | Spray Combustion and Atomization[END_REF], which describes at a mesoscopic level the evolution of the particle Number Density Function (NDF) in the phase space. However, the direct resolution of the WBKE being out of reach for practical applications (Murrone & Villedieu 2011), Eulerian methods rather solve transport equations derived from the WBKE for some particular moments of the NDF (see [START_REF] Emre | Eulerian moment methods for automotive sprays[END_REF] for a comprehensive review on Eulerian moment methods). An important part of the modelling is then how to deal with the discretization of the particle size distribution. Among the several approaches developed for this purpose in the framework of Eulerian moment methods, the sectional method [START_REF] Greenberg | On the origins of spray sectional conservation equations[END_REF][START_REF] Laurent | Multi-fluid modelling of laminar polydisperse spray flames: origin, assumptions and comparison of sectional and sampling methods[END_REF] has been retained in the strategy of [START_REF] Gaillard | Numerical simulation of cryogenic injection in rocket engine combustion chambers[END_REF] for its natural ability to account for phenomena that generate polydispersion, such as breakup. In the sectional approach, the size space is discretized into N s contiguous sections, within which the shape of the particle size distribution has to be postulated. To this end, the affine-TSM (Two Size Moment) reconstruction of [START_REF] Sibra | Simulation of reactive polydisperse sprays strongly coupled to unsteady flows in solid rocket motors: Efficient strategy using Eulerian Multi-Fluid methods[END_REF] was retained by [START_REF] Gaillard | Numerical simulation of cryogenic injection in rocket engine combustion chambers[END_REF], as it seems to offer the best compromise between accuracy and computational cost. Finally, N s systems of Eulerian equations are thus solved by the DPS, one for each section of the particle size distribution. They read in vector form (k denoting the number of a particular section):

∂ t q k + ∇ • (q k ⊗ v k ) = s k G + s k L + s k F + s k B , 1 k N s , (2.4)
with q k = (ρ k ρ k v k n k ) t the vector of conservative variables consisting of the bulk mass density, momentum and number density, with v k the mean velocity vector of droplets in the section. Note that the configuration studied in this paper is quasi isothermal.

Indeed, the gas and liquid phases are injected at the same temperature and temperature fluctuations due to acoustic perturbations are considered negligible. The temperature of droplets θ can then be considered constant and equal to the injection temperature of the liquid. As a result, no equation on the internal energy needs to be solved in the DPS. As for the SPS, the gravitational force is considered through a source term s k G = (0 ρ k g 0) t , while s k L and s k F are source terms referring to the coupling between the SPS and the DPS. The last term s k B stands for the exchange source term between the spray sections due to the fragmentation process.

Source terms

2.3.1. Coupling between the liquid phase and the spray (S L and s k L ) According to [START_REF] Gaillard | Numerical simulation of cryogenic injection in rocket engine combustion chambers[END_REF], primary atomization is modelled by a liquid mass transfer from the SPS to the DPS that reads:

S a = ρY l f a λ a , (2.5)
where f a is the characteristic frequency of the primary atomization process and λ a is an activation function. The atomization process downstream coaxial injectors mainly depends on the local velocity difference between the gas and the liquid. In the context of homogeneous mixture models, where the whole mixture is described by only one mean velocity u, the local velocity gradient can be used as an estimate. A simple choice is then to express the frequency of the primary atomization process as f a = f t . The shape of λ a in equation (2.5) is chosen in an empirical way to estimate the position of the interface but also to ensure numerical stability. This activation function thus reads λ a = 1-tanh(4Y 2 l ) (see figure 2). Indeed, it is important that the transfer operates on the gaseous side of the diffuse interface where the mass of liquid in the SPS is negligible. Otherwise, the low compressibility of the liquid would induce significant pressure oscillations because the volume of droplets is not considered by the SPS. Finally, to perform the liquid mass transfer toward the DPS, the initial mean diameter D a and velocity v a of the droplets thus created have to be defined. In absence of any information about the local velocity difference between both phases, a pragmatic approach is adopted in which these parameters are chosen constant and uniform and depend on the configuration to be treated (see section 3.4 for operating conditions simulated in this work). As a consequence, the mass transfer from the SPS is only operated toward one section of the DPS, namely k = k a .

In addition, a reverse liquid mass transfer from the DPS to the SPS accounts for the possible impingement of droplets onto the dense liquid core and ligaments. This mass transfer is expressed for each section k as:

s k c = ρ k f c λ c , (2.6)
with f c = 1/∆t the characteristic frequency, ∆t the integration time step, and λ c = 1tanh(-2 log α l ) the activation function empirically defined to meet with numerical issues (see figure 2). Note that f c is not based on a physical time scale but on a numerical one because impingement has to be treated as quasi instantaneous. According to the previous relations, the source terms S L and s k L can be expressed as:

S L =               0 -S a + Ns k=1 s k c -S a v a + Ns k=1 s k c v k -S a e 0 + 1 2 v a • v a + Ns k=1 s k c e 0 + 1 2 v k • v k               , (2.7)
and

s k L =         -s k c -s k c v k - s k c m k         k = k a , s ka L =       S a -s ka c S a v a -s ka c v ka S a m(D a ) - s ka c m ka       , (2.8)
with m k = ρ 0 πD k 3 /6 the mass of one droplet of diameter D k . Since the configuration studied in this paper is quasi-isothermal and pressure is not considered in the DPS, the density and specific internal energy of the pure liquid in this solver can be assumed constant and respectively equal to reference state values ρ 0 and e 0 . In the same way, m(D a ) = ρ 0 πD a 3 /6 is the mass of one droplet of diameter D a .

2.3.2.

Coupling between the gas phase and the spray (S F and s k F ) In the context of LRE applications, the density of the liquid phase is very high compared to the density of the gaseous phase and droplets may experience strong acoustic perturbations. Hence, in addition to gravity, the only forces acting on the dispersed phase are a priori drag and acoustic radiation forces.

The drag force experienced by a spherical isolated droplet of diameter D and velocity v in a medium of density ρ g and velocity u can be calculated from the standard expression:

F D = 1 8 πD 2 ρ g C D u -v (u -v).
(2.9)

The drag coefficient C D is estimated thanks to the correlation of [START_REF] Schiller | A drag coefficient correlation[END_REF] and depends on the droplet Reynolds number Re p defined as:

Re p = ρ g D u -v µ g , (2.10)
with µ g the dynamic viscosity of the gaseous phase. Note that in the sectional approach, the drag force needs to be integrated over the size interval of each section k to get the total drag force of the section F k D . In addition, liquid sprays submitted to a transverse acoustic modulation experience an acoustic radiation force. Indeed, as mentioned by [START_REF] King | On the Acoustic Radiation Pressure on Spheres[END_REF], any object located in an acoustic field acts as an obstacle for the acoustic wave, thus inducing a wave reflection at its surface. The acoustic radiation pressure distribution around the object arising from this mechanism may lead to a resulting acoustic radiation force F R . Regarding the SPS, the acoustic radiation pressure distribution around the liquid structures is naturally rendered thanks to Navier-Stokes equations. However, the coupling between the SPS and the DPS only intervene through source terms. Consequently, the ambient fluid does not take into account the volume occupied by the dilute spray and the local acoustic field around each spherical droplet is not simulated. Therefore, the acoustic radiation force has to be modelled in the DPS if one aims at taking into account its effect on the spray.

It is possible to evaluate a priori the deviation of droplets by F R in the configuration presented in this paper. To simplify, the acoustic field is supposed to be mono-harmonic. The time averaging of the Newton's second law applied to one isolated droplet of constant mass m projected in the acoustic modulation direction X leads to:

ma X = F R + F D, X , (2.11)
where a X is the droplet acceleration in this direction, while F D, X is the X-component of the drag force. Assuming that every droplet undergoes periodic oscillations around the injection axis under the effect of the acoustic velocity, and that the mean deviation due to the acoustic radiation force is small, F D, X can be considered null. Regarding F R , [START_REF] King | On the Acoustic Radiation Pressure on Spheres[END_REF] formulated the time averaged acoustic radiation force experienced by a small spherical droplet located in an intensity anti-node (IAN) of a standing acoustic field as:

F R = 3 4 m p ac ρ g c 2 k η 1 + 2 3 (1 -η) 2 + η , (2.12)
with p ac and f ac the amplitude and frequency of pressure fluctuations respectively, c the speed of sound in the medium, k = 2πf ac /c the acoustic wave number and η the density ratio between the gas and the liquid. Hence, for the operating conditions presented in this contribution, the mean deviation ∆X of an isolated droplet injected with a null transverse velocity in an IAN of a standing acoustic field can be deduced by equations 2.11 and 2.12. This leads to ∆X = F R /(2m) t 2 , where t is the time elapsed since the injection of the droplet. According to the axial velocity of droplets in the present configuration once accelerated by the gas, each one needs 6 ms in average to span the study area. Therefore, the spray is deviated up to 40 µm at a maximum, which can be neglected compared to the transverse dimension of the spray. As a result, the acoustic radiation force on droplets is not modelled in the DPS for the present contribution. Note that the evaluation of the mean deviation of droplets by F R is however case dependent. The modelling of acoustic radiation forces on droplets may then still be necessary when dealing with different operating conditions. As a result, it has been implemented in our numerical tool in view of future work. As a result, the source terms S F and s k F referring to the coupling mechanisms between the gas phase and the spray can be expressed as:

S F =              0 0 - Ns k=1 F k D - Ns k=1 F k D • v k              , s k F =        0 F k D 0        . (2.13) 2.3.3.
Exchanges between the sections of the spray (s k B ) The fragmentation of droplets due to the velocity differential with the gaseous phase is modelled at the scale of a single isolated droplet, and then formulated in the sectional formalism under the form of exchanges between the different sections. Indeed, droplets belonging to any section k are prone to break-up into smaller droplets whose size fall into the range of any inferior section. To estimate the physical properties of droplets after fragmentation, three key elements have to be evaluated: the fragmentation frequency ν fr of a single mother droplet of diameter D m and velocity v m , as well as D fr and v fr , the mean diameter and velocity of daughter droplets thus created.

Following the model proposed by [START_REF] Pilch | Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop[END_REF], the fragmentation frequency ν fr depends on different break-up regimes, and is expressed as a function of the local Weber number We = ρ g uv m 2 D m /σ. The model proposed by [START_REF] Wert | A rationally-based correlation of mean fragment size for drop secondary breakup[END_REF] is then used to provide D fr as the Sauter mean diameter of the daughter droplets, depending on D m , We and ν fr . Finally, a zero-dissipation fragmentation is assumed, so that the velocity of daughter droplets v fr equals to the velocity v m of the mother droplet.

Once these three elements are evaluated, the fragmentation operator may be formulated and integrated over the size interval of each section k. For the sake of readability, the resulting expression of the fragmentation source term s k B is not presented in this paper, but more details are given by [START_REF] Dufour | Étude d'un modèle de fragmentation secondaire pour les brouillards de gouttelettes[END_REF].

Geometry and numerical specifications

Computational domain and mesh

To investigate interaction mechanisms between acoustics and two-phase flows in injection conditions representative of those in LRE, the experimental configuration of [START_REF] Ficuciello | Investigation of airassisted sprays submitted to high frequency transverse acoustic fields: Droplet clustering[END_REF], from Coria laboratory, is simulated. The experimental test-rig, represented on figure 3, consists in a cavity equipped with four speakers at the bottom and a coaxial injector at the top. In order to limit numerical costs, the numerical domain is reduced to a square section box of dimensions 64.3 D l ×35.7 D l ×35.7 D l with a 1.79 D lhigh coaxial injector at the top, D l being the injection diameter of the central liquid flow at the exit of the injector. Figure 3 reveals that in this domain, represented as a grey rectangle, the experimental acoustic field can be considered as a plane standing acoustic wave. Note that in this configuration, the position of the injector can be set along the acoustic axis, mentioned as a.a. in the remainder of this paper. According to [START_REF] Baillot | Behaviour of an airassisted jet submitted to a transverse high-frequency acoustic field[END_REF], most acoustic effects on the liquid jet are maximum at velocity anti-nodes (VAN), such as the flattening of the liquid core, the modification of the spray angle and the transverse oscillation of droplets. However, at such location on the a.a., the deviation of the liquid jet by acoustic radiation forces is null. Therefore, a compromise is made in this paper in order to maximise the number of acoustic effects observed in the numerical simulations. The injector is thus located at an IAN of the second transverse mode of the cavity (2T), where all previously mentioned acoustic effects can intervene.

The computational mesh represented on figure 4 is made up with 28 million tetrahedra and prisms. To capture dense liquid structures and external vortices around the liquid jet, the mesh is refined near the injector in a 10.7 D l -high cylinder of diameter 2.5 D l with a constant element size of 1.0 δ lip , where δ lip is the injector lip thickness. A second 10.7 D l -high cylindrical box, of greater diameter 8.9 D l , is filled with 3.5 δ lip -size elements on average in order to assure the resolution of at least 80 % of the turbulent kinetic energy in the gas [START_REF] Pope | Turbulent flows[END_REF]. Based of the same criterion, the injector is filled with 0.06 D l -size tetrahedra and four layers of 0.03 D l -high prisms inside the central tube. Considering e as the thickness of the external coaxial channel, the numerical mesh is refined in the annular tube down to 0.1 e for tetrahedra and 0.05 e for prisms. This results in two tetrahedra within the injector lip thickness. In the rest of the domain, the mesh is coarsened up to an element size of 0.89 D l to minimize the dissipation of acoustic waves. As recommended by [START_REF] Troyes | Study of Impinging Supersonic Jet Noise with Aerodynamics and Acoustics Numerical Simulations[END_REF], this corresponds to more than 20 points for the discretization of the smallest acoustic wavelength to be simulated.

Operating and boundary conditions

The operating conditions investigated in the present contribution have been selected in collaboration with Coria in order to meet with both experimental and numerical requirements. As a result, the liquid mass flow rate has been increased compared to operating conditions studied by [START_REF] Ficuciello | Investigation of airassisted sprays submitted to high frequency transverse acoustic fields: Droplet clustering[END_REF] in order to limit computational costs. Water and air are both injected at atmospheric pressure and a temperature of 300 K. This gives rise to a momentum flux ratio J = 4 and a gaseous Weber number We g = 490, according to the following expressions:

J = ρ g U 2 g ρ l U 2 l ; We g = ρ g U 2 g D l σ , (3.1)
where U g and U l are the air and water velocities at the exit of the coaxial injector.

According to [START_REF] Lasheras | Liquid Jet Instability and Atomization in a Coaxial Gas Stream[END_REF], this corresponds to a fibre-type atomization regime representative of LRE operating under subcritical conditions. Constant and homogeneous injection velocities and temperatures are imposed at the inlet boundaries of the injector. Adiabatic no-slip conditions are prescribed at injector walls and lip, as well as at the top-wall of the cavity.

Acoustic forcing

Two cases are studied and compared to each other in experiments and numerical simulations: one with an acoustic modulation, namely the modulated case (MC), and the other one without any, namely the non modulated case (NMC).

Regarding the MC, a maximum sound pressure level of 166 dB has been prescribed experimentally by Coria to maximize the response from the two-phase flow to acoustic perturbations. The raw experimental pressure signal p measured at a Pressure Anti-Node (PAN) of the 2T mode near the top wall of the cavity is plotted on figure 5 i over one acoustic period. From this figure, it is clear that the acoustic field contains multiple harmonics causing the deformation of the pressure signal. More precisely, signal processing showed that 98.3 % of the acoustic energy of the raw experimental signal is comprised into the frequency range [0 Hz ; 3,000 Hz], and more precisely in the 2T, 4T and 6T modes of the cavity. Therefore, only these three modes are considered for the numerical reproduction of the experimental acoustic field. To this end, pressure fluctuations p (t) are prescribed on the lateral boundaries of the domain orthogonal to the a.a., following the form:

p (t) = 3 i=1 p ac, i sin (2πf ac, i t + ψ ac, i ), (3.2)
where f ac, i , p ac, i and ψ ac, i are the frequency, amplitude and phase of each harmonic i.

These parameters are evaluated thanks to a discrete Fourier transform on the experimental signal of figure 5 and are given in table 1. As concerns the lateral and bottom numerical boundaries parallel to the a.a., these correspond to free-flowing limits in the experiment.

It is thus chosen to prescribe on these boundaries an unsteady and heterogeneous pressure profile that complies with pressure fluctuations in the cavity along the a.a.. Note that it is necessary to evacuate the outgoing waves to establish a standing acoustic wave in the domain. To do so, non-reflecting conditions based on a relaxation approach are used for acoustic forcing and inlet boundaries.

In the NMC, the cavity boundaries orthogonal to the a.a. are replaced with adiabatic no-slip conditions in order to represent the lateral panels in the experiment, while the other acoustic forcing boundaries are set to constant pressure conditions.

Modelling parameters, numerical schemes and simulation strategy

As mentioned in section 2.3.1, the modelling of the primary atomization requires to postulate properties of the created droplets. To do so, the physical analysis of [START_REF] Marmottant | On spray formation[END_REF], linking the diameter and the velocity of created droplets to the injection velocities of gas and liquid flows, is used. For the present operating conditions, this approach gives rise to the following parameters: D a = 200 µm and v a = 3.5 m.s -1 . Note that coalescence between droplets is not taken into account in the present contribution. As a result, droplets issued from primary atomization are created into the spray section containing the biggest droplets, i.e. k a = N s . The spray size distribution is divided into N s = 3 sections for a good compromise between accuracy and numerical cost. For the sake of simplicity, these three sections are mentioned as small, medium and big sections in the present paper. Their upper bounds in terms of droplet diameter are respectively 50 µm, 150 µm and 300 µm.

As concerns numerical methods, both the SPS and the DPS are based on a cellcentred finite volume approach on general unstructured meshes. Numerical fluxes are computed by means of upwind numerical schemes in both solvers together with secondorder reconstructed values at the faces centroids thanks to the multislope MUSCL The first step of the numerical simulation has been to obtain an established solution of the two-phase flow without any acoustic perturbation. From this state, the computation is then carried out over a physical time equivalent to 18 and 23 acoustic periods for the NMC and the MC respectively. This corresponds to more than 20 and 0.4 convective times of the gaseous and liquid phases of the SPS respectively along the refined 10.7 D lhigh cylinder of figure 4. Even if this simulation time is not long enough to achieve average convergence for the liquid phase, it permits to establish the two-phase flow under acoustic modulation. Hence, a complete analysis of the transitory behaviour of the liquid from the beginning of the acoustic forcing is performed. To simplify the analysis of the numerical results, the time coordinate is referred to in terms of T ac , the number of acoustic periods of the cavity 2T mode, even in the NMC where no acoustic forcing is set.

The numerical simulation of the two cases required a total of 4.3 MhCPU on Intel Haswell processors, with typical runs involving 1,000 MPI processes.

Two-phase flow simulation without any modulation

Instantaneous shape of the two-phase flow

In order to get a first overview of the numerical strategy and thus facilitate upcoming interpretations, a qualitative comparison of experimental and numerical results without acoustic modulation is performed. Figure 6 represents snapshots of the two-phase flow near the exit of the injector for experiments and numerical simulations. Experimental images (see figure 6a) are acquired thanks to high-speed back-light visualizations at a frequency of 6.2 kHz. Regarding numerical snapshots (see figure 6b), the visualization of the dense liquid phase solved by the SPS and its diffuse interface model requires to fix a critical value of liquid volume fraction. Indeed, the interface between the gas (α l = 0) and the liquid (α l = 1) is spatially diffused over the mesh. Dense liquid structures are thus represented on figure 6b with a black 3D iso-surface of α l = 0.99, which will be denoted α l, 0.99 in the sequel. As regards the spray, the presence of droplets is represented thanks to a longitudinal slice of spray volume fraction α tot = k α k in grey scale.

On both images, the central liquid core can be observed as a long black structure that extends from the exit of the injector. Both snapshots reveal hydrodynamic instabilities growing at the surface of the liquid core under the action of shear stresses due to the coaxial gas flow. These instabilities then result in the detachment of large liquid structures at the tip of the central core which propagate along the injection axis. As regards the spray, figure 6a reveals a peeling of the liquid core near the exit of the injector by the gaseous flow which results in the creation of small droplets. This early atomization process seems to be correctly retrieved in the numerical simulation, with the presence of droplets in the mixing layer between the gas and liquid flows right at the exit plane. The droplets thus created are convected downward by the gaseous stream and expand radially to form a conical spray, both on experimental and numerical images. From these observations, it appears that the coupling strategy of section 2 is able to retrieve the macroscopic topology of the liquid flow.

Mean shape of the central liquid core

More quantitative comparisons are necessary to assess the numerical reproduction of the two-phase flow without acoustic modulation. Unfortunately, the presence of the spray in experimental back-light images impedes to perform any measurement of the liquid core length for the operating conditions simulated in this paper. However, experimental correlations obtained in the past can be used as reference. [START_REF] Villermaux | Mixing and Spray Formation in Coaxial Jets[END_REF] conducted experiments to characterized the time averaged probability of presence of the central liquid phase at each location downstream of the injector outlet thanks to a laser-induced fluorescence technique. A mean dye concentration C was thus measured along the jet axis and normalized by the initial mean concentration C 0 at the inner nozzle. For the dilution length L x , defined as the distance from the injector outlet to reach a given dilution level C/C 0 = x, the authors found the following correlations:

L 1.0 D l = 6 J -1 2 , L 0.5 D l = 12 J -1 2 . (4.1)
To compare with these correlations, the dilution level along the jet axis can be assessed in the NMC thanks to the time averaged liquid volume fraction α l . A planar cut of α l containing the injection axis is plotted on figure 7. The time average operation is performed over 18 acoustic periods. From this two-dimensional field, it is possible to extract the length of iso-lines α l, 1.0 and α l, 0.5 . For the present operating conditions where J = 4, the numerical results L 1.0 = 3.51 and L 0.5 = 5.95 are in good agreement with experimental correlations predicting L 1.0 = 3.0 and L 0.5 = 6.0. This contributes to the validation of the coupling strategy of section 2 for the simulation of two-phase flows without acoustic perturbation. The next step of this work is to evaluate the ability of the code to retrieve the effects of the acoustic modulation on the jet. First, the flattening of the liquid core by acoustic radiation pressure and the induced decrease of its length are analysed in section 5. Then, the impact of the acoustic modulation on the primary atomization process is studied in section 6. Finally, effects of acoustic perturbations on the shape and dynamics of the spray are examined in section 7.

5. Response of the liquid core to the acoustic forcing

Flattening of the liquid core

As mentioned above, the iso-surface α l, 0.99 can be used to represent liquid structures solved by the SPS. The structures thus obtained after 18 acoustic periods are represented Iso-lines of α l, 1.0 and α l, 0.5 . on figures 8a and 8b for the NMC and the MC respectively. In addition, transversal cuts of liquid volume fraction are displayed in grey scale at different distances from the injection plane in order to better visualise the evolution of the liquid core shape. By comparing these two figures it is clear that, from a certain distance downstream the injection plane, acoustic perturbations tend to flatten the liquid core in the X direction (i.e. the a.a.) and stretch it in the spanwise Y direction. This phenomenon is difficult to observe on experimental images for the present operating conditions due to the presence of multiple ligaments and droplets hiding the liquid core. However, it has already been observed by [START_REF] Ficuciello | Investigation of airassisted sprays submitted to high frequency transverse acoustic fields: Droplet clustering[END_REF] for other operating conditions with lower momentum flux ratio and gaseous Weber number. Back then, the authors explained this deformation by considering the acoustic radiation pressure distribution p rad (θ) around the liquid jet. Based on the work of [START_REF] King | On the Acoustic Radiation Pressure on Spheres[END_REF], [START_REF] Zhuk | Radiation force acting on a cylindrical particle in a sound field[END_REF] and [START_REF] Wu | Acoustic radiation pressure on a rigid cylinder: An analytical theory and experiments[END_REF][START_REF] Ficuciello | Analysis of high frequency/high amplitude acoustic field effects on coaxial injection: application to liquid rocket engines[END_REF] expressed the dimensionless acoustic radiation pressure distribution around the liquid in the case of a small cylindrical-like liquid core located in a mono-harmonic standing acoustic wave as:

p * rad (θ) = p rad (θ) P = cos 2 (kh) 4
(5.1)

+ sin 2 (kh) (1 + η) 2 η 2 cos 2 θ + (2η -1) sin 2 θ - kR sin(2kh) 2(1 + η) sin θ sin(2θ)(1 -η) + cos θ .
In this expression, R is the cylinder radius, h is the distance along the a.a. between the axis of the cylinder and any PAN and θ is the angular coordinate around the cylinder. Moreover, the acoustic radiation pressure is normalized by P = p 2 ac /(ρ g c 2 ), p ac being the acoustic pressure amplitude, and the small cylinder approximation is defined by kR 1. The theoretical distribution of p * rad (θ) is plotted on figure 9 for the particular case addressed in this paper (i.e. kh = π/4 for an injection at an IAN and kR = 4.9 × 10 -2 ). Since P depends on p 2 ac , the effect of harmonics can be neglected for simplification reasons. The acoustic radiation pressure is then only calculated for the 2T mode of the cavity. It can be seen on figure 9 that the acoustic radiation

f =    a b -1 if a < b 1 - b a otherwise. ( 5

Decrease of the intact liquid core length

In addition to the flattening of the liquid core, figure 8 seems to reveal a drastic decrease of its length once submitted to the acoustic modulation. The time evolution of the intact liquid core length L lc normalized by D l is plotted on figure 12 for both the NMC and the MC, with different α l threshold values. If we first consider the NMC, drops of L lc are noticed at low frequency, whatever the retained α l threshold value. This phenomenon has already been observed experimentally by [START_REF] Hardi | LOx Jet Atomization Under Transverse Acoustic Oscillations[END_REF] in transcritical operating conditions. For both sub-and transcritical conditions, this growthand-detachment behaviour results from the natural disintegration of the dense core into large dense structures, even if the two cases present strong differences of surface tension and density ratio between the two injected fluids. Note that L lc strongly depends on the value of α l . Indeed, since the interface is diffused in the SPS, a liquid fragment that appears to be just torn from the liquid core can still be considered as part of the jet for a smaller value of α l . This consequently complicates any quantitative comparison between instantaneous experimental and numerical measurements.

Regarding the response of the liquid core in the MC, the evolution of L lc can be split in two : a transitory and an established phase. The transitory phase starts at the beginning of the plot, when the acoustic forcing is set, and ends after approximately 10 to 13 acoustic periods, depending on the value of α l . During this time interval, the liquid core length keeps a similar evolution to that of the NMC, with only few differences at the end. After that, the length suddenly drops whatever the value of α l . From this point, the trend of L lc seems to be established, with a growth-and-detachment behaviour similar to that of the NMC. If we consider the threshold value α l, 0.99 , the drop of L lc occurs after 12 acoustic periods and reaches about 40 % of the initial length. Such behaviour has been noticed both in trans-and subcritical operating conditions by [START_REF] Hardi | LOx Jet Atomization Under Transverse Acoustic Oscillations[END_REF] and [START_REF] Baillot | Behaviour of an airassisted jet submitted to a transverse high-frequency acoustic field[END_REF] respectively. In the former contribution, this has been attributed to a transverse shedding mechanism of the dense phase due to gas oscillations. α l, 0.9 , α l, 0.99 α l, 1.0 (black: NMC ; red: MC).

For sufficiently high acoustic amplitudes, large pockets of oxygen were stripped away from the dense core in the direction of the a.a. by the transverse acoustic velocity, without giving the opportunity to large structures to detach from the tip of the dense core. However, in subcritical conditions, no such stripping mechanism has been observed in experiments and numerical simulations because of the higher density ratio between the gas and the liquid. To explain the decrease of L lc under acoustic perturbations, we thus refer to figure 13. It represents the liquid core for the NMC and the MC before (see figure 13a) and after (see figure 13b) a drastic decrease of the liquid core length in the MC (i.e. around t = 12 T ac ). In addition, transversal cuts of the liquid interface are plotted at Z = 4.5 D l , corresponding to the distance at which operates the cut of the liquid core. It can be noticed on figure 13a that the flattening effect analysed in the previous section tends to pinch the interface of the liquid core. This results in an early detachment of liquid structures at its tip on figure 13b. In the present case, the amount of time for this flattening mechanism to cause the cut of the liquid core actually corresponds to the time interval of the transitory phase, i.e. approximately 10 acoustic periods. After this amount of time, the thickness of the liquid core reaches a critical value causing the core break-up. Once it has been cut, the length of the liquid core keeps a lower value than in the NMC as long as acoustic perturbations are present.

The diffuse interface method in the SPS is thus well adapted to reproduce the flattening of the liquid core and the decrease of its length observed experimentally due to acoustic perturbations. In addition, the analysis of the transitory behaviour of the liquid flow gives the opportunity to advance understanding of the link between acoustic radiation pressure and the cut of the liquid core. This is of first importance since this may have a drastic impact on the flame length in LRE and thus participate to the thermoacoustic coupling.

6. Impact of the modulation on the primary atomization process

Response of the Separated Phases Solver

The main goal of coaxial injectors used in LRE is to assist the atomization of the liquid phase and maximize, for a given mass of propellant, the amount of interface shared with the gaseous phase. This is essential to promote the evaporation of the liquid and consequently the intensity of the combustion. Therefore, it is of great interest to investigate the evolution of the liquid interface area once submitted to acoustic perturbations. To characterize the interface, the interfacial area density Σ, defined as the amount of liquid surface per unit of volume, can be used. In the scope of diffuse interface models such as the one used in the SPS, [START_REF] Sun | Diffuse interface modeling of two-phase flows based on averaging: mass and momentum equations[END_REF] specified that the resolved interfacial area density may be locally defined as Σ SPS = ∇α l in the control volume V i . It is then natural to define the total interfacial area density in a given volume of interest Ω as Σ SPS, Ω = Vi∈Ω ∇α l V i /Ω. Based on this definition, the time evolution of Σ SPS, Ω is represented on figure 14 for both the NMC and the MC, with Ω corresponding to the refined zone of diameter 8.9 D l presented on figure 4b plus the volume of the injector. Since both cases start from an established solution, the value of Σ SPS, Ω at the beginning of the plot is not null. In the NMC, it can bee seen that Σ SPS, Ω does not depart much from its time averaged value. On the contrary, in the MC, Σ SPS, Ω starts to progressively increase after approximately one acoustic period, i.e. right after the establishment of the standing acoustic wave. Then, after approximately 13.5 acoustic periods, Σ SPS, Ω starts to decrease slowly. The increase of Σ SPS, Ω can be easily explained by the effect of the acoustic radiation pressure on the liquid. Indeed, the flattening of the liquid core results in the expansion of the liquid interface in the direction orthogonal to the a.a. and hence to the increase of the liquid core surface. In addition, after approximately 12 acoustic periods, the pinching of the liquid by acoustics promotes the break-up of the liquid core into smaller structures at its tip, hence increasing the amount of interface between both fluids. During this process, the liquid structures thus created themselves experience acoustic radiation pressure, which tends to promote their break-up as well. As a result, the amount of interface is supposed to keep increasing and then stabilize once the size of liquid fragments reaches an established value due to the equilibrium between acoustic radiation pressure and surface tension. However, for sufficiently high acoustic perturbation, the liquid core may break-up into structures too small to be correctly described by the local mesh refinement. In such a case, the volume fraction is spatially diffused which results in a local decrease of ∇α l . This can be observed on figure 15, representing slices of ∇α l at different instants for the NMC and the MC. By comparing figures 15b and 15c, a strong diffusion of the liquid structures is noticed at the tip of the liquid core, leading to the homogenization of ∇α l . This intervenes after the decrease of the liquid core length due to its stripping by acoustic radiation pressure and actually explains the drop of resolved interfacial area density observed in figure 14. It can also be seen, by comparing figure 15a to figures 15b and 15c, that the liquid expansion in the direction orthogonal to the a.a. is far greater in the MC than in the NMC. This is actually also due to the flattening of the liquid core and the promotion of its stripping by acoustic radiation pressure. Indeed, as observed by [START_REF] Ficuciello | Investigation of airassisted sprays submitted to high frequency transverse acoustic fields: Droplet clustering[END_REF], once the jet is flattened, small liquid ligaments are torn from the central core and ejected in the orthogonal direction to the a.a.. The SPS seems thus able to render the expansion of the liquid mass in this direction. Note that since the numerical mesh is too coarse in this region to discretize such small liquid structures, the liquid volume fraction is diffused. In such highly diffused zones, the loss of information about the interface shape and position as well as the fact that the velocities of the two separated phases are not known may deteriorate the evaluation of interactions between the gas and the liquid in the SPS. This is a well known limitation of the model. To quantify the amount of liquid that expands, let us define highly diffused zones by the double condition ∇α l < 100 and α l < 0.99. Based on this definition, it appears that approximately 15 % of the total mass of liquid in the SPS is found in these highly diffused zones in the MC, compared with only 3 % in the NMC. A way to gain accuracy in these zones without any additional mesh refinement would be to solve a transport equation of Σ SPS with specific source terms adapted to the subgrid modelling of the break-up process induced by shear stresses and the acoustic radiation pressure. Hence, the subgrid interfacial area density would bring useful information to enhance the modelling of the liquid behaviour. The use of an advanced modelling of the two fluids in the SPS giving access to their own velocities in the interface could also improve the overall strategy.

Response of the Dispersed Phase Solver

To complete the previous analysis, the amount of liquid interface in the DPS has to be measured as well. According to its definition, the interfacial area density of the spray in a given volume Ω reads Σ DPS, Ω = S Ω /Ω, with S Ω the total surface of droplets in the volume. By introducing the total volume of droplets V Ω , Σ DPS, Ω may be written as Σ DPS, Ω = α tot Ω S Ω /V Ω , with α tot Ω the volume fraction of droplets in Ω, all sections k considered. In addition, the introduction of the Sauter mean diameter of the spray D 32, Ω , which corresponds to the diameter of a single droplet with the same ratio S/V as the whole spray, leads to :

Σ DPS, Ω = 6 α tot Ω D 32, Ω . (6.1)
Based on this definition, the time evolution of Σ DPS, Ω is plotted on figure 16a for both the NMC and the MC. For both cases, an initial drop can be observed at the very beginning of the plot. This is actually due to a change in liquid properties operated at the end of the first simulation dedicated to the establishment of the flow, which is not presented in this paper. However, this modification has a minor influence on the simulation. Indeed, in the NMC, the trend of Σ DPS, Ω quickly establishes around a value of 3 m -1 , with very low fluctuations in time.

Regarding the MC, the plot reveals an increasing value of Σ DPS, Ω all along the simulation. Following the equation (6.1), this increase is either due to a growth of α tot Ω or to a drop of D 32, Ω . According to figures 16b and 16c, the volume fraction α tot Ω , and consequently the total mass of droplets, follow the same trend than Σ DPS, Ω while the Sauter mean diameter D 32, Ω remains quasi constant. Hence, this reflects an intensification of the mass transfer from the SPS toward the DPS once the jet is submitted to acoustic perturbations. This can be attributed to the increase of the overall interfacial area in the SPS (sum of the resolved and unresolved components) highlighted on figures 14 and 15. Since the primary atomization source term (equation (2.5)) depends on α l but not on ∇α l , the mass transfer from the SPS toward the DPS may increase even if the resolved component of the interfacial area density diminishes in the SPS.

For a better quantification of the phenomenon, let us define the net mass transfer rate per unit volume between the SPS and the DPS as:

S net = S a - Ns k=1 s k c , (6.2)
and the time averaged net mass transfer rate integrated over the whole volume Ω as:

ΩS net, Ω = Vi∈Ω S net V i . (6.3)
Computing the time average over 13 acoustic periods, it appears that the mean mass transfer rate is 2.24 g.s -1 under acoustic modulation, against 1.60 g.s -1 in the NMC. This increase of 40 % is significant and may have a drastic impact on the combustion stability in typical LRE configurations.

The previous analysis thus demonstrates the ability of the coupling strategy to render the increasing trend of the primary atomization process by acoustic radiation stresses, although it does not yet constitute a quantitative validation of the model.

Dynamics of the spray under acoustic modulation

Spray angle modification

An additional impact of acoustic perturbations on the two-phase flow is the drastic expansion of the spray angle in the direction orthogonal to the a.a., as already observed by [START_REF] Baillot | Behaviour of an airassisted jet submitted to a transverse high-frequency acoustic field[END_REF] and [START_REF] Ficuciello | Investigation of airassisted sprays submitted to high frequency transverse acoustic fields: Droplet clustering[END_REF]. In order to quantify this mechanism in the current simulations, the spray angle in the aforementioned direction has been defined as shown on figure 17a and reads:

γ + = Vi∈Λ+ γ α tot V i Vi∈Λ+ α tot V i , γ -= Vi∈Λ- γ α tot V i Vi∈Λ- α tot V i , (7.1)
where γ corresponds to the angle between the injection axis and the segment connecting the injector lip to the centre of the control volume V i , while α tot is the local volume fraction of droplets, all sections considered. Λ + and Λ -are domains of positive and negative Y-coordinates respectively, with the injection axis as reference. The spray angle is thus weighted by the local volume of droplets, and consequently its mass. Thanks to this definition, no threshold on the minimum volume fraction of droplets to measure is needed contrary to usual experimental measurements based on back-light images. The time evolution of the spray angle in the direction orthogonal to the a.a. is represented on figure 17b for the NMC and the MC. A significant increase of the angle is noticed on both directions once submitted to acoustic perturbations, while it remains quasi constant throughout the NMC. After approximately 10 acoustic periods, the angle establishes around ± 20 • which corresponds to twice the value without acoustic modulation.

To explain this phenomenon, figure 18 represents two-dimensional fields of S net for the NMC and the MC. The time average operation is performed over 13 acoustic periods. In addition, two-dimensional fields of the time averaged liquid mass fraction Y l are displayed, as well as iso-lines of Y l = 0.4. On this figure, the transverse expansion of the liquid mass fraction under the flattening of the liquid core is observed. This corresponds to ligaments being ripped from the liquid core under acoustic perturbations but too small to be correctly discretized by the mesh. Since the primary atomization source term (equation (2.5)) has been formulated to operate in the gaseous side of the interface (i.e. Y l → 0 on figure 2), the location of mass transfer toward the DPS consequently tends to depart from the injection axis, which in turn directly impacts the spray angle. This can be retrieved on figure 18, where zones of positive values of S net follow the expansion of the iso-line Y l = 0.4. Therefore, this demonstrates the ability of the coupling strategy between the SPS and the DPS to render the increasing trend of the spray angle once submitted to acoustic perturbations. However, comparisons with experimental data are needed to adjudicate on the modelling precision of the phenomenon. Indeed, the diffusion of the interface in the SPS might result in an overestimate of the volume taken by ligaments around the liquid core, and consequently distort the radial expansion of the spray.

Transverse oscillations

The last item of this study concerns the response of the spray to the acoustic velocity. Due to compression and decompression mechanisms, the propagation of acoustic waves results in both pressure and velocity fluctuations in the medium. In the case of the standing acoustic wave simulated in this paper, acoustic velocity reaches an amplitude of around 12 m.s -1 at the location of injection, i.e. at an IAN. This amplitude is not negligible compared to the gas injection velocity of around 70 m.s -1 and results in a transverse periodic movement of the gas. According to [START_REF] Méry | Experimental investigation of cryogenic flame dynamics under transverse acoustic modulations[END_REF], these velocity fluctuations may induce periodic oscillations of droplets around their injection axis and contribute to the coupling between flames and acoustics in a reactive case. Therefore, it is of first interest to be able to reproduce this transverse movement of the spray in numerical simulations. In order to measure the global response of the spray to the acoustic velocity, the spatial mean velocity of droplets in direction X (that of the a.a.) is calculated as:

< v k X > ρ = N cell i=1 v k X ρ k V i N cell i=1 ρ k V i (7.2)
where < > ρ refers to the spatial mean weighted by the local density of droplets ρ k , while v k X is the transverse velocity of droplets of section k in the control volume V i . Measurements are performed over two transverse slices located at a distance of 0.2 D l and 2.7 D l from the injection plane respectively, N cell being then the number of mesh cells intersected by a given slice.

Figure 19 represents Discrete Fourier Transforms (DFT) of the time evolution of < v k X > ρ over 18 acoustic periods for small, medium and big spray sections, for both transverse slices. At Z = 0.2 D l , droplets undergo transverse oscillations at the three modulation frequencies once submitted to the acoustic velocity, whatever the spray section. The amplitudes of oscillation for each mode depend on their inertia: the greater the inertia, the smaller the transverse movement amplitude. This response mechanism is retrieved thanks to the coupling strategy between the gas and the spray solvers, through the drag force term F k D (see section 2.3.2). Regarding the transverse slice at Z = 2.7 D l , it appears that droplets do not respond significantly to the 4T and 6T acoustic modes. In addition, the amplitude of oscillation at the main frequency f ac, 1 is greater compared to the one at Z = 0.2 D l , whatever the droplet size (+300 % for small droplets ; +100 % for medium and big droplets). This difference in the response amplitude is mainly caused by the fact that, due to their inertia, droplets do not reach instantly their maximum oscillation amplitude. Therefore, the oscillation of droplets seems to significantly depend on both their size and their position on the injection axis. Since the flame response in LRE configurations directly depends on the dynamics of the spray, the description of the droplet size polydispersion by the DPS appears essential in the framework of CFD analysis of combustion instabilities.

Conclusion

In this work, the unsteady simulation of a non-reactive two-phase flow submitted to a high-frequency transverse acoustic modulation has been performed. The methodology of [START_REF] Gaillard | Numerical simulation of cryogenic injection in rocket engine combustion chambers[END_REF] has been retained because it permits to describe the dynamics of the whole jet in an Eulerian framework, from injection to atomization processes. This methodology consists in the coupling between a diffuse interface method for the simulation of large liquid structures, and a kinetic-based Eulerian model for the description of droplets. The primary atomization source term on which is based the coupling between the solvers were formulated without any consideration for possible acoustic effects on atomization processes, while these may have a great impact on the spray properties. The aim of the present contribution was thus to evaluate the ability of such coupling strategy to render all the response mechanisms of two-phase flows submitted to transverse acoustic perturbations, typical of what actually occurs in liquid rocket engines under unstable operating conditions. The selected configuration corresponds to the one of [START_REF] Ficuciello | Investigation of airassisted sprays submitted to high frequency transverse acoustic fields: Droplet clustering[END_REF], which already brought much information on the response of two-phase flows under such acoustic forcing. Two simulations have been performed in this paper. The first one has been carried out without any acoustic modulation and used as reference, while the second one corresponds to the jet under transverse acoustic forcing. The results of these two simulations have thus been compared to one another as well as to experimental observations available in the literature.

From these simulations, it appears that the flattening of the liquid core by the heterogeneous profile of acoustic radiation pressure observed experimentally is well retrieved by the diffuse interface method. This flattening process gives rise to an early stripping of the liquid core associated to a release of small liquid structures at its tip and in the plane orthogonal to the acoustic axis. Due to the mesh refinement, this change in the primary atomization process results in a diffusion of the liquid volume fraction and may contribute to a loss of precision in the transport of the liquid structures thus created. To gain accuracy, the resolution of an additional transport equation of the interfacial area density could be adequate. However, even though information about the interface position and shape are lost, it is found that the current formulation of the mass transfer source term between both solvers is adapted to render the intensification of the primary atomization. In addition, the expansion of the liquid structures by acoustic radiation stresses directly results in the increase of the spray angle in the plane orthogonal to the acoustic axis. This phenomenon has already been observed in the literature and may have a great impact on the flame response in reactive cases. Finally, the transverse periodic motion of droplets under the action of the acoustic velocity is also captured thanks to the coupling between the gaseous phase and the spray through drag forces. Significant differences in the amplitudes of droplet oscillations have been noticed depending on both their size and their location on the injection axis. Since the flame response in liquid rocket engines is directly related to the dynamics of the spray, the description of the droplet size polydispersion by the spray solver turns out to be essential in the numerical study of combustion instabilities.

Future work is planned to complete the present contribution with quantitative comparisons to ongoing experimental measurements. In addition, the same kind of numerical simulation will be performed in a reactive framework in order to improve understanding of the driving mechanisms of HF combustion instabilities in subcritical operating conditions. These simulations are also expected to provide a basis for an accurate low-order modelling of phenomena.

Figure 1 :

 1 Figure 1: Schematic representation of the coupling strategy between the SPS and the DPS for the simulation of atomized two-phase flows.

Figure 2 :

 2 Figure 2: Activation functions for primary atomization and impingement source terms. λ a ; λ c .

Figure 3 :Figure 4 :

 34 Figure 3: Experimental test-rig of Ficuciello et al. (2017) coloured by normalized pressure amplitudes. (a) View containing the a.a. (b) View orthogonal to the a.a. Position of an IAN. Simulated part of the cavity.

Figure 5 :

 5 Figure 5: Pressure fluctuations at a PAN of the cavity.Raw experimental signal and numerical reconstruction for the acoustic forcing.

Figure 6 :

 6 Figure 6: (a) Experimental and (b) numerical snapshots of the two-phase flow near the injector exit (NMC).

Figure 7 :

 7 Figure 7: Longitudinal slice coloured by α l in grey scale (NMC).Iso-lines of α l, 1.0 and α l, 0.5 .

Figure 8 :

 8 Figure 8: Liquid core representation with blue iso-surfaces of α l, 0.99 and transversal cuts of α l in grey scale (black: α l = 1 ; white: α l = 0). (a) NMC and (b) MC at t = 18 T ac .

Figure 9 :Figure 10 :

 910 Figure 9: Theoretical p * rad distribution around a cylinder located at an IAN of a standing acoustic field. Red zones: compression ; blue zones: suction.

Figure 11 :

 11 Figure 11: Time evolution of the flattening parameter f at different distances Z from the injection plane. NMC ; MC.

Figure 12 :

 12 Figure 12: Time evolution of the intact liquid core length.α l, 0.9 , α l, 0.99 α l, 1.0 (black: NMC ; red: MC).

Figure 13 :

 13 Figure 13: Iso-surfaces of α l, 0.99 with transversal cuts of the liquid interface at Z = 4.5 D l ( NMC ; MC). (a) t = 11.5 T ac , (b) t = 12.5 T ac .

Figure 14 :

 14 Figure 14: Time evolution of the resolved interfacial area density in the SPS. NMC ; MC.

Figure 15 :

 15 Figure 15: Slices orthogonal to the a.a. coloured by ∇α l .Iso-lines of α l, 0.99 .

Figure 16 :

 16 Figure 16: Time evolution of (a) the interfacial area density in the DPS, (b) the total volume fraction of the spray and (c) the Sauter mean diameter. NMC ; MC.

Figure 17 :Figure 18 :

 1718 Figure 17: (a) Definition and (b) time evolution of the spray angle in the direction orthogonal to the a.a.. NMC ; MC.

Figure 19 :

 19 Figure 19: DFT of < v k X > ρ at two distances Z from the injection plane and for the three spray sections. NMC ; MC.

Table 1 :

 1 Mode fac, i [Hz] pac, i [Pa] ψac, i [rad] Characteristics of the acoustic signal at a PAN.

	1 2T	≈ 1000	5135	0
	2 4T	2 fac, 1	887	5.44
	3 6T	3 fac, 1	555	4.97

† Email address for correspondence: luc-henry.dorey@onera.fr

Acknowledgments

This work is supported and co-supervised by Cnes the French Space Agency and was granted access to the HPC resources of Cines under allocations 2016-c20162b7710 and 2017-A0022B07710 made by Genci. We gratefully thank Coria for sharing data as well as experience. The on-going work will be the occasion to maintain the joint effort between both labs.