EFFECT OF WATER DROPLETS ON EXTINCTION 177

plane several times, before being completely
vaporized. The most important finding of the
present work is that the flame extinction strain
rate shows nonmonotonic behavior for different
monodisperse droplet sizes considered for sev-
eral water mass loadings (1-3% in condensed
phase). Assuming that the droplets leave the air
stream at the same velocity as the gas, 15-um
droplets were scen to be the most effective. For
example, addition of 3% of water by mass in
condensed phase (in addition to 2.22% as satu-
rated vapor) was shown to reduce the extinction
rate to 134 s~', from about 400 s~! for the dry
case. The nonmonotonic effect of droplet size,
especially the superior flame extinction phe-
nomenon associated with 15-um droplets, is
attributed to the droplet dynamics in the coun-
terflow field and to the large mass evaporation
and heat sink observed near the oxygen con-
sumption layer. A careful analysis of the physi-
cal, thermal, and chemical contributions con-
trolling the phenomena predicted here is
currently under way and will be presented in a
separate paper.

Sensitivity analyses of the approximate ana-
lytical models used to simplify the numerical
integrations show that mass evaporation and
heat transfer terms need to be more accurately
modeled. Methods of improving the accuracy of
these models will be pursued in the future.

This work is supported by National Institute of
Standards and Technology, Gaithersburg, MD,
under the Grant No. 5D0098, with Dr. W. L.
Grosshandler as the technical monitor.

APPENDIX A

The Lagrangian equations describing the con-
servation of mass, momentum (both x and r
directions). and enthalpy of a single droplet can
be written as

dmy

di (Al)
d
Fatal = Ova t Fy (A2)
d
i (mtg) = ~ Qua + F, (A3)

& (nda) = -Qlhy +L) +H (A4)
where, m, is the droplet mass, Q the mass
evaporation rate, F, the axial component of the
external force acting on the droplet, F, the
radial component of the external force acting on
the droplet, f, the droplet specific enthalpy, H
the heat flux from the gas to the droplet, and L.
the average heat of evaporation at the droplet
surface.

The values of Q, F,, F,, and H depend in
general on the state of the droplet, its relative
motion with respect to the gas, the thermo-
chemical and transport properties of the gas,
and the presence of other droplets in the neigh-
borhood. To accurately determine these terms,
the regions around each droplet controlling the
mass and heat transfer should be resolved.
Fortunately, water vapor mass fraction, temper-
ature, velocity, and thermodynamic coefficients
of the gas around the droplet vary over dis-
tances greater than few diameters. Thus, aver-
age values of these gas-phase variables are used
to determine Q, F,, F,, and H.

Simplified Models for Evaporation, Heat
Transfer, and Forces Acting on the Droplet

When droplets are not at rest relative to the
surrounding gas phase, convection affects evap-
oration and heat transfer. This effect can be
assessed by estimating the droplet Reynolds
number defined as

(AS)

where p is the density averaged in the region
around the droplet, jt the viscosity averaged in
the region around the droplet, Vy, = (vy — v)?
+ (uy ~ u)?)'?, and d the diameter of the
droplet. When Re, > 20, there exists a bound-
ary-layer flow region around the front of the
droplet and a wake region behind it [11], and
therefore the droplets can interfere with cach
other. When Re, is of the order of 100 or
higher, the shear stress at the gas-liquid inter-
face can be large enough to induce internal
liquid-phase circulation [10, 11]. However, if the


EFFECT OF WATER DROPLETS ON EXTINCTION 171

suming thermal equilibrium condition given by
Eq. 24 (solid line), whereas other droplets show
a small thermal lag. However, this lag has
insignificant cffect on the flame structure and
extinction condition considering the small tem-
perature difference and the energy associated
with it. The variation of flux fraction ¥ as a
function of the axial location, shown in Fig. 8,
however, has a very important role on the
droplet source terms described by Eq. 17, there-
fore, the state of the flame. The reduction in ¥
seen as the droplet approaches the stagnation
plane is because of the flow divergence or
straining effect. If the droplet vaporization is
neglected, the curves shown for different drop-
let sizes shift somewhat, but the dramatic reduc-
tion of ¥ seen with droplet location does not
change. Similar results for axial liquid methanol
volume flux were observed experimentally by Li
et al. [5], with and without a flame.

SOURCE TERMS WITHOUT
INTERACTION

Once the solution of droplet variables is ob-
tained (i.¢., d, vy, Uy, Ty, and F), the source
terms contributing to the gas phase can be
evaluated through Eq. 17. Here, such contribu-
tions are shown, once again by turning “off” the
interaction of such source terms with the gas-
phase calculations to keep the flame structure
condition the same for different droplet sizes
considered.

Figures 9-11 show the variation of mass,
momentum, and energy source terms, i.e., S,,,,
Sy — US,,, and S, ~ h,S,,, along the axis of
symmetry for different droplet sizes. Figure 9
shows rather strange profiles, especially on the
air side before the flame, where evaporation
effects are negligible. It is scen that smaller
droplets exert a greater radial force per unit
volume on the gas than larger droplets. This is
however easily explained by the fact that larger
droplets shoot through that region faster than
smaller droplets, so that they have much shorter
residence time. A comparison of Figs. 10 and 11
indicates the profiles of S,, and S, — h,S,, to be
very similar (except for the change of the sign).
This suggests that heat transfer to the droplets
and evaporation are well correlated; the impli-

150-02

120-02

900-03

6.00-03

Sm (g/cm)

3.00-03

0.06400
-9.00-09 Lo dg
“0 0200 02 04 08

z (cm)

Fig. 9. Typical S,,, profiles of different droplet sizes, with Yo
= 0.02 and droplet source terms turned “off” in gas-phase
calculations. Also shown is the gas temperature (thick linc).

cations of such correlations are discussed fur-
ther in the next section.

EFFECT OF WATER DROPLETS ON
FLAME EXTINCTION

The most interesting question that can be
probed by the two-phase model developed here
is the effect of water droplets on flame extinc-
tion condition, when the source terms are in-
cluded in the gas-phase solution, in particular,
the role of water droplet size on the flame
extinction condition for a given mass loading.

900-01

806-01

706-01

600-01
$.00-01
~4.00-01

ODL

300-01
206-01

Sy =U x Sp (gm/em'/s)

~1.00-01

0106400

rey
100-0155 0.2 00 02 O4 06

Fig. 10. Typical Sy — US,, profiles of different droplet
sizes, with Y, = 0.02 and droplet source terms turned “off”
in gas-phase calculations. Also shown is the gas temperature
(thick line).



162

production by chemical reactions of species k,
h, the specific enthalpy of the species k, cp, the
specific heat at constant pressure of the species
k, V, the axial mass diffusion velocity of species
k, W, the molecular weight of species k, W the
average gas molecular weight, R° the universal
constant of gases, i the index of the evaporating
species, and 6, the Kronecker delta. The term
J = (1/r)(ap/ar) in Eq. 3 cannot be neglected
as its limiting value at the x-axis is not zero and,
in fact, is a constant [6, 12, 13].

The flow velocities, chemical compositions,
temperature, and pressure at the nozzle exits
are assumed given (as controllable or measur-
able quantities). Therefore, if x, and x, are the
axial positions of the lower and upper nozzles
shown in Fig. 1, the boundary conditions can be
written as

v(x) = Vp Vu) = Mu

U(x) = Up Ul%y) = Uns Yel) = Yay
¥,(x,) =¥y,k=1,...,N
T(x) = Ty Ty) = Tu

P(x) = Py) = Po (7)

Unless otherwise mentioned, the values used
here are T; = T,, = 300 K, py = 1 atm, x, =
—0.4 cm, x, = 0.6 cm, and U,, = U, = 0. The
composition of fuel stream is pure methane, so
that You, = 1, but that of air can be different,
depending on the water vapor or droplets used.
The values of v, and v, are varied to get
different flame extinction conditions. Because
only the first derivative of v appear in the
equation of mass, v cannot satisfy both the mass
equations and the two B.C.’s. Thus, the term J
acts as a constant pressure parameter, or an
cigen value, that allows the mass equation to be
satisfied, once v is specified at the boundaries.
The solution of the system of equations de-
scribed above is exact on the x-axis and is
approximate for small values of r compared to
the radius of the air nozzle. With the reduced
number of discretized equations, finite-rate
chemistry effects can be resolved accurately
using a detailed reaction mechanism in a highly
refined numerical grid along the axis of symme-

try.

A.M. LENTATI AND H. K. CHELLIAH

Singularities of Droplet Number Density
Equation

The dynamics of a single droplet can be deter-
mined by integrating in time its equations of
mass, momentum, and energy. Previous studies
[2, 3] have suggested that condensed-phase ef-
fects can be coupled to the gas phase by simply
multiplying the interaction of one droplet by the
local droplet number density, ny, that should
satisfy a steady-state equation. However, in the
case of the counterflow field, difficulties arise as
droplets can penetrate the stagnation plane and
reverse their trajectories [4]. If Eulerian formu-
lation is used for droplets, with self-similarity
formulation, the equation describing the con-
servation of droplet number density, ng, can be
written as (assuming droplets are not created or
destroyed)

d
(naa) + Inga = 0, @)

where Uy = Ug/rgo, with vq and uy being the
axial and radial velocity components of the
droplet and ryo the initial radial location of the
droplet (very close to the axis of symmetry) at
the air nozzle exit. The above equation assumes
that, for any point in the domain where the
equation is regular, n, and U, can be uniquely
defined. Assuming U, is assigned or already
calculated from single droplet equations, the
above equation indicates that dny/dx blows up
as vq — 0, and therefore ny.

Such a high peak in 7, is not expected in the
droplet flow, where the average distance be-
tween droplets can vary from about 20 to 35 or
more droplet diameters and their motion is
quite ordered. In fact, ng cannot even be strictly
defined as the number of droplets per unit
volume because the sample volumes are layers
perpendicular to the x-axis and their thickness
can be of the same order as the droplet diame-
ter. Under such conditions, as the droplet mo-
tion cannot be treated as a continuum, the
problem needs to be reformulated in a different
way, i.e., by using a Lagrangian description of
the droplet motion as described below.


178

droplets are very small, as considered here, such
effects can be neglected.

In the counterflow under consideration, the
estimated Re, is of the order of 1 or less (at
Jeast near the x-axis) so that the approximation
of droplets in quiescent atmosphere can be
adopted. To get to a simplified and explicit
analytical solution for Q and H, the following is
assumed: (1) the temperature and diffusion
fields are spherical-symmetric, (2) pressure is
constant, equal to po, (3) the droplets are
sufficiently far apart such that their evaporation
and heat transfer occurs as in a single isolated
droplet, (4) boundary conditions at infinity of
the surrounding gas of each single droplet are
given by the gas phase T and Y,’s at the droplet
location, (5) evaporation and heat transfer are
quasisteady (i.e., the droplet diameter, temper-
ature, etc. are assumed to vary slowly relative to
the transients of the gas properties around the
droplet), (6) the gas surrounding the droplet is
treated as a binary mixture of the evaporating
species and a gas whose thermodynamic prop-
erties are the same as a uniform mixture of the
other nonevaporating species present in the gas
phase at the droplet location, and (7) chemical
reactions in the regions around each droplet are
neglected (i.e., the evaporated species affects

the chemical reaction in the gas phase only
through the source term in Eq. 3). Integration
of equations for mass, diffusion, and enthalpy of
such single isolated droplets yield the following
expressions for mass evaporation and heat
transfer:

(A6)
a _ Qc, /21dd
H = dad\T.. » Texpt Q¢,/2mdX) - 1)
(AT)

where subscript i indicates the evaporating spe-
cies (in this case water), and s and © indicate the
conditions at the droplet surface and ~, respec-
tively. The gas properties with overbar are eval-
uated using the standard 1/3 law [10]. Here, Te;

= (Alpt, Dim is the Lewis number based on
averaged. properties, with Dj, being the overall
diffusivity of the evaporating species into the
mixture.

A. M. LENTATI AND H. K. CHELLIAH

The force components F, and F, include
viscous drag forces due to droplet motion rela-
tive to the gas phase and evaporation rate.
Gravity effects are neglected here because of
the high convective velocities. As Re, is of the
order of 1 or less, the force components are
approximated assuming Stokes drag coefficient

‘1p = 24/Rey, which yields

F, = 3adjulv — va) (A8)

F, = 3ndjitu — uy) (A9)

The equations above allow the existence of a
self-similar solution for the single droplet equa-
tions.

Self-Similar Formulation of Droplet Dynamics

For the gas-phase flow field written in the form
of Eq. 1, when Q, H, F,, and F, are given by
Eqs. A6-A9, the right-hand sides of Eqs. Al,
A2, and A4 depend only on xy, the droplet axial
coordinate, whereas the right-hand side of Eq.
A3 is proportional to ry, the droplet radial
coordinate. In Ref. [2], a new variable s =
rglug, Was introduced to cast the system of
equations in a self-similar form, with sy = 1/a
= constant for the potential flow field consid-
ered with a defined as the flow strain rate and
subscript 0 identifying the conditions at the air
inflow plane. For the general rotational flow
formulation considered here, with the flow
strain rate varying from zero at the air inflow
plane (i.e., for plug flow) to some finite value
near the mixing layer, assuming

ug = Udi)ra (A10)

the system (A1)-(A4) can be reduced to a
system of ODE’s independent from r4, so that
the trajectories of all the droplets can be deter-
mined with a single integration in time. With
these simplifications, the system of Eqs. Al-A4
together with the equations for the droplet
location can be written as

dxy _
Waa (All)

rg = ra Rt) (A12)


EFFECT OF WATER DROPLETS ON EXTINCTION 175

than 7,,. Thus, the assumption of uniform tem-
perature field within the droplet is expected to
be invalid. However, the error introduced in the
solution is not large because the fraction of
energy flux that goes into raising the tempera-
ture of the droplet is typically of an order of
magnitude less than that associated with latent
heat of evaporation, especially for the size of
droplets considered here. In addition, small
error in 7, (the droplet temperature at the
surface for nonuniform droplet temperature) is
expected to introduce negligible influence on
the mass evaporation rate as T,, => T,. This can
be better demonstrated by writing the mass
evaporation rate in terms of the transfer num-
ber based thermal driving force, Q', given by [8,

10]
roa yk GAT. — T,)
Q = and Saft 4 TST ] (31)

Therefore, the profiles of the source terms of
mass and enthalpy are hardly affected by the
above errors. Certainly, the comparison of char-
acteristic time scales indicates the need for
more-accurate modeling of the internal heat
transfer, especially for larger droplets such as
those encountered in water sprinklers.

Sensitivity of the Solutions to F,, Q, and H

Besides the deficiencies in droplet evaporation
and heat transfer modeling described above, the
estimation of the body force F, may involve
some uncertainties, primarily associated with
the neglect of gravity and low Reynolds number
approximation for the viscous drag force. For
the high-strain rate flows considered here, it is
relatively easy to establish the negligible role of
gravity based on the ratio of potential energy to
kinetic energy, given by 2g(x, — x)/vao ~
107°.

According to Refs. [2, 10, and 11] the models
developed assuming quiescent atmospheric con-
ditions underestimate evaporation and heat
transfer by ~25% and drag coefficients by
~15% for large droplets when Re, ~ 1. How-
ever, no mention is made about possible effects
of mass evaporation on drag or its interaction
with convection effects at high Re,. For the
conditions considered here, however, the maxi-

TABLE 3

The Sensitivity of @,,, to Variation of F,, Q, and H for
d= 30 pm, Yy = 3%

F= F143) 0=Q0 +a) H=HO4 A)

(%) %) OW) mac)
—20 —0.9 +87
—8.5,

+20 +41

mum observed Re, is about 2.5 ford = 50 um
and ¥, = 1% near extinction, whereas, for all
the other cases, it is equal to or less than 1. It
should also be pointed out that high Re, effects
cannot be directly included in the current mod-
eling approach, as Q, H, and F, would end up
depending on r, and the source terms them-
selves on r, invalidating the hypothesis used in
defining the self-similar solution. One way to
overcome such a difficulty is by assuming
Rey, > Reg, where Rey and Reg, are Reynolds
numbers based on the axial and radial velocity
components of the droplet, respectively. Then
F,, F,, Q, and H are modified by factors of the
kind 1 + a@Ref, where a and B are constants
<1 [4]. Such an assumption may be considered
valid only near the x-axis, when the gas-phase
flow is given as in Eq. 1.

Before undertaking the challenging task of
improving the submodels used, the sensitivities
of any possible errors associated with F,, Q, and
H are assessed here. For the case with d = 30
pm and Y) = 2%, the modification of F,, Q,
and H by +20% on the predicted extinction
strain rate is shown in Table 3. It is evident that
large errors on the axial force estimations do
not particularly affect a,,,, at least as long as the
force is strong enough to make the droplets
follow the gas well. The small increase in a,,, for
positive AF, is due to the fact that near the
flame the gas accelerates, dragging the droplets
with itself before slowing down again at the
stagnation flame as can be seen in Fig. 4. Errors
in both Q and H appear to affect the solution
more significantly. In the case of modified Q,
the lower variation on a,,, for A = +20% is
perhaps due to the coupled effect of Q on F,
and H, which reduces the droplet size rapidly to
the 15-um range. In case of modified H, the
variation of a@,,, is the same for both signs of A.
The a,,, is more sensitive to H not only because


EFFECT OF WATER DROPLETS ON EXTINCTION 167

Calculation of Source Terms

Once d, vy, Uy, Ty, and ¥ are determined from
Eqs. 18-23, the source terms arc first calculated
from Eq. 17 at the coarse and uniform grid
points x,,. If both the droplet locations x”) and
xy"! (given by Eq. 22) are in the same Ax,
centered i in x,, (ie. the droplet center of mass
remains in the same volume), then all the
contributions are assigned to the pointx,;, while
on all the other points x,, (k # j) the source
terms are not modified. Alternatively, x7*!
must be in the volume centered in one of xy" 's
two neighbors because of the small time steps
taken. In the latter case, the contributions to the
source terms are then shared by the two grid
points proportionately based on the fraction of
time 8f,, spent by the droplet in each of the
corresponding volumes. Subsequently, the val-
ues of S,,, Syj;, and Sij used at the points x;
G =i, » N, + 1) in the Eulerian formula-
tion are calculated interpolating from the values
given above,

The numerical approximations introduced
above in evaluating the terms S,,,;, Sy;, and S,,;
certainly can have some errors associated with
them. For example, the first-order temporal
integration scheme using the values of gas-
phase properties interpolated on a coarse grid
can be a source for such errors. However, the
procedure adopted here gives acceptable nu-
merical results when the gas-phase properties
that affect the droplet dynamics change on a
scale equal to or greater than Ax,, and the time
step 61, is taken small enough. Also, when a
droplet is crossing the boundary between two
volumes, the contributions that are not more
rigorously given to one of the volumes when
droplet mass center is, in fact, located in the
adjacent cell can be another source of an error.
Such errors are however minimized by taking
small time steps. Because a rigorous analysis of
the error is not possible due to the nonlinearity
of the problem, an attempt to quantify the
errors associated with the hypothesis intro-
duced, by using a sensitivity analysis, is dis-
cussed later in the paper.

Solution Algorithm

The equations describing the gas phase with the
droplet source terms can be written in the form

F(s, 8) =0 (25)

where s represents he matrix s;; (i = 1,
N+4,j=1,...,N, +1) of the (N +4) x
Wp + 1) unknowns , Yu, (kK = 1,..., .N),
Uj, vj, and J;, (j = 1,0... N, + 1) and § the
ect of the source terms Sinjs Sup and S,,;,
G = 1,...,.N, + 1). With a “good” initial
guess for s, these equations are integrated with
an algorithm that uses quasi-Newton steps.
After n + 1 steps, the approximation of the
solution is given by

"F(s", (1 — yy8" + y8""!)
(26)

stl = gt — (yl!

where s" and §" are the solution and droplet
source terms S after n steps, respectively, and
(J|s')" is the inverse of the Jacobian calculated
at step 2 keeping the source terms constant. The
Jacobian is in block-tridiagonal form and is
calculated numerically and then inverted, call-
ing a subroutine already used in previous work
[15-17]. Here, a relaxation parameter y < 1 is
introduced to speed up the convergence to
steady-state solution.

As shown in Eq. 17, the source terms depend
on the gas-phase values at the points that are
used to integrate the single droplet equations.
However, such dependence is not accounted in
the Jacobian, otherwise, the matrix would lose
its property of being block tridiagonal. This
means that J|g is only an approximation of the
real and less-sparse Jacobian. In general, con-
vergence is quite fast (no more than five or six
iterations), but it can become a problem when
Sj Su; and S,, peak values are large or
become very sensitive to the solution. For this
reason, in Eq. 26, F is sometimes calculated
introducing the relaxation factor y with $"~', to
compensate for very quick variation of S,,;, Syj,
and S,,; at each step.

The solution is assumed to have converged
when

srt! — 5

wow, (27)

>t 1)*
where § is the vector of the unknowns scaled by
some reference value of the solution variables,
and €q is a small prefixed value. If the solution
converges, a grid refinement step is performed


164

Fig. 2. An illustration showing the effect of flow strain on
areas formed by water droplet streams as a function of time.

droplets crossing a unit surface area perpendic-
ular to the x-axis at time ¢ and the number of
droplets crossing an equivalent area at time fy.
Because of the self-similarity approximation,
Eq. Al2 shows that, if the droplets are uni-
formly distributed at time ta, F(t) is indepen-
dent from the radial distance r, as the distance
from the x-axis of each droplet increases in time
with the same factor R(t). Therefore, F(t) may
be also defined as the ratio between Ay and A,,
as all the droplets crossing A, also cross A,, and
their distribution remains uniform over each
cross section. The number of droplets crossing
section Ag is NyqVyo. It is only near the turning
point that a droplet number density cannot be
satisfactorily defined. Therefore, the number of
droplets crossing a unitary surface perpendicu-
lar to the x-axis at time ¢ is ny vyoF(t). By
considering a differential volume, in the limit
time step dt — 0, it can be shown that the
function ¥ satisfies the following equation

dF

an = 2 FUa (14)
(Note: the above equation can also be obtained
by integrating Eq. 11 to find A,.) At time ty, ¥
= 1 by definition. Once the cquations for a
single droplet, described by Eqs. 9-12 and 14,
are integrated in time, values of my, Uy, Ay, and

A. M. LENTATI AND H. K. CHELLIAH,

| d<d,
i d

I
|
l
I t
|
|
!

yl

Fig. 3. An illustration of the axial motion of droplets in
time. for two extreme droplet sizes.

F are known along the entire droplet path, i.e.,
from the point on Ay to complete evaporation.
This information is next used to determine the
source terms contributing to each node in Eu-
lerian formulation of the gas-phase equations.

Consider the case where the counterflow re-
gion is divided into fixed layers whose thickness
is at least a few times larger than the droplet
diameter. As illustrated in Fig. 3, a droplet may
cross the boundaries of any volume dV; a few
times or go to complete vaporization in this
volume element. Such phenomena are very
likely to happen if the volume dV; is close to the
stagnation plane. Each droplet that moves in-
side dV; transfers mass, momentum, and en-
thalpy to the gas in a time interval dr by the
amount

8b = —b(0dt (15)

where ¢ represents the right-hand side terms in
Eqs. A13, A15, and Al6 (ic, -O, —QUy +
3adji(U — Uy) and —Q(hy + L) + H).
Because the number of droplets that cross the
plane Ay per unit area per unit time is nyVaos
the number of droplets reaching the element
dV, and contributing to Eq. 15 per unit area and
unit time is nav goF(t). Assuming that the path
of each droplet in dV; is divided into intervals
{ty, t, + dt}, k = 1,..., Kj, the mass,
momentum, and enthalpy transfer in the vol-
ume dx,dA are given by



EFFECT OF WATER DROPLETS ON EXTINCTION 163

Conservation Equations for the Condensed
Phase

The Lagrangian equations that describe the
conservation of mass, momentum, and enthalpy
of a single droplet in a steady counterflow field
are considered here. The droplet number den-
sity considered is small (corresponding to a
maximum water mass fraction of 3% in con-
densed phase), so that the average distance
between droplets is at least an order of magni-
tude greater than their diameters. Hence, drop-
let-droplet interactions, especially collisions,
can be neglected. The droplet is assumed to be
spherical and nonrotating. For simplicity, the
temperature of the droplet is assumed to be
uniform, but is allowed to vary in time from the
room temperature at the nozzle exit to the
equilibrium temperature determined by the
heat transfer and evaporation rate (see Eq. 24).

To integrate such equations describing the
droplet motion, simplified models are needed to
describe the viscous forces through Stokes drag
formulation and heat and mass transfer through
Spalding’s transfer number concepts (see Ap-
pendix A). The boundary conditions necessary
for such lumped models are averaged gas veloc-
ities, temperature, water vapor mass fraction,
and transport coefficients. All such variables are
assumed to vary along the flow stream lines with
characteristic lengths greater than the droplet
diameter. Under such assumptions, a self-simi-
lar solution for the droplet dynamics can also be
used. Introducing the subscript d to identify the
conditions of the droplet, the appropriate dy-
namical variables characterizing the state of the
droplet can be defined by velocity components
vy and Ug (= uy/rgo), Mass my (or diameter d,
knowing the density of water, p,), and temper-
ature 7,. The resulting Lagrangian equations
for mass, momentum, and energy to be solved
are (see Appendix A for further details)

d 8av
= (d?) = ~~~ 9
a” PL p ®)
dvg_ 18h)
dt pa Vd) (10)
dU, = —//2 18a =
Gen Ua 5 pl ~ Ua (11)

dTy__12A_[(T2= TY LY -
exp) -1 a, | ©)

dt Pac, >
where is the Spalding’s transfer number given

by

we (; - ¥) 13
=—In 3

Le; 1-Y,, (13)

The solution of the above equations yields the
transfer of mass, momentum, and energy from a
single droplet to the gas phase along the droplet
trajectory until it is completely vaporized. The
details of coupling the source terms arising from
a single droplet to the source terms appearing in
gas phase Eqs. 2-5 are described next.

Source Terms S,,, Sy, and S,

In principle, the source terms S,,, Sy, and S,,
contributing to the gas-phase equations can be
determined by integrating Eq. 8 for the droplet
number density, ny, and then multiplying Q,
QU, — Fira, and Q(hg + L) — H (ie,
negative of the right-hand side of Eqs. A13,
A15, and A16 by mq at each location along the
x-axis. As discussed earlier, the solution of such
a Eulerian equation for ny diverges if the drop-
let velocity approaches zero exactly at any node.
Consequently, the source terms will also blow
up when multiplied by , to include multiple
droplets in the flow. Instead of integrating the
droplet number density equation given by Eq. 8,
if a Lagrangian equation for droplet flux is
integrated in time, then the source terms con-
tributing to equations for mass, momentum, and
enthalpy will remain finite even at the droplet
turning point. The development of an appropri-
ate Lagrangian equation for droplet flux frac-
tion is described below.

In Appendix A, it is shown that, under some
simplifying but reasonable assumptions, my, Vy,
U4, and T,, are not functions of r,,, at least near
the axis of symmetry. This means that droplets
crossing the inlet section A, at the same time ty
then cross section A, together, at time f,, A, at
time fy, etc., with same mg, Vg, Ug, and Ty, as
shown in Fig. 2. Taking advantage of this prop-
erty, in the present calculation of source terms,
a droplet flux-fraction function F(t) is intro-
duced as the ratio between the number of


176
06 1.0
04 08
02 06
A o
w
0.0 04
droplet trajectory in time
2
0.2 0.2
“Bbevoo 100-02 200-02 308-02 400-02 5. 0-38
Time (s)

Fig. 18. Droplet trajectory as a function of time, ford = 50
pm and Y, = 0.03.

of the modification of Q by changing the equi-
librium temperature but also because of the
rapid increase in T,, before the droplet reaches
the flame.

Droplet-Droplet Collisions

A key assumption in the numerical model de-
veloped is that droplets are transported without
collision with others. This assumption allows for
a fairly simple description of the droplet flow
field and derivation of a self-similar solution for
both droplets and gas phase. It relies on the idea
that the droplets leaving the air nozzle are
monodisperse with the same exit velocity (vy),
hence, similar dynamics in the counterflow field.
Although such ideal conditions are difficult to
be realized in experimental situations, if the
initial droplet separation distance to droplet
diameter is large, as seen by Table 1, the
fraction of droplets colliding before evaporation
can be neglected.

The worst situation that can violate this as-
sumption arises when large droplets with high
droplet loading penetrate the stagnation plane
and reverse their trajectory, creating a region
where opposed flow of droplets exists near the
axis of symmetry. Such a region exists for 50-4m
droplets considered here, as seen in Fig. 5.
Figure 18 shows the calculated trajectory of
such 50-4m droplets along the x-axis as a func-
tion of time and the corresponding flux fraction

A. M. LENTATI AND H. K. CHELLIAH

function for Yy = 3%. With vj ~ 60 cm/s and
an average distance /, = 24d = 0.12 cm (sce
Table 1), each droplet is estimated to be pre-
ceded and followed by a droplet before and
after a time 1.8 X 107*s. Thus, if it is assumed
that droplets are Icaving the air nozzle at every
1.8 x 1074s, then the region between points |
and 2 (the region where each droplet group
would see a group of droplets flowing in the
opposing direction) in Fig. 18 consists of ~17
layers. By taking the average value of droplet
flux fraction, ¥, in region 1 —> 2 as 1/2 and in
region 2 — 3 as 1/4, the probability of droplet
collision in this worst case region can be roughly
estimated as

P. <1731(2) 0.4%
coll 24\,) 0%

ly 62)

Clearly, this estimate for collision probability
gets even smaller for smaller droplets because
of evaporation and the increased straining of
droplet flux. Hence, the assumption of negligi-
ble droplet collision is reasonable.

SUMMARY AND CONCLUDING REMARKS

The purpose of this work was to develop a new
two-phase numerical model to analyze the ef-
fectiveness of monodisperse condensed-phase
fire suppressants in extinguishing counterflow
flames. An important feature of the new model
developed is it can overcome the singularities
associated with solving the droplet number den-
sity equation in counterflow field accurately and
in a very robust manner. This was achieved by
using a hybrid Eulerian—Lagrangian formula-
tion for the gas—droplet flow, with the introduc-
tion of droplet flux fraction to describe the
droplet number density equation.

As part of this work, the model developed
was applied to a case where monodisperse water
droplets are introduced with the air stream to a
nonpremixed counterflow flame of methane and
air. Several cases with different droplet sizes,
ranging from 5-50 xm, were considered. Small
droplets (<20 jxm) tend to follow the gas closely
and go through rapid evaporation in the hot
mixing layer, and it was shown that they never
cross the stagnation plane. The 50-4m droplet
considered was shown to cross the stagnation


166

enough to allow analysis of different flow con-
ditions in a reasonable time and (3) most of all
robust, as the equations must be solved for
droplets of different sizes and different flow
field conditions. Because the variables control-
ling droplet equations vary slowly in space as
compared to some radical species or chemical
reaction rates, in the procedure employed here,
a new set of coarse grid points x,(j = 1,...,
N, + 1) is used for the integration of droplet
equations. To analyze the grid independence of
the solution, the coarse grid chosen consists of
equidistant points at intervals Ax, greater by
some factor than the initial droplet diameter.
On the selected coarse grid points, the values of
v, ¥,, U, and T are linearly interpolated. To be
consistent, these values are then used to deter-
mine the gas-phase properties using the
Chemkin and transport libraries [18, 19] as well
as liquid and liquid—-vapor equilibrium proper-
ties from thermodynamic table data [21].

The integration time step is chosen each time
so that the droplet will not travel more than a
fraction of the coarse grid interval. Denoting
each quantity calculated at time 7, with the
superscript 7, and at time ¢,,, = ft, + 6t, with
the superscript n + 1, the single droplet equa-
tions are approximated in the following form:

@eyr!

er __ situ

-—— (18)
Btn Pap
mt yn an
var tava BR ne
upd OM) (19)
Ue) = US myer 4 182"
uO Ua + pad)"
“(U"— UG") (20)
Ty) ~ Ty 12a"
ot, Pik pay"
(TL- Tv" LL"
1 — exp('¥”) (21)
Pe
¢ i, 2 = 0.5tvy + vy"! (22)
get _ gn .
ae UY + UG") (23)

5,

n

A.M. LENTATI AND H. K. CHELLIAH

The gas-phase properties appearing in these
equations are calculated only at time ¢,,, includ-
ing v and U, as implicit schemes can be very
time consuming. The terms in the right-hand
side of Eqs. 18-23 could be calculated at time
t,,, yielding a simple explicit, first-order method.
However, nonphysical oscillations might appear
in vy, Uy, and & (the instability of Ty is
discussed separately below) in case 8t,, is not
small enough, which is the case when mm, ap-
proaches zero (or d? — 0).

To improve the stability of the ordinary dif-
ferential equation (ODE) solver, vy, Uy, and ¥
appearing on the right-hand sides of Eqs. 18-23
are evaluated at time ¢,,, ;, with the exception of
the term U3, which is not fully implicit, ie.,
utu%*! instead of (U%*')?, to keep the equa-
tion in Uy*! linear. Also vy and U, in Eqs. 22
and 23 appear as averages between values at t,,
and f,,,, to improve the accuracy. The method
used in Eqs. 18-23, although not fully explicit or
fully implicit, is consistent (converges to the
ODE as 6t,, — 0), first order, and absolutely
stable. The error introduced can be made arbi-
trarily small by reducing 61, until it becomes
smaller than the error introduced calculating
the source terms on the coarse uniform grid
x(j = 1,.-., Ny + 1) instead of on the
original grid x,(j = 1,...,.N, + 1) where the
Eulerian equations are discretized.

As the droplet enters the hot mixing layer, its
temperature increases rapidly, and the numeri-
cal approximation for 7%*' can overshoot the
value determined by the thermal equilibrium,
ie., balance between the heat transfer to the
droplet and heat of evaporation. This equilib-
rium temperature of the droplet, Tj, where
these two terms balance is given by

L"
(Tz, ~ Tae) = = (1 — exp(¥"))
op

(24)

To eliminate spurious oscillations associated
with overshoots of droplet temperature, when
T,*' given by Eq. 21 is greater than T,,, the
condition T*' = Ty, is imposed. This essen-
tially implies that, as the characteristic time of
Tq becomes much smaller than &t,,, the asymp-
totic solution of thermal equilibrium is justifi-
able and is used here.


Dynamics of Water Droplets in a Counterflow Field and their
Effect on Flame Extinction

A. M. LENTATI and H. K. CHELLIAH*
Department of Mechanical, Aerospace and Nuclear Engincering, University of Virginia,
Charlotesville, VA 22903

‘The effect of fine-water droplets in extinguishing steady, laminar counterflow methane—air nonpremixed flames
is investigated here, using a numerical approach. A new two-phase model using a hybrid Eulerian-Lagrangian
formulation for the gas-droplet flow is developed as part of this work. A key feature of the model developed
is that it can avoid the singularity associated with the droplet number density equation in a consistent manner
by using a Lagrangian equation for droplet flux fraction. The gas phase is described by a detailed model
involving full chemical kinetics and transport, whereas droplet evaporation and heat transfer arc modeled
assuming quasisteady conditions. Application of the model to several monodisperse sizes of water droplets,
ranging from 5~50 ym, revealed an interesting nonmonotonic dependence of the flame extinction strain rate on
droplet size. This phenomenon is attributed to the droplet dynamics in the counterflow field considered here
and to the resulting nonmonotonic heat sink associated with mass evaporation observed at the oxygen-

consumption or radical production layer of the flame.

INTRODUCTION

Automatic water sprinklers, with a typical mass
median drop size of about 750 wm, are exten-
sively used as fire-suppressing devices in certain
types of fires (e.g., building fires). Recently,
however, fine-water droplets or water-mist sys-
tems with drop sizes less than 100 ym are being
considered as a replacement for recently
banned halon 1301 (bromotrifluoromethane—
CF3;Br, one of the most effective chemical
agents in suppressing gaseous flames but known
to cause ozone depletion in the upper atmo-
sphere [1]). One major advantage of a water-
mist system over a traditional sprinkler system is
that fine-water droplets with lower settling ve-
locity can suspend in air for much longer peri-
ods than larger droplets, hence, its applicability
as a gaseous fire suppressant. In addition, ther-
mal concepts are often used to distinguish the
effectiveness between the two systems. For ex-
ample, for the same mass of water, the latent
heat of vaporization is much greater than the
increase in sensible enthalpy from room tem-
perature to the boiling point. Thus, fine-water
droplets are believed to be more effective than
larger droplets because of the large surface area
to volume ratio leading to rapid evaporation
and faster cooling of the hot gases. However, in

*Corresponding author.

0010-2180/98/$19.00
PIT S0010-2180(97)00355-6

© 1998 by The Combustion Institute

flows associated with real fires, the two classes
of droplets can have very different droplet dy-
namics, and the above qualitative physical and
thermal reasoning may not hold. The urgency
for rapid transition from halon 1301 to alternate
systems, such as water-mist systems, has forced
researchers to rush into applications without
detailed understanding of the relative (1) phys-
ical cffects due to the water droplet size, (2)
thermal effects due to heat capacity and latent
heat of evaporation, and (3) chemical effects
due to enhanced overall three-body recombina-
tion reactions and shift in water-gas reactions.
The objective of this work is to understand the
role of such physical, thermal, and chemical
effects and accurately quantify their contribu-
tions using a simplified reacting flow field.
Laboratory-scale laminar flames established
within the mixing layer of counterflow of fuel
and air provide an excellent configuration to
analyze the dynamics of such water droplets in a
hot reacting environment and to understand the
flame extinction mechanism. Except for the
effects of turbulence (which can introduce many
undesired uncertainties in modelling transport
phenomena and turbulence-chemistry interac-
tions), this nonpremixed flow configuration has
many important features that are relevant for
understanding the extinction mechanism of real
fires. The most important feature is the ability
to accurately control the flow residence time to

COMBUSTION AND FLAME \15:158-179 (1998)
© 1998 by The Combustion Institute
Published by Elsevier Science Inc.


EFFECT OF WATER DROPLETS ON EXTINCTION 169

TABLE 1

Comparison of the Number Density and Droplet
Separation Lengths for Some Selected Water Mass

Loadings
a Ny ny hy

(um) (emTM*) (cmTM') (cm)

5 178408 56.3 0.0177

20 2787 14.1 O07

50 178 5.6 OAT78

2% 5 300457 2 0.0141

ld ~ 28 20 5632 178 0.0562

50 360 7 0.1406

3% 5 546259 81.7 00,0122

Nd ~= 24 20 8535 20.4 0.0485

50)

8.2 0.1223

of the same order. This implies that the neglect
of droplet-droplet interaction is reasonable for
all the cases considered here. It should be
pointed out that in actual water-mist systems, it
is highly unlikely that the water droplets gener-
ated are monodisperse. However, under very
low water mass flow rates, the generation of
monodisperse water droplets in laboratory ex-
periments is well established and is being pur-
sued by several groups, including our group at
the University of Virginia. The predictions per-
formed here with monodisperse droplets pro-
vides a better mechanism to analyze and under-
stand the basic droplet dynamics and
interactions with the gas-phase processes, which
can be easily extended to polydisperse droplet
flows in the future.

As mentioned previously for dilutely loaded
droplets. the gas displacement by water droplets
can be neglected because of the large density
ratio between the condensed water and the gas.
Under such conditions, the only interaction
between the droplets and the gas is because of
the source terms S,,,, Sy, and S,, appearing in
Eqs. 2-5. When these terms are included in the
numerical calculations, the flame structure and
extinction strain rate are affected dramatically,
depending on the amount of condensed water
added and the droplet sizes. Thus, for the
purpose of illustrating and comparing the tra-
jectories of various droplet sizes in the counter-
flow field and thcir associated source terms, a
fixed low-strain rate case (a = 130 s~') where
the droplet source terms are excluded (or the
source terms turned “off” in numerics) is con-

vu (m/s)

Lo —_
ary) 02 0.0 02 04 06
x (em)

Fig. 5. Comparison of the gas velocity and droplet velocity
of different sizes, with droplet source terms turned “off in
gas-phase calculations. Also shown is the gas temperature
(thick line).

sidered first. It should be mentioned here that
the profiles of droplet trajectories and source
terms presented below remain essentially the
same for high-strain rates as well. Figure 5
shows a comparison of the axial velocity along
the axis of symmetry of the gas and droplets of
different sizes for the selected low-strain rate
case. Although the coupling arising from the
droplet source terms on gas-phase structure is
neglected, in solving droplet equations for drop-
let trajectory shown in Fig. 5, the droplet vapor-
ization effects are included in Eqs. 18-23. In
Fig. 5S, the 5-zm droplets are seen to follow the
gas fairly closely, whereas large droplets deviate
considerably. The droplet lag scen becomes
slightly worse for high-strain rates, especially
near extinction conditions. The 5-m droplets
are also scen to be completely vaporized soon
after they enter the hot mixing layer, whereas
50-4m droplets penctrate through the flame
and also the stagnation plane and then reverse
their direction somewhere in the fuel stream. If
the droplet vaporization is suppressed, then the
penetration of 50-j4m droplets is much more
pronounced, with multiple crossings at the stag-
nation plane. These predictions are consistent
with those reported carlier by Chen et al. [4].
The droplet response to the variation in gas-
phase flow conditions can be characterized by
the Stokes number, S/. By defining the flow
residence time in the outer air flow as t, = Va,



168

2050

2000 oO.
a
cs

1950 *
x

1900

Trnaz (K)
a

Lo
15005 160 200 250-300 350 <400~—450

Strain Rate, a (s~?)

Fig. 4. Comparison of the maximum flame temperature,
Tyars ¥S. flaw strain rate, a, for different inflow water vapor

partial pressures.

to determine whether and where new grid
points must be added based on the evaluated
gradients and curvatures, and the iteration pro-
cess is repeated [17].

EFFECT OF WATER VAPOR ON FLAME
EXTINCTION

In a real-life fire situation, or in laboratory
experiments, when water droplets are intro-
duced to suppress flames, they are bound to
saturate the room air with water vapor. Conse-
quently, the fraction of oxygen displaced (or
diluted) by the water vapor can reduce the
flame strength, or the extinction strain rate
(@..,), Of the counterflow flame considered
here. Figure 4. for example, shows a comparison
of the predicted variation of maximum flame
temperature, T,,,,,, aS a function of the flow
strain rate (a), for some selected water vapor
concentrations. The line with symbol O corre-
sponds to a dry case, whereas the symbol *
correspond to a fully saturated case at a room
temperature of 300 K and atmospheric pressure
(i.e., partial pressure of water Py,0,sa: = 0.0351
atm or, in mass fractions, Yij.0,4 = 0.0222, with
Youu = 0.2278 and Yx,,, = 0.7500). The fuel
considered in these calculations is pure meth-
ane (Ycu,, = 1.0). Based on this figure, if the

A.M. LENTATI AND H. K. CHELLIAH

asymptotic value of highest strain rate for each
case is taken as the extinction strain rate, @,,,,
then the results indicate that saturated water
vapor alone can decrease the flame extinction
strain rate by 25%. This enhanced effect is not
solely due to dilution, but is partly duc to the
high thermal capacity of water.

DYNAMICS OF WATER DROPLETS
WITHOUT INTERACTION

If the amount of water being injected into the
air stream is greater than Pyy.o,sa, = 9.0351 atm
(or Yj1,9 > 0.0222), then the excess water will
be in condensed form, ic., in the form of
droplets. The mass fraction of water associated
with such droplets leaving the air nozzle can be
defined using

MaMa ao
Maa + Pura

where 7, is the initial mass of the droplet, 149
the initial number density, vj) the initial velocity
of the droplet, and p,, and v,, are the density and
velocity of air leaving the nozzle, respectively.
All the results shown here are obtained by
assuming that the droplet velocity leaving the
air nozzle is same as that of the gas, i.e. vyg =
v,» but other inflow conditions can be easily
analyzed with the model developed. One of the
major objectives of the present work is to un-
derstand the trajectories of condensed-phase,
monodisperse water droplets in a reacting coun-
terflow field considered here and investigate
their effect on the flame extinction condition.
In the present predictions, the effect of vari-
ous monodisperse droplet sizes for three differ-
ent fixed water mass fraction in the condensed
phase (Y, = 1, 2, and 3%) is considered, to
make fair comparisons. Table | shows the num-
ber density of water droplets leaving the air
nozzle, ny (#/cm’), for three different Y, and
for three different monodisperse droplet sizes
selected. The last column in Table 1 also lists
the characteristic separation length between
droplets, /, (cm). It is interesting to note that,
for all the combinations of water mass fraction
in condensed phase and droplet sizes, the ratio
of droplet separation length to the droplet
diameter, /,/d, remains fairly large (>20) and is

(28)

Yo=


EFFECT OF WATER DROPLETS ON EXTINCTION 165

K; tetdiy
Ab = fra > | waa
iy,

kat di,

(16)
In terms of per unit volume, as in Eqs. 2-5, Eq.
16 can be written as

Ky furan
Sy = — (dx nuva, 2 | F(thpdt
1

kal dy

(17)

In Eq. 17, the droplet mass center is assumed
to be in dV,. For finite sized droplets, when
their mass centers are near the boundary of
dV, part of the droplet surface is in the adja-
cent volume, so that the contributions to the
source terms should be shared accordingly.
However, the correction to Eq. 17 may be
considered negligibly small in case the droplet
has a diameter much smaller than dx; or the
source term contributions vary only slightly in a
distance of the order of the droplet diameter. In
case the latter condition is not satisfied, the
lumped description of the gas—droplet interac-
tion described in this work is not possible.

It should be pointed out that, in the present
formulation, the droplets leaving the air nozzle
are assumed be monodisperse. In actual exper-
iments, achieving a monodisperse droplet load-
ing is a challenging task. Although experiments
with water droplets have not been performed
yet at the present time, some droplet diameter
distribution is expected around a mean value.
The formulation described above, however, can
be easily generalized by approximating such a
continuous distribution of droplets as a collec-
tion of discrete sections with average droplet
diameters d; (1 = /, /). The single droplet
equations then can be solved, and Eq. 17 can be
evaluated for each section d; by replacing yy
with finy), Where f, is the fraction of droplets
having diameter d,. In the present study, the
above generalization is not undertaken until the
experimental results become available.

NUMERICAL APPROACH
Integration of Gas-Phase Equations

Equations 2-5 describe a set of nonlinear ordi-
nary differential equations for N + 4 un-

knowns, namely, T, ¥,(k = 1,..., N), U, v,
and J. Numerical methods have been used in
the past to obtain accurate solutions for the
system of equations, without the source terms
Siny Sy, and S$, [15-17]. These involve use of
first-order finite difference methods with
pseudo-Newton algorithm steps. Although
other higher-order methods are available, the
method used is unconditionally stable and re-
quires just a three-point stencil. Unconditional
stability is of paramount importance in the
solution of chemically reacting gases, where the
presence of sharp variations in mass fraction
profiles may lead to nonphysical solutions or no
convergence at all. The three-point stencil sim-
plifies the inversion of the Jacobian in the
pscudo-Newton steps, as the matrix turns out to
be a block tridiagonal, for which stable solution
codes are available. Convection terms are ap-
proximated with an upwind scheme. In particu-
lar, in calculating the enthalpy transport in Eq.
5, the upwind scheme is applied using the total
mass average velocity of each species k, by v +
V,. Adaptive grid refinement is also performed
after the solution converges. This allows the
first-order method to yield accurate results with-
out using too many grid points. For cases where
uniform grids are used, typical grid spacing is
about Ax; ~ 10 ym.

In solving Eqs. 2-5, the thermochemical and
transport properties, namely, p, W,, X,, @,
Coes Nyy Daglks @ = 1,..-,N), A, and p, are
calculated at each point x; using the Chemkin
and transport libraries [18, 19]. The chemical
reaction rates are modeled by 17 species in 39
elementary reaction steps [20]. For low Mach
number flow considered here, pressure py is
assumed given. Integration of droplet equations
and determination of the source terms S,,,, Sy,
and S,, are discussed below.

Integration of Droplet Equations

The single droplet equations must be solved at
the end of each pseudo—Newton step to update
the source term S,,,,, S,;, and S,,, contributing to
Egs. 2-5. These are performed at every iteration
step so that no droplet property needs to be
stored. The algorithm used for their solution
must be (1) accurate enough to produce error

not greater than the modeling itself, (2) fast



EFFECT OF WATER DROPLETS ON EXTINCTION 179

d

a cat
d -

gy (iva) = ~Qva + mdjlv vy) (ANA)

d
qi oma + mgU3, = —QUy + 3djpiu(U — Uy)
(A15)

S (ndg) = Olly +L) H (A16)

where R(t) = exp(fi, Uy(4dt) is a function of

hy

t only. The initial conditions are

Xalty) = Xu» Pato) = Tay

valto

dito

= Vay Ualto

(A17)

where 7, is the droplet temperature, assumed
to be uniform,

Equations A13-A16 can be further simplified
by introducing an averaged droplet specific
heat, c,g, and liquid density, py (independent of
T4), so that

ha = Cp, Ta — To) + hay

dy, Ta(To) = Ta,

(A18)

ma = F pat? (A19)
where /yg, is the droplet specific enthalpy at the
reference temperature Ty. Then, use of Eqs.
A6, A7, A18, and A19 in Eqs. A13-A16 yields
the simplified set of Eqs. 9-13 described previ-
ously.

REFERENCES

1. Grosshandler, W. L., Gann, R. G., and Pitts, W. M.,
(1994) Publ. SP-86/, National Institute of Standards
and Technology, Gaithersburg, MD.

2. Continillo, G., and Sirignano, W. A., Combust. Flame
81:325 (1990).

3. Lacas, F., Darabiha, N., Versaevel, P., Rolon, J. C.,
and Candel, S., (1992) Twenty-Fourth Symposium (In-
ternational) on Combustion, The Combustion Institute,
Pittsburgh, p. 1523.

4. Chen, N.-H., Rogg, B., and Bray, K. N. C., (1992)
Twenty-Fourth Symposium (International) on Combus-
tion, The Combustion Institute, Pittsburgh, p. 1513.

5. Li, $. C,, Libby, P. A., and Williams, F. A, Combust.
Flame 94:161 (1993).

6. Chelliah, H. K., Law, C. K., Ueda, T., Smooke, M. D.,
and Williams, F. A., (1990) Twenty-Third Symposium
(Intemational) on Combustion, The Combustion Insti-
tute, Pittsburgh, p. 503.

7. Seshadri, K., Combust. Flame 33:197 (1978).

8. Williams, F. A., Combustion Theory, 2nd ed., Addison
Wesley, Reading, MA, 1985.

9. Kuo, K. K., Principles of Combustion, Wiley-Inter-
science, New York, 1986.

10. Law, C. K,, Prog. Energy Combust. Sci. 8:171 (1982).

11. Abramzon, B., and Sirignano, W. A., Int. J. Heat Mass
Transfer 32(9):1605 (1989).

12, Kee, R. J., Miller, J. A., Evans, G. H., and Dixon-
Lewis, G., (1988) Twenty-Second Symposium (Interna
tional) on Combustion, The Combustion Institute,
Pittsburgh, p. 1479.

13. Smooke, M. D., Crump, J., Seshadri, K., and Giovan-
gigli, V., (1990) Rep. ME-100-90, Yale University, New
Haven, CT.

14, Seshadri, K., and Williams, F. A. Int. J. Heat Mass
Transfer 21:251 (1978).

15, Smooke, M. D., J. Comp. Phys. 48:72 (1982).

16. Miller, J. A, Kee, R. J., Smooke, M. D., and Grear,
J. F,, (1984) The Computation of the Structure and
Extinction Limit of a Methane-Air Stagnation Point
Diffusion Flame, The Combustion Institute, Pitts-
burgh:

17, Smooke, M. D., Puri, I. K., and Seshadri, K., (1986)
Twenty-First Symposium (International) on Combus-
tion, The Combustion Institute, Pittsburgh, p. 1783.

18. Kee, R. J., Rupley, F. M., and Miller, J. A., (1993)
Chemkin Il: A Fortran Chemical Kinetics Package for
the Analysis of Gas Phase Chemical Kinetics, Rep.
SAND89-8009B, Sandia National Laboratories, Liver-
more, CA.

19. Kee, R. J., Dixon-Lewis, G., Warnatz, J., Coltrin,
M. E., and Miller, J. A., (1992) A Fortran Computer
Code Package for the Evaluation of Gas-Phase Multi-
component Transport Properties, Rep. 86-8246, Sandia
National Laboratories, Livermore, CA.

20. Peters, N., and Rogg, B., Lecture Notes in Physics, Vol.
ml5, Springer-Verlag, New York, 1993.

21. Vargaftik, N. B., Tables on the Thermophysical Proper-
ties of Liquids and Gases, Hemisphere, New York,
1975.

22. Seshadri, K., and Peters, N., Combust. Flame 73:23
(1988).

23. Chelliah, H. K., and Williams, F. A., Combust. Flame
60:17 (1990).

24. Peters, N., Combust. Sci. Technol. 30:1 (1983).

25. Kim, J. S., and Williams, F. A., SIAM J. Appl. Math
53(6):1551 (1993).

Received 20 August 1997; revised 7 November 1997; accepted
21 November 1997



EFFECT OF WATER DROPLETS ON EXTINCTION 161

Simplified Conservation Equations for the
Gas Phase

The gas phase is described by a set of Eulerian
equations of conservation of mass, momentum,
and energy. The flow is assumed to be axisym-
mctric, i.e., no swirl, and steady. Because of high
flow strain rates considered here, the effects of
body forces both in the momentum and in the
energy equations are neglected. The gas is a
Newtonian fluid, but, due to low Mach number,
bulk viscosity effects are neglected. Diffusion
velocities are approximated using a simplified
multicomponent model. Viscous heating, ki-
netic energy, work due to gas compression,
radiative heat transfer, and Dufour effects are
neglected as compared to convective, diffusive,
and heat release terms. The condensed phase is
assumed to be monocomponent, so that only
one of the species equation is affected by its
evaporation. The condensed phase effects are
averaged as mass, momentum, and enthalpy
transfer from the condensed phase to the gas
phase per unit volume.

Because the flame established in the mixing
layer is planar, a self-similar solution that pre-
serves profiles of scalar variables along the
radial direction is used. Thus, we seek a solution
of the form [6, 12, 13]

ve=v(x),u =rU(x), Yy = ¥x(x),
T = T(x), p = plx.r) Q)

which becomes exact along the axis of symme-
try. Here, x is the axial coordinate, r the radial
coordinate, 1 the mass averaged gas velocity in x
direction, w the mass averaged gas velocity in
radial direction. Y, the mass fraction of species
k, T the temperature, and p the pressure. For
low Mach number flows considered here, pres-
sure is assumed to vary slowly along the axis of
symmetry, so that the gas density depends only
on the chemical composition and temperature,
and the momentum equation along such axis
can be eliminated. However, the dependence of
pressure in the radial direction cannot be ne-
glected a priori and is given by p = py(x) +
r°ép/ar? [14]. Hence, the equation for radial
momentum described below must be retained.

In the absence of droplets, such self-similar
formulations have been widely adopted in inves-

ligating counterflow gaseous flame structures
and flame extinction conditions, with compari-
sons to experimental data obtained along the
axis of symmetry. For low relative velocity of
droplets with respect to the gas velocity, i.c.,
small-droplet Reynolds number, it can be as-
sumed that the source terms S,,, S,,/r = Sy, and
S,, contributing to the gas-phase conservation
equations are functions of axial coordinate x, as
shown in Ref. [2] for outer potential flow and in
Appendix A for rotational flow with plug flow
boundary conditions. Here, S,, is the mass
transfer per unit volume from the droplets to
the gas, S,, the radial momentum transfer per
unit volume from the droplets to the gas, and S,,
the enthalpy transfer per unit volume from the
droplets to the gas.

With these source terms, the conservation
equations for mass, momentum in r direction,
species, and energy can be cast in a quasi—one-
dimensional form. These equations, together
with the equation of state specifying the gas
density, are listed below

d
Gi, (ov) + 2pU = Sy, @

du, d{ du
wt grass A (aie)
+Su-US,, (3)

dy, d
pv Tet ay (PVM a) = Wats
+ (Bu ~ YOSuv
k=1,...,N (4)
arid ( ar
m dx dx \" dx

N
pd wt hd Ye
kat

N
— Dd Wyergdty + Sy,

—hSiy (5)
Po
RWT (6)

where p is the gas density, » the gas viscosity, »
the gas thermal conductivity, @, the molar


172

4.00408

-2.00+08

x Sm (erg/em*)

1.00408

Sno

0.00400

a dy
rr | Ly
(cm)

1.00108
“0.

IL. Typical S,, — 4,8, profiles of different droplet
sizes, with ¥, = 0.02 and droplet source terms turned “off”
in gas-phase calculations. Also shown is the gas temperature
(thick line).

Such effects can be shown by plotting the vari-
ation of maximum flame temperature as a func-
tion of the flow strain rate, as shown in Fig. 12,
similar to that shown earlier in Fig. 4 for water
vapor effect. In addition to the saturated water
vapor at room conditions (solid line), when
condensed water 1% by mass is introduced, i.e.,
Y,) = 0.01, in the form of different droplct
sizes, different profiles for T,,,,. are obtained
and are plotted in Fig. 12 for comparison. This
figure indicates that 50-4m droplets are the
least effective with a relatively higher extinction

2000

1950

1900

Taz (K)

1650, —
100 150 200 250 300 350

Strain Rate, a (s7!)

Fig. 12. Comparison of T,. VS. @ for different droplet
sizes, with Y, = 0.01.

A.M. LENTATI AND H. K. CHELLIAH

2000 —
(oS eal vaporip_v=psa
Jo ad= Smcwns
1950 |& ~ Ad = 15 microns
fi “hd = 90 mirons
8 d= 50 microns
1900
*,
& 1850 a.
3 &
1800 ‘
1 .
1750 ‘
iN 4 \
4y
1700
1650 a —s
100 150 200 250 300 350

Strain Rate, a (s~!)

Fig. 13. Comparison of Tyg. VS. d, for different droplet
sizes, with Yy = 0.02.

strain rate (still below the saturated vapor case
with water vapor mass fraction of 0.022), while
15-~m droplets are the most effective with the
lowest extinction strain rate. As the size of
monodisperse droplets is decreased below 15
pm, they are seen to be less effective. The
lowest size considered here is 5 jzm, as the mass
and enthalpy source terms approach a 5-func-
tion as seen in Figs. 9 and 11, introducing severe
convergence problems of the numerical model.
In fact, droplets of 5 wm and below should
approach the limit where all the water consid-
ered is in vapor phase at the air side boundary
of the mixing layer, provided the inflow temper-
ature of air is not affected. Irrespective of this
difficulty, the most interesting nonmonotonic
behavior on droplet size is captured well by the
two-phase model developed here for steady,
laminar, counterflow, nonpremixed flames.

For higher water mass fractions in condensed
phase, similar profiles for 7,,,,. aS a function of
flow strain rate are obtained, as shown in
Figs. 13 and 14. However, the extinction strain
rate is seen to reduce significantly with increas-
ing water droplet loading. For various mass
fractions considered here, the nonmonotonic
effect of droplet size on the extinction strain can
be better illustrated as shown in Fig. 15. For all
the conditions considered here, 15- to 20-j4m
droplets are seen to be the most effective. This
supcrior effectiveness of 15- to 20-um droplets


174

cating a trend similar to that of a,,, shown in
Fig. 15.

The results presented here need to be vali-
dated with careful experimentation, and such
studies are currently under way. However, the
generation of truly monodisperse droplets is a
challenging task, and methods of extending the
present work to polydisperse droplets need to
be pursued in the future.

SENSITIVITY OF THE RESULTS TO THE
APPROXIMATIONS INTRODUCED

Several simplifications and approximations
were used in the development of an analytical
model for droplet equations. The implication of
some of these approximations introduced and
possible inaccuracies arising from these approx-
imations are addressed here separately.

Grid and Time Step Size

In solving the set of ODE’s described by Eqs.
2-5, an initial solution on a coarser grid is
traditionally used, with subsequent refinements
of the grid where the gradients and curvature
are large. This method works well when the
outer flow is approximated by a potential flow
field where U is independent of x in the outer
frozen flow region [16, 17]. However, for the
rotational flow field used here, with plug flow
boundary conditions, the grid refinement
scheme does not add additional points to refine
the outer flow region where the velocity gradi-
ent changes with x. To obtain accurate solution
for the flow strain rate, a, based on the velocity
gradient of the outer flow in the oxidizer side,
the domain outside the mixing layer is artifi-
cially refined until the strain rate solution be-
came independent of the grid selected. As men-
tioned previously, the droplet equations are
solved on a coarser grid x,; such that Ax, is at
least two times the largest droplet size consid-
ered, i.e., 100 wm. However, the final solution is
found to be invariant with the coarser grid
selection as long as the source terms contribut-
ing to each cell are appropriately partitioned,
depending on the time spent on each cell. In
addition to the grid selections, tests performed
with the Lagrangian time step modified with a

A.M. LENTATI AND H. K. CHELLIAH

4.00408

8.00408
‘g
7B 2.00108 1500
= 4g
f | z
% 1.06408 1000
<<
1
a

0.00400 }——-_——. ——. ee 500

1.004
on 86 4 02 00 02 O4 0. e

2 (cm)

Fig. 17. Typical energy source term profiles with different
values of f, for d = 5 ym and Y, = 0.02.

factor f are shown in Fig. 17. This figure indi-
cates no change in the flame structure and very
minor changes in S,,, with the f values selected.

Quasi-Steady Droplet Evaporation and Heat
Transfer Model

The expressions for Q and H described in
Appendix A are derived under the assumptions
that T, is uniform in the droplet and that the
quasisteady-state solutions of the surrounding
gas can accurately approximate the water vapor
concentration and temperature profiles. The
above relies on the hypothesis that the charac-
teristic times of mass diffusion tp and thermal
conduction in the droplet and in air, r,, and
Trgass Tespectively, are much smaller than the
characteristic flow residence time of the droplet,
ya, which depends on how fast the droplet is
traveling along the x-axis. Using the following
expressions for above time scales, tj = Ax/
Vamarr t.= PID, Tha = Pap Ral ar ANd Treas
= pc,d?/X, where overlined properties identify
averaged quantitics around the droplet (see
Appendix A). Here, Ax is a characteristic length
along the x-axis, where gas-phase temperature
and vapor mass fraction vary appreciably, A, is
the thermal conductivity of liquid water, and
Va,max iS the maximum value of velocity of the
droplet in its trajectory.

Estimation of the characteristic time scales
described above from the converged numerical
solution indicates that 7, is definitely larger


160
OXIDIZER + DROPLETS
| | NITROGEN
T ~
STAGNATION
PLANE
NITROGEN

Fig. 1. An illustration of the counterflow field with water
droplets. Also shown is the typical gas flow velocity profile,
and the definition of flow strain rate used.

transport data, The gas—droplet interaction in
the present model is described using hybrid
Eulerian—Lagrangian formulation, somewhat
similar to that employed in Refs. [2-4]. The
discretized gas-phase equations are solved using
an algorithm that include Newton and fixed-
point iteration steps. Unlike the previous work
(4], a rigorous method of handling the singulari-
ties associated with the droplet number density
equation at the stagnation point of droplets is
described. The resulting model is then used to
predict the trajectory of water droplets and the
modified extinction conditions for a range of
monodisperse droplet sizes.

FLOW CONFIGURATION

A schematic of the counterflow configuration
considered in the present study, very similar to
the ongoing experimental effort, is shown in Fig.
1. When ignited, a steady laminar flame can be
established within the mixing layer, as indicated
in Fig. 1. The exit velocity profiles of both air
and fuel streams are assumed to be of the plug
flow type. The droplets are introduced though
the upper nozzle at a steady rate with the air
stream and are assumed to be monodisperse. In
the present simulations, the droplet velocities at

A. M. LENTATI AND H. K. CHELLIAH,

the nozzle exit are assumed to be the same as
that of air. In actual experiments, these condi-
tions may differ slightly from those assumed
here but can be casily implemented with the
two-phase model developed here.

For the two-phase system considered here,
the gas velocity can differ from that of droplets,
The response of droplet trajectories for differ-
ent flow strain rates can be addressed by varying
the nozzle cxit velocities of methane and air
streams. For plug flow boundary conditions
used here (corresponding to viscous, rotational
flow formulation), the flow strain rate, a, is
defined as the negative of the axial velocity
gradient of the gas in the air side (i.c., —dv/dx),
just outside the mixing layer, as illustrated in
Fig. 1. At flame extinction, the critical value of
the flow strain rate is identified here as a,,,. In
experiments, the outer air flow is known to
deviate slightly from the plug flow boundary
conditions [6], but, once detailed laser Doppler
velocimetry data becomes available, such
boundary conditions can be easily implemented
in the present model.

FORMULATION

The Navier-Stokes equations, together with
multicomponent species and energy conserva-
tion equations, are well known and can be found
in any text on combustion [8, 9]. With the
presence of droplets in the flow, these equations
can become rather complicated, but somewhat
simplified two-phase conservation equations
have already been considered in the context of
spray combustion [2-4, 10, 11].

Because of the large density ratio between
the liquid and the gas and the low-number
density of droplets in the flow, effects of the
volume occupied by droplets on the convective
and diffusive terms can be easily neglected, with
errors introduced being less than 0.1%. As a
consequence, interactions between the droplets
and gas are only due to source terms associated
with transfer of mass, momentum, and energy.
Analyses of such couplings can be considerably
simplified by assuming that the droplets are
monodisperse, as considered here,


EFFECT OF WATER DROPLETS ON EXTINCTION 159

investigate the finite-rate chemistry effects. For
this reason alone, numerous theoretical and
experimental studies have been performed in
counterflow field containing homogeneous re-
actants and, more recently, heterogeneous reac-
tants in the context of fuel spray combustion.
For the first time, Continillo and Sirignano [2]
investigated the combustion of monodisperse
sized fuel droplets in a counterflow field numer-
ically for three different flow conditions. They
also investigated the conditions for validity of
the self-similarity approximation generally in-
troduced in solving the conservation equations
of the counterflow field. In this initial study of
Continillo and Sirignano [2], many limiting ap-
proximations on the chemical kinetic model,
thermochemical and transport data were used.
Hence, no attempt was made to predict the
flame extinction phenomena. In a subsequent
study, Lacas et al. [3] performed detailed nu-
merical simulations of a similar monodisperse
spray containing (1) n-heptane fuel droplets,
flowing against air, and (2) liquid oxygen drop-
lets, flowing against gaseous hydrogen. In both
these studies [2, 3], the droplets were trans-
ported in a carrier gas, and the sizes considered
were small enough that they were completely
vaporized before reaching the stagnation plane
established by the two opposed gaseous
streams. In another detailed numerical study,
Chen et al. [4] considered the motion of drop-
lets in a nonpremixed counterflow field. They
considered both n-heptane and water droplets.
For the droplet sizes considered, Chen et al. [4]
showed that a droplet can penetrate through the
stagnation plane established by the gaseous
streams and reverse its trajectory at some point
in the opposing side of the flow. This also raised
an interesting question on details of the numer-
ical technique that could overcome the singu-
larity associated with the droplet number den-
sity equation at the point where the droplet
trajectory is reversed (see comments by Li et al.
in Ref. [4]).

One major difference between the present
work and that of Continillo and Sirignano [2] is
the description of the outer frozen flow field. In
Ref. [2], the outer flow was described by irrota-
tional or potential flow field, assuming that the
fuel and air sources are located infinitely apart,
while in the present study a more general rota-

tional flow field is used to describe the finite
separation distance between fuel and air noz-
zles. Consequently, assumptions introduced to
simplify the Lagrangian equations for droplets
to a self-similar form are slightly different.
Recent experiments by Li et al. [5] with meth-
anol fuel sprays have indeed found that the
outer flow field lies somewhere between the
potential flow and plug flow cases, similar to
that found for methane-air gaseous flames [6].
Their experiments [5] using a phase-Doppler
particle analyzer (PDPA) have also shown that,
near the axis of symmetry, the axial velocity,
droplet size, and number density are indepen-
dent of the radial coordinate, as assumed in the
present study.

While the above theoretical and experimental
studies have addressed aspects of combustion of
fuel sprays in a counterflow field, the present
study is focused on addressing the flame extinc-
tion mechanism by fine-water droplets in a
similar flow field. To our knowledge, only Se-
shadri [7] performed experiments on extinction
condition of counterflow flames with water mist,
however, the water droplets generated in that
study were so small that they were completely
vaporized in the air stream. In parallel with the
present theoretical work, an experimental effort
is currently under way at the University of
Virginia on flame extinction by monodisperse
water droplets in the size range of 20-60 wm,
but here only the details of the two-phase
numerical model developed are reported. The
numerical model developed can provide de-
tailed description of the flame structure and an
improved understanding of the parameters con-
trolling the extinction phenomena. They can
also guide the supporting experiments. The
nonpremixed flame considered here is that es-
tablished by a counterflow of methane and air,
with monodisperse-sized fine-water droplets in-
troduced with the air stream. The water droplet
loading, or the number density, considered is
small (the mass fraction of water in condensed
phase is less than 3%) such that the ratio of
droplet separation distance to droplet diameter
is greater than 20. Thus, the collisional effects
between droplets are neglected. The gas-phase
chemical reactions are described using an ele-
mentary reaction mechanism having 17 species
in 39 steps, with variable thermochemical and


170

TABLE 2

Comparison of the Stokes Number (S¢) for Different

Droplet Sizes for a = 130s

a y ius

(4) w) a) St
5 0.7523 10° 4 0.7633 x 10 0.0098
10 0.3009 x 10° 0.7633 x 10° 0.0304
15 0.6770 « 10° 0.7633 x 10? 0.0887
20 0.1204 10°? 0.7633 10-2 (1.1877
25 0.1881 x 10°? 0.7633 x 1 > 0.2464
30 0.2708 x 10°* 0.7633 x 10 ? 00.3548
40 0.4814 x 10°? 0.7633 x 10° 0.6307
50 0.7523 x 10°? 0,7033 x 10° 0.9855
0 0.1083 x 10°! 0.7633 x 10 1.4191

and deriving an expression for the droplet re-

sponse time (f,) based on Stokes drag, the

following expression can be derived for St
_ 1 18m

© a pa?
For the strain rate of 130 s~' considered in Fig.
5, the variation of Sr for different droplet sizes
in the outer oxidizer flow is shown in Table 2. It
is seen that 5-~m droplets have truly small S¢
number, consistent with the results shown in
Fig. 5. As the droplet size approaches about 15
wm, the St number approaches about 0.1, indi-
cating slower response of the droplet to the gas,
as seen in Fig. 5.

Figures 6-8 show the variation of other vari-

(29)

Ons

Con O2~O«OSSC«
z (cm)

Fig. 6, Comparison of U of gas and U,, of different droplet
sizes, with droplet source terms turned “off in gas-phase
calculations. Also shown is the gas temperature (thick line)

A.M. LENTATI AND H. K. CHELLIAH

OD

2

B80 o
04 02 0.0 02 04 06
z (cm)

Fig. 7. Comparison of droplet temperature, 7,, of different
droplet sizes, with droplet source terms turned “off” in
gas-phase catculations. Also shown is the gas temperature
(thick line).

ables associated with the droplet, cvaluated
from Eqs. 20, 21, and 23, along the axis of
symmetry. These results are again obtained
under the same conditions as in Fig. 5. Figure 6,
for example, shows that U, is always less than
U, which is expected, as the gas is accelerating
radially from the axis of symmetry and is con-
sistent with experimental measurements of
methanol droplets in a counterflow by Li et al.
[5]. Figure 7 shows that the temperature of
5-um droplets (7) follows that obtained as-

or4

i)
“04

Fig, 8. Comparison of flux-fraction function, F, of different
droplet sizes, with droplet source terms turned “off” in
gas-phase calculations. Also shown is the gas temperature
(thick line).



EFFECT OF WATER DROPLETS ON EXTINCTION 173

2000 16.0
al vapor (PV = Psat)
1960 jo~-Od= 5 microns
te -4d=50micrns | 140 eae
1900
| 1850 - 120 ’
i y
BS 1800 4 °
B 100
1750
80
1700
1650, Fs tt ed 60 — J
100 150 200 250 300 350 00 100 200 300 400 500 600

Strain Rate, a (s~')

Fig. 14. Comparison of T,,,. vs. «, for different droplet
sizes, with Y, = 0.03.

can be explained by the location of the mass
source term seen in Fig. 9. In the case of 15-um
droplets, the peak value of S,,, and also S, —
hS,,, is scen to occur at the oxygen consump-
tion or radical production layer. The negative
effect of S;, — h,S,, on flame temperature can
result in lower radical production and hence
early flame extinction [22, 23]. The exact phys-
ical, thermal, and chemical contributions re-
sponsible for flame extinction condition in the
presence of water droplets are currently being

350

150

0.0 10.0 20.0 30.0 40.0 50.0 60.0
Droplet Diameter (jm)

Fig. 15, Comparison of a,,, vs. droplet size, for different
droplet mass fractions in condensed phase at inflow,

Droplet Diameter (um)

Fig. 16. Comparison of (x,)-x, vs. droplet size, for different
droplet mass fractions in condensed phase at inflow.

investigated separately, with comparison to that
of traditional chemical agents, such as halon
1301.

The concept of flow strain rate as defined
here, ie.,a = —dv/dx, is strictly applicable to
constant density or nonreactive mixing layer
flows. However, it is a useful fluid dynamical
parameter for comparison with near extinction
experimental velocity data obtained with nonin-
trusive laser-Doppler velocimetry techniques
[6]. For theoretical studies, a more useful fluid
dynamical parameter is the diffusion time scale
defined by the inverse of the scalar dissipation
rate [23, 24]

x = 2D|¥Z/* (30)

where D is an appropriate diffusion coefficient
(e.g., D = A/(pe,) for unity Lewis number) and
Z is the mixture fraction. For small stoichio-
metric mixture fractions, Z,, an analytical ex-
pression relating the stoichiometric scalar dissi-
pation rate, y,, with a taking into account the
thermal expansion effects, has been derived by
Kim and Williams [25]. Instead of using such an
analytical expression, the scalar dissipation rate
at the stoichiometric point at extinction is
directly computed here by solving for a con-
servation equation for mixture fraction. The
resulting variation of (y,).., as a function of
the initial droplet size and water mass fraction
in condensed phase is shown in Fig. 16, indi-