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

Frédérique Laurent-Nègre

Angelo Murrone
email: angelo.murrone@onera.fr

Joël Dupays
email: joel.dupays@onera.fr

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

Optimal Eulerian model for the simulation of dynamics and coalescence of alumina particles in solid propellant combustion

Keywords: Dilute Polydisperse Spray, Coalescence, Solid Propulsion, Alumina, Eulerian Multi-Fluid Model

The accurate simulation of polydisperse sprays encountering coalescence in unsteady gaseous flows is a crucial issue for solid rocket booster optimization. Indeed, the internal flow of the engine depends strongly on the alumina droplet size distribution, which spreads up with coalescence. Yet solving for unsteady two-phase flows with a high dimensional phase space is a challenge for both modelling and scientific computing. The usual Lagrangian approaches lead to a very high computational cost or to a low resolution level and they induce coupling difficulties to the Eulerian gaseous phase description. A wide range of Eulerian models have been recently developed to describe the disperse liquid phase at a lower cost and with an easier coupling to the carrier gaseous phase. Among these models, the Multi-Fluid model allows the detailed description of polydispersity and size/velocity correlations by separately solving fluids of size-sorted droplets, the so-called sections. On the one hand, the existing first order description of the size distribution in each section provides simple and fast resolution for coalescence. On the other hand, a second order method allows to reduce the number of sections required to capture accurately coalescence and to use elaborate droplet collision modelling, yet at the cost of heavier computation algorithms. This paper seeks to conclude on computational time and precision of both methods in order to choose the most efficient configuration for multi-dimensional unsteady rocket chamber simulations. Its objective is threefold: first, to validate the second order method by comparing simulations to reference solutions and dedicated experimental measurements conducted at ONERA, second to study the efficiency and robustness of both methods on coalescing size-conditioned dynamics problems, third, to draw some firm conclusions about the necessity to use first order or second order methods in order to capture the physics of solid propulsion configurations.

Introduction

Two-phase flows constituted of a gaseous phase and a disperse condensed phase play a key role in many industrial and scientific applications : fuel atomization and evaporation by Diesel injectors; fluidized beds; gas bubbles in oil or boiling water pipeworks; dynamics of planet formation in solar nebulae; etc. In all these applications the disperse phase is composed of particles of various sizes that can eventually coalesce or aggregate, break-up, evaporate and have their own inertia and sizeconditioned dynamics. So the importance of polydispersion is obvious for a full modelling of these phenomena.

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

Focusing on the dynamics of the alumina cloud in the combustion chamber is sufficient to evaluate specific impulse loss in the nozzle and slag material accumulation. Therefore the study takes place in the combustion chamber, where the principal physical processes that must be accounted for are : transport in real space, acceleration of droplets -due to drag conditioned by size-and coalescence, leading to polydispersity. These processes are remarkably sensitive to size distribution. We will therefore choose a model accurate as regards the size distribution and work with non-evaporating sprays throughout the paper, keeping in mind the broad application fields related to this study. By spray, we denote a disperse liquid phase constituted of droplets carried by a gaseous phase. We consider the specific case of dilute sprays i.e. where the liquid volume fraction is much smaller than one.

The retained approach called "mesoscopic" -or sometimes "kinetic" in reference to the kinetic theory of gases-describes the droplets as a cloud of point particles for which the exchanges of mass, momentum and heat are described using a statistical point of view, with eventual correlations : a finite set of global properties such as size of spherical droplets, velocity of the center of mass, temperature are modelled so that the total phase space is usually high-dimensional. More details about the droplets, such as angular momentum, can be predicted by increasing the size of the phase space : it is established that refined droplets models can be used as long as they do not include history terms. [START_REF] Williams | Spray combustion and atomization[END_REF][START_REF] Williams | Combustion Theory (Combustion Science and Engineering Series)[END_REF] proposed a transport equation based on kinetic theory that has proven to be useful for treating polydisperse, dilute and moderately dense liquid sprays.

Such an equation describes the evolution of the number density function or NDF of the spray due to the drag force of the gaseous phase and the droplet-droplet interaction of coalescence [START_REF] Kuentzmann | Aérothermochimie des suspensions[END_REF][START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF][START_REF] Laurent | 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 multifluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF].

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

As an alternative, the Eulerian Multi-Fluid model, furthered by [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 multifluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] from the ideas of [START_REF] Greenberg | On the origin of spray sectional conservation equations[END_REF], relies on the derivation of a semi-kinetic modelling from the Williams equation using a moment method for velocity, but keeping the continuous size distribution function [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]. This distribution function is then discretized using a "finite volume approach" in size phase space that yields conservation equations for mass, momentum (and eventually other properties such as number, energy) of droplets in fixed size intervals called sections, each of them constituting a different "fluid". Please note that integrating on a continuous size variable in each section is a key aspect : while other Eulerian approaches often consider discrete droplet sizes gathered into "classes" which cannot account for the new droplet sizes created by coalescence (except scarce examples such as in [START_REF] Vasenin | Two-phase flows in the nozzles of solid rocket motors[END_REF]), continuous size approaches such as the sectional method hereafter described are the only Eulerian methods handling coalescence naturally and rigorously. After integration on the sections, the resulting conservation equations are similar to those of the pressureless gas dynamics [START_REF] Bouchut | On zero pressure gas dynamics[END_REF] ; Zel ′ dovich (1970)) and lead to singu-lar behaviors such as delta shocks and vacuum zones. Well-suited numerical methods are thus required, some specific schemes being presented in De [START_REF] Chaisemartin | Modèles eulériens et simulation numérique de la dispersion turbulente de brouillards qui s'évaporent[END_REF]. The model finally requires closure relations for the phenomena accounted for by the Williams equation. We refer to [START_REF] Abramzon | Droplet vaporization model for spray combustion calculations[END_REF]; [START_REF] 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 multifluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] for detailed droplet models for which the Multi-Fluid model can be easily extended. These extensions do not have any impact on the conclusions of the present study.

The Eulerian Multi-Fluid model has proven its capability to simulate the size-conditioned dynamics of polydisperse sprays including coalescence with a first order resolution method of the size distribution in each section [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] even if the number of sections with the original method has to be large for accuracy purposes. In the context of evaporating sprays, several second order methods have been developed. A second order method using exponential size distributions in the sections is discussed in [START_REF] Laurent | Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF]. These two methods are in fact particular cases of a general method requiring n moments and introduced in Massot et al. (2010). The first order method can easily solve coalescence [START_REF] Laurent | Eulerian multifluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]) while the second order method has been extended to coalescence at the cost of some more elaborate algebra in [START_REF] Dufour | Modélisation multi-fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF]. It can yet easily include advanced models such as collision and coalescence efficiencies developed in [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].

For the purpose of this paper, the first and second order methods have been implemented in a research code solving dilute sprays in a pseudo 2D nozzle with one-way coupling to the gaseous phase. The physical features have yet been simplified : the drag force is modelled by a Stokes law; the temperature, composition, density and viscosity of the gaseous phase are assumed to be constant and uniform; as a consequence unsteady heating of the droplets is negligible so that their temperature and density are also constant. This configuration is used to compare the two methods solving the dynamics and size evolution of a lognormal distributed spray. Moreover it fully validates the second order method by providing detailed comparisons to lagrangian and Multi-Fluid reference solutions. The research code also allows to simulate the complete dynamics of an experiment on coalescence (D'Herbigny and Villedieu 2001) which had been initially designed to validate collision efficiency models with an averaged analytical formula. This work is thus a continuation of the experiment exploitation and provides precious validation for such two-phase flow models, which too often lack experimental back up. Moreover this test case confirms the robustness of the method towards sharp or monodisperse distributions. These validations globally show the compliance of the second order Multi-Fluid method to the features required by rocket booster simulations i.e. accuracy on polydispersity and dynamics, advanced coalescence models, robustness and fast computation. In addition, the second order Multi-Fluid model had previously been featured in CEDRE, a multiphysics 3D industrial code developed at ONERA, the French aerospace lab. CEDRE provides fully coupled aero-thermochemical resolution for energetics problems. It includes advanced models and can therefore compute drag force, evaporation and coalescence terms depending on all the gas parameters. The paper thus presents rocket booster simulations performed with an optimized version of the second order Multi-Fluid to achieve the proof of the efficiency of such models in a complex physical background.

The paper is organized as follows. Section one is dedicated to the derivation of the Multi-Fluid method as a theoretical framework : the origin and assumptions of the corresponding coalescence models are detailed, as well as the resulting coalescence formulas for first and second order methods. In section two, we evaluate the accuracy/speed compromise provided by the two Multi-Fluid methods implemented in the research code on a simple pseudo-2D validation test case used in [START_REF] Laurent | Eulerian multifluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. In section three, we validate the second order Multi-Fluid method and its collision and coalescence efficiency models on a 1D coalescence efficiency experiment conducted at ONERA and detailed in D' Herbigny and Villedieu (2001). In section four, and to illustrate the ability of the second order model to simulate in a reasonable time a complete solid propulsion test case, we give some numerical results in a modelled chamber and nozzle using the CEDRE code before concluding the study.

Mesoscopic Eulerian spray modelling : two Multi-Fluid methods

In this section, we introduce the framework of our study : the kinetic description of the disperse phase and the derivation of two Eulerian resolution methods. The separation between the two methods appears when presuming the size distributions on the sections. The origin and assumptions of the corresponding coalescence models are widely detailed.

A kinetic description : the Williams equation

Number density function When focusing on polydispersity, the size parameter ϕ of droplets is of capital importance but its natural expression depends on the ICMF 2010, Tampa, FL, May 30 -June 4, 2010 phenomena : volume v is relevant towards conservation of matter, surface S towards evaporation and radius r towards impact factor for instance. Since we assume spherical droplets, equivalence relation v = 4 3 πr 3 = 1 6 √ π S 3 2 allows the size to be expressed in this paper with the most comfortable notations.

Let us define the NDF function f ϕ of the spray, where f ϕ (t, x, u, ϕ)dxdϕdu denotes the averaged number of droplets (in a statistical sense), at time t, in a volume of size dx around a space location x, with a velocity in a du-neighborhood of u and with a size in a dϕneighborhood of ϕ. As for the size parameter conversions, we shall keep in mind that f r (r)dr = f S (S)dS = f v (v)dv and use the implicit notation f .

Williams equation

The evolution of the spray is described by the Williams transport equation [START_REF] Williams | Spray combustion and atomization[END_REF]. For a non evaporating case, it reads :

∂ t f + u • ∂ x f + ∂ u • (F f ) = Γ (1)
where F is the drag force per unit mass and Γ is the coalescence source term.

As an illustration, the Stokes law models drag force F due to the velocity difference with the gaseous phase and reads :

F(t, x, u, S) = U(t, x) -u τ p (S) , τ p (S) = ρ l S 18πµ g
where U is the gas velocity, µ g represents its dynamic viscosity and ρ l is the droplet (supposedly liquid) material density.

Coalescence operator

The kinetic modelling for the collision operator leading to coalescence is taken from [START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF]. We then assume :

[C1] We only take into account binary collisions (small volume fraction of the liquid phase)

[C2] The mean collision time is very small compared to the intercollision time

[C3] During coalescence, mass and momentum are preserved.

Thus, Γ = Q + -Q -where Q + and Q -respectively correspond to the quadratic integral operators associated with creation and destruction of droplets due to coalescence. These quadratic operators read [START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF][START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF]:

Q + = 1 2 ∫ ∫ u ⋆ ,v ⋆ ∈[0,v] f (t, x, u ⋄ (v, v ⋆ , u), v ⋄ (v, v ⋆ ))× f (t, x, u ⋆ , v ⋆ )B(|u ⋄ -u ⋆ |, v ⋄ , v ⋆ )Jdv ⋆ du ⋆ Q -= ∫ ∫ u ⋆ ,v ⋆ f (t, x, u, v)f (t, x, u ⋆ , v ⋆ )× B(|u -u ⋆ |, v, v ⋆ )dv ⋆ du ⋆ (2) where v ⋄ (v, v ⋆ ) = v-v ⋆ and u ⋄ = vu-v ⋆ u ⋆ v-v ⋆
are the precollisional parameters, J is the Jacobian of the transform (v, u) → (v ⋄ , u ⋄ ) : J = (v/v ⋄ ) d with d the dimension of the velocity phase space and B(|u -u ⋆ |, v, v ⋆ ) is the collision/coalescence probability kernel which reads :

B(|u-u ⋆ |, v, v ⋆ ) = E(|u-u ⋆ |, v, v ⋆ ) |u-u ⋆ | β(v, v ⋆ ) In this kernel, β(v, v ⋆ ) = π (r + r ⋆ )
2 is the impact parameter and E = E coll .E coal accounts for the collision efficiency E coll and the coalescence efficiency E coal . The modelling of these efficiencies is detailed in the following two paragraphs.

Collision efficiency

Collision efficiency E coll is a probability factor modelling the correlation of droplet velocities immediately before a collision. It is equivalent to consider that droplets may dodge each other due to the gas flow surrounding them. Collision efficiency laws thus require knowledge of local gas parameters such as density ρ g or viscosity µ g . In the case of unbalanced droplet sizes, D'Herbigny and Villedieu (2001) select two collision efficiency models : the Langmuir-Blodgett model, and the Beard-Grover model [START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF][START_REF] Achim | Simulation de collisions, coalescence et rupture de gouttes par une approche lagrangienne: application aux moteurs à propergol solide[END_REF]. These laws strongly depend on the bigger droplet Reynolds number and on a dimensionless number k which read, when taking r 1 the smaller radius and r 2 the bigger radius :

Re = 2ρgr2 µg | U -u 2 | , k = 2ρ l r 2 1 |u2-u1| 9µgr2
The number k is the ratio of the smaller droplet relaxation time τ r = 2ρ l r 2 1 9µg to its residence time in the bigger droplet influence zone τ h = r2 |u2-u1| .

In the case of low Reynolds numbers, Langmuir and Blodgett [START_REF] Langmuir | The production of rain by a chain reaction in cumulous clouds at temperatures above freezing[END_REF]) numerically get the following expression :

{ E 1 (k) = 0 if k ≤ 1.214 E 1 (k) = ( 1 + 3 ln(2k) 2(k-1.214)
) -2 otherwise whereas for high Reynolds, they get : ICMF 2010, Tampa, FL, May 30 -June 4, 2010 For intermediate cases they assume the following interpolation :

{ E 2 (k) = 0 if k ≤ 0.0833 E 2 (k) = k 2 (k+0.5) 2 otherwise 7 th International Conference on Multiphase Flow,
E LB coll (k, Re) = E 1 (k) 1 + Re/60 + (Re/60)E 2 (k) 1 + Re/60
(3) [START_REF] Beard | Numerical collision efficiencies for small raindrops colliding with micron size particles[END_REF] suggest to increase the accuracy of formula (3) which results from a simple interpolation between two limits. For this purpose, they use a numerical solution of the incompressible Navier-Stokes equations in order to determine the gaseous flow surrounding the bigger droplet depending on the Reynolds number. They can therefore evaluate precisely the forces on the smaller droplet and compute its trajectory. For Re ∈ [0, 400] and r 1 < r 2 , it yields:

E BG coll (k, Re) = 4 π 2 [arctan(max(H(k, Re), 0))] 2 (4)
where

H = 0.1465 + 1.302Z -0.607Z 2 + 0.293Z 3 Z = ln(k/k 0 ) k 0 = exp(-0.1007 -0.358 ln(Re) + 0.0261[ln(Re)] 2 )
Typical collision efficiency values with these two laws are in Fig. 1.

Coalescence efficiency

Coalescence efficiency E coal is the probability for two droplets to merge after a collision. Since the droplets can bounce on each other or separate by reflexion or stretching if the remaining kinetic energy of the new droplet is too high, coalescence efficiency depends a priori on the velocity parameters, the droplet material viscosity and surface tension, etc.

Most simple laws such as the Brazier-Smith model [START_REF] Achim | Simulation de collisions, coalescence et rupture de gouttes par une approche lagrangienne: application aux moteurs à propergol solide[END_REF][START_REF] Brazier-Smith | The interaction falling water drops : coalescence[END_REF]) are commonly used. Especially, the CEDRE code provides this model but it is not used in the following simulations. Therefore, the Brazier-Smith model is not detailed any further.

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

Semi-kinetic model

Velocity coherence assumption In a first step we reduce the size of the phase space, considering only the moments of order zero and one in the velocity variable at a given time, a given position and for a given droplet size : n = ∫ f du and ū = ∫ uf du/n which depends on (t, x, S). In order to close the system, the following assumptions are introduced :

[H1] For a given droplet size, at a given point (t, x), there is only one characteristic averaged velocity ū(t, x, S).

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

It is equivalent to presume the following NDF :

f (t, x, u, S) = n(t, x, S)δ(u -ū(t, x, S))
The monokinetic hypothesis [H1] is equivalent to reducing the velocity distribution support to a one dimensional submanifold parameterized by droplet size. It is fairly correct when the Stokes number is small enough for the droplets to follow the flow [START_REF] Massot | Eulerian Multi-fluid models for polydisperse evaporating sprays, chapter III of "Multi-Phase Reacting Flows: Modelling and Simulation[END_REF]). Yet it prevents droplet-crossing, which can occur in solid propulsion during ejection from fast vortices or high speed parietal injection.

The zero-dispersion hypothesis [H2] is justified when neither turbulence nor brownian phenomena interact with the disperse phase.

These hypotheses thus induce strong limitations for rocket booster simulation and will be discussed in the conclusion of this paper.

Equations This step leads to a system of conservation equations called the semi-kinetic model; it is a mono-kinetic model, saying that at a fixed position and time and for a fixed size, there is only one possible velocity. The semi-kinetic system reads :

∂ t n + ∂ x • (nū) = Q + n -Q - n ∂ t (nū) + ∂ x • (nū ⊗ ū) = n F + Q + u -Q - u (5)
Semi-kinetic coalescence operator In the semikinetic framework, the coalescence operator yields the evolution rate of the zeroth and first order moments of the velocity phase space. These two terms read, when omitting the (t, x) dependency :

Q + n = 1 2 ∫ v ⋆ ∈[0,v] n(v ⋄ (v, v ⋆ ))n(v ⋆ )β(v ⋄ , v ⋆ )I + n dv ⋆ Q - n = -n(v) ∫ v ⋆ ∈[0,+∞[ n(v ⋆ )β(v, v ⋆ )I - n dv ⋆ Q + u = 1 2 ∫ v ⋆ ∈[0,v] n(v ⋄ (v, v ⋆ ))n(v ⋆ )β(v ⋄ , v ⋆ )I + u dv ⋆ Q - u = -n(v) ∫ v ⋆ ∈[0,+∞[ n(v ⋆ )β(v, v ⋆ )I - u dv ⋆ (6
) where I - n , I + n , I - u and I + u are the partial collisional integrals. To preserve the monokinetic and zero-dispersion assumptions, these integrals are computed with average velocities and therefore read :

I + n =|ū(v ⋆ ) -ū(v ⋄ )| E(|ū(v ⋆ ) -ū(v ⋄ )|, v, v ⋆ ) I - n =|ū(v) -ū(v ⋆ )| E(|ū(v) -ū(v ⋆ )|, v, v ⋆ ) I + u = (v ⋄ ū(v ⋄ ) + v ⋆ ū(v ⋆ )) (v ⋆ + v ⋄ ) × |ū(v ⋆ ) -ū(v ⋄ )| E(|ū(v ⋆ ) -ū(v ⋄ )|, v, v ⋆ ) I - u = ū(v) |ū(v) -ū(v ⋆ )| E(|ū(v) -ū(v ⋆ )|, v, v ⋆ )

Multi-Fluid model

In a second step we choose a discretization 0 = S 0 < S 1 < • • • < S N for the droplet size phase space and we average the system of conservation laws over each fixed size interval [S k-1 , S k [, called section. The set of droplets in one section can be seen as a "fluid" for which conservation equation are written, the sections exchanging mass and momentum. In order to close the system, the following assumptions are introduced :

[H3] In one section, the characteristic averaged velocity does not depend on the size of the droplets.

[H4] The form of n as a function of S is supposed to be independent of t and x in a given section, thus decoupling the evolution of the mass concentration of droplets in a section from the repartition in terms of sizes.

These assumptions are equivalent to presuming the NDF in velocity and size inside each section k :

∀S ∈ [S k-1 , S k [ { ū(t, x, S) = ūk (t, x) n(t, x, S) = κ k (t, x, S)
Since we have chosen a constant value for the velocity distribution in each section, Multi-Fluid methods based on assumption [H3] feature a first order velocity convergence with the number of sections. It induces a convenient integration simplification but it can generate trajectory modes. For coarse size discretizations indeed, it assumes the same velocity for droplets of remarkably different sizes in a given section.

As for [H4], it allows to reduce the size distribution information of each section at (t, x) to a set of moments of S, the number of which depends on the choice of the (κ k ) k set of functions. The difference between the two methods we are about to introduce is yet the choice of their form.

In the so-called first order method, a constant 1 κ k function in each section yields a first order granulometry convergence with the number of sections. Whereas in the second order method, a two-coefficient 2 κ k function yields a second order granulometry convergence. The choice of refining the size distribution description is of course related to the need of a fine resolution of polydispersity, as told in introduction.

First order method

First order size assumption The first order Multi-Fluid method assumes, instead of [H4], the following notation in each section k :

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

∫ S k S k-1 1 κ k (S) ρ l 6 √ π S 3/2 dS = 1.
Such an approach only focuses on one moment of the distribution in the size variable : the moment in terms of mass is chosen because it is conserved by coalescence.

Please note that the distribution on the last section is a decreasing exponential with a fixed coefficient. This choice allows the final section to treat the bigger droplets but requires not to have a significant part of the mass. This is a major limitation compared to the second order method.

Equations

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

∂ t m k +∂ x •(m k ūk ) = 1 C m+ k -1 C m- k ∂ t (m k ūk )+∂ x •(m k ūk ⊗ ūk ) = m k Fk + 1 C mu+ k -1 C mu- k (7) 0 v ⋆ v ⋄ v i v i+1 v j v j+1
k s t r i p l o w e r l i m i t k s t r i p u p p e r l i m i t "(i, j, k) domain"

Figure 2: Coalescence terms are integrated on the "(i, j, k) domains"

where Fk is the "classical" averaged drag force per unit mass on a section [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] :

Fk = U -ūk τ k , 1 τ k = ∫ S k S k-1 ρ l 6 √ π 1 κ k (S)S 3/2 τ (S)
dS.

The total system thus counts twice as much equations as the number of sections.

First order coalescence terms

In each section equation, the creation coalescence terms result from a double integration : on the whole colliding partner size space at the kinetic level (eq. 2) and on the concerned section at the Multi-Fluid level. Yet the second dependency will not coincide with the section after mapping the natural variables of the two precursor colliding partners. Thus splitting the integration domain thanks to the section continuity yields elementary integrals Q ijk , triply indexed with the two precursor section numbers (i, j) and the destination section number (k). These domains are illustrated in Fig. 2.

In the particular case of the first order method, considering [H2] and assuming E coll = E coal = 1 [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF], the coalescence integrals Q ijk take the following form after factorizing the mass moments m i m j :

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

1 C m+ k , 1 C m- k , 1 C mu+ k and 1 C mu- k read : 1 C m+ k = I(k) ∑ i=1 m i N ∑ j=1 m j |ū i -ūk | ( 1 Q ⋆ ijk + 1 Q ⋄ ijk ) 1 C m- k = m k N ∑ j=1 m j |ū k -ūj | N ∑ i=1 1 Q ⋆ kji 1 C mu+ k = I(k) ∑ i=1 m i N ∑ j=1 m j |ū j -ūk | (ū i 1 Q ⋆ ijk + ūj 1 Q ⋄ ijk ) 1 C mu- k = m k ūk N ∑ j=1 m j |ū j -ūk | N ∑ i=1 1 Q ⋆ kji (8)
The integrands of the Q ijk integrals depend only on size parameters r,r ⋆ and r ⋄ thanks to [H4]. This allows the Q ijk integrals to be pre-calculated as soon as the section limits and the ( 1 κ k ) k are given, i.e. once and for all at the beginning of a simulation.

Second order method

Second order size assumption The second order Multi-Fluid model is based on a two-coefficient exponential approximation of the size distribution in each section. This means that [H4] reads, for

S ∈ [S k-1 , S k [ : 2 κ k (t, x, S) = a k (x, t) exp(-b k (x, t)S)
where

(a k (x, t), b k (x, t)) k ensures          ∫ S k S k-1 2 κ k (t, x, S)dS = n k (t, x) ∫ S k S k-1 2 κ k (t, x, S) ρ l 6 √ π S 3/2 dS = m k (t, x)
The choice of an exponential function ensures the positivity of the distribution function whatever scheme is chosen. It also aims at reducing the number of sections and is well suited for evaporation, which requires mass flux information at the section boundary. On the other hand, coalescence becomes more difficult to compute if one wants to respect realizability conditions on the sections i.e. the fact that (m k , n k ) k couples are conditioned by the section boundaries since they drive information on average droplet volume .

Equations

The conservation equations for the k th section now read :

∂ t n k +∂ x •(n k ūk ) = 2 C n+ k -2 C n- k ∂ t m k +∂ x •(m k ūk ) = 2 C m+ k -2 C m- k ∂ t (m k ūk )+∂ x •(m k ūk ⊗ ūk ) = m k Fk + 2 C mu+ k -2 C mu- k
(9) The total system thus counts three times as much equations as the number of sections.

Second order coalescence terms Since the time and space dependency of the size-distribution functions 2 κ k (t, x, S) is no longer factorizable as was m k in the first order method, the Q ijk integrals must be computed at each time step in each cell on its "(i, j, k) domain". So it is numerically interesting to integrate the collision/coalescence efficiency at the same time. That is why this method is "less inappropriate" for efficiency models. Let us take the notation :

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

2 Q n ijk = ∫ ∫ v ⋆ +v ⋄ ∈[v k ,v k+1 [ Ψ ij (t, x, r ⋆ , r ⋄ , |u ⋄ -u ⋆ |)dr ⋆ dr ⋄ 2 Q ⋆ ijk = ∫ ∫ v ⋆ +v ⋄ ∈[v k ,v k+1 [ Ψ ij (t, x, r ⋆ , r ⋄ , |u ⋄ -u ⋆ |) 4 3 πρ l r ⋆3 dr ⋆ dr ⋄ 2 Q ⋄ ijk = ∫ ∫ v ⋆ +v ⋄ ∈[v k ,v k+1 [ Ψ ij (t, x, r ⋆ , r ⋄ , |u ⋄ -u ⋆ |) 4 3 πρ l r ⋄3 dr ⋆ dr ⋄ The coalescence source terms 2 C n+ k , 2 C n- k , 2 C m+ k , 2 C m- k , 2 C mu+ k and 2 C mu- k
are still written as direct sums of the "(i, j, k) integrals" but therefore read :

2 C n+ k = I(k) ∑ i=1 N ∑ j=1 2 . 2 Q n ijk , 2 C n- k = N ∑ j=1 N ∑ i=1 2 Q n kji 2 C m+ k = I(k) ∑ i=1 N ∑ j=1 ( 2 Q ⋆ ijk + 2 Q ⋄ ijk ) , 2 C m- k = N ∑ j=1 N ∑ i=1 ( 2 Q ⋆ kji + 2 Q ⋄ kji ) 2 C mu+ k = i ∑ j=1 KH ∑ k=KL (ū j 2 Q ⋆ ijk + ūk 2 Q ⋄ ijk ) , 2 C mu- k =ū i . 2 C m- k (10)

Numerical strategy for coalescence

Towards computation, the main difference with the first order method is the time and space dependency of the in-

tegrals 2 Q n ijk , 2 Q ⋆ ijk and 2 Q ⋄ ijk .
The second order method requires indeed a costly numerical integration for each time step and space cell. We have chosen to perform these integrations with a 2D 5 point Newton-Cotes quadrature which requires

Ψ ij (t, x, r ⋆ , r ⋄ , |u ⋄ -u ⋆ |) to be computed 5 2 × N 2
sec × N cell times per time step. The accuracy of this method on such exponential functions has been tested but will not be discussed any further.

Yet performing this integration with an advanced collision or coalescence efficiency model is now natural and hardly more costly. That's why these efficiencies have only been included in the second order method. [START_REF] Laurent | Eulerian multifluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] introduces a general form for coalescence efficiency formulas to be suitable for a pre-calculated first order method but no collision efficiency models can be easily implemented since they require local gas parameters. That is why the first order method here excludes such collision efficiency models.

Dynamic Study : the Nozzle Test Case

As a main criterion, we want to compare the two Multi-Fluid methods on their ability to describe the dynamics of a coalescing cloud. We therefore need a wellsuited test case, inducing coalescence as well as sizeconditioned dynamics and difficult enough to highlight the limitations of the methods.

For that purpose, we have chosen the configuration which was used 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] to validate the first order method compared to a reference Lagrangian solver : a 2D axisymmetrical conical decelerating nozzle, designed in such a way that it admits, for one-way coupling spray dynamics a self-similar solution. We though inject a lognormal distribution instead of a bi-modal one.

Definition of configuration

The chosen configuration is stationary 2D axisymmetrical in space and 1D in droplet size. It is described in detail in [START_REF] Laurent | Eulerian multifluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. Hence, only its essential characteristics are given here.

For the problem to be 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 U z and the reduced radial velocity U r /r, for z ≥ z 0 :

U z = U (z) = z 2 0 U (z 0 ) z 2 , U r r = U z z = z 2 0 U (z 0 ) z 3
where z 0 > 0 is the coordinate of the nozzle entrance and the axial velocity U (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 case of coalescence, it is valid in the neighborhood of the centerline.

Because of the deceleration of the gas flow in the conical nozzle, droplets are slowing down too, 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 ) = -2U (z 0 )/z 0 ; it is chosen large enough so that the velocity difference developed by the various sizes of droplets is important. [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] chose rather large values, as well as strong deceleration, leading to extreme cases: U (z 0 ) = 5 m.s -1 , z 0 = 5 cm : these values generate a very strong coupling between coalescence and droplet dynamics. These severe conditions make the test cases under consideration very efficient tools for the numerical evaluation of the two Eulerian models. Finally please note that no efficiency law is used in this section so that the second order method remains comparable to the first order method, which lacks this sort of model. We thus assume E coll = E coal = 1 in this section.

Droplet initial distribution

The lognormal distribution is useful to characterize alumina particles in solid propergol rocket boosters. We therefore choose a lognormal distribution on the surface variable, without loss of generality. The droplets are constituted of liquid alumina, their initial velocity is the one of the gas, their initial temperature, fixed at the equilibrium temperature 3600 K (corresponding to an infinite conductivity model), does not change along the trajectories.

For the sake of the first order method accuracy, we do not want to transfer too much mass in the last section. The initial injected mass concentration is therefore taken as m 0 = 1.0609 kg.m -3 . For the same reason, the lognormal parameters are set to S LN = 1600µm 2 and δ LN = 1.5 which corresponds to a sharp distribution centered on a radius of 11.3µm. We shall yet keep in mind that the second order method can accurately treat cases where the final section hosts a significant portion of the mass whereas the first order method would need to increase the number of sections. The nozzle test case has thus been designed to meet optimal first order specifications i.e. where as much sections as possible get filled but the final section remains almost empty.

We now choose the tested numbers of sections : a 5 section test case illustrates what happens with a coarse discretization. A 13 section and a 25 section test cases prove the convergence of the method. They are compared to a 53 first order run which we use as a reference solution. The first order indeed has been fully validated compared to a lagrangian test case 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].

Distribution spread up results

The processing of polydispersion gives a first indication on the method's accuracy. We therefore compare the mass concentration distributions at the nozzle output to the 53 section reference test case in Fig. 3. With 5 sections, the first order strongly overestimates the size growth while the second order underestimates it. With 13 sections, the trend is the same but the error is smaller. Finally for 25 sections, we consider both methods are roughly converged. The growth overestimation of the first order method for poor size discretization brings about major consequences for the spray dynamics.

Spray dynamics results

Let us now appreciate the consequences of polydispersion treatment on dynamics. Fig. 4 shows the evolution of mass and number concentrations along the centerline of the nozzle computed with both Eulerian methods (5,13,25 sections) compared to the reference solutions (53 section first order method and lagrangian reference test case). For the first order case, number concentrations are computed considering section average droplet volumes given by integration of the ( 1 κ k ) k functions : data are therefore redundant with mass concentration data. We emphasize the fact that the convergence of the Table 1: Composition of the five section groupings two methods is achieved for 25 section simulations as we can see in Fig. 5, which is a zoom of Fig. 4.

To compare precisely the effect of polydispersion on dynamics, let us now consider five size intervals : they correspond to the sections in the 5 section case and are composed of section groupings in the other cases as illustrated in table 1. The evolution of the mass concentration of these groupings along the nozzle is given in Fig. 6. It is there obvious that the 5 section first order error on coalescence is severe, especially in the fifth and last grouping G 5 where no mass should be found consequently to the problem parametrization. Finally, the second order method reveals to have acceptable error with as few as 5 sections.

Computational time

To conclude this study, table 2 recaps the duration of the different runs. It is obvious that the first order method is much faster since it is highly optimized thanks to the pre-calculation. Even when considering the 5 section second order run -which is fairly as accurate as the 13 section first order-it is 6 time slower. Yet two restrictions must be added about the use of the first order with a thinner size space discretization. First, it is not possible to meet such a high computation speed with specific collision/coalescence laws. Second for typical solid propulsion simulations, each section requires the resolution of an unsteady 3D pressureless Euler-type system so that solving for coalescence is no longer the limiting step. Thus increasing the number of section becomes a costly operation, restraining the profit of the first order method.

Coalescence Experimentation : the D'Herbigny Test Case

In this section, we provide further validation for the second order Multi-Fluid method, by comparing it to a 1D experiment on collision/coalescence efficiency models. Moreover this configuration provides tough testing for Eulerian models because it deals with bimodal size dis-

The D'Herbigny collision efficiency experiment

In the D'Herbigny experiment, coalescence is studied through the growth of a bigger droplet falling through a fog of smaller droplets. Details about the experimental device and conditions can be found in D' Herbigny and Villedieu (2001) and in Fig. 7. In this experiment, the average injection radius of the bigger droplets is r 2 = 150µm thanks to piezoelectric injector which is fairly accurate. The smaller droplets, composing the fog, have a radius r 1 ∈ [2µm, 4µm]. The two types of droplets have an approximately constant velocity difference of 3m/s along the tunnel. One can therefore compute the values of the k parameter used in the Langmuir-Blodgett and Beard-Grover collision efficiency models. They are given in table 3.

Moreover, the collisional Reynolds number Re is around 60 which is comparable to the values in the Ariane 5 P230 booster. This range of values (moderate Re and low k) is supposed to be favorable to the Beard-Grover law.

Results

Droplet radius growth study D'Herbigny and Villedieu ( 2001) provide an analytic formula to estimate the Sauter mean radius growth ∆r 2 of the bigger droplets when assuming a small growth and a constant velocity difference :

∆r 2 = αE coll C v (11)
where C v is the volume fraction of the fog, E coll the average collision efficiency and α results from the integration of the impact parameter. As a first step, we thus perform the simulation with a constant velocity of 3m/s to fit the analytical result. Experimental, analytical and simulation results are given in Fig. 8. The first conclusion is that collision efficiency laws have a strong impact on radius growth, here undermining it up to 50%.

The second conclusion is a positive validation of both laws implemented in the second order Multi-Fluid Method. According to D' Herbigny and Villedieu (2001), this constant velocity study also concludes that Langmuir-Blodgett law is most suited for rocket chamber-like conditions.

Method numerical stability

The third conclusion of these simulations is that the method allows robust and accurate simulation of stiff distributions. Indeed the simulations are converged for a rough number of 100 iterations and provide a self-similar solution which means low diffusion in the size phase space.

Discussion on velocity

As a second step, let us revise the constant velocity assumption. On the one hand, droplets do reach a terminal velocity, i.e. the drag force compounds the weight but this limit velocity increases with the droplet size since the drag force increases like the droplet surface while the weight increases like the droplet volume. On the other hand, collisions with a static fog induce a momentum dilution i.e. a significant slow down. We can see the evolution of the effective velocity as calculated by the code when gravity, drag force and coalescence momentum transfer are enabled in Fig. 9. The droplet velocity tends to be slightly smaller than the terminal velocity because of momentum dilution. Yet the consequences of velocity on coalescence are not obvious. First let us neglect efficiency laws. The coalescence rate per unit time linearly increases with the velocity difference (eq. 2) but linearly decreases with the time spent in the tunnel so that it vanishes in the final growth formula (eq. 11). In fact radius growth for a fixed length travel only depends on the number of collisions i.e. fog concentration and travel length, assuming a fairly constant cross section. However when considering efficiency laws, the velocity difference plays a role since it controls the way smaller droplets can dodge the bigger ones. For instance eq. 11 depends linearly on an averaged collision efficiency factor strongly conditioned by velocity difference as shown in Fig. 1.

The new simulation results with full dynamics is given in Fig. 10. It shows that the Beard Grover efficiency correlation is not so bad while the average velocity might have been slightly overestimated in D' Herbigny and Villedieu (2001). As a final restriction, please note that no coalescence efficiency law has been used (E coal = 1) which suggests not to conclude on which collision efficiency law to use.

Conclusion on the model and methods

The previous results confirm that Eulerian Multi-Fluid methods can be used to simulate accurately the sizedistribution evolution of a coalescing spray and its sizeconditioned dynamics. The first order method provides, with a reasonable number of sections (between 10 and 20), extremely fast results since it is highly optimized.

The second order Multi-Fluid method provides good results with a very coarse size-space discretization (as few as 5 sections). If both methods can include coalescence efficiency models, only the second order method easily provides collision efficiency models. The corresponding laws can significantly reduce the particle growth rate which is experimentally observed so that these efficiency models are essential to capture the physics of coalescence.

7 th International Conference on Multiphase Flow, ICMF 2010, Tampa, FL, May 30 -June 4, 2010

Rocket booster simulation

In this section, we conclude the study by achieving a 2D complex simulation of a rocket booster chamber and nozzle which experiences the meeting of a parietal injection of gas and small alumina droplets resulting from the combustion of propergol and an axial two-phase flow carrying bigger droplets supposedly coming from priorly injected, burnt and coalesced particles. The upper part of the chamber generating these bigger droplets is not solved.

The CEDRE code

The CEDRE code has been developed at ONERA for several years to provide a multiphysics comprehensive tool for energetics problems in aeronautics. The main feature is the combination of a non structured spatial solver with models as varied as multiphase flows, multispecies chemistry, thermal conduction, radiation, wallfilm models, etc. The coupling can be one or twoway. As for the disperse phase, CEDRE includes a Lagrangian and two Eulerian methods i.e. a multi-class and the second order Multi-Fluid method studied here. Only the Lagrangian and the Multi-Fluid methods provide coalescence models.

Configuration

The 2D rocket booster simulation is performed on a simplified configuration, yet featuring the main difficulties of solid propulsion typical flows i.e. parietal injection and a supersonic nozzle. It takes place on a deformedstructured 1500 cell mesh (Fig. 11,top). When seeking a two-phase stationary flow solution, one usually takes a converged gas flow field as an initial state. With a total rate of flow of 10 kg.s -1 .m -2 from the wall (propellant combustion) and the head end (upstream flow), this gas flow involves extemely high velocity gradients in the nozzle (Fig. 11, middle), which will induce sizeconditioned droplet dynamics as studied in section two. The simulation strictly speaking starts when injecting the disperse phase : a monomodal 5µm radius wall injection represents the droplets resulting from recently burnt aluminum particles, directly expelled from the propellant grain, and another roughly monomodal distribution around 20µm is injected on the axis to model the previously coalesced droplets [START_REF] Fabignon | Combustion of aluminum particles in solid rocket motors[END_REF].

For the purpose of our simulation, the injected volume fractions (more or less 5.10 -5 ) approximate correctly typical rocket booster conditions and preserve the dilute spray hypothesis.

Results

For the disperse phase, we choose a 5 section discretization (0, 12.5, 25, 37.5 and 50 µm) as in section two. A 10 -2 s time interval is required to allow the first droplets to reach the end of the nozzle. With a 10 -6 s time step (10, 000 iterations), we perform a 1 h single processor AIX platform computation. The stationary volume fractions for the five sections are displayed in Fig. 12.

First, the disperse phase has an impact on the gas flow since we perform a two-way coupling simulation. Specific impulse loss can indeed be observed with the nozzle Mach number decrease in Fig. 11, bottom. Second, we do note bigger particle creation in Fig. 12 : coalescence occurs as soon as the two injected types of droplets meet since the third section gets filled with significant volume fractions (≈ 10 -5 ) in the chamber. Third, bigger droplets are created in the nozzle eventhough acceleration there induces a strong volume fraction dilution. In pursuance of section two, section 5 is left almost empty with a volume fraction around 10 -10 (the colormap is different from the other section ones); yet a different discretization could be used to activate this section thus increasing size accuracy.

Remaining Numerical difficulties in solid propulsion

The final simulation of this paper illustrates the ability of Multi-Fluid models to provide in a reasonable time fairly accurate information on polydisperse coalescing sprays in complex geometry gas flows. The above results however suggest several possible improvements, among which reducing computational time and considering more physical phenomena. Typical advanced rocket booster simulations are unsteady and performed on more complex meshes : they can include from 10, 000 cells to millions of cells for 3D meshes. Highly optimized code and multiprocessor computations (with several spatial domains for instance) are therefore required. In our constant time step simulation, the CFL condition is driven by the flow in the nozzle but the gap between the chamber and the nozzle stay times makes most of the early time steps too short and therefore uselessly numerous. An adaptative time step method considering both phases including operator splitting could prevent from ill-suited time steps.

Physically, high velocity gradients in the nozzle induce droplet secondary break-up which strongly undermines droplet growth. Since it is out of the scope of this paper, this phenomenon has been neglected but its modelling reaches two requirements : specific impulse computation is subject to the precise knowledge of droplet size after the nozzle throat and experimental data on alumina size distribution mainly consist of ejected material which has encountered complete break-up process. Developpement of fragmentation models for the Eulerian Multi-Fluid method has been achieved in [START_REF] Dufour | Modélisation multi-fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF].

Finally, the main drawback of such methods comes from the monokinetic hypothesis : it forbids droplet crossings which should occur on the centerline of our simulation and induces instead a droplet accumulation artifact due to the momentum averaging on the symmetry axis. On more advanced solid propulsion unsteady configurations, dedicated to hydrodynamic instability studies for instance, high-inertia droplets ejected from vortices need a multivelocity treatment which advanced high order moment methods can provide as developed in [START_REF] Kah | Eulerian quadraturebased moment models for dilute polydisperse evaporating sprays[END_REF].

Conclusions

In this paper, we provide a comprehensive validation of a Eulerian, second order in size, method for solving twophase polydisperse coalescing flows. The conclusion of this study is that the existing method is able to simulate accurately the dynamics of such flows with less sections than the first order method as demonstrated in section two yet in a slower way. It is however the only Eulerian Multi-Fluid method including validated advanced collision efficiency models, which are crucial for rocket booster studies. Finally, the method is practically efficient in rocket booster simulation contexts; moreover reducing the number of required sections can be profitable to the computational time of such simulations. The remaining limitations underlined in the final section require though further developments but the method can bear such improvements.

Figure 1 :

 1 Figure 1: Two collision laws selected for booster applications with r 2 = 150µm (red : r 1 = 2µm, green : r 1 = 3µm; blue : r 1 = 4µm)

Figure 3 :

 3 Figure 3: Mass concentration distribution at the nozzle's end computed with the first and second order methods (empty and filled symbols) and reference (black)

Figure 4 :

 4 Figure 4: Number and mass concentration along nozzle for 5, 13, 25 sections and reference solution (black)

Figure 6 :

 6 Figure 6: Evolution of mass concentration in the four finite section groupings and in the last section for 5, 13 and 25 sections

  Figure 7: D'Herbigny experimental device

Figure 8 :

 8 Figure 8: Radius growth ∆r 2 after 5m through a C v concentration fog (empty red : Analytical; filled green : Second order simulation with constant velocity; black : Experimental).

Figure 9 :

 9 Figure 9: Velocity along the tunnel with gravity, drag force and coalescence momentum transfer (blue triangles) and terminal velocity for the corresponding Sauter Mean Radius (black)

Figure 10 :

 10 Figure 10: Radius growth ∆r 2 with velocity resolution (filled blue : Second order simulation; black : Experimental)

Figure 11 :

 11 Figure 11: Top : Deformed-structured 1500 cell mesh (Arrows : injection zones); Middle : Gaseous Mach number without particle; Bottom : with particles

Figure 12 :

 12 Figure 12: Dispersed phase volume fractions

Table 2 :

 2 Computational time for 10, 000 iterations on a 2.66GHz Intel Core 2 Duo CPU

		1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06
	Mass concentration (kg/m 3 )	1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058 1.044 1.046 1.048 1.05 1.052 1.054 1.056 1.058
		1.042 1.042 1.042 1.042 1.042 1.042 1.042 1.042 1.042 1.042 1.042 1.042 1.042 1.042 1.042
		0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16	0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18	0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2	0.22 0.22 0.22 0.22 0.22 0.22 0.22 0.22 0.22 0.22 0.22 0.22 0.22 0.22 0.22	0.24 0.24 0.24 0.24 0.24 0.24 0.24 0.24 0.24 0.24 0.24 0.24 0.24 0.24 0.24
			Position (m)
	Figure 5: Zoom on the convergence of mass concen-
		tration at the nozzle's end (black : rough la-
		grangian validation)
			1st order MF 2nd order MF
		5 sections	1 s		10 s
		13 sections	1.5 s		50 s
		25 sections	2 s		180 s
		53 sections	5 s		n.a.

Table 3 :

 3 Collision parameter k values in the D'Herbigny experiment

th International Conference on Multiphase Flow, ICMF 2010, Tampa, FL, May 30 -June 4, 2010

Acknowledgments

The present research was done thanks to a Ph. D. Grant from DGA, Ministry of Defence (M. S. Amiet, Technical Monitor). The authors also thank G. Dufour and T. Fontfreyde for their preliminary contribution to the second order method in the CEDRE code.


Nicolas Le Sommer
email: nicolas.le-sommer@univ-ubs.fr

Salma Ben

Sassi Valoria

Location-based Service Discovery and Delivery in Opportunistic Networks

Keywords: Service Provision, Opportunistic Networking, Disconnected Mobile Ad Hoc Networks

Opportunistic networks are usually formed spontaneously by mobile devices equipped with short range wireless communication interfaces. Designing and implementing a routing protocol to support both service discovery and delivery in such kinds of networks is a challenging problem on account of frequent disconnections and topology changes. In these networks one of the most important issues relies on the selection of the best intermediate node(s) to forward the messages towards their destination(s).

In this paper, we propose a new location-based opportunistic routing protocol enabling a bandwidth-efficient service discovery and delivery in a wide area network composed of numerous mobile devices. This protocol implements self-pruning heuristics allowing mobile devices to decide whether they efficiently contribute in the delivery of the messages they receive from their neighbors. This protocol was validated through simulations, which proved its efficiency and cost-effectiveness.

I. INTRODUCTION

Mobile devices equipped with short range wireless communication interfaces (e.g., Wi-Fi, Bluetooth), and sometimes with a GPS receiver, are nowadays widespread and used daily by an increasing number of people. Netbooks, mobile internet devices or smartphones are some examples of such devices. Thanks to their short range wireless communication interfaces, these devices can spontaneously form a multi-hop disconnected mobile ad hoc network (DMANET). Designing a routing protocol supporting both service discovery and delivery in such kinds of networks is radically different than devising one for traditional infrastructure-based networks. Indeed due to the mobility of nodes and to the short communication range of wireless interfaces, the topology of DMANETs suffers from frequent and unpredictable changes, entailing an intermittent connectivity between nodes. In DMANETs, devices can thus communicate directly only when they are in range of one another. Otherwise, intermediate nodes can be used to relay a message from a source to a destination following the "store, carry and forward" principle. Since in DMANETs we do not have any knowledge about the network evolution, routes are computed dynamically at each hop while forwarding messages towards their destination(s). Each node receiving a message intended to another recipient, is expected to exploit its local knowledge in order to decide for instance which is the best next hop among its current neighbors to deliver the message. When no forwarding opportunity exists (e.g., no other nodes are in the transmission range, or the neighbors are evaluated not suitable to that communication) the node stores the message and waits for future contact opportunities with other devices to forward the message. Thanks to this principle, a message can be delivered even if the sender and the recipient are not present simultaneously in the network, or if they are not in the same network partition at the emission time. For instance in the disconnected mobile ad hoc network illustrated in Figure 1, Tom may discover and invoke the services offered by infostation I, even if he is not in the communication range of I. Indeed, the messages sent by Tom and by infostation I may be forwarded opportunistically by intermediate devices, such as those carried by Bob and Alice who are walking along the main street of the city.

In this paper, we propose a new message opportunistic forwarding protocol, called OLFServ, suited to the provision of location-aware software services in disconnected, partially connected or intermittently connected MANETs. This protocol is designed to perform an efficient and geographicallycontrolled broadcast of service advertisements and service discovery requests. This forwarding protocol is also suited to the invocation of remote service providers. Indeed, it enables service providers to specify in their advertisements their location and the geographical area where they can be discovered and invoked. It allows clients to define their location, the area where their messages must be disseminated, the location of the provider they want to invoke, as well as their moving direction and speed if they know them. Based on these pieces of information, intermediate nodes running OLFServ can decide themselves if they are "good" relays to deliver the messages they receive from their neighbors (i.e., if they contribute to bring a message closer to its destination). For this purpose, OLFServ implements self-pruning techniques in order to select a small subset of intermediate and relevant nodes that will have to broadcast the messages they receive in given geographical areas. In OLFServ, intermediate nodes can update location properties specified by clients and providers in service messages while forwarding these messages in order to refine the area where they must be disseminated progressively, until reaching the destination. In the current scenarios we focus on, service providers are fixed infostations deployed in a city, while clients are devices carried by humans (see Figure 1). Consequently, a client may not stay long enough at the location where he was when he sent its service invocation to receive at this same place the response returned by the provider it invoked. In order to address this issue, OLFServ enables service providers to compute, based on the properties exhibit by a client (location, speed, direction/destination), an "expected area" where the client is supposed to be when he will receive its response. Finally our protocol implements a random delayed message forwarding in order to avoid message collisions.

The remainder of the paper is organized as follows. Section II motivates our work by discussing some recent research works dealing with routing protocols in disconnected, partially connected or intermittently connected mobile ad hoc networks. Section III presents our new opportunistic and location-aware forwarding protocol. Section IV presents the results we obtained by running our protocol on a network simulator, and we compare these results with other solutions. Section V concludes this paper and gives some perspectives we plan to work along in the future.

II. RELATED WORK

Service provision relies on two main processes: service discovery and service invocation. In decentralized and highly changing environments such as those we consider, the service discovery process cannot rely on a centralized approach since no device is stable enough, or accessible permanently, to act as a service register. Therefore, each client is responsible for maintaining its own perception of the services offered in the network, and for discovering these services. The discovery process can be performed reactively, by processing the unsolicited service advertisements broadcast by service providers, or proactively, by broadcasting service discovery requests in the network and by processing the advertisements returned by providers in response. The service invocation, during which a given client actually interacts with the provider of a previously discovered service, is generally performed using a unicast and destination-based communication model. A selection process may precede the invocation, when the opportunity is given to the service client to choose among several providers.

Forwarding protocols designed for opportunistic networking and/or delay tolerant networking can be classified according to both the network knowledge required by mobile hosts in order to run them and the methods used to select the hosts that must forward the messages. Three kinds of methods can be identified: blind-based methods, capability-assessment-based methods, and neighbor-knowledge-based methods. The last category can be divided in two sub-classes: the designating methods and the self-pruning methods. In the first sub-class, the message forwarder (or the initial message sender) is responsible for selecting the next suitable intermediate node(s), whereas, in the second sub-class, intermediate nodes receiving a message from a neighbor node are expected to decide whether they must forward this message or not based on their local knowledge.

Blind-based protocols, such as the Epidemic Routing protocol [START_REF] Vahdat | Epidemic Routing for Partially Connected Ad Hoc Networks[END_REF], do not require any knowledge and do not make any attempt to reduce the number of forwarding nodes. The messages are blindly stored and forwarded to all neighboring nodes, generating a flood of messages. The drawback of such an epidemic dissemination lies in the very high number of message copies which are needed to obtain a successful delivery of messages. These protocols are thus not suitable to environments with high density regions, since they would generate too much network traffic and sometimes could even lead to network congestion. This problem is referred to as the broadcast storm problem [START_REF] Ni | The Broadcast Storm Problem in a Mobile Ad Hoc Network[END_REF].

This drawback is addressed by protocols implementing methods allowing to assess the capability of neighbor nodes to contribute in the message delivery. These methods traditionally use a probabilistic metric, often called delivery predictability, that reflects how a neighbor node will be able to deliver a message to its final recipient. Before forwarding (or sending) a message, a mobile host asks its neighbors to infer their own delivery probability for the considered message, and then compares these probabilities and selects the best next hop(s) among them. This computation can require a 1-hop, and sometimes a 2-hop, network knowledge. In Context-Aware Routing protocol (CAR) [START_REF] Musolesi | CAR: Context-Aware Adaptive Routing for Delay Tolerant Mobile Networks[END_REF] and GeOpps [START_REF] Leontiadis | GeOpps: Geographical Opportunistic Routing for Vehicular Networks[END_REF], the delivery probabilities are computed using both utility functions and Kalman filter prediction techniques. CAR allows each host to associate a utility function, representing the delivery probability, with every other host. These functions process contextual properties, and especially the change rate of connectivity and the probability of being located in the same cloud as the destination. Instead of using the current contextual information, CAR uses predicted future values of the context in order to have more realistic values. CAR assumes an underlying MANET routing protocol that connects nodes together in the same MANET cloud. To reach nodes outside the cloud, a sender looks for the node in its current cloud with the highest message delivery probability. In CAR, nodes compute delivery probabilities proactively, and disseminate them in their ad hoc cloud (they are exchanged with standard routing tables). Therefore, the contextual properties are exploited to evaluate probabilities just for those destinations that each node is aware of. GeOpps, which is a geographical delay tolerant routing algorithm, exploits information from the vehicles' navigation system to route messages to a specific location. To select the next packet carrier, neighbor vehicles that follow suggested routes to their driver's destination calculate the nearest point that they will get to the destination of the packet. Afterward, they use the nearest point and their map in a utility function that expresses the minimum estimated time that this packet would need in order to reach its destination. The vehicle that can deliver the packet quicker/closer to its destination becomes the next packet carrier. Like CAR, HiBOp [START_REF] Boldrini | Hibop: a history based routing protocol for opportunistic networks[END_REF] also exploits context information in order to compute delivery probabilities. However, HiBOp can be perceived as being more general than CAR since it does not require an underlying routing protocol, and because it is also able to exploit context for the destinations that nodes do not know. HiBOp uses history information in order to improve the delivery probability accuracy, and does not make predictions as CAR. HiBOp maintains an history table in order to estimate the probability of encountering a contextual property in the near future. Propicman [START_REF] Nguyen | Probabilistic routing protocol for intermittently connected mobile ad hoc network (propicman)[END_REF] also exploits context information and uses the probability of nodes to meet the destination, and infers from that the delivery probability, but in a different way. When a node wants to send a message to another node, it sends to its neighbor nodes the information it knows about the destination. Based on this information, the neighbor nodes compute their delivery probability and return it. The node that wants to send the message will send it only on the two-hop route(s) with the highest delivery probability, if this one is higher than its own. Propicman considers that mobile nodes carried by people are not likely to move around randomly, but in a predictable fashion based on repeating behavioral patterns at different timescales (day, week, month). If a node has visited a place several times before, it is likely to visit this location again in the future.

Prophet (Probabilistic Routing Protocol using History of Encounters and Transitivity) [START_REF] Lindgren | Probabilistic routing in intermittently connected networks[END_REF] has an approach similar to that of Propicman. However, according to our classification Prophet is in the third category (i.e., the protocols implementing neighbor-knowledge-based methods). Indeed, the selection of the best neighbor node is based on how frequently a node encounter another. When two nodes meet, they exchange their summary vectors, which contain their delivery predictability information. If two nodes do not meet for a while, the delivery predictability decreases. When a node wants to send a message to another node, it will look for the neighbor node that has the highest amount of time encountering the destination, meaning that has the highest delivery predictability to the destination. This property is further transitive.

All the above-mentioned protocols implement various strategies aiming to select the next best hop(s) to deliver a given message. These protocols could probably be used to perform service invocation, which traditionally relies on a unicast communication model, but not to achieve service discovery, which requires an efficient broadcast of service discovery requests and service advertisements. Indeed in order to avoid the broadcast storm problem and a network congestion, these messages must not be broadcast in a blindly epidemic manner, but instead in given geographical areas by protocols implementing both a delayed message forwarding and selfpruning heuristics, thus making it possible for mobile nodes to decide whether they should forward a message or not.

Over the last years, several broadcast protocols implementing such as features have been proposed for connected MANETs. Willians and Camp [START_REF] Williams | Comparison of Broadcasting Techniques for Mobile Ad Hoc Networks[END_REF], and later Stojmenovic and Wu [START_REF] Stojmenovic | Broadcasting and Activity-Scheduling in Ad Hoc Networks[END_REF], have proposed classifications of broadcast protocols for MANETs. They have identified four categories: simple flooding, probability-based methods, area-based methods and neighbor-knowledge-based methods. Like the epidemic forwarding protocol, the simple flooding [START_REF] Ni | The Broadcast Storm Problem in a Mobile Ad Hoc Network[END_REF] does not make any attempt to reduce the number of broadcasting nodes and does not require any knowledge. In probability-based methods, each node estimates a potential contribution to the broadcasting process. If this contribution is lower than a given threshold the message is not forwarded [START_REF] Wu | Broadcasting in ad hoc networks based on selfpruning[END_REF]. In Area-based and locationbased [START_REF] Ni | The Broadcast Storm Problem in a Mobile Ad Hoc Network[END_REF] methods, a message can be re-emitted only if the surface or the distance of this area is upper than a given threshold. Finally, neighbor-knowledge-based methods require the knowledge of the 1-hop or 2-hop neighborhood. DFCN (Delayed Flooding with Cumulative Neighborhood) [START_REF] Hogie | A bandwidthefficient broadcasting protocol for mobile multi-hop ad hoc networks[END_REF] is an example of such protocols. It proposes a bandwidth-efficient algorithm that introduces a benefit parameter allowing hosts to decide if a message must be broadcast or not according to their neighborhood status and to a random delay that aims at avoiding message collisions. These above-presented broadcast protocols are designed for connected MANETs and are not thus suitable to opportunistic networks that are intrinsically disconnected. Moreover, these protocols do not make it possible to restrict the propagation of messages geographically or to define the area where they must disseminated.

III. PROTOCOL OLFSERV

In this section we present OLFServ (Opportunistic and Location-based Forwarding protocol for Service provision), a protocol that supports both the discovery and the invocation of location-aware services in opportunistic networks, and especially in those formed by fixed infostations providing software services and by mobile devices carried by humans and behaving as service clients (see Figure 1). OLFServ addresses three main issues: 1) the efficient and geographicallyconstrained broadcast of service discovery messages, 2) the selection of relevant carriers to perform service invocation, and 3) the support of clients' mobility during the service invocation process. The remainder of this section presents the assumptions on which OLFServ relies, the heuristics it implements and how service discovery and invocation are performed using scenarios.

A. Requirements

Protocol OLFServ relies on 4 main assumptions: 1) Mobile hosts and fixed infostations are aware of their geographical location and able to compare their location with that of another host (for example with that of the recipient of a given message). 2) Mobile hosts are able to perceive their one-hop neighborhood. Such a piece of information is obtained using specific messages (beacon messages). 3) Service providers and clients are expected to include in messages the network address of the destination, but also

• the geographical location of the destination,

• their own location,

• a date of emission,

• a lifetime,

• a number of hops (equivalent to the TTL field of the IP protocol),

• the geographical area where messages can be disseminated,

• and their speed and direction (or destination) if they are mobile. Intermediate nodes (or at least the last forwarder) are expected to add their own location in messages before forwarding them. Messages must also include the identifiers of the 1-hop neighbor nodes of the mobile host that has originally sent the message in the network or that has forwarded it recently. 4) Each mobile host is able to store the messages it receives temporarily, and able to associate to each of them some information, and especially the identifiers of nodes that are known to have received them.

B. Overview of OLFServ

In OLFServ, communications are currently performed on an IP multicast group. In order to avoid message collisions and the broadcast storm problem when sending messages in the network, OLFServ implements a delayed message forwarding and location-based and neighbor-knowledge-based self-pruning heuristics.

The forwarding delay computed by the mobile hosts has two components: a component that is inversely proportional to the distance from the last forwarder and a random component. For this purpose, we have divided the wireless communication range of devices in several rings, and we have assigned to each ring a range of values so that the values of ring i are greater than those of ring i+1 (see Figure 2). When receiving a message, devices are expected to infer in which ring of the message sender they are located, and then to compute a random value included in the range of values of this ring. The farthest nodes from the sender are therefore expected to forward the message before the closest nodes, thus improving the geographical propagation of messages. Moreover, a node will cancel its message forwarding process if it perceives that all of its neighbors have already received the message it plans to forward, thus limiting the number of messages that roam in the network.

In OLFServ, the forwarding process is driven by four main events, namely the message reception, the expiration of the message delay forwarding, the detection of a new neighbor and the location changes. It is also constrained by location aspects. When receiving a message, hosts running OLFServ compute a new forwarding delay. When this delay has expired, the mobile hosts, which are in the area where this message can be disseminated, will forward this message only if some of their neighbor nodes have not received the message yet. Otherwise, it stores the message locally until it expires.

If the message must be broadcast in a whole given area, as is the case for both service advertisements and service discovery requests, all of the nodes composing the 1-hop neighborhood boundary of the message sender are thus likely to forward the message with the respect to above-presented constraints. In contrast, for service invocation requests and responses, for which the location of the recipient is known, only the nodes that are closer to the recipient than the previous hop are expected to forward the message towards the destination. Thus, we progressively refine the area where the message is broadcast until reaching the destination, and we drastically reduce the number of messages that are replicated in the network while having a good message delivery ratio. Mobile hosts that receive a response in reply to a service invocation request they forwarded in the past do not forward this request again in the future.

Finally, a mobile host will compute a new forwarding delay in order to send a message when its context changes. For instance, when it reaches the area where the message can be disseminated or when a new node appears in its vicinity. These properties are used in the invocation process in order to broadcast the response in the area where the client is supposed to be.

C. Service discovery

In order to understand how service discovery is achieved with OLFServ, let us consider the scenario of Figure 2. In this simple scenario, we consider a disconnected MANET composed of a fixed infostation I and a set of mobile devices carried by humans. Infostation I provides a service that is relevant only in a given geographical area, which is represented by the yellow circle in Figure 2. We assume that one of them is interested by the service proposed by I, namely node N 7 . We also assume that this network, which is composed of three distinct islands in Figure 2, evolves in an unpredictable way according to the nodes' mobility. We make no assumptions about the mobility of people. For instance, some of them can follow regular mobility patterns, while others can have a random mobility behavior or can be stationary.

In order to use the service offered by infostation I, N 7 must discover this service first either proactively or reactively. For the sake of illustration, let us consider this last discovery scheme and let us suppose that infostation I has injected in the network an advertisement A including its location, the geographical area where A must be disseminated, a date of emission, a lifetime, a number of hops and the set of nodes that are expected to receive this advertisement (i.e., I, N 1 , N 2 and N 3 ). Advertisement A will thus be received first by nodes N 1 , N 2 and N 3 . These nodes will store A locally, and then they will compute a forwarding delay in order not broadcast message A simultaneously. Since the forwarding delay is a random value that is inversely proportional to the distance from the last forwarder (in the present case infostation I), nodes N 1 and N 3 will compute a delay in a same range of values. This value will be greater than that computed by node N 2 .

Since infostation I has specified in advertisement A that hosts I, N 1 , N 2 and N 3 are supposed to have received this message, nodes N 1 and N 2 will not broadcast message A. Message A will thus be forwarded by nodes N 3 , N 4 and N 5 successively and will be received by N 6 . Thanks to the "store, carry and forward" principle implemented in OLFServ, and which exploits nodes' mobility and contact opportunities, this advertisement will be propagated in the whole area specified by the infostation, and only in this area. Indeed, the selfpruning heuristics implemented in our protocol enable mobile devices to forward messages only if they are in the areas where these messages are relevant. For instance, node N 6 , which left the island of the infostation in order to join that of the client (node N 7 ), will broadcast advertisement A in this new island. This message will be then received by node N 7 , which will also forward this message. Then, this message will be received by node N 8 . This node will not disseminate message A because it is outside of the area defined by infostation I. Message A will not be received by node N 9 because it is outside of the communication range of node N 7 .

D. Service Invocation

Once it has discovered the service offered by infostation I, node N 7 may invoke this service by sending an invocation request that will include especially the location of the provider, its own location, and its speed and direction. This invocation request will be received by nodes N 2 and N 6 , and will broadcast only by node N 2 , because it is closer to infostation I than node N 7 . This request will be then received by infostation I (see Figure 3).

Since N 7 has specified its location, its speed and its moving direction, infostation I can estimate the area where the client is expected to be when it should receive the response. So when returning its response, infostation I includes its location and the area where N 7 is supposed to be. For instance, the orange circle in Figure 3. The response will be then routed towards this area using a forwarding scheme comparable to that used for the invocation. When the response reaches this area, the precise location of N 7 is not known. Therefore, the response will be disseminated in the whole area following a broadcast scheme comparable to that used for the service discovery in order to improve the service delivery. When they receive a response for an invocation they have stored previously, the mobile devices remove the invocation message from their local cache and stop to forward it.

IV. EXPERIMENTS AND RESULTS

In order to evaluate our protocol, we have made a series of simulations using the Madhoc simulator1 , a metropolitan ad hoc network simulator that features the components required for both realistic and large-scale simulations, as well as the tools essential to an effective monitoring of the simulated applications. This simulator, which is written in Java, allow us to run our middleware platform on it.

A. Experiments

The simulation environment we consider is an area about 1km 2 , which is composed of 6 squares of 100m × 100m linked by predefined paths. An infostation is installed in each square. An infostation offers only one service. A same service is offered by two distinct infostations. The discovery and invocation area of each infostation is a circle with a radius of 300m. This environment is populated with 100 mobile devices equipped with Wi-Fi interfaces, which move along predefined paths in order to reach a given destination located in one of the predefined squares. The destinations of nodes are chosen randomly. When it has reached its destination, a mobile device strolls during 5 minutes in the square where it is located. Afterward, it chooses another destination. Each mobile host moves at a average speed between 0.5 and 2 m/s. The communication range of mobile devices and of infostations varies from 60 to 80m. Among the 100 mobile hosts, 70 act as clients of the services offered by the 6 infostations, whereas the 30 others only act as intermediate nodes. In our scenario, the infostations broadcast service advertisements every 30 seconds. Once they have discovered a relevant service, the clients invoke this service periodically (every 5 minutes) with a different request. All messages have a lifetime of 10 minutes and a number of hops of 8 hops. Messages are forwarded only in areas where they are relevant. Moreover, our middleware platform selects the best service provider when several service providers have been discovered (i.e., the closest provider) [START_REF] Sassi | Towards an Opportunistic and Location-Aware Service Provision in Disconnected Mobile Ad Hoc Networks[END_REF]. In OLFServ, beacon messages are sent every 20 seconds.

B. Results

In this section, we focus on a particular experiment whose objective was to measure the ability to satisfy the client invocation efficiently. For this purpose, we compared the service provision with OLFServ, with one implementing a 1hop discovery and invocation model, and with another one implementing a purely epidemic discovery and invocation model. In the 1-hop model, clients must be in the vicinity of a provider offering the service they require in order to discover and to invoke this service. In contrast in the epidemic model, messages are forwarded by mobile devices opportunistically and without any geographical constraints. In this model, service can be discovered and invoked from anywhere.

Figure 4 presents the network load for the simulations we have made, and in Table I we summarize the simulation results we have obtained. In contrast to the invocation with both the epidemic protocol and OLFServ protocol, the 1-hop service provision protocol has the advantage of offering a low network load, but to the detriment of the service provision since only 1,8 invocation requests have been made in average by clients during the simulation (i.e., 1 hour). As shown in Figure 4 and Table I, the OLFServ protocol has approximately the similar invocation success ratio as the epidemic protocol (0,94*4,7= 4,418 for OLFServ and 0,55*0,6=5,28 for the epidemic protocol), but with a lower network load. We also observe that OLFServ provides better results regarding the invocation delay than the epidemic routing protocol. Indeed, with the epidemic protocol the provider can be invoked from anywhere. Some clients can thus obtain some responses even if they are outside the areas specified by the providers, but with a more bigger delay introduces by the store, carry and forward principle. Moreover in the epidemic protocol, mobile devices broadcast their messages periodically every 20 seconds, whereas in OLFServ, messages can be forwarded by mobile devices as soon as they receive new messages. Moreover, our protocol makes it possible to support the mobility of clients, and thus to improve the delivery of responses.

V. CONCLUSION AND FUTURE WORKS

We presented in this paper a new opportunistic and locationbased forwarding protocol supporting service discovery and invocation in disconnected mobile ad hoc networks. This protocol addresses three main issues, namely the efficient and geographically-constrained broadcast of service discovery messages, the selection of relevant carriers to perform service invocation, and the support of clients' mobility during the service invocation process.

Despite the good preliminary results we obtained, we plan to investigate new directions with OLFServ. For instance, we plan to implement source routing mechanisms in order to help in the selection of intermediate nodes when forwarding service responses. Finally we would investigate how a nhop neighborhood perception can contribute to improve the self-pruning heuristics we have implemented in our protocol (typically with n=2). In the future, we also plan to perform more evaluation, and especially by varying the number of mobile devices in our simulation (for instance between 100 to 1000).

Figure 1 .

 1 Figure 1. Example of a disconnected MANET.

Figure 2 .

 2 Figure 2. Service discovery with OLFServ in disconnected MANETs.

Figure 3 .

 3 Figure 3. Service invocation with OLFServ in disconnected MANETs.

Figure 4 .

 4 Figure 4. Network load. 1-hop epidemic OLFServ Number of clients having discovered a provider 52/70 70/70 70/70 Average delay of success invocation(sec.) 3 122 39 Standard Deviation for delay of success invocation (sec.) 2 94 21 Average invocation success ratio 0,96 0,55 0,94 Average of number of invocations 1,8 9,6 4,7

Table I SIMULATION

 I RESULTS FOR SERVICE DISCOVERY AND INVOCATION.

http://agamemnon.uni.lu/~lhogie/madhoc/