Calculation of the Cellular ...

[ Pobierz całość w formacie PDF ]
Combustion, Explosion, and Shock Waves, Vol. 36, No. 6, 2000
Calculation of the Cellular Structure of Detonation
of Sprays in an H
2
{O
2
System
S. A. Zhdan
1
and E. S. Prokhorov
1
UDC 534.222.2
Translated from Fizika Goreniya i Vzryva, Vol. 36, No. 6, pp. 111{118, November{December, 2000.
Original article submitted August 7, 2000.
On the basis of the mathematical model of a two-phase two-velocity medium, deto-
nation of a cryogenic mixture (gaseous hydrogen{drops of liquid oxygen) was studied
numerically. The dynamics of formation and the special features of the structure of
the two-dimensional reaction zone of the detonation wave are discussed. The cellular
structure of detonation is modeled for the rst time for a cryogenic hydrogen{oxygen
spray.
Cellular structures of detonation in gases are
well known in experiments [1]. Cells in gas{droplet
systems were found only in sprays of small (d
0
5 m) droplets of decane [2]. The stability of det-
onation waves (DW) in homogeneous and hetero-
geneous explosives is studied theoretically in many
papers [3{9]. These waves may be described using
the model of mechanics of a one-velocity medium.
The stability of DW for multivelocity heterogeneous
reactive media, for example, sprays, has not yet
been studied analytically. This is primarily caused
by mathematical diculties arising even in the lin-
ear analysis of detonation stability in these me-
dia. Therefore, the only constructive tool for solv-
ing the problem of stability of a DW in a heteroge-
neous multivelocity medium is direct numerical simu-
lation. Thus, the one-dimensional instability of det-
onation in a cryogenic hydrogen{oxygen spray was
studied in [10], where, in particular, self-sustained
one-dimensional autooscillatory regimes of heteroge-
neous detonation were obtained. Based on the analy-
sis of these regimes, Zhdan [11] proposed a hypothesis
on the possibility of existence of cellular structures in
this spray in a wide range of droplet sizes.
Results of a numerical study of the two-
dimensional instability of detonation in a cryo-
genic hydrogen{oxygen spray are described in the
present paper. The dynamics of formation of a two-
dimensional reaction zone is discussed.
FORMULATION OF THE PROBLEM
Let a plane channel of width y
0
contain a mix-
ture of liquid oxygen droplets (with a diameter d
0
,
true density
2
, and volume concentration
20
) in
gaseous hydrogen with an initial pressure p
0
and tem-
perature T
10
. A DW propagates along the channel.
We have to nd its dynamics and structure depend-
ing on the channel width, droplet diameter, and ini-
tial composition of the mixture.
Within the framework of the model of mechan-
ics of a two-phase multivelocity medium [12], we con-
sider a two-dimensional unsteady motion of monodis-
persed oxidizer droplets in a gaseous fuel. A com-
bustion region is formed in the gas{droplet mixture
behind the shock wave (SW) after a chemical delay
of ignition. The rates of chemical reactions in this
region are much greater than the mass-transfer rate
between the phases.
We use the following assumptions in the math-
ematical model: 1) there is no energy release in the
zone of chemical induction; 2) chemical processes at
each point behind the ignition front proceed instan-
taneously, and the composition of the gas phase is
in chemical equilibrium; 3) the energy-release rate
in the gas{droplet mixture is limited by the phase-
transition velocity, and the thermal eect of chemical
reactions per unit mass of evaporated droplets and
the molar weight of the gas are variable quantities.
1
Lavrent'ev Institute of Hydrodynamics,
Siberian Division, Russian Academy of Sciences,
Novosibirsk 630090.
0010-5082/00/3606-0777 $25.00
c
2000
Plenum Publishing Corporation
777
778
Zhdan and Prokhorov
The equations of a two-dimensional unsteady
two-phase ow are written in the following form:
(
i
)
t
+ (
i
u
i
)
x
+ (
i
v
i
)
y
= (1)
i+1
j;
U
2;ch
= E
2
=
O
2
, and c
2
, l
2
, and E
2
are the heat ca-
pacity, the latent heat of vaporization, and the energy
of dissociation of molecules of the condensed phase,
respectively.
If the internal energy of the gas phase U
1
is reck-
oned from the maximally dissociated composition at
zero temperature, the following relations are valid in
the induction zone:
n
t
+ (nu
2
)
x
+ (nv
2
)
y
= 0;
(
i
u
i
)
t
+ (
i
u
i
)
x
+ (
i
u
i
v
i
)
y
+
i
p
x
= (1)
i
(f
x
ju
2
);
p
(
1
1)
1
U
1;th
=
;
(
i
v
i
)
t
+ (
i
u
i
v
i
)
x
+ (
i
v
i
)
y
+
i
p
y
(1)
(4)
= (1)
i
(f
y
jv
2
);
U
1;ch
=
E
2
z
O
2
E
1
(1 z)
:
H
2
(
i
E
i
)
t
+ [
i
u
i
(E
i
+ p=
i
)]
x
+ [
i
v
i
(E
i
+ p=
i
)]
y
= (1)
i
(f
x
u
2
+ f
y
v
2
jE
2
);
Here
1
is the ratio of specic heats in the induction
zone and E
1
and
H
2
are the dissociation energy and
the molar weight of hydrogen, respectively.
According to [14, 15], the thermodynamic and
chemical components of internal energy of the gas in
the region of chemical transformations have the form
Y
t
+ u
1
Y
x
+ v
1
Y
y
= 1=t
ind
;
z
t
+ u
1
z
x
+ v
1
z
y
= (1 z)j=
1
;
E
i
= U
i
+ (u
i
+ v
i
)=2;
i
=
i
i
;
1
+
2
= 1:
Here
i
and
i
and
i
are the volume concentration
and the mean and true densities of the phases, respec-
tively, u
i
and v
i
are the components of the velocity
vector u
i
, U
i
is the total internal energy of the ith
phase (i = 1 and 2), which is reckoned from the max-
imally dissociated composition at zero temperature,
p is the gas pressure, n is the number of droplets per
unit volume, f and j are the intensities of the force
and mass interactions between the phases, z is the
mass fraction of the evaporated substance of the sec-
ond phase in the gas regardless of the component it
enters, and Y is the fraction of the induction period
(Y = 1 at the SW front and Y = 0 at the ignition
front).
Based on the experimental data of [13], the
chemical delay of ignition for a hydrogen{oxygen
mixture behind the SW front, in the induction zone
(Y > 0), has the form
t
ind
=
K
a
O
2
1
z
h
min
1
2
U
1;th
=
+
min
=T
1
exp(=T
1
) 1
i
RT
1
;
+
+ 1
U
1;ch
= E
d
1
;
(5)
min
exp
"
a
(2)
where =
max
(
min
)=(
max
min
),
min
and
max
are the molar weights of the gas in the max-
imally dissociated (in this case, atomic) and max-
imally recombined states, respectively,
max
is the
molar fraction of triatomic molecules in the maxi-
mally recombined state, is the eective tempera-
ture of excitation of vibrational degrees of freedom
of molecules, and E
d
is the mean energy of disso-
ciation of the reaction products. All the parame-
ters are uniquely determined by the atomic composi-
tion of the mixture (by the mass fraction z of evap-
orated oxygen). For the hydrogen{oxygen mixture,
we have
max
= (1z)=
H
2
and
max
= 2z
max
=
O
2
for z < 8=9,
max
= z=
O
2
+ 0:5(1 z)=
H
2
and
max
= (1 z)
max
=
H
2
for z > 8=9,
min
=
2(z=
O
2
+ (1 z)=
H
2
), E
d
E
10
E
20
, and
= 3000 500
max
.
Since the rates of chemical reactions behind
the ignition front are much greater than the mass-
transfer rate between the phases, the composition of
the gas phase may be considered to be in chemical
equilibrium. Then behind the ignition front (in the
reaction zone), we may use the approximate equa-
tion for the shift in chemical equilibrium, which is
in agreement with the second principle of thermody-
namics [15]:
where "
a
is the activation energy, R is the universal
gas constant, K
a
is the preexponent, and
O
2
is the
molar weight of oxygen.
We supplement system (1) by the equations of
state of the phases
p =
1
RT
1
=;
2
= const;
(3)
U
i
= U
i;th
+ U
i;ch
;
where T
i
and U
i;th
and U
i;ch
are the temperature and
thermodynamic and chemical components of inter-
nal energy of the ith phase, is the current mo-
lar weight of the gas phase, U
2;th
= c
2
T
20
l
2
,
RT
1
;
 Calculation of the Cellular Structure of Detonation of Sprays in an H
2
{O
2
System
779
1
(1 =
max
)
2
=
min
1
exp
E
d
RT
1
i
=
i0
; u
i
= v
i
= 0;
U
i
= U
i0
(i = 1; 2); p = p
0
:
The boundary conditions are the no-slip condi-
tion of the gas phase (v
1
= 0) at the lower (y = 0) and
upper (y = y
0
) boundaries, transposition of param-
eters from the solution domain at the left boundary
moving with the mean SW velocity (this is correct
since the gas ow is supersonic relative to the moving
boundary), and relations on a strong discontinuity at
the right moving boundary [19].
Problem (1){(7) was solved numerically. To in-
tegrate the system of equations that describe the
gas ow, we used the second-order Godunov{Kolgan
scheme in moving grids [19, 20] with capturing of the
leading shock front and ignition front, where the ther-
mal decay of the discontinuity was calculated using
the procedure of [21]. The equations that describe
the motion of droplets were solved by the method of
coarse particles [22]. The accuracy of solving the dif-
ferential equations (1) was controlled by disbalance
in integral laws of conservation of mass and energy,
which was less than 1% in all calculations.
T
1
T
10
=2
h
T
1
i
= K
1
1 exp
: (6)
Here = 1 +
max
min
=(
max
min
) and K
1
is the
equilibrium constant.
The closing relations for the intensities of the
mass (J = j=
2
) and force (f) interactions between
the phases have the following form [12, 16, 17]:
8
<
6(1:5)
1=2
2
(
1
=
2
)
1=3
(
1
=
2
)
2=3
Re
1=2
=d
2
for We > 10;
J =
:
3K
vap
(1 + 0:27Pr
1=3
Re
1=2
)=d
2
for We < 10;
f = nd
2
1
C
D
ju
1
u
2
j(u
1
u
2
)=8; (7)
27Re
0:84
for Re < 80;
0:27Re
0:21
for 80 6 Re < 10
4
;
2 for Re > 10
4
:
Here d is the current diameter of droplets, Re =
1
dw=
1
is the Reynolds number, We = d
1
w
2
=
is the Weber number, Pr is the Prandtl number,
K
vap
is the vaporization constant, w = ju
1
u
2
j,
i
=
i
=
i
and
i
are the kinematic and dynamic
viscosities of the ith phase, respectively, the dynamic
viscosity of the gas has a power-law dependence on
temperature
1
=
0
(T
1
=300)
0:7
, and
0
is calculated
using Herning's method [1
8]:
0
= (
1a
)
H
2
+a
O
2
8
<
C
D
=
:
CALCULATION RESULTS
and a = z
p
O
2
=((1 z)
p
H
2
+ z
p
O
2
).
Note that the expressions for the mass-transfer
intensity J correspond to the model of droplet-mass
entrainment by the mechanism of boundary-layer
stripping, which transforms to the mechanism of
droplet vaporization for We 6 10, and are valid [16]
for small droplets (d
0
6 100 m). Exactly this range
of oxygen droplet diameters will be considered in the
present paper.
System (1){(7) is closed and completely de-
termines the two-dimensional unsteady two-velocity
motion of a reactive gas{droplet mixture with vari-
able heat release in the DW reaction zone.
The numerical study was performed for a
cryogenic hydrogen{oxygen mixture, which was a
monodispersed spray of oxygen droplets in gaseous
hydrogen (2H
2
+ O
2
) for the initial temperature of
the phases T
10
= T
20
= 80 K with the following
thermodynamic parameters [23]: (a) for gas phase,
H
2
= 2 kg/kmole,
H
2
= 0:9 10
5
kg/(msec),
O
2
= 2 10
5
kg/(msec), E
1
E
2
E
d
=
110 kcal/mole, K
1
= 1:81 10
2
kmole/m
3
,
1
= 1:4,
and Pr = 0:75; (b) for droplets of O
2
,
O
2
=
32 kg/kmole,
2
= 1135 kg/m
3
, c
2
= 1680 J/(kgK),
l
2
= 0:213 MJ/kg,
2
= 2:5 10
4
kg/(msec),
= 0:0131 N/m, and K
vap
= 2 10
6
m
2
/sec. The
constants in Eq. (2) for the chemical delay of igni-
tion were K
a
= 5:38 10
11
(molesec)/liter and
"
a
= 17:15 kcal/mole [13].
For these thermodynamic parameters of the
phases, the solution of system (1){(7) with initial
and boundary conditions depends on three governing
parameters: initial pressure of gaseous hydrogen p
0
and two scale factors | oxygen droplet diameter d
0
and channel width y
0
. For a xed value of p
0
, the
initial mean densities of gaseous hydrogen
10
and
oxygen droplets
20
=
20
2
are determined unam-
biguously. In particular, for p
0
= 1 atm, we have
10
= 0:305 kg/m
3
and
20
= 2:44 kg/m
3
; the ve-
locity of the ideal Chapman{Jouguet detonation is
D
CJ
= 2:97 km/sec.
INITIAL AND BOUNDARY CONDITIONS
The initial conditions in the solution domain
(0 < y < y
0
; x > 0) are as follows: a one-dimensional
steady solution in the reaction zone of ideal heteroge-
neous Chapman{Jouguet detonation for a monodis-
persed spray of oxygen droplets in gaseous hydrogen
for 0 < x < x
and constant initial values for x > x
:
 780
Zhdan and Prokhorov
Fig. 1. Velocity of the front of the one-
dimensional wave versus time for an autooscil-
latory detonation mode (p
0
= 1 atm and d
0
=
100 m).
One-Dimensional Instability of the DW.
It is known [10] that detonation of a cryogenic
hydrogen{oxygen spray in the one-dimensional for-
mulation (in contrast to heterogeneous detonation
of hydrocarbon{oxygen sprays [11]) is unstable and
propagates in a pulsating (autooscillatory) mode.
The unstable element in the wave structure, which
leads to periodic oscillations of parameters in the
reaction zone, is the ignition front. A typical de-
pendence of the wave-front velocity D on the time t
for p
0
= 1 atm and d
0
= 100 m is shown in
Fig. 1. It also follows from one-dimensional calcu-
lations that the time (t) and space (L
0
) periods of
longitudinal oscillations of the wave velocity, its am-
plitude, and detonation velocity averaged over the
period (D
0
= L
0
=t) depend on the initial diameter
of oxygen droplets. The qualitative information on
these quantities for a number of values of d
0
is given
in Table 1.
It is seen that the values of D
max
, D
min
, and
D
0
depend weakly on d
0
, whereas the period of self-
induced oscillations is an almost linear function of
the oxygen droplet diameter: L
0
d
1:1
Fig. 2. Longitudinal SW velocity versus time on the
upper (solid curves) and lower (dashed curves) walls of
the channel for y
0
= 26:6 (a), 22.15 (b), and 17.7 mm
(c).
0
. It should
be noted that, because of the nonmonotonic ther-
mal eect of chemical transformations in the reac-
tion zone [10], the mean velocity of one-dimensional
oscillating detonation D
0
is almost 8% greater than
the velocity of the steady Chapman{Jouguet detona-
tion (D
CJ
).
Two-Dimensional Instability of the DW.
The numerical algorithm and the code developed for
solving the two-dimensional unsteady problem (1){
(7) allow one to study the dynamics of propagation
of a heterogeneous DW in a plane channel with var-
ied channel width and oxygen droplet diameter. First
we consider the inuence of the channel width on the
dynamics of the detonation process for xed values
p
0
= 1 atm and d
0
= 100 m. Figure 2 shows the
calculated longitudinal velocity of the shock front D
versus the time t on the upper and lower walls for
several values of the channel width: y
0
= 26:6, 22.15,
and 17.7 mm. It is seen that the dependences D(t)
dier not only from the one-dimensional solution (see
Fig. 1) but also from one another. As the DW (ini-
tially one-dimensional) propagates, a crossow insta-
bility of the ignition front develops in the reaction
zone, which leads to transverse oscillations of gas-
dynamic parameters of the ow. The ow pattern
becomes two-dimensional, and the shock-front veloc-
ity experiences irregular oscillations (see Fig. 2c). By
varying y
0
, we show that oscillations of the SW ve-
locity can appear for t > 0:5 msec, which are regu-
lar with respect to the x axis and have dierent (see
TABLE 1
d
0
,
m
t,
sec
D
max
,
km/sec
D
min
,
km/sec
D
0
,
km/sec
L
0
,
mm
TABLE 2
d
0
,
m
D
max
,
km/sec
D
min
,
km/sec
hDi,
km/sec
a a=b
25
2.581
4.71
2.43
3.29
8.523
50
5.809
4.66
2.40
3.22
18.68
50
7.34
1.94
3.09
21.1 0.585
100 12.08
4.47
2.42
3.19
38.50
100
5.58
2.05
3.04
44.3 0.50
Calculation of the Cellular Structure of Detonation of Sprays in an H
2
{O
2
System
781
Fig. 3. Two-dimensional cellular structure of the DW in a heterogeneous mixture 2H
2
+ O
2
during
one period (d
0
= 100 m and a = 44:3 mm): the solid curve shows the calculated trajectory of the
triple point; (a) isobars of the gas phase; (b) isochores of the gas phase; (c) isochores of oxygen
droplets; .
Fig. 2a) and identical (see Fig. 2b) amplitudes on the
opposite walls of the channel. There exists a mini-
mum value of the channel width y
0
= y
= 22:15 mm
for which the SW velocity along the x axis at a dis-
tance x=y
0
> 70 in each longitudinal section changes
periodically (see Fig. 2b) with a time period t =
29:15 sec and a spatial period b = 88:6 mm. A
regular transverse perturbation is formed in the DW
structure, which is alternatively reected from the
channel walls so that the SW velocity and other gas-
dynamic parameters on the opposite walls change in
the same manner but in the opposite phase (shifted
by half a period). The calculated values of D
max
and
D
min
, and the wave velocity averaged over the pe-
riod hDi = b=t are listed in Table 2. The quantity
y
, which will be called the eigenvalue of problem
(1){(7), orders the two-dimensional instability of a
heterogeneous DW in the form of a regular periodic
solution (going a few steps forward, we note that y
corresponds to half of the transverse size of the deto-
nation cell a=2). The periodic solutions found in nu-
merical simulation are of special interest, and we will
dwell on them in more detail. We note that the longi-
tudinal velocity of the shock wave near the wall D(t)
is rather sensitive to variation of the channel width;
therefore, it is a convenient parameter for control-
ling the accuracy of determining y
. It follows from
the calculations that, for a deviation of the channel
width y
0
from y
by 5%, the velocity proles on the
opposite walls of the channel do not coincide, whereas
the shape of the leading front is still unchanged. The
fact that the SW velocity in a heterogeneous DW
has two maxima during one period (see Fig. 2b) is
worth noting. This fact will be explained below by
analyzing the structure of the reaction zone of the
heterogeneous DW.
The calculated periodic two-dimensional DW
structure (y
0
= y
) in a heterogeneous mixture
2H
2
+ O
2
during one period is shown in Fig. 3. The
isobars (P = p=p
0
) and isochores (R
1
=
1
=
20
) of
the gas phase are shown in Fig. 3a and b, and the iso-
chores of oxygen droplets (R
2
=
2
=
20
) are plotted
in Fig. 3c. To control the reliability of the periodic
solution found, we performed additional calculations,
where the channel width was doubled (y
0
= 2y
= a)
for a xed moment of time, and a mirror reection
of the distribution of gas-dynamic parameters from
the solution domain y
< y < 2y
was transferred
to the region y
< y < 2y
. The initial data ob-
tained in this way were used to calculate an unsteady
problem. It is shown that a self-sustained heteroge-
neous DW propagates over the channel. The reaction
zone of this DW has a two-dimensional periodic so-
lution with two opposing transverse waves, which is
symmetric about the channel centerline 0 < y < y
.
The internal structure of this DW coincides identi-
[ Pobierz całość w formacie PDF ]