F Doisneau

F Laurent

A Murrone
email: angelo.murrone@onera.fr

J Dupays
email: joel.dupays@onera.fr

M Massot
email: marc.massot@ecp.fr

Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion

Keywords: Polydisperse spray, High order Eulerian Multi-Fluid model, Adaptive quadrature for coalescence integrals, Solid propellant combustion, Aluminum oxide droplets, CEDRE code

The accurate simulation of polydisperse sprays undergoing coalescence in unsteady gaseous flows is a crucial issue. In solid rocket motors, the internal flow depends strongly on the alumina droplet size distribution, which spreads up with coalescence. Yet solving for unsteady two-phase flows with high accuracy on the droplet sizes is a challenge for both modeling and scientific computing. As an alternative to Lagrangian approaches, a wide range of Eulerian models have been recently developed to describe the disperse liquid phase at a lower cost, with an easier coupling to the gaseous phase and with massively parallel codes. Among these models, the Multi-Fluid model allows the detailed description of polydispersity and size/velocity correlations by separately solving fluids of size-sorted droplets, the so-called sections. The existing One Size Moment method, which describes the size distribution with one size moment per section, provides simple and fast resolution for coalescence. On the other hand, a Two Size Moment method has been suggested to reduce the number of sections but it lacks an efficient coalescence resolution method. After introducing a new strategy for two size moment coalescence, the two methods are compared on various configurations in a research code and an industrial-oriented code, in order to conclude on computational accuracy and cost. Then the paper aims at describing the most efficient approach for multi-dimensional unsteady and eventually coalescing rocket chamber simulations. Its objective is threefold : first, to validate the Two Size Moment method by comparing simulations to reference solutions and dedicated experimental measurements conducted at ONERA, second to study the efficiency and robustness of both methods, third, to draw some firm conclusions about the necessity to use the One Size Moment or Two Size Moment method to simulate solid propellant alumina sprays. We finally perform the first simulations of coalescence in realistic 2D boosters with the Two Size Moment method, implemented in the industrial-oriented code CEDRE.

Introduction

Two-phase flows constituted of a gaseous phase carrying a disperse condensed phase play a key role in many industrial and scientific applications : spray evaporation and combustion in Diesel engine combustion chambers, fluidized beds, dynamics of planet formation in solar nebulae, etc. In all these applications the disperse phase is composed of particles of various sizes that can eventually coalesce or aggregate, break-up, evaporate and have their own inertia and size-conditioned dynamics. So the importance of polydispersity is obvious for a comprehensive modeling of these phenomena.

In solid rocket motors (SRM), aluminum powder is frequently used as an additive to the solid propellant to increase the specific impulse of the engine. Unlike the other ingredients, aluminum particles can burn in a significant portion of the chamber and finally produce a high mass fraction of aluminum oxide -or alumina-in a condensed liquid disperse phase. This disperse phase encounters drag forces, coalescence and heat exchanges [START_REF] Dupays | Some aspects of two-phase flows in solidpropellant rocket motors[END_REF]. In the nozzle, the droplets accelerate suddenly and cool down with the gas, eventually becoming solid and breaking up in such velocity gradients [START_REF] Dukowicz | A particle-fluid numerical model for liquid sprays[END_REF][START_REF] O'rourke | Collective drop effects on vaporizing liquid sprays[END_REF][START_REF] Amsden | Kiva II, a computer program for chemically reactive flows with sprays[END_REF]. Thus, the disperse phase strongly interacts with the gaseous flow field during its way throughout the engine. It contributes to the performance loss of the SRM, first via a decrease in nozzle efficiency which penalizes specific impulse; second by modifying the acoustic response of the system and therefore having a complex influence on combustion and flow instability levels and frequencies; third some of the droplets, mainly the ones with high inertia, which end up in the eventual aft-dome region around the submerged nozzle, induce sloshing motion of this molten liquid slag and can lead to control problems and possible vehicle instability; fourth droplets are the source of slag material that may remain in the engine during firing, causing insulation erosion in high concentration zones. In such harsh conditions of pressure, temperature and velocity, solid propulsion experiments consume high technology materials and offer poor measurement output, especially on the disperse phase. The abundance of physical phenomena involved makes the models difficult to scale [START_REF] Tòth | Two-phase flow investigation in a cold-gas solid rocket motor model through the study of the slag accumulation process[END_REF]. Regarding the prohibitive cost of experiments, numerical simulation is the only available tool for optimizing rocket engines. Until now, complete 3D computations were achieved at the cost of drastic physics simplifications, allowing little self-reliance towards validation experiments. Yet time has come for comprehensive simulations, including advanced gas/droplet coupled models, to give predictive answers.

Focusing on the dynamics of the chemically inert alumina cloud in the motor is a first step to evaluate specific impulse loss in the nozzle, slag material accumulation and acoustic response of the system. Therefore solid propellant regression, aluminum particle combustion or fluid-structure interaction are here neglected. We work with non-evaporating sprays throughout the paper, keeping in mind the broad application fields related to the present study. By spray, we denote a disperse liquid phase constituted of droplets carried by a gaseous phase. We consider the specific case of what we refer to as moderately dense sprays, corresponding to hypotheses [START_REF] O'rourke | Collective drop effects on vaporizing liquid sprays[END_REF][START_REF] Mathiaud | Étude de systèmes de type gaz-particules[END_REF] that are realistic in a SRM : 1. the liquid mass fraction is high enough to generate strong drag and heat retro-coupling effects on the gas; 2. the liquid volume fraction is high enough to get significant effects from droplet collisions and coalescence; 3. but the liquid volume fraction is still much smaller than one so that the gas can be solved as a single phase flow with only source terms to ensure the coupling with the liquid phase. This level of coupling is sometimes referred to as four-way coupling to highlight the fact that coalescence yields particle-particle interactions but we rather call it two-way coupling hereafter to insist on the numerical peculiarities yielded by these particle retro-couplings on the gas. The study takes place in the combustion chamber, where the main physical processes that must be accounted for are : transport in real space, acceleration of droplets due to drag, conditioned by size, and coalescence, leading to polydispersity. These processes are remarkably sensitive to size distribution. On the one hand, size distribution impacts directly on droplet velocities through differences in the drag force relaxation times. On the other hand, velocity differences induce coalescence which impacts on polydispersity. So droplet size and velocity are intimately interlocked, which will be refered to as the "size/velocity coupling", illustrated in Fig. 1. We will therefore choose a model which is accurate as regards the size distribution. The retained approach called "kinetic", in reference to the kinetic theory of gases, describes the droplets as a cloud of point particles for which the exchanges of mass, momentum and heat are described using a statistical point of view, with eventual correlations : a finite set of global properties such as one size parameter (enough if droplets are spherical), the velocity of the center of mass, a temperature are explicitly included so that the total phase space is usually high-dimensional. More details about the droplets, such as angular momentum, non-sphericity factors, can be predicted by increasing the size of the phase space : it is established that refined droplet models can be used as long as they do not include history terms [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]. The Williams-Boltzmann equation, a transport equation based on kinetic theory, has proven to be useful for treating polydisperse, dilute and moderately dense liquid sprays [START_REF] Williams | Spray combustion and atomization[END_REF][START_REF] Williams | Combustion Theory (Combustion Science and Engineering Series)[END_REF]. Such an equation describes the evolution of the number density function or NDF of the spray due to the drag force of the gaseous phase and the droplet-droplet interaction of coalescence [START_REF] Kuentzmann | Aérothermochimie des suspensions[END_REF][START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF][START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. The statistical description of the disperse phase coupled to fluid equations such as Euler or Navier-Stokes for the carrier phase is a mesoscopic scale description, which will be the background of our contribution as a whole. The resolution of these fluid-kinetic coupled equations will be referred to as Direct Numerical Simulation (DNS), even if some level of modeling is involved and the droplet (microscopic) scales are not resolved. Besides, we will not introduce any turbulence modeling, such as turbulent dispersion models [START_REF] Simoes | Modélisation eulérienne de la phase dispersée dans les moteurs à propergol solide, avec prise en compte de la pression particulaire[END_REF] and sub-grid scale closure terms for both the gas and the liquid phases, as it is required in the context of Large Eddy simulation (LES) [START_REF] Wunsch | Theoretical and numerical study of collision and coalescence -Statistical modeling approaches in gasdroplet turbulent flows[END_REF][START_REF] Vié | Simulation aux grandes échelles d'écoulements diphasiques turbulents à phase liquide dispersée[END_REF][START_REF] Boileau | Robust numerical schemes for Eulerian spray DNS and LES in two-phase turbulent flows[END_REF][START_REF] Masi | Étude théorique et numérique de la modélisation instationnaire des écoulements turbulents anisothermes gaz-particules par une approche Euler-Euler[END_REF]]. Yet, efficient DNS is an incontrovertible starting point for disperse two-phase flow simulations while the adaptation of the model and methods to LES relies most of the time on theory and DNS tools. Thus, the framework of the present contribution will be the one of DNS of this fluid-kinetic coupled system of equations, even if the final simulations proposed in the context of the industrial-oriented code CEDRE will only be under-resolved and presented for feasibility purposes as a starting point for larger scale simulations and LES developments.

Size Distribution

There are several strategies in order to solve the liquid phase and the major challenge in numerical simulations is to account for the strong coupling between all the involved processes. A first choice is to approximate the NDF by a sample of discrete numerical parcels of particles of various sizes through a Lagrangian-Monte-Carlo approach [START_REF] Dukowicz | A particle-fluid numerical model for liquid sprays[END_REF][START_REF] O'rourke | Collective drop effects on vaporizing liquid sprays[END_REF][START_REF] Amsden | Kiva II, a computer program for chemically reactive flows with sprays[END_REF][START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF]. It is called Direct Simulation Monte-Carlo method (DSMC) in [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF] and is generally considered to be the most accurate for solving the Williams-Boltzmann equation; it is specially suited for direct numerical simulations (DNS) since it does not introduce any numerical diffusion, the particle trajectories in phase space being very accurately resolved. Its main drawback is the coupling of a Eulerian description for the gaseous phase to a Lagrangian description of the disperse phase, thus introducing spatial diffusion via projection steps on the Eulerian grid and offering limited possibilities of vectorization/parallelization or implicitation. Moreover for unsteady computations of polydisperse sprays, a large number of parcels in each cell of the computational domain is generally needed, thus yielding large memory requirement and CPU cost.

As an alternative, the Eulerian Multi-Fluid model, furthered 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] from the ideas developed in [START_REF] Greenberg | On the origin of spray sectional conservation equations[END_REF], relies on the derivation of a "semi-kinetic" model from the Williams equation using a moment method for velocity, but keeping the continuous size distribution function. This distribution function is then discretized using a "finite volume approach" in size phase space that yields conservation equations for mass, momentum (and eventually other properties such as number, energy) of droplets belonging to fixed size intervals called sections. Each of them constitutes a different "fluid". The two-way coupling, which occurs with moderately dense sprays, is natural to take into account in the Eulerian formalism, thanks to the similarity between the carrier fluid and the section systems of equations. It is well rendered as long as numerical methods are carefully designed, this question being tackled in the context of splitting methods in [START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF][START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF]. A total correlation of size and velocities is assumed at the semi-kinetic level and correlations are enforced at the Multi-Fluid level with a (usually constant) velocity closure in each section, which yet introduces little error if inertia and polydispersity are not too strong, these aspects being discussed in the present paper. Integrating on a continuous size variable in each section is a key aspect : while most Eulerian approaches, based on size sampling, consider discrete droplet sizes gathered into "classes" which cannot account for the new droplet sizes created by coalescence (except scarce examples [START_REF] Vasenin | Two-phase flows in the nozzles of solid rocket motors[END_REF]), continuous size approaches such as the Multi-Fluid method hereafter described are the only Eulerian methods handling coalescence naturally and rigorously. Some Eulerian approaches based on quadratures of the whole size support can capture the trend of the size distribution evolution due to coalescence [START_REF] Fox | Numerical simulation of spray coalescence in an Eulerian framework : direct quadrature method of moments and Multi-Fluid method[END_REF] but matters of size velocity correlations are not all solved yet. After integration on the sections, the resulting conservation equations for DNS are similar to those of the pressureless gas dynamics [START_REF] Bouchut | On zero pressure gas dynamics[END_REF][START_REF] Zel | Gravitational instability : an approximate theory for large density perturbations[END_REF] and lead to singular behaviors such as delta shocks and vacuum zones. Well-suited numerical methods are thus required, some specific schemes being presented in [START_REF] De Chaisemartin | Evaluation of eulerian Multi-Fluid versus lagrangian methods for ejection of polydisperse evaporating sprays by vortices[END_REF][START_REF] De Chaisemartin | Eulerian models for turbulent spray combustion with polydispersity and droplet crossing[END_REF][START_REF] De Chaisemartin | Polydisperse evaporating spray turbulent dispersion : Eulerian model and numerical simulation[END_REF]. The model finally requires closure relations for the phenomena accounted for by the Williams equation. We refer to [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][START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] for detailed droplet models for which the Multi-Fluid model can be extended. These extensions do not have any impact on the conclusions of the present study.

The Eulerian Multi-Fluid model has proven its capability to simulate the size-conditioned dynamics of polydisperse sprays including coalescence with the resolution of one moment of the size distribution in each section [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF], the so-called One Size Moment method, even if the number of sections has to be large for accuracy purposes as studied in [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] for coalescing cases. Methods based on two size moments are described in [START_REF] Laurent | Numerical analysis of Eulerian Multi-Fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF], that allow to reduce the number of sections. Their performances have been assessed for evaporating cases but not for coalescence. A two size moment method using exponential size distributions in the sections, introduced in [START_REF] Dufour | A second-order Multi-Fluid model for evaporating sprays[END_REF] for evaporation, is considered in this paper and referred to as the Two Size Moment method. One and Two Size Moment methods are in fact particular cases of a general n moment method described in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]. If the One Size Moment method can solve coalescence [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] and if an attempt has been made to include it in the Two Size Moment method [START_REF] Dufour | Modélisation Multi-Fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF], the purpose of this paper is to extend the Two Size Moment method to include coalescence, to validate the approach in general and to prove its performances in industrial simulations.

In this paper, the Two Size Moment Eulerian Multi-Fluid method is reformulated. It rigorously includes a droplet temperature and its treatment by coalescence. A detailed comparison to the One Size Moment method is thoroughly discussed. We then develop a complete numerical strategy, efficient for implementation and practical use in industrial codes. The Two Size Moment method specific algorithms that are size distribution coefficient inversion and coalescence integral computation are carefully designed and implemented to ensure robustness, accuracy and computational efficiency.

We then conduct a quantitative study of the validity and computational efficiency of the method. For the purpose of this paper, the One Size Moment and Two Size Moment methods have been implemented in a research code developed at EM2C Laboratory solving coalescing dilute sprays in a pseudo 2D nozzle1 with one-way coupling to the gaseous phase, a case which was introduced in [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] to verify the One Size Moment method. This configuration is used to compare the two methods solving the dynamics and size evolution of a lognormal distributed spray. It fully validates the Two Size Moment method in coalescing cases by providing detailed comparisons to Lagrangian and Multi-Fluid reference solutions. Then, an analytical solution in a bimodal case brings more quantitative validation, especially proving enhanced accuracy for the Two Size Moment method. This test-case confirms the excellent robustness of the Two Size Moment method towards steep and monodisperse distributions that are usually difficult to handle in a Eulerian formalism. The research code is then used to simulate the complete dynamics of an experiment on coalescence [START_REF] D'herbigny | Etude expérimentale et numérique pour la validation d'un modèle de coalescence[END_REF] which was conducted at ONERA to verify collision efficiency models. This work furthers the experiment data processing and provides a precious verification for such two-phase flow models, which too often lack experimental back up. These validations and verifications globally show the compliance of the Two Size Moment method to the features required for SRM simulations i.e. accuracy on polydispersity and dynamics, advanced coalescence models, robustness and fast computation.

In addition, we have extended the Two Size Moment Multi-Fluid model implemented in CEDRE, a multi-physics 3D industrial-oriented code developed at ONERA, to include coalescence with the efficient algorithms that are described, tested and validated in this paper. The CEDRE code, which provides fully coupled aero-thermochemical resolution for energetics problems, is used to perform two SRM simulations with the Two Size Moment method to achieve the validation and the proof Nozzle D'Herbigny TEPTEU (SRM) LP10 (SRM) of efficiency of such models in a complex physical background. The first case provides an ultimate validation to a reference Lagrangian simulation. The second case illustrates the ability of the Two Size Moment method to simulate in a reasonable time a complete solid propulsion unsteady test-case. The importance of coalescence in solid propulsion physics is brought to light through its influence on instability levels before concluding the study.

The paper is organized as follows. Section 1 is dedicated to the derivation of the Multi-Fluid method as a theoretical framework. In section 2, we evaluate the accuracy/cost compromise provided by the two Multi-Fluid methods implemented in the research code on a size/velocity coupling validation test-case. In section 3, we validate the Two Size Moment method by comparing it to a purposely derived analytical formula, in a 1D bimodal configuration. The difficulty lies in handling the bimodal size distributions, that are approximated by steep functions in our Eulerian formalism. In section 4, implemented collision and coalescence efficiency models are challenged to a 1D coalescence efficiency experiment that corresponds to the configuration simulated in the previous section. In section 5, we give 2D numerical results in two modeled engines, both composed of a chamber and a nozzle, using the CEDRE code. The different test-cases studied in sections 2, 3, 4 and 5, the corresponding codes and validation resources used are recapitulated in Tab. 1.

Mesoscopic Eulerian spray modeling : two Multi-Fluid methods

In this section, we introduce the framework of our study : the kinetic description of the disperse phase and the derivation of two Eulerian resolution methods. The separation between the two methods appears when presuming the size distributions inside the sections. The origin and assumptions of the corresponding coalescence models are detailed and an integration strategy is given to compute the coalescence terms.

A kinetic description : the Williams equation

Let us define the number density function (NDF) f φ of the spray, where f φ (t, x, u, T, φ)dxdφdT du denotes the average number of droplets (in a statistical sense), at time t, in a volume of size dx around a space location x, with a velocity in a du-neighborhood of u, with an internal temperature T in a dT -neighborhood of T and with a size in a dφ-neighborhood of φ. When modeling the droplet with an infinite conductivity hypothesis [START_REF] Abramzon | Droplet vaporization model for spray combustion calculations[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], we can choose to parametrize the droplet internal state with the only temperature variable. Its state function is then the enthalpy, which only depends on T and reads h(T ) = h ref + T T ref c p,l (θ)dθ, where c p,l (T ) is the droplet heat capacity and the index "ref" denotes a reference state.

When focusing on polydispersity, the size parameter φ of droplets is of primary importance but its natural expression depends on the phenomena : volume v is relevant towards conservation of matter, surface S towards evaporation and radius r towards impact parameters for instance. Since we assume spherical droplets, the equivalence relation v = 4 3

πr 3 = 1 6 √ π S 3
2 allows the size to be expressed in this paper with notations chosen to be most comfortable. As for the size parameter conversions, we shall keep in mind that f r (r)dr = f S (S)dS = f v (v)dv and use the implicit notation f . The evolution of the NDF is described by a Boltzmann-like equation, the Williams transport equation [START_REF] Williams | Spray combustion and atomization[END_REF][START_REF] Williams | Combustion Theory (Combustion Science and Engineering Series)[END_REF]. Considering local momentum and heat transfer with the gas due to drag forces and conduction-convection fluxes and considering collisions between droplets but no evaporation nor break-up, it reads :

∂ t f + u • ∂ x f + ∂ u • (Ff ) + ∂ T H c p,l (T ) f = Γ ( 1
)
where F is the drag force per unit mass, H is the heat transfer per unit mass and Γ is the collision source term. These terms require models which should take into account physical phenomena at the droplet scale. As an illustration, the Stokes law can model F when the particle Reynolds number defined in Eq. (A.1) is moderate. This force per unit mass is due to the velocity difference with the gaseous phase :

F(t, x, u, S) = u g (t, x) -u τ u (S) , τ u (S) = ρ l S 18πµ g
where u g (t, x) is the gas velocity, µ g its dynamic viscosity and ρ l is the supposedly liquid droplet material density. In the same conditions known as Stokes' Regime, the heat transfer H per unit mass reads [START_REF] Sirignano | Fluid Dynamics and Transport of Droplets and Sprays[END_REF] :

H(t, x, T, S) = c p,l (T ) T g (t, x) -T τ T (S) , τ T (S) = 3 2 c p,l (T ) c p,g (t, x) Pr τ u (S)
where T g (t, x) is the local gas temperature, c p,g (t, x) is the gas specific heat capacity and Pr is the Prandtl number, ratio of the gaseous kinematic viscosity and heat diffusivity. Finally, the kinetic modeling for the collision operator is taken from [START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF]. We then assume [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] :

[HC1] We only take binary collisions into account.

[HC2] The mean collision time is very small compared to the inter-collision time.

[HC3] Mass, momentum and droplet enthalpy are preserved during collisions.

[HC4] Every collision leads to coalescence.

Hypotheses [HC1] and [HC2] result from the small liquid phase volume fraction in the context of moderately dense sprays. Hypothesis [HC3] is a classic conservation hypothesis for most types of collisions. We assume that heat converted from relative kinetic energy by inelastic dissipation is neglected in the enthalpy balance. Hypothesis [HC4] is discussed in [START_REF] Brazier-Smith | The interaction falling water drops : coalescence[END_REF][START_REF] Achim | Simulation de collisions, coalescence et rupture de gouttes par une approche lagrangienne: application aux moteurs à propergol solide[END_REF][START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF] since colliding droplets can bounce on each other or separate by reflexion or stretching if the remaining internal kinetic energy of the new droplet is too high. In this context, the CEDRE code provides the Brazier-Smith model (depending on velocities, droplet material viscosity and surface tension, etc.) but it is not used in the following simulations (though adaptations are straightforward) and is therefore not detailed any further. The collision operator resulting from these hypotheses is developed as Γ = Q + -Q -where Q + and Q -respectively correspond to the quadratic integral operators associated with creation and destruction of droplets due to coalescence. Considering two precursor droplets of volumes v ⋆ and v ⋄ colliding to form a new droplet of volume v, the kinetic coalescence operators read [START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF][START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF]:

Q + = 1 2 u ⋆ ,T ⋆ v ⋆ ∈[0,v] f (t, x, u ⋄ , T ⋄ , v ⋄ )f (t, x, u ⋆ , T ⋆ , v ⋆ )B(|u ⋄ -u ⋆ |, v ⋄ , v ⋆ )Jdv ⋆ du ⋆ dT ⋆ Q -= u ⋆ ,T ⋆ v ⋆ f (t, x, u, T, v)f (t, x, u ⋆ , T ⋆ , v ⋆ )B(|u -u ⋆ |, v, v ⋆ )dv ⋆ du ⋆ dT ⋆ (2)
where

v ⋄ = v-v ⋆ , v ⋄ u ⋄ = vu-v ⋆ u ⋆ , v ⋄ h(T ⋄ ) = vh(T )-v ⋆ h(T ⋆ ) are pre-collisional parameters, J is the Jacobian of the mapping (v, T, u) → (v ⋄ , T ⋄ , u ⋄ ) : J = (v/v ⋄ ) n d +1 c p,l /c ⋄ p,l
with n d the dimension of the velocity phase space, c p,l = c p,l (T ) and c ⋄ p,l = c p,l (T ⋄ ) and B(|u -u ⋆ |, v, v ⋆ ) is the collision/coalescence probability kernel which reads :

B(|u -u ⋆ |, v, v ⋆ ) =|u -u ⋆ | β(v, v ⋆ )E(|u -u ⋆ |, v, v ⋆ ).
In this kernel, β(v, v ⋆ ) = π (r + r ⋆ ) 2 is the impact parameter or geometric cross-section and E is the collision efficiency. This efficiency accounts for the effects of the surrounding gas flow. Models can be found in [START_REF] Langmuir | The production of rain by a chain reaction in cumulous clouds at temperatures above freezing[END_REF][START_REF] Beard | Numerical collision efficiencies for small raindrops colliding with micron size particles[END_REF] and are recalled in Appendix A.

At this point, we can make quantitative remarks on the physics of coalescence that will be helpful to design relevant test cases in the following. We link coalescence intensity to the droplet rate of growth since coalescence is hereafter the only phenomenon to modify the droplet size. We define a local droplet growth time τ G (t, x, v) for a given size v :

τ G (t, x, v) -1 = v 0 v ⋆ v u ⋆ ,T ⋆ f (v ⋆ )B(v, v ⋆ )du ⋆ dT ⋆ dv ⋆ (3)
based on the coalescence rate with all the droplets that are smaller than v and weighted by the partner's relative volume v/v ⋆ . This is a local and instantaneous definition but it gives an insight into the intensity of coalescence : after a time τ G (t, x, v) and if all the conditions persist, droplets of size v which have not met any bigger droplet have encountered a volume growth of v i.e. their volume has doubled. The local growth time can be plotted as a chart to analyze the zones where coalescence is intense. Similarly, its evolution in time gives insight into the unsteadiness of coalescence.

In the following, we aim at characterizing the configurations for comparison purposes so a more general growth time is to be defined :

τ G (v) = m p (v) M βγ Ē (4)
where m p (v) is the mass of a particle of size v, M is an average total mass concentration, β is an average cross-section, γ is an average velocity difference, and Ē is an average collision efficiency. All these quantities are estimated to be representative on a volume and a period of interest, their order of magnitude being generally well known and not too variable. Except velocity differences that are notably responsive since they are here the reason for coalescence. We suggest the following estimation γ = ǭg ūg | τ u (v) -τ u | where ǭg and ūg are an average strain rate and velocity of the gas and τ u is an average dynamic time of the disperse phase. Finally the size v used for τ G (v) is chosen as the most represented at the entrance of the volume of interest. These growth times decrease -or equivalently coalescence intensifies-when either collision partner densities, velocity differences or collision efficiencies increase. We can define a coalescence Knudsen number Kn g as the ratio of the global growth time to a typical flow time : coalescence is intense when Kn g is close to one or smaller. This definition is analogous to the classical Knudsen number that measures the non-dimensionalized collision frequency [START_REF] Ferziger | Mathematical Theory of Transport Processes in Gases[END_REF][START_REF] Struchtrup | Macroscopic Transport Equations for Rarefied Gas Flows[END_REF] except that the velocity difference must be resolved in our balistic coalescence case instead of coming from microscopic agitation. The importance of coalescence is here measured by the non-dimensionalized volume growth according to the use of τ G . The characteristic flow time that is used to define Kn g can be a particle residence time τ R in the volume of interest : we then define coalescence intensity as the growth of a particle in the considered volume. In SRM nozzles, the gas accelerates from subsonic to supersonic regime on a very short distance, yielding high velocity differences but small residence times. In the decelerating nozzle that is considered in section 2, the residence time is much higher so that we take small velocity differences to keep a comparable coalescence intensity. The Kn g are estimated in the following to characterize the different configurations regarding coalescence.

Semi-kinetic model

The formalism and the associated assumptions needed to derive the Eulerian Multi-Fluid models are introduced in [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]. We shall now recall the two main steps which are the semi-kinetic derivation and the sectional integration in order to precisely introduce the coalescence terms.

In a first step we reduce the size of the phase space to the only droplet size variable. We only consider moments in the velocity and temperature variables at a given time, a given position and for a given droplet size that are the droplet number density n = f dudT , the average momentum nū = uf dudT and the average enthalpy n h = h(T )f dudT . They indeed depend only on (t, x, S). We can now define a notation for an effective temperature T so that the disperse phase formally uses the same relation to compute its enthalpy moment h = h( T ) than individual droplets2 . In order to close the system, the following assumptions are introduced :

[HV1] For a given size S and location (t, x), the only characteristic velocity is the average ū(t, x, S).

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

[HT1] For a given size S and location (t, x), the only characteristic temperature is the average T (t, x, S).

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

It is equivalent to presume the following form of NDF :

f (t, x, u, T, S) = n(t, x, S)δ(u -ū(t, x, S))δ(T -T (t, x, S)).
The set of hypotheses [HV1] and [HV2] known as the monokinetic hypothesis has been introduced in [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] and is equivalent to reducing the velocity distribution support to a one dimensional sub-manifold parameterized by droplet size. It is correct when τ u (S) is small compared to the gas characteristic dynamic time [START_REF] De Chaisemartin | Evaluation of eulerian Multi-Fluid versus lagrangian methods for ejection of polydisperse evaporating sprays by vortices[END_REF][START_REF] Massot | Eulerian Multi-Fluid models : modeling and numerical methods[END_REF][START_REF] De Chaisemartin | Eulerian models for turbulent spray combustion with polydispersity and droplet crossing[END_REF][START_REF] De Chaisemartin | Polydisperse evaporating spray turbulent dispersion : Eulerian model and numerical simulation[END_REF], the ratio of these two times forming the dynamic Stokes number St. These two hypotheses yield errors in the advection term if the Stokes number is close to one or higher but are correct below this limit as demonstrated extensively by comparisons to Lagrangian references in one-way coupled complex combustion cases [START_REF] Fréret | On the role of preferential segregation in flame dynamics in polydisperse evaporating sprays[END_REF][START_REF] Fréret | On the role of preferential segregation in flame dynamics in polydisperse evaporating sprays[END_REF]. The question of dense sprays of inertial particles remains tough and open. Hypotheses [HV1] and [HV2] can also yield errors in the source terms which depend directly on velocity such as coalescence or drag retroaction -in a twoway coupling context-if the Stokes number is too high. Hypotheses [HT1] and [HT2] are similarly introduced in [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] so that temperatures are locally correlated for droplets of a given size. They are correct when the thermal Stokes number formed, with τ T (S), is smaller than one, as well as the dynamic Stokes number. Otherwise, errors occur in temperature-dependent source terms such as evaporation or heating retroaction in a two-way coupling context. This step leads to a system of conservation equations called the semi-kinetic model, which reads :

   ∂ t n + ∂ x • (nū) = Q + n -Q - n ∂ t (nū) + ∂ x • (nū ⊗ ū) = n F + Q + u -Q - u ∂ t (n h) + ∂ x • (n hū) = n H + Q + h -Q - h (5)
where one gets the average transfer terms n F = Ff dudT and n H = Hf dudT .

In the semi-kinetic framework, the coalescence operator yields the evolution rate of the zeroth and first order moments of the velocity phase space. Coalescence phenomenon has no particular reason to preserve the velocity and temperature distributions as Dirac δ-functions. To preserve the monokinetic and zero-dispersion assumptions, these integrals have been evaluated with a formal dispersion which is constrained to be zero (this projection step is detailed in [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]). The same demonstration can be done for the temperature distribution. Both these projections are yet consistent with the size conditioned relaxations induced by drag and heat transfers. So these terms finally read, when omitting the (t, x) dependency :

Q + n = 1 2 v ⋆ ∈[0,v] n(v ⋄ )n(v ⋆ )β(v ⋄ , v ⋆ )I + n dv ⋆ Q - n = n(v) v ⋆ ∈[0,+∞[ n(v ⋆ )β(v, v ⋆ )I - n dv ⋆ Q + u = 1 2 v ⋆ ∈[0,v] n(v ⋆ )n(v ⋄ )β(v ⋆ , v ⋄ ) v ⋆ ū(v ⋆ ) + v ⋄ ū(v ⋄ ) v ⋆ + v ⋄ I + n dv ⋆ Q - u = n(v)ū(v) v ⋆ ∈[0,+∞[ n(v ⋆ )β(v, v ⋆ )I - n dv ⋆ Q + h = 1 2 v ⋆ ∈[0,v] n(v ⋆ )n(v ⋄ )β(v ⋆ , v ⋄ ) v ⋆ h( T (v ⋆ )) + v ⋄ h( T (v ⋄ )) v ⋆ + v ⋄ I + n dv ⋆ Q - h = n(v)h( T (v)) v ⋆ ∈[0,+∞[ n(v ⋆ )β(v, v ⋆ )I - n dv ⋆ (6)
where I - n and I + n are the partial collisional integrals, computed with the size correlated velocities :

I + n =|ū(v ⋆ ) -ū(v ⋄ )| E(|ū(v ⋆ ) -ū(v ⋄ )|, v ⋆ , v ⋄ ) I - n =|ū(v) -ū(v ⋆ )| E(|ū(v) -ū(v ⋆ )|, v, v ⋆ )
The phase space of the semi-kinetic system has too high a dimension to allow efficient deterministic resolution and requires therefore further modeling. We choose the Multi-Fluid model, a finite volume method on size moments that has proven to be efficient on polydisperse cases.

Multi-Fluid model

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

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

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

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

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

As for [HS1], it allows to reduce the size distribution information in each section at (t, x) to a set of moments of S, the number of which depends on the choice of the (κ k ) k set of size presumed form functions.

Let us now consider two methods based on different forms of presumed functions for [HS1]. In the One Size Moment method, a one parameter function in each section decouples size dependence 1 κ k (S) and space-time dependence m k (t, x), yielding a first order granulometry convergence with the number of sections. This is the classical Multi-Fluid method which has been developed and validated for evaporating cases [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] and in coalescing cases [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. We then introduce the Two Size Moment method using a two-coefficient 2 κ k (t, x, S) function which yields a second order granulometry convergence. The general method has been validated for evaporating cases [START_REF] Dufour | A second-order Multi-Fluid model for evaporating sprays[END_REF][START_REF] Laurent | Numerical analysis of Eulerian Multi-Fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF].

Although the two methods both have a first order reconstruction on size [HV3] and temperature [HT3], the One Size Moment method computes drag and heating with a fixed reconstruction while the Two Size Moment method computes these terms with a size distribution that is better described inside each section. Moreover the quality of the size reconstruction has a strong impact on the velocity and temperature errors since drag and heating times τ u and τ T usually depend on the square of the droplet radius. So the Two Size Moment method indirectly achieves a better approximation of velocity and temperature.

Both methods then solve for size moments calculated on the sections, the equations of which are derived by integrating the semi-kinetic system [START_REF] Beard | Numerical collision efficiencies for small raindrops colliding with micron size particles[END_REF]. Please note that the more parameters the presumed function has, the more moments one needs. The choice of refining the size distribution description with a high order moment method for a finer resolution of polydispersity therefore increases the computational cost. In the context of evaporation, the generalization to more moments can allow to reduce the number of sections to one [START_REF] Kah | Prise en compte des aspects polydispersés dans le contexte d'une approche couplée Eulérienne-Lagrangienne pour la modélisation d'un jet de carburant dans les moteurs à combustion interne d'écoulements diphasiques[END_REF][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]. Regarding coalescence, only the One Size Moment method has been validated [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] so far as explained in the introduction.

One size moment Multi-Fluid method

The One Size Moment Multi-Fluid method assumes a fixed size profile in each section, decoupling size and space dependency. The following notation therefore specifies [HS1] in each section k :

n(t, x, S) = m k (t, x) 1 κ k (S)
where m k is the mass concentration of droplets in the k th section, in such a way that :

S k S k-1 1 κ k (S) ρ l 6 √ π S 3/2 dS = 1.
Such an approach only focuses on one moment of the distribution in the size variable : the moment in terms of mass is chosen because it is conserved by coalescence. Please note that the distribution on the last section is a decreasing exponential with a fixed coefficient. This choice allows the final section to treat the bigger droplets but requires not to host a significant part of the mass for the sake of accuracy. It is well suited for evaporating sprays, the size support of which never extends [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]. This is though a major limitation when coalescence occurs since bigger droplets are generated. Comparisons to results provided by the Two Size Moment method when this final, unbounded section is "active" are discussed in section 2.

The conservation equations for the k th section result from the integration of the mass moment of the semi-kinetic system (5) in each section k and reads :

   ∂ t m k +∂ x •(m k ūk ) = 1 C m+ k -1 C m- k ∂ t (m k ūk )+∂ x •(m k ūk ⊗ ūk ) = m k 1 Fk + 1 C mu+ k -1 C mu- k ∂ t (m k hk )+∂ x •(m k hk ūk ) = m k 1 Hk + 1 C mh+ k -1 C mh- k (7)
where 1 Fk is the "classical" average drag force per unit mass on a section and 1 Hk is the average heat transfer per unit mass on a section :

1 Fk = u g -ūk 1 τ u k , 1 1 τ u k = S k S k-1 1 κ k (S)S 3/2 τ u (S) dS S k S k-1 1 κ k (S)S 3/2 dS , 1 Hk = ck T g -Tk 1 τ T k , 1 1 τ T k = S k S k-1 1 κ k (S)S 3/2 τ T (S) dS S k S k-1
1 κ k (S)S 3/2 dS .

The total system thus counts n d + 2 times more equations to be solved on the space domain than the number of sections with n d the dimension of the physical space. The gas-droplet coupling times 10

0 v ⋆ v ⋄ v i v i+1 v j v j+1
k s t r i p l o w e r l i m i t k s t r i p u p p e r l i m i t 1 τ u k and 1 τ T k are fixed as soon as the size discretization is given. An error on the disperse phase velocity and temperature dynamics occurs, especially if the droplet mass of the size distribution that is approximated in a section is supposed to be concentrated close to a section boundary.

v ⋆ + v ⋄ = v k ✲ v ⋆ + v ⋄ = v k+1 ✲ ◗ ◗ ◗ D ijk
For each section k, the coalescence creation terms result from a double integration : on the whole colliding partner size space at the kinetic level in Eq. ( 2) and on the concerned section at the Multi-Fluid level. Yet the second dependency will not coincide with the section after mapping the natural variables of the two precursor colliding partners i and j. When splitting the integration domain thanks to the continuity of section partitioning, one gets elementary integrals Q ijk , triply indexed with the two precursor section numbers i and j and the destination section number k. The corresponding elementary integration domains are illustrated in Fig. 2 and defined as :

D ijk = r ⋆ , r ⋄ r i-1 < r ⋆ < r i ; r j-1 < r ⋄ < r j ; r 3 k-1 < r ⋆3 + r ⋄3 < r 3 k
In the particular case of the One Size Moment method, considering [HV2] and assuming E = 1 [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF], the coalescence integrals Q ijk take the following form after factorizing the mass moments m i and m j :

1 Q ⋆ ijk = D ijk 1 κ i (r ⋆ ) 1 κ j (r ⋄ )π(r ⋆ + r ⋄ ) 2 4 3 πρ l r ⋆3 dr ⋆ dr ⋄ 1 Q ⋄ ijk = D ijk 1 κ i (r ⋆ ) 1 κ j (r ⋄ )π(r ⋆ + r ⋄ ) 2 4 3 πρ l r ⋄3 dr ⋆ dr ⋄ (8)
As for the disappearance terms, they can also be computed as sums of the elementary creation integrals and must be so to ensure the conservation of matter. After some algebra, the coalescence terms

1 C m+ k , 1 C m- k , 1 C mu+ k , 1 C mu- k , 1 C mh+ k and 1 C mh- k read [45] : 1 C m+ k = k i=1 m i m j i-1 j=1 |ū i -ūj | ( 1 Q ⋆ ijk + 1 Q ⋄ ijk ) 1 C m- k = m k N j=1 m j |ū j -ūk | N i=1 1 Q ⋆ kji 1 C mu+ k = k i=1 m i m j i-1 j=1 |ū i -ūj | (ū i 1 Q ⋆ ijk + ūj 1 Q ⋄ ijk ) 1 C mu- k = m k ūk N j=1 m j |ū j -ūk | N i=1 1 Q ⋆ kji 1 C mh+ k = k i=1 m i m j i-1 j=1 |ū i -ūj | ( hi 1 Q ⋆ ijk + hj 1 Q ⋄ ijk ) 1 C mh- k = m k hk N j=1 m j |ū j -ūk | N i=1 1 Q ⋆ kji (9)
The integrands of the Q ijk integrals depend only on size parameters r,r ⋆ and r ⋄ thanks to [HS1]. This allows the Q ijk integrals to be pre-calculated as soon as the section limits and the ( 1 κ k ) k are given, i.e. once and for all at the beginning of a simulation, the advantages and limits of which are discussed in section 1.6.

Two Size Moment Multi-Fluid method

The Two Size Moment method is based on a two-coefficient exponential approximation of the size distribution in each section [START_REF] Dufour | Modélisation Multi-Fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF][START_REF] Dufour | A second-order Multi-Fluid model for evaporating sprays[END_REF]. This means that [HS1] then reads, for S ∈ [S k-1 , S k [ :

2 κ k (t, x, S) = a k (t, x) exp(-b k (t, x)S) where (a k (t, x), b k (t, x)) k yields two moments per section (n k (t, x), m k (t, x)) k :          S k S k-1 2 κ k (t, x, S)dS = n k (t, x) S k S k-1 2 κ k (t, x, S) ρ l 6 √ π S 3/2 dS = m k (t, x)
The choice of an exponential function ensures the positivity of the distribution function. It also aims at reducing the number of sections and is well suited for evaporation, which requires mass flux information at the section boundary. On the other hand, a problem known as realizability emerges from the fact that the transported quantities are the (n k , m k ) k while all source terms are computed by integration of 2 κ k (t, x, S) on S, therefore requiring (a k , b k ) k . The inversion of the previous system as well as all PDE resolution and integration numerical methods must respect realizability conditions on the sections. Practically, (m k /n k ) k ratios are conditioned by the k th section boundaries since they drive information on average droplet volume. The (a k , b k ) k inversion algorithm was improved to reduce its cost and to increase its reliability. The details are provided in Appendix B. The conservation equations for the k th section now read :

       ∂ t n k +∂ x •(n k ūk ) = 2 C n+ k -2 C n- k ∂ t m k +∂ x •(m k ūk ) = 2 C m+ k -2 C m- k ∂ t (m k ūk )+∂ x •(m k ūk ⊗ ūk ) = m k 2 Fk + 2 C mu+ k -2 C mu- k ∂ t (m k hk )+∂ x •(m k hk ūk ) = m k 2 Hk + 2 C mh+ k -2 C mh- k (10)
where the heat and momentum source terms 2 Fk and 2 Hk are integrated similarly as for the One Size Moment method but with 2 κ k (S). The characteristic times 2 τ u k and 2 τ T k can now vary in space and time, depending on the reconstruction coefficients. The velocity and temperature dynamics are better evaluated, especially if the droplet mass of the size distribution that is approximated in a section is supposed to be concentrated close to a section boundary. The total system counts n d + 3 times more equations than the number of sections with n d the dimension of the physical space.

Since the time and space dependency of the size-distribution functions 2 κ k (t, x, S) is no longer factorizable as was m k in the One Size Moment method, the Q ijk integrals must be computed at each time step in each cell on D ijk . Let us introduce the notation :

Ψ ij (t, x, r ⋆ , r ⋄ , |u ⋄ -u ⋆ |) = 2 κ i (t, x, S(r ⋆ )) 2 κ j (t, x, S(r ⋄ ))π(r ⋆ + r ⋄ ) 2 |ū i -ūj | E(r ⋆ , r ⋄ , |u ⋄ -u ⋆ |)
The coalescence integrals have a different homogeneity than in the One Size Moment method. They now include the number or mass information and read :

2 Q n ijk = D ijk Ψ ij (t, x, r ⋆ , r ⋄ , |u ⋄ -u ⋆ |)dr ⋆ dr ⋄ 2 Q ⋆ ijk = D ijk Ψ ij (t, x, r ⋆ , r ⋄ , |u ⋄ -u ⋆ |) 4 3 πρ l r ⋆3 dr ⋆ dr ⋄ 2 Q ⋄ ijk = D ijk Ψ ij (t, x, r ⋆ , r ⋄ , |u ⋄ -u ⋆ |) 4 3 πρ l r ⋄3 dr ⋆ dr ⋄ (11)
The coalescence source terms

2 C n+ k , 2 C n- k , 2 C m+ k , 2 C m- k , 2 C mu+ k , 2 C mu- k , 2 C mh+ k and 2 C mh- k
are still written as direct sums of the "(i, j, k) integrals" and therefore read :

2 C n+ k = k i=1 i-1 j=1 2 Q n ijk 2 C n- k = N i=1 N j=1 2 Q n kji 2 C m+ k = k i=1 i-1 j=1 ( 2 Q ⋆ ijk + 2 Q ⋄ ijk ) 2 C m- k = N i=1 N j=1 ( 2 Q ⋆ kji + 2 Q ⋄ kji ) 2 C mu+ k = k i=1 i-1 j=1 (ū i 2 Q ⋆ ijk + ūj 2 Q ⋄ ijk ) 2 C mu- k = ūk . 2 C m- k 2 C mh+ k = k i=1 i-1 j=1 ( hi 2 Q ⋆ ijk + hj 2 Q ⋄ ijk ) 2 C mh- k = hk . 2 C m- k (12)
As shown in systems ( 7) and ( 10), drag, heat transfers and coalescence yield spatially located source terms in the Multi-Fluid equations, analogous to chemical source terms for instance. In the Two Size Moment method, these source terms must be evaluated by relevantly accurate integrations on the size distribution inside each section. A key issue for the Two Size Moment method is to compute the coalescence source terms for each section, in each cell and at each time step with a good cost over accuracy compromise.

Numerical strategy for the coalescence integral computation

The One Size Moment expressions given in ( 8) are double integrals on size variables and their integrands are space and time independent when assuming E = 1. The time/space dependency m i (t, x) has indeed been factorized. Pre-calculating the integrals 1 Q ⋆ ijk and 1 Q ⋄ ijk once for all allows to compute the coalescence terms as quadratic combinations of Eq. ( 9). The Two Size Moment expressions 2 Q n ijk , 2 Q ⋆ ijk and 2 Q ⋄ ijk given in [START_REF] Chalons | Beyond pressureless gas dynamics : Quadrature-based velocity moment models[END_REF] are size double integrals that can not be pre-calculated so they are more costly. Since the integration 2D domain D ijk is a non trivial polygon in volume variables (Fig. 2), we suggest two strategies based on point-wise evaluations of the integrand. The first strategy is an equidistributed quadrature using Newton-Cotes formula. The second strategy is an adaptive quadrature based on the moments of exponential kernels.

The 1D n-node Newton-Cotes formula approximates integrals as a sum of n point-wise evaluations of the integrand at n predefined abscissae. The method for 2D integrals, based on an n-node equidistributed quadrature and referred to as (NC n), imposes a cost mainly consisting in n 2 evaluations of the integrand Ψ ij per source term through the formula :

Ψ ij dS ⋆ dS ⋄ ≃ n-1 l=0 w l n-1 m=0 w m Ψ ij (S il , S jm ) ( 13
)
where

S il = S i-1 + l S i -S i-1 n-1
and the w l are tabulated coefficients [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]. This achieves a polynomial interpolation which is a priori not suited for exponential distributions.

An original solution is to use adaptive n-node quadratures (Ad n). They presumably better account for steep distributions i.e. where b is far from zero. Indeed, we aim at handling cases where droplets can be concentrated at an edge of the section so that computing the integrals with quadrature points evaluated at adapted abscissae is natural. We re-write the 2D integrand as Ψ ij = 2 κ i (S ⋆ ) 2 κ j (S ⋄ )Φ(S ⋆ , S ⋄ ) and we assume that the Φ part varies smoothly with the two size variables. The 2D quadrature formula keeps only the smooth part and reads :

Ψ ij dS ⋆ dS ⋄ ≃ n-1 l=0 w il n-1 m=0 w jm Φ(S il , S jm ) (14)
The exponential part is included in the quadrature kernel as follows. The moments of the i exponential kernel M p i =2 κ i (t, x, S)S p dS for p from 0 to 2n -1 are computed analytically since p is an integer. One then determines the abscissae S il and weights w il by inverting a linear system ensuring that the moments of order p can be reconstructed as M p i = n-1 l=0 w il S p il for p from 0 to 2n -1. The same calculations are performed for the S jm and w jm required to integrate on the second dimension. The quadrature is called adaptive since the abscissae and weights depend on 2 κ i and 2 κ j , used as weight functions, which depend on (t, x).

Regarding accuracy, each adaptive quadrature achieves a 2n th order integration so we can use less nodes compared to the equi-distributed quadrature. Regarding computational cost, it mainly consists in determining the eigenvectors of two linear systems of size 2n each to find the abscissae and weights. Yet analytical expressions for the 2 node case abscissae and weights [START_REF] Desjardins | A quadrature-based moment method for dilute fluid-particle flows[END_REF] We now perform a comparison of the suggested quadrature methods on their ability to compute the following integral I(β) :

I(β) = s 2 s 1 S 3
This function represents a typical coalescence partial 1D integral parametrized by β, which accounts for the exponential coefficient that causes integration stiffness. After normalizing I(β) by s 2 s 1 exp(-βS)dS and taking arbitrary numerical values for the integration bounds and the second droplet surface parameter, we show in Fig. 3 results compared to a reference integral computed with a 15-stage Romberg method. This proves the failure of polynomial quadratures for moderate to large values of the exponential coefficient. The adaptive quadratures behave extremely well, always achieving a precision better than 10 -10 with 2 nodes. In the following, we therefore use an adaptive quadrature with 2 nodes (Ad2).

Two Size Moment method validation on an academic nozzle configuration featuring size/velocity coupling

In this section, we validate the Two Size Moment method by comparing the two Multi-Fluid methods on their ability to describe the dynamics of a coalescing cloud in a size-velocity coupled case. Considering the size/velocity coupling (Fig. 1), we need a well-suited test-case, inducing coalescence as well as size-conditioned dynamics and difficult enough to highlight the limitations of the methods. We therefore use the nozzle test-case introduced in [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] to challenge the Two Size Moment method and determine for both methods :

• the convergence rate on mass density prediction regarding size discretization,

• the minimum number of sections to capture the physics of size/velocity coupling,

• the usefulness of the Two Size Moment method in the semi-infinite last section.

A decelerating nozzle for coalescence

But first, let us describe the test-case and the numerical strategy developed in the research code for its resolution. The geometry is a 2D axi-symmetrical diverging conical nozzle (Fig. 4). Numerically speaking, the chosen configuration is stationary, 2D axi-symmetrical in space and 1D in droplet size. It is described in detail in [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. Hence, only its essential characteristics are given here. For the problem to be simpler, an incompressible gas flow is considered. It is taken independent from the particles i.e. the particles are one-way coupled. We neglect the wall viscosity and all other sources of vorticity so the streamlines are straight with well chosen injection conditions. The expressions for the gaseous axial velocity u z and the reduced radial velocity u r /r come from the incompressibility and the self-similarity assumptions and they read for z ≥ z 0 :

u z (z) = z 2 0 u z (z 0 ) z 2 , u r r = u z z = z 2 0 u z (z 0 ) z 3
where z 0 > 0 is the coordinate of the nozzle entrance and u z (z 0 ) is the fixed axial velocity at the entrance. The velocity, which is plotted in Fig 4, decreases since the nozzle diverges. Remarkably, the gas residence time τ R = 1 3

z 3 max z 2
0 uz(z 0 ) can here be computed analytically, depending on the output abscissa z max since the the velocity is analytical and the configuration steady. This time can be used to estimate the particles' residence times though they will be slightly lower since the particle velocities are above u z (z) for all z, as shown in Fig. 4.

The particles are injected with velocities that are aligned to the gas streamlines so their trajectories are also straight lines. This assumption is only valid when no coalescence occurs. However, even in the case of coalescence, it remains valid in the neighborhood of the centerline [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. The droplets slow down because of the deceleration of the gas flow in the conical nozzle, however at a rate depending on their size and inertia, as illustrated in Fig. 4. This will induce velocity differences, as visible on the velocity graph, and coalescence. The injection velocity conditions the velocity differences and therefore the coalescence intensity though the residence time will decrease proportionally. But the deceleration rate at the entrance of the nozzle a(z 0 ) = -2u z (z 0 )/z 0 , which depends on the injection velocity because of the geometry, conditions the strength of the size-velocity coupling so the Kn g globally increases with u z (z 0 ). The droplets are made of alumina, their initial velocity is the one of the gas, their initial temperature, fixed at the one of the gas, is T k = 3600K and does not change along the trajectories. Since alumina is liquid at this temperature, droplets can coalesce. No efficiency law is used in this section (E = 1) since we want to validate the method with comparisons to other numerical methods. The lognormal size distribution is considered at injection since it is often used to characterize alumina particle size distributions in a SRM [START_REF] Salita | Survey of recent Al2O3 droplet size data in Solid Rocket Chambers, Nozzles and Plumes[END_REF]. When parametrized in droplet surface variable, it reads :

LN (S, S LN , σ LN ) = 1 S √ 2π log(σ LN ) exp - 1 2 log(S) -log(S LN ) log(σ LN ) 2 (15)
where S LN is called the geometric average surface and σ LN is called the geometric standard deviation. These parameters must not be mistaken for the usual moments of a distribution, they are of empirical use. We choose a lognormal distribution on the surface variable, without loss of generality. As an illustration, two distributions plotted on the radius are given in Fig. 5 that are used in the following. We now give the details of the numerical strategy that is used. The pseudo 2D steady equations describing this configuration, assuming a Stokes Law, are deduced from Eq. ( 7) for the One Size Moment method and read for section k :

     d z (m k ūk ) + 2m k ūk z = 1 C m+ k -1 C m- k d z (m k ū2 k ) + 2m k ū2 k z = m k u g (z) -ūk 1 τ u k + 1 C mu+ k -1 C mu- k (16)
and they are deduced from Eq. ( 10) for the Two Size Moment method :

             d z (n k ūk ) + 2n k ūk z = 2 C n+ k -2 C n- k d z (m k ūk ) + 2m k ūk z = 2 C m+ k -2 C m- k d z (m k ū2 k ) + 2m k ū2 k z = m k u g (z) -ūk 2 τ u k + 2 C mu+ k -2 C mu- k (17)
These systems are ordinary differential equation (ODE) systems which we integrate with an error smaller than 2.10 -4 . This error is defined with the L ∞ norm of the difference on mass concentration repartition to a reference solution.

Finally, we perform simulations on two different configurations :

• the first configuration studies the convergence and accuracy of the methods on coalescence and spray dynamics in a case that is comparable to the one of a SRM i.e. the granulometry relative shift is moderate. So Kn g is close to one to keep the average growth small compared to the average size;

• the second configuration assumes a strong average size increase to study the additional benefit of using two size moments in a final, unbounded section. The total injected mass, injection velocity and nozzle length are higher, making Kn g smaller and increasing coalescence intensity.

Asymptotic study

We here study the convergence of the One Size Moment and the Two Size Moment method on liquid mass concentration. We take for the nozzle z 0 = 0.05m, z max = 0.25m and u z (z 0 ) = 5m.s -1 so the residence time in the nozzle is τ R = 0.4s : these values generate a strong coupling between coalescence and droplet dynamics for the chosen droplet size, which means that the liquid mass density along the nozzle is completely different whether coalescence is considered or not.

The initial injected mass concentration is taken as m 0 = 1.06kg.m -3 , with lognormal parameters that are set to S LN = 1600µm 2 and σ LN = 1.5. This corresponds to a sharp distribution centered on a radius of 11.3µm as shown in Fig. 5. The growth time is estimated τ G = 0.5s for droplets close to the distribution maximum so that Kn g = 1.3. This corresponds to the coalescence intensity of a typical SRM case as computed in section 5, the velocity differences and residence times compounding each other as discussed in section 1.1. The growth is significant and the size velocity coupling is strong : these severe conditions also make the test-case under consideration a very efficient tool for numerical qualification of the two Eulerian methods.

For the sake of the One Size Moment method accuracy, we define a configuration where coalescence brings negligible mass in the last section so we always take the lower bound of the last section above 50µm.

We give in Fig. 6 the maximum error on mass density along the nozzle depending on the sectional discretization. This proves the convergence rate of both methods to be of order two. The One Size Moment method has a slightly better convergence constant.

From an applicative point of view, it is crucial to consider the performance of both methods for coarse size discretizations, say below 5 sections. With such a small number of sections, the One Size Moment method can yield error above 10%. Yet, the Two Size Moment method achieves remarkably accurate computations with always less than 3% relative error on total mass. This applicative advantage of the Two Size Moment method is studied deeply in the following and is the key issue for practical computations.

Minimizing the number of sections for practical applications

In the context of industrial codes, we want to use as few sections as possible. We therefore extensively compare both methods on coarse discretization cases regarding their accuracy on granulometry, velocity distribution prediction and mass concentration repartition. For the sake of the One Size Moment method accuracy, we still do not want to transfer too much mass in the last section so we keep the conditions described previously.

We choose the tested numbers of sections : a 5 section test-case illustrates what happens with a very coarse discretization. A 13 section and a 25 section test-cases illustrate the convergence of the method. They are compared to a 53 section One Size Moment run which we use here as a reference solution. These six simulations are successively compared to the reference solution regarding four data sets : the output size distribution at first, a velocity distribution in the nozzle at second, the total mass and number density along the nozzle axis at third and mass and number in 5 size intervals (or section groupings given in Tab. 2) along the nozzle at fourth. First, the correct treatment of the size distribution gives an indication on the method's accuracy. We therefore compare the mass concentration distributions at the nozzle output computed with the two methods and with 5, 13 and 25 sections and computed with the 53 section reference test-case in Fig. 7. With 5 sections, the One Size Moment method strongly overestimates the growth rate while the Two Size Moment method underestimates it. With 13 sections, the trend is the same but the error is smaller. Finally for 25 sections, we consider both methods are roughly converged. The growth rate overestimation by the One Size Moment method for coarse size discretization brings about major consequences for the spray dynamics since the gap of average Stokes number between the bounded sections and the final unbounded section. Second, we check in Fig. 8 the velocity distributions at an abscissa where the velocity differences and the effect of past coalescence are high. As expected, the smaller droplets have more decelerated as we can see on the reference response. The velocity distribution is observed to be monotonic, according to the monotonicity of the dynamic response time τ u (S) with size and the injection condition. But its behavior is not a simple law because some droplets, produced by coalescence, have brought their momentum with them to their new size, the velocity of the size then being average through the monokinetic hypotheses [HV1] and [HV2]. We now compare the results with the two approaches and coarse discretizations. In general the error on velocity increases with size, not because of [HV1] and [HV2] which are very correct for the entire Stokes number spectrum that is considered, but because coalescence is underestimated and so is inertia. Then, we observe that the constant velocity reconstruction in each section yields some error, especially in the 5 section case, but we still consider [HV3] as acceptable, the velocity spectrum to be captured in each section being not too wide as required. We finally compare the two approaches together. The velocity behavior is not very well rendered with 5 sections, though the One Size Method gets it right by chance. Then, the velocity distribution is shown to be better rendered with 13 sections by the Two Size Moment method than the One Size Moment method and [HV3] yields little error for both approaches. Finally no error is visible on the graph with the 25 section Two Size Moment method, while the One Size Moment method still has visible error. Third, let us quantify the influence of the treatment of size and velocity polydispersity on the global dynamics of the spray. Fig. 9 shows the evolution of mass and number total concentrations along the centerline of the nozzle computed with both Eulerian methods (5, 13, 25 sections) compared to the reference solutions (53 section One Size Moment method and Lagrangian reference test-case). For the One Size Moment case, number concentrations are computed considering section average droplet volumes given by integration of the ( 1 κ k ) k functions : data are therefore redundant with mass concentration data. We emphasize the fact that good convergence is achieved by the two methods Groups (µm) 5 section case 13 section case 25 section case G 1 =[0, 12.5 with 25 sections as we can see in Fig. 9, bottom, which is a zoom on the end of the nozzle and on a very small interval of the Y-coordinate. The convergence is proven by comparison to a 53 section case performed with the same research code and to a Lagrangian simulation result performed with another EM2C code.

Fourth and to compare precisely the effect of polydispersity on dynamics, let us finally consider the evolution of mass and number along the nozzle for five size intervals. These intervals do match with the sections in the 5 section case but correspond to section groupings in the other cases as illustrated in Tab. 2. The evolution of the mass concentration of these groupings along the nozzle is given in Fig. 10. It is there obvious that the 5 section One Size Moment error on coalescence is severe, especially in the fifth and last grouping G 5 where little mass should be found, as prescribed. The 5 section Two Size Moment method underestimation of size occurs in groups G 3 and G 4 and is more visible here than in Fig. 7. Yet it has moderate consequences on the total mass evolution as shown in Fig. 9 because these sections have Stokes numbers that are close enough.

As a conclusion this study validates the Two Size Moment Multi-Fluid method. It also shows that the Two Size Moment method converges to the reference solution with a second order slope as does the One Size Moment method but reveals to have, in this configuration, much smaller error when using few sections. Since the error on mass never exceeds 3% with the Two Size Moment method, this method yields acceptable dynamical response with as few as 2 sections which is most appreciated for industrial simulations.

Yet the test-case was limited to bounded sections in order to facilitate comparisons to the One Size Moment method. In the next paragraph, we evaluate the benefit of using two size moments in the final, unbounded section.

Two Size Moment method for the unbounded section

We shall now evaluate the impact of using the Two Size Moment method on an unbounded final section. We therefore change the injection conditions to intensify coalescence. The injection velocity is now U (z 0 ) = 500m/s. The origin of the injection is still z 0 = 0.05m but the length of the nozzle is increased and z max = 1.05m. The initial injected mass concentration is increased to m 0 = 2kg.m -3 with S LN = 1600µm 2 and σ LN = 1.3 that is a slightly narrower size distribution given in Fig. 5. In these new conditions, the growth time is now τ G = 0.5ms for droplets close to the distribution maximum and the residence time in the nozzle is τ R = 0.3s so that Kn g = 0.002 which is much smaller than the previous case. Coalescence is therefore more intense, as expected. Moreover, the particle inertia i.e. the Stokes number is higher so the velocity differences are higher and remain strong all along the nozzle. So coalescence is much more intense : the size distribution shifts to an average radius that is ten times bigger than previously as shown in Fig. 11. The lower bound of the final section is set at the same value for all the computations, and below the final average droplet size to assess the final section's particular role. The mass and Sauter Mean Radius (SMR) in the new conditions are given in Fig. 12 and commented in Tab. 3. They allow us to conclude on the versatility of the Two Size Moment method to capture an initially sharp distribution which shifts far from its initial average radius. The SMR failure occurs when too much mass is in the last section : the One Size Moment method requires less than a few percents mass in the last section while the Two Size Moment method allows the SMR to come quite close to the last section lower bound as shown by the SMR plateaus in Fig. 12 The conclusion of this additional study is that a two size moment method implemented in the final unbounded section strongly extends the allowed range for the size distribution which proves to be useful to reduce the number of sections. To complete the bounded section study with practical computational information, Tab. 4 recapitulates the duration of runs using different computation methods. The space integration cost is negligible compared to the coalescence source term computation so these costs uniquely arise from the coalescence integration method and we can use them to conclude on the efficiency of the different quadratures. First, (Ad2) reveals to be the fastest among the Two Size Moment quadrature methods. Moreover the quadrature cost has no reason to depend on the test-case. Since validations performed all along this paper show that (Ad2) brings satisfactory results, this quadrature method is naturally preferred for the Two Size Moment method.

Remarks on computational time

Second, it is obvious that the One Size Moment method is much faster on the nozzle test-case, thanks to coalescence integral pre-calculation. For the same number of sections, there are two orders of magnitude of difference on computational cost. For the same level of accuracy, consider the 5 section Two Size Moment simulation and compare it to the 25 section One Size Moment simulation which yields comparable error on mass repartition : the Two Size Moment method is 10 times slower than the One Size Moment method though having 5 times less sections. No firm conclusion regarding the cost in real applications can however be drawn from these numbers since the computational cost of spatial transport in multi-D configurations would then have to be taken into account. This cost is of course severely reduced when using less sections. The final discussion on cost is thus postponed to section 5.

Suitability of the Two Size Moment method for steep size distributions

We now perform a validation of the Two Size Moment Multi-Fluid method on a 1D configuration called the "D'Herbigny" case. This configuration considers the average radius growth rate of a big droplet falling in a fog of small droplets. The corresponding experiment is described in [START_REF] D'herbigny | Etude expérimentale et numérique pour la validation d'un modèle de coalescence[END_REF] and in section 4 and all the hypotheses are chosen to fit with the experiment conditions. We here consolidate the validation by comparing the computed big droplet distribution to a purposedly derived analytical distribution. We show that the Two Size Moment method achieves higher accuracy than the One Size Moment method.

This configuration provides tough testing for Eulerian models because it deals with bimodal size distributions, which are rather suited for the Lagrangian point of view since they are trivial to sample while their steepness remains troublesome in Eulerian modeling. The results provided here illustrate therefore the robustness of the Two Size Moment method and its ability to remain accurate when dealing with non-smooth size distributions.

Modal analytical resolution

We derive an original analytical approximation of the semi-kinetic system (5) in the particular D'Herbigny case in order to calculate the mean volume and volume standard deviation of the big droplet distribution. The main idea is to take advantage of the fog droplet monodispersity to foretell the big droplet discrete accessible sizes and to project their NDF on these sizes, thus reducing the phase space dimension by discretizing on size "modes".

Since the droplets belonging to the fog have a constant volume (v = v s ) and the bigger droplets can only grow by coalescence (v ≥ v b ), let us first split the semi-kinetic NDF n(t, z, v) regarding the size support as follows :

n s (t, z, v) = n(t, z, v = v s ) n b (t, z, v) = n(t, z, v ≥ v s )
where z is the 1D space coordinate. Of course, n = n s +n b . We now make the fundamental hypothesis :

[He1] n b is much smaller than n s .

This allows us to separate the different scales of the semi-kinetic system [START_REF] Beard | Numerical collision efficiencies for small raindrops colliding with micron size particles[END_REF]. Naturally, the coalescence terms involving (n b , n b ) (coalescence of big droplets) are too small to be taken into account. So the equations can be truncated as follows :

Order 0

∂ t n s + ∂ z • (n s ū) = ✘ ✘ ✘ ✘ ✘ ✘ Q + n (n s , n s ) -✘ ✘ ✘ ✘ ✘ ✘ Q - n (n s , n s ) ∂ t (n s ū) + ∂ z • (n s ū ⊗ ū) = n s F + ✘ ✘ ✘ ✘ ✘ ✘ Q + u (n s , n s ) -✘ ✘ ✘ ✘ ✘ ✘ Q - u (n s , n s ) (18)
Order 1

∂ t n b + ∂ z • (n b ū) = Q + n (n s , n b ) + Q + n (n b , n s ) -Q - n (n s , n b ) -Q - n (n b , n s ) ∂ t (n b ū) + ∂ z • (n b ū ⊗ ū) = n b F + Q + u (n s , n b ) + Q + u (n b , n s ) -Q - u (n s , n b ) -Q - u (n b , n s ) ( 19
) Please note that coalescence terms involving (n s , n s ) vanish since small droplets all have the same velocity, thanks to the monokinetic hypothesis [HV1]. The term F stands here for the sum of the forces on droplets i.e. drag and gravity. Considering negligible forces in the zeroth order momentum equation i.e. neglecting gravity for the fog which sediments too slowly, the small droplet velocity is constant. Choosing the fog as the reference frame, we have ū(v s ) = 0. Finally the zeroth order mass equation shows that small droplets have a constant concentration which makes the problem steady as long as the injection conditions are. We shall therefore consider :

n s (t, z, v) = α s δ(v -v s )
The One Size Moment equations describe the big droplet evolution. The term Q - n (n s , n b ) vanishes on the big droplet size support. Moreover,

Q + n (n b , n s ) = Q + n (n s , n b ) and Q + u (n b , n s ) = Q + u (n s , n b ) so that : ∂ z • (n b ū) = 2.Q + n (n s , n b ) -Q - n (n b , n s ) ∂ z • (n b ū ⊗ ū) = n b F + 2.Q + u (n s , n b ) -Q - u (n b , n s ) ( 20
)
With the previous hypotheses on the fog features, one can easily calculate the remaining coalescence terms, where we have considered E = 1 :

Q + n (n s , n b ) = α s n b (z, v -v s )β(v -v s , v s ) |ū(v -v s )| /2 Q - n (n b , n s ) = α s n b (z, v)β(v, v s ) |ū(v)| Q + u (n s , n b ) = α s n b (z, v -v s )β(v -v s , v s )ū(v -v s ) |ū(v -v s )| (v -v s )/(2v) Q - u (n b , n s ) = α s n b (z, v)β(v, v s )ū(v) |ū(v)|
We now make a continuity assumption on velocities so that, considering the very small size of the fog droplets (v s ≪ v b ), we have :

[He2] ū(v -v s ) ≈ ū(v) for v ≥ v s .
Replacing the momentum equation by a linear combination of the two equations in [START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF], the big droplet evolution system is finally composed of number and velocity equations that are decoupled :

∂ z n b (z, v) = α s [β(v -v s , v s )n b (z, v -v s ) -β(v, v s )n b (z, v)] ✭ ✭ ✭ ✭ ✭ ✭ ✭ ✭ ✭ n b ū(z, v)∂ z ū(z, v) = n b (z, v) F -α s β(v -v s , v s )ū(z, v) |ū(z, v)| n b (z, v -v s ) v s v (21)
with the boundary condition :

n b (z = 0, v) = α b δ(v -v b ) ū(z = 0, v = v b ) = ū0
where α b is the big droplet number density and ū0 their velocity at the boundary. Let us rewrite the decoupled momentum equation in [START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF] to characterize the bigger droplet behavior in the D'Herbigny conditions :

n b (z, v) F = α s β(v -v s , v s )ū(z, v) |ū(z, v)| n b (z, v -v s ) v s v
The rhs models the effect of coalescence with the static fog droplets. Absorbing static mass induces a momentum dilution i.e. a slow down. When neglecting this effect, this equation reduces to a classic dynamic balance ( F = 0) which, in our case, means that the drag force compounds the weight. Therefore, droplets do reach a terminal velocity but this limit velocity increases with the droplet size since the Stokes number depends on the size. To illustrate this fact, Fig. 13 shows the evolution of the big droplet velocity in the D'Herbigny configuration computed with the Two Size Moment method, supposing that droplets are injected with an initial velocity of 3m/s. The black curve is the dynamic equilibrium solution ( F = 0) for the corresponding droplet Sauter mean radius : they both increase. Finally, the slight overestimation yielded by the terminal velocity approach comes from neglecting the so-called momentum dilution. This difference is no longer negligible when coalescence becomes intense. See the case in Fig. 13, right of a C v = α s v s = 60ppm fog where the terminal velocity is 45% overestimated.

C v = 10ppm C v = 60ppm
C v = 10ppm; Right : C v = 60ppm.
Now that we have discussed the evolution of velocity, confirming that the big droplet distribution remains monokinetic and slowly increases with size, let us solve the mass equation by projecting it on size modes. A remarkable point in the bimodal limit condition is that it yields size modes since the volume increment is quantified by the small droplet volume v s . Therefore the big droplet volume v can only take the following values :

v k = v b + kv s ; k ∈ N
This observation allows us to exactly project the mass equation from system [START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF] on the size modes.

Let us rewrite the big droplet size distribution :

n b (z, v) = +∞ k=0 n k (z)δ(v -v k ).
The mass equation becomes a system of ODEs with the boundary conditions n 0 (0

) = α b and n k (0) = 0, k ≥ 1. Defining ǫ = vs v b , b k = [(1 + kǫ) 1 3 + ǫ 1 3 ] 2 3 and z-1 = α s π 3v b 4π 2 3
, it reads :

               d z n 0 (z) = - 1 z b 0 n 0 (z) d z n 1 (z) = 1 z [b 0 n 0 (z) -b 1 n 1 (z)] ... d z n k (z) = 1 z [b k-1 n k-1 (z) -b k n k (z)] (22)
We highlight that the z/b k are lengths that define the typical coalescence lengths or mean coalescing free paths for droplets n k . Consider now another length criterion, that evaluates the importance of volume growth rate for the initial class of big droplets :

L g = v b v s z b 0 .
If z/b 0 is the typical travel length of a big droplet between two coalescing collisions with a fog droplet that is a mean free path regarding coalescence, we need v b /v s times more distance for a big droplet to absorb enough matter to double her own volume. Length L g gives the length after which the big droplet average volume has significantly changed because of coalescence. This length is defined as was the growth time τ G in section 1.1. This means that the non-dimensionalized number that describes coalescence intensity in the D'Herbigny configuration is a coalescence Knudsen number Kn g = L g /z based on a length, for instance the length of the tunnel, in accordance with the steadiness of the problem.

To solve system [START_REF] Dufour | Modélisation Multi-Fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF], we consider two modelings whether z is small or big compared to L g which allows us to approximate β constant or not.

Approached and exact analytical formulae for bimodal coalescence

As a first estimation, we assume that b i = b 0 and the collision lengths are all equal to z 0 . We formalize this temporary hypothesis :

[He3] the collision lengths z/b k are all equal to z 0 = z/b 0 . This is equivalent to assuming the cross-section β(v k , v s ) = β(v b , v s ) = β constant so [He3] is valid only when the big droplet size does not vary too much i.e. z ≪ L g . System [START_REF] Dufour | Modélisation Multi-Fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF] can now be integrated recursively, yielding the following modal densities :

n k (z) = α b 1 k! z z 0 k exp(-z/z 0 ) (23)
In fact, the modal densities follow Poisson's law with z/z 0 being the law parameter. This means that after traveling a z = kz 0 distance in the fog, the most probable droplet size is v k = v b + kv s which corresponds exactly to k coalescence events. We now want to compare formula [START_REF] Dufour | A second-order Multi-Fluid model for evaporating sprays[END_REF] to simulations. In all the following bimodal simulations i.e. in sections 3 and 4 of the paper, we enable the size discretization of the big droplet distribution to be as fine as possible so we choose non equidistributed section bounds : the first section encloses the entire small droplet distribution; a dummy section accounts for the gap; the big droplet distribution lives in all the remaining sections, that are of constant width in radius. The number of sections given to specify the discretization of the following simulations is therefore the one of the big droplets that is sections above r b . A comparison of formula [START_REF] Dufour | A second-order Multi-Fluid model for evaporating sprays[END_REF] with a fine section grid Two Size Moment simulation in the case of few collisions is given in Fig. 14 to illustrate the relevance of the analytical model. In this configuration, [He3] is perfectly valid since ǫz/z 0 = v b /v s z/z 0 = (50) -3 z/z 0 ≪ 1. We also note that Poisson's law can be approximated with a normal distribution with very good match as soon as z/z 0 > 5 which is clear in Fig. 14, bottom. This fact eases computations for high z. Let us add that the modal i.e. discrete nature of the exact distribution is obvious in Fig. 14 and is correctly captured in our Two Size Moment simulations. Yet, the number of modes increases with z/z 0 and the corresponding number of sections required becomes unsustainable. To explore higher z, we now consider simulations that are modally unresolved because of "coarse" section discretizations, the big droplet distribution being then computed as continuous. In this first modeling, with β constant, we can analytically evaluate the moments of the distribution

n b (z, v) = n k (z)δ(v -v k
) by direct summation of the analytical number density of Eq. ( 23). The average volume reads :

v(z) = vn b (z, v)dv n b (z, v)dv = v b + v s z z 0 (24)
and the granulometry standard deviation reads :

σ 2 (z) = (v -v) 2 n b (z, v)dv n b (z, v)dv = v 2 s z z 0 (25)
The model validity domain is limited to small deviations of z/b k from z 0 i.e. v(z) ≈ v b . This corresponds to a coalescence number criterion kv s ≪ v b which retrieves the condition on coalescing length :

z ≪ L g
This first model yields an increasing size dispersion for the big droplets. This dispersion is intrinsic to the way droplets grow, independently of any velocity or cross-section variations. Moreover, the resulting distribution is, as said previously, fairly gaussian as it is a Poisson law with a high parameter. This invalidates any monodisperse approach for the D'Herbigny configuration.

The results of the β constant model are a first step that confirms the validity of our approach. Yet in the D'Herbigny conditions, one has z/z 0 ≈ 2.10 5 when C v = 60ppm while v b vs ≈ 10 5 so that the coalescence number criterion is not fulfilled. This discrepancy is obvious when comparing the analytical and simulation-derived growth rate laws. Indeed the analytical formula [START_REF] Dukowicz | A particle-fluid numerical model for liquid sprays[END_REF] gives a linear increase of average volume with C v , which does not match the linear radius increase obtained numerically and experimentally and showed in Fig. 17. The constant β model is not satisfactory in the D'Herbigny conditions.

As a second and broader modeling, we now solve system [START_REF] Dufour | Modélisation Multi-Fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF] with the exact i.e. variable β(v k , v s ) coefficients. We thus remove hypothesis [He3]. We still have n 0 (z) = α b exp(-s 0 z z ) and we can prove by induction that :

n k (z) = α b k-1 i=0 b k-1 b k -b i λ ik-1 [exp(-b i z z ) -exp(-b k z z )] with λ ik = k-1 j=0 b j k j=0 j =i (b j -b i ) (26)
No further analytical expression has been derived from this formula, which has to be computed directly or at least approximated as done in the following.

Confrontation of simulations to the analytical model

We now consider a particular case where the fog volume concentration goes from 10ppm to 60ppm and the travelled distance is z = 2m. We therefore have z ≫ L g . Let us compare the different model and simulation results. The One Size Moment simulations are performed with 200 sections while using 80 sections for the Two Size Moment ones. The different NDF are given in Fig. 15 for the One Size Moment, Two Size Moment, constant β and variable β models.

First, we note that the predicted average radii are close except with the constant β model. The constant β model fails, as expected, because z ≫ L g invalidates [He3] which is equivalent to saying that the big droplet radius variation is too important and requires to reevaluate the impact parameter. The other three models well agree on the radius growth. They predict an average radius growth ∆r b that is linearly dependent on C v and is confirmed in Fig. 17 in section 4. In contrast, the constant β model predicts a volume growth linearly dependent on fog concentration, which is natural since coalescence intensity is in this case only dependent on C v .

Second, there is a strong difference in Fig. 15 between a refined One Size Moment simulation where the size distribution is extremely widespread, and a coarser Two Size Moment simulation where the size support remains narrow. The peak widening effect featured by the One Size Moment method increases with fog concentration as visible in Fig. 15. This widening effect increases also with z. In general, steep distributions are abnormally widespread along their way in the size phase space : this is comparable to numerical diffusion in the size phase space. We refer to this effect as numerical diffusion though it results from numerical errors on source terms and not on fluxes. It is remarkable that these test-cases are performed with a high number of sections so the One Size Moment numerical diffusion is unacceptable since it will be even more penalizing in practical cases where size discretizations are coarse. In contrast, the Two Size Moment method brings moderate diffusion in size phase space as shown in Fig. 15 since widening is smaller to that due to physical size dispersion quantified in our analytical models. The Two Size Moment method handles correctly the D'Herbigny case regarding average size and size dispersion. As a conclusion, the Two Size Moment Multi-Fluid method, here proven to be robust and validated, is as accurate and far less diffusive than the One Size Moment method when dealing with coalescing steep size distributions. This feature is relevant for SRM applications where the alumina cloud is sometimes modeled as bimodal [START_REF] Price | Combustion of Metallized Propellants[END_REF].

C v = 10ppm C v = 30ppm C v =

An experimental validation showing the importance of collision efficiency modeling

In this section, we validate the Two Size Moment Multi-Fluid method by referring to a 1D coalescence experiment which was conducted at ONERA by F.-X. D'Herbigny [START_REF] D'herbigny | Etude expérimentale et numérique pour la validation d'un modèle de coalescence[END_REF]. This experiment, referred to as the D'Herbigny experiment can be simulated with the configuration described in section 3. The importance of modeling collision local phenomena is here highlighted and the Two Size Moment method proves to capture the physics of such a configuration.

The D'Herbigny collision efficiency experiment

The D'Herbigny experiment consists in the growth of a big droplet falling through a fog of smaller droplets. The main features are summed up in Fig. 16. Details about the experimental device can be found in [START_REF] D'herbigny | Etude expérimentale et numérique pour la validation d'un modèle de coalescence[END_REF]. We only specify that the average injection radius of the bigger droplets is r b = 150µm and their initial velocity u s = 3m/s thanks to a piezoelectric injector. The fog droplet radius r s ∈ [2µm, 4µm] is much smaller. These conditions bring the efficiency law parameters Re p and k, defined in Eq. (A.1), close to typical SRM ones. The growth time for the big droplets is τ G = 1.1s when assuming E = 1 and C v = 60ppm for the small droplet concentration. The residence time in the tunnel is τ R = 1.7s for the big droplets so that Kn g = 0.6.

We perform simulations with the Two Size Moment Multi-Fluid method in the experiment conditions, though testing six fog droplet concentration values, with C v from 10ppm to 60ppm. We use the three collision efficiency laws described in Appendix A i.e. E = 1, the Langmuir-Blodgett law from Eq. (A.2) and the Beard-Grover law from Eq. (A.3). The measurements in the D'Herbigny experiment and the computations that we performed of the average big droplet radius after 5m in the fog, depending on its concentration, are given in Fig. 17. The experiment and the simulations show a linear dependency of radius growth on fog concentration. This allows us to determine an average collision efficiency coefficient Ē given in Tab. 5 for the two empirical laws and for the D'Herbigny experiment. The simulation results finally show that collision efficiency laws have a dramatic effect on radius growth, undermining it up to 72% in the Langmuir-Blodgett case for instance.

As for the experiment results, they corroborate the linear dependency of radius growth on fog concentration that was proven in section 3. They also confirm the dramatic undermining of coalescence by the collision efficiency factor, tending to designate the Langmuir-Blodgett law even-though the

❄ ❄ ❄ g ✛ Transparent tunnel (5m) ✛ Monodisperse injector ❜ ❜ ❜ ❝ ❝ ✟ ✟ ✟ ✙ Granulometry measurement ❍ ❍ • Laminarization Grid
Fog generator (r s , C v , u s ) experimental average collision efficiency is even smaller than the model forecast. We therefore conclude that coalescence simulations in SRM conditions should take a collision efficiency law into account. As a conclusion, the D'Herbigny experiment underlines the importance of collision efficiency model. Moreover, though the experimental results are shadowed by uncertainty on collision efficiencies, local fog concentrations and size measurements, their comparison to simulations still provides an ultimate validation of the Two Size Moment method, regarding the linear dependence of growth on fog concentration.

Conclusion on the model and methods

The previous results in sections 2, 3 and in the present section confirm that Eulerian Multi-Fluid methods can be used to simulate accurately the size distribution evolution of a coalescing spray and its size-conditioned dynamics. When collision efficiency is taken equal to one, the One Size Moment method provides, with a reasonable number of sections -at least 10-extremely fast results since the coalescence term computations are reduced to quadratic combinations of section masses. The Two Size Moment Multi-Fluid method provides good results with a very coarse size-space discretization (as few as 5 sections) and is far less "diffusive" in size phase space for steep size distributions.

But section 4 proves the importance of modeling the collision efficiency factor. The Two Size Moment method features any such modeling for very little additional cost. A general form for coalescence efficiency formulae to be suitable for a pre-calculated One Size Moment method are suggested in [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] but no classical collision efficiency models can be easily implemented since they require local gas parameters. The One Size Moment method can therefore no longer benefit from pre-calculation when including such models.

We finally highlight that the quantitative conclusions that are drawn regarding the performance of the Multi-Fluid section discretization on all these cases totally depends on the respective configurations. First because the Stokes number local spectra, which condition the treatment of velocity and temperature polydispersity in the sections, varies from a case to another. Second because the coalescence intensity, which conditions the size phase space refinement requirement, varies from a case to another. Since the nozzle test case has levels of polydispersity and coalescence that are comparable to the ones of a SRM, the previous conclusions on the number of sections and on coalescence intensity quantitatively hold for the motor simulations of the following section. Likewise the D'Herbigny experiment has comparable levels of coalescence parameters so the conclusions on coalescence efficiency impact hold. Of course, the qualitative conclusions on the advantages of the Two Size Moment method hold in any case.

Validation and feasibility on SRM configurations

In this section, we conclude the study by achieving two 2D simulations with an industrial-oriented code in which we have implemented the Two Size Moment Multi-Fluid method featuring efficiency models :

1. a simplified SRM steady configuration allows us to compare our Eulerian Multi-Fluid coalescence method to a Lagrangian reference method 2. a realistic geometry of a SRM sub-scale model allows us to challenge the new method on a complex unsteady case by capturing acoustic chamber modes and vortex shedding.

Both simulations are two-way coupled. We recall that sections 2, 3 and 4 have provided detailed comparisons of the Two Size Moment method with Lagrangian and experimental data on configurations that are representative of the peculiarities of many two-phase flows, especially in a SRM. Therefore, the role of the following section is to assess the feasibility of the approach in real configurations, especially its numerical cost, and highlighting the reasons why our approach is efficient. So the Langmuir-Blodgett collision efficiency law is enabled for feasibility demonstration purposes, in accordance with the satisfactory results of section 4. Finally these simulations illustrate the impact of coalescence in solid propulsion.

The CEDRE code

The CEDRE code is a multi-physics platform on general unstructured grids, for both research and industrial applications, 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 follows a multi-domain, multi-solver approach. Solvers are considered for each physical system : gas phase, dispersed phase, thermal fields in solids and radiation. These solvers share the CEDRE architecture and libraries, and can be coupled to perform a multi-physics computation or be operated alone. They tackle subjects as varied as multiphase flows, multispecies chemistry, thermal conduction, radiation, wall-film models, advanced thermodynamics, etc. The coupling can be one-way or two-way [START_REF] Errera | Multiphysics coupling approaches for aerospace numerical simulations[END_REF]. As for the disperse phase, CEDRE includes a Lagrangian and two Eulerian methods i.e. a multi-class or sampling method and the Two Size Moment Multi-Fluid method studied and extended here [START_REF] Murrone | Numerical modeling of dispersed two-phase flows[END_REF]. Only the Lagrangian and the Multi-Fluid methods provide coalescence resolution, the latter being coded and validated in this paper.

A simplified steady validation : the TEPTEU case

The first SRM simulation is performed on a 2D simplified configuration, yet featuring the main difficulties of SRM typical flows i.e. parietal injection and a supersonic nozzle. Parietal injection refers here to two-phase mixture injected normally to the wall : it is a classical modeling for the gas and inclusions produced and ejected by the solid propellant combustion. A flow that is injected from a cylindrical wall may have a peculiar structure with for instance Taylor-Culick instabilities, widely studied in simple cases.

The configuration has two injection zones as shown in Fig. 18 : a parietal injection corresponding to the flow coming from the wall that carries small alumina droplets, typical of immediately burnt solid propellant and an axial injection corresponding to the flow coming from the rear, that carries bigger droplets to model priorly injected, burnt and coalesced particles. The upper part of the chamber generating these bigger droplets is not solved so that the mesh is significantly shortened. The meeting of these two types of particles may generate coalescence.

The simulation takes place on a deformed-structured 1500 cell mesh (Fig. 18). When seeking a two-phase stationary flow solution, one usually takes a converged gas flow field as an initial state. With a total flow rate of 10kg.s -1 .m -2 from the wall (propellant combustion) and the head end (upstream flow), this gas flow involves extremely high velocity gradients in the nozzle (Fig. 19, top), which will induce size-conditioned droplet dynamics as studied in section 2. In the chamber, the growth time is about τ G = 24ms and the residence time τ R = 10ms so that Kn g = 2.4. In the nozzle, these times are τ G = 0.5ms and τ R = 0.2ms so that Kn g = 2.5. The two-phase simulation strictly speaking starts when injecting the disperse phase : a monomodal 5µm radius wall injection represents the droplets resulting from recently burnt aluminum particles, directly expelled from the propellant grain, and another roughly monomodal distribution around 20µm is injected on the axis to model the previously coalesced droplets [START_REF] Fabignon | Combustion of aluminum particles in solid rocket motors[END_REF]. For the purpose of our simulation, the injected volume fractions approximate correctly typical SRM chamber conditions and preserve the dilute spray hypothesis. For the disperse phase, we choose a 5 section discretization (0, 12.5, 21, 25 and 30 µm) as in section 2 so that only the first or the second section contains droplets at the wall. A 10 -2 s time interval is required to allow the first droplets to reach the end of the nozzle. With a 10 -6 s time step (10, 000 iterations), we perform a 1h single processor AIX platform computation which corresponds to 360ms per cycle (about 50% more than without coalescence). The gaseous Mach number field, once the permanent two-phase regime is reached, is given in Fig. 19, bottom; the stationary volume fractions for the five sections are displayed in Fig. 20. The first section droplets do come from the wall, they fill the combustion chamber and the nozzle, where their volume fraction decreases due to the strong acceleration, which we refer to as dilution. The second section droplets do come from the rear and they are pushed towards the axis by the parietal injection. The three other section droplets are not present close to the wall since they are not injected.

wall injection (5µm) ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ 20µm ✲ ✲ ✲ ✲ ✲
First, we observe in Fig. 19 that the disperse phase has an impact on the gas flow, in accordance with the two-way coupling enabled in the simulation. Specific impulse loss can therefore be observed with the nozzle Mach number decrease in the presence of droplets. Second, we note in Fig. 20 the creation of bigger particles : coalescence yields droplets in section three as soon as the two injected types of droplets meet; moreover the volume fraction of section three droplets in the chamber is significant which shows the importance of coalescence. Third, even bigger droplets, especially in section four, are created in the nozzle even though acceleration induces there a strong volume fraction dilution. The Lagrangian reference computation, also performed with the CEDRE code, takes 6h on the same platform. Fig. 21 shows a mean diameter (d 30 ) with both simulations : its evolution is correctly predicted by the Eulerian Multi-Fluid method as implemented in CEDRE. The main difference occurs on the axis where the Eulerian model undermines coalescence. Indeed momentum is averaged on each section because of the monokinetic hypothesis [HV3] and thus the Eulerian model is unable to account for particle trajectory crossing (PTC) among the sections. This effect is also visible close to the wall where streamlines become tangential. On the contrary, the Lagrangian framework naturally handles PTC and can feature any wall rebound model, which allows the particle streamlines to avoid the wall.

Another difference is visible in the nozzle : there are no particles far from the axis in the Lagrangian simulation (see the first section volume fraction in Fig. 20, dark blue convention in Fig. 21). On the contrary, particles are present almost everywhere in the nozzle with the Eulerian simulation. This reveals that the Eulerian transport scheme is spatially diffusive. This also appears on the size gradient at the rear of the chamber, which is smoothened. Yet spatial diffusion is rather a matter of numerical method and can be easily improved [START_REF] De Chaisemartin | Polydisperse evaporating spray turbulent dispersion : Eulerian model and numerical simulation[END_REF] but this is not the point of this paper. Fig. 22 shows the relative error on the Mach Number field compared to the reference Lagrangian simulation. Mach number is important to quantify the total specific impulse loss in the nozzle. The defects of the Eulerian Multi-Fluid method, previously commented on the d 30 evolution, appear here too. Size and velocity prediction errors on the axis imply an overestimation of the Mach number at the rear end. In the nozzle, particle trajectories are axi-parallel and close to the axis, still because of Eulerian momentum averaging, so that the Mach number is overestimated up to 10% in the near-axis zone. Finally, small particle Eulerian diffusion increases flow inertia close to the wall so that the Mach number is slightly underestimated (5%).

To conclude this 2D complex validation, we observe on the one hand that computational time is dramatically cut down with the Eulerian method, which is six time faster and on the other hand, that monokinetic hypothesis and numerical diffusion generate small discrepancies of the Eulerian momentum and volume fraction, having up to 10% repercussions on the Mach number in our current configuration.

Complex unsteady configuration : the LP10 case

The second SRM simulation is performed on a 2D representation of a real sub-scale motor called LP10. To be more precise, the simulation takes place at a time when one has observed the most intense instabilities. The simulated time is long enough to perform a harmonic study of the instabilities. At that time, the star shaped grain has totally burnt so that no flow comes from the front end. Droplets are only injected from the two segment walls (Fig. 23) with the same lognormal distribution : coalescence occurs only because of the droplet size dispersion3 with an intensity that is comparable to the one of the nozzle test case.

The configuration still features parietal injection and a supersonic nozzle but has a more realistic geometry and we aim at capturing flow instability. The test-case is therefore more complex and numerically costly. The simulation now takes place on a deformed-structured 27000 cell mesh showed in Fig. 23. The case is fully solved on a 300ms time interval with a 1µs time step in order to capture unsteady pressure signals and perform modal frequency analysis. The disperse phase granulometry is solved with a 3 section Two Size Moment method where the section lower bounds are : 3.33, 6.66 and 10µm. The 300, 000 iterations are performed with a 32 CPU NEHALEM platform in 10 hours. This corresponds to 800ms per cycle (25% more than without coalescence). Regarding coalescence intensity in the chamber, the growth time is about τ G = 24ms though the velocity difference can be higher in vortices compared to the average value we take. The residence time in the chamber is τ R = 20ms so that Kn g = 1.2. In the nozzle, these times are close to those of the TEPTEU case with τ G = 0.5ms and τ R = 0.2ms so that Kn g = 2.5. We can evaluate the growth of particles from section three : they have a d 30 time average diameter of 14.1µm when crossing the nozzle throat on the centerline. Since they were all injected with a 12.83µm diameter, they encounter an increase of 1.3µm due to coalescence in the chamber. Results not presented here show that the diameter increase reaches 1.8µm when the collision efficiency is set to one. As expected, not using collision efficiency modeling results in an overestimation of coalescence intensity. So it is demonstrated here on a practical case that advanced collision efficiency modeling plays a sensitive role. To illustrate the configuration unsteadiness, let us consider the gaseous vorticity field in Fig. 24 which is typical of Parietal Vortex Shedding (PVS) instabilities [START_REF] Lupoglazoff | Parietal vortex shedding as a cause of intstability for long solid propellant motors. numerical simulations and comparisons with firing tests[END_REF]. These instabilities occur in long, segmented motors and are strongly coupled to the disperse phase. The instantaneous volume fractions for the three sections and the third section mean diameter (d 30 ) are displayed in Fig. 25. As expected, the three droplet sizes react differently to the gaseous vortices described in Fig. 24. The third section is mostly segregated, since its characteristic time is close to the vortices' one ( St = 1). One also notices that the aft-dome cavity filling depends strongly on droplet size. These phenomena, along with specific impulse loss conditioned by size (illustrated in Fig. 22), are among the main issues motivating polydisperse two-phase SRM studies as exposed in our introduction. Therefore, the importance of accurately capturing polydispersity becomes obvious. Pressure signal spectrum is displayed as a log-log diagram in Fig. 26. The first harmonic has a level 20% lower than in a polydisperse non-coalescing simulation. This gives quantitative evidence that coalescence modeling is a crucial issue in SRM simulation to capture instability levels. A numerical study on instability levels depending on polydispersity and coalescence can be found in [START_REF] Doisneau | Two-way coupled simulation of acoustic waves in polydispersed coalescing two-phase flows : application to Solid Rocket Motor instabilities[END_REF].

wall injection (lognormal) ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ wall injection (lognormal) ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄

Computational cost of coalescence with the Two Size Moment method of CEDRE

The numerical features of the 2D SRM simulations and the computational times with and without coalescence are gathered in Tab. 6 and compared to the Lagrangian computational time. Accounting for coalescence is nonetheless proven to be necessary, but we show here that it is numerically accessible in the Eulerian Multi-Fluid framework, yielding at most a doubling of the numerical cost for 5 sections with the Two Size Moment method. The dramatic increase of the cost of the Two Size Moment method with the number of sections on a complex, practical case is visible here. So why have we not used the One Size Moment method since we have proven in section 2 that the One Size Moment method is much cheaper on the nozzle test-case (Tab. 4)? First, it is not possible to meet the high computational speed of the One Size Moment method on the nozzle test-case if we abandon pre-calculation to include specific collision/coalescence laws, which are indeed required (section 3). Second, the Two Size Moment coalescence integral computations, which are local, are significantly accelerated on parallel architectures with many-core CPUs so these architectures are well adapted to the Two Size Moment method. Third, for a given level of accuracy, the One Size Moment method requires more sections than the Two Size Moment method. For typical SRM simulations, each section requires the resolution of an unsteady 3D pressureless Euler-type system. So that increasing the number of sections is a costly operation, thus limiting the profit yielded by the One Size Moment method. This fact can be seen on the computational times of Tab. 6 where solving for coalescence represents half of the total cost while it represented almost the whole cost for the nozzle case (Tab. 4).

Conclusion on the feasibility of SRM simulations

The necessity to use two size moments to capture coalescence with Eulerian methods is proven through the accuracy and feasibility studies of this paper. The ability of the Two Size Moment Eulerian Multi-Fluid method as implemented in CEDRE to simulate unsteady polydisperse sprays encountering coalescence in a complex geometry at a reasonable computational cost is proven by the final simulations. It is the only method that is able to capture size/velocity coupling with such a coarse size discretization, this being discussed for two-phase acoustics in [START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF][START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF]. Higher order size moment methods with advanced velocity reconstruction are promising regarding more inertial droplet treatment or computational cost-cutting via the reduction of the number of sections [START_REF] Vié | Size-velocity correlations in high order moment methods for polydisperse evaporating sprays : modelling and numerical issues[END_REF] but their compatibility with coalescence modeling has not been addressed up to now.

Our approach is suitable and ready for tackling the first DNS issues in solid propellant combustion among four considered in the introduction : hydrodynamic and acoustic instabilities. The quantitative study of this problematic on the LP10 sub-scale motor is fully discussed in [START_REF] Doisneau | Two-way coupled simulation of acoustic waves in polydispersed coalescing two-phase flows : application to Solid Rocket Motor instabilities[END_REF]. For other studies involving smaller droplets (typically micronic), the modeling is compatible with two-way coupling but the numerical strategy is to be improved to prevent the use of very small time steps : this question is tackled in [START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF][START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF] where a highly tunable splitting strategy adapted to the stiffness of moderately dense polydisperse spray Eulerian systems is described and quantitatively evaluated. This coupling approach is used in the case of small droplets with two-way coupling and coalescence in [START_REF] Doisneau | Eulerian VS Lagrangian simulation of unsteady two-way coupled coalescing two-phase flows in solid propellant combustion[END_REF][START_REF] Doisneau | Eulerian VS Lagrangian simulation of unsteady two-way coupled coalescing two-phase flows in solid propellant combustion[END_REF]. Finally, the approach scales towards practical DNS configurations as expected and proven with the computational cost study of Tab. 6. It is moreover a required basis before including further modeling such as LES for instance.

The extension of spray modeling to LES configurations requires to remove hypothesis [HV2], which is incompatible with sub-grid scale effects [START_REF] Boileau | Robust numerical schemes for Eulerian spray DNS and LES in two-phase turbulent flows[END_REF][START_REF] Chalons | Beyond pressureless gas dynamics : Quadrature-based velocity moment models[END_REF]. Other velocity closures (e.g. gaussian) yield turbulent transport terms [START_REF] Simoes | Modélisation eulérienne de la phase dispersée dans les moteurs à propergol solide, avec prise en compte de la pression particulaire[END_REF][START_REF] Masi | Étude théorique et numérique de la modélisation instationnaire des écoulements turbulents anisothermes gaz-particules par une approche Euler-Euler[END_REF], turbulent collision kernels [START_REF] Fede | Numerical study of the subgrid fluid turbulence effects on the statistics of heavy colliding particles[END_REF][START_REF] Wunsch | Theoretical and numerical study of collision and coalescence -Statistical modeling approaches in gasdroplet turbulent flows[END_REF] and turbulent coalescence kernels [START_REF] Zaichik | Statistical models for predicting the effect of bidisperse particle collisions on particle velocities and stresses in homogeneous anisotropic turbulent flows[END_REF]. The context of sub-micronic droplets also brings sub-grid scale closures into play. Classical modeling [START_REF] Friedlander | Smoke, Dust and Haze, Fundamentals of Aerosol Dynamics[END_REF] has to be adapted to account for Brownian agglomeration and out-of-equilibrium transport.

For most droplet size ranges in a SRM, drag is strong enough to locally correlate velocities to size, allowing to consider our method as valid, especially hypothesis [HV1]. In specific cases however, the Stokes number is high enough for droplets not to track the flow [START_REF] De Chaisemartin | Evaluation of eulerian Multi-Fluid versus lagrangian methods for ejection of polydisperse evaporating sprays by vortices[END_REF][START_REF] Massot | Eulerian Multi-Fluid models : modeling and numerical methods[END_REF][START_REF] De Chaisemartin | Eulerian models for turbulent spray combustion with polydispersity and droplet crossing[END_REF][START_REF] De Chaisemartin | Polydisperse evaporating spray turbulent dispersion : Eulerian model and numerical simulation[END_REF] and PTC can occur. Hypothesis [HV1] impedes to account for same size droplet PTC since it locally averages velocities for a given size. For instance droplet crossings should occur on the centerline of our simulation but we have instead a droplet accumulation artifact due to the momentum averaging on the symmetry axis. High-inertia droplets ejected from vortices and from the solid propellant surface need a specific multi-velocity treatment such as quadrature methods [START_REF] Desjardins | A quadrature-based moment method for dilute fluid-particle flows[END_REF], and high order velocity moment methods [START_REF] Chalons | A multi-gaussian quadrature method of moments for gasparticle flows in a les framework[END_REF][START_REF] Kah | Eulerian quadrature-based moment models for dilute polydisperse evaporating sprays[END_REF]. Multi-velocity approaches are yet rarely extended to coalescence except in [START_REF] Fréret | Turbulent combustion of polydisperse evaporating sprays with droplet crossing : Eulerian modeling of collisions at finite Knudsen and validation[END_REF].

Other issues in SRM modeling are close to being accessible to Eulerian Multi-Fluid simulation in CEDRE. First, to refine the specific impulse loss quantification in the nozzle, secondary break up modeling has to be added. Indeed high velocity gradients in the nozzle induce droplet secondary break-up which strongly undermines droplet growth. Models that are compatible with the Eulerian Multi-Fluid approach [START_REF] Dufour | Modélisation Multi-Fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF] have recently been implemented in the CEDRE code. Second, aluminum combustion, moving meshes, bi-species droplets and evaporation modeling are still to be implemented in CEDRE to tackle combustion instabilities in burst regime. They are also compatible with the present Eulerian formalism.

Conclusion

We provide a development and a comprehensive validation of a two size moment Eulerian Multi-Fluid method based on exponential size reconstructions for solving two-phase polydisperse coalescing flows. The conclusion of this study is that the new method is able to simulate accurately the dynamics of such flows with very few sections. The compromise on capturing the size-velocity correlations is acceptable in solid rocket motor configurations and reducing the number of required sections is extremely profitable to reduce the computational cost of complex simulations. The Two Size Moment method is the only Eulerian Multi-Fluid method that can include validated advanced collision efficiency models, which are crucial for solid rocket motor studies. Finally, the method is computationally accessible and practically efficient in the industrial-oriented code CEDRE, to predict for instance the level of hydrodynamic and acoustic instabilities in a solid rocket motor, which is one of the main issues in this field.

Appendix A. Collision efficiency models

Collision efficiency E has been introduced in the coalescence terms of Eq. ( 2). It is a probability factor modeling the correlation of droplet velocities immediately before a collision. It is equivalent to consider that droplets may dodge each other due to the gas flow surrounding them. Collision efficiency laws thus require knowledge of local gas parameters such as density ρ g or viscosity µ g . The simplest approach is to neglect this effect and to consider E = 1. This simplifies the computation but a small error on growth rate is amplified when coalescence is intense and droplet size growth significant. In the case of unbalanced droplet sizes, two collision efficiency models are suggested on empirical basis [START_REF] D'herbigny | Etude expérimentale et numérique pour la validation d'un modèle de coalescence[END_REF] : the Langmuir-Blodgett model, and the Beard-Grover model [START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF][START_REF] Achim | Simulation de collisions, coalescence et rupture de gouttes par une approche lagrangienne: application aux moteurs à propergol solide[END_REF]. These laws strongly depend on the bigger droplet particle Reynolds number Re p and on a non-dimensionalized number k which read, when taking r s the smaller radius and r b the bigger radius : The number k is the ratio of the smaller droplet relaxation time τ u (r s ) to its residence time in the bigger droplet influence zone τ R = r b |u b -us| . In the case of low particle Reynolds numbers, Langmuir and Blodgett [START_REF] Langmuir | The production of rain by a chain reaction in cumulous clouds at temperatures above freezing[END_REF] numerically get the following expression :

E 1 (k) = 0, if k ≤ 1.214 E 1 (k) = 1 + 3 ln(2k) 2(k-1.214) -2
, otherwise whereas for high particle Reynolds numbers, they get : Beard and Grover [START_REF] Beard | Numerical collision efficiencies for small raindrops colliding with micron size particles[END_REF] suggest to increase the accuracy of formula (A.2) which results from a simple interpolation between two limits. For this purpose, they use a numerical solution of the incompressible Navier-Stokes equations in order to determine the gaseous flow surrounding the bigger droplet depending on the Reynolds number. They can therefore evaluate precisely the forces on the smaller droplet and compute its trajectory. For Re p ∈ [0, 400] and r b < r s , it yields:

E 2 (k) = 0, if k ≤ 0.0833 E 2 (k) =

Appendix B. Reconstruction algorithm for two size moments

The two size moment Eulerian Multi-Fluid method transports two size moments for section k but the properties that have to be averaged on section [S No analytical expression can be derived for its reciprocal function so that g k is inverted with Ridders' method. This requires to compute g k several times, usually around ten times for the required accuracy. The g k function is explicitly computable with an expression featuring the error function or Dawson function. We consider this method as too costful so we use the adaptive quadrature introduced in 1.6 to compute the 3/2 order moment integration. The accuracy of the adaptive quadrature in this case is presented in Fig. B.28 and is better than 10 -4 for (Ad2).

Figure 1 :

 1 Figure 1: Size/velocity coupling mechanism for a liquid spray.

Figure 2 :

 2 Figure 2: Domains D ijk on which coalescence terms are integrated.

Figure 3 :

 3 Figure 3: Comparison of quadrature errors (left : absolute; right : relative) on the test integral I(β) with S 1 = 1, S 2 = 2 and S ⋆ = 0.75 -: NC5; • : NC9; ∆ : Ad2; ∇ : Ad3; Black : Reference.

Figure 4 :

 4 Figure 4: Geometry of a conical 2D self-similar diverging nozzle (Left) and typical velocities in the nozzle (Right) -(Solid) : Gas u z (z); Dashed : Particles of radius r 1 = 15µm corresponding to a Stokes number St 1 ; Dotted : Particles with r 2 = 30µm and St 2 > St 1 .

Figure 5 :

 5 Figure 5: Lognormal droplet mass distributions with S LN = 1600µm 2 used as boundary condition in the nozzle test case -Solid : m 0 = 1.06kg.m -3 and σ LN = 1.5; Dashed : m 0 = 2kg.m -3 and σ LN = 1.3.

  e r o n e o r d e r t h r e e

Figure 6 :

 6 Figure6: Maximum error on total mass concentration compared to a reference case (One Size Moment 1, 000 sections) -Empty : One Size Moment method; Solid : Two Size Moment method with (Ad2) quadrature.

Figure 7 :

 7 Figure 7: Mass concentration distribution at the nozzle's end (z = 0.25m) computed with the One Size Moment (empty symbols) and Two Size Moment method (solid symbols) with 5, 13 and 25 sections and reference (solid line).

Figure 8 :

 8 Figure 8: Average velocity distribution at z = 0.08m computed with the One Size Moment (empty symbols) and Two Size Moment method (solid symbols) with 5, 13 and 25 sections and reference (solid line).

Figure 9 :

 9 Figure 9: Top : Mass and number concentrations along nozzle -Bottom : zoom on the nozzle's end.

Figure 10 :

 10 Figure 10: Evolution of mass concentration in the four bounded section groupings and in the last section for 5, 13 and 25 sections.

Figure 11 :

 11 Figure 11: Mass concentration distribution at different nozzle abscissae -Solid : 13 sections, One Size Moment; Dashed : 7 sections, Two Size Moment; Histogram : reference).

Figure 12 :

 12 Figure 12: Mass concentration and Sauter Mean Radius along the nozzle.

Figure 13 :

 13 Figure 13: Velocity and Sauter Mean Radius along the tunnel computed with the Two Size Moment method with E = 1 (∆) and terminal velocity (black) corresponding to the computed Sauter Mean Radius -Left : C v = 10ppm; Right : C v = 60ppm.

4 kg/m 3 )Figure 14 :

 4314 Figure 14: Poisson's law (+) compared to Normal law (line) and 200 section Eulerian Multi-Fluid mass density distribution (histogram) and total mass per mode (∆) for z/z 0 = 3.8 (top) and z/z 0 = 19.2 (bottom). For analytical validation purposes, this case specifically features a strictly monodisperse fog with r s = 10µm.

Figure 15 :

 15 Figure 15: NDF after 2m in a 10, 30 and 60ppm fog -Histogram : 200 sections, One Size Moment MF method; ∆ : 80 sections, Two Size Moment MF method; : constant β model; • : variable β model.

Figure 16 :

 16 Figure 16: D'Herbigny experimental device (ONERA).

Figure 17 :

 17 Figure 17: Radius growth ∆r b after 5m depending on the fog concentration C v -Solid : Two Size Moment simulation with E = 1; : Two Size Moment simulation with Langmuir-Blodgett model; • : Two Size Moment simulation with Beard-Grover model; Errorbars : Experimental. Two Size Moment Multi-Fluid simulation Experiment Langmuir-Blodgett Beard-Grover Average efficiency ĒLB = 0.28 ĒBG = 0.46 ĒXP = 0.21

Figure 18 :

 18 Figure 18: TEPTEU : a deformed-structured 1500 cell mesh (Arrows : injection zones).

Figure 19 :

 19 Figure 19: Gaseous Mach number without droplets (top) and with a disperse phase computed with the Two Size Moment method (bottom).

Figure 20 :

 20 Figure 20: Dispersed phase volume fractions per section.

Figure 22 :

 22 Figure 22: Eulerian relative difference to the Lagrangian simulation on the Mach number field (%).

Figure 23 :

 23 Figure 23: LP10 : Deformed-structured 27000 cell mesh (Arrows : injection zones).

Figure 24 :

 24 Figure 24: Gaseous phase vorticity module (s -1 ).

Figure 25 :

 25 Figure 25: Top to bottom : Volume fraction of section 1, 2 and 3 and d 30 mean diameter (µm) of section 3.

Figure 26 :

 26 Figure 26: FFT of pressure sensor signals (log-log diagram).

= 2ρ l r 2 s

 2 Re p = 2ρgr b µg | u g -u b | , k |u b -us| 9µg r b . (A.1)

k 2 (

 2 k+0.5) 2 , otherwise.For intermediate cases they assume the following interpolation :E LB (k, Re p ) = E 1 (k)1 + Re p /60 + (Re p /60)E 2 (k) 1 + Re p /60 . difference |ū b -ūs | (m/s) Langmuir-Blodgett ❝ Beard-Grover

Figure A. 27 :

 27 Figure A.27: Two collision laws selected for SRM applications with r b = 150µm (Solid : r s = 2µm, Dashed : r s = 3µm; Dotted : r s = 4µm)

E

  BG (k, Re p ) = 4 π 2 [arctan(max(H(k, Re p ), 0))] 2 (A.3)whereH = 0.1465 + 1.302Z -0.607Z 2 + 0.293Z 3 Z = ln(k/k 0 ) k 0 =exp(-0.1007 -0.358 ln(Re p ) + 0.0261[ln(Re p )] 2 ). Typical collision efficiency values with these two laws are given in Fig. A.27.

k 3 2

 3 ] are computed assuming a two coefficient exponential reconstruction for (n k , m k ). This requires to compute the reconstruction coefficients as follows :n k b k exp(-b k S k-1 )-exp(-b k S k ) if b k = 0 n k S k -S k-1 if b k = 0 (B.1) with a unique solution because g k (b k ) = m k /n k is a bijective function from R to ρ l 6 depending only on b k : g k (b k ) = ρ l 6 √ π S k S k-1 S exp(-b k S)dS S k S k-1 exp(-b k S)dS

Figure B. 28 :

 28 Figure B.28: Left : The g k inverse function with its limiters -Right : Relative error on g k (b) -∆ : Ad2; ∇ : Ad3; Reference : analytic To avoid overflow, a constant slope limiter is used in the exponential part of the NDF reconstructions that imposes | b k ∆S |< K lim . Its effect on the accuracy of size distribution reconstructions has been evaluated as very moderate. The shape of the reciprocal function once limited is presented in Fig. B.28.

Table 1 :

 1 Summary of the cases studied in this paper and corresponding result media.

	Time dependence	Steady	Steady	Steady	Unsteady
	Physical Space	Pseudo 2D	1D	2D	2D
	Phase space (size)	1D	1D	1D	1D
	One Size Moment MF

  allow to significantly reduce this cost.

	-100 -100 -100 -100 -100 -100 -100	-50 -50 -50 -50 -50 -50 -50	0 0 0 0 0 0 0	50 50 50 50 50 50 50	100 100 100 100 100 100 100	-100 -100 -100 -100 -100 -100 -100 -100	-50 -50 -50 -50 -50 -50 -50 -50	0 0 0 0 0 0 0 0	50 50 50 50 50 50 50 50	100 100 100 100 100 100 100 100

Table 2 :

 2 Composition of the five section groupings

				]			1		1 to 3			1 to 6
			G 2 =[12.5, 25]			2		4 to 6			7 to 12
			G 3 =[25, 37.5]			3		7 to 9			13 to 18
			G 4 =[37.5, 50]			4		10 to 12			19 to 24
			G 5 (> 50µm)			5		13			25
	0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5					0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8					0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
						0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7
	0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
						0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6					0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15
	0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3					0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
						0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4					0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
	0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2					0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
	0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1	G 1				0.1 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1 0.2 0.1 0.2	G 2				0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05			G 3
	0 0 0 0 0 0 0 0 0 0 0 0 0 0					0 0 0 0 0 0 0 0 0 0 0 0 0					0 0 0 0 0 0 0 0 0 0 0 0 0
	0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05	0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1	0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15	0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2	0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25	0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05	0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1	0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15	0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2	0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25	0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05	0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1	0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15	0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2	0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25
	0.1 0.12 0.1 0.12 0.1 0.12 0.1 0.12 0.1 0.12 0.1 0.12 0.1 0.12 0.1 0.12 0.1 0.12 0.1 0.12 0.1 0.12 0.1 0.12 0.1 0.12	G 4				0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25	Last section				Empty 1 size moment
	0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08					0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2						Solid 2 size moment
	0.04 0.06 0.04 0.06 0.04 0.06 0.04 0.06 0.04 0.06 0.04 0.06 0.04 0.06 0.04 0.06 0.04 0.06 0.04 0.06 0.04 0.06 0.04 0.06 0.04 0.06					0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1						∆ ❝ 25 sections 5 sections 13 sections
	0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02					0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05						-Reference
	0 0 0 0 0 0 0 0 0 0 0 0 0					0 0 0 0 0 0 0 0 0 0 0 0 0
	0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05	0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1	0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15	0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2	0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25	0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05	0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1	0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15	0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2	0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25

  .

	Legend	Size moment(s) Sections r(n sec ) Mass failure absc. SMR failure absc.
	Solid -	1	101	200µm	reference	reference
	Empty ∆	1	13	100µm	0.20m	0.20m
	Solid	2	7	80µm	0.50m	0.70m
	Solid ∆	2	13	80µm	0.60m	0.70m
	Dashed --	2	23	80µm	0.80m	0.70m

Table 3 :

 3 Nozzle configurations and result comments.

Table 4 :

 4 Computational time on a 2.66GHz Intel Core 2 Duo CPU.

Table 5 :

 5 Average collision efficiency values in the D'Herbigny configuration

  Overall average droplet diameter (d 30 ) and particle streamlines in Eulerian and Lagrangian computations.

	EULERIAN
	D 1 0 _ f i n a l
	20.0
	17.8
	15.5
	13.2
	11.0
	8.8
	6.5
	4.2
	2.0
	LAGRANGIAN
	Figure 21:

Table 6 :

 6 Computational cost per iteration (time per cycle × number of CPUs).

The configuration can almost be considered as 3D since it involves two space dimensions and one droplet size dimension which have to be resolved.

But T = n -1 R T f dudT a priori

( √ S ⋆ + √ S) 2 exp(-βS)dS.

Droplet size dispersion can be induced by the propellant combustion for instance

Acknowledgments

The present research was done thanks to a Ph. D. Grant from DGA, Ministry of Defense (M. S. Amiet, Technical Monitor) and an ONERA research project (PEA Nano2). The authors also thank G. Dufour and T. Fontfreyde for their preliminary contribution to the Two Size Moment method in the CEDRE code.