Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: P1.16 - Optimal Eulerian model for the simulation of dynamics and coalescence of alumina particles in solid propellant combustion
Full Citation
Permanent Link:
 Material Information
Title: P1.16 - Optimal Eulerian model for the simulation of dynamics and coalescence of alumina particles in solid propellant combustion Computational Techniques for Multiphase Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Doisneau, F.
Laurent-Nègre, F.
Murrone, A.
Dupays, J.
Massot, M.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
Subject: dilute polydisperse spray
solid propulsion
Eulerian multi-fluid model
Abstract: 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.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00439
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: P116-Doisneau-ICMF2010.pdf

Full Text

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

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

Frangois Doisneau*t, Fr6d6rique Laurent-Negret, Angelo Murrone*,

Joel Dupays* and Marc Massott
Department of Fundamental and Applied Energetics, ONERA, 91120 Palaiseau, FRANCE
SLaboratoire EM2C, Ecole Centrale Paris and CNRS, 92295 Chatenay-Malabry, FRANCE,,, and
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.


Two-phase flows constituted of a gaseous phase and a
disperse condensed phase play a key role in many in-
dustrial 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 ap-
plications the disperse phase is composed of particles of
various sizes that can eventually coalesce or aggregate,
break-up, evaporate and have their own inertia and size-
conditioned dynamics. So the importance of polydisper-
sion is obvious for a full modelling of these phenomena.
In rocket boosters, aluminum powder is frequently

used as solid propellant additive to increase specific im-
pulse. Unlike the other ingredients, aluminum parti-
cles can burn in a sigmlik.iii portion of the chamber
and produce a condensed liquid disperse phase of alu-
mina. This disperse phase encounters drag forces, coa-
lescence and heat exchanges (Dupays et al. 2000). In the
nozzle, the droplets accelerate suddenly and cool down
with the gas, becoming solid particles which break-up
in such velocity gradients (Dukowicz 1980; O'Rourke
1981; Amsden et al. 1989). 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 insula-
tion 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 ve-
hicle instability. In such harsh conditions of pressure,
temperature and velocity, solid propulsion experiments
consume high technology materials and offer poor mea-
surement output, especially on the disperse phase. The
abundance of physical phenomena involved makes the
models difficult to scale (T6th 2008). 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 accumula-
tion. 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 some-
times "kinetic" in reference to the kinetic theory of
gases- describes the droplets as a cloud of point parti-
cles 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 cen-
ter 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. Williams
(1958, 1985) proposed a transport equation based on
kinetic theory that has proven to be useful for treating
polydisperse, dilute and moderately dense liquid sprays.

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

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 in-
teraction of coalescence (Kuentzmann 1973; Hylkema
1999; Laurent and Massot 2001; Laurent et al. 2004).
There are several strategies in order to solve the liq-
uid phase and the major challenge in numerical simula-
tions 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 (Dukowicz 1980; O'Rourke 1981; Ams-
den et al. 1989; Hylkema 1999; Rilger et al. 2000). It is
called Direct Simulation Monte-Carlo method (DSMC)
by Bird (1994) and is generally considered to be the
most accurate for solving Williams equation; it is spe-
cially 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 dis-
perse phase, thus offering limited possibilities of vec-
torization/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 de-
scription. Moreover for unsteady computations of poly-
disperse 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 Laurent and Massot (2001); Laurent et al.
(2004) from the ideas of Greenberg et al. (1993), 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
(Laurent and Massot 2001). 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 Eu-
lerian approaches often consider discrete droplet sizes
gathered into .L,,,s, which cannot account for the
new droplet sizes created by coalescence (except scarce
examples such as in Vasenin et al. (1995)), 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
(Bouchut (1994); 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 ( I.iiciii.illiiii
(2009). The model finally requires closure relations for
the phenomena accounted for by the Williams equation.
We refer to Abramzon and Sirignano (1989); Laurent
and Massot (2001); Laurent et al. (2004) 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 capa-
bility to simulate the size-conditioned dynamics of poly-
disperse sprays including coalescence with a first order
resolution method of the size distribution in each sec-
tion (Laurent and Massot 2001) even if the number of
sections with the original method has to be large for ac-
curacy 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 Laurent (2006). 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 (Laurent et al. 2004) while the second or-
der method has been extended to coalescence at the cost
of some more elaborate algebra in Dufour (2005). It can
yet easily include advanced models such as collision and
coalescence efficiencies developed in Hylkema (1999);
Hylkema and Villedieu (1998).
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, compo-
sition, 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

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

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 pro-
vides fully coupled aero-thermochemical resolution
for energetic 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 ded-
icated 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 Laurent
et al. (2004). In section three, we validate the second or-
der Multi-Fluid method and its collision and coalescence
efficiency models on a 1D coalescence efficiency exper-
iment 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 reason-
able time a complete solid propulsion test case, we give
some numerical results in a modelled chamber and noz-
zle using the CEDRE code before concluding the study.

1 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 pre-
suming the size distributions on the sections. The origin
and assumptions of the corresponding coalescence mod-
els are widely detailed.

1.1 A kinetic description : the Williams

Number density function When focusing on poly-
dispersity, the size parameter p of droplets is of capi-
tal importance but its natural expression depends on the

7th International Conference on Multiphase Flow,
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 = 4ir3
1 S3
6, S 2 allows the size to be expressed in this paper with
the most comfortable notations.
Let us define the NDF function f9 of the spray, where
f (t, x, u, y)dxdydu 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 do-
neighborhood of p. As for the size parameter conver-
sions, we shall keep in mind that f (r)dr fs(S)dS
f" (v)dv and use the implicit notation f.

Williams equation The evolution of the spray is de-
scribed by the Williams transport equation (Williams
1958). For a non evaporating case, it reads :

Jf + u Oxf + 0 (F f) = (1)

where F is the drag force per unit mass and F 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

r (S)

187 p9

where U is the gas velocity, pg represents its dynamic
viscosity and pi is the droplet (supposedly liquid) mate-
rial density.

Coalescence operator The kinetic modelling for
the collision operator leading to coalescence is taken
from (Hylkema and Villedieu 1998). We then assume :

[Cl] 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 pre-

Thus, r Q+ Q where Q+ and Q respectively
correspond to the quadratic integral operators associated
with creation and destruction of droplets due to coales-
cence. These quadratic operators read (Hylkema and

Villedieu 1998; Hylkema 1999):

Q + i f (t, X, U (, *, U), V (, V*))x
2 J1u*,v*n[o,v]
f(t, x, u*, v*) (|u u u* , v*)Jdv*du*

Q- j(, x, u, v)f(t, x, u*, w) x
T(|u u*|, v, v*)dv*du*
where v(v, v*) v-v* and u = vu** arethepre-
collisional 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 :
T(|u-u*|, V, J*) -= &(|u-u* ,v,*) iu-u*| B(v, *)

In this kernel, B(v, v*) = (r + r*)2 is the impact pa-
rameter and = con.coal accounts for the collision
efficiency (col and the coalescence efficiency (coal. The
modelling of these efficiencies is detailed in the follow-
ing two paragraphs.
Collision efficiency Collision efficiency Ecol 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 p, or viscosity pg.
In the case of unbalanced droplet sizes, D'Herbigny
and Villedieu (2001) select two collision efficiency mod-
els : the Langmuir-Blodgett model, and the Beard-
Grover model (Hylkema 1999; Achim 1999). These
laws strongly depend on the bigger droplet Reynolds
number and on a dimensionless number k which read,
when taking rl the smaller radius and r2 the bigger ra-
dius :
Re -2P2 U U2 k pl u2 ul
sg 9'Pgr2
The number k is the ratio of the smaller droplet relax-
ation time T, = 2pj to its residence time in the bigger
droplet influence zone T7 = .r

In the case of low Reynolds numbers, Langmuir and
Blodgett (Langmuir 1948) numerically get the following
expression :

Er (k)
E (k)

31n(2k ) 2
+ 2(k 1.214))

if k < 1.214

whereas for high Reynolds, they get:

SE (k)
E2 (k)

0 if k < 0.0833
(k+o5)2 otherwise

For intermediate cases they assume the following inter-
polation :

FLB El(k) (Re/60))E (k)
Col(k, Re) 1 + Re/60 + 1 + Re/60

Beard and Grover (1974) suggest to increase the ac-
curacy of formula (3) which results from a simple inter-
polation between two limits. For this purpose, they use a
numerical solution of the incompressible Navier-Stokes
equations in order to determine the gaseous flow sur-
rounding 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 E [0,400] and ri < r2, it yields:

colj(k, Re) = [arctan(max(H(k, Re),0))] (4)


H = 0.1465 + 1.302Z 0.607Z2 + 0.293Z3
Z ln(k/ko)
ko= exp(-0.1007 0.358 n(Re) + 0.0261[ln(Re)]2)

Typical collision efficiency values with these two laws
are in Fig. 1.





0 05 1 15 2 25 3
Velocity difference |u2 ui (m/s)

Figure 1: Two collision laws selected for booster appli-
cations with r2 150pm (red : rl = 2/m,
green: rl = 3pm; blue : rl 4pm)

Coalescence efficiency Coalescence efficiency
(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 ki-
netic 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
(Achim 1999; Brazier-Smith et al. 1972) are commonly
used. Especially, the CEDRE code provides this model

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

but it is not used in the following simulations. There-
fore, 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 Laurent and Massot (2001). 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.

1.2 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 = fdu and u f ufdu/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
u(t, x, S).

[H2] The dispersion around the averaged velocity
u(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)6(u u(t, x, S))

The monokinetic ]i p' ,lle',i [H1] is equivalent to reduc-
ing 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 (Massot 2007). Yet it pre-
vents droplet-crossing, which can occur in solid propul-
sion during ejection from fast vortices or high speed
parietal injection.
The zero-dispersion h] p, llwc'is [H2] is justified when
neither turbulence nor brownian phenomena interact
with the disperse phase.
These ]h.lp'lllw, 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 conser-
vation 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 ve-
locity. The semi-kinetic system reads :

3tn + x, (nu) = cQ Q,
9t(nu) + 9x (nu (F u) = nF + Q'

Semi-kinetic coalescence operator In the semi-
kinetic 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(o (,, -*))n(V*)B3(o I*) do*
n 2 J v* [0,]
Qn -n(v) n(o*)B3(, o*)I do*

J2 v* [O ,+oo

Q.n(O) n(v*)B(vo, v*)Iu d*[
*C[ o,+ oo[
where I, I,, I and 1+ are the partial collisional inte-
grals. To preserve the monokinetic and zero-dispersion
assumptions, these integrals are computed with average
velocities and therefore read :

I+ =lu(v*) u( )| (I(V*) u( )|, v, v*)
I =u(o) U()| &(u|() U(o)|, o, o)

Su(<) (|w) Vu(o*) (|u() -u( )|, o, o*)
fU(ou) U(W)1 CO(W) U(O), V, ow)
4U- U(V) U(V) U(*)| C(IUM U(V*)|, V, V*)

1.3 Multi-Fluid model

In a second step we choose a discretization 0 = So <
S1 < ... < SN for the droplet size phase space and
we average the system of conservation laws over each
fixed size interval [Sk 1, Sk [, called section. The set of
droplets in one section can be seen as a "fluid" for which
conservation equation are written, the sections exchang-
ing 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 :

SCr [u(t, X ix, S)
VS [ -1, S [ n(t, x, S)

Kk(t, x, S)

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

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 conver-
gence with the number of sections. It induces a conve-
nient integration simplification but it can generate tra-
jectory 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) 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 IK
function in each section yields a first order granulome-
try convergence with the number of sections. Whereas in
the second order method, a two-coefficient 2Kk 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 polydis-
persity, as told in introduction.

1.4 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) = k (t, x) k(S)

where mk, is the mass concentration of droplets in the
kth section, in such a way that


k(S) Pi S/dS = 1.
6 v

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 ,igiilik.iIl part of the mass.
This is a major limitation compared to the second order
Equations The conservation equations for the kth
section result from the integration of the mass moment
of the semi-kinetic system (5) in each section k and
reads :

dtmuk+9x-(mkUk) u=k
9t (m kUk )+ x- (mk Uk 0gUk)

- M C



vj+: \ -----

0 vi Vi+01 v*
Figure 2: Coalescence terms are integrated on the
"(i, j, k) doii.nmiii

where Fe is the "classical" averaged drag force per unit
mass on a section (Greenberg et al. 1993; Laurent and
Massot 2001):

SU a 1 s 1 l (S)S3/2
F T T Sk 61/ I (S)

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 depen-
dency will not coincide with the section after mapping
the natural variables of the two precursor colliding part-
ners. Thus splitting the integration domain thanks to
the section continuity yields elementary integrals Qijk,
triply indexed with the two precursor section numbers
(i, j) and the destination section number (k). These do-
mains are illustrated in Fig. 2.
In the particular case of the first order method, con-
sidering [H2] and assuming co,, = coal 1 (Laurent
and Massot 2001), the coalescence integrals Qijk take
the following form after factorizing the mass moments

Qi*JI J ( JI 2}+ r < 7,. ,,ddr*dr+
Q dv*<+v *[vV2+l4 3
f k= 1 K(r*)wlI (r )w(r*+ r )2 pir dr*dr
S +v*+v'C [vk,vk+l[

As for the disappearance terms, they can also be com-
puted as sums of the elementary creation integrals and
must be so to ensure the conservation of matter and en-
ergy. After some algebra, the coalescence terms 1C',

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

Ct 1C+ and 'C" read:

1(k) N
crn Zm Zm luj (1Q^I,'Q^I,)
i=1 j=1

j-l i-I
1(k) N
1CrU+ i= i j *.j I (U dQ+ + IQk)
ik1 : 1

1Cku- = meUfk E mj u,



The integrands of the Qsjk integrals depend only on size
parameters r,r* and r' thanks to [H4]. This allows the
Qijk integrals to be pre-calculated as soon as the section
limits and the (IKsk) are given, i.e. once and for all at
the beginning of a simulation.

1.5 Second order method

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

,((t, x, S) a (x, t)exp(-bke(x, t)S)

where (a (x, t), bk (x, t)) I ensures

S kl 2kI(t, x, S)dS = nk(t, x)

Sk 2(t, S) S3/2dS- = m(t, x)
S _k-1 6

The choice of an exponential function ensures the pos-
itivity 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 sec-
tions i.e. the fact that (ma, nk) couples are conditioned
by the section boundaries since they drive information
on average droplet volume.

Equations The conservation equations for the kth
section now read:


Ot(mkUk)+9x. (mu uku )= m C+)U)C "-+- 2C"T
The total system thus counts three times as much equa-
tions as the number of sections.

dtnl+d,(ncul) 2Cn+
09mfe+89x -(mUk) k

Second order coalescence terms Since the time
and space dependency of the size-distribution functions
Kk(t, x, S) is no longer factorizable as was mk in the
first order method, the Qijk 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 colli-
sion/coalescence efficiency at the same time. That is
why this method is "less inappropriate" for efficiency
models. Let us take the notation :

i (t, x, r*r u -u*|) = (t, x, S*)l (t, x, S )x
(r* + r )2((r*, ru, u_ U* ) i Uj

The coalescence integrals have a different homogeneity
than in the first order method. They now include the
number, mass or momentum information and read :

Q ijt J 'x(t, r, |u,
iQ Jj*+V[Vk,Vk+1
2Q* [4 x r[
J Jv*+v'>[vk,vk+1[
11 ij (t, x, r*, r ,
J Jv*+v'[- vk,vk+ [


u*) ,. ddr*dr

U*) w4Tplr dr*dr

The coalescence source terms 2+', k2 2C(', ,
2C" +and 2C are still written as direct sums of the
"(i, j, ) iinici.iI, but therefore read:

I(k) N N N
k 11 2. ZZk, ki kji
i lj 1 j i= 1
I(k) N N N
2Cr-E EZQ 2Q ZZ EE1(Q,+ t Q<;)
i=i"=1 j =i 1
i KH
2cmu *Zu 2, 2cm-
kc 1 ZZ (Uj2QIti~uUk6,2CTUk k k
j=l kKL

1.6 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 >,, 2Q1 and Qij. The second order method
requires indeed a costly numerical integration for each
time step and space cell. We have chosen to per-
form these integration with a 2D 5 point Newton-Cotes
quadrature which requires i3 (t, x, r*, r |u u*|) to
be computed 52 x N,2 x Ne,1 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 col-
lision 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. Laurent
et al. (2004) introduces a general form for coalescence

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

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 pa-
rameters. That is why the first order method here ex-
cludes such collision efficiency models.

2 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 dynam-
ics of a coalescing cloud. We therefore need a well-
suited test case, inducing coalescence as well as size-
conditioned dynamics and difficult enough to highlight
the limitations of the methods.
For that purpose, we have chosen the configuration
which was used in Laurent and Massot (2001) to vali-
date the first order method compared to a reference La-
grangian solver : a 2D axisymmetrical conical decel-
erating nozzle, designed in such a way that it admits,
for one-way coupling spray dynamics a self-similar so-
lution. We though inject a lognormal distribution instead
of a bi-modal one.

2.1 Definition of configuration

The chosen configuration is stationary 2D axisymmetri-
cal in space and 1D in droplet size. It is described in
detail in Laurent et al. (2004). Hence, only its essential
characteristics are given here.
For the problem to be one-dimensional in space, con-
ditions for straight trajectories are used and are compat-
ible with the assumption of an incompressible gas flow.
This leads to the following expression for the gaseous
axial velocity U, and the reduced radial velocity U,/r,
for z _> z :

U U(z)

S1(Zo) U, U,


r Z Z

where zo > 0 is the coordinate of the nozzle entrance
and the axial velocity U(zo) 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 con-
ical 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 noz-
zle is taken as a(zo) = -2U(zo)/zo; it is chosen large
enough so that the velocity difference developed by the
various sizes of droplets is important. Laurent and Mas-
sot (2001) chose rather large values, as well as strong de-
celeration, leading to extreme cases: U(zo) 5 m.s-1,

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

zo = 5 cm : these values generate a very strong cou-
pling between coalescence and droplet dynamics. These
severe conditions make the test cases under considera-
tion 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 sec-
ond order method remains comparable to the first order
method, which lacks this sort of model. We thus assume
Ecoi coal 1 in this section.

2.2 Droplet initial distribution

The lognormal distribution is useful to characterize alu-
mina 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 equi-
librium temperature 3600 K (corresponding to an infinite
conductivity model), does not change along the trajecto-
For the sake of the first order method accuracy, we
do not want to transfer too much mass in the last sec-
tion. The initial injected mass concentration is therefore
taken as mo 1.0609 kg.m 3. For the same reason, the
lognormal parameters are set to SLN 1600pm/2 and
8LN = 1.5 which corresponds to a sharp distribution
centered on a radius of 11.3pm. We shall yet keep in
mind that the second order method can accurately treat
cases where the final section hosts a ,igniliik.i 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 specifica-
tions 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 com-
pared 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 Laurent and Mas-
sot (2001).

2.3 Distribution spread up results

The processing of polydispersion gives a first indica-
tion 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

1 5 sections


a o

.o 008
00 006

" 004

U 002

- 008

I 13 sections |

I 25 sections I

0 10 20 30 40 50 60
Radius (pm)
Figure 3: Mass concentration distribution at the noz-
zle's end computed with the first and second
order methods (empty and filled symbols) and
reference (black)

roughly converged. The growth overestimation of the
first order method for poor size discretization brings
about major consequences for the spray dynamics.

2.4 Spray dynamics results

Let us now appreciate the consequences of polydisper-
sion treatment on dynamics. Fig. 4 shows the evolution
of mass and number concentrations along the center-
line 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 concentra-
tions are computed considering section average droplet
volumes given by integration of the (1k)k functions :
data are therefore redundant with mass concentration
data. We emphasize the fact that the convergence of the

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



o 02


1 05

5Y 1

1 05

- 09

0 85 L

01 015
Position (m)

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

Groups (im) 5 sections 13 sections 25 sections
Gi=[0, 12.5] 1 1 to 3 1 to 6
G2-[12.5, 25] 2 4 to 6 7 to 12
G3-[25, 37.5] 3 7 to 9 13 to 18
G4=[37.5, 50] 4 10 to 12 19 to 24
Gs (> 50p/m) 5 13 25

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 il-
lustrated in table 1. The evolution of the mass concen-
tration 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 G5 where no mass should be found conse-
quently to the problem parametrization. Finally, the sec-
ond order method reveals to have acceptable error with
as few as 5 sections.

016 018 02 022 024
Position (m)

Figure 5: Zoom on the convergence of mass concen-
tration at the nozzle's end (black : rough la-
grangian validation)

5 sections
13 sections
25 sections
53 sections

I 1st order MF 2nd order MF
Is 10s
1.5 s 50 s
2s 180 s
5 s n.a.

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

2.5 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 typ-
ical 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 be-
comes a costly operation, restraining the profit of the
first order method.

3 Coalescence Experimentation : the
D'Herbigny Test Case

In this section, we provide further validation for the sec-
ond 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-

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

Empty First order
SFilled Second order
0 / A 5 sections
S 13 sections
S25 sections
005 / Reference

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

tributions. Indeed, modal distributions (i.e. monodis-
perse or discrete distributions) are rather suited for the
Lagrangian point of view while they introduce stiff size
distributions in Eulerian modelling.

3.1 The D'Herbigny collision efficiency

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 r2 150pm thanks to piezoelectric
injector which is fairly accurate. The smaller droplets,
composing the fog, have a radius ri e [2/m, 4pm].
The two types of droplets have an approximately con-
stant 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 ef-
ficiency models. They are given in table 3.
Moreover, the collisional Reynolds number Re is
around 60 which is comparable to the values in the Ari-
ane 5 P230 booster. This range of values (moderate Re
and low k) is supposed to be favorable to the Beard-
Grover law.

Figure 7: D'Herbigny experimental device

rl 2pm 2, 5pm 3pm 3, 5pm 4pm
k 0.78 1.22 1.75 2.39 3.12

Table 3: Collision parameter k values in the
D'Herbigny experiment

3.2 Results

Droplet radius growth study D'Herbigny and
Villedieu (2001) provide an analytic formula to esti-
mate the Sauter mean radius growth Ar2 of the bigger
droplets when assuming a small growth and a constant
velocity difference :

Ar2 = Q&collCv

where C, is the volume fraction of the fog, coll the av-
erage collision efficiency and a results from the integra-
tion 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 sim-
ulations are converged for a rough number of 100 itera-
tions and provide a self-similar solution which means
low diffusion in the size phase space.

M-onodisperse injector

o --- tunnel (5m)


| Fog machine (r, Ci, ,

Tranquilizing Grid


7 i~8a

n Langmuir-Blodgett
25 O Beard-Grover


. 15

0 10 20 30 40
Fog concentration C, (ppm)

50 60

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

3.3 Discussion on velocity

As a second step, let us revise the constant velocity as-
sumption. 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 momen-
tum 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 momen-
tum transfer are enabled in Fig. 9. The droplet velocity
tends to be slightly smaller than the terminal velocity
because of momentum dilution.

2 25




0 1 2 3 4 5
Position (m)
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)

Yet the consequences of velocity on coalescence are
not obvious. First let us neglect efficiency laws. The co-
alescence rate per unit time linearly increases with the

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

velocity difference (eq. 2) but linearly decreases with
the time spent in the tunnel so that it vanishes in the fi-
nal growth formula (eq. 11). In fact radius growth for a
fixed length travel only depends on the number of colli-
sions i.e. fog concentration and travel length, assuming
a fairly constant cross section. However when consider-
ing 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 aver-
age velocity might have been slightly overestimated in
D'Herbigny and Villedieu (2001). As a final restric-
tion, please note that no coalescence efficiency law has
been used (coal 1) which suggests not to conclude on
which collision efficiency law to use.

o Langmuir-Blodgett
25 Beard-Grover

S 10 20 30 40
Fog concentration C, (ppm)

50 60

Figure 10: Radius growth Ar2 with velocity resolu-
tion (filled blue : Second order simulation;
black: Experimental)

3.4 Conclusion on the model and methods

The previous results confirm that Eulerian Multi-Fluid
methods can be used to simulate accurately the size-
distribution evolution of a coalescing spray and its size-
conditioned dynamics. 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 re-
sults 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 correspond-
ing laws can significantly reduce the particle growth rate
which is experimentally observed so that these efficiency
models are essential to capture the physics of coales-

Y s-L~~

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

4 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 in-
jection of gas and small alumina droplets resulting from
the combustion of propergol and an axial two-phase flow
carrying bigger droplets supposedly coming from pri-
orly injected, burnt and coalesced particles. The upper
part of the chamber generating these bigger droplets is
not solved.

4.1 The CEDRE code

The CEDRE code has been developed at ONERA for
several years to provide a multiphysics comprehensive
tool for energetic problems in aeronautics. The main
feature is the combination of a non structured spatial
solver with models as varied as multiphase flows, multi-
species chemistry, thermal conduction, radiation, wall-
film models, etc. The coupling can be one or two-
way. As for the disperse phase, CEDRE includes a La-
grangian 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 co-
alescence models.

4.2 Configuration

The 2D rocket booster simulation is performed on a sim-
plified configuration, yet featuring the main difficulties
of solid propulsion typical flows i.e. parietal injection
and a supersonic nozzle. It takes place on a deformed-
structured 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 to-
tal rate of flow of 10 kg.s 1.m 2 from the wall (pro-
pellant combustion) and the head end (upstream flow),
this gas flow involves extremely high velocity gradients
in the nozzle (Fig. 11, middle), which will induce size-
conditioned droplet dynamics as studied in section two.
The simulation strictly speaking starts when injecting
the disperse phase : a monomodal 5om radius wall in-
jection represents the droplets resulting from recently
burnt aluminum particles, directly expelled from the pro-
pellant grain, and another roughly monomodal distri-
bution around 20/m is injected on the axis to model
the previously coalesced droplets (Fabignon et al. 2003).
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 p 1 ,liL'wNi

wall injection (5p/m)


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

4.3 Results

For the disperse phase, we choose a 5 section discretiza-
tion (0, 12.5, 25, 37.5 and 50 pm) 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 frac-
tions 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. Spe-
cific 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 ignilik.lI vol-
ume fractions ( 10 5) in the chamber. Third, bigger
droplets are created in the nozzle eventhough accelera-
tion 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.

4.4 Remaining Numerical difficulties in solid

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 re-
sults however suggest several possible improvements,
among which reducing computational time and consid-
ering more physical phenomena.
Typical advanced rocket booster simulations are un-
steady and performed on more complex meshes : they
can include from 10, 000 cells to millions of cells for

I Soctinn 1 [0 19 lii m I

Figure 12: Dispersed phase volume fractions

3D meshes. Highly optimized code and multiprocessor
computations (with several spatial domains for instance)
are therefore required. In our constant time step sim-
ulation, 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 in-
duce droplet secondary break-up which strongly under-
mines droplet growth. Since it is out of the scope of this
paper, this phenomenon has been neglected but its mod-
elling reaches two requirements : specific impulse com-
putation is subject to the precise knowledge of droplet
size after the nozzle throat and experimental data on alu-
mina size distribution mainly consist of ejected material
which has encountered complete break-up process. De-
veloppement of fragmentation models for the Eulerian
Multi-Fluid method has been achieved in Dufour (2005).
Finally, the main drawback of such methods comes
from the monokinetic ]ihp ,liLsis : 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 sym-
metry axis. On more advanced solid propulsion un-
steady configurations, dedicated to hydrodynamic insta-
bility studies for instance, high-inertia droplets ejected
from vortices need a multivelocity treatment which ad-
vanced high order moment methods can provide as de-
veloped in Kah et al. (2010).

1 0x10
1 4x10'
5 2x10
1 9x10
7 2x10(
2 7x10
1 0x10
1 4x10
5 2x10
1 9x10
7 2x10

1 4x10
7 19x10
10x lO
1 9x10
7 2x10

SSection 2: [12.5, 3-]i |.

Section 3: [25, 37.5]/um

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

5 Conclusions

In this paper, we provide a comprehensive validation of
a Eulerian, second order in size, method for solving two-
phase 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 Eule-
rian Multi-Fluid method including validated advanced
collision efficiency models, which are crucial for rocket
booster studies. Finally, the method is practically effi-
cient in rocket booster simulation contexts; moreover re-
ducing the number of required sections can be profitable
to the computational time of such simulations. The re-
maining limitations underlined in the final section re-
quire though further developments but the method can
bear such improvements.

6 Acknowledgments

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


B. Abramzon and W.A. Sirignano. Droplet vaporization
model for spray combustion calculations. Int. J. Heat
Mass Transfer, 32:1605-1618, 1989.

P. Achim. Simulation de collisions, coalescence et rup-
ture de gouttes par une approche lagrangienne: appli-
cation aux moteurs a propergol solide. PhD thesis, Uni-
versit6 de Rouen, 1999.

A. A. Amsden, P. J. O'Rourke, and T. D. Butler. Kiva ii,
a computer program for chemically reactive flows with
sprays. Technical Report LA-11560-MS, Los Alamos
National Laboratory, 1989.

K. V. Beard and S. N. Grover. Numerical collision ef-
ficiencies for small raindrops colliding with micron size
particles. J. of atmospheric sciences, 31:543-550, 1974.

G. A. Bird. Molecular gas dynamics and the direct sim-
ulation of gas flows. Oxford Science Publications, 42,

F Bouchut. On zero pressure gas dynamics. In Advances
in kinetic theory and computing, pages 171-190. World
Sci. Publishing, River Edge, NJ, 1994.

I Q- -rtinn -, f,,m I

PR. Brazier-Smith, S.G. Jennings, and J. Latham. The
interaction falling water drops : coalescence. Proceed-
ings of the Royal Society, 326:393-408, 1972.

S. De Chaisemartin. Moddles euleriens et simulation
numerique de la dispersion turbulente de brouillards qui
s'evaporent. PhD thesis, Ecole Centrale de Paris, 03

F. X. D'Herbigny and P. Villedieu. Etude exp6rimentale
et num6rique pour la validation d'un module de coales-
cence. Technical Report RF1/05166 DMAE, ONERA,

G. Dufour. Modelisation multi-fluide eulerienne pour les
ecoulements diphasiques a inclusions dispersees. PhD
thesis, Universit6 Paul Sabatier Toulouse III, 2005.

J. K. Dukowicz. A particle-fluid numerical model for
liquid sprays. J. Comput. Phys., 35(2):229-253, 1980.

J. Dupays, Y. Fabignon, P. Villedieu, G. Lavergne, and
J. L. Estivalezes. Some aspects of two-phase flows
in solid-propellant rocket motors. In Solid Propellant
Chemistry, Combustion, and Motor Interior Ballistics,
volume 185 of Progress in Astronautics and Aeronau-
tics. AIAA, 2000.

Y. Fabignon, O. Orlandi, J.F. Trubert, D. Lambert, and
J. Dupays. Combustion of aluminum particles in solid
rocket motors. AIAA Paper 2003-4807, July 20-23 2003.
In 39th AIAA/ASME/SAE/ASEE Joint Propulsion Con-
ference and Exhibit, Huntsville, Tx.

J. B. Greenberg, I. Silverman, and Y. Tambour. On the
origin of spray sectional conservation equations. Com-
bustion and Flame, 93, 1993.

J. J. Hylkema. Modelisation cinetique et simulation
numerique d'un brouillard dense de gouttelettes. Appli-
cation auxpropulseurs a poudre. PhD thesis, Ecole Nat.
Sup6rieure de l'A6ronautique et de l'Espace, 1999.

J. J. Hylkema and P. Villedieu. A random particle
method to simulate coalescence phenomena in dense liq-
uid sprays. In Lecture Notes in Physics, volume 515,
pages 488-493, Arcachon, France, 1998. Proc. 16th Int.
Conf. on Num. Meth. in Fluid Dyn.

D Kah, F Laurent, L Freret, S De Chaisemartin, R
Fox, J Reveillon, and M Massot. Eulerian quadrature-
based moment models for dilute polydisperse evap-
orating sprays. Flow Turbulence and Combustion,
pages 1-26, 04 2010.

P. Kuentzmann. Aerothermochimie des suspensions.
Gautier-Villars Editeur, 1973.

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

I. Langmuir. The production of rain by a chain reaction
in cumulous clouds at temperatures above freezing. J.
Meteor, 5:17-192, 1948.

F. Laurent. Numerical analysis of Eulerian multi-fluid
models in the context of kinetic formulations for dilute
evaporating sprays. M2AN Math. Model. Numer. Anal.,
40(3):431-468, 2006.

F. Laurent and M. Massot. Multi-fluid modeling of lam-
inar poly-dispersed spray flames: origin, assumptions
and comparison of the sectional and sampling methods.
Comb. Theory and Modelling, 5:537-572, 2001.

F. Laurent, M. Massot, and P. Villedieu. Eulerian multi-
fluid modeling for the numerical simulation of coales-
cence in polydisperse dense liquid sprays. J. Comp.
Phys., 194:505-543, 2004.

M. Massot. Eulerian Multi-fluid modelsforpolydisperse
evaporating sprays, chapter III of "Multi-Phase React-
ing Flows: Modelling and Simulation", pages 79-123.
CISM Courses and Lectures No. 492. Springer, 2007.
Editors D.L. Marchisio and R. O. Fox,.

M. Massot, F. Laurent, D. Kah, and S. De Chaise-
martin. A robust moment method for evaluation of the
disappearance rate of evaporating sprays. to appear
in SIAM J. of App. Math., 2010. http://hal.archives-

P. O'Rourke. Collective drop ifh. i on vaporizing liquid
sprays. PhD thesis, Los Alamos National Laboratory
87545, University of Princeton, 1981.

M. Rilger, S. Hohmann, M. Sommerfeld, and
G. Kohnen. Euler/Lagrange calculations of turbulent
sprays: the effect of droplet collisions and coalescence.
Atomization and Sprays, 10(1):47-81, 2000.

B. T6th. Two-phase flow investigation in a cold-gas
solid rocket motor model through the study of the slag
accumulation process. PhD thesis, Universit6 Libre de
Bruxelles, 2008.

I. M. Vasenin, R. K. Narimanov, A. A. Glazunov, N. E.
Kuvshinov, and V. A. Ivanov. Two-phase flows in the
nozzles of solid rocket motors. J. Propulsion and Power,
11(4):583-592, 1995.

F A. Williams. Spray combustion and atomization.
Phys. Fluids, 1:541-545, 1958.

F A. Williams. Combustion Theory (Combustion Sci-
ence and Engineering Series). ed F A Williams (Read-
ing, MA: Addison-Wesley), 1985.

Y. B. Zel'dovich. Gravitational instability : an approxi-
mate theory for large density perturbations. Astronomy
and Astrophysics, 5:84-89, 1970.

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs