François Doisneau
email: francois.doisneau@onera.fr.

Alaric Sibra
email: alaric.sibra@onera.fr

Joel Dupays
email: joel.dupays@onera.fr.4

Angelo Murrone
email: angelo.murrone@onera.fr.

Frédérique Laurent
email: frederique.laurent@ecp.fr.

Marc Massot
email: marc.massot@ecp.fr.

Numerical strategy for unsteady two-way coupling in polydisperse sprays: application to Solid Rocket Motor instabilities

come

I. Introduction

The accurate simulation of two-phase flows constituted of a moderately dense polydisperse liquid or solid phase in a carrier gaseous flow field is crucial in many applications and especially for aluminized solid propellant combustion. The importance of polydispersity in such two-phase flows has been demonstrated in general cases [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] and more specifically in solid propulsion cases [START_REF] Doisneau | Two-way coupled simulation of acoustic waves in polydispersed coalescing two-phase flows : application to Solid Rocket Motor instabilities[END_REF][START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF]. In a solid rocket motor (SRM), droplet mass loadings of the aluminum combustion residuals are high enough to induce a significant two-way interaction between the phases in terms of momentum and heat. Besides, the typical density and sizes of the droplets are such that this level of mass loading can also lead to collisions or coalescence. However, the volume occupied by the liquid or solid phase remains a very small fraction of the total volume and does not need to be taken into 1 account explicitly in the carrier phase modeling. The resulting liquid or solid phase is then disperse and we will describe it as a spray. Such a spray will then be referred to as moderately dense since (1) the volume fraction does not play a role, and (2) droplet-droplet interactions, even if present, will not be the main phenomena driving the trajectories, as opposed to dense sprays or fluidized beds.

The strong coupling that occurs in moderately dense sprays between the carrier gas and the spray requires well-suited numerical methods, all the more for unsteady flows. Among unsteady cases, two-phase acoustics is a reference but a tough problem since one needs to capture accurately: (i) dissipation and dispersion of acoustic waves in two-phase media, (ii) acoustic disturbance yielded by droplet transport and segregation in the fluid, (iii) interaction of the gaseous flow non-linearities with the dynamics of droplets. Many two-phase approaches from the literature encounter difficulties in order to account for (i), which is though an important feature. Some aspects of these methods are described and the reasons for failure are exposed. We then focus on this point in the more general case of polydispersity, which we refer to as polydisperse acoustics. The impact of polydispersity on acoustics-spray interaction is presumably complex but few studies are available and many authors suggest equivalent diameter approaches such as in [START_REF] Peirano | A mean diameter for numerical computations of polydispersed gas-solid suspensions in fluidization[END_REF], which is not suitable when the physical phenomena involve very different levels and time scales of interaction between the phases depending on the size of the droplets.

The Williams-Boltzmann equation, a transport equation based on kinetic theory, has proven to be useful for treating dilute and moderately dense polydisperse phases. Such an equation describes the evolution of the number density function (NDF) and it is coupled through transfer terms to a fluid model for the gas, i.e. conservation equations such as Euler or Navier-Stokes systems. The statistical description of the disperse phase coupled to the fluid description of the carrier phase yields a mesoscopic scale description, which will be the background of our contribution as a whole. The resolution of these fluid-kinetic coupled equations will be referred to as Direct Numerical Simulation (DNS), even if some level of modeling is involved and the droplet (microscopic) scales are not resolved. Besides, we will not introduce any turbulence modeling, such as turbulent dispersion models [START_REF] Simoes | Modélisation eulérienne de la phase dispersée dans les moteurs à propergol solide, avec prise en compte de la pression particulaire[END_REF] and sub-grid scale closure terms for both the gas and the droplet phases, as it is required in the context of Large Eddy Simulation (LES) [START_REF] Boileau | Robust numerical schemes for Eulerian spray DNS and LES in two-phase turbulent flows[END_REF][START_REF] Masi | Étude théorique et numérique de la modélisation instationnaire des écoulements turbulents anisothermes gaz-particules par une approche Euler-Euler[END_REF][START_REF] Vié | Simulation aux grandes échelles d'écoulements diphasiques turbulents à phase liquide dispersée[END_REF][START_REF] Wunsch | DNS/DPS of inertial droplet coalescence in homogeneous isotropic turbulence and comparison with pdf model predictions using the direct quadrature method of moments[END_REF]]. Yet, efficient DNS is an incontrovertible starting point for disperse two-phase flow simulations while the adaptation of the model and methods to LES relies most of the time on theory and DNS tools. Thus, the framework of the present contribution will be the one of DNS of this fluid-kinetic coupled system of equations, even if the final simulations proposed in the context of the industrial-oriented code CEDRE will only be under-resolved and presented for feasibility purposes as a starting point for larger scale simulations and LES developments.

As an alternative to Lagrangian approaches, which are difficult to couple to the Eulerian representation of the gas phase [START_REF] Capecelatro | An Euler-Lagrange strategy for simulating particle-laden flows[END_REF], the Eulerian Multi-Fluid model (MF) [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] discretizes the NDF using a finite volume approach of the size variable but conditions velocity by size in order to reduce the phase space dimension. After integration on size intervals, the so-called "sections", the resulting conservation equations are similar to those of the pressureless gas dynamics with source terms accounting for gas-droplet interactions. The reconstruction of the size distribution allows to treat rigorously the size dependency of all the terms of a given problem since the kinetic equation is integrated over the sections, performing a partitioning of the whole size space. Properties, such as the conservation of a set of local moments, are then enforced independently on each section and the combination of the reconstructed NDFs converges towards the exact NDF. This allows the MF method to deal with phenomena such as coalescence [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] -which modeling is highly needed in SRM simulations [START_REF] Doisneau | Eulerian VS Lagrangian simulation of unsteady two-way coupled coalescing two-phase flows in solid propellant combustion[END_REF]-contrary to so-called class methods [START_REF] Bendakhlia | Multiradii modeling of counterflow spray diffusion flames[END_REF][START_REF] Konan | Reactive multiphase flow simulation of Uranium Hexafluoride conversion reactor[END_REF][START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]. Moreover the continuous treatment of the size variable is done in a way that is consistent at the section level to the treatment of other phenomena (drag, heating), which makes a strong difference with the classical methods described for Population Balance Equations (PBEs) in [START_REF] Kostoglou | Evaluation of zero order methods for simulating particle coagulation[END_REF] and in [START_REF] Kumar | On the solution of population balance equations by discretization -II. A moving pivot technique[END_REF], which conserve a few global properties for a given phenomenon. The MF model treats particle inertia independently in each section so it can capture, with respect to polydispersity, acoustics and unsteady hydrodynamics thanks to the fact that it is spatially resolved.

So the MF method achieves an accurate treatment of size distributions for a large size range including inertial particles/droplets in a CFD framework: it is more accurate and rigorous than other Eulerian methods as regards polydispersity when multiple phenomena including inertia occur.

But computing transport in physical space for both the gas and the droplet sections, and eventually capturing the coupling at all time scales, is computationally expensive. So we seek a MF method:

• that features an accurate treatment in the sections in order to reduce their number,

• that allows to avoid taking a small time step when small, low inertia droplets are involved.

For size accuracy high order MF methods exist with reconstructed size distributions converging towards the NDF at order one [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF], two [START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF], or even four [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF]. For coupling and time accuracy efforts must be done to allow the use of larger time steps. In the context of dense sprays, the stiffness of small droplets has been identified and occurs similarly by means of short velocity and heat relaxation times [START_REF] Burman | Numerical analysis of two operator splitting methods for an hyperbolic system of conservation laws with stiff relaxation terms[END_REF]. Numerical approaches based on modifications at the Riemann solver level allow to capture accurately all the waves in the asymptotic limit of small droplets [START_REF] Béreux | A Roe-type Riemann solver for hyperbolic systems with relaxation based on time-dependent wave decomposition[END_REF] but it is complex to design methods suited and accurate for all the size range [START_REF] Béreux | Zero-relaxation limit versus operator splitting for two-phase fluid flow computations[END_REF].

We therefore introduce, for moderately dense sprays, a time integration method that has an extended stability domain and we assess its accuracy with the objective of designing a tunable accuracy/cost method. Moreover the method is based on operator splitting, the flexibility of which makes it comfortable to implement, opened to explicit and implicit time integration of each operator and can be generalized to other types of couplings. Regarding the latter point, the coupling among the sections yielded by coalescence is indeed investigated in [START_REF] Doisneau | Eulerian VS Lagrangian simulation of unsteady two-way coupled coalescing two-phase flows in solid propellant combustion[END_REF] but a two-way coupling with an external phase can also be considered e.g. radiation or three phase flows. The paper first addresses the specificities of polydisperse two-phase acoustics and the conditions for a Multi-Fluid-type method to account for them. The coupling difficulties are then identified and analyzed in order to conceive a numerical method based on sound arguments. Controlling accuracy appears crucial for methods to be used in both research and industrial contexts. The developed strategy is implemented in a 1D research code called SAP1 for verification purposes on two-phase acoustics cases and in an industrial-oriented code for validation on general cases and feasibility demonstration. The code, called CEDRE and developed at ONERA, is used in energetics and aeronautics, and particularly for SRM simulation [START_REF] Refloch | CEDRE software[END_REF]. A complex, unsteady SRM featuring instabilities is therefore computed and exploited, thus showing the potential of the method.

The paper is organized as follows: § II recalls the derivation of the Multi-Fluid model, its hypotheses, and the way it captures polydispersity. The different time scales featured by the model and the level of coupling between all the phases are analyzed, as a guide to an efficient resolution. The model's ability to account for two-phase dissipation and dispersion is discussed. § III describes an analytical approach, which has been used in the literature to solve linear acoustic propagation in non-rotational sprays: it is extended to polydispersity and is then used as a reference solution. The physics captured by this theory is proven to be well approached by the Multi-Fluid model, the convergence criteria being discussed. § IV reviews the numerical requirements and peculiarities of spray-acoustics simulations. A new numerical strategy called ACS, based on operator splitting, is deduced from what precedes to solve strong coupling with a tunable accuracy/cost trade off.

§ V provides additional quantitative arguments to parametrize the new ACS algorithm. The new strategy is implemented in a research code to assess the limits and requirements for spray-acoustics interactions. It is then implemented in the industrial-oriented code CEDRE and validated on acoustic test cases. § VI shows the results of an unstable SRM simulated with the new strategy as implemented in CEDRE, which prove the method to be robust and accurate. The numerical strategy for SRM simulations is assessed before concluding.

II. Modeling unsteady moderately dense polydisperse flows

In this section, the derivation of the Multi-Fluid model as a description of the disperse phase and the coupling to the gas equations are recalled. The dissipation of the global system is defined. The ability of this modeling to handle acoustic waves in polydisperse sprays is then discussed in detail as well as the required size discretizations to do it accurately. The numerical peculiarities of moderately dense two-phase coupling are finally investigated and benchmark criteria for research and industrial codes are provided.

A. The Eulerian Multi-Fluid model

We model the spray with a system of coupled conservation equations for the disperse phase and for the gas. The gas carrier phase is described with the Euler equations with ideal fluid equation of state. These equations, generally used at high Reynolds number, are nonetheless sufficient to capture two-phase acoustics but the numerical methods that are hereafter developed can be generalized straightforwardly to Navier-Stokes equations, which are used for the final SRM simulations. The droplet phase is described with the Eulerian Multi-Fluid (MF) model, which discretizes the size phase space using sections, i.e. size intervals, and then amounts to a finite volume method [START_REF] Laurent | Numerical analysis of Eulerian Multi-Fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF][START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF].

Kinetic description of the disperse phase

We now recall its derivation from the kinetic level, where droplets are point particles for which the interactions with the carrier gas are modeled. Let us define the number density function (NDF) f of the disperse phase, where f (t, x, u, T, S)dxdudT dS denotes the average number of droplets (in a statistical sense), at time t, in a volume of size dx around a space location x, in a dudT dSneighborhood of velocity u, internal temperature T and droplet surface S. The evolution of the NDF is described by the Williams-Boltzmann equation. Considering local momentum and heat transfer with the gas due to drag forces and conducto-convective fluxes and considering collisions between droplets but no evaporation nor break-up, it reads:

∂ t f + u • ∂ x f + ∂ u • (Ff ) + ∂ T H c p,l f = Γ (1)
where F and H are the drag force and heat transfer per unit mass, c p,l is the pressure specific heat capacity and Γ is the collision source term. Though collisions and coalescence are important features of moderately dense sprays in a SRM [START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF], the method is exposed without them, keeping in mind that the adaptation is straightforward. The transfer terms require models to take into account physics at the droplet scale. They are typically modeled by Stokes' laws with eventually Schiller-Naumann or Ranz-Marshall corrections. Whatever models are used, characteristic dynamic and thermal times (or relaxation times) τ u (S) and τ T (S) can be defined from F and H, that strongly depend on droplet size.

The semi-kinetic level

In a first step, we reduce the size of the phase space to the only droplet size variable. We therefore consider moments in velocity and temperature variables of order zero and one, conditioned by size: the droplet number density n = f dudT , the average velocity u = n -1 uf dudT and the average enthalpy h = n -1 h(T )f dudT . They depend only on (t, x, S). We can now define the effective temperature T so that h = h(T ). In order to close the system, the following assumptions are introduced:

[HV1] For each droplet size S, at a given point (t, x), the only characteristic velocity is the average u(t, x, S).

[HV2] The velocity dispersion around u(t, x, S) is zero in each direction, whatever the point (t, x, S).

[HT1] For each droplet size S, at a given point (t, x), the only characteristic temperature is the average T (t, x, S).

[HT2] The temperature dispersion around T (t, x, S) is zero whatever the point (t, x, S).

It is equivalent to presume the following form of NDF:

f (t, x, u, T, S) = n(t, x, S)δ(u -u(t, x, S))δ(T -T (t, x, S)). (2)
The set of hypotheses [HV1] and [HV2] known as the monokinetic hypothesis has been introduced in [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] and is equivalent to reducing the velocity distribution support to a one dimensional sub-manifold parameterized by droplet size. It is correct when τ u (S) is small compared to a gas characteristic convective time [START_REF] De Chaisemartin | Polydisperse evaporating spray turbulent dispersion : Eulerian model and numerical simulation[END_REF][START_REF] De Chaisemartin | Eulerian models for turbulent spray combustion with polydispersity and droplet crossing[END_REF][START_REF] Massot | Eulerian Multi-Fluid models: modeling and numerical methods[END_REF], the ratio of these two times forming the dynamic Stokes number St c .

These two hypotheses render the droplet convection and the velocity dependent source terms satisfactorily for Stokes numbers below one as demonstrated extensively by comparisons to Lagrangian references in one-way coupled complex combustion cases [START_REF] Fréret | On the ability of the Eulerian multifluid model to predict preferential segregation and flame dynamics in polydisperse evaporating sprays[END_REF][START_REF] Fréret | On the role of preferential segregation in flame dynamics in polydisperse evaporating sprays[END_REF]. The question of moderately dense sprays of more inertial droplets remains tough and open [START_REF] Doisneau | Eulerian modeling and simulation of small scale trajectory crossing and coalescence for moderate-Stokes-number spray flows[END_REF][START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial particle-laden flows[END_REF][START_REF] Vié | On the direct numerical simulation of moderate-Stokesnumber turbulent particulate flows using algebraic-closure-based and kinetic-based moment methods[END_REF]. Hypotheses [HT1] and [HT2] are similarly introduced in [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] so that temperatures are locally correlated for droplets of a given size. They are correct when the thermal Stokes number, formed with τ T (S), is below one, as well as the dynamic Stokes number. Once again for moderately inertial to inertial droplets, errors would occur in temperature-dependent source terms such as evaporation or heating retroaction in a two-way coupling context. This step leads to a system of conservation equations called the semi-kinetic model, which reads:

   ∂ t n + ∂ x • (nu) = 0 ∂ t (nu) + ∂ x • (nu ⊗ u) = nF ∂ t (nh) + ∂ x • (nhu) = nH (3)
where one gets the average transfer terms nF = Ff dudT and nH = Hf dudT . The phase space of the semi-kinetic system has still too high a dimension to allow direct resolution and requires further modeling.

The Multi-Fluid models

The second modeling step is a finite volume discretization of the size phase space that leads to the Eulerian Multi-Fluid model. Multi-Fluid models rely on the choice of a discretization 0 = S 0 < S 1 < • • • < S N = ∞ for the droplet size phase space and the averaging of the conservation law system over each fixed size interval [S k-1 , S k [, called section [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]. The set of droplets in one section can be seen as a "fluid" for which conservation equations are written, the sections exchanging mass, momentum and enthalpy. In order to close the system, the following assumptions are introduced:

[HV3] In each section, the velocity does not depend on the size of the droplets.

[HT3] In each section, the temperature does not depend on the size of the droplets.

[HS1] In each section, the form of n as a function of S is presumed.

We choose for assumption [HV3] the notation u(t, x, S) = u k (t, x) to designate the constant velocity distribution in section k. Hypothesis [HT3] is similarly introduced and is equivalent to presuming h(T (t, x, S)) = h k (t, x) = h(T k (t, x)) as a constant enthalpy distribution in section k, and the corresponding effective temperature T k (t, x) which allows to define the sectional specific heat capacity c p,k = c p,l (T k ). The validity of these assumptions is linked to the strength of polydispersity in each section, which is quantified in a section by comparing the smallest to the biggest dynamic Stokes numbers. If the dynamic Stokes number spectrum is too wide, [HV3] and [HT3] no longer hold while only [HT3] fails for a wide thermal Stokes number spectrum. In both cases, the discretization must then be refined, as discussed in the following in the particular case of acoustics in polydisperse two-phase flows, or high order methods such as the Coupled Size-Velocity Moment (CSVM) method [START_REF] Vié | Size-velocity correlations in hybrid high order moment/multi-fluid methods for polydisperse evaporating sprays: Modeling and numerical issues[END_REF] must be used.

As for [HS1], it is expressed with size presumed form functions κ k : this assumption allows to reduce the size distribution information in each section at (t, x) to a set of moments of S, the number of which depends on the choice of the (κ k ) k set. Let us consider two methods based on different forms of presumed functions for [HS1]. In the One Size Moment method (OSM), a one parameter function in each section κ k , decouples size dependence 1 κ k (S), noted with the left-exponent 1, and space-time dependence m k (t, x), homogeneous to a mass concentration:

n(t, x, S) ≈ k m k (t, x) 1 κ k (S)I [S k-1 ,S k [ (S) (4)
where

I [S k-1 ,S k [ (S)
is the characteristic function, equal to one in section k and to zero elsewhere.

With the constraint

S k S k-1 κ k (S)dS = 1
, this yields a first order size convergence with the number of sections and corresponds to the classical Multi-Fluid method which has been developed and validated in evaporating cases [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] and in coalescing cases [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. The Two Size Moment method (TSM) uses a two-coefficient reconstruction in each section, noted with a 2 left-exponent as follows 2 κ k (t, x, S), which yields a second order size convergence with the reconstruction:

n(t, x, S) ≈ k 2 κ k (S)I [S k-1 ,S k [ (S). (5)
and with the two constraints on the size moments of order zero and 3/2:

n k = S k S k-1 2 κ k (S)dS m k = S k S k-1 ρ l S 3/2 6 √ π 2 κ k (S)dS (6)
where n k is the number concentration and m k the mass concentration. The general method has been validated for evaporating cases [START_REF] Laurent | Numerical analysis of Eulerian Multi-Fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF] and for coalescence, especially in solid rocket motors [START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF].

Although the two methods both have a first order reconstruction in velocity [HV3] and temperature [HT3], OSM computes drag and heating with a fixed reconstruction while TSM computes these terms with a size distribution that is better described inside each section. The TSM method achieves indirectly a better approximation of velocity and temperature by estimating the characteristic times more accurately. Moreover the quality of the size reconstruction has a strong impact on the velocity and temperature errors since the characteristic dynamic and thermal times τ u and τ T usually depend on the square of the droplet radius.

Both methods then solve for size moments calculated in the sections, the equations of which are derived by integrating the semi-kinetic System (3). Please note that the more parameters the presumed function has, the more moments one needs. The choice of refining the size distribution description with a high order moment method for a better resolution of polydispersity therefore increases the computational cost.

Equations of the One Size Moment Multi-Fluid method

The TSM method is implemented in the CEDRE code. But when no size modifying phenomena occurs at any time and location, TSM is equivalent to a purposely parameterized OSM method so the latter is the only method for which the equations are presented in this paper, for the sake of legibility. The OSM Multi-Fluid model leads to the following definition of the bulk density of droplets of surface S ∈ [S k-1 , S k [:

m k (t, x) = ρ l 6 √ π S k S k-1 S 3 2 1 κ k (t, x, S)dS (7)
where ρ l is the supposedly liquid droplet material density. When adding the gas equations on density ρ g , velocity u g and internal energy e g = c v,g T g to the N sec section equations, we obtain the following coupled system for moderately dense sprays:

                           ∂ t ρ g + ∂ x •(ρ g u g ) = 0 ∂ t (ρ g u g ) + ∂ x •(ρ g u g ⊗u g ) = -∂ x p - Nsec k=1 m k F k ∂ t (ρ g e g ) + ∂ x •(ρ g e g u g ) = -p∂ x •u g - Nsec k=1 m k H k + Nsec k=1 m k F k (u g -u k ) ∂ t m k + ∂ x •(m k u k ) = 0 ∂ t (m k u k ) + ∂ x •(m k u k ⊗u k ) = m k F k ∂ t (m k h k ) + ∂ x •(m k h k u k ) = m k H k    k = 1, N sec (8)
where F k and H k come from F and H, averaged on a section. The total droplet system thus counts n d + 2 times more equations than the number of sections with n d the dimension of the physical space. The gas source terms have been written so that System (8) satisfies the conservation of momentum and total energy for the mixture of gas and droplets.

In the following, we will assume constant heat capacities c p,g and c p,k for the gas and the sections and we will assume Stokes' laws for both F k and H k without loss of generality on the presented results. In Stokes' regime, the drag force per unit mass, due to the velocity differences with the gaseous phase, and the heat transfer, due to temperature differences, read for section k:

F k = u g -u k τ u k , H k = c p,k T g -T k τ T k (9)
where u g (t, x) and T g (t, x) are the gas local velocity and temperature. The section characteristic times for these transfers read, once integrated on each section k:

τ u k = ρ l (d 31 k ) 2 18ν g , τ T k = 3 2 c p,k c v,g Prτ u k ( 10
)
where ν g is the gas dynamic viscosity, Pr is the Prandtl number, ratio of the gaseous kinematic viscosity and heat diffusivity and d

31

k is a diameter averaged on section k from the general formula

( √ πd pq k ) p-q = S k S k-1 S p 2 1 κ k (t, x, S)dS S k S k-1 S q 2 1 κ k (t, x, S)dS . ( 11
)
Regarding further studies, we note diameters averaged on the whole size distribution d pq . For SRM applications, the gas-droplet transfer can occur, in the nozzle for instance, at high particular Reynolds numbers Re p = d p | u g -u k | /ν g so that Stokes' regime is not valid. In this case, we use Schiller-Naumann and Ranz-Marshall convective corrections [START_REF] Ranz | Evaporation from drops[END_REF][START_REF] Schiller | A drag coefficient correlation[END_REF] for drag and heat transfer respectively, which result in a transfer rate increase when velocity differences are strong.

B. Stiffness of polydisperse MF systems

The physics of two-phase flows is complex and multi-scale. Transport, drag and heat transfer are considered and they experience a large spectrum of time scales which yields stiffness, that must be accounted for when designing the numerical strategy.

The kinetic and semi-kinetic systems have a continuous range of time scales for convection through u(S) and for relaxation through τ u (S) and τ T (S). In the MF method, sizes are discretized so that velocities, temperatures and time scales are average for a group of droplets in the size interval. The discrete spectrum still spans a large range of frequencies so the system remains stiff. Yet we highlight the fact that discretization gathers the scales of the smallest and biggest droplets into the first and last sections respectively. This narrows the time scale spectrum. The level of size discretization that is required to capture the physics, and the treatment of stiffness are discussed in § III B. An overview of two-phase flow spatial discretization for resolution purposes is provided in § IV A 1. Once a space discretization is chosen, a time scale appears for the gas, that conditions stability and accuracy of the numerical methods. The most relevant scale for the gas, noted τ g , is defined as the time required for the fastest perturbation to cross the length of a mesh cell ∆x:

τ g = ∆x c 0 + | u g | . ( 12
)
This characteristic time is compared to the integration time step to build the CFL stability criterion, crucial for most numerical transport schemes. But it also has a physical meaning: the gas variables in a cell -which are coupled in a compressible approach-cannot change faster than τ g , if the fields are smooth.

A convective time scale can be defined similarly for each section, based on the average droplet velocity u k : This time scale only drives the droplet convection and is often much bigger than τ g . It is also the basis of the stability/accuracy criterion of the droplet phase transport scheme. In the following, we will use τ c k to refer to the convection scale of the section with the highest velocity, which will be the most problematic, keeping in mind that all these scales are generally slower than the gas scale τ g , except for supersonic flows.

τ c k = ∆x | u k | . ( 13
Each section has also a dynamic and a thermal relaxation time, which come from the drag and heat transfer laws. These times noted τ u k and τ T k are generally close so we consider the smaller time as the reference time scale of a section τ k :

τ k = min{τ u k , τ T k }. (14)
It is useful to distinguish droplets on their ability to react to acoustic disturbances so we call:

• big droplets those which relax slower than the gas varies in a cell (τ k > τ g );

• small droplets those which relax faster than the gas varies in a cell (τ k < τ g ).

For other applications, modeling droplet-droplet interactions, evaporation, etc. yield additional time scales. The example of growth times due to coalescence is given in [START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF]. They also add coupling links between the sections: collisions lead to momentum transfer, evaporation, coalescence or breakup lead to mass, momentum and heat transfer, etc. but these phenomena are not considered here.

In System [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF], the conservation equations of all the phases are coupled through the gas phase. This is schematically summed up in Figure 1.

So the problem features 3N sec + 1 time scales and is tightly coupled. This requires a careful strategy to render the two-phase physics. For polydisperse sprays, we define the minimum relaxation time scale:

τ min = min{τ k } (15)
so that only three time scales are relevant to characterize the coupling stiffness: the gaseous and fastest droplet convection scales τ g and τ c k and the relaxation scale τ min . The evolution of all these time scales depending on droplet size is given in Figure 8 for typical SRM. In these cases, meshes are defined with nozzle cell sizes that typically yield τ g > 10 -6 s. With usual alumina particles, the τ min of most resulting droplets generally remains above τ g but it is troublesome for micronic droplets and may fall below 10 -8 s for droplets coming from propellants with nano-aluminum particles. To sum up, stiffness emerges from the fast scales that lie in acoustics with τ g , and in polydispersity with τ min .

C. Dissipativity of two-phase systems

System (8) has been written to ensure the conservation of total momentum and total energy. Yet, drag and heat transfer between the droplets and the gas yields local dissipation regarding velocity and temperature differences, so that the system tends to relax towards a equilibrium state. This fact needs to be quantitatively rendered by the chosen two-phase model and method, which we discuss in the following sections.

Two-phase entropy equation and local equilibrium

We consider the gas specific entropy s g and the section specific entropies s k to define a mixture entropy ρs = ρ g s g + Nsec k=1 m k s k and its flux F s = ρ g s g u g + Nsec k=1 m k s k u k , the phase entropies being built analogously to a single phase entropy [START_REF] Woods | The Thermodynamics of Fluid Systems[END_REF]. The use of a global entropy for disperse phase flows is inspired from the work of [START_REF] Nigmatulin | Dynamics of multiphase media[END_REF]. We derive from Gibbs' relations for single phases the following entropy budget for System [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF], assuming smooth solutions:

∂ t ρs + ∂ x F s = Nsec k=1 m k T g τ u k (u g -u k ) 2 + Nsec k=1 m k c p,k τ T k (T g -T k ) 2 T g T k . (16)
The entropy Eq. ( 16) has two-phase source terms corresponding to the irreversibility of heat transfer and viscous effects in the droplet vicinity. These transfers at an unresolved level have a significant effect on the two-phase dynamics and have been modeled through drag and heating correlations.

Besides, the transfers within the bulk gas may not be significant compared to convective transfers, allowing to use Euler equations depending on the Reynolds number of the flow.

The entropy source terms are positive and cancel at dynamic and thermal equilibrium, i.e. the state where the velocities of all phases are equal to u ∞ and the temperatures are equal to T ∞ , which read:

u ∞ ρ eq = ρ g u g + k m k u k T ∞ c v,eq = ρ g c v,g T g + k m k c p,k T k ( 17
)
where ρ eq = ρ g + k m k and c v,eq = ρ g c v,g + k m k c p,k . The relaxation of the variables of System (8) towards this local equilibrium state is exact in 0D. This trend to reduce the velocity and thermal differences is responsible for the attenuation of acoustic waves which is the first specificity of two-phase acoustics, referred to as dissipation. The second specificity is dispersion and is due to velocity and temperature lags. Numerical methods which pretend solving two-phase acoustics should therefore capture the relaxation towards equilibrium efficiently. Typical numerical approaches are analyzed regarding this problematic in the following. In polydisperse sprays, the equilibrium state depends on the gas variables and all the disperse phase ones: this is a consequence of the full coupling among all phases, as explained in System (II B).

Failure of equilibrium models to capture dissipation

We discuss the necessity to solve the gas and droplet velocities and temperatures to properly capture dissipation and the limits of relaxed, or equilibrium, approaches such as the Equivalent Gas (EG) or the Eulerian Equilibrium Model (EEM). We also introduce a hybrid method referred to as Partial Equivalent Gas (PEG), where the smallest droplets only are relaxed.

A full equilibrium model, i.e. assuming a homogeneous density of droplets constantly at dynamic and thermal equilibrium with the gas, is a classical first approach for two-phase flows and is referred to as Equivalent Gas (EG). It practically consists in modifying the density and heat capacity of the gas for its equations to account for the mixture [START_REF] Rudinger | Fundamentals and applications of gas-particle flow[END_REF]. The model has also been described as the limit of dense two-phase flow model where the pressures and normal velocities are equal at interfaces [START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF] and is also known as homogeneous equilibrium model [START_REF] Kah | Taking into account polydispersity for the modeling of liquid fuel injection in internal combustion engines[END_REF]. No additional equations are then required for the disperse phase. EG is satisfactory only if all the Stokes numbers are very small compared to one. Because velocities and temperatures are assumed at equilibrium at all times, dissipation is not considered. In other terms, the entropy source terms of Eq. ( 16) depend on velocity and temperature differences, which are no longer resolved in the EG equations. In the particular case of acoustic waves -which is fully studied in the following section-no attenuation occurs and the speed of sound is constantly equal to c eq that is the speed of sound in a mixture of density ρ eq and heat capacity c v,eq . Since these facts are true only for very small droplets, EG fails to capture crucial two-phase acoustics features such as dissipation for most droplet sizes.

So the velocity and temperature differences are crucial variables for dissipation. When the droplet velocity field is not resolved, it can still be reconstructed based on a Taylor expansion of the drag force. The truncation error then increases with the Stokes number [START_REF] Ferry | A fast Eulerian method for disperse two-phase flow[END_REF]. In a similar way, the droplet temperature field can be estimated [START_REF] Ferry | Equilibrium Eulerian approach for predicting the thermal field of a dispersion of small particles[END_REF] so that data is available to approximate the dissipation due to velocity and temperature differences. This type of modeling, referred to as EEM in the previous references, is said to be extendable to two-way coupling [START_REF] Ferry | A fast Eulerian method for disperse two-phase flow[END_REF], where the velocities and temperatures are formally assessed with equivalent mixture density and specific heat: the approach is similar to EG so it departs as the Stokes numbers increase and the velocity and temperature differences that are then assessed still suffer a truncation error. Other expansions were performed by directly considering two-way coupling [START_REF] Druzhinin | Concentration waves and flow modification in a particle-laden circular vortex[END_REF][START_REF] Druzhinin | On the two-way interaction in two-dimensional particle laden flows: the accumulation of particles and flow modification[END_REF] but with a partial knowledge of the flow so they are available for simple configurations only. Anyway, the extension to polydispersity or droplet-droplet collisions was never performed. As a consequence of these remarks, near-equilibrium approaches are potentially not suited for capturing dissipation in industrial configurations.

Up to now, the two-way coupling is properly rendered only when both droplet and gas variables are solved such as in the dusty gas modeling [START_REF] Marble | Dynamics of dusty gases[END_REF], valid for monodisperse cases. The Multi-Fluid models have this feature and they are valid for polydisperse cases as well as for collision-coalescence.

An idea to avoid stiffness of small droplets is to combine full resolution strategies, such as MF, to EG in order to force the small droplets at equilibrium. Such techniques are referred to as Partially Equivalent Gas (PEG). An approach referred to as PEG-MF is considered in § VI where a part of the disperse phase mass at injection, corresponding to the droplets with St 1 everywhere, is transferred to the gas according to the EG technique and the rest is considered with its polydispersity with a MF technique. The same limitations exist as for EG i.e. the "unresolved" disperse phase should effectively be close to equilibrium at all time. We then highlight that small droplet evaporation, chemical reaction, coalescence or radiation can no longer be accounted for in this case since they are perfectly mixed and indistinguishable from the gas. This technique is classically used in non-reactive SRM simulations to get rid of the small droplet stiffness.

III. Linear acoustics in monodisperse and polydisperse sprays

An analytical model is considered for some simple cases of two-phase acoustics in order to investigate the relevance of MF methods in such contexts and validate them. The following linear acoustic case is a classical approach for waves in monodisperse sprays. It has been used to test numerical schemes [START_REF] Dupays | Contribution à l'étude du rôle de la phase condensée dans la stabilité d'un propulseur à propergol solide pour lanceur spatial[END_REF][START_REF] Morfouace | Two-phase flow analysis of instabilities driven by vortex-shedding in solid rocket motors[END_REF] in a monodisperse case.

It is here extended to polydisperse sprays in order to be the basis of quantitative evaluation of the new numerical strategy -this evaluation being performed in § V. Then, the ability of MF methods to account for polydisperse acoustics is shown and the effect of size discretization is studied to draw some firm conclusions on size discretization requirements.

A. Non rotational linear acoustic theory: a reference

We consider a monochromatic wave with a pulsation ω propagating in a spray that is homogeneous when at rest. We can determine analytically the solution by assuming a decaying plane wave structure. We define the speed of sound in the gas alone c 0 = √ γrT 0 where γ is the isentropic coefficient and r is the specific constant of the gas. The corresponding wave number is k 0 = ω/c 0 and the wave period τ ω = 2π/ω. Defining the complex wave number k = k 1 + ik 2 , we are interested in decaying harmonic solutions with the following form:

u g (t, x) = sin(ωt -k 1 x) exp(-k 2 x)u 0 g . (18)

Dispersion relation for a monodisperse spray

Temkin and Dobbins [START_REF] Temkin | Attenuation and dispersion of sound by particulate-relaxation processes[END_REF] have considered conservation equations for the gas and for a phase of monodisperse droplets, which are formally equivalent to a 1-section MF system, both interacting through drag and heat transfer. A dispersion relation for decaying harmonic waves can be derived for droplets of surface S: [START_REF] De Chaisemartin | Eulerian models for turbulent spray combustion with polydispersity and droplet crossing[END_REF] with C = m 1 /ρ g the droplet mass loading i.e. the ratio of droplet bulk density to gas density and κ = c p,l /c p,g the pressure heat capacity ratio. We note that k 1 /k 0 = c 0 /c(ω) with c(ω) the effective speed of sound in the spray. Temkin and Dobbins also suggest to define dimensionless coefficients to quantify attenuation and dispersion:

k k 0 2 = 1 + C 1 -iωτ u (S) 1 + Cκ 1-iωτ T (S) 1 + Cγκ 1-iωτ T (S)
       α = 2k 2 Ck 0 β = 1 C k 0 k 1 2 -1 (20)
where α is a dimensionless dissipation, independent of C for moderate mass loadings and β is a dimensionless dispersion, null for asymptotically big droplets. From this equation, two dimensionless groupings St ω = ωτ u (S) and St θ ω = ωτ T (S) naturally appear that are the acoustic Stokes number St ω and the thermal acoustic Stokes number St θ ω , which is generally close to St ω . We plot in Figure 4 the dissipation and dispersion coefficients as functions of the wave frequency, which is proportional to the dynamic Stokes number. At low St ω , the speed of sound converges to c eq that is the speed of sound in a homogeneous mixture of gas and droplets at equilibrium. It can be computed with the equivalent mixture density and the gas compressibility. It reaches a maximum for α , when St ω ∼ 1 for a monodisperse spray and at a slightly higher St ω , based on a d 31 , for the polydisperse sprays retained as examples and presented in Figure 2. Another local maximum, usually smaller, appears when the thermal and the dynamic times are separated by more than a decade. In all the following conditions, 25 µm diameter droplets yield St ω = 1 for a wave frequency around 700 Hz. Finally at very high Stokes number, the speed of sound is the one of the gas alone c 0 and the particulate phase is decoupled from the gas.

Dispersion relation for a polydisperse spray

An acoustic wave propagating in a quiet monodisperse or polydisperse spray yields drag force and heat transfer that are local and only depend on droplet size. The initial condition being the same for droplets of the same size, their velocity and temperature will evolve in the same way so the monokinetic assumptions [HV1], [HT1], [HT1] and [HT2] are exact regarding polydisperse acoustics. And so is the semi-kinetic model. We can therefore derive an extension of the formula for polydisperse sprays in a MF context:

k k 0 2 = 1 + k C k 1 -iωτ u k 1 + k C k κ 1-iωτ T k 1 + k C k γκ 1-iωτ T k (21)
where

C k = m k /ρ g .
This formula converges towards a continuous formulation:

k k 0 2 = 1 + R + ρ l ρ g S 3 2 6 √ π n(S) 1 -iωτ u (S) dS 1 + R + ρ l ρg S 3 2 6 √ π κn(S) 1-iωτ T (S) dS 1 + R + ρ l ρg S 3 2 6 √ π γκn(S) 1-iωτ T (S) dS (22)
which is intrinsic and does not require the MF formalism and hypotheses. It is however subjugated to the semi-kinetic hypotheses.

The above relations can be studied in the limit of distributions of very small or very big droplets. We define St ω,min and St ω,max the acoustic Stokes number of the smallest and biggest droplet of the spray respectively. In the particular case of small mass loadings C 1, we recover two limits of attenuation stated in [START_REF] Harrje | Liquid propellant rocket combustion instability[END_REF]:

α =    ωτ u 53 + (γ -1)κωτ T 53 St ω,max 1 1 ωτ u 31 + (γ -1)κ 1 ωτ T 31 St ω,min 1 (23)
which depend each on a specific moment of the droplet size distribution d pq through the relaxation times noted e.g. τ u pq = ρ l (d pq ) 2 /(18ν g ). In general cases yet, approaches based on such global moments will fail as said previously.

Response of some typical sprays

We consider the three size distributions shown in Figure 2, a monodisperse distribution, a lognormal distribution -the form of which is classical in solid propulsion [START_REF] Price | Combustion of metallized propellants[END_REF]-, and a window distribution. These distributions are chosen to have the same average diameter d 31 = 6.53 µm. The exact attenuation and dispersion induced by these sprays on acoustic waves are shown in Figure 4, as computed from the analytical relation of Eq. [START_REF] Doisneau | Eulerian VS Lagrangian simulation of unsteady two-way coupled coalescing two-phase flows in solid propellant combustion[END_REF]. Compared to the classical response of a monodisperse spray, polydisperse sprays have a wider attenuation range but a lower attenuation maximum. Moreover the response of the lognormal spray is more biased for low frequencies than the one of the window spray as a consequence of the lognormal size distribution tail. As a conclusion, polydispersity has an impact on acoustic waves as soon as it has a significant mass of droplets with Stokes numbers separated by at least one decade, which we refer to as "genuine polydispersity".

Conclusions on the linear theory

The linear theory gives the exact solution of acoustics in a two-phase, quasi-uniform medium that is attenuation and dispersion. The theory has been extended to polydispersity. However it does not handle the segregation of droplets by the flow, droplet inhomogeneities and structures having significant effects on acoustic waves. Neither does it handle non linear wave propagation. So these interesting results are inadequate for many applications, especially solid rocket motors, where inhomogeneities and non-linearities are important.

The theory is however mostly useful for code validation purpose: two-phase acoustics brings into play all the couplings between the phases of a two-phase flow. Challenging a code on such a simple case as wave propagation in a quasi-uniform medium allows to evaluate its ability to capture accurately all the gas-droplet couplings.

In the particular context of the Multi-Fluid method, a sectional accuracy study is conducted to assess the requirements on size discretization. These requirements are given hereafter for sprays featuring some typical size distributions, the response of which are also computed with the analytical approach. These responses will also be used as reference solutions in § V.

B. Ability of MF methods to handle polydisperse acoustics

Polydisperse acoustics is completely and exactly captured by the semi-kinetic model. But at the Multi-Fluid level, the droplet size distribution is discretized and a unique velocity is assumed in each size interval so that the acoustic response of a polydisperse spray is only approximated. Convergence is gained if the section width in terms of Stokes number is not too large which is the reciprocal property of the "genuine polydispersity" criterion given in § III. Because OSM and TSM both rely on a unique velocity moment to capture droplet dynamics in the section, they should converge at the same order. Yet TSM provides a better estimation of the average relaxation times, compared to OSM where they are fixed by the section bounds. Especially TSM is remarkably better with low numbers of sections.

Sectional accuracy study

We now assess some aspects of the convergence of the MF method towards the two-phase acoustic response: we focus on the ability of the method to render accurately the shape of the dispersion relation with a low number of sections. There is a wide variety of options to parameterize the MF, e.g. the number of conserved moments and the form of the reconstruction We first consider the window distribution from Figure 2 and we choose OSM with a 1 κ k (r) = cst. The size phase space is discretized with intervals of constant width as regards the variable r, which is a simple and natural first approach. The convergence of the acoustic response with such a discretization is shown in Figure 5 and is satisfactory with as few as 4 sections.

We now consider the lognormal distribution from Figure 2: this distribution is wider in terms of mass repartition (MDF) as shown in Figure 3. We here choose a discretization that guarantees the same amount of mass within each section. In addition we assess the acoustic response of the MF spray with the exact d 31 , computed from the analytic distribution on the section. This discretization strategy is more complex than the previous one and results are given in Figure 6. The acoustic response is well estimated with four sections.

As a conclusion, MF methods are efficient to treat polydisperse two-phase acoustics, even with a low number of sections: the discretization technique must be well-chosen, depending on the width and complexity of the spray's size distribution. Appendix A details the results that are obtained for the response of the lognormal spray, depending on the discretization technique.

Stiffness quantification of acoustics with MF

As explained in § II B, the MF method narrows the time scale spectrum. The previous sectional accuracy study shows that a rather coarse discretization can capture the features of polydisperse acoustics. If we refer to the "genuine polydispersity" criterion, a decade or more -depending on the size distribution smoothness-can be spared on the smallest time scale through the size averaging effect of an appropriate size discretization. The MF method intrinsically reduces the stiffness of two-phase Eulerian modeling.

Conclusion on polydisperse acoustics with MF

As a conclusion, the linear acoustic theory, here extended to polydispersity, gives the key to polydisperse acoustics and the way to handle it with MF methods. Yet, the non rotational linear approach is limited to cases with uniform initial droplet repartition. On the contrary, the Multi-Fluid method -which can capture accurately, with few sections, polydisperse acoustics and which reduces stiffness-is spatially resolved and no linearization is performed so it naturally captures the physics of rotational flows and non-linear waves. The Multi-Fluid method is therefore well adapted for any acoustic or dynamic study.

IV. Numerical methods and strategy for two-phase flows

In this section, we define the common numerical features required to solve two-phase flows with a fully Eulerian strategy. A classical time integration technique referred to as BWC is detailed, that easily couples the two-phases but has two drawbacks. First it fails for strong couplings i.e. for high mass loadings. Second, it requires a numerical time step to integrate all the disperse phase equations that is shorter than the smallest relaxation time, therefore yielding a high computational cost when low inertia droplets are present. So a new strategy is designed to remain accurate and robust for high mass loadings and for low inertia droplets. It is designed to require numerical time steps, only driven by the gas CFL constraints instead of spray stiffness. It is based on a two-stage operator splitting to decouple droplet convective phenomena from two-phase acoustics. We refer to this strategy as the acoustic-convection splitting (ACS). We finally describe the implementation of ACS in a research code and in the CEDRE code.

A. Common numerical features of two-phase flow approaches 1. Space discretization

In fluid mechanics, finite volume approaches are classical as they ensure a good mass conservation and an intuitive form for the source terms. Non-structured meshes are mandatory in industrial applications but structured meshes are used in some codes, for a wider range of schemes is available, achieving some specific properties or being more robust and accurate. The numerical focal point of a flow simulation is then the convective part: high order reconstructions of the flux are required to increase accuracy on the gradients, because space repartition is nonetheless of interest but drives many phenomena through source terms of turbulence, reaction, etc. Yet the order is often limited to three or four for stability and compactness reasons.

In laminar regimes, the monokinetic closure [HV1] and [HV2] models the evolution of dilute to moderately dense disperse phases as pressureless gas dynamics (PGD) [START_REF] Bouchut | On zero pressure gas dynamics[END_REF][START_REF] Zel Dovich | Gravitational instability : an approximate theory for large density perturbations[END_REF], with strong gradients and ratios of neighbors' densities, even vacuum zones and δ-shocks, which numerically result in high amounts of mass within a cell. The main issue in this context is to preserve positivity, so little oscillation is admitted and to cope with strong gradients: these issues directly condition robustness. Second order schemes, based on the seminal work of Bouchut [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF], have been developed for the PGD of dilute to moderately dense sprays on structured meshes [START_REF] De Chaisemartin | Polydisperse evaporating spray turbulent dispersion : Eulerian model and numerical simulation[END_REF][START_REF] De Chaisemartin | Eulerian models for turbulent spray combustion with polydispersity and droplet crossing[END_REF]. For industrial applications, the use of non-structured meshes is required for sprays [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF], which is possible in the CEDRE code with fully generic meshes [START_REF] Courbet | Space discretization methods[END_REF]. The disperse phase schemes, based on MUSCL reconstructions inspired from the gas solver schemes [START_REF] Murrone | Numerical modeling of dispersed two-phase flows[END_REF], reach order two [START_REF] Le Touze | Numerical methods on unstructured meshes for the simulation of sprays within liquid rocket engines[END_REF].

Polykinetic closures always yield a particle pressure term that can either render deterministic (resolved) crossings [START_REF] Doisneau | Eulerian modeling and simulation of small scale trajectory crossing and coalescence for moderate-Stokes-number spray flows[END_REF][START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial particle-laden flows[END_REF], the so-called Particle Trajectory Crossing (PTC) or unknown sub-grid scale crossings [START_REF] Fevrier | Partitioning of particle velocities in a gas-solid turbulent flow into a continuous field and a spatially-uncorrelated random distribution: theoretical formalism and numerical study[END_REF][START_REF] Kah | Eulerian quadrature-based moment models for dilute polydisperse evaporating sprays[END_REF][START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF], yielded by unresolved turbulence in LES for instance, though the validity of these closures has not been fully demonstrated yet. Recent works have been carried out to qualify and compare some methods [START_REF] Vié | On the direct numerical simulation of moderate-Stokesnumber turbulent particulate flows using algebraic-closure-based and kinetic-based moment methods[END_REF]. The structure of the convective part of the disperse phase equations is then similar to Euler equations but the pressure is much weaker so peculiarities close to those of PGD occur and the schemes of standard Gas Dynamics (GD) are often unstable. This situation is referred to as Hypercompressible Gas Dynamics (HGD). For instance, the Mesoscopic Eulerian Formalism (MEF) [START_REF] Fevrier | Partitioning of particle velocities in a gas-solid turbulent flow into a continuous field and a spatially-uncorrelated random distribution: theoretical formalism and numerical study[END_REF], that is implemented in the AVBP code [START_REF] Kaufmann | Vers la simulation des grandes échelles en formulation Euler-Euler des écoulements réactifs diphasiques[END_REF] and has been used for SRM simulation [START_REF] Simoes | Modélisation eulérienne de la phase dispersée dans les moteurs à propergol solide, avec prise en compte de la pression particulaire[END_REF], renders unresolved turbulence through particle pressure and viscosity. The code uses a cell-vertex formulation and Two-Step Taylor Galerkin approaches [START_REF] Colin | Development of high-order Taylor-Galerkin schemes for LES[END_REF] that ensure a third or fourth order space-time convergence on cartesian discretizations. These schemes are however unstable and require artificial viscosity that is applied with a sensor designed for the disperse phase [START_REF] Martinez | Simulation aux grandes échelles de l'injection de carburant liquide dans les moteurs à combustion interne[END_REF][START_REF] Vié | Simulation aux grandes échelles d'écoulements diphasiques turbulents à phase liquide dispersée[END_REF].

For monokinetic closures however, the numerical peculiarities of PGD can bring in unphysical singularities [START_REF] De Chaisemartin | Polydisperse evaporating spray turbulent dispersion : Eulerian model and numerical simulation[END_REF][START_REF] De Chaisemartin | Eulerian models for turbulent spray combustion with polydispersity and droplet crossing[END_REF][START_REF] Kah | Eulerian quadrature-based moment models for dilute polydisperse evaporating sprays[END_REF]: in a two-way coupling context, these singularities would generate a singular retrocoupling towards the gas -or an arbitrarily strong retrocoupling in a spatially discretized context-as highlighted in [START_REF] Doisneau | Eulerian VS Lagrangian simulation of unsteady two-way coupled coalescing two-phase flows in solid propellant combustion[END_REF]. This requires at least robust schemes if the impact of the singularities remains limited, whereas additional modeling is needed if the two-way coupling is significant in such zones. To cope with singularities, polykinetic closures can indeed be introduced [START_REF] Doisneau | Eulerian modeling and simulation of small scale trajectory crossing and coalescence for moderate-Stokes-number spray flows[END_REF][START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial particle-laden flows[END_REF][START_REF] Vié | On the direct numerical simulation of moderate-Stokesnumber turbulent particulate flows using algebraic-closure-based and kinetic-based moment methods[END_REF]: these methods are promising for two-way coupling of inertial (crossing) sprays but the effect of the disperse phase on the gas has never been assessed yet. So we develop a method that ensures two-way coupling between the gas and a disperse phase, here modeled with PGD though discrepancies can be expected when δ-shocks form. Our numerical approach has however been designed to remain relevant for new modelings, e.g. polykinetic closures.

Time integration

The first approach that achieves two-way coupling is a simultaneous time integration of the droplet and gas operators, space fluxes being considered as time dependent source terms. When using explicit time integrations, the time step must be smaller than the smallest characteristic time. As said in § II B, the two-phase system is stiff so that there are several orders of magnitude between this small scale and the biggest scales so the approach is not cost-effective. An overall implicitation can improve stability at higher time steps, and so reduce the cost but the numerical task is harsh because the number of variables and their cross-dependencies increases the scheme's complexity. Moreover, the overall time integration must meet various requirements linked to the discretization of the operators: the time integration should be performed at the same order than space discretization to achieve time consistency [START_REF] Allaire | Analyse numérique et optimisation: Une introduction à la modélisation mathématique et à la simulation numérique[END_REF][START_REF] Duarte | Adaptive numerical methods in time and space for the simulation of multi-scale reaction fronts[END_REF] and high order space schemes often need high order time integration for stability reasons e.g. WENO schemes [START_REF] Wang | Linear instability of the fifth-order WENO method[END_REF]. The phase transport operators must achieve positivity on density, maximum principle on velocity and temperature etc.

Operator splitting techniques can be used for unsteady problems with multiple time scales. The principle of splitting techniques is to define operators in a system of equations and to integrate them consecutively. The benefit is twofold compared to an all-at-once integration. First the numerical methods used for an operator can meet their theoretical specifications (stability, positivity, maximum principle etc.). Second time steps that are internal to each numerical operator are subjected to the only numerical stability conditions of their respective numerical schemes. The major drawback is that the splitting time step ∆t c has to be chosen specifically for the configuration that is solved, keeping in mind that the splitting error should remain limited. Adaptive time stepping could also be used in order to achieve a parametrized and optimal level of error [START_REF] Duarte | Adaptive numerical methods in time and space for the simulation of multi-scale reaction fronts[END_REF][START_REF] Duarte | A new numerical strategy with space-time adaptivity and error control for multi-scale streamer discharge simulations[END_REF]. This overall time step has to be based on physical time scales of the coupling phenomena. Splitting the operators often yields pure convection problems and source term problems. The first ones can be integrated classically with an order compliant method and the CFL requirements apply. The sources are mostly local problems (except for radiation) and they are eventually stiff.

We now consider two approaches that are more or less based on operator splitting. The first one is classical and widely used [START_REF] Murrone | Numerical modeling of dispersed two-phase flows[END_REF] but not efficient for strong two-way coupling. The second one is here developed to yield an efficient and accurate method for all levels of coupling and all sizes of droplets.

B. The Balanced Weak Coupling: a classical but limited approach 1. Definition A convenient technique for small mass loadings of droplets is referred to as the Balanced Weak Coupling (BWC). On a time step ∆t BWC , it consists in computing the droplet phase alone, assuming that the gas is frozen as in one way coupling approaches. But source terms are then reinjected in the gas equations to ensure total momentum and energy conservation on the time step. In this sense, the technique can be considered as a splitting technique with modified source terms. The key parameter of this technique is the ratio of the integration time step to the droplet relaxation time τ min :

K BWC = ∆t BWC τ min (24)
which drives the quality of the coupling.

Study in 0D

For 0D cases, this method ensures correct convergence to the equilibrium values whatever K BWC as shown in Figure 7. For small mass loadings C, the frozen gas hypothesis is correct on all time steps as shown in the top left graph in Figure 7 because droplets weakly influence the gas so u g The 0D study illustrated in Figure 7 shows that BWC oscillates and converges difficultly for high mass loading 0D cases. But in higher dimension cases, transport in physical space is performed at each time step with the computed transient velocities, which are poorly evaluated. So BWC is proven to be ill-suited for moderately dense disperse phase flows. Moreover, high mass loading zones can also be created from dilute average conditions e.g. at the edge of vortices so that BWC is also inaccurate for sprays that are dilute in an average sense as soon as the flow is likely to generate droplet segregation and therefore zones of strong mass loading.

Conclusion on the BWC

It is shown that because of stiffness and coupling, a first approach for solving System (8) would be to integrate the whole set of equations with a time step shorter than the smallest time scale of the spray flow [START_REF] Billet | Impact of Volume Viscosity on a Shock/Hydrogen Bubble interaction[END_REF]. For big droplets, this approach introduces no additional constraints compared to a single-phase case since the fastest time scale is the one of the gas. With an explicit time integration for instance, recommended for unsteady solution determination and parallel computation, the integration time step has anyway to be smaller than τ g for CFL stability requirement. With an implicit integration, CFL should remain moderate, for accuracy purpose. But for small droplets, the integration time step can be severely constrained by τ min which decreases like the square of the mean droplet diameter in the section. So we suggest in the following a new numerical strategy where the time step can be chosen to reduce computational cost while keeping robustness and accuracy.

C. The Acoustic-Convection Splitting: a sound and adaptable strategy A new strategy is designed, based on a two-stage operator splitting to decouple droplet convective phenomena from two-phase acoustics. We aim at increasing computational efficiency while preserving two-phase acoustics features such as dissipation and solving accurately for convection. We refer to this technique as the acoustic-convection splitting (ACS).

We consider System (8) as composed of transport operators T g and the T k 's and a relaxation operator R. These operators apply on U(t) ∈ R 3Nsec+3 where:

U(t) = [U g , (U k ) k=1,Nsec ] t (t) = ρ g , ρ g u g , ρ g e g , (m k , m k u k , m k h k ) k=1,Nsec t (t). (25)
The gas transport operator T g accounts for full-Mach gaseous dynamics (convection and acoustics) and the N sec transport operators T k account for the droplet convection of the sections. They do not modify the other phases' variables so that transport occurs independently for all "fluids". The gas transport operator reads:

T g    ∂ t ρ g + ∂ x •(ρ g u g ) = 0 ∂ t (ρ g u g ) + ∂ x •(ρ g u g ⊗u g ) = -∂ x p ∂ t (ρ g e g ) + ∂ x •(ρ g e g u g ) = -p∂ x •u g (26)
and the transport operator for each section k reads:

T k    ∂ t m k + ∂ x •(m k u k ) = 0 ∂ t (m k u k ) + ∂ x •(m k u k ⊗u k ) = 0 ∂ t (m k h k ) + ∂ x •(m k h k u k ) = 0 (27)
The relaxation operator, which yields all the coupling, reads:

R                            ∂ t ρ g = 0 ∂ t (ρ g u g ) = - Nsec k=1 m k F k ∂ t (ρ g e g ) = - Nsec k=1 m k H k + Nsec k=1 m k F k (u g -u k ) ∂ t m k = 0 ∂ t (m k u k ) = m k F k ∂ t (m k h k ) = m k H k    k = 1, N sec (28)
This relaxation ensures the coupling between all the "fluids" but this coupling is local in space since convection is not considered on the time step. It has two major features:

• conservation of local mass, momentum and total energy of the spray.

• evolution towards an equilibrium state described by formula [START_REF] Courbet | Space discretization methods[END_REF].

So we can exactly rewrite System (8) between t and t + ∆t c :

U(t + ∆t c ) = R + T g + Nsec k=1 T k (∆t c )U(t). ( 29
)

Droplet phase convection splitting

We now introduce the structure of the first stage of ACS. It can be performed according to two different types of Lie splittings:

U Lie1 (t + ∆t c ) = R + T g (∆t c ) Nsec k=1 T k (∆t c )U(t) (30)
or

U Lie2 (t + ∆t c ) = Nsec k=1 T k (∆t c ) R + T g (∆t c )U(t) (31)
or according to a Strang splitting:

U Strang (t + ∆t c ) = R + T g ∆t c 2 Nsec k=1 T k (∆t c ) R + T g ∆t c 2 U(t). (32)
If each operator is solved with a numerical strategy that is at least second order in time, the suggested Strang splitting [START_REF] Strang | On the construction and comparison of difference schemes[END_REF] ensures a second order convergence with ∆t c as long as the solution is regular. Large scale droplet structures can have a strong impact on the overall flow so they have to be accurately captured by the numerical method: the droplet phase transport scheme has to have a low level of dissipation. Excessively refined time steps can introduce numerical diffusion. We suggest a two-stage splitting based on the following time scale hierarchy: droplet phase convection generally yields the slowest scales while the fastest scales are gas transport scales since they propagate acoustic waves. The relaxation scales are to be linked to the acoustic scales in order to account for two-phase acoustics.

In the following, we use the Strang splitting with the stiffer operator, i.e. acoustics, outside in order to end a cycle close to the acoustics equilibrium manifold [START_REF] Descombes | Operator splitting for nonlinear reaction-diffusion systems with an entropic structure: singular perturbation and order reduction[END_REF]. Moreover, the acoustics operator R + T g is the only one that features an equilibrium manifold so that the corresponding operator should always be applied last for accuracy purpose. The convective splitting time step is chosen to allow droplet phase convection to occur at its physical time scale which matches numerical criterion of the transport scheme:

∆t c ∼ τ c k . (33)
In addition, the droplet transport solver is called less often than in a fully coupled strategy which spares CPU resources.

Acoustic splitting

We now tackle the operator R + T g , which contains the physics of the so-called two-phase acoustics. We approximate these terms with a Strang splitting on a time step ∆t a . This defines the "acoustic" operator A:

A(∆t a ) = R ∆t a 2 T g (∆t a )R ∆t a 2 . ( 34
)
The two-phase acoustic operator A performs the gas transport and the relaxation, assuming that the droplets encounter acoustic effects locally and do not need to be convected at this time scale. This is proven to be satisfactory in the following.

Time step prescription strategy

Let us now consider the requirements on the splitting time step ∆t a , which ensures the level of coupling between the gas and the droplets and is therefore referred to as the coupling time step. It is set to follow the largest time scale (corresponding to the slowest phenomenon) between gas-droplet relaxation and gas acoustics:

∆t a = max K p τ min ; K g τ g ( 35
)
where K g and K p are security coefficients, which values are addressed later. It is indeed not necessary to ensure a coupling at so large a frequency that either the gas or the smallest droplets cannot react. So τ min is used for sprays involving only big droplets since it is larger but τ g is used as soon as small droplets are present. In both cases, CPU resource can be spared because the coupling time step is never smaller than the one of the gas and it can become very large when droplets have a large inertia.

There is however an upper bound for the coupling time step, imposed by the wave period τ ω for stability reasons. The relaxation operator should capture the average value of the acoustic wave so a minimum number of coupling instants are imposed in a period:

∆t a < τ ω 4 (36)
which avoids offsetting the equilibrium value of the two-phase system. The evolution of the prescription for ∆t a depending on the smallest droplet size is illustrated in Figure 8 with typical values for K g and K p .

Summary of the features of ACS

The ACS strategy operates within a time step defined with the smallest droplet convective scale ∆t c :

U ACS (t + ∆t c ) = A ∆t c 2 Nsec k=1 T k (∆t c )A ∆t c 2 U(t) (37)
where we decompose the acoustic operator

A(∆t c ) = [A(∆t a )]
∆tc/∆ta [START_REF] Fréret | On the role of preferential segregation in flame dynamics in polydisperse evaporating sprays[END_REF] in order to solve it with respect to a two-phase acoustic time scale ∆t a :

A(∆t a ) = R ∆t a 2 T g (∆t a )R ∆t a 2 . ( 39
)
This method is appropriate for smooth gaseous fields but discontinuous states can introduce strong variations in short times that are not compatible with the time steps suggested splitting techniques. This requires to check that the droplet velocities have not varied too much during A, which for instance can be done for smooth fields with an embedded error control method [START_REF] Descombes | Adaptive time splitting method for multi-scale evolutionary partial differential equations[END_REF][START_REF] Duarte | Adaptive numerical methods in time and space for the simulation of multi-scale reaction fronts[END_REF][START_REF] Duarte | A new numerical strategy with space-time adaptivity and error control for multi-scale streamer discharge simulations[END_REF]. Such adaptive methods are not discussed since we work with smooth acoustic cases. In general, a strategy based on operator splitting and relatively large time steps is valid in the absence of discontinuities e.g. even for choked nozzles. Furthermore to treat discontinuities such as shocks, time steps smaller than all the scales are required [START_REF] Billet | Impact of Volume Viscosity on a Shock/Hydrogen Bubble interaction[END_REF]. In these cases, which are out of the scope of this paper, the mesoscopic treatment of multiphase flows with shocks requires additional modeling.

The ACS strategy has been designed to decouple the stability and accuracy issues for unsteady strongly coupled two-phase problems. It yields the following advantages:

• the overall coupling time step ∆t c is prescribed only regarding accuracy on droplet transport;

• the acoustic time step ∆t a is prescribed only regarding two-phase dissipation/dispersion accuracy;

• all the operators are implemented as modules having their own adapted time integration methods.

Finally, the ACS strategy is flexible enough to consider other physical phenomena through operator splitting such as coalescence, break-up or evaporation. Coalescence has been for instance implemented and feasibility has been assessed on a scale 1 Ariane 5 motor [START_REF] Doisneau | Eulerian VS Lagrangian simulation of unsteady two-way coupled coalescing two-phase flows in solid propellant combustion[END_REF].

D. Implementation of ACS

The ACS is implemented in a research code referred to as SAP1 in order to validate the strategy and evaluate its performance on a two-phase acoustic test case. The ACS is also implemented in CEDRE to demonstrate the feasibility of the method in complex cases.

The SAP1 code is a 1D code that solves a disperse phase, with the one size moment Eulerian MF approach, coupled to a compressible gas. The code has been designed to test splitting errors alone. It therefore features high quality numerical methods that are very well-suited to the requirements of acoustics, no matter their cost. The gas transport scheme is based on WENO5 space fluxes integrated with a third order explicit Runge-Kutta method (RK3) to generate very low numerical diffusion. This scheme can handle shock waves but this feature is not used in the present paper. For the droplet transport, which is decoupled as regards each section, a second order in space and time Bouchut scheme is chosen, according to strategy for PGD described in § IV A 1. The Bouchut scheme ideally operates at CFLs that are close to one to minimize numerical diffusion. Thanks to the two-stage splitting of the ACS strategy, we can increase ∆t c to solve the section convection with a CFL k chosen to limit numerical diffusion in physical space, as shown in § V C. The relaxation step is based on RADAU5, which is a fixed error implicit RK5 solver for ODEs with adaptive time step [START_REF] Hairer | Solving ordinary differential equations[END_REF]. This method is A-stable and L-stable and is thus specifically designed for stiff ODEs. Besides, the Geneva group around E. Hairer has provided a very efficient implementation of such implicit RK methods, which we use here.

The CEDRE code can solve a disperse phase, with a two size moment Eulerian MF method, coupled with a compressible and eventually viscous gas. But the code is designed for research and industrial simulations of complex multi-physics problems so the numerical methods are optimized regarding cost and accuracy. The transport schemes for both gas and droplets are based on second order MUSCL space fluxes inspired from the Bouchut scheme, integrated with a third order implicit Runge-Kutta method (RKI3) to extend the stability domain above CFL= 1. As for relaxation numerical methods, though an analytical solution exists for monodisperse sprays in Stokes regime, we will use a more general numerical method for the relaxation step in order to easily generalize to polydisperse sprays or high Re p regimes. We insist on the fact that, even with multiple sections, operator R is a system of ODEs and can be easily solved, though the different droplet size characteristic time scales can be very different, therefore yielding a stiff problem. So we use in the following a simple numerical method which is sufficient when K p is small enough.

V. Validation and limits of ACS

Accounting for two-phase acoustics is a demanding ability as it requires to capture fully coupled, unsteady phenomena. So it is a good criterion to evaluate the quality of the coupling that is achieved between the gas and the disperse phase. We therefore use the classical acoustic theory [START_REF] De Chaisemartin | Eulerian models for turbulent spray combustion with polydispersity and droplet crossing[END_REF] as well as the polydisperse dispersion relation derived in § III A 3 to assess the performance of the suggested numerical strategy regarding the couplings, as explained in § III A 4. We consider in this section two test cases for ACS: the first one is based on a damped propagation and shows the ability of the method to reproduce the dissipation and dispersion of an acoustic wave according to the formulae of § III. It also allows to assess the limits of the security coefficient. The second one is a sine profile -therefore propagating like the previous wave-but in a moving frame. This second case proves the validity of decoupling droplet transport from two-phase acoustics, by reaching the same level of accuracy in a convective context. The validation is then performed on the industrial-oriented code CEDRE to show that the ACS strategy is efficient for general cases and realistic geometries.

A. Acoustic case description

Let us now describe the acoustic wave and the numerical parameters used for both configurations. We consider a 1D acoustic wave of a given frequency of 1000 Hz, with a small relative amplitude of 10 -4 so that its propagation remains in the linear regime and Stokes' laws are valid. In SRM conditions of 3600 K and 50 bar, the speed of sound in a burnt gas mixture can reach c 0 = 1275 m/s. The mixture has a dynamic viscosity ν g = 8.855.10 -5 Pa.s. The droplet density is 1766 kg/m 3 and the mass loading is set to a high value of C = 0.394. As an illustration, the instantaneous 1D fields for both phases are given in Figure 9, where it can be seen that the oscillations are damped as the wave penetrates the two-phase medium.

The domain length for both codes is set to 25.51 m which allows 20 periods in average for the propagative case and the 2000 cell mesh has a spatial resolution of 100 points per wavelength (ppw). Using 100 ppw ensures for SAP1 as well as for CEDRE numerical attenuation of the acoustic waves that is lower than the two-phase attenuation for all cases and second to increase the accuracy on the post-processing of attenuation and dispersion. We insist on the fact that the post-processing is done carefully since attenuation and dispersion are very sensitive: a simple approach based on local extremum tracking can be used to estimate amplitude decrease and wavelength, respectively linked to α and β but errors due to discretization must be minimized e.g. by interpolation. A fitting of the best decaying sinus to the signal, e.g. with a least square method, can be used to increase accuracy on the coefficients.

For monodisperse distributions, the response can be plotted versus the acoustic Stokes number and we can change the droplet diameter in order to sweep the acoustic Stokes number range without changing the frequency and the mesh. But for polydisperse distributions, the distribution is fixed and the mesh is changed; the response is plotted versus the excitation wave frequency.

B. Study of the ACS splitting error with SAP1

The quality of the two-phase acoustics depends on the splitting time steps that are chosen for ACS (cf. § IV C). The level of accuracy will depend on the choice of K g if there are small droplets or K p otherwise. We hereafter make an empirical study of the splitting error with SAP1. The chosen test case is a pure acoustics one where the two-phase dispersion relation (( 21)) applies. •: Kg = 1; ∆: Kg = 0.9; 2: Kg = 0.5; ∇: Kg = 0.2.

Splitting error for small droplets

We show in Figure 10 the attenuation α and dispersion β for different values of the K g coefficient. First, the computation can give arbitrary results when K g > 1 for very small droplets as shown on both curves of Figure 10, this corresponds to the small droplet "no go area" illustrated in Figure 8. Second, an erroneous constant attenuation is yielded below a size depending on the value of K g , as shown on the left curve in Figure 10. This overdissipation seems acceptable as soon as K g < 0.9 and the analytic result is satisfactorily recovered for K g = 0.2. An analytical study of the acoustic splitting error which is not detailed here confirms that this overdissipation plateau for small droplets is naturally yielded by the splitting strategy itself.

Splitting error for big droplets

The margin on the relaxation time, K p , depends on the numerical method used to compute the relaxation operator. We have conducted a numerical study with a simple relaxation approach that suggests to set K p < 0.25 to ensure an accuracy on dissipation and dispersion better than 1% but the details are not provided for the sake of legibility. We consider that the computation can fail for most methods when K p > 1 which corresponds to the "no go area" illustrated in Figure 8. Finally the prescription of K p is not a key issue since it only allows to spare calling the relaxation module, which is not costly.

Capturing two-phase acoustics

Results for SAP1 with typical security coefficients (K p = K g = 0.5) are compared to the analytic solution. The results are given together with those of a further study in Figure 14. There is an excellent agreement of SAP1 results. The high order numerical methods in SAP1 allow an error smaller than 1% on both α and β which is very good, regarding the sensitivity of these coefficients.

We then consider a polydisperse case: the window size distribution that is presented in Figure 2. Contrary to the monodisperse case where we could sweep the droplet size, we can only sweep a frequency range, so a mesh has to be specifically created for each frequency, for the number of periods in the domain to remain constant. Moreover the test case is too costly to perform with the BWC approach because of the presence of very small droplets introducing small time scales. The polydisperse case is performed only for SAP1 and CEDRE (ACS) and the results are given in Figure 15, using 10 sections which is enough to capture polydisperse acoustics with good accuracy as discussed in § III B 1. The MF method handles naturally polydisperse cases as linear combinations of monodisperse cases so this test case is just a formal check of the code, the numerical performances being already evaluated in the monodisperse case. The results on dissipation and dispersion for SAP1 and CEDRE are correct, which validates the approach for polydisperse sprays too. These test cases prove the efficiency and accuracy of the splitting technique to treat the acoustic operator A. This approach is more general than ACS itself, and it does not require the splitting of the droplet convection part to work. The authors assume that the previous results could have been obtained without separating the droplet convection. The aspects of treating droplet convection separately from gas transport are evaluated hereafter.

C. Capturing convection with ACS

As exposed previously, the convective splitting allows the use of time steps that are adapted to the droplet convection operator. But droplet convection has small impact in our previous test case since the wave amplitude is small: droplet velocity remains below 10 -4 c 0 . We now test the ACS strategy on a case featuring significant droplet convection.

Benefit of ACS for droplet phase convection schemes

As stated previously, the droplet transport scheme can have an ill-suited time step imposed by acoustics if a coupled strategy is chosen. This generates an important numerical diffusion with most transport schemes. With ACS, the droplet convection time step can be taken as close to one as desired which solves the droplet convection accuracy issue.

We perform a numerical test case with SAP1 on the evolution of a window-shaped packet of droplets when purely convected. We compare its convection by a second order scheme either when the time step is conditioned by acoustics or when it can be chosen freely thanks to ACS. The results are presented in Figure 11. The packet that is convected at CFL k = 0.081 -due to a time step that is constrained by gas physics-is smeared: the Bouchut transport scheme is ill-suited for such low CFL numbers and it operates at best with a CFL close to one, which is a classical requirement for this class of kinetic schemes [START_REF] De Chaisemartin | Polydisperse evaporating spray turbulent dispersion : Eulerian model and numerical simulation[END_REF][START_REF] Massot | Eulerian Multi-Fluid models: modeling and numerical methods[END_REF]. The packet convected with the same scheme but in ACS time integration framework is much better preserved, simply because the scheme now operates at CFL k = 0.94, at the user's choice. With ACS indeed, the best CFL can be chosen for transport, as the strategy decouples the operators.

Accuracy of ACS on acoustics with droplet convection

We now consider an acoustic disturbance propagating in a spray with ensemble motion. This disturbance has the same characteristics as the previous propagating wave but it has a shorter space-time extension. Convection occurs for both droplets and gas which is equivalent to changing the reference frame. The test case is thus analytically trivial and yields the same dispersion relation as Eq. ( 19). Yet for numerical codes, the space translation corresponding to convection is far from trivial as shown previously in Figure 11. The acoustic disturbance is a 3-period sine profile that propagates in a spray with ensemble motion at velocity u g = u p = 150m/s. We choose a 3-period profile to ease result post-processing. The initial conditions for the gas are shown in Figure 12 while the spray is initially uniform.

There is no need to study polydisperse convection because sections are naturally decoupled regarding transport in physical space. So the test case is monodisperse and the results are plotted versus the acoustic Stokes number. The test case is not performed with CEDRE because of the complexity of the initial condition. Figure 13 presents the results for SAP1 with ACS only for dissipation, since no post-processing has been developed by the authors for dispersion for complexity reasons. The results on dissipation that are presented match very well with the non-moving analytical attenuation, which validates the principle of two-stage splitting technique introduced with ACS.

D. Evaluation of ACS in an industrial code

We now evaluate the performances of ACS with an industrial-oriented code: CEDRE. The numerical methods are designed to be less costly and to work on non-structured meshes so they introduce small errors in two-phase acoustics. The difference in the number of moments in size between SAP1 and CEDRE has no incidence on the following monodisperse cases and small incidence on the polydisperse cases for the reasons given in § III B. Finally, the conditions of the acoustic studies make the comparison possible between SAP1 and CEDRE. Results on acoustic wave attenuation and dispersion for the propagative case (Figure 9) in a monodisperse medium with CEDRE with BWC and CEDRE with ACS are compared to SAP1 and analytical results in Figure 14. The ACS in CEDRE yields very good results so its acoustic part is validated. The BWC in CEDRE doesn't yield satisfactory results below a given acoustic Stokes number. The limit is in this case St ω < 2, this particular value being linked to the integration time step that is chosen for our test case: the range of over-dissipative droplets corresponds exactly to the small droplets. So the BWC is over-dissipative when the droplet relaxation time is close to the integration time step or smaller, as expected regarding the errors on instantaneous velocity differences in the 0D test case given in Figure 7.

Results on acoustic wave attenuation and dispersion in a polydisperse medium with CEDRE with ACS are compared to SAP1 and analytical results in Figure 15. The distribution that is used is the window size distribution that is presented in Figure 2. The ACS in CEDRE still yields very good results which is expected as explained previously. Test cases have been performed down to 5 Hz, which is a very low frequency that requires many iterations, considering the wide scale separation between the wave pulsation ω and the chosen time step.

VI. Polydisperse two-way coupled Solid Rocket Motor simulation

In this final section, we validate the ACS time strategy and we illustrate its robustness and accuracy on a 2D unsteady polydisperse moderately dense case with a complex geometry and hydrodynamic structures. This test case consists in a SRM seeded with polydisperse particles ranging from large, slightly inertial particles to very small particles behaving like tracers and therefore belonging to the "Equivalent Gas" assumption validity domain. This is a tough case with various time scales linked to flow scales and droplet polydispersity. The study of the space discretization errors is beyond the scope of the present contribution, accordingly to the methodology used for similar studies on splitting methods [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multi resolution and dedicated high order implicit/explicit time integrators[END_REF]. In a number of applications indeed, the maximal level of space refinement is set by the computational capabilities so we investigate the new numerical strategy solely for the resolution of a semi-discretized problem on a fixed, "most refined" mesh: in our case, this mesh is far from that of a DNS which is intractable anyway, but it is still a mesh that is used in practical SRM simulations. Once the space (mesh) and size (MF) discretizations are chosen, the time resolution is to be done, for a stiff dynamical system associated to the semi-discretized problem. The abilities of BWC and ACS to deal with the multi-scale features of the problem are assessed. The methodology to prove the time accuracy of the method consists in a comparison between ACS and BWC: the latter is indeed a validated approach and it is considered with a time step belonging to the domain of its convergence so BWC can be considered as a reference.

Considering a time step adapted to the fastest flow time scales i.e. the supersonic nozzle flow, we perform two types of simulations: validation and efficiency test cases. In the validation test case, used as a reference we get rid of the stiffness of very small droplets with the MF-PEG method -since they feature St 1 everywhere-and we prove ACS to match BWC results. In the efficiency test case, the smaller particles are fully and "explicitly" solved within a fourth section by a MF approach -which is therefore stiffer-and ACS retrieves the previous results while BWC departs. Another efficiency test case is finally performed: a MF simulation with even smaller (nano) droplets allows ACS to retrieve the reference results, as expected since very small and nano particles both behave like tracers. On the same nano simulation with the same time step, BWC crashes and requires, as expected, an overall time step smaller than the nano droplet relaxation time. This validates ACS and proves its accuracy and robustness, with larger time steps than BWC for very small (or nano) droplets.

A. Case description and numerical configuration

We consider a SRM case featuring a coupling between acoustic and hydrodynamic instabilities. The considered motor, called LP10, was fired with solid particles at a temperature that prevents them from encountering size evolution so that their distribution is well known and the two-phase effect can be accurately studied. The flow in the chamber is subject to Parietal Vortex Shedding, a hydrodynamic instability coupled to the acoustic mode of the chamber as discovered in [START_REF] Lupoglazoff | Parietal Vortex Shedding as a cause of instability for long solid propellant motors. Numerical simulations and comparisons with firing tests[END_REF], theorized in [START_REF] Casalis | Spatial instability of planar channel flow with fluid injection through porous walls[END_REF], and actively studied up to now [START_REF] Cai | Two-phase vorticoacoustic flow interactions in solid-propellant rocket motors[END_REF][START_REF] Casalis | Some recent advances in the instabilities occurring in long solid rocket motors[END_REF][START_REF] Chedevergne | Direct Numerical Simulation and biglobal stability investigations of the gaseous motion in solid rocket motors[END_REF][START_REF] Chu | Premixed flame response to acoustic waves in a porous-walled chamber with surface mass injection[END_REF][START_REF] Fabignon | Instabilities and pressure oscillations in solid rocket motors[END_REF]. For this specific LP10 case, a driving effect of inert droplets on instability levels has been observed experimentally [START_REF] Prévost | Overview on the representativity of small scale motor tests[END_REF] and numerically [START_REF] Dupays | Inert condensed phase driving effect of combustion instabilities in Solid Rocket Motor[END_REF]. The strong impact of polydispersity on acoustic levels has been demonstrated with CEDRE [START_REF] Doisneau | Two-way coupled simulation of acoustic waves in polydispersed coalescing two-phase flows : application to Solid Rocket Motor instabilities[END_REF]. No other simulation of this case has been led up to now. The simulation of the LP10 instabilities requires capturing flow and spray spatial structures as well as their coupling, these issues being the ones we aim at tackling all along the present contribution, this with a more efficient and still accurate method.

The LP10 simulation is performed with CEDRE on a 2D axisymmetrical representation of an LP10 motor fired during an experimental campaign. A multi-species Navier-Stokes formulation is used, for which the generalization of the splitting strategy can be done. Efficient splitting strategies are indeed qualified for convection-reaction-diffusion equations in any regimes in [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multi resolution and dedicated high order implicit/explicit time integrators[END_REF]; but in our case diffusion is slower than convection since the Reynolds numbers are high and the boundary layers are coarsely meshed and no reactions occur so we choose to time-integrate the viscous operator simultaneously to the carrier gas convection. This approach is satisfactory in our particular regime and it is convenient to implement as it complies to the legacy of the CEDRE code. With the present mesh, representing the geometry of the propellant grain, assumed steady at a given time of the firing, the instability spectrum has well-defined modes based on a monochromatic fundamental. The simulation takes place at a time when the most intense instabilities are observed. The simulated time is long enough to perform a harmonic study of the instabilities. At this time, the forward cylindrical grain has totally burnt so that no flow comes from the front end. The simulation takes place on the deformed-structured 27, 000 cell mesh already used in previous studies [START_REF] Doisneau | Two-way coupled simulation of acoustic waves in polydispersed coalescing two-phase flows : application to Solid Rocket Motor instabilities[END_REF] and shown in Figure 16. No propellant combustion is modeled: gas and droplets are injected at a fixed surface flow rate of 16.01 kg/m 2 /s and a fixed temperature of 2255 K from the two walls representing the two propellant segments. The boundary condition ensures that the normal pressure gradient is null, which determines the gas density and velocity. Particles are chemically inert, non coalescing since they are solid in the experimental firing, and they are at equilibrium with the burnt gas flowing from the propellant. The same distribution, based on experimental droplet data, is injected from both segments to render polydispersity: this distribution has a fairly lognormal shape. In the MF method, the size phase space is discretized in four sections with the diameter bounds 0, 3.33, 6.66 and 10 µm, which is enough to capture polydispersity for a lognormal distribution according to the sectional accuracy study performed in § III B 1 on analytical acoustic solutions. We refer to the corresponding droplets as very small VS, small S, medium M and large L droplets. With the MF-PEG method, the size phase space is discretized in three sections with the diameter bounds 3.33, 6.66 and 10 µm while the surface flow rate at the wall corresponding to VS droplets is accounted for with an EG technique. BWC and ACS strategies are considered for MF and MF-PEG cases. Finally, a MF computation is performed with nano N droplets instead of VS droplets. The details on injection and Stokes numbers for the different droplet categories are given in Table 1.

The case is very stiff since (i) the chamber flow is slow while the nozzle flow is supersonic, and (ii) the disperse phase is strongly polydisperse. In the chamber, the characteristic flow times are associated to hydrodynamic instabilities and more precisely vortex revolution times but also to the fundamental acoustic wave period ω. In our case in fact, the two phenomena are locked so that their characteristic times are linked: we can qualify droplet behavior in the chamber with an acoustic Stokes number St ω,k . In the nozzle, the characteristic time is linked to the gas acceleration and is defined as the inverse of an average strain rate = (u g (x ) -u g (x))/(x -x). We consider its maximum value max = 3.10 4 s -1 and the droplet behavior in the nozzle is qualified with the corresponding nozzle Stokes number St ,k . The gas has a viscosity ν g = 7.057.10 -5 Pa.s and the droplets are made of zirconium oxide with ρ l = 6100kg/m 3 . We can sort the problem's characteristic times as given in Table 2. We insist on the fact that the droplet time scales are estimated with dynamic Stokes times τ u k while Schiller-Nauman corrections are used for strong slip velocities and can be linked to shorter values for the time scales. Moreover the thermal times can be smaller, especially with the Ranz-Marshall correction. So the hierarchy given in Table 2 is only indicative.

The two-phase simulation starts from a steady state of the single phase gaseous flow field. With the chosen mesh, the most restrictive CFL time is τ g = 6.10 -7 s. and the splitting time steps are ∆t c = ∆t a = 2.10 -7 s, being equal because of the coupling between acoustic and convective aspects in the supersonic nozzle. This yields the gas security coefficient K g = 1/3. With the MF-PEG method, the smallest section time is τ min = τ u S = 1.2.10 -4 s which yields a droplet security coefficient K p = 1/600. With the MF method, the smallest section time is τ min = 1.7.10 -6 s which yields a droplet security coefficient K p = 1/9. Once again, these security coefficients are based on Stokes time but the relaxation is faster when Re p is large due to drag and heat corrections.

When the droplets are injected, the pressure builds up to a new value of average pressure. We give instantaneous droplet volume fraction fields in Figure 17 for the PEG-MF ACS as an illustration. Their structure matches the one of the BWC results, at a point that differences are insensible to the naked eye, so that the BWC results are not reproduced here. As expected in classical analysis, the three droplet sizes behave differently in the gaseous vortices: the smaller droplets featuring a uniform repartition while more inert droplets from section three are ejected from the vortices. This behavior complies to the droplets' respective acoustic Stokes numbers St ω,k (cf. Table 1) because the vortices are, in this particular case, tuned on the acoustic instability so their revolution time is equal to the wave period.

Finally, the computation with nano droplets N instead of VS yields the same results for ACS, regarding chamber acoustics and nozzle dynamics. On the contrary, the nano droplet N computation crashes with BWC. This shows the improvement on robustness brought by the new ACS strategy when extreme stiffness is introduced by very low inertia droplets.

B. Dynamic study in the nozzle

We perform a quantitative comparison between the four computations MF-PEG and MF on the instantaneous hydrodynamic data along the three cuts in the nozzle given in Figure 16.

We first compare the two MF-PEG methods. Regarding droplet mass concentrations, given in Figure 18, the MF-PEG results match perfectly between BWC and ACS. Along the A-A cut, concentrations increase and decrease after the throat. The B-B and C-C cuts show that the droplet concentrations are maximum on the centerline and decrease as we get closer to the wall. Regarding velocities and temperatures, given in Figure 19, the MF-PEG results also match perfectly between BWC and ACS. Along the A-A cut, the gas velocity increases as the temperature decreases, and the droplets follow partially this behavior, but with differences that increase with the droplet size. We recall that the three resolved sections S, M and L feature St ,k 1 so that this behavior, i.e. wall depletion, velocity and temperature lag, is as characteristic of very inertial droplets in the nozzle. The B-B and C-C cuts also prove ACS and BWC to match perfectly. Finally the ACS reproduces perfectly the BWC results when droplets have a relaxation time that is large compared to the numerical time step i.e. K p 1: so the validity and accuracy of ACS is proven to time integrate space and size semi-discretized SRM problems.

We now compare the two MF methods and we keep a MF-PEG result as a reference. Regarding droplet mass concentrations, given in Figure 20, the ACS results match very well to the reference while BWC ones depart. Along the A-A cut, concentrations increase and decrease after the throat especially for the VS droplets which dilute much more in the diverging part of nozzle. The BWC overestimates S, M, and L concentrations in the diverging part. The B-B and C-C cuts now show that the VS droplet concentration is very uniform in nozzle sections, compared to inertial droplets (S, M, L). Regarding velocities and temperatures, given in Figure 21, the ACS results also match the reference well and the BWC results still depart. Along the A-A cut, the ACS gas velocity and temperature match the reference, contrary to BWC. The VS droplet temperature as resolved by ACS sticks perfectly to the gas and the ACS VS velocity matches well the gas, according to the analysis of its St ,k ; the inertial sections in ACS match perfectly the corresponding sections of the reference, which is computed in a PEG approach: this is natural as soon as the gas fields are similar in MF and MF-PEG. On the contrary, the BWC which has not well resolved the gas yields discrepancies on all the sections. The B-B and C-C retrieve these results: ACS matches the reference perfectly except for the VS velocities at throat which have a slight lag; the BWC has velocity and temperature discrepancies at the output, close to the centerline. Finally the ACS reproduces very well the reference results with a MF method, even with security coefficient K p = 1/9 close to one; the slight difference at the throat is presumably due to the fact that the PEG hypothesis is at its limit of validity. On the contrary, the BWC results for the MF method depart from the PEG reference so BWC is proven to be inaccurate when time steps are close both to droplet relaxation times and flow characteristic times, which is not satisfactory in a perspective of cutting computational costs.

C. Harmonic study

To study the instabilities, data is gathered during the steady regime and a harmonic analysis is performed on three pressure signals, measured at the front, the inter-segment and at the rear of the chamber. The spectra are shown in Figure 22 for the four computations. The structure of the modes is similar and typical of a longitudinal steady mode: harmonics are integer multiples of a unique fundamental wave and they decrease quasi-linearly for sensors 1 and 3. They match very well, with differences on the first two peaks below 0.5%. So the use of ACS or BWC in this case has no influence on the harmonic behavior of the problem. We link this result to the fact that the chosen time step is far below the acoustic time scales of the problem so that BWC or ACS are not challenged. While proofs of ACS accuracy on two-phase acoustics have been given in § V, the ACS strategy is now qualified to capture unsteady two-way coupled hydrodynamics and acoustics.

D. Conclusion on ACS for SRM simulations

The ACS method is proven to be well suited for SRM simulations. It solves accurately the two-way coupled convection problem in the nozzle while BWC needs a smaller time step to meet the same level of accuracy. It is robust and accurate with a splitting time step based on the flow convective scale, whatever scales are introduced by the spray; in comparison, BWC needs to base its time step on the smallest time scale, which may be imposed by the disperse phase. With such a time step, constrained by the supersonic nozzle, and such a discretization, the acoustic in the chamber is over-resolved so no difference is observed between BWC and ACS; but the monochromatic wave studies have proven ACS to be more accurate than BWC on two-phase acoustics too, when time steps are close to the acoustic time. So this allows ACS to spare cost in further SRM configurations.

General conclusion

Theoretical issues of two-phase acoustics have been tackled: an analytical formula for attenuation and dispersion by polydisperse suspensions is derived and matches analytical asymptotic rules from Temkin and Dobbins' theory. The discretization requirements for MF methods to capture twophase acoustics are investigated and the time scales and couplings of the corresponding system are described. A numerical strategy with adaptable cost/accuracy trade off is developed for unsteady polydisperse moderately dense sprays. It is quantitatively studied and validated on acoustic cases. It is available in an advanced CFD platform and a solid rocket motor case proves the feasibility, the accuracy and the robustness of the method.

Since moderately dense sprays also encounter collisions and coalescence, ACS is extended to the corresponding physics by considering an appropriate splitting of the coalescence operator so that simulations of more complete industrial cases can be achieved with CEDRE. The ACS two-way coupling integration strategy is promising to complete the physical phenomena accounted for in the disperse phase such as evaporation, combustion or break-up but also to couple multi-physics approaches such as droplet radiation, given that Eulerian radiation solvers are available, which is the case in CEDRE.

Interdisciplinary Studies (P. Moin/M. Massot), a Ph. D. Grant from ONERA and an internship grant from EM2C.

APPENDIX A: SYSTEMATIC SIZE DISCRETIZATIONS FOR POLYDISPERSE ACOUSTICS

The ability of the MF method to render polydisperse two-phase acoustics with a low number of sections has been presented in § III B 1. The discretization method had to be carefully chosen, especially with a wide size distribution such as the lognormal one presented in Figure 2, which tail corresponds to a significant widening of the mass distribution (see Figure 3). The MF convergence study is here furthered and four types of systematic discretization techniques are assessed on the acoustic response of the lognormal spray:

• the so-called Constant ∆r discretization method defines equal radius intervals (and has been successfully used in § III B 1 for the window distribution),

• the so-called Constant ∆S discretization method defines equal Surface intervals,

• the so-called Geometric discretization method has bounds defined by r k = r 0 R k with R > 0,

• and the so-called Fixed mass ∆r discretization method has its bounds (r k ) k defined to enclose a constant amount of mass in each section k:

m k = r k+1 r k r 3 f (r)dr = m tot /N sec .
The Fixed mass ∆r discretization has been successfully used in § III B 1 for the lognormal distribution.

These four approaches are illustrated in Figure 23 in the case of a 4 section discretization of the lognormal distribution of Figure 2. The Constant ∆S discretization method makes the first section very wide while the three other sections host very little mass; the Constant ∆r discretization method also suffers from the fact that the lognormal distribution is skewed on the right. The Geometric and the Fixed mass ∆r methods discretize more finely the small sizes: they adapt better to the lognormal skewness, the Fixed mass ∆r discretization having a balanced amount of mass in each section by construction.

The convergences of the acoustic responses with the four discretizations for the lognormal spray are given in Figure 25, 24, 26, and 27 respectively. Compliantly to the observations made before about the quality of the discretizations, the constant ∆r and constant ∆S methods yield poorer results than the Geometric and the Fixed mass ∆r methods. Finally, we consider the approximate response with a Fixed mass ∆r discretization and the exact d 31 computed from the lognormal distribution on the section in Figure 6: the accuracy at low numbers of sections is remarkably good.

As a conclusion, the Geometric discretization method and the Fixed mass ∆r discretization method, as well as the use of an accurately reconstructed d 31 strongly improve the accuracy of MF, which provides a satisfactory response with as few as 4 sections.

Fig. 1

 1 Fig.1Coupling of the gas and section equations with their time scales.

Fig. 2 Fig. 3

 23 Fig. 2 Different number density functions (NDF) with the same total bulk density and average diameter d31 -Left: Dirac size distribution (monodisperse); Middle: lognormal size distribution; Right: window size distribution.

Fig. 4

 4 Fig. 4 Attenuation α and dispersion β of a monochromatic acoustic wave in a polydisperse spray -Solid: Dirac size distribution (monodisperse); Dotted: window size distribution; Dashed: lognormal size distribution.

Fig. 5 Fig. 6

 56 Fig. 5 Attenuation α and dispersion β of a monochromatic acoustic wave in a polydisperse window spray -Dotted: MF with 2, 3, 4, and 10 sections with constant ∆r; Solid: reference.

16 CFig. 7

 167 Fig.7BWC computation of a 0D gas-droplet relaxation case for different mass loadings and different K BWC -: gas; •: droplets; Line: exact solution for each phase.

Fig. 8

 8 Fig. 8 Typical acoustics time scales plotted versus droplet diameter -Black: gas time τg; Dotted: typical droplet relaxation time τ min (BWC timestep must be smaller); Red: prescribed acoustic time step ∆ta for ACS with Kg = Kp = 0.5.

m k (kg/m 3 )Fig. 9

 39 Fig. 9 Propagation test case for Stω = 20 with C = 0.106 -Gas variables (thin line) and droplet variables (+) once the wave has completely penetrated.

Fig. 10

 10 Fig. 10 Attenuation and dispersion of a monochromatic wave in a monodisperse spray with C = 0.106 for different coupling time steps (SAP1 code, Kp = 0.25) -Line: analytic; +: Kg = 2; •: Kg = 1; ∆: Kg = 0.9; 2: Kg = 0.5; ∇: Kg = 0.2.

3 )Fig. 11

 311 Fig.11Convection of a droplet structure with a time step driven by acoustics (window packet convected at 30m/s, second order Bouchut scheme, SAP1 code) -+ droplet phase transport at CFL k = 0.081; × droplet phase transport at CFL k = 0.94; -: initial state.

3 )Fig. 12 Fig. 13

 31213 Fig.[START_REF] Casalis | Spatial instability of planar channel flow with fluid injection through porous walls[END_REF] Convective test case in a spray with C = 0.106 -Gas variables at three different times (Initial state: centered 3-period sine profile).

Fig. 14 Fig. 15

 1415 Fig. 14 Attenuation and dispersion of a monochromatic acoustic wave in a monodisperse spray with C = 0.106 -Line: analytic; ∆: CEDRE (BWC technique); 2: SAP1 (ACS splitting); •: CEDRE (ACS splitting).

  Fig. 16 Deformed-structured 27, 000 cell mesh of the LP10 geometry (Arrow: injection; Line: cut; +: pressure sensor).

1. 0 7 Fig. 17

 0717 Fig. 17 Instantaneous volume fraction per section in LP10 with CEDRE, MF-PEG method with ACS (top to bottom: S, M and L droplets).

Fig.

  Fig. Droplet mass concentrations along cuts (left: cut A-A, middle: cut B-B, right: cut C-C) -Solid: MF-PEG BWC (reference); •: MF-PEG ACS.

Fig.

  Fig. Velocities and temperatures along cuts (left: cut A-A, middle: cut B-B, right: cut C-C) -Solid: MF-PEG BWC (reference); •: MF-PEG ACS -Empty symbols: disperse phase (per section); Filled symbols: gas.

Fig.

  Fig. Droplet mass concentrations along cuts (left: cut A-A, middle: cut B-B, right: cut C-C) -Solid: MF-PEG reference (BWC, 1 curve for each of the 3 sections); ∆: MF BWC; 2: MF ACS. Compared to MF-PEG, the MF solutions have 1 more curve for the section of VS droplets.

  0

Fig.

  Fig. Velocities and temperatures along cuts (left: cut A-A, middle: cut B-B, right: cut C-C) -Solid: MF-PEG reference (BWC); ∆: MF BWC; 2: MF ACS -Empty symbols:disperse phase; Filled symbols: gas. Compared to MF-PEG, the MF solutions have 1 more curve for the section of VS droplets, which is almost superimposed to that of the gas, given the fast time scales of VS droplets.

Fig. 22

 22 Fig. 22 Spectra of pressure signal non dimensioned by the ACS highest peak -Left: front-end sensor (1); Middle: middle sensor (2); Right: rear-end sensor (3) -Solid: MF-PEG BWC (reference); •: MF-PEG ACS; ∆: MF BWC; 2: MF ACS.

Radius 10 MDF (m - 1 )Fig. 23 FourFig. 24 Fig. 25 Fig. 26 Fig. 27

 1012324252627 Fig.[START_REF] Doisneau | Two-way coupled simulation of acoustic waves in polydispersed coalescing two-phase flows : application to Solid Rocket Motor instabilities[END_REF] Four methods of discretization/reconstruction for the lognormal distribution using 4 sections; ---section bounds with area proportional to mass; ×: d31 computed with κ(r) = cst; 2: exact d31 -Top Left: Constant ∆S discretization method; Top Right: Constant ∆r discretization method; Bottom Left: Geometric discretization method; Bottom Right: Fixed mass ∆r discretization method.

Table 1

 1 Particle injection conditions Category Section bounds (µm) d 30 k diameter (µm) St ω,k St ,k Flow rate (kg/m 2 /s)

	N	[0, 3.33[	0.20	2.10 -4 0.3	0.15
	VS	[0, 3.33[	0.60	2.10 -3 2.8	0.15
	S	[3.33, 6.66[	5.02	0.51 8.10 2	0.17
	M	[6.66, 10[	8.14	1.34 2.10 3	0.30
	L	> 10	12.83	3.33 4.10 3	0.61
		Table 2 Timescale hierarchy (µs), exhibiting stiffness.
		τ u N Stokes time		0.1
		Time step ∆ta		0.2
		Time step ∆tc		0.2
		Nozzle CFL time (Full Mach)	0.6
		τ u V S Stokes time		1.7
		Nozzle strain-rate-based time	30
		τ u S Stokes time		120
		τ u M Stokes time		320
		Eddy revolution time		750
		τ u L Stokes time		780
		First acoustic mode (L/2)	1500
		Typical simulation time	500, 000

Acknowledgments

The present research was done thanks to a Ph. D. Grant from DGA, Ministry of Defence (M. S. Amiet, Technical Monitor), a collaborative project grant from the France-Stanford Center for