R O Fox
email: rofox@iastate.edu

F Laurent

M Massot
email: marc.massot@em2c.ecp.fr

Numerical simulation of spray coalescence in an Eulerian framework: direct quadrature method of moments and multi-fluid method

Keywords: Liquid Sprays, Coalescence, Direct Quadrature Method of Moments, Multi-fluid Model, Spray Equation, Number Density Function

 to the Williams spray equation. Both the multi-fluid method and DQMOM yield systems of Eulerian conservation equations with complicated interaction terms representing coalescence. In order to focus on the difficulties associated with treating size-dependent coalescence and to avoid numerical uncertainty issues associated with two-way coupling, only one-way coupling between the droplets and a given gas velocity field is considered. In order to validate and compare these approaches, the chosen configuration is a self-similar 2D axisymmetrical decelerating nozzle with sprays having various size distributions, ranging from smooth ones up to Dirac delta functions. The second key feature of the paper is a thorough comparison of the two approaches for various test-cases to a reference solution obtained through a classical stochastic Lagrangian solver. Both Eulerian models prove to describe adequately spray coalescence and yield a very interesting alternative to the Lagrangian solver. The third key point of the study is a detailed description of the limitations associated with each method, thus giving criteria for their use as well as for their respective efficiency.

Introduction

In many industrial combustion applications such as Diesel engines, fuel is stocked in condensed form and burned as a dispersed liquid phase carried by a gaseous flow. Two phase effects as well as the polydisperse character of the droplet size distribution (since the droplets dynamics depend on their inertia and are conditioned by size) can significantly influence flame structure. Size distribution effects are also encountered in a crucial way in solid propellant rocket boosters, where the cloud of alumina particles experiences coalescence and become polydisperse in size, thus determining their global dynamical behavior [START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF][START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF]. The coupling of dynamics, conditioned on particle size, with coalescence or aggregation as well as with evaporation can also be found in the study of fluidized beds [START_REF] Tsuji | Cluster patterns in circulating fluidized beds predicted by numerical simulation (discrete particle versus two-fluid models[END_REF] and planet formation in solar nebulae [START_REF] Bracco | Particle aggregation in a turbulent keplerian flow[END_REF][START_REF] Chavanis | Trapping of dust by coherent vortices in the solar nebula[END_REF]. Consequently, it is important to have reliable models and numerical methods in order to be able to describe precisely the physics of two-phase flows where the dispersed phase is constituted of a cloud of particles of various sizes that can evaporate, coalesce or aggregate and also have their own inertia and size-conditioned dynamics. Since our main area of interest is combustion, we will work with sprays throughout this paper, keeping in mind the broad application fields related to this study.

Generally speaking, two approaches for treating liquid sprays corresponding to two levels of description can be distinguished. The first, associated with a full direct numerical simulation of the process, provides a model for the dynamics of the interface between the gas and liquid, as well as the exchanges of heat and mass between the two phases using various techniques such as the Volume Of Fluids (VOF) or Level Set methods [START_REF] Aulisa | A geometrical area-preserving volumeof-fluid method[END_REF][START_REF] Herrmann | A Eulerian level set/vortex sheet method for two-phase interface dynamics[END_REF][START_REF] Josserand | Droplet impact on a dry surface: triggering the splash with a small obstacle[END_REF][START_REF] Tanguy | Development of a level set method for interface tracking: application to droplet collisions[END_REF]. This "microscopic" point of view is very rich in information on the detailed properties at a more local level concerning, for example, the resulting drag exerted on one droplet depending on its surroundings. The second approach, based on a more global point of view, describes the droplets as a cloud of point particles for which the exchanges of mass, momentum and heat are described globally, using eventually correlations, and the details of the interface behavior, angular momentum of droplets, detailed internal temperature distribution inside the droplet, etc., are not predicted. Instead, a finite set of global properties such as mass, momentum, temperature are modeled. Because it is the only one for which numerical simulations at the scale of a combustion chamber or in a free jet can be conducted, this "mesoscopic" point of view will be adopted in the present paper. Furthermore, we are interested in sprays where droplet interactions (e.g., coalescence) have to be taken into account, which corresponds to liquid volume fractions between 0.1% and 1%. O'Rourke [START_REF] O'rourke | Collective drop effects on vaporizing liquid sprays[END_REF] classified the various regimes from the "very thin spray", which are transported by the gaseous carrier phase without influencing the gaseous phase, through the "thin spray" regime, for which there is two-way coupling through the momentum equation between the two phases, up to the "thick spray" regime for which the volume fraction of liquid is high enough so that droplet-droplet interactions have to be taken into account, but is still low enough so that the liquid volume fraction is negligible as compared to the gaseous one. Because our primary focus is on the ability of Eulerian methods to capture droplet coalescence, our study is limited here to the "thick spray" regime. By restricting our attention to one-way coupling, we can avoid difficulties (e.g., grid convergence) associated with using Lagrangian methods with two-way coupling, and it will thus be possible to make detailed comparisons between Eulerian and Lagrangian simulation results.

In the mesoscopic framework, there exists considerable interest in the development of numerical methods for simulating sprays [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][START_REF] Miller | Direct numerical simulation of a confined three-dimensional gas mixing layer with one evaporating hydrocarbon-droplet-laden stream[END_REF][START_REF] Miller | Direct numerical simulation and subgrid analysis of a transitional droplet laden mixing layer[END_REF][START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF][START_REF] Réveillon | Eulerian analysis of the dispersion of evaporating polydispersed sprays in a statistically stationary turbulent flow[END_REF]. The principal physical processes that must be accounted for are [START_REF] Alexander | The direct simulation Monte Carlo method[END_REF] transport in real space, (2) droplet evaporation, (3) acceleration of droplets due to drag, and (4) coalescence of droplets leading to polydispersity. The major challenge in numerical simulations is to account for the strong coupling between these processes. Williams [START_REF] Williams | Spray combustion and atomization[END_REF] proposed a relatively simple transport equation based on kinetic theory that has proven to be a useful starting point for testing novel numerical methods for treating coalescing liquid sprays. In the context of one-way coupling, the Lagrangian Monte-Carlo approach [START_REF] Dukowicz | A particle-fluid numerical model for liquid sprays[END_REF], called Direct Simulation Monte-Carlo method (DSMC) by Bird [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF], is generally considered to be more accurate than Eulerian methods for solving Williams equation. However, its computational cost is high, especially in unsteady configurations. Moreover, in applications with two-way coupling, Lagrangian methods are difficult to couple accurately with Eulerian descriptions of the gas phase. There is thus considerable impetus to develop Eulerian methods for describing sprays. In this paper, we limit our attention to one-way coupling with a given (laminar) gas velocity field (i.e., one-way coupling with a given gas velocity field.) Thus no turbulence models are required to close the spray equation.

In a recent paper Laurent et al. [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] have demonstrated the capability of an Eulerian multifluid model to capture the physics of polydisperse evaporating sprays with one-way coupling. This approach 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 that yields conservation equations for mass, momentum (and eventually other properties such as temperature) of droplets in fixed size intervals called "sections" extending the original work of Tambour, Greenberg and collaborators [START_REF] Greenberg | An opposed jet quasi-monodisperse spray diffusion flame[END_REF][START_REF] Greenberg | On the origin of spray sectional conservation equations[END_REF]. Even though this approach has recently been extended to higher order by Laurent [START_REF] Laurent | Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF] and Dufour [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], the necessity to discretize the size phase space can be a stumbling block in some applications. Moment methods, on the other hand, do not encounter this limitation.

In this work, we apply the recently developed direct quadrature method of moments (DQ-MOM) [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF] to treat Williams equation in a Eulerian framework. As its name implies, DQ-MOM is a moment method that closes the non-linear terms (e.g., droplet coalescence) using weighted quadrature points (abscissas) in phase space. Such a closure relates to the construction of an approximated number density function from a set of moments under the form of a sum of Dirac delta functions, the support of which corresponds to the abscissas. However, it is important to make a clear difference between such an Eulerian approach and the corresponding Lagrangian approach, for which the number density is approximated by a large number of numerical "parcels". The evolution of abscissas and the corresponding weights are governed by the dynamics of a few moments, whereas the evolution of the parcels are governed by the Williams equation since they are a stochastic discretization of this equation. Consequently, the DQMOM usually involves a very restricted number of unknowns on a Eulerian mesh, whereas the Lagrangian method involves a very large number of unknowns that are followed along their trajectories in phase space.

The DQMOM method distinguishes itself from other quadrature methods (e.g., QMOM [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF][START_REF] Marchisio | Quadrature method of moments for aggregationbreakage processes[END_REF]) by solving transport equations for the weights and abscissas directly (instead of transport equations for the moments). The source terms for the transport equations depend on the physical processes involved. For Williams equation, we show in Section 2 that laminar transport and drag result in source terms that are independent of the choice of moments and, in fact, are equivalent to those used in Lagrangian formulations. When evaporation does not lead to the disappearance of droplets in finite time, this is also true for the evaporation process. On the other hand, coalescence leads to a linear system for the source terms for which the coefficient matrix depends on the choice of moments. The applicability of DQMOM to Williams equation thus depends on whether or not a particular choice(s) of moments can be found that leads to a non-singular linear system. When the evaporation law allows the disappearance of droplets in finite time the equations for the moments of the number density function not only involve unclosed integral terms, but also the flux of disappearing droplets, i.e. the pointwise value of the number density function at zero size. This quantity has then to be closed since it has a strong influence on the dynamics of the whole set of moments; it leads to a significant difficulty since it corresponds to the reconstruction of a pointwise value of the number density function from a set of its moments. In this study, we propose a solution to this difficult issue. Note that because spatial transport is treated explicitly, it suffices to tackle the flux problem in the homogeneous case. We will see that a key point is to provide a flux closure that yields stable moment dynamics and a non-singular linear system in the DQMOM framework.

Let us also underline that the transport terms in the systems of conservation equations for both Eulerian models are the same and given by pressureless gas dynamics. The structure of these transport terms and the associated difficulties have been the subject of several studies and there are numerical methods designed in order to treat the resulting singularities as shown in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. The question of the computational efficiency of such Eulerian approaches (especially in coalescing systems) is a key question since these methods are intended to be used in more realistic unsteady configurations as an alternative to the too costly Lagrangian methods for polydisperse sprays. We have already studied this question in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] where the Eulerian multi-fluid approach was shown to offer a good precision with a relative low cost [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. Because of the similarity of the transport terms for both Eulerian approaches, the conclusions about the computational efficiency presented in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] are also valid for the DQMOM method. Consequently we focus our study and comparisons on stationary configurations for which we are sure to have a reference solution at our disposal and from which we can obtain firm conclusions about the capabilities of the various approaches.

In Section 3, we present the chosen test configuration, which is a self-similar 2D axisymmetrical decelerating nozzle and sprays with two inlet distributions: a smooth monomodal function and Dirac delta functions. We also discuss in detail the reasons (e.g., significant coalescence rates) for the choice of the test cases, and why they are particularly challenging for the various numerical methods. Finally the Lagrangian solver, the numerical subtleties for obtaining the associated reference solution, as well as the multi-fluid method are then presented. In Section 4, we consider the results for the various test cases including combinations of coalescence, linear evaporation in terms of volume (since it conserves the number of droplets and thereby eliminates the need to model the evaporative flux) and the usual nonlinear evaporation law (for which the evaporative flux must be modeled.) We present results for the most difficult test cases, designed to highlight the challenges one would encounter in more realistic cases. The results are compared to a reference solution obtained through a Lagrangian stochastic algorithm [START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF]. The advantages and limitations of the Eulerian methods are then analyzed in detail in terms of precision and efficiency. It is shown that the DQMOM method offers very interesting features in a number of situations (e.g., strongly coalescing droplets), and is a good candidate for more complex configurations.

DQMOM for Williams equation

The Williams transport equation [START_REF] Williams | Spray combustion and atomization[END_REF] for the joint volume, velocity number density function f (v, u; x, t) is

∂ t f + u • ∂ x f + ∂ v (R v f ) + ∂ u • (Ff ) = Γ, (1)
where R v is the evaporation rate, F is the drag force acting on the droplet, and Γ is the coalescence term. Note that specific forms for the evaporation rate and drag law are not required for DQMOM. However, in this work we will consider one-way coupling with a given gas velocity that appears in F. Using standard assumptions [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF], we can write the coalescence term in two parts:

Γ = Q - coll + Q + coll where Q - coll = - ∞ 0 B(|u -u * |, v, v * )f (v, u)f (v * , u * ) dv * du * , (2)
Q + coll = 1 2 v 0 B(|u ⋄ -u * |, v ⋄ , v * )f (v ⋄ , u ⋄ )f (v * , u * )J dv * du * , (3)
v ⋄ = v -v * , u ⋄ = (vu -v * u * )/(v -v * ), and J = (v/v ⋄ ) 3 is the Jacobian of the transform (v, u) → (v ⋄ , u ⋄ ) with fixed (v * , u * ).
The collision frequency function B is defined by

B(|u -u * |, v, v * ) = E coal (|u -u * |, v, v * )β(v, v * )|u -u * |, (4)
where E coal is the coalescence efficiency probability, which, based upon the size of droplets and the relative velocity, discriminates between rebound and coalescence, and

β(v, v * ) = π 3v 4π 1/3 + 3v * 4π 1/3 2 . ( 5
)
For simplicity, we will take E coal = 0 (no coalescence) or E coal = 1; however, any other functional form could be used in the derivation that follows. A more general version of the spray equation would include the droplet temperature and molecular composition. For simplicity, we consider only the volume and velocity in this work. Finally note that adding spatial diffusion terms in Eq. ( 1) would generate additional terms in DQMOM [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF].

One of the principal mathematical difficulties when developing Eulerian solvers for Eq. ( 1) is the accurate treatment of the coalescence term. Indeed, the integral form of Γ leads to highly non-local and non-linear interactions in volume-velocity phase space. A "direct" Eulerian solver would require discretization of the high-dimensional phase space (in addition to real space), and would thus be computationally intractable. In contrast, multi-fluid models discretize only the volume phase space and use the average velocity conditioned on droplet size (i.e., the mono-kinetic assumption [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]), while moment methods (such as DQMOM) provide closures based on a finite set of moments. Before applying DQMOM to Eq. ( 1), we should note that the coalescence term is defined such that the moments representing mass and momentum are conserved:

ρvΓ(v, u) dv du = 0 (6)
and

ρvuΓ(v, u) dv du = 0, (7)
where the liquid density ρ is assumed to be constant. These conservation properties must be retained in numerical approximations used to treat Eq. ( 1) (as we shall see is the case with DQMOM).

The DQMOM approximates the density function by weighted delta functions in volumevelocity phase space [START_REF] Fox | Computational Models for Turbulent Reacting Flows[END_REF][START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF]:

f (v, u) = N n=1 w n δ(v -v n )δ(u -u n ) (8)
where

δ(u -u n ) ≡ δ(u 1 -u 1,n )δ(u 2 -u 2,n )δ(u 3 -u 3,n ).
Note that in this formulation, the weights w n and abscissas (v n , u n ) are Eulerian fields. Application of DQMOM results in closed transport equations for the number density, mass density, and momentum density, respectively, of the form:

∂ t w n + ∂ x • (w n u n ) = a n , (9)
∂ t (w n ρv n ) + ∂ x • (w n ρv n u n ) = ρb n , (10)
and

∂ t (w n ρv n u n ) + ∂ x • (w n ρv n u n u n ) = ρc n , (11)
where a n , b n , and c n are source terms that are found from the right-hand side of Eq. ( 1) as described below. These equations can be solved with appropriate initial and boundary conditions to find the fields w n (x, t) and (v n (x, t), u n (x, t)) appearing in Eq. [START_REF] Dufour | A second-order multi-fluid model for evaporating sprays[END_REF]. Note that Eqs. [START_REF] Dukowicz | A particle-fluid numerical model for liquid sprays[END_REF][START_REF] Dupays | Some aspects of twophase flows in solid-propellant rocket motors[END_REF][START_REF] Fox | Computational Models for Turbulent Reacting Flows[END_REF] are equivalent to an Eulerian multi-fluid model [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF], but with the source terms on the right-hand side determined using DQMOM.

The DQMOM approximation for the moments of the number density function are found directly from Eq. ( 8):

v k u l 1 u m 2 u p 3 ≡ v k u l 1 u m 2 u p 3 f (v, u) dv du = N n=1 w n v k n u l 1,n u m 2,n u p 3,n . (12)
The fundamental idea behind DQMOM is that we should choose the weights and abscissas such that as many moments as possible are determined by the moment transport equations found from Eq. (1). Note that there are a total of N weights, N volume abscissas, and 3N velocity abscissas and (equivalently) 5N unknown source terms in Eqs. [START_REF] Dukowicz | A particle-fluid numerical model for liquid sprays[END_REF][START_REF] Dupays | Some aspects of twophase flows in solid-propellant rocket motors[END_REF][START_REF] Fox | Computational Models for Turbulent Reacting Flows[END_REF]. We will thus need to choose 5N independent moments to determine the source terms. We will return to the subject of how to choose the moments in Section 2.4. The procedure for using these moments to find the source terms is described next.

Space and time derivatives

The space and time derivatives in Eq. ( 1) generate the corresponding terms in Eqs. [START_REF] Dukowicz | A particle-fluid numerical model for liquid sprays[END_REF][START_REF] Dupays | Some aspects of twophase flows in solid-propellant rocket motors[END_REF][START_REF] Fox | Computational Models for Turbulent Reacting Flows[END_REF]. These are found by formally inserting Eq. ( 8), and differentiating:

∂ t f + u • ∂ x f = N n=1 δ(v -v n )δ(u -u n ) [∂ t w n + ∂ x • (u n w n )] - N n=1 w n δ (1) (v -v n )δ(u -u n ) [∂ t v n + u n • ∂ x v n ] - N n=1 w n δ(v -v n )δ (1) (u -u n ) • [∂ t u n + u n • ∂ x u n ] (13)
where δ (1) (ψ) = dδ(ψ)/dψ and δ (1) (ψ) is a vector with components δ

(1)

1 (ψ) = δ (1) (ψ 1 )δ(ψ 2 )δ(ψ 3 ), δ (1)
2 (ψ) = δ(ψ 1 )δ (1) (ψ 2 )δ(ψ 3 ), and δ

(1) 3 (ψ) = δ(ψ 1 )δ(ψ 2 )δ (1) (ψ 3 ). Using the definitions of the source terms, Eq. ( 13) can be rewritten as

∂ t f + u • ∂ x f = N n=1 δ(v -v n )δ(u -u n ) + v n δ (1) (v -v n )δ(u -u n ) a n - N n=1 δ (1) (v -v n )δ(u -u n ) -v -1 n δ(v -v n )δ (1) (u -u n ) • u n b n - N n=1 v -1 n δ(v -v n )δ (1) (u -u n ) • c n . (14)
Note that this expression is linear in the source terms (a n , b n , c n ).

The next step is to apply the moment transform to Eq. ( 14). Formally, this yields

v k u l 1 u m 2 u p 3 (∂ t f + u • ∂ x f ) dv du = N n=1 (1 -k)v k n u l 1,n u m 2,n u p 3,n a n + N n=1 (k -l -m -p)v k-1 n u l 1,n u m 2,n u p 3,n b n + N n=1 v k-1 n u l 1,n u m 2,n u p 3,n lu -1 1,n c 1,n + mu -1 2,n c 2,n + pu -1 3,n c 3,n , (15)
where, unless otherwise noted, the definite integrals cover all of phase space. The next step is to consider the terms in Eq. ( 1) that correspond to transport in volume-velocity phase space.

Phase-space transport

We begin by rewriting Eq. (1) as

∂ t f + u • ∂ x f = P, (16)
where the phase-space transport terms are defined by

P ≡ -∂ v (R v f ) -∂ u • (Ff ) + Γ. (17)
We can then define the moment transform of the phase-space terms by

P (k, l, m, p) ≡ v k u l 1 u m 2 u p 3 P dv du. ( 18
)
Note that if the moments P (k, l, m, p) are known, Eq. ( 15) forms a linear system that can be solved to find the unknown source terms. We can compute the phase-space moments using Eq. ( 17):

P (k, l, m, p) = -v k u l 1 u m 2 u p 3 [∂ v (R v f ) + ∂ u • (Ff ) -Γ] dv du. ( 19
)
As shown next, the integrals on the right-hand side can be expressed in terms of the weights and abscissas, and a flux term corresponding to disappearance of droplets due to evaporation.

Starting with the evaporation term in Eq. ( 19), we can use integration by parts to find

∞ 0 v k ∂ v (R v f ) dv = -δ k0 R v (0, u)f (0, u) - ∞ 0 kv k-1 R v (v, u)f dv, (20)
where δ k0 is the Kronecker delta. Using Eq. ( 8) in the final integral, we find

v k u l 1 u m 2 u p 3 ∂ v (R v f ) dv du = -δ k0 ψ(t)u l f 1 u m f 2 u p f 3 - N n=1 kw n v k-1 n u l 1,n u m 2,n u p 3,n R v (v n , u n ), (21)
where ψ(t) is the evaporative flux of droplets at zero size and u f is the velocity of droplets with zero volume (which will normally correspond to the fluid velocity). Note that the first term on the right-hand side of Eq. ( 21) will be non-zero only for k = 0, and corresponds to the loss of droplets due to evaporation. A fundamental question when applying DQMOM to evaporation problems is how to determine ψ(t) from the weights and abscissas. The value of ψ(t) corresponds to the value of the number density function at zero size, and in the case of the d 2 evaporation law, it is precisely the value of the number density as a function of droplet surface, which has no reason to be zero in general. Determining the value of ψ(t), a pointwise information, from the values of moments is clearly a difficult task, for which we will propose a solution in the next subsection. On the other hand, the second term on the right-hand side of Eq. ( 21) is non-zero only for k > 0, and appears in closed form.

Turning next to the drag-force term in Eq. ( 19), we can use integration by parts to find

u l j ∂ u j (F j f ) du = -lu l-1 j F j f du for j = 1, 2, 3. (22)
Thus, the drag-force term becomes

v k u l 1 u m 2 u p 3 ∂ u • (Ff ) dv du = - N n=1 w n v k n u l 1,n u m 2,n u p 3,n lu -1 1,n F 1 (v n , u n ) + mu -1 2,n F 2 (v n , u n ) + pu -1 3,n F 3 (v n , u n ) . ( 23
)
Note that this term appears in closed form.

Turning now to the coalescence term, we will treat each of the two parts Q - coll and Q + coll separately. The first part yields in a straightforward manner

v k u l 1 u m 2 u p 3 Q - coll dv du = - N n=1 N q=1 w n w q v k n u l 1,n u m 2,n u p 3,n B(|u n -u q |, v n , v q ). (24)
The second part requires a change in the order of integration, and a change of variables:

v 0 h(v, u)B(|u ⋄ -u * |, v ⋄ , v * )f (v ⋄ , u ⋄ )f (v * , u * )J dv * dv du * du = ∞ v * h(v, u)B(|u ⋄ -u * |, v ⋄ , v * )f (v ⋄ , u ⋄ )f (v * , u * )J dv dv * du * du = h v * + v ⋄ , v * u * + v ⋄ u ⋄ v * + v ⋄ B(|u ⋄ -u * |, v ⋄ , v * )f (v ⋄ , u ⋄ )f (v * , u * ) dv * dv ⋄ du * du ⋄ , ( 25
)
where h is an arbitrary function of v and u. It then follows that

v k u l 1 u m 2 u p 3 Q + coll dv du = 1 2 N n=1 N q=1 w n w q (v n + v q ) k v n u 1,n + v q u 1,q v n + v q l × v n u 2,n + v q u 2,q v n + v q m v n u 3,n + v q u 3,q v n + v q p B(|u n -u q |, v n , v q ). ( 26
)
Note that the right-hand side of this expression is in closed form.

Collecting together all of the terms, the moments appearing on the right-hand sides of Eqs. [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF][START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF][START_REF] Laurent | Accurate treatment of size distribution effects in polydispersed spray diffusion flames: multi-fluid modeling, computations and experiments[END_REF][START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF][START_REF] Marchisio | Quadrature method of moments for aggregationbreakage processes[END_REF][START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF] become

P (k, l, m, p) = δ k0 ψ(t)u l f 1 u m f 2 u p f 3 + N n=1 kw n v k-1 n u l 1,n u m 2,n u p 3,n R v (v n , u n ) + N n=1 w n v k n u l 1,n u m 2,n u p 3,n lu -1 1,n F 1 (v n , u n ) + mu -1 2,n F 2 (v n , u n ) + pu -1 3,n F 3 (v n , u n ) + 1 2 N n=1 N q=1 w n w q   (v n + v q ) k v n u 1,n + v q u 1,q v n + v q l v n u 2,n + v q u 2,q v n + v q m v n u 3,n + v q u 3,q v n + v q p -v k n u l 1,n u m 2,n u p 3,n -v k q u l 1,q u m 2,q u p 3,q   B(|u n -u q |, v n , v q ). ( 27
)
Note that due to the form of the coalescence term, the moments conserve mass (P (1, 0, 0, 0) = 0) and momentum (P (1, 1, 0, 0) = P (1, 0, 1, 0) = P (1, 0, 0, 1) = 0) when evaporation and drag are null. Thus, the weights and abscissas in the DQMOM representation will keep the same conservation properties as the original model (i.e., as Eq. ( 1)).

Comparing the terms in Eqs. ( 15) and ( 27), we can note that the evaporation and drag terms in the DQMOM representation can be solved for explicitly. Thus, the source terms can be written as

b n = b * n + w n R v (v n , u n ), (28)
c n = c * n + w n u n R v (v n , u n ) + w n v n F(v n , u n ), (29)
where source terms a n , b * n and c * n in the transport equations are found by solving the linear system

N n=1 (1 -k)v k n u l 1,n u m 2,n u p 3,n a n + N n=1 (k -l -m -p)v k-1 n u l 1,n u m 2,n u p 3,n b * n + N n=1 v k-1 n u l 1,n u m 2,n u p 3,n lu -1 1,n c * 1,n + mu -1 2,n c * 2,n + pu -1 3,n c * 3,n = P * (k, l, m, p), (30)
with the right-hand side given by

P * (k, l, m, p) = -δ k0 ψu l f 1 u m f 2 u p f 3 + 1 2 N n=1 N q=1 w n w q   (v n + v q ) k v n u 1,n + v q u 1,q v n + v q l v n u 2,n + v q u 2,q v n + v q m v n u 3,n + v q u 3,q v n + v q p -v k n u l 1,n u m 2,n u p 3,n -v k q u l 1,q u m 2,q u p 3,q   B(|u n -u q |, v n , v q ). ( 31
)
The expression for the source terms (Eq. 30) completes the derivation of the DQMOM transport equations for the Williams spray equation.

In the absence of coalescence, Eq. ( 31) is particularly simple. Thus, the pure evaporation case is an interesting limit case for which a n , b * n , and c * n will be non-zero only if the evaporative flux ψ is non-zero. However, the evaporative flux cannot be determined by moment constraints alone (see Section 2.3). If the evaporative flux is assumed to be null, the zeroorder moment will remain unchanged in the absence of coalescence as long as some abscissa crosses the zero size limit and yields a pointwise singular and infinite flux as in Lagrangian methods when some parcels reach the zero size limit. However, as mentioned in the Introduction, since there are only a few abscissas that describe the moment dynamics, such a singular behavior is not ideal for smooth number density functions (whereas it is the correct one if the number density function is a sum of Dirac delta function from the beginning as in the bimodal case that will be studied later). Consequently we need an evaluation of this flux function that guarantees a smooth flux as a function of time for smooth distribution functions. Even when coalescence is included, the moments may be poorly estimated if the evaporative flux is neglected. An example of such behavior can be found in the work of Mossa [START_REF] Mossa | Extension polydisperse pour la description euler-euler des écoulements diphasiques réactifs[END_REF] where the droplet size distribution is presumed to be log-normal and where the evaporative flux at zero size is neglected, leading to numerical difficulties and a poor prediction of the second moment. Thus, we will use a separate procedure, described next, to approximate the contribution due to the evaporative flux that yields a continuous in time flux, as well as a guarantee that the abscissas never cross the zero size limit.

Evaporative flux

The source terms cannot be computed directly from the moment constraints in Eq. ( 31) because the evaporative flux is unknown. We must therefore apply additional (or different) constraints to determine all of the unknowns. Considering only evaporation and setting drag and coalescence to zero in the right-hand side of Eq. ( 30), we obtain the following linear system:

N n=1 (1 -k)v k n u l 1,n u m 2,n u p 3,n a n + N n=1 (k -l -m -p)v k-1 n u l 1,n u m 2,n u p 3,n b * n + N n=1 v k-1 n u l 1,n u m 2,n u p 3,n lu -1 1,n c * 1,n + mu -1 2,n c * 2,n + pu -1 3,n c * 3,n + δ k0 u l f 1 u m f 2 u p f 3 ψ = 0 (32) with 5N + 1 unknowns a n , b * n , c * n and ψ.
Note that because the right-hand side is null, only trivial solutions can be found using moment constraints. We will therefore introduce ratio constraints of the form

D Dt w n w n+1 evap = 0, D Dt v n v n+1 evap = 0 and D Dt u jn u jn+1 evap = 0,
which are applied only for the changes due to evaporation. These constraints are motivated by the behavior of the weights and abscissas corresponding to sufficiently smooth and continuous density functions. For example, if the surface density function is exponential and the evaporation rate is proportional to the surface area of a droplet, then the abscissas remain constant and the weights decrease monotonely. On the other hand, for singular density functions (e.g., composed of delta functions), the ratio constraints are expected to perform poorly. We will look more closely at this issue in Section 4.

It can be observed that the choice of k = l = m = p = 0 in Eq. ( 32) leads to

ψ = - N n=1 a n . (33)
Thus, the evaporative flux depends only on a n . Note that physically ψ ≥ 0. Hence, if the value computed for ψ from Eq. ( 33) is negative (which is possible for very general evaporation rates), then a n , b * n , c * n and ψ are set equal to zero. However, for the evaporation rate considered in this paper, it can be shown that the flux will be non-negative.

Conservation of mass (k = 1 and l = m = p = 0 in Eq. ( 32)) leads to

n b * n = 0. ( 34
)
Applying the ratio constraint for the abscissas yields

w n+1 v n+1 b * n -w n v n b * n+1 = E n for n = 1, . . . , N -1; (35)
where the right-hand side is defined by

E n = w n w n+1 [v n R v (v n+1 ) -v n+1 R v (v n )] . (36)
Note that in order for there to be an evaporative flux, we will normally have

E n ≥ 0 for all n (assuming that v 1 < v 2 < . . . < v N ). The case where E n = 0 occurs when R v (v) is proportional to -v (i.e.
, the evaporation rate is proportional to the droplet volume). The more common case where E n > 0 occurs when R v (v) is proportional to -v 1/3 (i.e., the droplet surface area decreases linearly). In general, R v (v) ∝ -v γ with γ < 1 will lead to positive E n . The physical interpretation for this difference is that for γ < 1 the droplets will disappear due to evaporation in a finite time, while for γ ≥ 1 the disappearance time is infinite. The linear system formed from Eqs. ( 34) and ( 35) can be solved separately to find b * n .

Conservation of momentum (k = 1 and l, m, or p = 1 in Eq. ( 32)) leads to

N n=1 c * n = 0. (37)
Likewise, the ratio constraint for each component of the velocity yields

w n+1 v n+1 u jn+1 c * jn -w n v n u jn c * jn+1 = u jn u jn+1 E n for n = 1, . . . , N -1. (38)
Together with Eq. ( 37), this equation can be solved separately for each component (j = 1, 2, 3) to find c * n .

The ratio constraint for the weights yields N -1 equations for a n :

w n+1 a n -w n a n+1 = 0 for n = 1, . . . , N -1. (39)
Note that this constraint is satisfied by a n = αw n where α is unknown. We must therefore choose one independent moment in Eq. ( 32) in order to solve for α. Since b * n and c * n are known, we can rearrange Eq. ( 32) as

α N n=1 (k -1)v k n u l 1,n u m 2,n u p 3,n + δ k0 u l f 1 u m f 2 u p f 3 w n = N n=1 (k -l -m -p)v k-1 n u l 1,n u m 2,n u p 3,n b * n + N n=1 v k-1 n u l 1,n u m 2,n u p 3,n lu -1 1,n c * 1,n + mu -1 2,n c * 2,n + pu -1 3,n c * 3,n , (40)
which can be solved with k = 1 to find α. If we choose, for example, k = 2 and l = m = p = 0 as the independent moment, then the constraint becomes

α = 2 N n=1 v n b * n N n=1 v 2 n w n (41)
and α depends only on b * n . However, if we choose k = 2 and l = m = p = 1, then the constraint becomes

α = N n=1 v n (u 2n u 3n c * 1n + u 1n u 3n c * 2n + u 1n u 2n c * 3n -u 1n u 2n u 3n b * n ) N n=1 v 2 n u 1n u 2n u 3n w n . (42)
For this choice, α is independent of u f . A choice that leads to a fully coupled system is k = 2, l = 2, m = p = 0, which yields

α = 2 N n=1 v n u 1n c * 1n N n=1 v 2 n u 2 1,n w n ( 43
)
or k = m = p = 0 and l = 1, which yields

α = N n=1 v -1 n u 1,n b * n -c * 1,n N n=1 (u 1,n -u f 1 ) w n . (44)
Note that when v n → 0, we have u 1,n → u f 1 and c * 1,n → u f 1 b * n ; hence, this last constraint is consistent with this limiting behavior. These choices are asymmetric in the velocity components, and thus do not treat all components the same. A "symmetric" choice with similar properties is k = 2 and l = m = p = 2 or k = 0 and l = m = p = 1, which lead to a more complicated constraint. The "best" choice will most likely be problem dependent. In our test cases, the choices with k = 2 give similar results, better than the ones with k = 0. The calculations are thus done with the value of α given in Eq. (41): this value is the simplest and can be shown to be non-positive as soon as E n ≥ 0, at least for the case N = 2.

In summary, the contribution due to evaporation is estimated by first solving separate linear systems for b * n and c * n . The estimate for a n = αw n is found using an independent moment constraint from Eq. ( 40) to find α. Finally, the evaporative flux ψ is computed from Eq. ( 33), and should be non-negative. If ψ is negative (or equivalently if α is positive), then the contribution due to evaporation is null. The contribution due to coalescence is found by solving a linear system of the form of Eq. ( 30) where the right-hand side is given by Eq. ( 31) with ψ = 0. As described below, the final source terms (a n , b * n , c * n ) are found simply by adding together the contributions from the evaporative flux and coalescence.

DQMOM linear system

The DQMOM representation of Williams spray equation is given by the transport equations for the weights and abscissas (Eqs. [START_REF] Dukowicz | A particle-fluid numerical model for liquid sprays[END_REF][START_REF] Dupays | Some aspects of twophase flows in solid-propellant rocket motors[END_REF][START_REF] Fox | Computational Models for Turbulent Reacting Flows[END_REF]). The source terms for these equations are found by solving the linear system as described above. The exact form of the linear system depends on the choice of moments. This choice, in turn, will determine if the system is well defined in the sense that the coefficient matrix is non-singular. A choice of moments that is consistent with the mono-kinetic assumption used in the multi-fluid model is to consider only moments of orders zero and one in the velocity components (i.e., l, m, p ∈ {0, 1}). In this work, in order to make direct comparisons with the multi-fluid model, we will limit our consideration to such moments. In general, this choice of moments should allow for the best possible description of droplet coalescence, while at the same time ensuring that droplet mass and momentum are conserved.

A choice of 5N moments that has been found to be non-singular is

k = (i -1)/3 i ∈ {1, . . . , 2N} with l = m = p = 0 k = i i ∈ {1, . . . , N} with l = 1, m = p = 0 k = i i ∈ {1, . . . , N} with m = 1, l = p = 0 k = i i ∈ {1, . . . , N} with p = 1, l = m = 0. (45)
For N ≥ 2, this choice of moments includes the surface area and the volume of the droplets, which are important variables for evaporating spray, as well as their momentum. The linear system can then be written in matrix form (showing only non-zero components) as

              A 1 A 2 E 1 E 2 E 3 A 3 A 4 B 1 C 1 D 1 B 2 C 2 D 2 B 3 C 3 D 3                             a b * c * 1 c * 2 c * 3               =               P a P b P 1 P 2 P 3               , (46)
where the matrices A j , B j , C j , D j and E j are all N×N, and a, b * and c * j are column vectors formed from the components a n , b * n and c * jn , respectively. In general, the exact definitions of the other matrices will depend on which constraints are used to define the system, i.e., Eq. [START_REF] O'rourke | Collective drop effects on vaporizing liquid sprays[END_REF] or those described in Section 2.3. Nevertheless, the form of the linear system is the same in all cases. As noted earlier, the linear system is solved twice at each time step. First with the matrices for the evaporative flux without coalescence (i.e., A 3 = B j = C j = 0), and second with the matrices for coalescence without evaporation (i.e., E j = 0). The unknowns a, . . ., c * 3 are found by adding the two solutions.

As discussed earlier, for the evaporative step the linear system can be decomposed into five N×N systems that can be solved sequentially. Likewise, for the coalescence step a and b * can be found separately by solving a 2N×2N system:

   A 1 A 2 A 3 A 4       a b *    =    P a P b    , (47)
and then each of the vectors c j can be found separately:

D j c * j = P j -B j a -C j b * , (48)
where (for coalescence) D j = V is a Vandermonde matrix [START_REF] Press | Numerical Recipes in Fortran 77[END_REF] formed from the volume abscissas v n . Other choices of moments have also been found to be numerically stable. For example, another valid choice is

k = (i -1)/3 i ∈ {1, . . . , 2N} with l = m = p = 0 k = (2i -1)/3 i ∈ {1, . . . , N} with l = 1, m = p = 0 k = (2i -1)/3 i ∈ {1, . . . , N} with m = 1, l = p = 0 k = (2i -1)/3 i ∈ {1, . . . , N} with p = 1, l = m = 0. ( 49
)
This choice can be found to give more accurate results for our test cases and still includes the surface area and the volume of the droplets, as well as their momentum. Thus, it will be used for the computations in Section 4. We should note that for a given value of N, the simulation results found with the moment sets in Eqs. ( 45) and ( 49) were nearly identical. The choice between these two systems was thus made based on ease of solution of the linear system.

In general, when moments involving the velocity are limited to first order, the matrices that must be inverted will be non-singular as long as the volume abscissas are distinct. The numerical treatment of the singularities associated with Eq. (47) has been discussed elsewhere [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF]. The coalescence operator will normally force the v n to remain distinct if they have distinct velocities. However, if due to initial conditions two or more of the volume abscissas are equal, it suffices to perturbate the values of v n enough to allow for the coefficient matrix in Eq. ( 47) to be invertible. We should also note that for cases dominated by coalescence (e.g., without evaporation) the volume abscissas grow rapidly, leading to matrices that are more and more ill-conditioned as the abscissas increase. Thus, even though they are strictly nonsingular, such matrices lead to severe numerical difficulties. However ill-conditioning can be almost completely alleviated by using iterative improvements of the linear solver as described in Section 2.5 of Press et al. [START_REF] Press | Numerical Recipes in Fortran 77[END_REF] after rescaling Eq. ( 30). The latter is done by defining positive scaling factors v s and u s , and dividing both sides of Eq. ( 30) by v k s u l+m+p s . Note that the abscissas and unknown source terms are rescaled in a consistent manner:

v n → v n /v s , u n → u n /u s , a → a, b * → b * /v s , and c * j → c * j /(v s u s ).
The evaporative flux constraints (Eqs. [START_REF] Tanguy | Development of a level set method for interface tracking: application to droplet collisions[END_REF], (38) and (40)) can be rescaled in a similar manner by introducing a positive scaling factor w s for the weights: w n → w n /w s . In this work, we use the following scaling factors: v s = max n v n , u s = max n |u n | and w s = n w n . We find that using the scaled variables and at most three iterative improvements of the linear solver are enough to completely eliminate round-off error in the solution to the DQMOM linear system. Moreover, because round-off error leads to poor performance of the differential equation solver, the overall computational cost using the iterative improvements can be significantly reduced.

Nozzle Test Problem

In order to validate the proposed DQMOM approach for strongly coalescing sprays, and to compare this method to both a reference Lagrangian solver solution as well as the solution obtained with the multi-fluid model, we need a well-suited test problem that is difficult enough to highlight the limitations of the methods under consideration. For that purpose, we have chosen for the gas phase a 2D axisymmetrical conical decelerating nozzle, designed in such a way that it admits, for one-way coupling spray dynamics a self-similar solution. After presenting the details of this configuration, we will provide the set of DQMOM equations to be solved in this framework. We have selected six representative test cases, combining coalescence/no coalescence with evaporation/no evaporation, which are then presented. Next we give an overview of the Lagrangian solver that provides the reference solution for the various test cases. Because the problems under consideration can be difficult to solve numerically, we must be very careful as far as this reference solution in concerned and thus we provide the details of the Lagrangian numerical integration in the limit of one-way coupling with the gas phase. Finally, before discussing the results in Section 4, we review the fundamentals of the Eulerian multi-fluid model for the sake of self-consistency of the paper.

Definition of configuration

The chosen configuration is stationary 2D axisymmetrical in space and 1D in droplet size. It is described in detail, along with the Lagrangian solver, 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.

A spray of pure heptane fuel is carried by a gaseous mixture of heptane and nitrogen into a conical diverging nozzle of axis (0 < z). At the entrance, z = z 0 , 99% of the mass of the fuel is in the liquid phase, whereas 1% is in the gaseous mixture. The temperature (400 K) as well as the composition of the gas mixture (mass fraction is 2.9% for heptane and 97.1% for nitrogen) is fixed during the entire calculation. The gas density is then 0.871733 mg.cm -3 . The influence of the evaporation process on the gas characteristics is not taken into account in our one-way coupled calculation. It is clear that the evaporation process is going to change the composition of the gas phase and then of the evaporation itself. However, we do not attempt to achieve a fully coupled calculation, but only to compare two ways of evaluating the coupling of the dynamics, evaporation and coalescence of the droplets. It has to be emphasized that it is not restrictive in the framework of this study, which is focused on the numerical validation of Eulerian solvers for the liquid phase under conditions of strong coalescence.

For the problem to be one-dimensional in space, conditions for straight trajectories are used and are compatible with the assumption of an incompressible gas flow. This leads to the following expression for the gaseous axial velocity v z and the reduced radial velocity v r /r:

v z = V (z) = z 2 0 V (z 0 ) z 2 , v r r = U(z) = V (z) z = z 2 0 V (z 0 ) z 3 for z ≥ z 0 ( 50
)
where z 0 > 0 is the coordinate of the nozzle entrance and the axial velocity V (z 0 ) at the entrance is fixed. The trajectories of the droplets are also assumed straight since their injection velocity is co-linear to the one of the gas. This assumption is only valid when no coalescence occurs. However, even in the presence of coalescence, it is valid in the neighborhood of the centerline.

Let us finally consider two droplet distribution functions. The first one, called monomodal, is composed of droplets with radii between 0 and 35 µm, with a mean radius of 12 µm, a variance of 5 µm and a Sauter mean radius of 15.6 µm. It is represented in Fig. 1 and is typical of the experimental conditions reported in [START_REF] Laurent | Accurate treatment of size distribution effects in polydispersed spray diffusion flames: multi-fluid modeling, computations and experiments[END_REF]. The droplets are constituted of liquid heptane, their initial velocity is the one of the gas, their initial temperature, fixed at the equilibrium temperature 325.4 K (corresponding to an infinite conductivity model), does not change along the trajectories. The second distribution is called bimodal since it involves only two groups of radii, respectively, 10 and 30 microns with equal mass density. This bimodal distribution function is typical of alumina particles in solid propergol rocket boosters [START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF], and is represented in Fig. 1.

The initial injected mass density is taken as m 0 = 3.609 mg.cm -3 so that the volume fraction occupied by the liquid phase is 0.57%. Because of the deceleration of the gas flow in the conical nozzle, droplets are going to decelerate, however at a rate depending on their size and inertia. This will induce coalescence. The deceleration at the entrance of the nozzle is taken as a(z 0 ) = -2V (z 0 )/z 0 ; it is chosen large enough so that the velocity difference developed by the various sizes of droplets is important. We have chosen rather large values, as well as strong deceleration, leading to extreme cases: V (z 0 ) = 5 m.s -1 , z 0 = 10 cm for the monomodal case and V (z 0 ) = 5 m.s -1 , z 0 = 5 cm for the bimodal case. These values generate a very strong coupling between coalescence, evaporation and droplet dynamics. These severe conditions as well as the two types of size distributions make the test cases under consideration very efficient tools for the numerical evaluation of the two Eulerian models.

DQMOM model equations in nozzle configuration

For the nozzle test case, Eqs. [START_REF] Dukowicz | A particle-fluid numerical model for liquid sprays[END_REF][START_REF] Dupays | Some aspects of twophase flows in solid-propellant rocket motors[END_REF][START_REF] Fox | Computational Models for Turbulent Reacting Flows[END_REF] reduce to a set of ordinary differential equations (ODEs) defined on the interval z ∈ [z 0 , ∞) for the variables w n , v n , η n (z) = u r /r and ξ n (z) = u z :

2w n η n + ∂ z (w n ξ n ) = a n , (51) 2w n v n η n + ∂ z (w n v n ξ n ) = b n , (52) 3w n v n η 2 n + ∂ z (w n v n η n ξ n ) = c rn /r, (53)
2w n v n η n ξ n + ∂ z w n v n ξ 2 n = c zn , (54)
where u z = ξ(z) and u r = rη(z) are the axial and radial components of the spray velocity, respectively. The corresponding fluid velocities are given in Eq. ( 50). The terms on the righthand side of Eqs. (52-54) are given by

b n = b * n + w n R v (v n ), (55)
c rn /r = c * rn + w n η n R v (v n ) + w n v n F r (v n , η n )/r, (56)
c zn = c * zn + w n ξ n R v (v n ) + w n v n F z (v n , ξ n ), (57)
where the drag model is

F r (v, η)/r = α 4π 3v 2/3 (U -η), (58)
F z (v, ξ) = α 4π 3v 2/3 (V -ξ) (59)
with α = 1.566 × 10 -7 m 2 .s -1 .

From the form of the governing equations, it is straightforward to show that if η n = ξ n /z at z = z 0 , then this relation will hold for all z and the droplet trajectories are straight lines.

The system of DQMOM model equations can thus be reduced to three nonlinear ODEs for w * n = w n (z/z 0 ) 2 , v n , and ξ n by eliminating Eq. (53):

∂ z (w * n ξ n ) = a n , ( 60
) ∂ z (w * n v n ξ n ) = b * n + w * n R v (v n ) (61)
and

∂ z w * n v n ξ 2 n = c * zn + w * n ξ n R v (v n ) + αw * n v n 4π 3v n 2/3 (V -ξ n ). ( 62
)
The terms on the left-hand side represent changes in the weights and abscissas due to transport. The terms on the right-hand side represent, respectively, the changes due to coalescence, evaporation, and drag. The coalescence terms are found by solving

(1 -k) N n=1 v k n ξ m n a n + (k -m) N n=1 v k-1 n ξ m n b * n + m N n=1 v k-1 n ξ m-1 n c * zn = 1 2 z 0 z 2 N n=1 N q=1 w * n w * q B(|ξ n -ξ q |, v n , v q ) (v n + v q ) k v n ξ n + v q ξ q v n + v q m -v k n ξ m n -v k q ξ m q . ( 63
)
Note the presence of the scaling factor (z 0 /z) 2 in the coalescence rate. As discussed in Section 2.4, we will use moments given in Eq. ( 49) that decouple Eq. ( 63) into two smaller systems.

Test cases

For evaporation, we will consider three cases described below: (i) no evaporation (R v = 0), (ii) linear evaporation (R v ∝ v), and (iii) non-linear evaporation (R v ∝ v 1/3 ). For each case, we will consider two sub-cases: without coalescence (E coal = 0) and with coalescence (E coal = 1). The two evaporation laws correspond to the two cases described in Section 2.3, for which droplets disappear either in infinite time (ii), thus leading to a evaporative flux at zero size, or in finite time (iii) for which the evaporative flux depends on the structure of the number density function in size phase space.

No evaporation

For the special case of no evaporation and no drag, the right-hand sides of Eqs. (52-54) are null. This special case has an analytical solution with v n , w * n , and ξ n constant. In the opposite limit of no evaporation and infinite drag, ξ n = V and w * n ∝ (z/z 0 ) 2 .

For non-evaporating droplets, R v = 0. In the absence of coalescence, a n = b * n = c * zn = 0. The DQMOM model reduces to v n and w * n constant, and

ξ n ∂ z ξ n = α 4π 3v n 2/3 (V -ξ n ). ( 64
)
This result is consistent with our earlier remark concerning the cases of zero and infinite drag. Finally, we should note that even with coalescence the momentum is conserved (k = m = 1) so that c * zn = 0. Thus, we can expect w * n ξ n to be approximately constant for all values of drag. For this case we expect excellent agreement between DQMOM and the Lagrangian solver in the absence of coalescence since the corresponding transport equations are identical (i.e., each DQMOM abscissa behaves like a Lagrangian particle). On the other hand, with coalescence the droplets grow very large and we expect differences due to how the coalescence term is treated in each method. This test case will, however, be very difficult for the multifluid model, since it was especially designed to tackle the problem of evaporation. In the presence of strong growth of droplet size, the number of sections that must be used in order to reproduce the physics with the multi-fluid model will also dramatically increase. Consequently, this test case will allow us to both test the capability of the DQMOM to capture the coupling of dynamics and coalescence at low cost in comparison to the Lagrangian solution, and to see if the multi-fluid model can provide good results, even if it is not competitive in terms of computational efficiency.

Linear evaporation

For evaporating droplets with linear evaporation, we take

R v (v n ) = -E v v n , (65)
with E v = 7.1262 s -1 for the monomodal case and E v = 14.2524 s -1 for the bimodal case. For this case, the evaporative flux ψ is zero. The coalescence terms are again found by solving Eq. ( 63). In the absence of coalescence, we have a n = b * n = c * zn = 0 and the DQMOM model reduces to w * n ξ n constant, Eq. ( 64), and

ξ n ∂ z v n = R v (v n ). (66)
Thus the volume v n and velocity ξ n are coupled through evaporation and drag, but are independent of w * n in the absence of coalescence. For this case we again expect excellent agreement between DQMOM and the Lagrangian solver in the absence of coalescence since the corresponding transport equations are identical. On the other hand, with coalescence there is a competition between growth and evaporation leading to smaller droplets than in the nonevaporating case. This is a very interesting test case, since it will allow us to compare both methods in an evaporative configuration, but without getting into the difficulty of modeling the droplet disappearance with the DQMOM approach.

Non-linear evaporation

With non-linear evaporation we will use

R v (v n ) = - E s 2 3v n 4π 1/3 (67)
with E s = 1.99 × 10 -7 m 2 /s. For this case the evaporative flux ψ will generally be non-zero, and is found using the method described in Section 2.3 with w * n in place of w n . However, we will also compare predictions for the bimodal initial distribution found by setting ψ = 0. As for the previous cases, we will investigate the effect of the flux model with and without coalescence. From a practical standpoint, the behavior of DQMOM with non-linear evaporation is of great interest and it is a configuration with which the comparison of both Eulerian models will be of practical relevance.

Reference Lagrangian solution

Euler-Lagrange numerical methods are commonly used for the calculation of polydisperse sprays in various application fields (see for example [START_REF] O'rourke | Collective drop effects on vaporizing liquid sprays[END_REF][START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF][START_REF] Miller | Direct numerical simulation of a confined three-dimensional gas mixing layer with one evaporating hydrocarbon-droplet-laden stream[END_REF][START_REF] Miller | Direct numerical simulation and subgrid analysis of a transitional droplet laden mixing layer[END_REF][START_REF] Rüger | Euler/Lagrange calculations of turbulent sprays: the effect of droplet collisions and coalescence[END_REF][START_REF] Dupays | Some aspects of twophase flows in solid-propellant rocket motors[END_REF] and the references therein). In this kind of approach, the gas phase is generally computed using a deterministic Eulerian solver, while the dispersed phase is treated in a Lagrangian way. The influence of the droplets on the gas flow is taken into account by the presence of source terms in the system of gas conservation equations. Two Lagrangian methods can be used as far as the dispersed phase is concerned depending on the level at which the physical processes are modeled. The first one is a Discrete Particle Simulation, which describes the evolution of numerical particles, each one representing one or several droplets. The physical processes such as transport, evaporation, collisions are then described by Liouville equations and the Eulerian fields usually recovered through ensemble averages. However, in the present study, we have preferred the Williams governing equation and thus a statistical description of the coalescence process. We then coherently use a Direct Simulation Monte-Carlo method (DSMC), the second kind of approach. It can be seen as the uncoupling, over a small time step, of the droplet transport in phase space (dynamics and evaporation), described by a particle method, and the collisions described by a Monte-Carlo method.

A complete exposition on the derivation and implementation of this method is outside the scope of this paper. We refer the reader, for example, to [START_REF] Amsden | Kiva II, a computer program for chemically reactive flows with sprays[END_REF][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] for more details. Here, for the sake of completeness, we present only the main features of the numerical method that we used in order to provide a "reference numerical solution" for the test cases.

Lagrangian solver

The Lagrangian solver can be roughly interpreted as a stochastic representation of the kinetic equation [START_REF] Alexander | The direct simulation Monte Carlo method[END_REF]. In other words, in the limit of a sufficiently large number of stochastic particles and a sufficiently fine computational grid (at least in the case of one-way coupling), the statistical estimates for the moments found from the particles should converge to those computed from the Eq. ( 1). In the Lagrangian solver, at each time step k, the droplet distribution function f (t k ) is approximated by a finite weighted sum of Dirac masses, f(t k ), which reads

f(t k ) = N k i=1 n k i δ z k i ,u k i ,v k i . (68)
Each weighted Dirac mass is generally called a "parcel" and can be physically interpreted as an aggregated number of droplets (the weight n k i ), located around the same point, z k i , with about the same velocity, u k i and about the same volume, v k i . N k denotes the total number of parcels in the computational domain at time t k . In all our calculations, the weights n k i were chosen in such a way that each parcel represents the same volume of liquid (

n k i v k i = Const).
Each time step of the particle method is divided in two stages. The first is devoted to discretization of the left-hand side of the kinetic equation ( 1), modeling the motion and evaporation of the droplets. In our code, the new position, velocity and volume of each parcel are calculated according to the following numerical scheme:

u k+1 i = u k i exp(-∆t/τ k i ) + V (z k i ) 1 -exp(-∆t/τ k i ) (69)
v k+1 i v k i dv/R v (v) = ∆t (70)
z k+1 i = z k i + ∆t u k+1 i = z k i + ∆t V (z k i ) + ∆t (u k i -V (z k i )) exp(-∆t/τ k i ) (71)
where V denotes the axial gas velocity, z k i (u k i ) corresponds to the axial coordinate of the position (velocity) of the parcel i at time t k , R v is the evaporation rate (independent of t and x since the gas composition and temperature are assumed constant in the domain). Eq. ( 70) is resolved analytically and depends on the chosen evaporation model. For linear evaporation, it can be written as

v k+1 i = v k i exp(-E v ∆t) (72)
and for non-linear evaporation, it is written as

s k+1 i = s k i -E s ∆t (73)
where s k i is the parcel surface area. The parcel relaxation time τ k i is defined as

τ k i = 2ρ(r k i ) 2 9µ g , (74)
with r k i being the parcel radius, ρ the liquid density and µ g the gas viscosity.

In Eqs. (69-71), the parcel radial coordinate is not calculated because the trajectories of the parcels are straight lines. Besides, as mentioned above, the influence of the droplets on the gas flow is not taken into account. Hence, Eq. ( 50) is used to calculate the gas velocity, V (z k i ), at the parcel location.

The second stage of a time step is devoted to the discretization of the collision operator. Several Monte-Carlo algorithms have been proposed in the literature for the treatment of droplet collisions [START_REF] O'rourke | Collective drop effects on vaporizing liquid sprays[END_REF][START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF][START_REF] Rüger | Euler/Lagrange calculations of turbulent sprays: the effect of droplet collisions and coalescence[END_REF][START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF][START_REF] Schmidt | A new droplet collision algorithm[END_REF]. They are all inspired by the methods used in molecular gas dynamics [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF] and, in particular, they suppose that the computational domain is divided into cells, or control volumes, which are small enough to consider that, within them, the droplet distribution function is almost uniform.

The algorithm used in our reference Lagrangian solver is close to the one proposed by O'Rourke. It consists of the following 3 steps (see also [START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF] for more details):

1. For each computational cell C J , containing N J parcels, we choose randomly, with a uniform distribution law, N J /2 pairs of parcels, (N J -1)/2 if N J is odd. 2. For each pair p, let p 1 and p 2 denote the two corresponding parcels with the convention n 1 ≥ n 2 , where n 1 and n 2 denote the parcel numerical weights. Then for each pair p of the cell C J , we choose randomly an integer ν p , according to the Poisson distribution law:

P (ν) = λ 12 ν! exp (-λ 12 ), with λ 12 = π n 1 (N J -1)∆t vol(C J ) (R 1 + R 2 ) 2 |u 1 -u 2 |
with vol(C J ) being the volume of cell C J , which is proportional to (z J /z 0 ) 2 for the nozzle test case, and R 1 , R 2 being the radii of the parcels p 1 , p 2 . The coefficient λ 12 represents the mean number of collision, during (N j -1) time steps, between a given droplet of parcel p 2 and any droplet of parcel p 1 . Note that a given pair of parcels is chosen, on average, every (N j -1) time steps.

3. If ν p = 0, no collisions occur during this time step between the parcels p 1 and p 2 . Otherwise, if ν p > 0, parcel p 1 undergoes ν p coalescence with parcel p 2 and the outcome of a collision is treated as follows. First the weight n 1 of the parcel p 1 is replaced by n ′ 1 = n 1 -ν p n 2 and its other characteristics are left unchanged. If n ′ 1 ≤ 0, parcel p 1 is removed from the calculation. Secondly, the velocity u 2 and the volume v 2 of parcel p 2 are replaced by

v ′ 2 = v 2 + ν p v 1 , u ′ 2 = v 2 u 2 + ν p v 1 u 1 v 2 + ν p v 1 ,
and its weight, n 2 , is left unchanged.

Let us mention that, for each time step and each control volume C J , the computational cost of this algorithm scales like O(N J ). This is a great compares to the O'Rourke method, which scales like O(N 2 J ). Another algorithm, with the same scaling features, has been introduced by Schmidt and Rutland in [START_REF] Schmidt | A new droplet collision algorithm[END_REF].

To obtain good accuracy, the time step, ∆t, must be chosen small enough to ensure that the number of collisions between two given parcels, p 1 and p 2 , is such that for almost every time: ν p n 2 ≤ n 1 . The average value of ν p being λ 12 , this constraint is equivalent to the condition

n 2 N J ∆t vol(C J ) π(R 1 + R 2 ) 2 |u 1 -u 2 | ≪ 1. (75)
For the nozzle test case described above, this constraint reveals to be less restrictive than the "CFL" condition ∀i = 1, ....N,

|u i |∆t ∆z ≪ 1, (76)
with ∆z being the mesh size. This condition is necessary to compute accurately the droplet movement and in particular to avoid that a parcel goes through several control volumes during the same time step. This is essential in order to have a good representation of the droplet distribution function in each mesh cell.

Reference solution

The Lagrangian solver just described is used to provide reference solutions in stationary cases with and without coalescence. In order to obtain a converged solution, particular attention must be devoted to the choice of the number of parcels, the size of the cells, and the time step.

For cases without coalescence, the computational cells are only used to have spatially averaged quantities to compare with Eulerian results. Moreover, the stationary aspect of the problem allows averaging in time in order to obtain smooth solutions. For these reasons, the conditions on the number of parcels and on the size of the computational cells are not very restrictive in the absence of coalescence. The time step is only limited by the CFL-like condition (76) needed for the convergence, with a low value. This last condition is the most restrictive since the scheme used for the transport of the particles is first order. For our test cases, the time step must be 10 -6 s or smaller. For cases with coalescence, there are additional restrictions. First, the algorithm used for coalescence assumes that the droplet distribution function of the spray is nearly uniform over each computational cell. However, in the region with high gradients of the gas velocity, that is to say at the entrance of the nozzle, this distribution can change quickly and the size of the cells must be small enough to avoid numerical errors. Moreover, in order to properly describe the coalescence phenomenon in each cell with the stochastic algorithm, a sufficient number of parcels must be present in each cell, typically on the order of 50, with a minimum of 20 [START_REF] Alexander | The direct simulation Monte Carlo method[END_REF]. The smaller are the cells, the larger must be the number of parcels in the computational domain. The required size of the cells is evaluated for the case where the size distribution changes the most rapidly (the case without evaporation). We then employ a nonuniform space discretization with 130 cells, with smaller cells near the entrance of the nozzle defined using a uniform discretization for the variable z 3/10 . The number of parcels injected per second is given in Tables 1 and2.

Eulerian multi-fluid solver

Eulerian multi-fluid methods were developed as an alternative to Lagrangian methods for the simulation of polydisperse evaporating sprays. A complete derivation of such methods from the kinetic model is performed 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] for dilute sprays and in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] for sprays with coalescence. The principle of the method is quite different than the one used in DQMOM. Indeed, it can be considered as a finite-volume discretization in the droplet size phase space for moments of order 0 and 1 of the velocity distribution conditioned on size.

In laminar flows, it can be proven rigorously that it is sufficient to work with only these two moments as long as the velocity distribution conditioned on droplet size is mono-kinetic [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] (i.e. all droplets with the same volume have identical velocities so that the size-conditioned velocity variance is null.) By construction, the nozzle test problem will be mono-kinetic for non-coalescing droplets. However, with coalescence the size-conditioned velocity variance can be nonzero. Comparisons between the Lagrangian and Eulerian results in the presence of coalescence will therefore allow us to quantify the magnitude of the error caused by invoking the mono-kinetic assumption in the Eulerian models. (Recall that the choice of moments used in the DQMOM linear system is equivalent to the mono-kinetic assumption in the multi-fluid model.) In this section, we provide only the main points of the derivation of the multi-fluid model, as well as the underlying assumptions that are implied, and the resulting system of equations that will be solved.

The first step consists of writing equations for the two moments in velocity. This leads to the closed semi-kinetic model if the following assumption is made concerning the structure of f :

f (v, u; x, t) = n(v; x, t)δ(u -ū(v; x, t)).
In other words, the droplet velocity conditioned on the size is assumed to be Dirac delta function. In the case of a coalescing spray, the compatibility of such a condition with droplet coalescence is far from obvious; however, the semi-kinetic system of conservation equations can be obtained by using an asymptotic limit as presented in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF].

The second step consists of discretizing n(v) in sections [v (j-1) , v (j) ) and in integrating the semi-kinetic model over each section. This leads to a multi-fluid model (by using a presumed distribution κ (j) (v) in each section), thereby yielding a conservation equation on the moment associated with the mass density n(v; x, t) = m (j) (x, t)κ (j) (v) where

v (j) v (j-1)
ρvκ (j) (v)dv = 1.

In addition, only the averaged velocity is considered in each section, i.e. ū(v; x, t) = ū(j) (x, t) if v (j-1) ≤ v < v (j) . The resulting system can be found in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. It can be rewritten and simplified in the stationary, self-similar, 2D axisymmetrical configuration we are considering.

The resulting set of equations is

2m (j) u z (j) z + ∂ z (m (j) u z (j) ) = -(E (j) 1 + E (j) 2 )m (j) + E (j+1) 1 m (j+1) + C (j) m (77) 2m (j) u z (j) z 2 + ∂ z (m (j) (u z (j) ) 2 ) = -(E (j) 1 + E (j) 2 )m (j) u z (j) + E (j+1) 1 m (j+1) u z (j+1) + m (j) F (j) z + C (j) muz (78)
where u z (j) is the axial velocity, which only depends on z, and ru z (j) /z is the radial velocity, since the trajectories are straight lines. Moreover, E

1 and E (j) 2 are the "classical" pre-calculated vaporization coefficients [START_REF] Greenberg | On the origin of spray sectional conservation equations[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]:

E (j) 1 = -ρ v (j-1) R v (v (j-1) ) κ (j) (v (j-1) ) and E (j) 2 = - v (j) v (j-1) ρ R v (v) κ (j) (v) dv,
and

F (j) z (v (j) u , u z (j)
) is the axial component of the "classical" pre-calculated drag force [START_REF] Greenberg | On the origin of spray sectional conservation equations[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]:

F (j) z = v (j) v (j-1) ρ v F z (v, u z (j) ) κ (j) (v) dv where v (j) u =   v (j) v (j-1) v κ (j) (v) dv v (j) v (j-1) v 1/3 κ (j) (v) dv   3/2
.

The source terms associated with coalescence phenomenon in the mass and momentum equation, respectively, of the jth section read

C (j) m = -m (j) N k=1 m (k) V jk Q jk + I (j) i=1 m (o ⋄ ji ) m (o * ji ) V o ⋄ ji o * ji (Q ⋄ ji + Q * ji ), (79)
C (j) muz = -m (j) u z (j) N k=1 m (k) V jk Q jk + I (j) i=1 m (o ⋄ ji ) m (o * ji ) V o ⋄ ji o * ji u z (o ⋄ ji ) Q ⋄ ji + u z (o * ji ) Q * ji , (80)
where

V jk = |u z (j) -u z (k) | and the collision integrals Q jk , Q ⋄ ji and Q * ji do not depend on z. The disappearance integrals Q jk are evaluated on rectangular domains L jk = [v (j-1) , v (j) ] × [v (k-1) , v (k) ], whereas the appearance integrals, Q ⋄ ji and Q * ji , are evaluated on the diagonal strips D ⋄ * j = {(v ⋄ , v * ), v (j-1) ≤ v ⋄ + v * ≤ v (j) }/ ∪ N k=1
L kk , which are symmetric strips with respect to the axis v ⋄ = v * . These strips D ⋄ * j are divided into domains, denoted by X ji and the symmetric one, X sym ji , where the velocity of the partners is constant. The domains X ji and X sym ji are the intersection of D ⋄ * j with L kl , k > l, and L kl , k < l, respectively; their index is denoted i ∈ [1, I (j) ] and we define two pointers that indicate the collision partners for coalescence, at fixed i: o ⋄ ji = k and o * ji = l.

The coefficients used in the model, either for the vaporization process or the drag force E (j)

1 , E (j) 2 and F (j) z , j = [1, N] in Eqs. (77-78), or for the coalescence:

Q jk , j = [1, N], k = [1, N], k = j, Q ⋄ ji , Q * ji , j = [2, N], i = [1, I (j)
] in Eqs. (79-80) can be pre-evaluated from the choice of κ (j) in each section. The algorithms for the evaluation of this coefficients are provided in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. The distribution function is chosen constant as a function of the radius in the sections 1 to N and exponentially decreasing as a function of the surface in the last section, as done in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF].

Because only the one-way coupled equations are solved and since the structure of the gas velocity field is prescribed and stationary, we only have to solve the 1D ordinary differential Eqs. (77,78) for each section. The problem is then reduced to the integration of a stiff initial value problem from the inlet where the droplets are injected until the point where 99.9% of the mass has evaporated. The integration is performed using LSODE for stiff ordinary differential equations from the ODEPACK library [START_REF] Hindmarsh | ODEPACK, a systematized collection of ODE solvers[END_REF]. It is based on BDF methods [START_REF] Hairer | Solving ordinary differential equations[END_REF] (Backward Differentiation Formulae) where the space step is evaluated at each iteration, given relative and absolute error tolerances [START_REF] Hindmarsh | ODEPACK, a systematized collection of ODE solvers[END_REF]. The relative tolerance, for the solutions presented in the following, are taken to be 10 -4 , and the absolute tolerance are related to the initial amount of mass in the various sections, since it can vary of several orders of magnitude. Repeated calculations with smaller tolerances have proved to provide essentially the same solutions.

Results and Discussion

Simulations for the cases presented in the previous section were carried out with the Lagrangian method, the multi-fluid model, and DQMOM. Except for the cases without evaporation for which the multi-fluid method is not well suited (it requires a large number of sections and is only presented for comparison purposes), the Eulerian methods were solved using a initial-value solver for ODEs and required very short computational times (i.e., CPU secs) on a desktop computer. It is interesting to note that in the case without evaporation, the small computational cost still holds for the DQMOM approach.

In contrast, the time-and z-dependent Lagrangian simulations required several CPU hrs for each case. Because the DQMOM and multi-fluid results do not depend on time, it is not appropriate to compare the computing times directly. Nevertheless, it will generally be the case that using Eulerian methods will result in a substantial reduction in the computing time for solving the spray equation. Such a statement was studied in details in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] for unsteady calculations and the conclusions drawn from that paper are applicable to the two Eulerian methods presented here. Thus, the principal open question is whether or not the DQMOM results are of comparable accuracy to the multi-fluid model and to the more costly Lagrangian simulations. We will compare predictions for selected statistics from the three solution methods in order to answer this question. For the monomodal distribution and DQMOM resolution, several initial conditions will be used in the following an are presented in Table 3. The representative moments used to compare the three solution methods are the number density m 0 , the mass density m 1 , the average axial velocity difference between droplets and the gas phase u d , and the Sauter mean radius r 32 . They are defined by

Monomodal distribution Vol. moments, N =4 Rad. moments, N =4 Rad. moments, N =6 Rad. moments, N =8 n w n /N 0 r n w n /N 0 r n w n /N 0 r n w n /N 0 r n 1
m 0 = f (v, u) dv du, m 1 = ρvf (v, u) dv du, u d = ρvu z f (v, u) dv du m 1 -V, r 32 = (4π/3) 1/3 vf (v, u) dv du v 2/3 f (v, u) dv du .
With the DQMOM approach, these quantities are written

m 0 = N n=1 w n , m 1 = N n=1 w n ρv n , u d = N n=1 w n v n (ξ n -V )/m 1 , r 32 = (4π/3) 1/3 N n=1 w n v n N n=1 w n v 2/3 n .
And with the multi-fluid method, they are

m 0 = N j=1 m (j) v (j) v (j-1) κ (j) (v)dv, m 1 = N j=1 m (j) , u d = N j=1 m (j) (u z (j) -V )/m 1 , r 32 = (4π/3) 1/3 N j=1 m (j) v (j) v (j-1) vκ (j) (v)dv N j=1 m (j) v (j) v (j-1) v 2/3 κ (j) (v)dv .
Note that in practical applications, the mass density is a key quantity because it represents the total mass of liquid contained in the droplets. In the nozzle test case, the rate of coalescence is strongly dependent on the velocity difference between droplets, which we find to be strongly correlated with the average axial velocity difference. Indeed, if u d is not accurately captured, then we find that the predictions for all moments will degrade accordingly. In addition to the moments, we will also compare the mean droplet velocity conditioned on the radius u z |r at selected downstream locations, as well as the mass distribution function (ρvf ). For the DQ-MOM, the scaled weights will be used to represent the mass distribution function. Obviously, since the sum of the weights equals the area under the mass distribution function, the absolute value of the heights of the scaled weights is arbitrary. Nevertheless, the relative heights and the locations provide insight into how well the quadrature points represent the distribution function.

We should note that for the monomodal cases without coalescence, the results with no evaporation were essentially identical for all three solution methods. The results presented below for the monomodal case with linear evaporation are representative of the quality of the predictions for all cases without coalescence and no evaporation. Likewise, for the bimodal case without coalescence and with linear or no evaporation, DQMOM and the Lagrangian method were essentially identical. The multi-fluid method also yielded very good results for these cases if the number of sections was chosen large enough to mitigate the numerical diffusion in the size phase space associated to the description of evaporation that leads to broadening of the peaks. Nonetheless, because none of these cases revealed any unanticipated problems for any of the simulation methods, we will not discuss them further. Instead, we will primarily focus on cases that present particular challenges to one or more of the solution methods.

Monomodal case: linear evaporation without coalescence

We begin with a representative case where all three solution methods yield essentially identical results for all statistics. As noted in the discussion of the methods, for linear evap- oration without coalescence the DQMOM equation for each node is the same as the Lagrangian model. Thus, the only difference between the two solution methods is that the Lagrangian method uses many more particles to represent the spray than the DQMOM method (N = 4). For the monomodal distribution, the multi-fluid model does not require many sections (N = 10) to accurately capture cases with linear evaporation without coalescence. The simulation results for the three methods are shown in Fig. 2. It can be observed that the mass density and velocity difference predicted by the three methods are nearly identical. From the plot of the mass distribution function at z = 16 cm, we can see that the multi-fluid model with ten sections does a good job of capturing the Lagrangian mass distribution function.

Likewise, the DQMOM weights and abscissas follow the general shape of the Lagrangian mass distribution function. Finally, for the conditional velocity u z |r we see that all three methods produce the same curve. We should note that for cases without coalescence the Lagrangian simulations predict essentially no velocity dispersion about the conditional value.

In other words, conditional velocity fluctuations defined by u ′ (r) ≡ (u z -u z |r ) 2 |r 1/2 are null. This is exactly one of the necessary conditions evoked when deriving the multi-fluid model, which would explain why its predictions for this case are in excellent agreement with the Lagrangian method.

Monomodal case: nonlinear evaporation without coalescence

Cases with nonlinear evaporation result in a loss of droplets in finite time, which translates into a nonzero flux ψ(t) in DQMOM. For the monomodal case without coalescence, we expect the flux term to be a smooth function of t, and thus it cannot be neglected. In the multi-fluid model, the flux is computed directly from the shape of the first section (i.e., the section near the origin) and does not yield any difficulty.

In our multi-fluid simulations, we use the "optimal" choice of sections with N = 12 shown in Fig.

3 [START_REF] Laurent | Accurate treatment of size distribution effects in polydispersed spray diffusion flames: multi-fluid modeling, computations and experiments[END_REF]. Obviously, a finer discretization (larger N) could be used in the multi-fluid model to attain closer agreement with the Lagrangian method, but this would increase the computational cost. Note that the first section is represented by a constant slope, which corresponds to a constant flux level at each time step. For the DQMOM, we use N = 4 and the evaporative flux is computed using the ratio constraints introduced in Section 2.3. It can be noticed that the increase of N do not imply an increase of the number of conserved moments during the evaporation step since the number of ratio constraints is also increasing in the same way. The value of N is then conditioned by the capacity of the method to follow the dynamics of droplets of different sizes. Representative results for the three solution methods are shown in Fig. 4. In general, all three methods produce very similar predictions. From the number density, we can observe that DQMOM with the ratio constraints does a good job of predicting the loss of droplets due to evaporation. Likewise, the mass densities found from all methods are very close. We should note that for z > 20 cm the number of remaining droplets is very small and the statistics computed from the Lagrangian method are subject to statistical errors. Comparing the Sauter mean radii predicted by the three methods, we can observe that the agreement is generally satisfactory up to z = 20 cm. The DQMOM shows the largest deviation from the Lagrangian method at z = 20 cm due to errors in the flux model, but the agreement is still acceptable. The differences in the Sauter mean radius are reflected in the predictions of the velocity difference. In general, droplets with a larger radius will have a higher velocity difference. Thus, we see that initially the Sauter mean radius predicted by DQMOM is larger than that from the Lagrangian method, resulting in a slightly higher velocity difference at z = 12 cm. Later on (z > 15 cm) this trend is reversed. Finally, we can note that neglecting the flux in DQMOM yields poor predictions of number density since we can observe the artificial jumps in the number density related to the singular fluxes associated to one abscissa crossing the zero size limit, as well as the oscillating dynamics of the Sauter mean radius for this case.

Bimodal case: nonlinear evaporation without coalescence

By changing from the monomodal to the bimodal distribution, we change the nature of the initial distribution function and thus the nature of the numerical difficulties. For the multi-fluid model, the bimodal case is difficult because a relatively large number of sections (N = 30) is needed to capture the two peaks with acceptable numerical diffusion. The use of a secondorder method developed in [START_REF] Laurent | Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF] would reduce this number to around 10; however, it would still be difficult to describe Dirac delta function by a finite-volume approximation. On the other hand, this case is "optimal" for DQMOM because only two (N = 2) abscissas are required (one for each peak) and the flux is null, expect when a peak passes the origin. In Fig. 5 results from the three simulation methods are shown and it is clear that DQMOM performs extremely well for this case by setting ψ = 0. For example, the number density function shows step changes at z = 7.2 cm and 13.8 cm (i.e., when a peak passes the origin), and DQMOM exactly reproduces this behavior. With N = 30, the multi-fluid model does a good job of predicting the mass density. However, from the plots of number density and Sauter mean radius, we can observe the negative effects of numerical diffusion, which tends to smooth out the peaks in the distribution (the method has been shown to be first order in the droplet size discretization step in [START_REF] Laurent | Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF], where some higher-order methods have been proposed). Nevertheless, all three methods yield reasonable predictions for all of the cases without coalescence. We should note, however, that for more complicated initial distributions (e.g., delta functions combined with smooth functions) specifying the evaporative flux in DQMOM may be problematic.

Monomodal case: linear evaporation with coalescence

We now turn to the more difficult cases that include coalescence. As mentioned earlier, the coalescence of droplets with different volumes (and velocities) will lead to velocity dispersion (u ′ (r) > 0). Physically, this implies that two droplets with the same volume will have a nonzero probability of colliding (due to the difference in velocity). Thus, the rate of coalescence when u ′ (r) > 0 will be larger than when u ′ (r) = 0. Numerical approximations (such as the multi-fluid model) that assume u ′ (r) = 0 should therefore predict smaller droplets than the Lagrangian method. In Fig. 6 we present results for the three methods for linear evaporation (ψ = 0) with coalescence. ¿From the velocity difference we can observe that coalescence leads to a slower relaxation to the gas velocity due to formation of larger droplets than without coalescence. Note that in general all three methods predict similar results for the velocity difference. However, due the differences in the predictions of the shape of the mass distribution function, the multi-fluid model predicts slightly slower relaxation and the DQMOM slightly faster than is found with the Lagrangian method. Comparing with Fig. 1, we can observe that coalescence leads to much larger droplets than are present in the initial distribution function.

In general, the multi-fluid model predicts a slightly larger number of droplets above 80 µm than the Lagrangian method. Nevertheless, the predictions are in reasonably good agreement.

The predictions for the conditional velocity u z |r are also good. Finally, note that we used N = 6 with DQMOM, the reason for which will be discussed for a more difficult case in Section 4.7.

Bimodal case: linear evaporation with coalescence

We now consider a more difficult case where the initial distribution is bimodal. As discussed previously, the peaks in the distribution are difficult to resolve accurately in the multi-fluid model with a limited number of sections . When combined with coalescence, this has important consequences because numerical diffusion can lead to spurious coalescence of droplets with slightly different volumes (and hence velocities) as observed in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. For example, with the bimodal distribution with droplets of radii 10 and 30 µm, coalescence cannot produce droplets below 30 µm. However, spurious coalescence between droplets of radii near 10 µm leads to droplets in the range below 30 µm. We overcome this difficultly by using a large number of sections (N = 500) in the multi-fluid model. This number could also be reduced by using a second-order method for the evaporation such as the one of [START_REF] Laurent | Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF] but this is not the point we want to make with this configuration. Note that the same problem arises in the Lagrangian method when the spatial cell size ∆z is too large. While DQMOM does not suffer from spurious coalescence, the bimodal case is still difficult because the initially two-peak distribution will quickly form multiple peaks due to pair-wise collisions. With N = 6 in DQ-MOM, it is at best possible to represent six peaks. Results for the three methods are shown in Fig. 7 where it can be seen that the mass density and the velocity difference are reasonably well predicted by the multi-fluid model and DQMOM. From the mass distribution function at z = 11 cm, the multi-peak structure due to coalescence is quite evident, as is the slight numerical diffusion in the multi-fluid model (even with N = 500, but this is expected since this is a first-order method). Note that DQMOM with N = 6 has two abscissas at points corresponding to the major peaks (10 and 33 µm), and the remaining abscissas at points without major peaks. Comparisons of the conditional velocity predicted by the three methods are also quite favorable for this difficult case.

Monomodal case: nonlinear evaporation with coalescence

We will now consider the more physically relevant case of nonlinear evaporation. As discussed earlier, the evaporative flux for this case is nonzero, so we will need to model it in DQMOM. Here, we consider two models for ψ: (a) ratio constraints and (b) ψ = 0. Because the initial distribution is monomodal, we might expect that using ratio constraints is always a better choice. On the other hand, if coalescence is much faster than evaporation, it might happen that droplets grow faster than they disappear so that the evaporative flux is closer to zero. For the multi-fluid model, we use N = 15 sections. Results for the three methods are shown in Fig. 8. The number density illustrates the effect of the choice of ψ in DQMOM. With ψ = 0, the number density changes discontinuously whenever an abscissa passes the origin. However, DQMOM with ratio constraints yields predictions very similar to the other two methods. Likewise, the mass density is predicted to be very similar for all three methods; however, using zero flux with DQMOM is slightly worse. The predictions for the Sauter mean radius show opposing trends. In general, the multi-fluid model overpredicts the mean radius (i.e. predicts too much coalescence), while DQMOM underpredicts it. As before, for the DQMOM predictions, the results with the ratio constraints are best. The predictions for the velocity difference follows the same trend. As discussed in the next example, the differences observed between the Lagrangian method and the two Eulerian methods is likely due to the latter's inability to capture velocity dispersion. Moreover, we have used N = 6 with DQMOM since, as shown in Section 4.7, it is adequate to describe coalescence phenomenon for this particular set of moments.

Monomodal case: coalescence with no evaporation

In order to highlight the role of coalescence on determining the evolution of the number density function, we now consider a case with no evaporation. For this case, droplets will grow continuously due to coalescence, and velocity dispersion will enhance the collision rate and lead to even larger droplets. Because the multi-fluid model uses fixed sections, it is necessary to fix the maximum radius at 200 µm with N = 500 in order to capture the largest droplets at z = 30 cm. In contrast, the abscissas in DQMOM move in phase space to accommodate growth. Nevertheless, we can anticipate that the number of abscissas will affect the accuracy of the DQMOM predictions. In Fig. 10 it can be observed that when the number of moments increases, the accuracy of the DQMOM solution increases, from something almost ignoring the coalescence phenomenon with N = 2 to a saturation of the accuracy for N ≥ 8. Indeed, the accuracy of the DQMOM for the description of the coalescence is related to the accuracy of the approximation of the coalescence operator by the quadrature formula ( 24) and ( 26), which increases with N. Since the results are quite good and at a low cost (and the linear system is reasonably well conditioned), we will use N = 6 for comparisons with the other two methods.

As mentioned earlier, without coalescence or evaporation all three methods predict essentially identical results. In Fig. 9 the results for the pure coalescence case are shown. Notice that the mass density does not decrease to zero because there is no evaporation; however, it does change due to transport. From the velocity difference, we can see that the multi-fluid model and DQMOM overpredict the relaxation rate. As discussed previously, this is due to both methods underpredicting the mean droplet size. From the mass distribution functions at z = 22 cm, we can observe that the Lagrangian method has more droplets with radii above 80 µm than the multi-fluid model, which is consistent with the observed trend in the velocity difference. In order to explore the link between velocity dispersion and coalescence, we have computed 50% probability intervals for the conditional velocity. These are defined to be the values of v for which the conditional velocity PDF f (v|r) is fifty percent of its peak value. Note that the absence of velocity dispersion f (v|r) is a delta function centered at u z |r , so the width of the intervals is a measure of dispersion. From the plot of conditional velocity, we can note that for large droplets the velocity dispersion is significant. We can also note that using DQMOM essentially results in points along the curve u z |r , i.e., increasing N with the same choice of moments does not capture the velocity dispersion.

Conclusions

In this work, we have implemented DQMOM to treat the Williams spray equation that describes evaporation, acceleration and coalescence of liquid droplets in a laminar gas flow. The derivation of the DQMOM equations was shown to be a straightforward task, and resulted in a linear system for the source terms. The right-hand side of this linear system is non-zero only in the presence of coalescence or non-linear evaporation. The coefficient matrix depends on the choice of moments used in DQMOM.

We have compared this method, as well as the solution obtained with another Eulerian method: the multi-fluid model, to the reference solution produced by a classical Lagrangian solver. As far as coalescence phenomena are concerned, the efficiency of DQMOM has been shown to be better than the multi-fluid model due to its limited numerical diffusion in the size phase space, especially for the bimodal distribution function. However, as far as the evaporation process is concerned, it is comparable to the multi-fluid model, but still needs a further study in order to fully understand how to treat optimally the issue of the evaporative flux due to droplet disappearance. Although this issue has been so far neglected in the literature on moment methods, our study illustrates that it has an important effect on the moment dynamics.

The principal conclusion from this study is that DQMOM is numerically robust and straightforward to implement for the Williams spray equation and that it will be a very good candidate for more complex two-phase combustion applications once the issue of the evaporative flux

Figure 1 .

 1 Figure 1. Initial number density functions for droplet radius. Left: Monomodal distribution. Right: Bimodal distribution.

Figure 2 .

 2 Figure 2. Monomodal case with linear evaporation (R v ≈ 7.126v). Top left: mass density. Top right: velocity difference. Bottom left: mass distribution function (z = 16 cm). Bottom right: conditional velocity (z = 16 cm).

Figure 3 .

 3 Figure 3. Monomodal distribution function with optimal sections.

Figure 4 .

 4 Figure 4. Monomodal case with nonlinear evaporation. Top left: number density. Top right: mass density. Bottom left: Sauter mean radius. Bottom right: velocity difference.

Figure 5 .

 5 Figure 5. Bimodal case with nonlinear evaporation. Top left: number density. Top right: mass density. Bottom left: Sauter mean radius. Bottom right: velocity difference.

Figure 6 .

 6 Figure 6. Monomodal case with linear evaporation (R v ≈ 7.126v) and coalescence. Top left: mass density. Top right: velocity difference. Bottom left: mass distribution function (z = 22 cm). Bottom right: conditional velocity (z = 22 cm).

Figure 7 .

 7 Figure 7. Bimodal case with linear evaporation (R v ≈ 14.252v) and coalescence. Top left: mass density. Top right: velocity difference. Bottom left: mass distribution function (z = 11 cm). Bottom right: conditional velocity (z = 11 cm).

Figure 8 .

 8 Figure 8. Monomodal case with nonlinear evaporation and coalescence. Top left: number density. Top right: mass density. Bottom left: Sauter mean radius. Bottom right: velocity difference.

Figure 9 .

 9 Figure 9. Monomodal case with coalescence but no evaporation. Top left: mass density. Top right: velocity difference. Bottom left: mass distribution (z = 22 cm). Bottom right: conditional velocity (z = 22 cm).

Figure 10 .

 10 Figure 10. Monomodal case with coalescence but no evaporation. Top left: mass density. Top right: error on the mass density. Bottom left: velocity difference. Bottom right: error on the velocity difference.

Table 2

 2 Number of parcels for the Lagrangian simulations for the cases with coalescence.

	distribution evaporation No. of parcels No. of parcels inj./s
	monomodal	linear		41,560	100,000
	monomodal	nonlinear	20,440	1,000,000
	bimodal	nonlinear	6,320	200,000
	Table 1
	Number of parcels for the Lagrangian simulations for the cases without coalescence.
	distribution evaporation No. of parcels No. of parcels inj./s min. No. of
						parcels /cell
	monomodal	linear	160,000	200,000	40
	bimodal	linear	126,000	560,000	50
	monomodal	nonlinear	35,000		1,300,000	260
	monomodal	no		44,200		300,000	65

Table 3

 3 Initial conditions for weights and abscissas found using QMOM.

		0.7323	9.9955	0.1845	4.4079 8.5573E-2 3.3423 4.6445E-2 2.8465
	2	0.2545	18.5282	0.5397	11.0409	0.2779	7.5262	0.1488	5.5373
	3 1.288E-2 27.5630	0.2635	18.2840 5.5339E-2 12.9743	0.3089	9.6916
	4 2.279E-4 36.0142 1.212E-2 28.3910 4.9778E-3 18.8823	0.3438	14.2697
	5					3.1137E-4 26.3693	0.12931	19.2986
	6					1.6671E-5 34.7171 2.0905E-2 25.2866
	7							1.6982E-3 31.5808
	8							6.5627E-5 37.5149

Acknowledgments

We would like to thank P. Villedieu for providing valuable help for the Lagrangian solver and several helpful discussions. The present research was done thanks to a Young Investigator Award from the French Ministry of Research (New Interfaces of Mathematics -M. Massot, 2003-2006), the support of the French Ministry of Research (Direction of the Technology) in the program: "Recherche Aéronautique sur le Supersonique" (Project coordinator: M. Massot 2003-2006) and an ANR (National Research Agency -France) Young Investigator Award (M. Massot, 2006-2009). One of the author (R. O. Fox) was partially supported by the U. S. National Science Foundation (CTS-0403864). The strong support of the Ecole Centrale Paris through the invitation of R. O. Fox to Visiting Professor positions in 2005 and 2006 is gratefully acknowledged.

, supported by

is further improved.