A Sibra

J Dupays

A Murrone

F Laurent

M Massot

Simulation of reactive polydisperse sprays strongly coupled to unsteady flows in solid rocket motors: Efficient strategy using Eulerian Multi-Fluid methods

Keywords: Reactive polydisperse sprays, Eulerian Two Size Moment Multi-Fluid methods, Time operator splitting, Dedicated scheme for complex evaporation laws, Solid propellant combustion, CEDRE code

In this paper, we tackle the issue of the accurate simulation of evaporating and reactive polydisperse sprays strongly coupled to unsteady gaseous flows. In solid propulsion, aluminum particles are included in the propellant to improve the global performances but the distributed combustion of these droplets in the chamber is suspected to be a driving mechanism of hydrodynamic and acoustic instabilities. The faithful prediction of two-phase interactions is a determining step for future solid rocket motor optimization. When looking at saving computational ressources as required for industrial applications, performing reliable simulations of two-phase flow instabilities appears as a challenge for both modeling and scientific computing. The size polydispersity, which conditions the droplet dynamics, is a key parameter that has to be accounted for. For moderately dense sprays, a kinetic approach based on a statistical point of view is particularly appropriate. The spray is described by a number density function and its evolution follows a Williams-Boltzmann transport equation. To solve it, we use Eulerian Multi-Fluid methods, based on a continuous discretization of the size phase space into sections, which offer an accurate treatment of the polydispersion. The objective of this paper is threefold: first to derive a new Two Size Moment Multi-Fluid model that is able to tackle evaporating polydisperse sprays at low cost while accurately describing the main driving mechanisms, second to develop a dedicated evaporation scheme to treat simultaneously mass, moment and energy exchanges with the gas and between the sections. Finally, to design a time splitting operator strategy respecting both reactive two-phase flow physics and cost/accuracy ratio required for industrial computations. Using a research code, we provide 0D validations of the new scheme before assessing the splitting technique's ability on a reference two-phase flow acoustic case. Implemented in the industrial-oriented CEDRE code, all developments allow to simulate realistic solid rocket motor configurations featuring the first polydisperse reactive computations with a fully Eulerian method.

Introduction

In solid rocket motors (SRM), aluminum particles are included in the propellant with a significant mass fraction to improve the global performances by increasing the chamber temperature. The distributed combustion of aluminum droplets in a portion of the chamber and the presence of polydisperse Al 2 O 3 residues in the entire volume that can be described as a spray can have a significant impact on the motor behavior. As one of the potential driving mechanisms of these instabilities, interactions between the spray combustion, hydrodynamic instabilities and the acoustics of the chamber could be a key issue. To improve solid rocket motor design, performing unsteady reactive two-phase flow computations with a reliable prediction of spray-acoustics interactions is a determining step to estimate flow instability levels [START_REF] Gallier | Aluminum combustion driven instabilities in solid rocket motors[END_REF][START_REF] Radenac | Fluctuating energy balance to post-process multiphase flow computations[END_REF][START_REF] Dupays | Two-phase unsteady flow in solid rocket motors[END_REF]. In the context of industrial computations, constraints such as an optimization of the cost/accuracy ratio have to be taken into account. An efficient numerical coupling strategy, which can be easily implemented in an industrial code is a challenge for scientific computing. In a solid rocket motor, droplet mass loadings are high enough so that the particles can appreciably modify the flow. This is called two-way coupling or strong coupling since the disperse phase and the gaseous carrier phase are strongly coupled through mass, momentum and energy exchanges. In this paper, we focus our study on the accurate description of the coupling dynamics between a moderately dense polydisperse spray and a carrier gaseous phase. We define a moderately dense spray as a dispersed phase constituted of spherical droplets of various size, the volume of which does not occupy more than one per cent of the total volume, thus leading to potential high mass loading without to take into account the liquid volume fraction in the gaseous phase description and for the derivation of the two phase exchanges. The following developments take into account a multi-species gaseous phase, representative evaporation and simplified droplet combustion models.

In order to treat polydisperse moderately dense sprays, we choose a kinetic approach in reference to the kinetic theory of gases. Using a statistical point of view, the spray is described as point particles for which global properties are established as the velocity of the center of mass, the temperature and as we focus on spherical particles, only one parameter for size. A number density function (NDF) describes the spray and its evolution, due to the drag force, the mass and heat transfers with the gaseous phase and the droplet-droplet interactions. The Williams-Boltzmann transport equation appears particularly relevant to follow the evolution of the NDF when dealing with polydisperse sprays [START_REF] Williams | Spray combustion and atomization[END_REF][START_REF] Williams | Combustion Theory (Combustion Science and Engineering Series)[END_REF]. To solve the liquid phase, two approaches are possible: Lagrangian methods solve accurately the NDF itself through a direct discretization of the transport equation. The repartition of the mass, momentum and heat source terms at the droplet locations onto the gas Eulerian grid was for a long time a challenging issue due to numerical diffusion but recent works [START_REF] Capecelatro | An Euler-Lagrange strategy for simulating particle-laden flows[END_REF][START_REF] Zamansky | Radiation induced turbulence in particle-laden fluids[END_REF] provide a new volume filtering approach that ensures interphase coupling convergence under mesh refinement. This technique is currently limited to structured meshes and the accuracy of the two-way coupling between both phases is still an issue especially when using high resolution methods for the gaseous phase. These recent Lagrangian approaches prove to be more competitive to scale up on parallel architectures. Nevertheless the ratio of accuracy versus computational expense is still to be investigated. Indeed, such a model can be costly when aiming at describing properly the interactions between a polydisperse phase coupled to unsteady gaseous flows. An alternative would be to use a Eulerian approach for non-inertial particles (smokes) following the work of Ferry and Balachandar and a Lagrangian description of more inertial particles [START_REF] Najjar | Simulations of solid-propellant rockets: Effects of aluminum droplet size distribution[END_REF][START_REF] Attili | Numerical simulation of multiphase flows in solid rocket motors[END_REF]. Moreover, it suffers from the same drawbacks as the usual Lagrangian approach. The fully Eulerian description of the spray provides an interesting alternative when dealing with parallel architectures and with two-way coupling issues.

The Eulerian approaches are based on moment methods of the NDF conditioned on droplet size leading to a "semi-kinetic" system for which size/velocity and size/temperature correlations are described. This enables to reduce the phase space dimension. Among Eulerian models, the Multi-Fluid methods (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] are derived from the "semikinetic" model by discretizing the size phase space using a finite volume approach and integrating continuously the size variable on fixed droplet size intervals. Systems of conservation equations are established for mass, momentum and energy on each size intervals called sections. These sectional systems can be seen as different "fluid" that are twoway coupled through source terms with a carrier fluid. At the Multi-Fluid level, size/velocity and size/temperature correlation resolutions are ensured considering velocity and temperature closures in each section. Other Eulerian methods, such as the so-called sampling model, adopt a discrete approach of the size variable using a sum of Dirac delta functions at given sizes. Systems of conservation equations are obtained for each droplet size and coupled to the carrier phase through source terms accounting for the gas-particle interactions. This sampling approach proves its efficiency in numerous reactive two-phase flow computations [START_REF] Lupoglazoff | Numerical simulations of the unsteady flow inside segmented solid-propellant motors with burning aluminum particles[END_REF][START_REF] Gallier | Aluminum combustion driven instabilities in solid rocket motors[END_REF][START_REF] Murrone | Numerical modeling of dispersed two-phase flows[END_REF] featuring complex evaporation and combustion models. We refer to [START_REF] Doisneau | Eulerian modeling and simulation of polydisperse moderately dense coalescing spray flows with nanometric-to-inertial droplets: application to Solid Rocket Motors[END_REF] for a detailed bibliography on the position of the sampling method among the other Eulerian models as well as for a precise connexion with what are sometimes called pivot and class methods. Nevertheless, this model, which does not rely on a continuous representation of the size variable, can hardly account for inter-particle interactions. A major advantage of the Multi-Fluid methods is to tackle particle-particle interactions such as secondary break-up or coalescence [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF]. Indeed, the MF formalism permits flux exchanges between sections so that new droplet sizes are easily handled. Moreover the MF methods appear as a flexible approach in order to describe accurately any size distribution of the NDF. The sampling method is appropriate to solve the evolution of size distributions for which the particles can be represented by discrete size values. Nevertheless, the method is less efficient to describe continuous size distributions. MF methods can follow complex evolutions of the size variable, which are non-local when dealing with coalescence or break-up and which are continuous for evaporating cases.

For Multi-Fluid methods, the droplet size-conditioned dynamics mainly results from the level of description of the size distribution function in the sections. A hierarchy of MF methods has been developed and it is based on the resolution of a certain number of size moments of the distribution: the One Size Moment method (OSM) was derived in [START_REF] Laurent | Numerical analysis of Eulerian Multi-Fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF], the Two Size Moment Multi-Fluid approach (TSM) has been introduced in [START_REF] Dufour | Modélisation Multi-Fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF] and the four size moment MF method (EMSM for Eulerian Multi Size Moment) has been proposed in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF][START_REF] Vié | Size-velocity correlations in high order moment methods for polydisperse evaporating sprays : modeling and numerical issues[END_REF]. They respectively yield to a first, second and fourth order convergence for the NDF description with the number of sections. Working with several size moments is not without constraints since realizability conditions have to be respected, which means that the moment vector should stay inside the moment space at anytime. This systematically leads to the development of realizable numerical methods and algorithms, which do not degrade the computational efficiency. When dealing with cost saving requirements, the TSM-MF methods offer an interesting compromise to use a limited number of sections while keeping size accuracy in the sections for complex dynamics. The TSM-MF approach and specific numerical methods for the coalescence source term computation were widely studied in [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF]. For evaporating cases, studies [START_REF] Dufour | Modélisation Multi-Fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF][START_REF] Laurent | Numerical analysis of Eulerian Multi-Fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF] were limited to simplified mass transfer models for which the evolution rate is the same for all the size range. For industrial simulations, advanced models featuring strong interactions between the gas environment and the size-conditioned droplet dynamics are required. So far piecewise exponential reconstructions are used as basis form functions in the sections; essential properties such as the positivity of the distribution are guaranteed. Nevertheless, the inversion step to determine the exponential function parameters from the size moments is costly and limits the competitiveness of the method. In this context, the TSM-MF method seems to be the best candidate for industrial simulations but it requires new develoments, extensive validations and well-suited numerical schemes need to be developed for accurately solving the NDF evolution in the size phase space. Once a method is identified, which is supposed to describe properly the physics for unsteady two-way coupled simulations, the choice of an efficient time integration strategy is crucial. Two-way interactions such as the coupling between the gaseous flow non-linearities or the acoustics of the chamber and the disperse phase dynamics underline the multi-time scale complexity of SRM computations. A first approach can be a phase-per-phase strategy, which integrates separately the gas and liquid systems; coupling source terms are modified and integrated at each phase integration step considering that the other phase is frozen. It has been shown that this time integration is not robust to solve stiff problems that are featured by a large range of numerical and physical time evolution rates and when having high mass loadings [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF]. Another numerical strategy based on IMEX methods [START_REF] Crouzeix | Une méthode multipas implicite-explicite pour l'approximation des équations d'évolution paraboliques[END_REF][START_REF] Shampine | IRKC: An IMEX solver for stiff diffusion-reaction PDEs[END_REF], which associate explicit/implicit schemes for the solution of stiff problems should be interesting. Nevertheless, the development of a such an approach depends on the mathematical bakcground of the governing equations. When modifying the physics of the problem, one has to rethink about the whole method and its implementation. In terms of code legacy, the flexibility of the implementation is a key condition, which appears as a constraint for an IMEX strategy. As an alternative and in order to optimize the cost/accuracy ratio of industrial-oriented codes, a time integration method based on an operator splitting seems particularly appropriate. It consists in decoupling physical phenomena and treating them in several operators which are integrated separately in time. This splitting technique has many advantages but the main one is the possibility to solve each problem with numerical dedicated methods for which the time step can be chosen in function of the physical time scales and the splitting time step is thus not constrained by numerical stability conditions [START_REF] Duarte | New resolution strategy for multi-scale reaction wave using time operator splitting, space adaptative multiresolution and dedicated high order implicit/explicit time integrators[END_REF]. By using the splitting strategy, we control two main sources of numerical errors: the first error is the so-called splitting error that is introduced because of the separate time solving of each time-dependent subproblem; such an error has been discussed in [START_REF] Hundsdorfer | Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations[END_REF][START_REF] Duarte | Adaptive numerical methods in time and space for the simulation of multi-scale reaction fronts[END_REF] for advection-reaction-diffusion systems. The second error results from the time and space discretization issues and associated numerical methods used for the time integration of the inner subproblems. In the context of two-phase flow computations, the splitting technique has been validated for inert particles in [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF] and studied for coalescing droplets in [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF][START_REF] Doisneau | Eulerian modeling and simulation of polydisperse moderately dense coalescing spray flows with nanometric-to-inertial droplets: application to Solid Rocket Motors[END_REF]. For reactive sprays, a new splitting strategy has to be designed looking for accuracy and robustness for two-phase interactions. The development of such a time integration is based on the identification and the analysis of the numerical peculiarities of capturing two-way coupling. Moreover, evolving time scales due the evolution of the polydispersity have to be taken into account. In order to perform competitive industrial computations, both flexibility and computational efficiency of the method are crucial so they have to be discussed and evaluated.

In this paper, we present a new TSM-MF model based on a relevant size reconstruction approach where the function is continuous and affine on a part of the section (zero otherwise). The accuracy of this new reconstruction proves to be the same as the exponential one but with numerous advantages especially in terms of CPU cost. In order to capture faithfully the evaporation dynamics, an accurate resolution of the mass, momentum and heat transfers with the gas is required. In addition, a well-suited numerical method is needed to account for complex mass transfer models in the context of MF methods. Indeed, the computation of flux exchanges between successive sections needs a specific numerical treatment, which preserves the realizability condition due to the use of two size moments. To ensure both robustness for the sectional exchanges and an efficient computation of the gas-particle interactions, we develop a dedicated scheme namely TSM-SRST scheme (for Simultaneous Reactive Source Term scheme). In a preliminary 0D verification, we do not consider any gas-particle couplings in order to focus on the evaluation of the exchanges between sections for different evaporation laws.

We then develop an efficient splitting strategy based on a time integration that separates gas and droplet convections from the coupling phenomena between reactive sprays and multi-species gaseous flows. The splitting technique is implemented in the research code SAP1 (Simulation d'Acoustique Polydisperse 1D) developed at EM2C. We propose a two-phase acoustic study that offers a relevant validation framework in order to evaluate two-phase flow solvers that aim at solving unsteady phenomena. The coupling between the gas and the spray expressed through the mass, momentum and energy exchanges is responsible for the dispersion and the dissipation of acoustic waves. Ensuring the physical propagation of acoustic disturbances appears as a demanding test case to demonstrate that a time splitting integration is relevant to account for unsteady two-phase interactions. First, the ability of MF methods to deal with evaporating two-phase acoustics is investigated, where some interesting new conclusions are drawn on the importance of polydispersity in comparison with related monodisperse studies [START_REF] Dupays | Propagation of acoustic waves in a two-phase vaporizing mixture[END_REF]. Then, simulations validate the time splitting strategy for different time step values. We also underline the accuracy of the integration. With this acoustic test case, we complete the validation of the new SRST scheme as regards the integration of the phase coupling source terms.

Finally, we implement and assess the affine TSM-MF method in the industrial-oriented code CEDRE (Calcul d'Ecoulements Diphasiques Réactifs pour l'Energétique). The CEDRE code, developed at ONERA, aims at doing research and industrial simulations of complex multi-physics problems. We propose reactive two-way coupled computations on two 2D solid rocket motor configurations. With the first configuration called TEU, we challenge the new affine reconstruction method to the exponential one and we highlight its advantage to save computational resources. Moreover, we implement and adapt the SRST scheme in order to match up to an industrial multi-solver architecture. With the CEDRE codre, we provide an additional validation of the SRST scheme comparing to a reference simulation using the sampling method. The second computation performed on the P230 configuration, which represents a full-scale solid rocket motor, features unsteady mechanisms such as acoustics or hydrodynamic structures coupled to droplet aluminum combustion. This feasibility test case illustrates the potential of both affine TSM-MF method and new SRST scheme to account for the impact of burning polydisperse spray on motor instability levels whereas it was limited to a monodisperse spray in previous studies [START_REF] Lupoglazoff | Numerical simulations of the unsteady flow inside segmented solid-propellant motors with burning aluminum particles[END_REF][START_REF] Gallier | Aluminum combustion driven instabilities in solid rocket motors[END_REF].

The paper is organized as follows: section 2 presents the derivation of the Eulerian Two Size Moment Multi-fluid model and its main hypotheses. The new affine basis form function is introduced and compared to the exponential one. We also present models at the level of a single particle, especially various heat/mass transfer models that we use in all the different numerical validations of this paper. In section 3, we design the time splitting integration strategy in the context of multiple time scale problems, which feature unsteady polydisperse two-phase flows. Moreover, it matches up to the numerical requirements of industrial two-phase flow solvers which are to have an acceptable cost/accuracy ratio. We also describe the numerical convection schemes used in the research code SAP1 and in the CEDRE code. In section 4, we derive the dedicated SRST scheme to integrate simultaneously the gas and the sectional source terms. Then, we propose a 0D validation of the TSM-SRST scheme using arbitrary evaporation laws in order to only validate the flux exchanges between sections. Section 5 provides numerical validations of our splitting strategy on a reference test case in the framework of linear evaporating polydisperse two-phase flow acoustics using the SAP1 code. This study also proves the SRST scheme's efficiency to tackle the gas-droplet interactions. Section 6 is composed of CEDRE computations on solid rocket motor configurations, which provide validations and illustrations of our developments in a complex physical background and using an industrial CFD platform.

Eulerian polydisperse two-phase flow modeling

Kinetic description of the disperse phase

The choice of a statistical description is particularly appropriate when studying a disperse phase that is considered as a set of point particles, which exchange mass, momentum and energy with a carrier continuous gaseous phase. Assuming that the spray is constituted of droplets of various sizes that are characterized by a finite set of global properties such as velocity and temperature, the polydisperse phase is described by a number density function (NDF) f where the quantity f (t, x, u, T, S )dxdudT dS denotes the average number of droplets, at a time t, in a volume of size dx around a space location x, in a dudT dS -neighborhood of velocity u, internal temperature T and droplet surface S . A Williams-Boltzmann-type equation, based on kinetic theory, is used to describe the evolution of the NDF since it has demonstrated to be useful for treating dilute and moderately dense polydisperse sprays [START_REF] Williams | Spray combustion and atomization[END_REF][START_REF] Williams | Combustion Theory (Combustion Science and Engineering Series)[END_REF][START_REF] Doisneau | Eulerian modeling and simulation of polydisperse moderately dense coalescing spray flows with nanometric-to-inertial droplets: application to Solid Rocket Motors[END_REF]. It is coupled through exchange terms to an Eulerian gas description such as Euler or Navier-Stokes approaches. Neglecting particles interactions such as collisions, coalescence and secondary break-up phenomena and focusing on external force, heat transfer and evaporation/condensation, the transport equation reads:

∂ t f + u.∂ x f + ∂ u . F p m p f + ∂ T H p m p c p,l f -∂ S K p f = 0 (1)
where F p and H p are the drag force and heat transfer and m p = (ρ l S 3/2 )/(6 √ π) is the mass of the particle of surface S ; the term ρ l is the liquid droplet material density, which is considered constant. Since we consider spherical droplets, the droplet size can be expressed indifferently through the surface variable or the diameter variable d or the radius variable r, that are linked by the relation S = πd 2 = 4πr 2 . We define c p,l as the pressure specific heat capacity of the liquid droplet. In this work, we neglect the dependance of c p,l on the droplet temperature. The term K p represents the size variation rate of the droplets. Without dilatation, only evaporation and condensation processes are modeled through the term K p . As a convention, we take K p > 0 when evaporation occurs. More complex models could be envisioned [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 can be handled within the proposed framework, however the chosen level of modeling will be adequate for the purpose of SRM reliable simulations. Models for the momentum, heat and mass transfers are presented in Appendix A. In the context of the development of reference numerical methods for Eulerian modeling of polydisperse sprays, studying complex mass transfer models aims at extending and at generalizing previous works [START_REF] Dufour | Modélisation Multi-Fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF][START_REF] Laurent | Numerical analysis of Eulerian Multi-Fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF] that focus on the d 2 -law for which the size variation rate K p is constant for all the droplet sizes. Detailed in Appendix A.2, Appendix A.3 and Appendix A.4, we describe several heat transfer/evaporation models, which are representative of physical models that are encountered in industrial two-phase simulations. In these models, the heat transfer is differently coupled to the mass exchange; so it offers a variety of interactions, which are relevant in order to evaluate numerical coupling strategies. Indeed, we offer:

• A two-phase acoustic validation framework for two-phase flow solvers in Section 5 using a thermal conductivitycontrolled evaporation model (c.f.Appendix A.2). Only the heat transfer drives the mass exchange between a particle and a single-species gaseous phase.

• A validation of the following new numerical methods and schemes in Section 6.1 using a representative heat/mass transfer model for the droplet exchanges mass with a multi-species gaseous phase. The evaporation dynamics is controlled by the diffusion of the species as described in Appendix A.3.

• An analysis of the potential of the presented developments in the context of solid propulsion simulations in Section 6.2. A simplified model for aluminum droplet combustion is given in Appendix A.4.

The semi-kinetic level

To follow the evolution of the NDF, we choose Eulerian moment methods that are based on moment methods of the NDF, which leads to the derivation of a set of conservation equations. The first step of Eulerian models relies on the reduction of the phase space dimension. Considering presumed NDF in velocity and temperature conditioned on size thus resolving size/velocity and size/temperature correlations, the phase space only depends on the droplet size variable. Moments in the velocity and enthalpy variable conditioned in size are then defined using the following formalism: ∀φ, φ = n -1 φ f dudT where n = f dudT is the droplet number density. Here the enthalpy is directly linked to the temperature such as h = h(T ). The closure of the system is obtained through the so-called monokinetic assumptions introduced in [START_REF] De Chaisemartin | Eulerian models and numerical simulation of turbulent dispersion for polydisperse 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]:

[HV1] At (t, x) and for a given droplet size S , there is only one velocity u(t, x, S ) and the velocity dispersion is zero in each direction.

[HT1] At (t, x) and for a given droplet size S , there is only one temperature T (t, x, S ) and the temperature dispersion is zero.

The following NDF is now only conditionned by droplet size:

f (t, x, u, T, S ) = n(t, x, S )δ(u -u(t, x, S ))δ(T -T (t, x, S )) (2)
This reduces the support of the NDF to a one dimensional sub-manifold parametrized by droplet size. Finally, this step leads to a system of conservation laws called the semi-kinetic model:

           ∂ t n + ∂ x (nu) = ∂ S (nK) ∂ t (nu) + ∂ x (nu ⊗ u) = ∂ S (nKu) + nF ∂ t (nh) + ∂ x (nhu) = ∂ S (nKh) + nH (3)

Multi-Fluid model hierarchy and related assumptions

The second step to derive the Eulerian Multi-Fluid model (MF) is based on a finite volume discretization of the size phase space such as 0 = S 0 < ... < S N sec = ∞ for the droplet surface variable. A system of conservation equations is obtained for each fixed size interval [S k-1 , S k [ called section. 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. The final system is written on a set of moments in size of n,nu and nh; so as to close the system, three other assumptions are necessary. The first one is:

[HS1] In each section, the form of n can be presumed and expressed as a function of S . This fundamental assumption leads to the choice of a size basis form function κ (k) in each section at (t, x):

n(t, x, S ) ≈ κ (k) (t, x, S )1 [S k-1 ,S k [ (S ) (4)
At space-time location (t, x), the size distribution profile in a section is reduced to a set of moments of S. In the classical approach developed 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], referred to as the one size moment Multi-Fluid model (OSM-MF), a constant form for the κ (k) function is used such that κ (k) = m (k) (t, x)κ (k) (S ) where m (k) is the bulk mass density of droplets. This one parameter function yields a first order size convergence with the number of sections [START_REF] Laurent | Numerical analysis of Eulerian Multi-Fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF]. When studying evaporation or coalescence, a high number of sections may be needed to get an accurate size evolution of the distribution function and a better description of the polydispersity of the spray. So higher order size moment methods have been developed in order to improve the quality of the size distribution profile in each section by using multi-parameter basis functions.

The number of coefficients determines the number of size moments to be solved and we look for a family of profiles for which there is a one to one correspondance between the chosen moments and the kinetic distribution. A Two Size Moment model (TSM-MF) was developed in [START_REF] Dufour | Modélisation Multi-Fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF] and studied in [START_REF] Doisneau | Eulerian modeling and simulation of polydisperse moderately dense coalescing spray flows with nanometric-to-inertial droplets: application to Solid Rocket Motors[END_REF][START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF] for coalescing droplets. A four size moment MF method (EMSM for Eulerian Multi Size Moment) has been proposed in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF][START_REF] Vié | Size-velocity correlations in high order moment methods for polydisperse evaporating sprays : modeling and numerical issues[END_REF] to reach a high level of precision in the size phase space treatment. Two assumptions are required on the velocity and the temperature to close the MF model:

[HV2] In each section, the velocity u(t, x, S ) = u (k) (t, x) does not depend on the droplet size.

[HT2] In each section, the temperature T (t, x, S ) = T (k) (t, x) does not depend on the droplet size. Assumption [HT2] is equivalent to consider a constant enthalpy distribution in section k, namely h(t, x, S ) = h (k) (t, x) = h(T (k) (t, x)). The validity of hypotheses [HV2] and [HT2] relies on the range of droplet characteristic times in a given section. In other words, when comparing the Stokes numbers of the smallest droplets to the biggest ones encountered in a section, the assumptions of constant velocity and temperature distributions may be no longer valid. One of the objectives of high order moment methods such as the EMSM model is to use a very low number of sections or a unique section especially when dealing with industrial applications. In solid propulsion applications, the size polydispersity is quite important, the spectra of Stokes numbers in the sections are wide so that too few sections yields poor representativity of the droplet dynamics, heating and combustion. Using a high order size moment model with a relative high number of sections is not relevant in terms of computational cost. Nevertheless, another family of high order methods such as Coupled-Size-Velocity Moment method [START_REF] Vié | Size-velocity correlations in high order moment methods for polydisperse evaporating sprays : modeling and numerical issues[END_REF] are well-suited to tackle non constant velocity profiles in the sections. When dealing with high-inertia particles compared to the time scales of the flow, particle trajectory crossing (PTC) can occur. At a given location and given range of particle size, the velocity distribution is multi-valued. Because of the monokinetic hypothesis [HV1], the Eulerian MF method fail to describe the physics of PTC. A new method of moments for particle-laden flows based on a Gaussian velocity closure and a second order MUSCL/HLL transport scheme are proposed in [START_REF] Vié | On the anisotropic gaussian velocity closure for inertial-particle laden flows[END_REF]. In specific SRM applications where coalescence occurs, inertial droplets with large Stokes numbers are created and lead to PTC in a specific portion of the chamber, which is the nozzle. In the present contribution, we do not study coalescing droplets, so that the occurrence of droplet trajectory crossings is very limited and the assumptions [HV1] and [HT1] are reasonable. Moreover, accounting for turbulent two-phase flows needs a peculiar modeling effort. In [START_REF] Vié | Accounting for Polydispersion in the Eulerian Large Eddy Simulation of the Two-Phase Flow in an Aeronautical-type Burner[END_REF], a new approach based on a Multi-Fluid Mesoscopic Eulerian Formalism has been proposed. It consists in capturing polydispersity with size-conditioned turbulent dynamics and it demonstrates its efficiency for Large Eddy simulation of two-phase flows. LES computations and associated developments are not treated in this work since turbulence is barely developed in the presented SRM configurations. In this contribution, we choose Two Size Moment MF models. They offer an interesting compromise between a precise resolution of the size evolution of the spray with a second order size convergence method and the industrial computational requirements by using a moderate number of sections.

Two Size Moment Multi-Fluid models

The chosen size moments of TSM-MF methods are the moment of order 0, which is the number concentration of droplets and the moment of order 3/2, which matches to the bulk mass density of droplets:

n (k) (t, x) = S k S k-1 κ (k) (t, x, S )dS m (k) (t, x) = ρ l 6 √ π S k S k-1 S 3/2 κ (k) (t, x, S )dS (5)
Those models yield two main issues: the first problem is to choose a family of size basis functions κ (k) such as to ensure that the density function remains positive. The other constraint when using TSM-MF methods consists in preserving the moment space which means that the moments n (k) and m (k) in each section should respect the realizability condition given by the relation:

S 3/2 k-1 n (k) < 6 √ π ρ l m (k) < S 3/2 k n (k) (6)
To ensure the positivity of the density function, a relevant size basis distribution function was proposed in [START_REF] Dufour | Modélisation Multi-Fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF][START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF] by using a piecewise exponential form:

κ (k) (t, x, S ) = a k (t, x) exp(-b k (t, x)S ), S ∈ [S k-1 , S k [, a k = Ψ a (n (k) , m (k) ), b k = Ψ b (n (k) , m (k) ) (7)
The numerical procedure, which consists in determining (a k , b k ) k 's from the size moments is called the inversion step. We use an algorithm based on a Ridder's method that is proven to be reliable but its computational cost is significant. Moreover a numerical limitation is needed for the exponential parts b k S 's of the NDF for certain steep size distributions. In these cases, the accuracy of the inversion method is reduced and computational limits can be reached for calculators. We propose a new reconstruction where the function inside the section is continuous and affine on a part of the section (zero otherwise) as represented in Fig. 2. This reconstruction uses four parameters that guarantee the positivity condition:

κ (k) (t, x, S ) = 0, S ∈ [S k-1 , S a ] or S ∈ [S b , S k [ κ (k) (t, x, S ) = α k + (β k -α k ) (S -S a ) (S b -S a ) , S ∈ [S a , S b ] (8) 7 f a (S) 6 Z Z Z Z Z Z a k S k 1 =S (k) a b k = 0 S (k) b S k f a (S) 6 a k b k S k 1 =S (k) a S k =S (k) b f a (S) 6 ⇢ ⇢ ⇢ ⇢ ⇢ ⇢ b k S k 1 a k = 0 S (k) a S k =S (k) b Figure 1:
The three types of a ne reconstructions.

The reconstructed function f a is then taken continuous inside each section ]S k 1 ,S k [ and affine on the whole section or on a part of it, the function being zero elsewhere. This function can then be written:

f a (S) = 8 > > < > > : 0 if S k 1  S  S (k) a a k +(b k a k ) S S (k) a S (k) b S (k) a if S (k) a  S  S (k) b 0 if S (k) b  S  S k i 2 {1,...,N}, (3.4)
with only three considered cases illustrated in Figure 3:

[C.1] S (k) b 2]S k 1 ,S k [, S (k)
a = S k 1 and b k = 0. [C.2] S (k) a = S k 1 and S (k) b = S k . [C.3] S (k) a 2]S k 1 ,S k [, S (k)
b = S k and a k = 0. Let us denote R k 1 = p S k 1 , R b = q S (k) b and µ (k) in f = 2 DS 2 k Z S k S k 1 (S k S)S 3/2 dS, µ (k) sup = 2 DS 2 k Z S k S k 1 (S S k 1 )S 3/2 dS.
These two last scalars belong to ]S 3/2 k 1 ,S

Z S k S k 1 f a (S)dS = n k , Z S k S k 1 S 3/2 f a (S)dS = m k . (3.5)
The function's parameters are given by: Among the four parameters α k , β k , S a and S b used for the reconstruction, two are always fixed, depending on the three possible cases which are illustrated in Fig. 2. The choice of the case depends on the value of the ratio m (k) /n (k) . More details on the reconstruction and on its mathematical properties are given in [START_REF] Laurent | Two-size moment multi-fluid model: a robust and high-fidelity description of polydisperse moderately dense evaporating sprays[END_REF].

The affine reconstruction is an interesting alternative to the exponential technique as it offers a significant computational gain when processing the inversion step. When looking at a relevant accuracy/cost ratio, an improvement on such a numerical procedure is a key development. Moreover, this affine formalism avoids computation issues as those encountered by the exponential form. Multiple advantages of the affine size basis function are illustrated in the following Sections such as the source term computation. CPU cost comparisons between the exponential and affine approaches are provided in section 6.1.

Finally, the two-phase coupled system of conservation equations based on a Two Size Moment Multi-Fluid model reads:

                                                                 ∂ t ρ g + ∂ x . ρ g u g = N sec k=1 M (k) ∂ t ρ g Y i + ∂ x . ρ g Y i u g = 0, i = 1, N spec -1 Y N spec = Y vap ∂ t ρ g u g + ∂ x . ρ g u g ⊗ u g + p = - N sec k=1 F (k) -M (k) u (k) ∂ t ρ g E g + ∂ x . (ρ g E g + p) u g = - N sec k=1 H (k) + F (k) .u (k) -M (k) h (k) + L (k) v + 1 2 u (k) 2 ∂ t n (k) + ∂ x .(n (k) u (k) ) = N (k+1) -N (k) ∂ t m (k) + ∂ x .(m (k) u (k) ) = E (k+1) -E (k) -M (k) ∂ t (m (k) u (k) ) + ∂ x .(m (k) u (k) ⊗ u (k) ) = E (k+1) u (k+1) -E k u (k) -M (k) u (k) + F (k) ∂ t (m (k) h (k) ) + ∂ x .(m (k) h (k) u (k) ) = E (k+1) h (k+1) -E (k) h (k) -M (k) h (k) + H (k)                  k = 1, N sec (9)
The N sec fluid systems are coupled to N spec + 2 Euler gas equations through their source terms. Here N spec is the number of species in the gaseous phase and in the rest of the work the term Y vap represents the vapor coming from the droplet evaporation. The gas system is the classical of transport equations on the local density ρ g , the velocity u g and the total energy

E g = E re f + T T re f c v,g dT + 1 2 u g 2 with c v,g
is the specific volume heat capacity and E re f is a reference energy calculated at a reference temperature T re f . For a sake of legibility, we choose the classical Euler system but the Navier-Stokes approach is also possible and, besides, it is used in Sections 6.1 and 6.2 with the CEDRE code. Note that we have restricted the study to moderate dense sprays so the volume occupied by the disperse phase remain a very small fraction of the total volume and can be neglected in the carrier phase modeling. The gas source terms that are the contribution of all the sections satisfy the conservation of mass, momentum and energy of the two-phase mixture.

For N sec fluid systems written in the Two Size Moment formalism and considering n d dimensions for the physical space, we have N sec (n d + 3) conservation equations. For each liquid system, mass, momentum and energy source terms are obtained by averaging K, F and H on the sections. We finally obtain the MF source terms:

N (k) = K(S k-1 ) κ (k) (S k-1 ), E (k) = ρ l S 3/2 k-1 6 √ π N (k) M (k) = ρ l 6 √ π S k S k-1 d S (S 3/2 ) K(S ) κ (k) (S )dS , κ (k) (S ) = κ (k) (S ) n (k) , m (k) F (k) = ρ l 6 √ π S k S k-1 S 3/2 F(S ) κ (k) (S )dS , H (k) = ρ l 6 √ π S k S k-1 S 3/2 H(S ) κ (k) (S )dS (10)
MF systems are featured by two types of source terms. The first ones represent the exchange terms between the gaseous phase and a given section through average exchange terms given by M (k) , F (k) and H (k) ; this is the so-called average sectional exchange terms. The other source terms are flux exchange terms between successive sections. The number of droplets transfered from the k th section to the (k -1) th section is expressed through the N (k) term. Similarly, the mass exchange term between the sections are given by E (k) .

Another Eulerian polydisperse model: the sampling model based on a discrete size approach

Another method to get a polydisperse size description of sprays consists in realizing a sampling of the distribution but keeping an Eulerian description. The NDF is approximated by a sum of mass Dirac delta functions as shown in Fig. [START_REF] Attili | Numerical simulation of multiphase flows in solid rocket motors[END_REF]. This discretization can be seen as N sample samples at given size locations {S 1 , ..., S N sample } in the size phase space for which mean velocities u (k) 's and mean temperatures T (k) 's are considered:

f (t, x, u, T, S ) = N sample k=1 n (k) (t, x) δ(S -S k (t, x)) δ(u -u (k) (t, x)) δ(T -T (k) (t, x)) (11)
For each sample, we define the number concentration of droplets n (k) and the bulk mass density m (k) = ρ l n (k) S 3/2 k / (6 √ π).

The liquid system of conservation equations for the k th sample reads:

               ∂ t n (k) + ∂ x .(n (k) u (k) ) = 0 ∂ t m (k) + ∂ x .(m (k) u (k) ) = -n (k) ṁ(k) p ∂ t (m (k) u (k) ) + ∂ x .(m (k) u (k) ⊗ u (k) ) = -n (k) ṁ(k) p u (k) + n (k) F (k) p ∂ t (m (k) h (k) ) + ∂ x .(m (k) h (k) u (k) ) = -n (k) ṁ(k) p h (k) + n (k) H (k) p ( 12
)
where ṁ(k) p , F (k) p and H (k) p are respectively the mass transfer term, the drag force and the heat exchange term for a single particle of a given droplet size S k . As presented in Eq.( 9), the N sample liquid systems are coupled to the gas system through the contributions of all the source terms.

In solid and liquid propulsion, this model has been used for reactive two-phase flow computations in the context of motor instability studies [START_REF] Lupoglazoff | Numerical simulations of the unsteady flow inside segmented solid-propellant motors with burning aluminum particles[END_REF][START_REF] Gallier | Aluminum combustion driven instabilities in solid rocket motors[END_REF] and more recently in the CEDRE code [START_REF] Murrone | Numerical modeling of dispersed two-phase flows[END_REF][START_REF] Le Touze | Numerical coupling strategies for the separated-to-dispersed transition within the liquid phase of cryogenic rocket engines[END_REF]. As a reference approach for solid rocket motor simulations featuring complex mass transfers and combustion, the sampling model is presented here in order to validate the TSM-MF method for which it is the first evaluation. As an additional remark, the sampling method is accurate for solving the evolution of discrete mass distributions but it proves to be less efficient when continuous distributions are considered [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]. In solid propulsion, lognormal size distributions and smooth bimodal size distributions are commonly encountered.

Unsteady two-way coupling strategy and dedicated numerical methods

In this section, we first present a relevant time operator strategy designed to take into account the main constraints of unsteady two-phase flow simulations. The proposed time splitting technique consists in successively solving gas and liquid convection and phase coupling subproblems with associated time independent operators. Then, we describe the time integrations and the numerical methods used in the transport operators and their implementation in the research code SAP1 and in the CEDRE code. After that, we focus on the source term operator by exposing the numerical peculiarities of TSM-MF models to account for the computation of the average sectional terms. We also highlight the difficulty to correctly deal with the flux exchanges between sections, which require a specific treatment in order to respect MF realizability conditions. As a main development of this paper, we finally provide a dedicated scheme that aims at satisfying those conditions and at simultaneously integrating mass, momentun and energy exchanges with the gas, especially in the context of complex evaporation models.

Choice of a time splitting integration

One major issue when dealing with unsteady two-way coupling problems featuring complex phenomenon interactions is the time integration strategy. A time coupling technique based on an operator splitting seems appropriate to treat unsteady two-phase flow problems. The use of specific numerical methods for each operator and the possibility to have internal time steps respecting the stability conditions lead to the accurate resolution of each problem. As explained before, we highlight that the two-phase coupling has to be treated in the same subsytem and has to be solved in a dedicated solver, which is necessary in order to obtain the asymptotic equilibrium for this subproblem. Moreover, two-phase exchanges feature commonly large ranges of time evolution rates, this implies to be stiffly robust and to use accurate numerical methods [START_REF] Hairer | Solving ordinary differential equations II, Stiff and differential-algebraic problems[END_REF]. A preliminary version of the scheme has been introduced in [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF] for inert two-phase coupling; in Subsection 4.3 of this paper we propose a dedicated scheme that extends two-phase exchanges to mass transfers.

Splitting techniques are also interesting to reach an acceptable cost/accuracy ratio. As long as the numerical stability is guaranteed in each solver of the split subproblems, the global stability of the splitting method is ensured. As a consequence, large time steps can be chosen for certain operators. The advantage of splitting methods to use large time step can be also a main drawback. Indeed a significant splitting error appears when the overall time step is greater than the physical coupling time scales. In any case, splitting strategies introduce a decoupling error but it is acceptable if the splitting time step is below or close to the smallest physical coupling time scale. Moreover the numerical errors introduced by the schemes to solve each subproblem are far from being negligible, we do not want them to be at the same order as the splitting error. It is crucial to have a limited amount of numerical error when integrating in time the subsystems. In this case, the global error on the numerical solutions resulting from the splitting technique mainly results from the splitting error, which can be controlled by choosing the suitable splitting approach as explained in Section 3.4.

SRM time scale hierarchy and industrial requirements

In order to get an acceptable CPU time cost, we are interested in using large splitting time steps, which means that the splitting time step shoud not be constrained by the fastest time scales of the problem. Four time scales characterize typically reactive two-phase flow systems: the gaseous and liquid convection scales τ c g and τ c liq , the droplet relaxation scales τ min = min(τ u , τ T ) and the evaporation one τ evap . SRM instabilities result from a coupling between hydrodynamic structures and acoustics; it is crucial to solve accurately these time scales. For both accuracy and CPU cost gain requirements, we choose to set the splitting time step to solve the convection and acoustic time scales. Classical SRM meshes are designed to avoid harsh CFL values especially in the nozzle where the flow is supersonic, which yields classical values (τ c g , τ c liq ) ≈ 10 -6 s. The droplet evaporation time scales are not troublesome since they are usually around 10 -4 s. Most of the droplet relaxation time scales are above the convection ones but stiffness emerges from micronic and nanometric particles with time scales that fall far below 10 -6 s. Taking into account the whole size range of time scales cannot be envisioned for an industrial approach. When using ∆t ≈ τ c g , the coupling between the gas and the smallest droplets is not correctly solved. As explained in Subsection 3.1, it introduces a splitting error that may impact the accuracy level expected for industrial simulations. Nevertheless for SRM computations, the bulk of the smallest droplets feature St << 1 and behave like tracers. They do not significantly modify the gas flowfield and the local two-phase equilibrium as explained in [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF]. Despite a splitting time step limited by the convection and acoustic time scales, most of the gas/particle interactions are captured. In Subsection 5.5 and 6.1.2, we illustrate the impact of several splitting time steps on the numerical solution.

As an additional remark, it may be crucial to choose specific methods that can handle discontinuities such as a shock. Splitting methods are still possible to treat discontinuities and singularities. Nevertheless, as explained in [START_REF] Billet | Impact of volume viscosity on a shock/hydrogen bubble interactions[END_REF] these techniques lead to the use of time steps smaller than all the time scales of the problem. It is important to specify that our study only concerns choked nozzles presenting smooth flows.

Splitting decomposition in subsystems

We choose to decouple the convective liquid and gas phenomena from the phase coupling. We consider system (9) as composed of transport operators T g and T k 's and a source term operator S. These operators apply on Q(t, x) ∈ R N sec (n d +3)+N spec +n d +1 where we recall that n d is the number of dimensions of the physical space. One can write:

Q(t, .) = [ρ g , ρ g Y 1 , ..., ρ g Y N spec -1 , ρ g u g , ρ g E g , (n (k) , m (k) , m (k) u (k) , m (k) h (k) ) k=1,N sec ] t (t, .)
The gas transport semigroup T g accounts for full-Mach gaseous dynamics i.e. convection and acoustics:

T g                  ∂ t ρ g Y i + ∂ x .(ρ g Y i u g ) = 0 i = 1, N spec -1 ∂ t ρ g + ∂ x .(ρ g u g ) = 0 ∂ t ρ g u g + ∂ x . ρ g u g ⊗ u g + p = 0 ∂ t ρ g E g + ∂ x . (ρ g E g + p) u g = 0 (13)
Note that T ∆t g Q 0 is the solution after a time step ∆t from the initial condition Q(0, .) = Q 0 ; the liquid variables are unchanged. As presented in Eq.( 9), the Euler system is chosen for the gas description. Nevertheless, the generalization to Navier-Stokes equations taking into account viscous effects is possible. The use of splitting techniques for reaction/convection/diffusion problems has already been discussed and validated in [START_REF] Duarte | New resolution strategy for multi-scale reaction wave using time operator splitting, space adaptative multiresolution and dedicated high order implicit/explicit time integrators[END_REF][START_REF] Duarte | Adaptive numerical methods in time and space for the simulation of multi-scale reaction fronts[END_REF]. In the scope of our solid propulsion simulations, Reynold numbers are particularly high and the meshes are not refined enough at the boundary layers. Since the diffusion time scale is slower than the convection one, we choose to integrate the viscous term simultaneously in the gas convection operator. The liquid convection occurs independently for all sections so we have N sec transport operators:

T k                    ∂ t n (k) + ∂ x . n (k) u (k) = 0 ∂ 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 (14)
Note that T ∆t k Q 0 is the solution after a time step ∆t from the initial condition Q 0 ; the gas variables are unchanged. All the "fluid" transport operators for both gas and droplets evolve at their own advection velocity without any coupling. As a consequence, different CFL values determine their stability conditions. When choosing a splitting time step that is superior to the stability time step of a subproblem, it could be interesting to optimise the internal time steps independently in each operator in order to introduce a limited amount of numerical diffusion and respect their stability domains. Moreover, it offers a gain of computational time since the stability condition of certain operator do not appear as a constraint for the others that use larger time steps. Such a decoupling method finally contributes to get an interesting cost/accuracy ratio. The source term subsystem S ensures all the coupling interactions between the gas and the polydisperse phase: [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF] Note that S ∆t Q 0 is the solution after a time step ∆t from the initial condition Q 0 . This operator represents a collection of ODE systems for all the gas and liquid variables parametrized by the space location x. For each system, the gas and liquid variables relax to an equilibrium state at a rate that is a function of the different droplet characteristic times. Note that the influence of the polydisperse phase on the gas is represented through the sum of the average exchange terms of all sections, which ensures the global conservation of mass, momentum and energy of the two-phase system. In the context of evaporating droplets, sections are also strongly coupled since the mass transfers written in the context of Multi-Fluid methods introduce exchange source terms between successive sections. As a consequence, the choice of dedicated time integration methods is essential to ensure robustness and an accurate coupling between all "fluids".

S                                                                  ∂ t ρ g = N sec k=1 M (k) ∂ t ρ g Y i = 0, i = 1, N spec -1 Y N spec = Y vap ∂ t ρ g u g = - N sec k=1 F (k) -M (k) u (k) ∂ t ρ g E g + ∂ x . (ρ g E g + p) u g = - N sec k=1 H (k) + F (k) .u (k) -M (k) h (k) + L (k) v + 1 2 u (k) 2 ∂ t n (k) = N (k+1) -N (k) ∂ t m (k) = E (k+1) -E (k) -M (k) ∂ t (m (k) u (k) ) = E (k+1) u (k+1) -E k u (k) -M (k) u (k) + F (k) ∂ t (m (k) h (k) ) = E (k+1) h (k+1) -E (k) h (k) -M (k) h (k) + H (k)                  k = 1, N sec

Splitting structure

Various splitting methods are possible such as the first order Lie splitting scheme [START_REF] Trotter | On the product of semi-groups of operators[END_REF] and the second order Strang one [START_REF] Strang | On the construction and comparison of difference schemes[END_REF]. The first scheme can be seen as a one-step technique for which the subproblems are successively solved with its associated operator and the same splitting time step. The Strang approach is a symmetric composition of one-step methods using half-sized time steps to solve each subsystem. High order methods have been developed [START_REF] Schatzman | Towards non commutative numerical analysis: High order integration in time[END_REF][START_REF] Thalhammer | High-order exponential operator splitting methods for time-dependent schrodinger equations[END_REF][START_REF] Castella | Splitting methods with complex times for parabolic equations[END_REF][START_REF] Hansen | High order splitting methods for analytic semi-groups exist[END_REF] but they imply significant stability restrictions and they are difficult to implement. We insist on the fact that it is judicious to promote, for each operator, numerical schemes featuring a higher time accuracy than the splitting method. It would be useless and costly to choose a high order splitting integration technique with inferior order methods. Indeed the computational time accuracy is constrained by the lowest order scheme.

Convection/Source splitting

The key point is to choose appropriately how to organize the source term operator and the transport ones. Several splitting compositions offering different time orders are possible: the Strang splitting scheme ensures a second order convergence in time as long as there is no shock in the solution; the Lie splitting is first order in time. For both schemes, the way the operators are sequenced has an impact on the level of the splitting error. It has been demonstrated in [START_REF] Descombes | Operator splitting for nonlinear reaction-diffusion systems with an entropic structure: singular perturbation and order reduction[END_REF] that operators that have the fastest time scales have to be resolved at the end of the splitting for a maximal order and minimal impact of fast scales on splitting errors. As seen in Appendix A and Section 3.2, the fastest time scales are mostly due to small droplets in SRM applications. So we suggest two order schemes that both end up with the source term subsystem:

Q Strang = S ∆t/2 [ N sec k=1 T k + T g ] ∆t S ∆t/2 Q 0 ( 16
) Q Lie = S ∆t [ N sec k=1 T k + T g ] ∆t Q 0 ( 17
)
We define Q Strang and Q Lie as the solutions after a splitting time step ∆t, which result from the exact resolution of each subsystem that are successively integrated from the initial condition.

Two-stage splitting: acoustic operator

For specific problems dealing with large physical time scale spectrum, one or several operators such as the transport ones do not need to be called as often as the others. Larger time steps can be used in order to have a low level of numerical dissipation and to spare CPU resources. To do so, multi-stage splitting technique are appropriate. As a relevant example, the acoustic study presented in Section 5 features slow droplet convection compared to the acoustics and two-phase coupling time scales. So we use a two-stage splitting approach to call the liquid transport operators less often than the gas convection and the source term ones. In order to gain computational ressources, we choose the Lie splitting formalism, it is written:

Q acou = A ∆t c [ N sec k=1 T k ] ∆t c Q 0 ( 18
)
At this step, the operators are exactly solved. Nevertheless when applying the space discretization and when choosing the numerical schemes, the splitting time step ∆t c will be chosen to be the most restrictive droplet convective time step that guarantees the liquid CFL conditions. The so-called acoustic operator A performs gas convection and two-phase interactions. Note that the acoustic operator ends up the upper stage of the splitting because we consider it as the stiffer subproblem compared to the liquid convection operator. In order to be stiffly accurate on the two-phase coupling, we use the Strang splitting scheme to solve the acoustic operator:

A ∆t a = S ∆t a /2 T ∆t a g S ∆t a /2 , A ∆t c = [A ∆t a ] ∆t c /∆t a , ∆t c /∆t a ∈ N (19)
As before, the splitting time ∆t a will be chosen when dealing with the time integration of the subproblems with respect to two-phase acoustic characteristic time scales.

3.5. Numerical methods for the transport operators and implementation in the SAP1 and CEDRE codes Previously, we described relevant time splitting strategies to deal with two-phase problems by decomposing it into subsystems. We are now interested in the time and space discretization issues and the choice of the time integration methods for each operator. We aim at presenting the space discretizations and the numerical methods that are implemented in the two-phase CFD codes, which we use to perform the following numerical validation studies. We are particularly interested in the schemes of the T g and the T k 's operators, which account for the gas and liquid convection.

The SAP1 code, which is a research code, has been designed to evaluate Eulerian Multi-Fluid models describing polydisperse sprays using new numerical schemes and different time splitting integration strategies for a high-fidelity two-way coupling computations and acoustic studies. Note that we deal with a single-species carrier phase. For the space discretization, we choose structured 1D meshes that are particularly relevant to capture two-phase flow acoustics. To do so, the gas convection is solved with a scheme based on a fifth order WENO space fluxes [START_REF] Shu | Essentially-non-oscillatory and weighted essentially-non-oscillatory schemes for hyperbolic conservation laws, advanced numerical approximation of nonlinear hyperbolic equations[END_REF] that is time integrated with an explicit RK3 scheme, which is necessary to ensure stability as specified in [START_REF] Wang | Linear instability of the fifth-order weno method[END_REF]. For the liquid transport, a second order space and time explicit kinetic scheme [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF][START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF] has been implemented.

The CEDRE code is an industrial-oriented multi-physics platform that is mainly used in the fields of energetics and propulsion [START_REF] Refloch | CEDRE software[END_REF][START_REF] Scherrer | Recent CEDRE applications[END_REF]. The software architecture is based on a multi-solver strategy for which several solvers are considered to tackle specific physical problems. Spatial parallelized computations can be performed using a multidomain approach and general unstructured meshes. The gas solver accounts for multi-species gas in the context of reactive compressible flowfields. The gas transport scheme that is used in this paper is based on second order MUSCL space fluxes; both implicit and explicit time integrations are possible and reach a third order in time. The Eulerian disperse phase solver tackles evaporating polydisperse sprays with the Eulerian sampling model or with exponential and affine MF models. The liquid transport is a second order numerical method dedicated to pressureless gas dynamics but it is here extended to general unstructured meshes [START_REF] Le Touze | Multislope MUSCL method for general unstructured meshes[END_REF].

Computation of the gas/particle coupling source terms in the context of TSM-MF models

Now that we have detailed the numerical methods for the transport in the physical space, we are concerned with the treatment of the two-phase interactions especially with the mass exchanges in the context of TSM-MF methods. So we focus on the intra sectional reconstructions and on the time integration of the source term operator S.

We divide this section into four parts: the first one concerns the computation strategy to compute the average exchange terms in the sections, then we highlight the numerical peculiarities to deal with the mass flux exchanges between sections. In the third Subsection, we present a new scheme in order to solve those issues and to integrate simultaneously all the two-phase exchange terms. Finally, we propose 0D validations and a numerical analysis of the new scheme.

Intra sectional reconstruction and related source term evaluation

The specificity of MF methods is to compute averaged source terms using a basis function κ (k) in each section. For computational efficiency, a crucial point concerns the numerical methods that are chosen to operate the integrations over the sections. For the TSM-MF models, we point out the advantages of the affine reconstruction compared to the exponential one.

At a given time and location, sectional hypotheses [HV2], [HT2] and [HS1] guarantee that the integrands of source terms such as F (k) and H (k) are only functions of the size variable S . As given in Eq.( 10), the computation of the drag force and heat transfer source terms can be expressed as fractional moments of κ (k) . The same conclusion on the mass exchange terms can be applied as expressed in Subsection 2.4. When choosing the exponential reconstruction, it appears that the integrals are difficult to compute because the integrand has the following form:

I p/q = α k S k S k-1 S p/q e -β k S dS (20)
When q = 1 or 2, the integral can be integrated analytically but it can be particularly costly; for higher values of q numerical strategies have to be used. Quadrature methods can be chosen such as Gauss-Legendre based on fixed N-node distribution as illustrated by the relation:

I p/q ≈ N i=0 W GL i f (S k -S k-1 ) 2 Y GL i + 1 + S k-1 (21)
where Y GL i and W GL i are respectively the abscissas and the weights of the N quadrature nodes. The main problem of fixed node approaches is that they are not well suited for steep size distributions. Indeed repartitions where droplets are concentrated at an edge of the section cannot be efficiently evaluated since most of the nodes are not properly located. An alternative strategy [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF] is to choose adaptative node quadratures for which the abscissas are distributed as a function of the integrand profile.

When choosing the affine reconstruction, the source term computation is only done over the interval [S a , S b ] where the basis size function is positive as defined in Eq.( 8). The source terms have the generic form:

I p/q = S b S a S p/q (β k -α k ) (S b -S a ) S - (β k -α k ) (S b -S a ) S a + α k dS (22)
Thanks to the polynomial form of the integrand, an exact analytical integration is always possible and it is written:

I p/q = 1 P α k - (β k -α k ) (S b -S a ) S a (S P b -S P a ) + 1 Q (β k -α k ) (S b -S a ) (S Q b -S Q a ) (23)
where we define:

P = 1 + p q et Q = 2 + p q .

Sectional mass transfer numerical peculiarities

Eulerian MF models rely on the choice of a continuous discretization of the size phase space and on the integration of the semi-kinetic system (3) using a finite volume method on the sections. The evaporation process results in a time evolution of spray distribution function in the size phase space is consequently treated as exchange terms between the successive sections and exchange terms with the gaseous phase and the spray. In a one size moment MF approach, there is only one equation of mass concentration per section, which means that only mass fluxes have to be taken into account. Classical ODE integration schemes can be used to compute these mass exchange terms without any significant numerical issues. For the TSM-MF method, which ensure an improved flux information at the section boundaries, the moments n (k) and m (k) in each section should respect the realizability relation given by Eq.( 6). When solving the ODE system at a given location x as expressed in Eq.( 15), it appears that classical explicit time integrations do not preserve the realizability condition. Implicit methods applied on the mass transfer terms are still possible but they are complex and costly because of the important number of variable dependencies.

So we investigate a time integration method that accounts robustly for the mass transfers and its associated heat transfer. As highlighted in the description of the complex evaporation models in Appendix A, the heat transfer and the mass exchanges are strongly coupled to each other. In order to get a high-fidelity restitution of the gasspray couplings especially when combustion-thermoacoustic interactions are studied, all the gas and particle coupling source terms should preferably be computed using an all-at-once time integration scheme. Moreover we also want to solve accurately the momentum exchange. As a consequence, we propose time integration strategy for which the mass transfers, the drag and the heat exchange terms are solved simultaneously. A new scheme has been developed in order to deal efficiently with evaporation considering number and mass transport in the size phase space and including the drag and heat exchange source term resolution. The main issue is to correctly evaluate the number and mass concentrations that are respectively transfered from section (k + 1) to section (k) and from section (k) to section (k + 1) when evaporation or condensation occurs. For the sake of legibility, this new scheme is presented here in its upwind evaporating version for which exchanges occur from upper sections to the lower ones. The derivation of the scheme follows three steps:

1. Reconstruction of the distribution in each section at time t n 2. Computation of moments (n out , m out , n source , m source ) on "source" and "out" sub-sections at time t n ; the moment abscissas are determined based on a one point quadrature strategy.

Transport of the abscissas following the evaporation rate and computation of the moments at time t n+1

The scheme is an explicit time integration. It is developed in the spirit of other developpements presented in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] where the spray is described using moment methods for which size-dependent evaporation laws can be treated robustly. More details on the quadrature based method dedicated to kinetic equation are given in [START_REF] Laurent | Two-size moment multi-fluid model: a robust and high-fidelity description of polydisperse moderately dense evaporating sprays[END_REF]. As an asset of the new scheme, we present the strategy to integrate simultaneously the average sectional terms, which are the drag force and the heat transfer.

Derivation of the scheme

Several assumptions are done: we suppose an averaged size variation rate K (k) for each section and we consider it constant during a time step. It is deduced from Eq.(A.4), which links the mass transfer term to the droplet size variation rate:

K (k) = 4 √ π ρ l S k S k-1 S -1/2 ṁp (S ) κ (k) (S )dS S k S k-1 κ (k) (S )dS (24)
At time t n and for each section [S k-1 , S k ], we determine the size basis form function κ (k) of the distribution from the moment n (k) n and m (k) n . Then, as shown in Fig. (3), we compute fluxes that represent the number and mass concentration quantities that leave the k th section. For a sake of legibility, we define α = ρ l 6 √ π .

n (k) out = S k-1 +K (k) ∆t S k-1 κ (k) (S )dS m (k) out = α S k-1 +K (k) ∆t S k-1 S 3/2 κ (k) (S )dS (25)
Note that those integrals have the same form as the ones presented in Eqs. [START_REF] Dupays | Effects of particulate phase on vortex shedding driven oscillations in solid rocket motors[END_REF][START_REF] Gallier | Aluminum combustion driven instabilities in solid rocket motors[END_REF]. The mass quantities are difficult to integrate analytically when using the exponential reconstruction, as explained in Subsection 4.1. The affine approach offers the possibility to determine exact fluxes, which limits the numerical errors that can be introduced. Based on a quadrature strategy, we provide one quadrature point over an "out" sub-section [S k-1 , S k-1 + K (k) ∆t] and a "source" sub-section [S k-1 + K (k) ∆t, S k ]. The quadrature abscissas s (k) source and s (k+1) out that can be considered as averaged surfaces on the sub-sections are:

s (k) source =       m (k) source α n (k) source       2/3 s (k+1) out =        m (k+1) out α n (k+1) out        2/3 (26)
They are obtained from the moments (n (k) source , m (k) source ) and (n (k+1) out , m (k+1) out ) of the two sub-sections:

               n (k) source = n (k) n -n (k) out m (k) source = m (k) n -m (k) out n (k+1) out = S k +K (k+1) ∆t S k κ (k+1) (S )dS m (k+1) out = α S k +K (k+1) ∆t S k S 3/2 κ (k+1) (S )dS (27)
The moments n (k) n+1 and m (k) n+1 at time t n+1 = t n + ∆t are computed by:

           n (k) n+1 = n (k) source + n (k+1) out m (k) n+1 = m (k) source -ṁ(k) source ∆t + m (k+1) out -ṁ(k+1) out ∆t (28)
where ṁ(k) source and ṁ(k+1) out are the mass evaporation rates evaluated on the sub-sections. The following proposition gives the expression of the evaporation rates that ensures the realizability conditions given in Eq.( 6) for the moments n (k) n+1 and m (k) n+1 . Note that when k = N sec , all quantities namely "N sec + 1" do not exist.

Proposition 4.1. Let consider that the two size moments (n

(k) n ,m (k) n ) k∈[1,...,N sec ] for all size intervals [S k-1 , S k [∈ R + respect
the realizability conditions of Eq.( 6) at time t n . The moments (n (k) n+1 ,m (k) n+1 ) k∈[1,...,N sec ] computed at time t n+1 from Eq.( 28) ensure the realizibility relations:

α n (k) n S 3/2 k-1 < m (k) n+1 < α n (k) n S 3/2 k , k ∈ [1, ..., N sec ] ( 29
)
if the two following conditions are respected: (i) the CFL stability criterion: CFL evap ≡ max{CFL (k) evap , k = 1, ..., N sec } ≤ 1, with:

CFL (1) evap = K (1) ∆t (S 1 -S 0 ) ; CFL (k) evap = K (k) ∆t min ((S k-1 -S k-2 ), (S k -S k-1 )) , k ≥ 2 (30)
(ii) the ṁ(k) source and ṁ(k+1) out terms expressed in Eq.( 28) and defined on the sub-intervals introduced in Eq.( 27) are computed such that:

               ṁ(k) source ∆t = α n (k) source s (k) source 3/2 -s (k) source -K (k) ∆t 3/2 ṁ(k+1) out ∆t = α n (k+1) out s (k+1) out 3/2 -s (k+1) out -K (k+1) ∆t 3/2 (31)
Proof. Using Eq.( 28) and Eq.( 31), one can write:

m (k) n+1 α n (k) n+1 = n (k+1) out s (k+1) out -K (k+1) ∆t 3/2 + n (k) source s (k) source -K (k) ∆t 3/2 n (k+1) out + n (k) source (32)
Using Eq.( 26) and Eq.( 27) and when respecting the CFL condition K (k) ∆t ≤ (S k -S k-1 ), we have:

s (k) source -K (k) ∆t ∈]S k-1 , S k -K (k) ∆t[ ]S k-1 , S k -K (k) ∆t[ ∈ ]S k-1 , S k [ (33)
In the same way, respecting the CFL criterion K (k+1) ∆t < (S k -S k-1 ), we have:

s (k+1) out -K (k+1) ∆t ∈]S k -K (k+1) ∆t, S k [ ]S k -K (k+1) ∆t, S k [ ∈ ]S k-1 , S k [ (34)
Since n (k+1) out ≥ 0 and n (k) source ≥ 0, it finally proves the realizability criterion presented in Eq.( 29).

The drag and heat exchange source terms are evaluated in their corresponding "out" and "source" sub-sections:

                     F (k) out = α S k-1 +K (k) ∆t S k-1 F(S )S 3/2 κ (k) (S )dS F (k) source = α S k S k-1 +K (k) ∆t F(S )S 3/2 κ (k) (S )dS H (k) out = α S k-1 +K (k) ∆t S k-1 H(S )S 3/2 κ (k) (S )dS H (k) source = α S k S k-1 +K (k) ∆t H(S )S 3/2 κ (k) (S )dS (35)
Finally, we obtain:

                                                                                                                   ρ g Y vap n+1 = ρ g Y vap n + ∆t N sec -1 k=1 ṁ(k) source + ṁ(k+1) out ρ g u g n+1 = ρ g u g n -∆t N sec -1 k=1 F (k) source -ṁ(k) source u (k) n + F (k+1) out -ṁ(k+1) out u (k+1) n ρ g E g n+1 = ρ g E g n -∆t N sec k=1 H (k) source + F (k) source .u (k) n -ṁ(k) source h (k) source + L (k) v + 1 2 u (k) n 2 -∆t N sec -1 k=1 H (k+1) out + F (k+1) out .u (k+1) n -ṁ(k+1) out h (k+1) out + L (k+1) v + 1 2 u (k+1) n 2 n (k) n+1 = n (k) source + n (k+1) out m (k) n+1 = m (k) source -ṁ(k) source ∆t + m (k+1) out -ṁ(k+1) out ∆t m (k) n+1 u (k) n+1 = m (k) n u (k) n -m (k) out u (k) n + m (k+1) out u (k+1) n + ∆t F (k) source -ṁ(k) source u (k) n +∆t F (k+1) out -ṁ(k+1) out u (k+1) n m (k) n+1 e (k) n+1 = m (k) n e (k) n -m (k) out e (k) n + m (k+1) out e (k+1) n + ∆t H (k) source -ṁ(k) source e (k) n +∆t H (k+1) out -ṁ(k+1) out e (k+1) n                                                          k = 1, N sec (36)
The new scheme is presented in its one step explicit version, which matches to a first order scheme for the average sectional exchange terms such as the drag force one. In order to improve the accuracy resolution of those terms, high order explicit methods are possible using the same principle as the TVD Runge-Kutta schemes developed in [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF][START_REF] Gottlieb | Total variation diminishing runge-kutta schemes[END_REF]. We will use the second order or third order RK-TVD method as in [START_REF] Dupays | Effects of particulate phase on vortex shedding driven oscillations in solid rocket motors[END_REF]:

Q 1 = Q n + ∆t S (Q n ) Q n+1 = 1 2 Q n + 1 2 Q 1 + ∆t 2 S (Q 1 ) Q 1 = Q n + ∆t S (Q n ) Q 2 = 3 4 Q n + 1 4 Q 1 + ∆t 4 S (Q 1 ) Q n+1 = 1 3 Q n + 2 3 Q 2 + ∆t 2 S (Q 2 ) (37)
The main constraint is still to ensure the realizability conditions for the quantities (n (k) , m (k) ) k 's. As TVD-RK methods are linear combinations of integration steps that respect those conditions for each section, the global scheme is realizable. At each stage of the RK schemes, the size moments have changed, which means that updating the affine or exponential reconstructions in each section is required. So the computational gain obtained for the inversion step of the affine function becomes a key advantage as the inversion procedure is often called.

0D validation of the TSM-SRST scheme

In this subsection, our goal is to verify the new scheme for various evaporation laws. We propose to study the evolution of a smooth polydisperse size distribution, namely L. We represent the evolution of the NDF through its size moments of order 0 and 3/2. In this study, we choose a one-way approach for which the phase coupling terms with the gas are not considered and the evaporation laws do not depend on the evolution of the gaseous phase. So we use the Euler explicit version of the SRST scheme; the RK-TVD versions have been developed to increase the time order integration of the average drag force and heat transfer terms.

Analytical background for the study of evaporating polydisperse distributions

We first propose to test the scheme considering a constant evaporation law K p (S ) = cst. An analytic solution L(t, S ) can be established from the initial NDF L 0 (S ) = L(0, S ); this solution is a translation of L 0 at the constant speed K p in the surface space phase: L(t, S ) = L 0 (S + K p t).

With a constant evaporation law study, we are interested in the choice of the most appropriate size discretization in the same way that it was proposed in [START_REF] Laurent | Accurate treatment of size distribution effects in polydisperse spray diffusion flames: multi-fluid modelling, computations and experiments[END_REF] for the OSM-MF method. As presented in Eq.( 30), the so-called evaporating CFL evap criteria that are computed in each section have to be lower than 1 to ensure the stability of the scheme. For the specific case where all the CFL evap criteria are equal to one, the solution is exact. Nevertheless, as explained in Subsection 3.2, the time step used in SRM computations is chosen as regards convection and acoustic time scales, which are usually much faster than the evaporation ones. So low CFL evap values are commonly encountered. In this study, we focus on the impact on the numerical solution of the use of different time discretizations. Then, we provide results with a fixed time step and several number of sections. To perform this analysis, we focus on the temporal evolutions of the errors on the total bulk mass M tot and number concentration n tot , both quantities are defined in B.3. The error calculation relies on the difference between the total number concentration n tot (respectively the total bulk mass M tot ) at a time t computed with the new scheme to the ones obtained from the exact solution. Both differences are then divided by their initial values. So the errors are expressed in percentage. Note that in Appendix B, we perform a similar analysis with a square root evaporation law in order to extend the presented study and to validate the scheme for non constant size variation rates per section.

Numerical results for a constant evaporation law

We take K p = 1.10 -4 m 2 /s and a simulated period of 1.5 10 -4 s. Simulations are performed with two different size discretization strategies: a constant ∆S approach and a constant ∆r one. For each case, we choose a fixed number of 5 sections and we study the accuracy of the solution for three time steps: ∆t 1 = 2.10 -4 s, ∆t 2 = 2.10 -5 s and ∆t 3 = 2.10 -7 s, which correspond exactly for the constant ∆S approach to CFL evap values equal to CFL evap = 0.1, CFL evap = 0.01 and CFL evap = 1.10 -4 .

When comparing results of the two size discretization approaches, as illustrated in Fig.( 4) and Fig. [START_REF] Capecelatro | An Euler-Lagrange strategy for simulating particle-laden flows[END_REF], the constant ∆r strategy provides better accuracy on both bulk mass density and number concentration. It is particularly relevant for the number concentration for which the error at ∆t 2 is around 4% for the constant ∆r approach whereas it is about 7% for the constant ∆S one. The advantage of the constant radius strategy mainly stems from a more refined size discretization for the small droplets. This point will be further discussed for the square root study in Appendix B.

The accuracy level for both size moment evolutions decreases when using time steps featuring small evaporating CFL evap criteria. Nevertheless, we observe an error saturation, which means that the scheme converges to a solution for the small time steps. It is noticeable that the bulk density and number time evolutions are well respected with an interesting number of 5 sections. Considering the smallest time step ∆t 3 for both size discretizations, the total bulk density error is still inferior to 1% and the error on the total number concentration does exceed about 8%, which is acceptable. Indeed, having correct results on the bulk density evolution leads to a satisfying mass exchange dynamics with a gaseous phase. The errors presented for a fixed size discretization interval in Fig. [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF] and Fig. [START_REF] Capecelatro | An Euler-Lagrange strategy for simulating particle-laden flows[END_REF], with a time step tending to zero, thus a CFL number based on the evaporation rate going to zero, and compared to the exact solution, thus exhibiting a convergence toward the intrinsic error related to the size discretization step. Hence, it is natural to converge toward this intrinsic error, which is the one related to the chosen size discretization. In the fully unsteady simulations provided in Section 6, there is little hope to obtain a precise evaluation of the influence of these errors on the simulations. We have conducted simulations with a smaller size step and considered the simulations as correct when the influence of the size discretization become small enough. We finally perform three computations with 5, 10 and 20 sections and a constant ∆r discretization. We choose to use the ∆t 2 time step for which we show that the scheme is converged; the use of a smaller time step is not relevant. As presented in Fig. [START_REF] Castella | Splitting methods with complex times for parabolic equations[END_REF], we obtain better results on the number concentration with an error around 0.5% with 20 sections; the error on the bulk density is comparable.

As a conclusion, the new scheme provides satisfying results even when using small time steps. As expected, we have better time size moment evolutions when using a high number of sections. For the same number of sections, we improve the accuracy level of the solution when considering the ∆r size discretization, especially on the number concentration. Another validation performed with a square root evaporation law confirms the accuracy of the scheme and that the constant ∆r discretization is more appropriate to treat accurately the evaporation of the small droplets.

Linear evaporating two-phase acoustic study

In this section, we first introduce the context of linear evaporating two-phase flow acoustics, which offers an interesting validation framework for two-phase flow solvers that aim at studying unsteady two-way coupled phenomena. The presence of monodisperse evaporating sprays may have a driving impact on acoustic distrubances as explained in [START_REF] Dupays | Propagation of acoustic waves in a two-phase vaporizing mixture[END_REF]. In this context, we highlight the importance of taking into account size polydispersity on the acoustic response. Then, we derive for MF models a polydisperse dispersion relation that feature the behavior of an acoustic wave in a two-phase medium. Then, we evaluate the ability of MF models to tackle two-phase acoustics. Finally, we validate the time splitting strategy to account for polydisperse acoustics using relatively large time steps and we complete the validation of the SRST scheme as regards the mass, momentum and energy exchanges with the gaseous phase.

Analytical theory for evaporating two-phase flow acoustics

In this reference case, we consider a two-phase flow medium composed of a polydisperse spray that can exchange mass, momentum and energy with its own vapor. When at rest, the gas and liquid phases are at the same uniform temperature that is the saturated state of the droplets. Initially, the gas and the particles have no velocities, which finally means that there are no interactions between both phases. We are interested in finding an analytical solution to describe the propagation of plane wave disturbance in such a medium.

We investigate the case of a monochromatic wave of a small amplitude with a pulsation ω propagating in such a mixture. This unsteady acoustic perturbation creates a thermal disequilibrium between both phases, which results in alternate evaporating/condensing regimes about the saturation state put at T 0 . The mass and heat exchange model that is chosen to describe the two-way interactions is based on the thermal conductivity-controlled evaporation model presented in Appendix A.2. The decaying or amplifing harmonic solutions of the wave are presented with the following form on the gas variable:

u g (t, x) = sin(ωt -k 1 x) exp(-k 2 x)u 0 g ( 38
)
where we define the complex wave number k = k 1 +ik 2 and the one-phase wave number k 0 = ω/c 0 . We aim at deriving a dispersion relation that gives the spatial attenuation α and dispersion β coefficients that feature the acoustic wave behavior. We have α = 2k 2 and β = k 2 1 -(ω/c 0 ) 2 . The speed of sound in the gas alone is defined as c 0 = (γrT 0 ) 1/2 where γ is the isentropic coefficient and r is the specific constant of the gas.

Evaporating polydisperse dispersion relation for MF models

The liquid system describing the spray was based on the sampling approach presented in Section 2.5 for which a unique size of particles is used. Such a dispersion equation is here extended to polydisperse sprays using the continuous size approach that features the MF models. We start from the set of conservation equations of the semikinetic model developed in Eq.( 3) but written in a 1D context. Applying the several assumptions of MF models [HV1,2], [HT1,2] and then [HS1], we derive in Appendix C an expression of the complex wave number only function of the gas and liquid properties. The polydisperse dispersion relation derived from the MF formalism is written:

kc 0 ω 2 =         1 -γ c p,l (1 -λ(φ + (γ -1)/γ)) c p,g N sec k=1 µ k iωτ T k                 1 - c p,l (1 -φλ) c p,g N sec k=1 µ k iωτ T k                 1 + N sec k=1 µ k (1 -iωτ u k )         (39)
where we have defined:

λ = c p,g T 0 L v , φ = - ρ g,0 (γ -1) γρ l (40)
We have also introduced the mass ratio µ k = m (k) /ρ g,0 comparing the bulk mass density of the k th section to gas density i.e. the mass loading of each section. The term µ tot has been chosen to represent the total mass loading of the spray. The terms T 0 and ρ g,0 are respectively the temperature and the bulk density of the gas at the saturated state.

Differences between the OSM or TSM Multi-Fluid approaches are expressed through the droplet relaxation times that are computed for each section as a function of the size distribution. In Stokes' regime, we have:

τ u k = ρ l S k 18 π µ g τ T k = 3c p,l 2c p,g Pr τ u ( S k ), S k 2 = S k S k-1 S 3/2 κ (k) 0 (S ) dS S k S k-1 S 1/2 κ (k) 0 (S ) dS ( 41
)
where the terms κ (k) 0 represent the initial reconstruction functions in the sections. For the OSM-MF model, the basis size distributions are constant in the sections which means that the relaxation times only depend on the section boundaries. They are better evaluated in the TSM-MF method for which the basis size functions enable an accurate approximation of the mass repartition. In Section 5.4 and Appendix D, we provide qualitative comparisons between MF approaches, especially when having a limited number of sections.

We can also derive a dispersion relation that is continuous in size from the linearized semi-kinetic step given in Eq.(C.3). The polydisperse formula developed in the context of MF models converges for a large number of sections towards this continuous relation:

kc 0 ω 2 = 1 -γ c p,l (1 -λ(φ + (γ -1)/γ)) c p,g ρ l 6 ρ g,0 √ π R + S 3/2 n(S ) iωτ T (S ) dS 1 - c p,l (1 -φλ) c p,g ρ l 6 ρ g,0 √ π R + S 3/2 n(S ) iωτ T (S ) dS 1 + ρ l 6 ρ g,0 √ π R + S 3/2 n(S ) (1 -iωτ u (S )) dS (42)
For the rest of the study, we define the dimensionless attenuation α * and dispersion β * terms: α * = c 0 ω α and β * = c 0 ω β. Note that when α * > 0, the acoustic wave is attenuated. On the contrary, the perturbation is amplified when it has negative values. The dispersion term β * features the speed of the sound that decreases for positive values.

Importance of size polydispersity on the acoustic response

First of all, we want to determine the impact of size polydispersity on the acoustic response. To do so, we get the exact attenuation and dispersion from Eq.( 42) for two size distributions: a lognormal distribution and an equivalent monodisperse one that has the same average radius r 31 , which is the suitable average size to choose in Stokes' regime as explained in [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF]. Contrary to the lognormal function defined in Section 4.4, we choose to express the NDF here as a function of radius; it is namely here L r . Its radius r 31 is computed such as:

L r (r) = 1 r ln(σ ln ) √ 2π exp        -0.5 (ln r -ln(µ ln )) ln(σ ln ) 2        r 31 =                +∞ 0 r 3 L r (r) dr +∞ 0 r L r (r) dr                1/2 ( 43
)
We choose the following parameters for the lognormal function: σ ln = 1.9, µ ln = 3.10 -6 . This size distribution is studied on the finite interval r = [0 µm; 100 µm]. Both number and mass distributions are shown in Fig. [START_REF] Chedevergne | Instabilités intrinsèques des moteurs à propergol solide[END_REF]. The gas and spray properties are given in Tab. [START_REF] Abramzon | Droplet vaporization model for spray combustion calculations[END_REF]. As an additional study, we present results for two latent heats of vaporization L v,1 and L v,2 . Indeed, it has been demonstrated in [START_REF] Dupays | Propagation of acoustic waves in a two-phase vaporizing mixture[END_REF] that the amplification or the attenuation of the acoustic disturbance is very sensitive to this parameter. In Fig. [START_REF] Crouzeix | Une méthode multipas implicite-explicite pour l'approximation des équations d'évolution paraboliques[END_REF], we confirm that a sufficient low value for the latent heat of vaporization could amplify the acoustic wave, as illustrated by the negative α * coefficient for the L v,1 case. This shows that the mass transfer acts as a driving mechanism with this particular modeling. Focusing on the influence of size polydispersity, it is significant that the lognormal distribution and its monodisperse equivalent present numerous differences on the acoustic response. For the L v,1 case, both dispersion and attenuation curves of the lognormal polydispersion feature lower extrema that are around 40% inferior to the monodisperse distribution, which is quite important. Moreover, these extrema are located at different frequency values, especially for the L v,2 case for which the acoustic attenuation of the lognormal repartition reaches its maximum one decade lower than the monodisperse case. 5.0 10 6 J/kg L v,2 10.0 10 6 J/kg

Discussion on the ability of MF models to capture two-phase acoustics

We are here concerned about the ability of MF models to render the polydisperse nature of sprays in terms of acoustic behavior. Several parameters impact the way the MF method converges to the exact acoustic response of a given polydisperse distribution. There are two main criteria: the discretization of the size phase space and the choice of the size basis function in the sections. The size phase space can be discretized with constant radius or surface intervals. It is also relevant to share the size space into sections that contain the same amount of mass. Two discretization studies done in an inert context in [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF] and extended here for evaporating sprays in Appendix D, show that the constant mass discretization is well-suited. To highlight the impact of having an appropriate size reconstruction function in the sections, we compare the influence of using relaxation droplet times computed from a constant basis size function κ(r) = cst (OSM-MF model) to the exact characteristic times obtained from the exact d (k) 31 's of the exact size distribution function. It underlines that the droplet characteristic times are better estimated in the second case. So we can use a limited number of sections to have excellent results on the acoustic response. As illustrated in Fig. (D.25), we conclude that a constant mass discretization and the use of the exact size distribution are the best setting combination. As provided in Appendix D, 4 sections are enough whereas 8 sections are needed to obtain an accurate acoustic response with the constant size function κ(r) = cst. In pratical cases, the exact distribution approach is not possible so that we use the TSM-MF approach that ensures a better description than the constant size function κ(r) = cst. In order to validate the splitting technique with the SAP1 code, we use several time steps that are chosen as regards the acoustic time scale as it is expected for industrial applications and we illustrate the impact on the quality of the numerical solution. We propose the following unsteady numerical configuration: a 1D acoustic wave of a given frequency is injected in a non moving and uniform two-phase flow medium. The wave amplitude is sufficiently small in order to remain in the linear regime. When a sufficient number of period has propagated in the domain, we provide a post-processing of the instaneous 1D gaseous field to determine the attenuation α * and the dispersion β * of the acoustic disturbance. Note that the post-processing based on a local extremum tracking method is a crucial step since non negligeable errors may be introduced. So techniques such as interpolation have been chosen to decrease the errors due to the spatial discretization of the wave. The choice of the spatial discretization is done considering a spatial resolution of 70 points per wavelength (ppw), which ensures a limited amount of numerical diffusion and a better accuracy for the post-processing. So the domain length is set in function of the frequency of the injected acoustic wave and a minimum of 8 periods in the domain. In a preliminary approach, we perform a computation in a pure gaseous phase medium for which the theoretical attenuation and dispersion are zero. The preliminary acoustic single-phase simulation guarantees values of the coefficients α * and β * inferior to 10 -5 , which is quite satisfying. Then we study the case of a two-phase mixture; we consider a polydisperse distribution, which is the size lognormal size repartition presented in the previous Section 5.3. We propose two series of computations for the two different latent heats of vaporization. For each physical case, we perform several simulations for different frequencies. The physical parameters for the gas and the spray are the same as the ones used in the previous analytical study given in Tab. [START_REF] Abramzon | Droplet vaporization model for spray combustion calculations[END_REF].

We choose the two-stage splitting strategy introduced in Section 3.4.2 where the liquid convection is called less often than the so-called acoustic operator for which the gas convection and the source terms are integrated using a Strang splitting and a time step ∆t a . We use the RK2-TVD version of the SRST scheme and the disperse phase is described with the affine MF model. We choose to take 8 sections and a constant mass discretization of the size phase space as justified in the previous Section. Indeed 8 sections was enough to obtain a precise acoustic response with the OSM-MF model; so we expect that for the same number of sections, the TSM-MF model provides better results. Note that a ∆r = cst size discretization is strongly advised when dealing with droplet evaporation as explained in Section 4.4. This configuration is quite different in this study since the NDF distribution oscillates around its initial state depending on the evaporation or the condensation regimes. The simulated period T simu depends on the frequency case that is studied. For the 10Hz case, the domain length is 1500 m and the simulated time is 1.5s whereas the 10000Hz case implies a 1.5m domain and the simulated period of 1.5 10 -3 s. Finally, the splitting time step ∆t c chosen for the liquid convection is chosen equal to T simu /10. For the computations featuring a frequency superior to 500 Hz, the gas CFL criterion is very restrictive due to the small mesh cells, which imply to choose smaller time steps than for the inferior frequency cases. For the lower frequency values, we can test the splitting strategy with two time steps ∆t a,1 = 5.10 -5 and ∆t a,2 = 8.10 -5 s. For the ∆t a,1 cases, droplets with a radius inferior to 6 µm have relaxation times faster than the time step, which means that the gas-particle interactions are not correctly treated. For the ∆t a,2 case, the unsolved range of droplet sizes concerns the radius inferior to 7.53 µm. For a sake of legibility, one can find the time steps used for each simulation in Tab. [START_REF] Attili | Numerical simulation of multiphase flows in solid rocket motors[END_REF].

Numerical SAP1 results on dissipation and dispersion, as illustrated in Fig. [START_REF] Daniel | Eulerian approach for unsteady two-phase solid rocket flows with aluminum particles[END_REF] and Fig. [START_REF] De Chaisemartin | Eulerian models and numerical simulation of turbulent dispersion for polydisperse evaporating sprays[END_REF] for both L v,1 and L v,2 cases, are compared with the analytical solutions of the exact polydisperse acoustic responses. Depending on the value of the latent heat of vaporization, the opposite behaviors of the attenuation and the amplication of the acoustic wave are well restituted in both cases. For ∆t a,1 cases, quite satisfying agreement of the SAP1 results are to be underlined since the errors are smaller than 5% on both α * and β * relative to the extremum values.

When using the larger time step ∆t a,2 , we provide also very correct results, especially for the β * coefficient for which we obtain similar findings as the ∆t a,1 simulations. The main errors are noticed on the α * curve. Indeed, for the peak of amplitude of the L v,1 configuration at 50 Hz and for the maximum of attenuation of the L v,2 case that occurs around 100 Hz, there is an over-attenuation of the acoustic wave. The errors on the α * relative to the extremum values do not exceed 15%. Such an overdissipation has already been observed in [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF] for inert sprays when the splitting time step is larger than a certain range of droplet relaxation times.

This numerical study validates three main developments: we confirm the ability of TSM-MF model to handle evaporating polydisperse acoustic behavior, we extend the validation of the SRST scheme to momentum and energy exchanges with a gaseous phase whereas it was restricted only to the study of size moment evolutions in Part 4. [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF]. We finally highlight the potential of the splitting to treat unsteady reactive two-phase flows even using time steps that are superior to the fastest time scales of the smallest droplets. Depending on the accuracy level expected for industrial SRM simulations, the splitting technique confirms that it is an appropriate time strategy to control the cost/accuracy ratio. In the next Section, we propose to validate those methods on SRM configurations where we can take into account liquid convection, hydrodynamic structures and multi-dimensional aspects. We also test the developments with advanced droplet mass transfer models used in industrial computations.

Validation and feasibility on reactive two-coupled SRM simulations

In this section, we highlight the efficiency of both TSM-MF method and time splitting strategy by achieving 2D simulations of two SRM configurations with the industrial-oriented code CEDRE. Both computations are two-way coupled. The first is the so-called TEU built around a simplified cylindrical chamber geometry and a convergingdiverging nozzle. We first provide validation of the TSM-SRST scheme for evaporating droplets in an unsteady and multi-species gaseous flow. To do so, we compare quantitatively the results between the sampling approach and both TSM-MF models in order to highlight the advantage of the affine reconstruction compared to the exponential one in terms of computational cost saving. The second SRM configuration is the P230 motor, which presents unsteadiness flowfield with various types of hydrodynamic instabilities and various types of coupling between these instabilities, the acoustic modes of the chamber and the aluminum droplet combustion and with the inert aluminum oxide residues.

A previous SRM study [START_REF] Lupoglazoff | Numerical simulations of the unsteady flow inside segmented solid-propellant motors with burning aluminum particles[END_REF] with burning aluminum droplets described by a sampling spray approach was proposed but it was restricted to a monodisperse distribution. The smallest droplets, which were not inertial particles, were solved considering an equivalent gas approach as explained in [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF]. We extend and improve this study to a polydisperse case and we validate the affine MF method comparing to the sampling method. For this full-scale motor, the relevance of the splitting strategy is emphasized since the range of the time scales is wide and the droplet mass loading is high. In a recent two-phase SRM simulation [START_REF] Najjar | Simulations of solid-propellant rockets: Effects of aluminum droplet size distribution[END_REF][START_REF] Attili | Numerical simulation of multiphase flows in solid rocket motors[END_REF], a Lagrangian approach is proposed to describe the burning aluminum droplets and a convective-diffusive equation using an Eulerian formalism is chosen to treat the alumina residues that are restricted to the case of smoke particles. The following P230 computation reinforces the idea that a fully Eulerian formalism based on a MF method for the whole size spectrum is appropriate to follow complex burning droplet dynamics.

TEU test case for the splitting technique evaluation and Eulerian model comparisons

TEU numerical configuration

The TEU test case is a 2D axisymmetric simplified SRM configuration [START_REF] Dupays | Two-phase unsteady flow in solid rocket motors[END_REF]. It has been designed as a first approach to validate reactive two-phase flow solvers. This motor configuration is a pure numerical test case for which we choose to vaporize water droplets using the diffusion-controlled evaporation model presented in Appendix A.3. The simulation takes place on a deformed-structured 8 640 cell mesh. The domain is initialized with a pre-established two-phase flow computation featuring a steady state. A two-species gas is injected axially at 1000 K and at a constant surface flow rate of 39 kg/m 2 /s from the end head as seen in Fig. [START_REF] Descombes | Operator splitting for nonlinear reaction-diffusion systems with an entropic structure: singular perturbation and order reduction[END_REF]. The gas is composed of a dry mixture: Y air = 1.0 and Y vap = 0.0. Particles are also axially injected at the temperature of 300 K. The polydispersity is rendered with a lognormal size distribution L s , as introduced in Section 4.4. The NDF is defined over the finite interval d ∈ [ 0 µm ; 20 µm ]. We take the following parameters for the injected NDF distribution: σ ln = 0.4, µ ln = 2.79 10 -10 . The total droplet surface flow rate is 1.0 kg/m 2 /s, which then represents 2.5% of the total flow rate. Computations are performed with the sampling approach and with the affine and exponential MF models. We choose 5 sections/samples to offer a correct size refinement of the given size spectrum. A constant ∆r size discretization is chosen as advised in Section 4.4. The injection details are given in Tab. [START_REF] Billet | Impact of volume viscosity on a shock/hydrogen bubble interactions[END_REF]. We propose to study an unsteady state by pulsating the gaseous mass flow rate at a given frequency of 1500 Hz and an amplitude of 10% around the mean value: D g (t) = 39 (1 + 0.05 sin(2π f t)) kg/m 2 /s. To evaluate the accuracy level of the splitting strategy, we perform a numerical reference simulation for which most of all the droplet relaxation time scales are solved. For the sampling and affine MF model, two simulations are proposed using two time steps: the first computation called "reference" uses ∆t re f = 1.10 -7 s. The second computation called "fast" is performed with the time step ∆t f ast = 1.10 -6 s, which is superior to the acoustic time scale. So for all the computations, the gas convection is solved using an implicit second order Runge-Kutta scheme (RKI2) with a second order space resolution. The liquid transport is solved using an explicit second order RK (RK2) time integration method with a first order space scheme. For the source term treatment, a third time order integration scheme is chosen; the new SRST scheme is used for the MF methods as presented in Eq.( 37) for the MF methods. In order to integrate the different operators, we choose the Strang splitting introduced in Eq.( 16).

Eulerian model comparisons and splitting validations

We compare the three Eulerian models by verifying the evaporation dynamics of the polydisperse spray along the motor chamber. As shown in Fig. [START_REF] Descombes | Operator splitting for nonlinear reaction-diffusion systems with an entropic structure: singular perturbation and order reduction[END_REF], the enthalpy of the gas decreases significantly along the motor due to the important amount of energy required for the droplet evaporation. All methods feature very similar fields. In accordance with the pulsating gas injection, we observe sinusoidal flow structures for which the extrema have the same locations and the same amplitudes. As presented in Fig. [START_REF] Descombes | Operator splitting for nonlinear reaction-diffusion systems with an entropic structure: singular perturbation and order reduction[END_REF], a pressure sensor S has been used in the chamber in order to check that the frequency remains equal to 1500 Hz. We verify that the mean pressures are the same for all the methods; we have P mean = 2.21 10 5 Pa. The spray volume fraction fields are also quite comparable. This underlines that choosing 5 sections/samples is sufficient to obtain an equivalent evaporation dynamics. Differences are noticeable at the end of the evaporation, which occurs at the middle of the chamber. The exponential MF method presents a smoother gradient for the droplet volume fraction compared to the other approaches. Indeed, the exponential reconstruction leads to a more important diffusion in the size phase space during evaporation. On the contrary, the affine approach enables to have a sharp end of the queue of the distribution. In the end of the evaporation process, these differences are minor since the volume fraction value is low, which means that the retrocoupling on the gas is negligeable. As presented in Fig. [START_REF] Desjardin | Mechanistic model for aluminum particle ignition and combustion in air[END_REF], we also provide comparisons using two instantaneous cuts: an axial cut A-A and a radial cut R-R. We focus on the axial one, which highlights more precisely the gas temperature and volume fraction evolutions in the chamber. It presents minor differences between the sampling and the affine MF method. So it finally validate both new spray model and SRST scheme.

We are now interested in comparing the computational efficiency of the affine and exponential reconstructions. Comparing the wall clock times, the MF models are more expensive than the sampling approach. It is obvious since numerous intra-sectional and flux computations between sections have to be done. It is more relevant to investigate the competitiveness between MF methods. The exponential approach is about 4.1 times more expensive than the sampling method whereas it is less 2.3 times for the affine method. The difference is remarkably important since the WTC represents the computational time for both gas and liquid phase resolutions whereas the MF methods only differ in the choice of the sectional size description. In this study, we also show that the splitting technique reaches the same level of accuracy even when some time scales are not resolved. For the "ref" case, droplets with a diameter superior to 500 nm have relaxation times slower than the time step, which means that the gas-particle interactions are correctly treated. In the "fast" configuration, the time step ∆t f ast enables to only solve time scales for droplets with the diameter superior to 1.5 µm. We perform other computations using larger time steps but it finally leads to stability issues. They have been identified and they stem from the explicit resolution of the source terms. Now, regarding the results between the "fast" and "ref" simulations for both methods, we obtain very similar results as shown in Fig. [START_REF] Desjardin | Mechanistic model for aluminum particle ignition and combustion in air[END_REF]. The cut A-A provides no noticeable differences for the chamber conditions, which demonstrates that the evaporation dynamics is not impacted by the use of relative large time step. The radial cut R-R done in the nozzle features fast convection time scales and flow conditions resulting from the two-phase interactions that occur in the chamber. Differences on the gas temperature and axial velocity are observable but they are minor since they do not exceed 1%.

To conclude this unsteady SRM configuration, we show on the one hand that the splitting technique offers relevant accuracy levels for time steps that are larger than a certain range of droplet relaxation time scales and on the other hand, that the MF methods and the SRST scheme provide the suitable evaporation dynamics. The computational time is dramatically cut down with the affine approach, which is almost two times faster than the exponential one.

P230 configuration with burning aluminum particles

The P230 configuration is based on a full-scale internal motor geometry for which the propellant is divided into two grains containing aluminum particles. The injected size distribution is typically bimodal and burns near the propellant surface as the droplet sizes are small compared to the chamber size. The released energy in the combustion layer impacts the global acoustics of the chamber [START_REF] Radenac | Fluctuating energy balance to post-process multiphase flow computations[END_REF]. After burning, it only remains inert aluminum oxide residues, which results in a bimodal size distribution (see Appendix A.4) that can couple to hydrodynamic instabilities and the acoustics. Indeed, the aluminum oxide particles interact with two kinds of flow structures: the Parietal Vortex Shedding (VSP) commonly encountered in long motors with a high length-to-diameter ratio [START_REF] Chedevergne | Instabilités intrinsèques des moteurs à propergol solide[END_REF] and the Obstacle Vortex Shedding (VSO) due the presence of protruding thermal protections between the two propellant segments. For reasons of confidentiality, some values are not given and others are dimensionless.

Numerical configuration

The simulation takes place on a 2D axisymmetric deformed-structured 68, 500 cell mesh and the configuration is parallelized on 32 domains. At the simulated time, the forward star-shaped grain has totally burnt. The geometry is here considered fixed at a time when intense instabilities are observed. The simulation starts from a steady state of single-species gaseous flow field. Aluminum particles are injected from the two walls representing the two propellant grains as shown in Fig. [START_REF] Doisneau | Eulerian modeling and simulation of polydisperse moderately dense coalescing spray flows with nanometric-to-inertial droplets: application to Solid Rocket Motors[END_REF]. The droplet mass flow rate represents 18% of the total mass flow rate. As it is favorable to the sampling method, which was so far the unique spray description used for reactive computations in the CEDRE code, the distribution is represented here using two Dirac delta functions located at the diameter d 1 and d 2 with the following mass repartition: 12% and 6%. We provide two computations using the sampling approach and the affine MF model; all the other numerical parameters are the same. For the MF method, the SRST scheme has been modified in order to take into account the specificity of the aluminum droplet combustion model. Indeed, aluminum droplets burn until they reach a predetermined diameter d res as explained in Appendix A. It implies that the initial dirac size distribution evolves in a Lagrangian manner in the size phase space to another dirac distribution. This configuration appears as a tough validation test case for MF methods because we have to deal with discrete size distributions. So this study underlines the flexibility of the sectional spray description. For the SRST scheme, the fluxes between sections have been limited so that the final diameters d 1,res and d 2,res can be obtained. So we choose an appropriate size phase space discretization using 4 sections as given in Fig. [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF]. Note that this is a schematical illustration; we use the distributions, which are steeper and closer to the upper section boundaries in order to describe discrete repartitions. For the sampling method, we use two delta Dirac functions.

The simulation is solved on 3s time period, which is a sufficiently long period to provide a relevant frequency analysis featuring a 1/3 Hz precision to be compared, for instance with the 20Hz frequency of the first acoustic mode. We choose the time step equals to 1 µs, which permits to have an explicit time integration to solve the acoustics with respect to gas convection CFL criterion and to the given mesh. A RK2 scheme for gas convection is chosen with a second order space resolution. The liquid phase is also solved with an explicit RK2 method for the time integration and a second order space approach. For the source term computation, we use the RK2-TVD schemes. Total volume fraction Y Figure 16: Comparisons between the sampling and the affine MF method for the radial cut R-R in the combustion layer (Y=0 is the propellant surface). Results for the sampling method: --and for the affine MF method: -•-

Feasibility of the affine Multi-Fluid method

In order to evaluate the MF methods, we study the combustion dynamics and the impact of the aluminum oxide disperse phase on the gaseous field and on the instability levels. Note that the distributed combustion occurs in the very limited portion of the chamber near to the propellant called the combustion layer. So we propose a radial cut R-R as presented in Fig. [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF] in order to observe and to compare more precisely the sampling and MF methods. As illustrated in Fig. [START_REF] Duarte | Adaptive numerical methods in time and space for the simulation of multi-scale reaction fronts[END_REF], we obtain the same final temperature, which validates the modified SRST scheme as regards the energy release by the droplet combustion. Differences on the combustion dynamics can be noticed. When using the MF approach, the combustion layer is about 30% thicker than with the sampling method. We propose two main explanations: we describe with a continuous size approach a strictly bi-disperse distribution and we use only one section to follow the evolution of each Dirac delta function as presented in Fig. [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF]. Consequently, both lead to an important diffusion in the size phase space. Nevertheless, the differences on the gas temperature between both Eulerian methods remain acceptable since they do not exceed 15%, which is satisfying. One may use more sections to describe the size evolution of this discrete distribution but it will increase the computational cost. This configuration is self-sufficient in order to demonstrate the potential and the flexibility of the MF method and associated numerical developments in the context of burning polydisperse sprays. To be more representative of realistic cases, continuous bi-lognormal distributions should be used; in this context the MF methods are more competitive. Such a distribution is not proposed here since it requires another combustion model, which cannot be parametrized only by the initial and residual droplet sizes. Such a complex model is beyond the scope of the present contribution but developed in [START_REF] Sibra | Modeling and study of the evaporation and the combustion of droplets in solid rocket motors using an Eulerian Multi-Fluid approach[END_REF]. As an added bonus compared to the sampling approach, the MF methods offer the possibility to treat coalescence of the aluminum oxide residues, which may occur in the chamber.

The aluminum oxide residues are located in different portions of the motor as regards their sizes and their relaxation times. For the smallest droplets d 1,res , they follow quickly the gas disturbances whereas the d 2,res are more inertial so they go further in the chamber. Indeed, the size ratio d 2,res /d 1,res is equal to 20, which leads to a ratio of 400 between their relaxation times. One can notice that both residues are located where most of the hydrodynamic structures are present. As illustrated in Fig. [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF], the instantaneous volume fraction fields highlight differences on the droplet locations within the hydrodynamic vortices. As a consequence, the two droplet sizes impact the gaseous phase differently, which underlines the importance of capturing polydispersity.

To complete the validation of the MF method, we provide a qualitative study of the motor instabilities. We use two pressure sensors, S1 put at the front and S2 at the rear of the chamber as given in Fig. [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF]. The results on the mean pressures are quite satisfying. Considering the sampling method as a reference, the differences do not exceed 3%. We get good agreements on the instability levels, which are slightly superior for the sampling approach. As presented in Fig. [START_REF] Duarte | New resolution strategy for multi-scale reaction wave using time operator splitting, space adaptative multiresolution and dedicated high order implicit/explicit time integrators[END_REF], the main frequencies between both simulations are the same. Due to a lot of noise, the frequency peaks cannot be clearly distinguished, especially for the first modes. In a previous P230 computation [START_REF] Lupoglazoff | Numerical simulations of the unsteady flow inside segmented solid-propellant motors with burning aluminum particles[END_REF], the first acoustic mode was coupled to the VSP and VSO instabilities, which tends to reinforce the first harmonic level. It does not occur in this present simulation for which the numerical parameters differ in two points. A monodisperse distribution describing only the d 2 droplets was considered in [START_REF] Lupoglazoff | Numerical simulations of the unsteady flow inside segmented solid-propellant motors with burning aluminum particles[END_REF]. Moreover, the aluminum mass loading was limited to 10%. This finally demonstrates the impact of the disperse phase on instability frequencies and levels and the importance of well describing its size-conditioned spray dynamics when dealing with complex interactions between acoustics and various types of hydrodynamic structures.

Conclusion

In order to describe reactive polydisperse sprays, we provide a relevant Two Size Moment MF method based on affine size reconstructions. This fully Eulerian approach is particularly appropriate to accurately follow the size distribution evolution with a reasonable number of sections. The efficiency for the coupling source term computation is improved, which leads to a significant reduction of the CPU cost for complex two-phase flow simulations as highlighted by the TEU computation. The 0D validation demonstrates that the new SRST scheme is a robust method to solve sectional mass transfer peculiarities inherent of the TSM-MF methods. Both evaporating two-phase acoustic study and SRM computations validate widely the scheme, which is developed to solve simulaneously the gas-particle interactions and the sectional exchanges in the context of complex evaporation models. Both affine TSM-MF method and SRST scheme prove to be suitable candidates to predict precisely the physics of the aluminum combustion zone, which is conditioned by very specific mass/heat transfer models. A numerical coupling strategy based on a splitting technique is designed to solve unsteady two-phase flow problems for two-way coupled SRM simulations featuring Results for the sampling method ---: and for the MF method: -large range of time scales. As a flexible approach, it is successfully implemented in the advanced CFD platform CEDRE and in the research SAP1 code. The splitting strategy is robust and proves that time steps superior than the fastest time scales of the problem can be used, while preserving global accuracy of the resolution. So it offers an interesting cost/accuracy ratio that is considered as a key requirement for industrial computations. The splitting integration exhibits both efficiency and ease of use for the P230 simulation by capturing faithfully the interactions between droplet combustion, acoustics and flow structures. As a perspective, it is interesting to reach high fidelity description of the Al/Al 2 O 3 droplet combustion. In the spirit of the presented developments, a dedicated MF method for bi-component sprays is introduced in [START_REF] Sibra | A new eulerian multi-fluid model for bi-component polydisperse sprays: an essential approach to evaluate the impact of aluminum combustion on solid rocket motor instabilities[END_REF] but it needs to be implemented and validated in an industrial platform and for representative SRM configurations.

S . An analytic solution for the NDF resulting from a square root evaporation law can be obtained assuming a smooth initial size repartition. As presented in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF], the solution at time t is:

L(t, S ) = L 0 (Φ(t; 0, S )) K p (Φ(t; 0, S )) K p (S ) (B.1)
where Φ(t; 0, S ) is the solution at time t of Eq: d t Φ = K p (Φ) with Φ(0) = S . When taking the square root evaporation law K p (S ) = √ a max + S /S max where a max and S max are constants that we define later, one can get the solution :

Φ(t; 0, S ) = S max             t 2S max -a max + S S max       -a max       (B.2)
In this study, we consider an initial NDF that is a lognormal distribution function of the surface S , namely L s , over the finite interval S ∈ [ 0 ; S max ]:

M tot = n tot ρ l 6 √ π S max 0 S 3/2 L s (S ) dS , L s (S ) = 1 S ln(σ ln ) √ 2π exp        -0.5 (ln S -ln(µ ln )) ln(σ ln ) 2        (B.3)
where M tot is the total bulk mass density of the distribution and n tot represents the total number concentration of the distribution. In order to initialize the simulation, we provide a discretization into sections over the finite interval [ 0 ; S max ] and we determine in each section the two size moments from the exact size distribution L 0 s such that:

n (k) = n tot S k S k-1 L s (S ) dS ; m (k) = n 0 ρ l 6 √ π S k S k-1 S 3/2 L s (S ) dS (B.4)
Then we compare the number and mass moment time evolutions computed with the new scheme to the ones obtained from the exact solutions for both the constant evaporation law and the square root one Eq.(B.2). For both evaporation law studies, we choose the same parameters for the initial NDF: σ ln = 0.

        1 + N sec k=1 µ k γ h (1 -γφ) (1 -iωτ T k )         = c 0 (γ -1) ω k u (C.7)
For convenience, we have introduced: λ = c p,g T 0 L v , φ = -ρ g,0 (γ -1) γρ l We finally solve the previous analytical system to express the complex wave number as a function of the gas and liquid properties. It leads to the polydisperse dispersion relation:

kc 0 ω 2 =         1 -γ c p,l (1 -λ(φ + (γ -1)/γ)) c p,g N sec k=1 µ k iωτ T k                 1 - c p,l (1 -φλ) c p,g N sec k=1 µ k iωτ T k                 1 + N sec k=1 µ k (1 -iωτ u k )         (C.8)

Appendix D. SIZE DISCRETIZATION AND MF RECONSTRUCTION STRATEGIES TO RENDER EVAP-ORATING POLYDISPERSE ACOUSTICS

In order to reach an interesting cost/accuracy ratio, it is asked for MF methods to render accurately physical phenomena with a minimum number of sections. This is particularly observable when dealing with large size distributions featuring wide spectra in terms of Stokes numbers. Two parameters are significant to account for polydispersity with MF approaches: first the discretization technique that is applied on the size phase space. Second, the way the size distribution is presumed in each section. Such a study has been performed in the context of inert polydisperse acoustics in [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF]; it is here extended for two-phase acoustic interactions taking into account mass exchanges. As a typical size distribution, we choose the lognormal size repartition presented in 5.3 to discuss on the size discretization strategies and to demonstrate the advantage of having a high level description for the size repartition in each section. First, we propose to compare two size discretizations:

• The constant ∆r discretization technique which is featured by the choice of equal radius intervals.

• The constant ∆m discretization strategy for which the Multi-Fluid boundaries S k = 4πr 2 k are determined in order to have equal m (k) 's in each section: Then, we suggest two strategies to describe the size repartition on each section and to compute the average radii and the droplet characteristic times:

m (k) =
• The constant κ(r) distribution function: As it is done for the one size moment MF model, we describe the size repartition by choosing a constant size function κ (k) (r) on each section. As a consequence, the droplet relaxation times in the sections are computed from average radii r (k) 31 's obtained that are only a function of the fixed section boundaries: r (k) 31 = 0.5 r 2 k + r 2 k-1

1/2 .
• The exact size distribution function: We consider the exact expression of the size distribution function, which is a lognormal one for this study, to determine the exact average radii and relaxation times of each interval by using the relations Eq [START_REF] Schatzman | Towards non commutative numerical analysis: High order integration in time[END_REF]. The size discretization has been performed on a large interval r = [0 µm; 100 µm] which is quite sufficient to integrate faithfully the lognormal distribution. As a correction, the amount of mass that is contained on the interval r = [100 µm; +∞[ is added to the section of the biggest droplets. To highlight the influence of each parameter on sections with a ∆r = cst discretization strategy and using exact r (k) 31 in each section the two-phase acoustic response, we perform four results that cover all the combinations of the size discretization and description strategies. In each case, we provide a MF convergence study by comparing the acoustic response of 2, 4, 8 and 12 sections to the analytical polydisperse solution computed from Eq.( 42). First, we compare the results of the constant ∆r discretization to the constant ∆m discretization considering the constant κ(r) distribution in the sections; they are respectively given in Fig. D.22 and D.24. The constant ∆r strategy yields poor convergence since it still mislocates the extrema of dispersion and attenuation curves even with 12 sections whereas the constant ∆m technique is quite good at low numbers of sections. The amplitude of both curves are also overestimated for the constant ∆r discretization. The accuracy of the ∆m = cst approach is satisfying for a moderate number of sections, 8 intervals are sufficient to evaluate correctly the extremum values. As a conclusion, the constant ∆m appears appropriate to tackle polydisperse distributions that span wide spectra of Stokes number. Indeed it fashions better to the skewness of the mass repartition of the lognormal for the smallest droplets as represented in Fig. [START_REF] Chedevergne | Instabilités intrinsèques des moteurs à propergol solide[END_REF]. Second, we highlight the advantage of a high level of description in each section using the same discretization technique. For the constant ∆r discretization, the results of the κ(r) = cst distribution and those obtained with the exact profile of the size distribution function are shown in Figure D.22 and D.23. The exact approach demonstrates that a better evaluation of the relaxation times τ u k and τ T k in the sections decreases significally the error done the extremum locations. A slight overestimation of the extremum values is still observable at 12 sections but it is quite good compared to the κ(r) = cst distribution strategy. Then, we propose to compare again the two size distribution techniques but using the constant ∆m discretization. When comparing Figure D.24 and D.25, both strategies provide remarkably fast convergence in sections which is mainly due to the appropriate choice for the size discretization. Nevertheless the exact computation of the droplet relaxation times in the sections ensures a better accuracy on both dispersion and amplication responses with an interesting low number of sections. Using at least 4 sections seems very satisfying to render the polydisperse acoustic response of the lognormal distribution. As a conclusion, we highly recommand to use the constant ∆m discretization especially when the spectrum of Stokes number is wide and the mass is concentrated in particular intervals of the size phase space. Moreover it is noticeable that a relevant estimation of the mass repartition within the sections provides satisfying results on the droplet characteristic time computations. It leads to a faster convergence to the analytical polydisperse acoustic response for an equal number of sections. In this study, we compute the exact r (k) 31 's of the sections because the size distribution is known. In more general cases, the size distribution function is undetermined so that the Two Size Moment MF model is suitable to better describe the sectionnal mass repartition compared the one size moment MF approach presented here using a constant κ(r) strategy.

Figure 1 :

 1 Figure 1: Left to right: Sketch of size distribution representations for the sampling, exponential and affine Multi-Fluid models

Figure 2 :

 2 Figure 2: Various configurations for the affine reconstruction as regards the realizability conditions

4. 3 .

 3 New scheme for a Simultaneous Reactive Source Term computation: the TSM-SRST scheme 4.3.1. Quadrature based strategy

Figure 3 :

 3 Figure 3: Principle of the Simultaneous Reactive Source Term scheme for TSM-MF models (with affine reconstruction)

Figure 4 :Figure 5 :

 45 Figure 4: Time evolution of the errors on the total bulk density and the total number concentration relative to initial values. Configuration for 5 sections and a constant ∆S discretization: --: ∆t 1 = 2.10 -4 s or CFL = 0.1 •-•: ∆t 2 = 2.10 -5 s or CFL = 0.01 -: ∆t 3 = 2.10 -7 s or CFL = 1.10 -4

Figure 6 :

 6 Figure 6: Time evolution of the errors on the total bulk density and the total number concentration relative to initial values. Configuration for a constant ∆r discretization and ∆t 2 = 2.10 -5 s: -: 5 sections; •-• 10 sections; --: 20 sections

Figure 7 :

 7 Figure 7: LEFT: Number density functions for the lognormal distribution and the monodisperse average radius r 31 / RIGHT: Mass density function (MDF=4/3πr 3 NDF) for the lognormal size distribution: Dashed line;

Figure 8 : 2 5. 5 .

 825 Figure 8: Analytical dispersion β * (top) and attenuation α * (bottom) of an acoustic wave in two-phase media: Solid: Lognormal size distribution and Dashed: Equivalent r 31 monodisperse distribution. Results for two different latent heat of vaporizations: Left: L v,1 , Right: L v,2

Table 2 :

 2 Time steps used for the acoustic operator corresponding to the different frequency simulations 5 Hz ⇒ 200Hz ∆t a,1 = 5.10 -5 s ∆t a,2 = 8.10 -5 s 500 Hz ∆t a,1 = ∆t a,2 = 2.8 10 -5 s 1000 Hz ∆t a,1 = ∆t a,2 = 1.4 10 -5 s 5000 Hz ∆t a,1 = ∆t a,2 = 2.8 10 -6 s 10000 Hz ∆t a,1 = ∆t a,2 = 1.4 10 -6 s

Figure 9 :Figure 10 :

 910 Figure 9: Dispersion β * and attenuation α * of an acoustic wave in a polydisperse lognormal spray. Configuration for 8 sections and the L v,1 latent heat of vaporization: --: Analytical; ∆t a,1 = 5.10 -5 s; : ∆t a,2 = 8.10 -5 s; •: ∆t a,1 = ∆t a,2 as given in Tab.(2)

Figure 11 :Figure 12 :

 1112 Figure 11: Comparisons between the Eulerian models for the "fast" computations: Sampling (top); Affine MF (middle); Exponential MF (bottom). Instantaneous fields for the total enthalpy of the gaz (left) and for total droplet volume fraction (right).

Figure 13 :

 13 Figure 13: P230 geometry with 32 domain decomposition; gas and droplet parietal injections (arrows)

Figure 14 :

 14 Figure 14: Size phase space discretization for the MF method

Figure 15 :

 15 Figure 15: Results for the affine MF method: Instantaneous volume fraction of aluminum oxide residues for d 1,res (top) and for d 2,res (bottom)

Figure 17 :

 17 Figure 17: FFT (log-log diagram) of pressure signals given by the sensors S1 (left) and S2 (right). Results for the sampling method ---: and for the MF method: --

3 ,Figure B. 18 :Figure B. 19 :Figure B. 20 :

 3181920 Figure B.18: LEFT: Square root law discretization with 5 sections and a constant ∆S discretization / RIGHT: NDF time evolution at t = 0 s and t = 2.210 -3 s with the affine reconstruction: Dashed line; Analytical solution: Solid line

Figure D. 22 :

 22 Figure D.22: Dispersion β * and attenuation α * of an acoustic wave in evaporating polydisperse lognormal spray using parameters given in Tab.1 with L v = 5.10 -6 J/kg. Solid: Exact polydisperse response for the NDF-lognormal distribution and Dashed: MF response with 2, 4, 8 and 12 sections with a ∆r = cst discretization strategy and using r (k) 31 's computed with κ(r) = cst

Figure D. 23 :

 23 Figure D.23: Dispersion β * and attenuation α * of an acoustic wave in evaporating polydisperse lognormal spray using parameters given in Tab.1 with L v = 5.10 -6 J/kg. Solid: Exact polydisperse response for the NDF-lognormal distribution and Dashed: MF response with 2, 4, 8 and 12 sections with a ∆r = cst discretization strategy and using exact r (k) 31 in each section

Figure D. 24 :Figure D. 25 :

 2425 Figure D.24: Dispersion β * and attenuation α * of an acoustic wave in evaporating polydisperse lognormal spray using parameters given in Tab.1 with L v = 5.10 -6 J/kg. Solid: Exact polydisperse response for the NDF-lognormal distribution and Dashed: MF response with 2, 4, 8 and 12 sections with a ∆m = cst discretization strategy and using r (k) 31 's computed with κ(r) = cst

Table 1 :

 1 Gas and droplet properties

	Property Value	Property Value
	T 0	3498.8 K	µ tot	0.396
	ρ 0	3.78 kg/m 3	ρ l	1766 kg/m 3
	c p,g	2021.8 J/kg/K	c p,l	1375 J/kg/K
	µ g	8.855 10 -5 kg/m/s	Pr	0.8
	γ	1.23
	L v,1

Table 3 :

 3 Particle injection conditions for the TEU configuration

		Section ranges (µm) d 30 diameter (µm) Surface flow rate (kg/m 2 /s)
	Sample/Section 1	[0.0, 4.0[	3.84	5. 10 -7
	Sample/Section 2	[4.0, 8.0[	7.23	≈ 7.88 10 -2
	Sample/Section 3	[8.0, 12.0[	9.85	≈ 6.512 10 -1
	Sample/Section 4	[12.0, 16.0[	13.20	≈ 2.512 10 -1
	Sample/Section 5	[16.0, 20.0[	17.01	≈ 1.95 10 -2

Table 4 :

 4 Comparison between Eulerian methods on instability amplitudes and dimensionless mean pressures

		Sensor 1 Sensor 2
	P mean sampling	100.00	100.00
	P mean affine MF	99.97	99.98
	∆P/P mean sampling	2.4%	3.9%
	∆P/P mean affine MF	2.0%	3.7%

Acknowledgments

The present research was done thanks to a Ph. D. Grant from ONERA. This work and its associated developments come from a collaboration work between the applied mathematic team at EM2C laboratory and the Fundamental and Applied Energetics -DEFA department at ONERA.

Appendix A. DRAG, HEAT AND MASS TRANSFER MODELS AT THE LEVEL OF A SINGLE DROPLET

In this Appendix, we give details of the drag, heat and mass models used in this study. We particularly focus on different heat transfer/evaporation models. In these models, the heat transfer is differently coupled to the mass exchange. Indeed the two first models rely respectively on a thermal-conductivity approach and on a species diffusion-controlled process that account for mass exchanges with single-species and multi-species gas flowfields. Both models do not take into account internal conduction as it is done in [START_REF] Abramzon | Droplet vaporization model for spray combustion calculations[END_REF][START_REF] Sirignano | Fluid Dynamics and Transport of Droplets and Sprays[END_REF] since having a temperature profile within the droplet is not compatible with the kinetic description of the spray and hard to deal with in terms of a statistical approach as explained 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] where a comprehensive study is provided. In order to take into account the droplet heating, an infinite liquid-conductivity assumption is made for which the droplet temperature is uniform inside the particle, equal to the temperature at its surface. In solid propulsion, this hypothesis is valid since we study relative small droplets and because liquid aluminum has a high thermal conductivity. Finally, we present a simplifed combustion model for Al/Al 2 O 3 droplets, which we will be used in the unsteady SRM computation of Section 6.2.

Appendix A.1. Drag force model

The classical relation for the drag force applied on a spherical particle is:

where ρ g is the local gas density and C D is the drag coefficient which is expressed differently in function of the particle Reynolds number Re p defined by Re p = ρ g d u gu p /µ g where µ g is the gas dynamic viscosity. When Re p < 1, the Stokes' regime is considered. At higher Reynolds numbers such as in solid rocket motor nozzles, one can consider Schiller-Naumann convective corrections [START_REF] Schiller | A drag coefficient correlation[END_REF]. Whatever regimes are considered, a dynamic characteristic time τ u can be defined in order to represent the time required for a particle to respond to a change in gas velocity. In Stokes' regime, the drag force and the droplet relaxation time are written:

One can distinguish droplets comparing their relaxation times to a gas characteristic time [START_REF] De Chaisemartin | Eulerian models and numerical simulation of turbulent dispersion for polydisperse evaporating sprays[END_REF][START_REF] Massot | Eulerian Multi-Fluid models: modeling and numerical methods[END_REF]. This ratio defines the dynamic Stokes number St. When St > 1, we consider inertial droplets. The particles are respectively moderately inertial and non-inertial when the Stokes number is close to one and much below one. Note that the droplet dynamics is strongly conditionned by its size. Indeed, as it is expressed in Stokes' regime, the relaxation time τ u depends on the square of the droplet diameter, which means that a large spectrum of droplet sizes leads to a wide spectrum of Stokes numbers. This case is often encountered in solid propulsion since the polydispersion of droplets in terms of size is important due to evaporation or other interaction phenomena such as coalescence. Note that for higher Reynolds numbers, the convective correction increases the drag transfer rate and the dynamic droplet times becomes smaller.

Appendix A.2. Heat transfer and thermal conductivity-controlled evaporation model When studying droplets that exchange mass with a single-species gaseous phase, the evaporation process is not controlled by the diffusion of species. Only the heat transfer drives the mass exchange between the particle and the gas. As presented in [START_REF] Williams | Combustion Theory (Combustion Science and Engineering Series)[END_REF], the mass transfer results from a thermal disequilibrium between phases. Indeed no evaporation occurs when the two phases are at the same temperature nor when the droplet has a temperature lower than its saturation one T sat . In this particular case, the saturation temperature corresponds exactly to the boiling point of the liquid. In this model, the mass transfer rate depends on the two phase temperature differences expressed through the thermal Spalding number B T . The mass transfer term ṁp for a single droplet reads:

where c p,g is the specific pressure heat capacity of the gas, λ g is the thermal conductivity of the gas, L v is the latent heat of vaporization and Pr is the Prandtl number. The term Nu c is the Nusselt number and it represents the convective correction term. The size variation rate K p that is used in Eq.( 1) can be directly deduced from the mass transfer term rate since it matches to mass flux ṁp released by the droplet surface. We have:

The heat exchange term is composed of two contributions: the energy needed to the phase change ṁp L v and the heat absorbed by the liquid core Q p . The heat transfer term is written:

Note that, when mass transfer occurs, the two contributions balance each other. All the heat received Q p is used for phase change ṁp L v and the particle temperature remains constant equal to T sat . As done before for the drag response time, a thermal characteristic time τ T can be obtained. In Stokes' regime, the heat received by the droplet is written:

where c p,p is the specific heat capacity of the liquid core. Note that the dynamic and thermal relaxation times have the same order of magnitude.

Appendix A.3. Diffusion-controlled evaporation model

In more general cases, the droplet exchanges mass with a multi-species gaseous phase. The evaporation dynamics is controlled by the diffusion of the species, which is represented by the mass Spalding number B M . As presented in [START_REF] Sirignano | Fluid Dynamics and Transport of Droplets and Sprays[END_REF], the mass transfer term is written:

where Y vap sur f and Y vap ∞ are respectively the vapor mass fractions at the droplet surface and far from the droplet. Compared to the single-species model, the evaporation occurs even if the droplet has not reached its saturation temperature. In this regime, the terms Q p and ṁp L v are not equal and the heat exchange H p defined in Eq.(A.5) is positive; this enables to increase the droplet temperature. In this model, Q p and the thermal Spalding number B T are functions of the evaporation rate that leads to a strong interaction between the mass transfers and the droplet heating:

The droplet stops heating when it is at its saturation temperature. Note that T sat is a function of the pressure and of the composition of the gas mixture that surrounds the droplet.

Appendix A. [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF]. Aluminum droplet combustion model with a single-species gas We propose to use the simplified model for aluminum combustion presented in [START_REF] Lupoglazoff | Numerical simulations of the unsteady flow inside segmented solid-propellant motors with burning aluminum particles[END_REF][START_REF] Gallier | Aluminum combustion driven instabilities in solid rocket motors[END_REF] that is supposed to follow a d 2 law. The released single-species by the burning droplet is the same as the surrounding gas, which is considered as a gas product. So the mass transfer term ṁp has an expression close to the one of the thermal-conductivity-controlled model but it takes into account the heat per unit mass Q r released by the droplet combustion:

The heat transfer term H p is the same as expressed in Eq.(A.5). As explained in the thermal conductivity-controlled evaporation model, there is no mass transfer until the droplet reaches its saturation temperature. Then, when droplet starts burning, its temperature remains constant. Finally, the gas temperature increases during the combustion thanks to an appropriate definition of the enthalpies of the liquid aluminum and of the gas product ensuring a correct amount of energy to the gas. So we determine the heat released by the droplet combustion Q r as a constant for the whole computation.

As an additional remark, we want to take into account the bi-component nature of the droplet [START_REF] Yetter | Metal particle combustion and nanotechnology[END_REF][START_REF] Sibra | A new eulerian multi-fluid model for bi-component polydisperse sprays: an essential approach to evaluate the impact of aluminum combustion on solid rocket motor instabilities[END_REF]. We consider the presence of nonreacting aluminum oxide that is present in the droplet when it leaves the propellant surface. In this model, the combustion starts when the droplet is at its saturation temperature T sat . The mass of aluminum oxide remains constant during the droplet burning. The combustion is finally stopped when the droplet size reaches a predetermined diameter d res that corresponds to the remaining mass of inert aluminum oxide. More representative combustion models [START_REF] Daniel | Eulerian approach for unsteady two-phase solid rocket flows with aluminum particles[END_REF][START_REF] Desjardin | Mechanistic model for aluminum particle ignition and combustion in air[END_REF] can be found to treat the ignition of the droplets and the deposition of condensed aluminum oxides on the droplet surface. Nevertheless, the presented model is quite relevant to validate our developments accounting for a fully Eulerian description of burning polydisperse sprays.

Appendix B. VALIDATION OF THE TSM-SRST SCHEME: Square root evaporation law

Appendix B.1. Analytical background for the square root law Then we focus on an evaporation law that is representative of physical evaporating models. Indeed the size evaporation rate has a regular evolution in term of droplet size through various evolutions of the droplet velocity and temperature. So we consider a smooth evaporation law with a size evolution rate K p that depends on the size variable

Appendix B.2. Numerical results

We choose a max = 1.10 -8 and S max = 0.6; the average size variation rates of the sections are determined such that:

For both size discretizations, we perform several computations with the same time step equal to ∆t = 2.10 -5 s and different number of sections (5, 10 and 20 sections). We first want to validate one of the main assumption of the SRST scheme that is a constant evaporation rate per section.

In ), the evaporating rates between the different sections features a ratio: K max /K min > 1.5, which is high enough to test harshly the SRST scheme. The biggest droplets evaporate faster than the smallest ones and we expect that the parameter σ ln of the distribution tends to decrease. For both size discretization approaches, the new scheme provides the appropriate time evolution of NDF, which tends to prove its ability to deal with complex evaporation law.

Comparing .20) where we propose an instantaneous view of the NDF evolution at t = 2.2 10 -3 s. With the ∆r = cst approach, the NDF is better discretized between [0; 2.10 -7 m 2 ] with 2 sections whereas only one section is used for the other size discretization. Moreover, the size evaporation rate is overestimated in the first section of the constant ∆S discretization. Focusing on the number of sections, we obtain very satisfying results for both discretizations with 10 and 20 section simulations that feature errors inferior to 5% for the number concentration and around 0.5% for the bulk density. Nevertheless, the results of the ∆r = cst approach are better for the bulk density time evolution, especially for the 5 section case. Finally, the use of the constant ∆r discretization is more appropriate to treat accurately the evaporation of the small droplets.

Appendix C. DERIVATION OF THE POLYDISPERSE DISPERSION RELATION FOR MF MODELS

The two-phase interaction terms are modeled in Stokes' regime such as: where τ u , τ T and τ m are respectively the dynamical, thermal and evaporating droplet relaxation times. The liquid properties c p,l and L v are assumed constant whatever the droplet size. Because of the linear acoustic framework, we linearize the different variables around their constant values:

We get the linearized semi-kinetic system:

Now we proceed the continuous discretization of the size phase space using sections defined as [S k-1 , S k [ and on which we suppose that the droplet velocity and temperature are constant: u p (S ) = u (k) and T p (S ) = T (k) . When considering a Two Size Moment approach, one can write the same linearized formalism for the number concentration of droplets and the bulk mass density per section:

The linearized system of conservation equations for the k th section is:

The averaged surface S k computed on each section k is derived from the general formula: S k p-q = S k S k-1 S p/2 n 0 (S ) dS

S q/2 n 0 (S ) dS In Stokes' regime, it can be shown that p = 3 and q = 1 are relevant. For a sake of legibility, we use the following formalism for the droplet relaxation times of each section:

The same strategy of linearization is applied on the Euler gas system of conservation equations. The gas source terms are established by summing the contribution of all the sections:

where ρ g,0 represents the bulk density of the gas at the saturated state. Then, we take adimensionnalized variables such as: ρ = ρ g /ρ g,0 , u = u g /c 0 , T = T g /T 0 , u (k) = u (k) /c 0 , T (k) = T (k) /T 0 . We also consider the mass ratio µ k = m (k) /ρ g,0 comparing to the bulk mass density of the k th section to gas density i.e. the mass loading of each section. We also introduce the term µ tot as the total mass loading of the spray. Assuming a harmonic regime Φ = Φ 0 exp(i(kx -ωt)), a system of analytical relations can be obtained by suppressing the temporal and spatial variations: