Combustion, Explosion, and Shock Waves, Vol. 37, No. 3, pp. 315 320, 2001
Simulation of the Erosion Burning of a Granular Propellant
M. M. Gorokhov1 and I. G. Rusyak1 UDC 536.46
Translated from Fizika Goreniya i Vzryva, Vol. 37, No. 3, pp. 76 82, May June, 2001.
Original article submitted February 17, 2000; revision submitted June 2, 2000.
For combustion of axisymmetric propellant grains under blowing conditions, a math-
ematical model is proposed and numerical simulation is performed. The effects of
incoming-flow parameters (velocity, pressure, and temperature) and surface dimen-
sions and geometry on grain-burning rate are studied. Physical patterns of flow
around burning propellant grains are presented.
Granular propellant particles (or powder particles) An analysis of the literature on the problem of ero-
are used in many technological processes and equip- sion burning of propellants showed that all studies were
ment. Combustion of such particles, as a rule, occurs performed for objects of simple geometry such as a plate
under conditions of intense flow of an external jet or or a channel. For granular propellants and grained pow-
combustion products. Examples are the combustion of ders there are no literature data. In the present paper,
coal particles in pressurized furnaces, ignition of a solid we first undertake an attempt to calculate the erosion
rocket motor (SRM) charge by a multiphase jet con- burning rate for granular propellants.
taining igniter grains, and combustion of granular pro- In many cases of practical importance, the relax-
pellants in artillery systems during shot. Obviously, for ation time in gas and condensed combustion zones is
more accurate calculations of parameters of such pro- much smaller than the characteristic time of change
cesses and phenomena, it is necessary to increase the of internal-ballistic parameters in SRM or artillery sys-
accuracy in predicting the grain-burning rate with blow- tems. In such cases, propellant combustion can be stud-
ing. ied in a quasistationary formulation. In addition, in
The effect of increase in propellant-burning rate reality, the propellant burning rate and the blowing ve-
due to increase in incoming-flow velocity was first re- locity differ by several orders of magnitude. Because
vealed by Leipunskii [1] in experimental studies of of this, the dynamics of change in grain size can be ig-
double-base propellant combustion under blowing con- nored, although for grain sizes comparable to the heated
ditions. Physical explanations of the nature of this condensed-phase (c-phase) layer, this assumption is not
phenomenon are given by Zel dovich [2] and V. I. Vi- valid. In addition, the plane formulation of the com-
lyunov [3]. bustion problem should be corrected in determining the
As applied to SRM problems, a coupled problem c-phase burning rate. Below, we study the combustion
of chemical kinetics and hydrodynamics (thermal prob- of grains with characteristic size much larger than the
lem of erosion burning of propellants) in a plane hydro- thickness of the heated c-phase layer.
dynamical formulation was considered in [3]. However, The proposed formulation of the problem includes
the analytical approach used in [3] required that a num- simultaneous solution of the equations of hydrodynam-
ber of serious assumptions be introduced. In particular, ics and chemical kinetics for propellant gasification
the assumption of an asymptotic flow regime was used, products near axisymmetric bodies of various configura-
which is justified only for large blowing velocities. The tions. The c-phase burning rate is determined from its
effect of negative erosion was revealed experimentally surface temperature using the well-known Merzhanov
by Vilyunov and Dvoryashin [4]. Bulgakov and Lipanov Dubovitskii analytical solution [6]. This approach al-
[5] proposed a detailed, physically founded mathemati- lows one to study the mutual effect of the external flow
cal model for propellant combustion with blowing and and propellant-combustion processes over a wide range
explained the mechanism of negative erosion. of velocity, pressure, and incoming-flow velocity.
1
Izhevsk State Technical University, Izhevsk 426069.
0010-5082/01/3703-0315 $25.00 © 2001 Plenum Publishing Corporation 315
316 Gorokhov and Rusyak
We consider stationary, substantially subsonic (at
Mach numbers M < 0.3) flow around a single fixed par-
ticle of axisymmetric shape in the range of Reynolds
numbers constructed from the midsection diameter
Re = 10 2 · 105. In this range up to the separation
point, the body boundary layer is laminar [7] and tur-
bulence can be ignored in the calculations. Gas-phase
combustion can be represented as the result of two suc-
cessive overall reactions [8].
Fig. 1. Object of investigation, coordinate system, and
The formulation of the problem of homogeneous
domain of integration.
solid propellant combustion in flow of a viscous incom-
where p is the pressure and R is the gas constant of the
pressible heat-conducting gas thinks into account the
two-stage nature of the chemical transformation of gasi- combustion products.
Our investigations showed that the traditional
fication products includes:
spherical coordinate system is ineffective in this case be-
 the equation of motion
cause it cannot be used to study the burning of particles
" · Ávv = " · P , (1)
that deviate considerably from the spherical shape (see,
 the continuity equation
e.g., Fig. 2). In addition, in the case of blowing, the
flow is deformed over the stream, and it is impossible
" · Áv = 0, (2)
to perform equally accurate calculations over the entire
 the heat-transfer equation
particle combustion zone in a region with spherical sym-
µ
metry. In numerical analysis of flow past burning bodies
" · ÁvT = " · "T
Pr
of complex shape, it is reasonable to use finite-difference
grids attached to the body surface and domains of inte-
+ Q1f1(C1, T ) + Q2f2(C2, T ), (3)
gration adapted to the flow field. This approach applies
 the equations of chemical kinetics
to calculations for axisymmetric bodies of various geom-
µ
etry by a unified scheme and simplifies the formulation
" · ÁvC1 = " · "C1 - f1(C1, T ), (4)
Sc
and numerical implementation of boundary conditions.
µ In addition, it simplifies the solution of the problem re-
" · ÁvC2 = "· "C2 +f1(C1, T ) - f2(C2, T ). (5)
lated to the occurrence of scheme viscosity and allows
Sc
grid lines to be made finer near a solid surface [9, 10].
Here Á is the gas density, v is the flow-velocity vector,
The last circumstance is important in this case because
P is the stress tensor, T is the temperature, µ is the
the spatial gas-dynamic scale and combustion-zone di-
molecular viscosity, Pr is the Prandtl number, Sc is the
mensions depend on flow conditions and can differ by
Schmidt number, and Qi is the heat effect of the first
several orders of magnitude.
(i = 1) and second (i = 2) stages of the chemical reac-
The system of equations and boundary conditions
tion, Ci is the concentration of reaction products at the
were written in the orthogonal axisymmetric coordinate
corresponding stage (i = 1, 2), and R0 is the univer-
system shown in Fig. 1. Numerical integration of sys-
sal gas constant, fi(Ci, T ) is the chemical-reaction rate
tem (1) (6) was performed in the region bounded by
for the corresponding stage (i = 1, 2), defined by the
the symmetry line AD, the body contour (burning sur-
formula
face) FE and the entry (AB), exit (CD), and outer (BC)
Ei
½i
i
fi(Ci, T ) = Á½ Ci Zi exp - ,
boundaries located in the external flow.
R0T
The orthogonal curvilinear finite-difference grid is
where Ei, Zi, and ½i are the activation energy, preexpo-
constructed by the complex method of boundary ele-
nent, and order of the overall reaction of the first and
ments [9 11]. Examples of construction of such grids
second stages.
are given in Fig. 2.
In the calculations, the Prandtl and Schmidt num-
On the burning surface, the boundary conditions
bers are assumed to be constant. The effect of tempera-
are
ture on the molecular viscosity of combustion products
u = 0, Áv = Ácvc,
is taken into account by the known dependence pro-
posed by Sutherland [7].
c("·Tc) = ("·T ) - (cp - cc)ÁcvcTs, (7)
System (1) (5) is closed by the equation of state
p
µ µ
= RT, (6)
("·C1) = Ácvc(C1 - 1), ("·C2) = ÁcvcC2,
Á
Sc Sc
Simulation of the Erosion Burning of a Granular Propellant 317
u· = v = T· = (C1)· = (C2)· = 0.
On the entry and upper boundaries, the flow parameters
take the following values:
p = p", u = u", v = 0,
T = T", C1 = C2 = 0.
Here p", u", and T" are, respectively, the pressure,
velocity, and temperature of the incoming flow. On
the exit boundary, we specify unperturbed-flow pres-
sure and  soft boundary conditions for the remaining
variables:
p = p",
u¾ = v¾ = T¾ = (C1)¾ = (C2)¾ = 0.
The flow field is calculated by the Simple method [12].
The asymptotic convergence of the method was es-
tablished by making the difference grid finer and vary-
ing the dimensions of the domain of numerical inte-
gration. Numerical calculations were carried out on a
finite-difference grid with the following parameters: 160
grid cells were specified in the ¾ direction; 40 of them
Fig. 2. Examples of construction of finite-difference
were placed on the body surface and 60 cells were lo-
grids for calculation of flow parameters about axisym-
metric bodies: (a) ellipsoid shape; (b) dumbbell shape.
cated in the · direction. Further increase in the number
of nodes in the ¾ and · directions is inexpedient because
it does not lead to a significant change in the burning
where u and v are the projections of the velocity vector
rate averaged over the surface (a two-fold increase in the
v onto the ¾ and · axes, respectively, vc, c, cc, Ác, and
number of nodes changes the burning rate by 0.1%). In
Tc are the stationary burning rate, thermal conductiv-
the ¾ direction, the space ahead of the body was var-
ity, specific heat, density, and temperature of the solid
ied in the range (3 15)Rm, and the space behind the
propellant, respectively, Ts is the burning surface tem-
body was varied in the range (15 30)Rm, where Rm is
perature of the solid propellant, and  and cp are the
the midsection radius of the body. All changes depend
thermal conductivity and specific heat of combustion
on the Reynolds number, whose increase leads to a de-
products in constant volume. The stationary burning
crease in space before and above the body and to an
rate of the solid propellant is defined by the formula [6]
increase in space behind the body. The grid lines were
2
made finer toward the body surface so that a sufficient
2(c/ccÁc)ZcR0Ts Ec
vc = exp , (8)
(for the required accuracy) number of cells were in the
Ec(2Ts - Ts,1 - Tin) 2R0Ts
regions of considerable velocity gradients in the bound-
where Tin is the initial propellant temperature, Ts,1 =
ary layer and considerable temperature gradients in the
Tin + Qc/cc is the propellant surface temperature for
combustion zone.
flashless combustion, Qc, Ec, and Zc are the heat ef-
The problem of flow around a sphere was consid-
fect, activation energy, and preexponent of the overall
ered as the test problem. We compared the position of
condensed-phase reaction.
the separation point of the boundary layer, the depen-
The temperature gradient on the condensed-phase
dences of drag coefficients (see [9]) and heat-exchange
surface in the direction toward the interior of the solid
coefficients of the sphere [13] and the length of the cir-
propellant is given by
culation zone behind the sphere [14] on the Reynolds
c("·Tc) = ccÁcvc(Ts - Ts,1). (9) number. In all cases, calculations agree with experi-
ment.
If the surface does not burn, the boundary condi-
The combustion of grains (Rm = 5 · 10-3 m) was
tions on the body surface become
studied numerically in the range of u" = 2 200 m/sec.
u = 0, v = 0, T = Ts, C1 = C2 = 0.
Most of the calculations were carried out at p" =
On the symmetry line, we have 10 MPa. The incoming-flow temperature was equal to
318 Gorokhov and Rusyak
TABLE 1
Heat Effects of Overall Reactions
of Double-Based Propellant N at Various Pressures
p, MPa
Parameters,
kJ/kg
1 2 5 10
Qc 314 377 469 523
Q1 1360 1403 1507 1507
Q2 0 925 1327 1340
The normal grain burning rate without blowing vk,0
(u" = 0) cannot be determined by the developed cal-
culation procedure. Therefore, as the normal burning
rate, we used the average-integral value of the grain
surface burning rate at u" = 2 m/sec because at this
value, the effect of the external flow is practically ruled
Fig. 3. Distribution of the erosion coefficient on the
out (the particle surface temperature becomes uniform
surface of a spherical grain at various blowing velocities.
and takes a value of Ts = 717 K).
The variation in the local erosion coefficient µ =
the temperature of combustion products of the tested
vc/vc,0 along the generator surface of a spherical grain
propellant at constant pressure unless otherwise speci-
is presented in Fig. 3. Calculations show that the ero-
fied. The calculations were performed for double-based
sion effect is most pronounced on the frontal and rear
propellant N, for which thermal characteristics and for-
surfaces of the burning grain. This is due to the fact
mally kinetic constants of the condensed and gas phases
that in the frontal part there is dynamic pressing of
of the combustion zone are explored in [8, 15, 16].
the gas-phase combustion zone to the grain surface. In
In the present work, we used the following thermal
the rear part of the grain there is intense vortex for-
and formally kinetic characteristics:
mation and, as a consequence, the rate of heat transfer
for the c-phase,
to the burning surface increases. On the upper part of
Ác = 1600 kg/m2, cc = 1465 J/(kg · K),
the sphere, The effect of negative erosion is observed at
u" < 200 m/sec.
c = 0.235 J/(m · sec · K),
Figure 4 gives a curve of the integral erosion co-
efficient µ averaged over the surface of the spheri-
Ec = 82.9 kJ/(g · mole), Zc = 6 · 109 1/sec;
cal grain versus incoming-flow velocity. It is evident
for the gas phase,
that the erosion effect is strongly affected by the blow-
cp = 1466.5 J/(kg · K), R = 330 J/(kg · K),
ing velocity. In this range of incoming-flow veloci-
ties, grain combustion is characterized by negative (at
Pr = 1, Sc = 1,
2 < u" 75 m/sec) and positive (at u" 75 m/sec)
E1 = 29.3 kJ/g-mole,
Z1 = 0.29 · 106 m0.6/(kg0.2· sec), ½1 = 1.2,
E2 = 167.4 kJ/g-mole,
Z2 = 0.40 · 108 m3/(kg · sec), ½2 = 2.0,
1.5
T TC + T0
µ = µ0Pr .
T0 TC + T
Here µ0 = 0.175 · 10-4 nsec/m2, T0 = 273 K, and TC =
628 K is the Sutherland constant.
Table 1 gives the heat effects of the overall reactions
Fig. 4. Average-integral (over the spherical grain sur-
for the condensed and gas phases Qc, Q1, and Q2 of face) value of the erosion coefficient versus incoming-
flow velocity.
double-based propellant N at various pressures.
Simulation of the Erosion Burning of a Granular Propellant 319
Fig. 6. Streamlines for flow around burning nonspher-
ical grains: (a) ellipsoid; (b) dumbbell.
Fig. 5. Distribution of the local erosion coefficient over
grain surfaces of ellipsoid (1) and a dumbbell (2) shapes.
ultimately, leading to positive values of the erosive ef-
fect. Here it should be noted that the additional normal
component of the gas velocity in the chemical-reaction
erosion. The minimum (0.79) and maximum (1.08) val- zone can be smaller, commensurable or larger than the
ues of the average-integral erosion coefficient correspond velocity of gas outflow from the grain surface due only
to values of u" = 10 and 200 m/sec. Additional calcu- to combustion (without blowing). Its values depend
lations showed that the nonmonotonic behavior of the weakly on the burning rate and are about 1 m/sec.
curve of µ (u") is explained by competition of heat The proposed mathematical model, which takes
fluxes from the chemical-reaction zone and the exter- into account the interaction of the external flow and
nal flow. In the normal combustion regime, there is processes in the combustion zone near the propellant
no convective heat transfer from the external flow and surface, allows one to estimate the effect of flow pressure
the combustion is determined by the heat flux from the and temperature on the change in grain-burning rate.
chemical-reaction zone. With increase in the flow ve- Calculations showed that a rise in pressure diminishes
locity, the boundary-layer thickness decreases but an the erosion effect. Thus, variation in the pressure p" in
additional dynamic normal component of the gas out- the range of 1 10 MPa, (u" = 100 m/sec) corresponds
flow velocity appears near the burning surface, which to variation in the coefficient µ from 1.39 to 1.03. At
results in additional displacement of chemical-reaction the same time, with increase in incoming-flow temper-
products from the burning surface. As a result, at this ature, the erosion effect enhances. As T" varies from
stage, the decrease in heat flux to the burning surface 1100 to 3200 K (u" = 100 m/sec), the coefficient µ
from the chemical-reaction zone is more considerable varies in the range 0.95 1.16. The results obtained are
than the increase in heat flux from the external flow in qualitative agreement with the corresponding data
region. As a result, the burning rate decreases and for flow around a planar propellant surface.
the negative erosion effect takes place. With further The developed numerical procedure can be used to
increase in the external-flow velocity, the dynamic dis- study the effect of surface geometry and grain geometry
placement of gaseous products of propellant decomposi- and size on burning rate in the case of blowing. Thus,
tion from the grain surface increases monotonically, and for example, to establish the dependence of the erosion
the heat transfer from the chemical-reaction zone to the effect on surface geometry, we performed calculations
burning surface tends asymptotically to zero. Simulta- of model grains, whose results are presented in Fig. 5.
neously, the convective component of the heat flux in- At u" = 100 m/sec, the average-integral erosion co-
creases (practically linearly), making a major contribu- efficient is equal to 0.97 and 1.02 for grains 1 and 2,
tion to the development of the combustion process, and, respectively. The flow patterns around these grains are
320 Gorokhov and Rusyak
given in Fig. 6, from which it is evident that in the re- 6. A. G. Merzhanov and F. I. Dubovitskii,  Theory of
gion of  cutout of grain 2, a vortex zone forms, leading steady-state propellant combustion, Dokl. Akad. Nauk
to a considerable increase in burning rate in this region. SSSR, 129, No. 1, 153 157 (1959).
7. L. G. Loitsyanskii, Mechanics of Liquids and Gases, Perg-
Investigation of the dependence of the erosion effect
amon Press, Oxford New York (1966).
on grain geometry and dimensions showed that these
8. A. A. Zenin,  Processes in combustion zones of double-
factors influence significantly grain-burning character-
base propellants, in: Physical Processes in Combustion
istics. As the longitudinal dimension of the ellipse (at
and Explosion [in Russian], Atomizdat, Moscow (1980),
constant midsection radius Rm = 5 · 10-3 m) increases
pp. 68 104.
in the range of (0.3 5)Rm, the average-integral value
9. M. M. Gorokhhov, I. G. Rusyak, and V. A. Tenenev,
of µ decreases from 1.08 to 0.88. With decrease in the
 Numerical studies of flow past axisymmetric bodies with
midsection radius in the range (1 0.2)Rm (at the con-
blowing from the surface, Izv. Vyssh. Uchebn. Zaved.,
stant longitudinal dimension Rm = 5 · 10-3 m), the
Mekh. Zhidk. Gaza, No. 4, 162 166 (1996).
average-integral value of µ increases from 1.03 to 1.10.
10. V. A. Tenenev, I. G. Rusyak, and M. M. Gorokhhov,
A two-fold decrease in the radius of a spherical particle
 Numerical studies of aluminum particle combustion in
led to a change in the value of µ from 1.03 to 1.39.
a two-phase flow, Mat. Model., 9, No. 5, 87 96 (1997).
11. T. V. Hromadka and C. Lai, The Complex Vari-
able Boundary Element Method in Engineering Analysis,
CONCLUSIONS
Springer-Verlag, New York (1987).
12. S. V. Patankar, Numerical Heat Transfer and Fluid
" The numerical studies of combustion of spherical
Flow, Hemisphere McGraw-Hill, Washington New York
grains of propellant N at external-flow velocities of 2
(1980).
200 m/sec showed that at this incoming-flow velocities,
13. E. R. Ekkert and R. M. Drake, Theory of Heat and Mass
grain combustion is characterized by both negative (at
Transfer [Russian translation], Gosenérgoizdat, Moscow
u" < 75 m/sec), and positive (at u" 75 m/sec)
Leningrad (1961).
erosion effects.
14. O. M. Belotserkovskii,  Numerical Simulation in Con-
" The effect of the pressure and temperature of
tinuum Mechanics [in Russian], Nauka, Moscow (1984).
the blowing flow on the change in grain-burning rate
15. V. E. Zarko, V. F. Mikheev, S. V. Orlov, et al.,  Ignition
is determined. It is shown that with rise in pressure,
of propellant by a hot gas, in: Combustion and Explo-
the erosion burning effect is diminished and with rise in
sion [in Russian], Nauka, Moscow (1972), pp. 34 37.
temperature, it is enhanced.
16. A. A. Zenin,  Formal kinetic characteristics of reac-
" The effect of surface geometry and grain geome-
tions involved in propellant combustion, Fiz. Goreniya
try and sizes on the rate of propellant combustion with
Vzryva, 2, No. 2, 28 32 (1966).
blowing is studied.
The authors thanks V. A. Tenenev for assistance
in choosing the finite-difference grid and computational
method.
REFERENCES
1. O. I. Leipunskii, Doct. Dissertation, Inst. of Chem. Phys.,
Acad. of Sci. of the USSR, Moscow (1945).
2. Ya. B. Zel dovich,  On the theory of propellant combus-
tion in a gas flow, Fiz. Goreniya Vzryva, 7, No. 4, 163
176 (1971).
3. V. N. Vilyunov,  On the theory of erosion burning of
propellants, Dokl. Akad. Nauk SSSR, 136, No. 2, 381
384 (1961).
4. V. I. Vilyunov and A. A. Dvoryashin,  Combustion
of double-based propellant N in a gas stream, Fiz.
Goreniya Vzryva, 7, No. 1, 45 51 (1971).
5. V. K. Bulgakov and A. M. Lipanov,  On the theory of
combustion of condensed materials with blowing, Fiz.
Goreniya Vzryva, 19, No. 3, 32 41 (1983).