Burning-Rate Models and Their Successors: 

A Personal Perspective 

by Martin S. Miller 

ARL-TR-2996 June 2003

Approved for public release; distribution is unlimited.
NOTICES 

Disclaimers 

The findings in this report are not to be construed as an official Department of the Army position unless 

so designated by other authorized documents. 

Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the

use thereof. 

Destroy this report when it is no longer needed. Do not return it to the originator. 
Army Research Laboratory 

Aberdeen Proving Ground, MD 21005-5066 

ARL-TR-2996 June 2003

Burning-Rate Models and Their Successors: 

A Personal Perspective 

Martin S. Miller 

Weapons and Materials Research Directorate, ARL 

Approved for public release; distribution is unlimited.
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering

and maintaining the data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, 

including suggestions for reducing the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson 

Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to

comply with a collection of information if it does not display a currently valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

June 2003

2. REPORT TYPE 

Final 

3. DATES COVERED (From - To) 

1 September 2002–25 March 2003

5a. CONTRACT NUMBER

5b. GRANT NUMBER

4. TITLE AND SUBTITLE 

Burning-Rate Models and Their Successors: A Personal Perspective 

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER 

622618.H80

5e. TASK NUMBER

6. AUTHOR(S) 

Martin S. Miller 

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

U.S. Army Research Laboratory 

ATTN: AMSRL-WM-BD 

Aberdeen Proving Ground, MD 21005-5066 

8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

ARL-TR-2996 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S)

11. SPONSOR/MONITOR'S REPORT NUMBER(S)

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Approved for public release; distribution is unlimited. 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

The development of models to compute the burning rate of energetic materials has undergone a paradigm shift in the last 15

years from treating the gas-phase reactions as a few global reactions to attempting to describe all of the elementary reactions

involved. This shift toward chemical specificity has brought the promise of true predictive capability in such models. 

Thwarting this promise, however, has been the unattainability of treating the condensed-phase and surface processes with 

commensurate detail and assurance. In this work, it is shown that limited predictability has recently been achieved by a semiempirical finesse of the condensed-phase reactions and surface gasification through the agency of a universal pyrolysis law for

certain families of propellant ingredients. It is also argued that the long-range prognosis for generalized predictability arising 

from a first-principles three-phase model lies with developing submodels for the condensed-phase and surface phenomena 

based on molecular dynamics (MD). However, since a continuum description of the gas-phase processes will always be

preferable, the mating of these MD submodels with the continuum gas-phase description will be a non-trivial research issue. 

To facilitate the coming dialogue between current continuum modelers and future MD modelers, this work describes the 

conceptual and mathematical basis of the continuum models in considerable detail. It is suggested that the mating of these 

disparate descriptions might be accomplished by developing simple, idealized continuum submodels whose parameters and 

assumptions might be justified by the MD submodels. An example of such a continuum submodel is developed for the process 

of multicomponent evaporation as a surface-gasification mechanism.

15. SUBJECT TERMS

propellant, burning rate, model, molecular dynamics 

16. SECURITY CLASSIFICATION OF: 19a. NAME OF RESPONSIBLE PERSON

Martin S. Miller a. REPORT 

UNCLASSIFIED 

b. ABSTRACT

UNCLASSIFIED

c. THIS PAGE 

UNCLASSIFIED

17. LIMITATION 

OF ABSTRACT 

UL 

18. NUMBER 

OF PAGES 

60 

19b. TELEPHONE NUMBER (Include area code) 

410-306-0718 

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18

ii
Contents 

List of Figures iv

List of Tables vi

1. Introduction 1

2. Phenomena 2

3. Concepts 3

3.1 Intraphase Conservation Equations.................................................................................5

3.1.1 Solid Phase ............................................................................................................5

3.1.2 Liquid Phase ..........................................................................................................6 

3.1.3 Gas Phase ..............................................................................................................7

3.2 Phase-Matching Continuity Conditions ..........................................................................7 

3.2.1 Species-Flux Continuity at Solid/Liquid Boundary, Liq x = −x .............................8 

3.2.2 Energy-Flux Continuity at Solid/Liquid Boundary, Liq x = −x ..............................9 

3.2.3 Species-Flux Continuity at Liquid/Gas Boundary, x = 0 .....................................9 

3.2.4 Energy-Flux Continuity at Liquid/Gas Boundary, x = 0 ......................................9 

3.3 Surface Regression Mechanism ......................................................................................9 

3.3.1 Single-Component Evaporation Mechanism.........................................................9 

3.3.2 Multicomponent Evaporation Mechanism..........................................................10 

3.4 Mathematical Closure of the 3-Phase Problem.............................................................11 

4. Models 12

4.1 Frozen Ozone ................................................................................................................14 

4.2 Deficiencies in the Idealization .....................................................................................17 

4.2.1 Multicomponent Evaporation..............................................................................17 

4.2.2 Liquid-Phase Diffusion .......................................................................................19 

4.2.3 Real-Gas Equation of State .................................................................................20 

4.2.4 Phase Separation..................................................................................................20 

4.3 RDX...............................................................................................................................20 

4.4 Multi-Ingredient Propellant Mixtures: A New Approach ............................................24 

iii
5. Challenges and Opportunities 38

5.1 “Molecular” Continuum Model of Multicomponent Evaporation................................38 

5.2 MD Simulations of the Condensed Phase .....................................................................44 

6. Conclusions 45

7. References 47

Appendix. Propellant Ingredients and Formulations 51

List of Figures 

Figure 1. Photograph of an (a) RDX-composite propellant (M43) deflagrating at 1.6 MPa 

and (b) a caricature of the 3-phase molecular processes involved.............................................4 

Figure 2. Energy fluxes at the phase boundaries arising from convection, molecular 

diffusion, and thermal conduction. ............................................................................................8

Figure 3. Single-component evaporative surface-regression mechanism.....................................10 

Figure 4. Flow chart illustrating the logic for determining the eigenvalues for the 3-phase 

problem for single-component evaporation.............................................................................13 

Figure 5. Computed eigenvalues and profiles for the steady-state deflagration of frozen 

ozone at 0.1 MPa and an initial temperature of 40 K. .............................................................15 

Figure 6. Rates of production of each species in a pure frozen-ozone problem source at 

0.1 MPa and 40 K initial temperature. The f and b suffixes indicate that the reaction 

identified is predominately proceeding in the forward and backward directions, 

respectively. The whole flame is divided into zones based on the predominant 

reactions there. .........................................................................................................................17 

Figure 7. Muticomponent mixture of O3 and O2 illustrating how the differing rates of

evaporation of the two components result in depletion of O2 mole fraction and 

enrichment of the O3 mole fraction at the surface necessitating the consideration of 

liquid-phase molecular diffusion to ensure species continuity................................................19 

Figure 8. Several 3-phase models compared to experimental RDX burning rates.......................23 

Figure 9. Zenin’s universal pyrolysis law with sampling of his experimental data (23) for 

a wide range of double-base propellant ingredient proportions...............................................25

Figure 10. Zenin’s data for a wide variety of nitramine/binder combinations (31–33). 

The pyrolysis law identified as NTRB is a least-squares fit to all of the RDX/ and 

HMX/binder data. The RDXBA pyrolysis law is a fit to the RDX/BAMO-AMMO 

data alone and the CL20PUNE pyrolysis law is the least-squares-fit to the 

CL20/PUNE data alone............................................................................................................27 

iv
Figure 11. Pyrolysis laws for a number of neat energetic materials. The RDX Miller 

calculation was performed with the 3-phase model previously used for frozen 

ozone using the same input data for RDX as Liau and Yang (with the exception of 

no condensed-phase reactions). All data and calculations are for an initial 

temperature of 293 K except for frozen ozone, which was at 40 K.........................................28 

Figure 12. Triply nitrated NC repeat unit. ....................................................................................30 

Figure 13. Monte Carlo model for distribution of repeat units among cellulose tri-, di-, 

and mononitrates for an NC specimen of given percent nitrogen. Comparison is made 

of the model predictions with the NMR data of Todd and Glasser (35)..................................31

Figure 14. Conceptual deconstruction of a propellant containing NC into subingredients 

and then into net condensed-phase decomposition products entering the gas phase...............32 

Figure 15. Comparison of CYCLOPS-code calculations of burning rate of M10 

propellant with experimental data............................................................................................33

Figure 16. Comparison of CYCLOPS-code calculations of burning rate of JA2 

propellant with experimental data............................................................................................34

Figure 17. Comparison of CYCLOPS-code predictions of dark-zone thermal structure 

compared with the experimental data of Vanderhoff et al. (42)..............................................35 

Figure 18. Comparison of CYCLOPS-code calculations of the burning rate of and 

RDX/thermoplastic elastomer (TPE) propellant using two different pyrolysis laws 

with the experimental data of Zenin (32).................................................................................36 

Figure 19. Comparison of CYCLOPS-code predictions of dark-zone thermal structure 

for an RDX/TPE propellant with the microthermocouple data of Zenin (32).........................36 

Figure 20. Model for continuous-phase molecular forces experienced by a molecule 

evaporating from a liquid surface. ...........................................................................................39

Figure 21. Accuracy of the simple heat-of-vaporization theory using 61 L-J parameter 

sets for both polar and nonpolar molecules. ............................................................................41 

Figure 22. Accuracy of the simple vapor-pressure model using 61 L-J parameter sets for 

both polar and nonpolar molecules..........................................................................................43 

Figure 23. One view of a relaxed configuration of JA2 propellant computed by a MD 

simulation (36) consisting of 2 chains of 15 monomers representing NC, 16 molecules 

of diethylene glycol dinitrate (DEGDN), and 10 molecules of NG. Oxygen atoms are 

in red, nitrogen in blue, carbon in grey, and hydrogen in white. A number of the 

component molecules can be identified as indicated. This is believed to be the first 

computed molecular representation of a real propellant formulation......................................44 

v
List of Tables 

Table 1. Condensed-phase reaction-rate coefficient parameters assumed by different 

models for RDX in the Arrhenius form. ..................................................................................22 

Table 2. Comparisons among RDX models of selected features at 1 atm....................................24 

Table 3. Zenin’s nitramine/binder test-material formulations (31–33). .......................................26 

Table 4. Parameters for various pyrolysis laws, equation 36. ......................................................28 

Table 5. Comparison of major species mole fractions in the dark zone of double-base 

propellant with various experimental measurements...............................................................35 

vi
1. Introduction 

The capability to calculate the burning rate of propellants from their ingredients has long been 

recognized as desirable though, until recently, not realizable. Even today, this capability is limited 

to classes of energetic materials that have been fairly extensively studied experimentally. So far, 

the goal of computing all the properties of propellants with ingredients that have never actually 

been synthesized is still beyond our reach. However, considerable progress towards obtaining 

properties of notional materials such as heats of formation, density, and detonation sensitivity has 

been made in the last 10 years. Over the same period of time, fledgling ability to compute the 

burning rate for unstudied formulations within known classes of propellants has emerged. 

What is it that makes the burning-rate calculations so difficult? The combustion of energetic 

materials involves coupling physical and chemical phenomena in the condensed phase, the gas 

phase, and at the interface between the two. Our greatest knowledge and experience is centered in

the gas phase. We now know that many dozens of chemical species reacting by many hundreds of 

elementary reactions are at play there in concert with the physical processes of molecular 

diffusion, convection, and thermal conduction. Even by itself, the gas phase presents us with a 

daunting array of nonlinear phenomena to understand. Fortunately, general scientific progress 

over the last several decades has armed us with the experimental and theoretical tools to approach 

this problem in a systematic way. Though many uncertainties remain, the conceptual means, if not 

always the resources, are available to resolve them. Such is not the case in dealing with the 

condensed-phase and interfacial processes. In these cases, we truly do not know what we do not 

know. 

It would seem that the best hope for breaking this impasse lies in molecular dynamics (MD) 

simulation using reactive potentials, but this approach is fairly said to be in its infancy. Thus far, 

only continuum-mechanics models have been developed to describe the 2- or 3-phase combustion 

process, and arguments are easily mounted to suggest that one will always want to use a continuum

model to describe the gas phase. Yet, I believe that mating these two descriptions will not be a 

trivial task. Anticipation of the coming need for communication between practitioners of these 

two approaches is what motivates this discourse. Accordingly, emphasis is placed on the 

mathematical framework and concepts behind the continuum-mechanics approach, the history of 

model development, their current status, and unsolved problems. 

1
2. Phenomena 

Before approaching the issues surrounding the mathematical modeling of propellant burning rates, 

one should become familiar with the range of phenomenological behaviors associated with the 

combustion of propellants and propellant ingredients. Solid propellants are of two basic types: 

homogeneous and composite. Single- and double-base propellants are examples of homogeneous 

propellants. Homogeneous propellants are considered to be well mixed on a molecular scale, 

although this may not be strictly true because of the practical constraints on mixing. Single-base 

propellant consists almost exclusively of nitrocellulose (NC), which becomes increasingly less 

soluble in the manufacturing solvents as its nitration level increases above ~13%N, leading to a 

very viscous liquid during mixing and extrusion. My first experience in measuring burning rates 

was with 6-in-long strands of M10.*

 Based on the measurement technique and degree of control 

over known variables such as pressure, I was expecting a standard error of about 1%. Instead, the 

standard deviation in measured burning rates using dozens of 1-in-long specimens was as large as 

25%. These large deviations arose from inhomogeneities in the propellant material itself due to 

imperfect mixing of the highly viscous feedstock material. This was a most inauspicious 

introduction to my new field! Fortunately, this would become the worst case I would ever 

encounter in my career. Most other homogeneous propellants benefit from the plasticizing 

properties of nitroglycerin (NG) and burn very reproducibly. However, one always must be aware, 

particularly with experimental propellants, that manufacturing procedures may influence the 

combustion properties in unexpected and nonreproducible ways. 

Composite propellants comprise a heterogeneous mixture of crystalline oxidizers and 

polymeric binders. Double-base propellant with crystalline nitroguanidine (NQ) added 

(M30) is an example. Only in the last 15 years have composite propellants involving 

cyclotrimethylenetrinitramine (RDX) come into use for guns. The use of such a secondary 

explosive as a major ingredient in gun propellant required considerable testing to pass safety 

criteria for use in the field. Composites utilizing RDX, cyclotetramethylenetetranitramine (HMX), 

and 2,4,6,8,10,12-hexanitrohexaazaisowurtzitane (CL20) coupled with various energetic polymeric 

binders are currently under active investigation. 

The propellant type can have important consequences for modeling the burning rate. In many 

cases, it is found that the oxidizer particle size strongly influences the burning rate. If this is so, a 

one-dimensional (1-D) model may not be suitable. On the other hand, observations of burning 

surfaces of nitramine-composite propellants (1) at pressures on the order of 1 MPa indicate that a 

melt layer exists at the surface. This melt layer provides an opportunity for the solid ingredients to 

 *

 Propellant formulations are given in the Appendix. 

2
become intimately mixed prior to gasification, potentially restoring a 1-D character to the 

combustion phenomena. 

Despite their uniform morphology, one cannot necessarily assume that homogeneous propellants 

burn one-dimensionally. The following description of double-base propellant combustion is 

enough to give pause to the most intrepid model builder: “It can be observed visually that the 

burning surface exhibits a wave-like mode of consumption. It appears as if glowing filaments of 

carbonaceous material periodically move over the surface consuming a thin layer of propellant. 

Ciné photography of the propellant surface shows that a smooth area on this surface appears to 

darken and to roughen. This area is then consumed by a wave of combustion which moves across 

the surface and leaves behind a network of carbon filaments which are blown off the surface by the 

steady evolution of gas. ... The consumption of the reacting surface layer leaves a smooth surface 

which subsequently repeats this sequence. It seems probable that the overall, average rate of 

burning is the result of two components: a steady rate analogous to the burning of a liquid 

propellant and a surface or condensed phase reactive wave moving laterally across the surface” 

(2). It is clear that any 1-D model of homogeneous or composite propellant combustion is to be 

understood as an idealization of the phenomena and not as an exact description. 

Despite all of the aforementioned complexities, propellants burn, at least on a macroscopic scale,

in “parallel layers” (i.e., in a direction normal to the local surface curvature). This is the property 

that allows for interior-ballistic control of the net gasification rate through intricate and ingenious 

propellant-grain geometries. It is also the property that gives the model builder some basis for 

hope. Another source of encouragement is that most propellants burn with a very simple powerlaw pressure dependence, despite radical changes in the appearance of the visible flame attached to

the burning surface. At pressures below ~1 MPa, there is no visible flame above the burning 

surface of homogeneous and many composite propellants. As the pressure increases, a weak flame 

appears ≥1 cm distant from the surface. With further pressure increases, the visible flame resides 

closer and closer to the surface, under steady-state, constant-pressure conditions. At ≥10 MPa, the 

flame appears attached directly to the surface. In most cases, this evolving flame behavior does 

not perturb the power-law pressure dependence of the burning rate. The nonluminous zone 

between the visible flame and the surface is known as the dark zone, which is now known to be a 

consequence of the relatively slow reduction of NO to N2 (and slow reactions involving HCN in 

the case of nitramines). All this suggests that burning-rate control resides in a thin gas layer close 

to the surface, an inference that further encourages the 1-D idealization. 

3. Concepts 

To date, there have been no MD models of multiphase combustion. This may change in the 

distant future; however, a nearer-term prospect is that MD submodels of the condensed-phase 

and/or interface processes will be developed and that these submodels will be used to supplement 

3
the continuum-mechanics combustion models. Anticipating the need to merge these two 

viewpoints, we develop in this section the mathematical formulation of the continuum paradigm. 

Studying Figure 1 will provide the reader with some appreciation of the range of phenomena that 

are involved in the combustion of propellants. No model exists which has treated all of the 

processes suggested in the figure; the list is intended as a conceptual transition from nature to 

mathematics. Some propellants are known to exhibit a liquid layer at the burning surface and 

some may not. I shall assume that, in the general case, a 3-phase (solid, liquid, gas) problem must 

be solved. The photograph in Figure 1 is of an RDX composite propellant (M43) burning in the 

steady state at a pressure of 1.6 MPa. In this case, the dark zone is clearly defined, and one can see 

that the idealization of a 1-D semi-infinite solid does not appear to be unreasonable. 

(b) (a)

Figure 1. Photograph of an (a) RDX-composite propellant (M43) deflagrating at 1.6 MPa and (b) a caricature of the 

3-phase molecular processes involved. 

4
Getting now to the formalism, I consider a semi-infinite solid combusting in the steady state at 

constant pressure. The spatial coordinate, taken as x, extends from −∞ deep in the solid at an 

initial temperature of T0, to the solid/liquid interface at − xliq, to the liquid/gas surface at x = 0, and 

finally to the region of equilibrium gas products at the adiabatic flame temperature Tf at x = ∞ . 

Equations conserving mass, atomic species, and energy must be solved in each phase subject to the 

boundary conditions for that phase. The intraphase solutions must also satisfy the equations of 

continuity at each phase boundary and whatever additional constraints are imposed by the surfaceregression mechanism.

3.1 Intraphase Conservation Equations 

The conservation equations within each phase are developed and discussed elsewhere (3), but will 

be summarized for convenience and completeness here. The + and − superscripts on the location 

superscripts indicate the side of the boundary where the values are taken, e.g., means

evaluated at the liquid side of the interface at the coordinate 

+ − Liq x

Liq x = −x . 

3.1.1 Solid Phase 

Species conservation: 

W k N

dx

dY

m k k

k & = ω& = 1,2,..., , (1)

assuming molecular diffusion is negligible on a combustion timescale, where is the total mass 

flux and an eigenvalue for the problem, Y is the mass fraction of the k

m&

k

th species, ωk & is the net 

rate of production of species k due to chemical reactions, Wk is the molecular weight of the kth

species, and N is the total number of distinct species in all three phases. These equations are 

subject to the domain boundary conditions 

= = = −∞ . (2) −∞ Y Y k 1,2,...,N at x k k

The set of mass fractions { } −∞ Yk reflect the composition of the unreacted propellant. 

Energy conservation: 

∑=

 − − = 

 



 NSol

k

Sol Sol k k k W h

dx

dT

mc

dx

dT

dx

d

1

λ & ω& 0 . (3)

Both the thermal conductivity λSol and the average specific heat of the solid mixture Sol c are, in 

general, functions of the independent variable x. W is the molecular weight and is the 

enthalpy of species k. This equation is subject to the domain boundary conditions: 

k hk

T = T0 at x = −∞ (4)

and 

5
m Liq T = T at x = − x , (5)

where T0 is the initial temperature of the propellant, Tm is the melting point of the propellant 

compound, and −xLiq is the coordinate of the solid/liquid boundary, which will be an eigenvalue of

the complete combustion problem. It should be noted here that these conservation equations for 

the solid phase and liquid phase in the next section embody precepts deriving from long experience 

with low-pressure gas-phase processes and may not be as general as presumed. In the gas phase at 

sufficiently low densities each reaction event takes place in essential isolation. Thermal reaction 

coefficients, obtained by averaging over velocities and reaction cross sections, can therefore be 

used to characterize reaction events anywhere and everywhere. In the condensed phase, on the 

other hand, reactive events do not take place in isolation but in a dielectric field produced by the 

close proximity of many “spectator” molecules. This dielectric field will evolve as a function of 

the instantaneous configuration of all the molecules in the vicinity of the reacting molecules during

the course of the reaction. The potential therefore exists that there will be too many variables to 

make the notion of thermal rate coefficients an adequate and useful idealization. Of course, at 

present the issue is moot in the case of energetic materials as the reactions are virtually unknown. 

An MD approach, in principle, would handle these complexities in a natural way. 

3.1.2 Liquid Phase 

All models to date have neglected molecular diffusion in the liquid phase as well, but I retain it in 

the general formulation here because I later show that it may be of some importance. 

Species conservation in liquid phase: 

( ) W k N

dx

dY Y V m

dx

d

k k

k

Liq k k ρ + & −ω& = 0 = 1,2,..., , (6) 

where ρ Liq is the liquid density and V is the diffusion velocity of species k. These equations are 

subject to the boundary conditions: 

k

 (7) Liq

x

k

x

k

x

Liq k

x

k mY Y V mY k 1,2,...,N at x x Liq Liq Liq Liq + = = = − + + + − − − − − & ρ &

and 

Y Y k 1,2,...,N at x 0 . (8) 0

k = k = = −

Energy conservation in liquid phase: 

∑ ∑=

 − − − = 

 



 N

k

k k k

N

k

k

Liq k k p

Liq

Liq p Y V c W h

dx

dT

mc

dx

dT

dx

d

1

λ & ρ ω& 0 , (9) 

subject to the boundary conditions: 

m Liq T = T at x = −x (10)

6
and 

 0 T = Ts at x = . (11)

3.1.3 Gas Phase 

Species conservation in gas phase: 

( ) W 0 k 1,2,..., N. dx

dY Y V m

dx

d

k k

k ρ Gas k k + & −ω& = = (12) 

subject to the boundary conditions: 

mY Y V mY Y V k 1,2,...,N at x 0 (13) 0

k

0

Gas k

0

k

0

k

0

Liq k

0

k + = + = = − − − + + +

& ρ & ρ

and 

= 0 k = 1,2,..., N at x = ∞. dx

dYk (14)

Energy conservation in gas phase: 

0, 1

∑ ∑=

 − − − = 

 



 N

k

k k k

N

k

k

Gas k k p

Gas

Gas p Y V c W h

dx

dT

mc

dx

dT

dx

d λ & ρ ω& (15) 

subject to the boundary conditions: 

T = Ts at x = 0 (16)

and 

= 0 at x = ∞. dx

dT (17)

There is, finally, the mass conservation equation

= 0, dx

dm& (18)

which has the trivial solution 

tan , (19) 0m cons t r = = ρ s &

through all three phases. In particular, at T0 , the solid-phase density is and the linear burning 

rate is r. 

0 ρ s

3.2 Phase-Matching Continuity Conditions 

Some of the above boundary conditions are expressed in quantities that are unknown at the outset 

and coupled to solutions in adjacent domains. These initially unknown quantities, m& , T , , s Liq x

7
and { }− 0 Yk are the eigenvalues for the complete problem. The final solutions for the temperature 

and mass fractions in each domain must satisfy the conservation equations with their respective 

boundary conditions in each phase for unique values of the eigenvalues. In order to accomplish 

this, one needs further constraints on the problem. These are to be found in the inter-phase 

equations of continuity, discussed in this section, and in the surface-regression mechanism,

discussed in the next section. Figure 2 shows the contributions to the energy fluxes across the 

phase boundaries. The meaning of the subscripts and superscripts will be clear from the figure 

contexts. The continuity conditions at each phase interface are constructed by equating the 

summed contributions on each side of a particular interface. In the figure, I have neglected 

contributions from the kinetic energy, as these are very small for typical propellant burning rates 

(of order 0.01% of the starting enthalpies), and also contributions from molecular diffusion in the 

solid phase, because it will be too slow on the time scale of importance to combustion. The 

boundary conditions are shown in Figure 2. 

Conduction 

Diffusion 

Convection 

Conduction 

Diffusion 

Convection 

Figure 2. Energy fluxes at the phase boundaries arising from convection, molecular diffusion, and thermal

conduction. 

3.2.1 Species-Flux Continuity at Solid/Liquid Boundary, Liq x = −x

This relation is the same as equation 7 and thus is not a new constraint. 

8
3.2.2 Energy-Flux Continuity at Solid/Liquid Boundary, Liq x = −x

∑ ∑ ∑ + + + + +

+

− −

−

− − − − −

−

− −

−

 + + 

 



  + = − 

 



 −

N

k

x

k

x

k

x

Liq k

N

k

x

k

x

k

x

Liq

N

k

x

k

x

k

x

Sol

Liq Liq Liq Liq Liq

Liq

Liq Liq

Liq

m Y h Y V h

dx

dT

m Y h

dx

dT λ & λ & ρ . (20) 

3.2.3 Species-Flux Continuity at Liquid/Gas Boundary, x = 0

This relation is the same as equation 13 and thus is not a new constraint. 

3.2.4 Energy-Flux Continuity at Liquid/Gas Boundary, x = 0

∑ ∑ ∑ ∑ + + + + +

+

− − − − −

−

 + + 

 



  + + = − 

 



 −

N

k

0

k

0

k

0

Gas k

N

k

0

k

0

k

0

Gas

N

k

0

k

0

k

0

Liq k

N

k

0

k

0

k

0

Liq m Y h Y V h

dx

dT

m Y h Y V h

dx

dT λ & ρ λ & ρ . (21) 

3.3 Surface Regression Mechanism 

In general, one might consider that many mechanisms contribute to surface regression during 

combustion: evaporation or desorption of surface species without change of identity, reaction of 

surface species to produce different gas-phase species, reaction of gas-phase species with surface 

species to produce other gas-phase species, and, possibly, physical ejection of molecular 

aggregates due to explosive reaction in the surface. It is also possible that some combination of 

these mechanisms might be active in the condensed phase below the surface, creating bubbles 

which are subsequently convected to the surface and released. Obviously, it will be a very great 

challenge to include all of these processes in a full combustion model, although attempts have been 

made to model RDX using a mix of evaporation and bubble formation (4, 5). In these treatments, 

one-dimensionality was preserved through the artifice of continuous porosity, though with 

unknown accuracy. In my opinion, the only surface-regression mechanism that has been treated 

with reasonable rigor is evaporation, and even this has been approximate with unknown accuracy. 

The remainder of this section will be devoted to a discussion of the approach taken to an 

evaporative mechanism and its potential shortcomings. 

3.3.1 Single-Component Evaporation Mechanism 

To date, the only evaporation mechanism used in published burning-rate model developments is 

based on the following reasoning. Consider a pure liquid substance in equilibrium with its vapor at 

some temperature Ts. I will refer to the multiphase presence of a single chemical species as “single 

component.” At equilibrium, the mass flux impinging on the surface is given by the kinetic-theory

result , 4

1

nvW where n is the molar number density in the vapor, v is the mean molecular velocity 

at Ts, and W is the molecular weight of the species being considered. If α is the fraction of 

impinging molecules that are absorbed into the surface (i.e., the accommodation coefficient), then 

the mass flux actually absorbed into the surface is nvW

4

1

α . Under equilibrium conditions, the

mass flux escaping the surface equals the mass flux being absorbed. Thus, the escaping mass flux 

9
can be expressed in terms of the number density at equilibrium or, through the equation of state, 

the equilibrium vapor pressure pe

. Under combustion conditions, the vapor pressure at the surface 

is less than the equilibrium value due to depletion by reactions and molecular diffusion. Thus, one 

can express the net regression rate of the surface during combustion in terms of the equilibrium

vapor pressure at Ts, the total pressure, and the mole fraction of the mother-liquor species on the 

gas-phase side of the surface, + 0 X . These arguments are summarized graphically in Figure 3 and 

in equation 22, which expresses the net flux of evaporating molecules as the difference between

the gross escaping flux and the gross condensing flux, i.e., 

X ,Ts

0+

( ) [p ( ) T X p . 2 RT

W

m total

0

s

e 2

1

s

+

− 













 = π & α ] (22) 

Figure 3. Single-component evaporative surface-regression mechanism.

3.3.2 Multicomponent Evaporation Mechanism 

On the face of it, there is no problem with the rigor of the arguments just presented. The difficulty

comes when applying it to even the simplest combustion problem. Take the self-deflagration of 

pure frozen ozone, for instance. This case will be developed more fully later, but it serves to 

illustrate the problem with single-component evaporation. Under combustion conditions, by 

definition, there will be chemical species other than ozone in the gas phase. These will be products 

and intermediates of ozone decomposition, i.e., O2 and O. These other species will also impinge 

upon and be absorbed to some extent in the surface. Thus, the surface will have a multicomponent 

nature, and the equilibrium vapor pressure of pure ozone will no longer determine the mass flux of 

ozone escaping the surface but also be a function of the mole fractions and molecular properties of

the other molecules present. In addition, the mass flux leaving the surface will include

contributions from the surface-absorbed products and intermediates originating from the gas-phase 

10
reactions. In principle, then, the single-component evaporation mechanism is intrinsically 

inappropriate for a problem involving combustion. Of course, as an approximation, it may turn out 

to be useful, but its accuracy is impossible to assess a priori.

3.4 Mathematical Closure of the 3-Phase Problem 

To illustrate the final posing of the 3-phase mathematical problem, I adopt single-component 

evaporation as the surface-regression mechanism and assume that molecular diffusion in the liquid 

phase can be neglected. Neglect of molecular diffusion means that only one boundary condition in 

the liquid-phase domain is required and I choose this to be 

. (23) Liq

x

k k Y Y k 1,2,...,N at x x Liq = = = − −

Note that in the absence of diffusion, Y . This single boundary condition replaces 

equations 7 and 8. Neglect of diffusion, through replacement of equation 8, also greatly simplifies 

the problem by reducing the number of eigenvalues to three, , T , and . One begins by 

providing starting estimates for these three eigenvalues of the coupled boundary-value problems; 

call them , T , and . The solid phase conservation equations are then solved using the 

estimated values and . This solution will provide values of

− + − − = Liq Liq x

k

x

k Y

Liq x s m&

Liq ˆ

x s

ˆ

ˆ

mˆ&

Liq x mˆ& { } Liq −x Yk needed by the liquidphase boundary condition equation 23. After solving the liquid-phase conservation equations test 

to see if

[ ( ) ( )] m

Sol

m k

Liq

k

N

k

x

k

x x

Liq

x x

Sol m Y h T h T

dx

dT

dx

dT Liq

Liq Liq

 + − 

 



 



 



 ∑ −

=− =−

ˆ

ˆ ˆ

λ ? λ & . (24) 

This is equation 20 rewritten for the neglect of liquid-phase diffusion. The quantity in brackets is 

the latent heat of fusion for simple substances. If the equality is not met at some chosen level of 

approximation, then choose another value of and repeat until equation 24 is satisfied. When 

equation 24 is satisfied, solve the gas-phase conservation equations using the liquid-phase 

solutions just obtained for

Liq x

{ }− 0 Yk , then test to see if the surface regression mechanism is satisfied, 

i.e., if

( ) [ ( ) s total

e

s

s p T X p

RT

W

m X T + +− 













 0

2

1

0 ˆ ˆ 2

? ˆ , ˆ

π & α ]. (25) 

Normally, equation 25 will not be satisfied at this point. Choose another value of mˆ& and solve the 

solid, liquid, and gas phases again, iterating until both equations 24 and 25 are satisfied. Then test 

to see if

∑( ) ∑ + + − − + + +

− +

 + − + 

 



  −  



 −

N

k

Gas k k k

N

k

Liq Gas k k k k m Y h Y h Y V h

dx

dT

dx

dT 0 0 0 0 0 0 0

0 0

ˆ λ ? λ & ρ . (26) 

11
This is equation 21 rewritten for the neglect of liquid-phase diffusion. Normally, equation 26 will 

not be satisfied at this point, so choose another value of Ts

ˆ and iterate all of the above until 

equations 24–26 are all satisfied to the required degree of accuracy. At this point, one will have 

determined the eigenvalues of the problem, xLiq, , and . The whole process can be more easily 

visualized as a logic flow chart in Figure 4. Actually, the process is straightforward and successive 

guesses for the eigenvalues can be based on physical reasoning depending on the results of each 

energy-flux equality test so that convergence of each eigenvalue can be approached with 

monotonic decreasing error. The solution process is obviously more complicated if molecular 

diffusion in the liquid phase is not neglected. 

Ts m&

4. Models 

The quantitative modeling of the burning rate of solid energetic materials really began in the dark 

hours of World War II with, predominantly, the efforts of Parr and Crawford (6) at the University

of Minnesota and Rice and Ginell (7) at the University of North Carolina. These and related 

wartime works were published in a special issue of the (then) Journal of Physical and Colloid 

Chemistry. At that time, virtually all propellants, gun and rocket, used some combination of NC 

and NG as their energetic ingredient. Since that time, composite propellants consisting of 

ammonium perchlorate in rubber binders (e.g., in the space shuttle boosters) have dominated the 

solid rocket propellant applications with solid composites based on HMX in polymeric binders 

also in use. Fielded gun propellants are still dominated by nitrate-ester propellants although 

composites of these conventional propellants with NQ are common. The variety of propellant 

ingredients exacerbates the difficulties faced by combustion modelers both because of their 

different mechanisms and because of the paucity of detailed experimental data available for many 

of them. Recent trends are toward even more rapid proliferation of new chemical ingredients such

as oxetanes, with their functional-group tailorability, and azides, with their attractive 

environmental advantages. In addition to enhanced performance and safety, the new constraints of

minimal environmental impact in manufacture, use, and demilitarization are now driving concerns 

in the development of new propellants. Coincident with the emergence of many promising new 

energetic materials, with attendant dilution of experimental characterization, is the growing 

urgency for theoretical guidance in the formulation of propellants incorporating these materials. 

It is generally true that the higher the performance required of the weapon system, the smaller the 

margin of safety will be in the functioning of all the components of the system including the 

propellant. If having a burning-rate model was deemed important 50 years ago, it is considerably 

more so with today’s new mix of developmental constraints and advanced ingredients. 

In this section, my aim is not to give a comprehensive history of burning-rate modeling but, rather, 

to provide a sense of the conceptual development of modeling approaches. This background is 

12
Converged eigenvalues

xLiq, Ts, m 

.

Gas/liquid B.C.

satisfied? 

Surface regression

mech. satisfied? 

Solve conservation 

eqns. in gas phase 

Liquid/solid B.C.

satisfied? 

Solve conservation eqns.

in solid & liquid phases

Choose trial values

of xLiq, Ts, and m 

.

Figure 4. Flow chart illustrating the logic for determining the eigenvalues 

for the 3-phase problem for single-component evaporation.

essential to an assessment of future avenues of progress in the field. Three subtopics will be 

addressed: frozen ozone, RDX, and multi-ingredient propellants. This progression of increasing 

system complexity allows us to illustrate some of the detailed mechanistic challenges facing the 

model builder and some new approaches to realizing a workable tool for the propellant formulator. 

13
4.1 Frozen Ozone 

Frozen ozone is the simplest chemical system falling within the scope of 3-phase self-sustained 

deflagration. Though simple from a theoretical point of view, it is anything but a straightforward 

subject for experimental investigation. Its propensity to detonate is legendary and the attendant 

dangers have undoubtedly inhibited the kind of extensive measurements of burning rate that one 

would like for comparison with model outputs. On the other hand, a wealth of high-quality 

experimental data has been obtained on thermophysical properties such as specific heat, thermal 

conductivity, melting and boiling points, latent heats, reaction paths and rates, and equations of 

state. This comprehensiveness and reliability of the input database on frozen ozone, coupled with 

its simplicity, makes it an attractive subject for modeling despite the paucity of burning-rate data. 

Its conceptual simplicity encourages and enables a more thorough study of mechanisms than with 

any other energetic material. 

Frozen ozone melts at about 80 K and has a normal boiling point of 161.3 K. The rate of reaction 

in the condensed phase is known to be very slow compared to the time scale of self-sustained 

deflagration. Thus, all of the uncertainties of describing condensed-phase reactions are 

conveniently (and legitimately!) sidestepped. The gas-phase reaction mechanism is known with 

good confidence to consist of the following three reactions: 

 O3 + M ↔ O2 + O + M H 25.65 kcal / mol ; (I) 0 ∆ 298.15K = +

 O3 + O ↔ O2 + O2 H 93.41 kcal / mol ; (II) 0 ∆ 298.15K = −

O + O + M ↔ O2 + M H 119.06 kcal / mol . (III) 0 ∆ 298.15K = −

If one is considering a pure ozone feedstock, then equation 22 describes the mass flux resulting 

from evaporation. The only measurement of the linear deflagration rate of condensed-phase ozone 

was made by Streng (8) for a liquid mixture of 90% O3 with 10% O2. The total mass flux leaving 

the surface in this multicomponent case is 

( ) 3 3

0 m m 1 YO mO & = & − + & −

, (27)

an equation first used by Ben-Reuven and Caveny (9) in connection with an RDX evaporation 

model. m is the mass flux of O O3

& 3 alone. Assuming that Y (i.e., no liquid-phase reactions 

or molecular diffusion), this expression reduces to 

−∞ ≈

−

3 O3

0

O Y

. Y

m

m

3

3

O

O

−∞ ≈ & & (28)

This mass-fraction-equivalency assumption is not strictly true for a number of reasons to be 

discussed subsequently, but it may be an adequate approximation. may be computed using a 

modification of equation 22 in which the equilibrium vapor pressure term is replaced by the O

O3 m&

e p 3

14
partial pressure, which may be approximated by , an expression of Raoult’s law for ideal 

solutions. is the O

e

O

0

O3 3 X p

−

− 0

O3 X 3 mole fraction at the surface, which is assumed to be approximated by 

, consistent with equation 28. Of course, Raoult’s law is only an approximation that can 

sometimes result in considerable error. 

−∞

O3 X

An example of a calculation (10) of the burning rate of pure frozen ozone at 0.1 MPa and initial 

temperature 40 K is given in Figure 5. As indicated in the figure, the computed linear burning rate 

is 0.25 cm/s at 0.1 MPa and an initial temperature of 40 K. By comparison, pure RDX burns 

~1 order of magnitude more slowly at a sevenfold higher initial temperature! It is interesting that 

the surface temperature is ~3 K lower than the boiling point at this pressure. An earlier model (11) 

assumed that the surface temperature is equal to the boiling point; evidently, in this case, it is not a

bad assumption. Unlike any of the burning-rate calculations for more complex systems, such as 

RDX, the calculations by Miller (10) include the effects of thermal diffusion in the gas phase, 

although the burning rate is changed by only a few percent if this process is neglected. 

Figure 5. Computed eigenvalues and profiles for the steady-state deflagration 

of frozen ozone at 0.1 MPa and an initial temperature of 40 K. 

-50 -40 -30 -20 -10 0 10 20 30 40 50 60 70

Temperature ( K )

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

(K) = 158.06

FROZEN OZONE COMBUSTION AT 1 ATM & 40 K

r (cm/s) = 0.2526 x (micron) = 24.07 T liq s

1.0

0.9 O2

0.8 T

0.7

0.6

M

ole Fractio

n

GAS

0.5

SOLID

LIQUID

0.4

0.3 O3

0.2

0.1 O

0.0

Distance from Surface ( microns )

Source: (10) 

15
Another mechanism investigated by Miller (10) and unique to this ozone study, is the possibility of

a heterogeneous reaction in which an oxygen atom from the gas phase reacts with a surface ozone 

molecule resulting in two gas-phase O2 molecules: 

O(gas) + O3(liq) → 2 O2(gas) H 93.41 kcal / mol . (IV) 0 ∆ 298.15K = −

This reaction apparently has never been measured but seems not only plausible but probable. It 

was assumed that the probability is unity in order to determine the maximum effect. Surprisingly, 

the burning rate increased only 1% at 0.1 and 2.0 MPa, despite this highly exothermic reaction 

occurring at the surface, where it adds directly to the heat feedback. In addition to the enhanced 

heat feedback, this reaction contributes to the destruction of ozone on the surface and should 

thereby contribute directly to the regression rate. It turns out that the O atoms near the surface that 

are consumed by the heterogeneous reaction would be consumed anyway in the reaction zone 

between 1 and 10 μm from the surface in the gas phase, and reactions this close to the surface 

contribute their heat with high efficiency to the heat feedback (Figure 6). The analysis leading to 

this conclusion is based on an important concept in steady-state combustion that can be gained 

from a study of the previously mentioned conservation equations. Starting with equation 15 and 

assuming that the mixture thermal conductivity is constant and that the specific heats for all 

species are equal and constant, one can show (12) that 

( ) ∫

∞ −

=

 = 

 





+ 0 0

q x e dx

dx

dT x mc

x

Gas Gas

p

λ λ

&

, (29)

where q(x) is the net volumetric rate of heat release in the gas phase. This expression indicates that 

heat released within a characteristic transport distance, 

p

Gas

m& c

λ , will contribute to the heat feedback 

with high efficiency. Equation 29 is valuable in very complex reaction mechanisms to sort out 

which reactions are materially affecting the burning rate. The characteristic distance is ~10 μm in

Figures 5 and 6. 

The ozone problem is simple enough to achieve a complete and unambiguous analysis of the 

chemistry in each part of the flame. The monotonic behavior of the temperature profile in Figure 5 

belies the underlying complexity of parallel and sequential reactions, which switch directions 

during the course of the flame, as seen in Figure 6. A curiosity indicated in Figure 6 is that 

spatially the first reaction step in the ozone decomposition flame is the production of ozone by 

reaction I running in reverse. In more complex reaction mechanisms with dozens of species and 

hundreds of reactions, it is impossible to completely analyze the chemistry. In those cases sorting 

tools such PREAD (13), ChemPlot (14), and ELEMAP (15) coupled with considerable experience 

are essential to gain insights into the mechanisms affecting both the burning rate and the flame 

features such as the dark zone. 

16
Distance from Surface (microns)

0.1 1 10 100 1000

-4

-3

-2

-1

0

1

2

3

4

5

R(O3)

R(O)

R(O2)

ZONE A ZONE B ZONE C ZONE D

P = 1 atm T0 = 40 K

Gas/Surface Reaction Probability = 0

 Rxns

Ib & IIIf

Rxn IIf Rxn If Rxn IIIf

Source: (10) 

Figure 6. Rates of production of each species in a pure frozen-ozone

problem source at 0.1 MPa and 40 K initial temperature. 

The f and b suffixes indicate that the reaction identified is 

predominately proceeding in the forward and backward

directions, respectively. The whole flame is divided into

zones based on the predominant reactions there. 

The lone measurement of a condensed-phase burning rate involving ozone, to which I previously 

alluded, was for a liquid mixture of 90% O3 and 10% O2. The result was ~0.4 cm/s; our 

calculation, utilizing equation 28 and Raoult’s law, produced 0.30 cm/s. Sandri (11) has 

speculated that Streng’s measured rate is too high because of heating of the liquid by the flame via 

conduction through the containing-vessel walls and by radiation from the flame. Plausible as they 

might be, however, the degree to which errors from these sources affected the burning-rate 

measurement is not possible to determine. Thus, there may remain a significant error in our 

burning-rate model, and it is worthwhile to look for shortcomings in the idealization which might 

account for the discrepancy. If such shortcomings matter to the ozone case, they may well matter 

in more complex systems as well. 

4.2 Deficiencies in the Idealization 

4.2.1 Multicomponent Evaporation 

As pointed out in section 3.3.2, any combustion problem driven by an evaporative surfaceregression mechanism is multicomponent in essence since product species will coexist in the liquid 

17
surface either as a result of the presence of more than one ingredient, condensed-phase reaction or

molecular diffusion of gas-phase reaction products back to and adsorption onto the surface. This 

issue is most easily discussed within the context of the O3/O2 mixture problem. Let us assume that 

there will be no atomic oxygen in the liquid surface due to its high reactivity. The attractive 

intermolecular forces dominate the heats of vaporization, and the strength of that interaction is 

ordered as follows: O3-O3 > O3-O2 > O2-O2. The heat of vaporization of an O3 molecule from a 

surface of both O3 and O2 molecules will therefore be less than that from a surface of pure O3

because the latent heat is just a measure of the work required to remove an O3 from the surface. 

On the other hand, the heat of vaporization of an O2 molecule from the same surface mixture will 

be greater than that for O2 from a pure O2 surface. These modifications, in turn, alter the 

equilibrium partial pressures of both O3 and O2 outside the mixture surface, and the equilibrium

vapor pressure is a key determinant of the surface-regression rate by equation 22. Thus, O3 will 

escape at a faster rate and O2 at a slower rate from a surface mixture of the two substances than 

would be expected using the single-component evaporation formulation coupled with equation 28. 

Beyond the greater difficulty of computing the rates of escape of each of the molecules in a multicomponent description, there are now fundamental consequences to the posing of the mathematical 

problem. Though the mass fluxes of O3 and O2 leaving the surface are different, in order to 

maintain a steady state, the linear rates of regression of both species must be equal, i.e., 

O3 O2 r = r . (30)

Since 

(31) 3 3 3

0

O O O m Y ρr −

& =

and 

, (32) 2 2 2

0

O O O m Y ρr −

& =

then 

. 0 0

2

2

3

3

− = −

O

o

O

o

Y

m

Y

m& & (33)

I have assumed that the liquid surface is composed only of O3 and O2. This means that 

; (34) − −= − 0

O

0

O2 3 Y 1 Y

O3 m& and can be obtained from equation 22, modified as previously described, using the proper 

multicomponent equilibrium partial pressures for each species, of course. Ideally, one would also 

dispense with Raoult’s law by computing the partial pressures based, for example, on model 

potential-energy functions as described in section 5.1. However, I am now left with a new 

O2 m&

18
unknown or eigenvalue of the problem, Y , with equation 33 as a new constraint. I now must add 

a new nested loop to the flow chart in Figure 4. While this is a serious increase in computational 

complexity for the ozone problem, it becomes an overwhelming increase in complexity for a 

compound like RDX which has dozens of potential species in the surface. Taking the general case 

where there are n chemical species in the surface, the number of eigenvalues, and hence nested 

loops, will be (n + 2). Solution of a complex set of equations like this may be possible, but it will 

undoubtedly require a different approach than nested loops, possibly a global optimization scheme

in which trial values of all of the variables are changed after each iteration. 

− 0

O3

4.2.2 Liquid-Phase Diffusion 

In the O3/O2 mixture case, the inclusion of multicomponent evaporation will mean that the mass 

fractions of O3 and O2 at the liquid side of the surface are different from the proportions in the 

feedstock at −∞. This is illustrated in Figure 7. Since there are no reactions in the liquid phase (by 

justifiable assumption) the discontinuity in feedstock/surface mass fractions will be resolved in 

nature by molecular diffusion, which must be considered as an attendant consequence of treating 

multicomponent evaporation. Because diffusion is generally slow in liquids, it will manifest itself

in this case as very steep gradients in the species just below the surface. Obviously, a whole new 

class of supportive data will be required. 

O3

O2

LOW RATE 

OF ESCAPE

HIGH RATE

OF ESCAPE

Diffusional

Flow of O3

Diffusional

Flow of O2

LIQUID GAS

Figure 7. Muticomponent mixture of O3 and O2 illustrating how the differing

rates of evaporation of the two components result in depletion of O2

mole fraction and enrichment of the O3 mole fraction at the surface 

necessitating the consideration of liquid-phase molecular diffusion to

ensure species continuity. 

19
4.2.3 Real-Gas Equation of State 

Since our calculations of the burning rate of frozen ozone and the liquid O3/O2 mixtures assumed

an ideal-gas equation of state, it is worth considering if real-gas effects play a role. The critical 

pressure and temperature of ozone are 5.46 MPa and 261.1 K, so the shortfall in computed burning

rate at 0.1 MPa cannot be attributed to this source. However, calculations by Miller (10) for frozen 

ozone were performed for pressures up to 2 MPa. At 2 MPa and 40 K initial temperature, the 

computed surface temperature is 217 K. Using standard tables (16), one finds that ozone gas at the 

surface has a compressibility factor ( RT

pv , where v is the molar volume) of about 0.75. At these 

conditions, the attractive molecular forces are dominating and number densities are higher than 

ideal, speeding up the reactions. At the same time, for pressures near the critical point, the heat of

vaporization is decreasing rapidly; this means that a given amount of heat feedback will vaporize 

more surface molecules. Thus, it is possible that use of a real-gas equation of state could have 

significant impact on a burning-rate calculation. To our knowledge, there are no published 

burning-rate calculations using a real-gas equation of state. 

4.2.4 Phase Separation

Below ~93 K, mixtures of O3 and O2 separate into two phases, the upper one being O2 rich and the 

lower one O3 rich. The proportion in each phase depends on the temperature. At the boiling point 

of liquid O2, 90.2 K, for example, a starting mixture of 50% O3 by weight, O3 in the upper layer 

ultimately settles to 12% by weight and 38% by weight in the lower phase. Streng’s measurement 

of the burning rate of a 90% O3/10% O2 mixture was conducted in a 9-mm Pyrex* tube cooled on 

the outside by liquid oxygen. No mention is made of the presence of two phases during 

combustion of the mixture, but if the temperature were actually 90.2 K, two liquid phases may 

have been present. It is possible that this phenomenon could explain the discrepancy between his 

measured rate and our computed one, since the O2-rich upper layer might be expected to evaporate 

faster. On the other hand, this would lead to a progressively richer mixture as the O2 selectively 

escaped; in turn, the specimen would not burn at a steady rate. It also might have been the case 

that the mixture surface regressed too fast to establish this phase separation. Finally, the presence 

of the flame inside the tube may have warmed the liquid mixture a few Kelvins, reaching the point 

where O3 and O2 are miscible in all proportions. 

4.3 RDX 

By far, the most intensive efforts to compute burning rate based on elementary gas-phase reactions 

have been directed towards neat RDX. It is known that RDX chemically decomposes both in the 

solid and liquid phases, but the solid-phase reaction is slow and probably not relevant to 

combustion time scales (17). RDX also has a relatively high vapor pressure and was hypothesized 

to gasify during steady combustion by evaporation first by Ben-Reuven and Caveny (9), whose 

 *

 Pyrex is a registered trademark of Corning Inc. 

20
model employed global reactions in the gas phase. The first use of the evaporation mechanism in a 

model with elementary gas-phase reactions was by Melius (18) ,who also assumed the single 

condensed-phase reaction, 

 RDX(Liq) → 3 CH2O + 3 N2O. (V)

His calculated burning rates at 0.1 and 2 MPa are in excellent agreement with experimental data.

Other elementary gas-phase reaction RDX models combining evaporation and liquid-phase 

reactions include those of Liau and Yang (4), Davidson and Beckstead (5), and Prasad et al. (19). 

The last of these works, though sharing most of the same condensed- and gas-phase reactions as 

the first two, is a semi-empirical model, requiring the experimental value of the surface 

temperature to obtain the burning rate, and will not be discussed here. Both Liau and Yang and 

Davidson and Beckstead assume that reaction V occurs and one other liquid-phase RDX 

decomposition reaction. In the Liau and Yang model, this other reaction is 

 RDX(Liq) → 3 HCN +3/2 NO + 3/2 NO2 + 3/2 H2O, (VI) 

whereas in the Davidson and Beckstead model, it is 

 RDX(Liq) → 3 H2CN +3 NO2. (VII)

Both models assume that RDX can form bubbles in the liquid layer and that the following 

secondary reaction takes place in the bubbles: 

 NO2 + CH2O → NO + CO + H2O. (VIII)

The reaction parameters chosen by each of the authors are given in Table 1. 

Both models describe the melt region as 2-phase, consisting of liquid and bubbles; however, since 

both models are 1-D, this 2-phase character is treated in the conservation equations as a continuous 

“porosity,” φ, defined as the ratio of the cross-sectional area occupied by the gas to the total crosssectional area. The rate of evaporation into “bubbles,” or more accurately, into the porosity, is 

computed using equation 22 at the local subsurface temperature instead of the surface temperature 

and the “surface area” of bubbles defined by Liau and Yang (4) as 

( )

( ) ( ) 2

1 A 36 n 1 ,

2

1 A 36 n ,

1 / 3 2 / 3

bubbles

1 / 3 2 / 3

bubbles

= − >

= <

π φ φ

π φ φ

.

(35)

where n is the bubble number density. These equations were given without explanation, except to 

say that n was to be determined empirically. Davidson and Beckstead (5), who used these same

equations, evidently understood that n would be a constant, independent of pressure, and they 

chose 1 × 10 (13) as the appropriate value without further explanation. Liau and Yang provide no 

value for n used in their calculations.

21
Table 1. Condensed-phase reaction-rate coefficient parameters assumed by different models for RDX in

the Arrhenius form. 

Reaction V VI VII VIII

A 4.66 × 1018 — — —

B 0 — — — Melius (18) 

E 47.8 kcal/mol — — — 

A 6.0 × 1013 — 1.6 × 1017 802 

B 0 — 0 2.77

Liau and Yang (4) 

E 36 kcal/mol — 45 kcal/mol 13.73 kcal/mol 

A 4.88 × 1011 6.5 × 1014 — 802

B 0 0 — 2.77

Davidson and Beckstead (5) 

E 36 kcal/mol 45 kcal/mol — 13.73 kcal/mol 

Notes: 

Burning-rate calculations based on these models are given in Figure 8. 

The Miller (10) 3-phase model, assuming no condensed-phase reactions, is applied to RDX, and its results are also 

shown in Figure 8. 

A is in appropriate centimeter, mole, Kelvin, and seconds units. 

The consequences of these varied assumptions and data on the computed burning rate are shown in 

Figure 8. For comparison, I applied our 3-phase code used in the ozone case (10) to RDX using 

the same input data and reaction mechanism employed by Liau and Yang with the exception that 

no condensed-phase reactions were considered. This model thus assumes that surface regression is 

by evaporation at the surface only, driven by the heat feedback from reactions in the gas phase 

whereas in the Liau and Yang model some RDX decomposes in the condensed phase and some

RDX evaporates into the liquid-phase porosity (representing bubbles). The close agreement 

between all of these models and the experimental burning rates (20–26) is nothing short of 

astonishing. Given the varied inputs and assumptions, one is tempted to conclude that the single 

common element is the evaporative mechanism and that the condensed-phase reactions are not 

playing a significant role in the burning rate. However, Davidson and Beckstead use a 

substantially lower vapor pressure (vapor over solid) than that used by Liau and Yang (vapor over 

liquid) and find a low sensitivity of the burning rate to vapor pressure. Moreover, Davidson and 

Beckstead predict a fairly large amount of the RDX decomposes in the liquid state, unlike Liau 

and Yang (Table 2). Unfortunately, Melius did not report the vapor-pressure expression used in 

his model, and his results are not otherwise directly comparable to the other models because his 

model used an earlier gas-phase reaction mechanism. All the other models used a mechanism

developed by Yetter et al. (27) based on Melius’s (18) but extended with new reactions thought 

important. 

22
Pressure (atm)

1 10 100 1000

Burnin

g Ra

t

e (

cm/

s) 0.01

0.1

1

10

Zenin (1995)

Miller (1988)

Homan (1997)

Atwood (1998)

Ulas et al. (1995)

Zimmer-Galler (1968)

Glaskova (1974)

Melius (1990)

Liau & Yang (1995)

Davidson & Beckstead (1997)

Miller Ozone Code with Liau & Yang Inputs

RDX

Figure 8. Several 3-phase models compared to experimental RDX burning

rates. 

In the final analysis, no definitive comparison is possible because of the incomplete reportage of

input data used in the Melius (18), Liau and Yang (4), and Davidson and Beckstead (5) models. 

Although each of the models based their solution of the gas-phase conservation equations on the 

PREMIX code (28) developed at Sandia National Laboratories, they each implemented the 3-

phase numerical solution in different ways. None of the authors described these methods in 

sufficient detail to enable others to reconstruct them unambiguously. One is left with the 

inescapable conclusion that the close agreement between models and experiment in Figure 8 is not 

as astonishing as it first seems but the result of finding a combination of input parameters that 

works. In the end, reconciliation of these various model assumptions, input data, and numerical 

approaches is probably not productive because of the fundamental and, so far, irreducible 

uncertainties associated with the condensed-phase processes, principally, but not limited to, the 

23
Table 2. Comparisons among RDX models of selected features at 1 atm.

Model Ts 

(K) 

Vapor Pressure 

(atm) at 600 K 

Sensitivity of 

Burning Rate to

Vapor Pressure 

RDX Fraction

Decomposed in

Liquid Phase

(%)

Melius (18) 549 Not reported Not reported 4.7 

Liau and Yang (4) 573 0.69 High <1

Davidson and Beckstead (5) 595 0.17 Low 25

Miller (10) 633 0.69 High 0

chemical reactions. Modeling the condensed phase in the kind of detail attempted by Liau and 

Yang and Davidson and Beckstead is not currently supported by the availability of adequate 

experimental measurements and may never be. 

4.4 Multi-Ingredient Propellant Mixtures: A New Approach 

In view of the situation previously described, how can one go beyond a relatively simple 

propellant ingredient like RDX, which has been the subject of the most extensive modern research 

to date, to treat multi-ingredient propellants? I asked this question of myself 6 years ago and in 

seeking an answer had an idea which became the basis for launching a new approach. The idea 

was to construct a hybrid-rigor model in which the gas phase would be treated in full elementaryreaction detail but the condensed phase and surface gasification would be treated in semi-empirical 

fashion using the Arrhenius-like expression first used by Rice and Ginell (7) and known in the 

propellant community as the “pyrolysis law”: 

m( ) T A e , s

s

RT

E

s s

−

& = (36)

coupled with a set of hypothesized decomposition products from each propellant ingredient. 

The decomposition products would conform to proper chemical balance and be selected according 

to available experimental results and theoretical reasoning. Furthermore, each propellant 

ingredient would be assumed to decompose without interference from other ingredients (i.e., 

noninteractively). A full development of this model for nitrate-ester propellants is given elsewhere 

(29, 30), but the basic concepts will be discussed here. 

Of course, if one had to measure a pyrolysis law for every combination of ingredients, there would 

be no hope of predictive capability from such a model. However, Anatoli Zenin, having spent the 

last 40+ years refining the technique of measuring surface temperatures of numerous propellant 

mixtures and single propellant ingredients with microthermocouples, discovered that a single, 

universal pyrolysis law (23) (Figure 9) holds for a wide range of double-base propellant 

24
Distance from Surface (microns)

0.1 1 10 100 1000

Species Production Rat

e (moles/

cc-s) -4

-3

-2

-1

0

1

2

3

4

5

R(O3)

R(O)

R(O2)

ZONE A ZONE B ZONE C ZONE D

P = 1 atm T0 = 40 K

Gas/Surface Reaction Probability = 0

Gas-Phase Species-Production Rates

 in Frozen-Ozone Combustion

 Rxns

Ib & IIIf

Rxn IIf Rxn If Rxn IIIf

Figure 9. Zenin’s universal pyrolysis law with sampling of his experimental data

(23) for a wide range of double-base propellant ingredient proportions. 

formulations with NG content from 0% to 50% and NC with percent nitration from about 

11.5% to 13.5%. As seen in Figure 9, even extreme initial temperature data are reasonably well 

accommodated by the same pyrolysis law. Furthermore, the same double-base pyrolysis law also 

works for additions of HMX (23). I wondered if other classes of ingredients (e.g., nitraminebinder systems) might display this same universality. If so, this idea might lead to a predictive 

tool, at least for members of these classes. To test this hypothesis, Zenin (23) applied his 

embedded microthermocouple technique to a wide range of nitramine-oxidizer/polymeric-binder 

combinations (Table 3). Figure 10 shows the results of this work (31–33) to date. First of all, 

notice that there is little difference between RDX and HMX materials with the same binder. Also, 

most binders with RDX and HMX oxidizers group closely about the pyrolysis law labeled 

RDXBA, which was obtained as a least-squares fit to the RDX/poly 3,3-bis(azidomethyl) oxetane 

(BAMO)-poly 3-azidomethyl-3-methyl oxetane (AMMO) data alone. It is of great interest that the

CL20/PUNE results, while closely adhering to the form of the pyrolysis law, occupies a very 

different region on the graph. HMX/PUNE, on the other hand, is in excellent accord with the 

25
Table 3. Zenin’s nitramine/binder test-material formulations (31–33). 

Test-Mixture Designation 

Binder-Ingredient 

Proportions 

Oxidizer-Binder 

Proportions 

RDX (or HMX)

a

/BAMO-AMMO 50:50 80:20

RDX (or HMX)

a

/BAMO-THF 50:50 80:20

RDX (or HMX)/CBIH — 80:20 

RDX (or HMX)/GAP1U — 80:20 

RDX (or HMX)/GAP2 — 80:20 

CL20 (or HMX)/PUNE 20:80 70:30 

a

 50% by weight of particle sizes < 50 μs, 50% in range 150–300 μm, 99% purity of RDX and 

HMX.

Notes: 

THF = tetrahydrofuran. 

CBIH = copolymer of butadiene and isoprene with hydroxyl terminated groups. 

GAP1U = glycidyl azide-polyurethane copolymer. 

GAP2 = glycidyl azide polymer (molecular mass of 2000). 

RDXBA pyrolysis law. The reason for this difference is not known; perhaps, the rate-limiting step 

in the decomposition of CL20 is significantly different than for RDX and HMX. The tentative 

conclusion from these measurements is that propellants composed of RDX or HMX in both inert 

and energetic polymeric binders do indeed follow a universal pyrolysis law. CL20 clearly follows 

such a law for one energetic binder and may for other binders as well. Other binders are currently

being investigated. The parameters for these pyrolysis laws along with Zenin’s double-base law 

are given in Table 4. It would be of considerable practical interest to explore the theoretical basis 

for the existence of these universal pyrolysis laws; they seen to suggest a common rate-limiting 

step associated with the condensed phase and/or interfacial region for each family. An 

understanding of this behavior could well lead to a predictive capability for pyrolysis laws, thereby 

removing a large empirical component in the present model. 

By way of reinforcing the general applicability of the form of the pyrolysis law, equation 36, it is 

worth noting that neat RDX, HMX, and ammonium dinitramide (ADN) all may be well described 

by this same functional relationship, as can be seen in Figure 11. In addition to the experimental 

data, I have plotted in Figure 11 the results of a least-squares fitting of the pyrolysis law to results 

calculated using my 3-phase code for both frozen ozone and neat RDX. It can be seen that the 

relation between the calculated surface temperatures and burning rates are in excellent accord with 

the form of equation 36. In view of all the foregoing evidence, there is considerable support for 

the general applicability of the pyrolysis law, and this should encourage efforts to calculate the 

parameters from first principles, an important step in removing the empiricism from this modeling 

approach and opening the door to much wider prediction capability. 

26
1/Ts

(K)

0.0010 0.0012 0.0014 0.0016 0.0018 0.0020 0.0022

Mass Regression Rate (g/cm2-s) 0.01

0.1

1

10

RDX/BAMO-AMMO @ 20C

RDX/BAMO-AMMO @ -100C

RDX/BAMO-AMMO @ 90C

HMX/BAMO-AMMO @ 20C

RDX/BAMO-THF @ 20C

HMX/BAMO-THF @ 20C

RDX/CBIH @ 20C

HMX/CBIH @ 20C

RDX/GAP1U @ 20C

HMX/GAP1U @ 20C

RDX/GAP2 @ 20C

HMX/GAP2 @ 20C

HMX/PUNE @ 20C

CL20/PUNE @ 20C

CL20/PUNE Pyrolysis Law (CL20PUNE)

RDX/BAMO-AMMO Pyrolysis Law (RDXBA)

Nitramine/Binder Pyrolysis Law (NTRB)

RDXBA CL20PUNE

NTRB

Figure 10. Zenin’s data for a wide variety of nitramine/binder combinations (31–33). The pyrolysis law

identified as NTRB is a least-squares fit to all of the RDX/ and HMX/binder data. The 

RDXBA pyrolysis law is a fit to the RDX/BAMO-AMMO data alone, and the CL20PUNE 

pyrolysis law is the least-squares-fit to the CL20/PUNE data alone. 

27
Table 4. Parameters for various pyrolysis laws, equation 36. 

Designation As 

(g/cm2-s)

Es 

(cal/mol) 

DB 1800 9935

RDXBA 10470 1371

NTRB 2004 11827

CL20PUNE 1.868 × 106 16032 

1/Ts(K)

0.0010 0.0020 0.0030 0.0040 0.0050 0.0060

M

ass Regression Rate (g/

c

m2-s)

0.01

0.1

1

10

100

1000

10000

RDX Zenin Experimental Data 

RDX Zenin Data Fit

RDX Miller 3-Phase Model

RDX Miller 3-Phase Model Fit

HMX Zenin Experimental Data

HMX Zenin Data Fit

ADN Zenin Experimental Data

ADN Zenin Data Fit

O3 Miller 3-Phase Model

O3 Miller 3-Phase Model Fit

ADN Exp.

HMX Exp.

RDX Exp.

RDX Model

O3 Model

Figure 11. Pyrolysis laws for a number of neat energetic materials. The RDX Miller calculation was performed

with the 3-phase model previously used for frozen ozone using the same input data for RDX as Liau and 

Yang (with the exception of no condensed-phase reactions). All data and calculations are for an initial 

temperature of 293 K except for frozen ozone, which was at 40 K.

As previously indicated, the model requires the identity and mole fractions of the chemical species 

resulting from the condensed-phase processes (i.e., the gases leaving the surface). These species 

will be far from the equilibrium distribution, and, therefore, one has only chemical balance to 

constrain the set of surface products absolutely. With ingredients as complex as those for 

propellants, this means that there are usually very many possible product sets to consider. 

Experimental thermal decomposition studies may be useful guides to selecting an appropriate 

28
product set but not perfect since those experiments may not accurately mimic the conditions of the 

burning surface (temperature and heating rate) and probably all suffer from significant unknown 

changes in the surface products as a result of very fast reactions close to the surface. Use of 

general theoretical reasoning can also help limit the number of probable product sets. In the end, 

however, an important selection criterion must be whether the proposed product set results in a 

calculated burning rate close to the experimental one. Use of this criterion, of course, assumes that 

the pyrolysis law, the gas-phase reaction mechanism, and all the other input data the model needs

are accurate. The predictive ability of the model for the whole propellant is redeemed by insisting 

that the decomposition products from each ingredient, “calibrated” ideally for that ingredient 

alone, remain consistent for all propellant formulations using that ingredient. This strategy is 

probably the best-compromise approach, but one should recognize that the assumption of noninteractivity in decomposition might be a better approximation for some systems than others. 

Since all fielded gun propellants include at least some NC and certain propellants, such as M10, 

have 98% NC, any burning-rate model of practical importance will have to deal with NC as an 

ingredient. Not only is NC a complex long-chain polymer, but the repeat units vary among four 

different types depending on whether the glucose ring is triply, doubly, singly nitrated or 

unnitrated altogether. Figure 12 shows a triply nitrated repeat unit. NC is characterized for 

propellant use by its average percent nitrogen (%N) by weight, and this quantity varies typically 

from ~11.5 to ~13.5 %N. Pure cellulose trinitrate, cellulose dinitrate, and cellulose mononitrate 

have 14.14, 11.11, and 6.76 %N, respectively. The work of Leider and Seaton (34) showed that a 

12.3 %N NC has all three nitration states present but no unnitrated sites. We developed a Monte 

Carlo model (29, 30) to determine for an NC specimen of given percent nitrogen what is the 

distribution among these three nitration states. Predictions of this model are compared to the 

nuclear magnetic resonance (NMR) data of Todd and Glasser (35) in Figure 13. The burning-rate 

model then assumes that NC consists of three separate, hypothetically pure ingredients 

corresponding to cellulose tri-, di-, and mononitrate in the proportion indicated by the Monte Carlo 

model. The densities of these hypothetically pure nitrates are obtained from a MD calculation (36) 

using the COMPASS force field. 

The model handling of multi-ingredient propellants is summarized graphically for JA2 in 

Figure 14. This figure does not represent the current decomposition-product sets but serves to 

illustrate the principles behind the method. The CYCLOPS (29, 30) code, named for the mythical 

race of creatures who forged thunderbolts for Zeus, is the computer-code embodiment of the 

burning-rate model under discussion. CYCLOPS accepts the weight percents of each ingredient 

and checks to see if NC is among them. If so, then the Monte-Carlo subroutine is called and the 

nitrate state distribution determined. The ingredient list is then expanded to include the three NC 

subingredients (cellulose tri-, di-, and mononitrate). The condensed-phase decomposition 

29
O C

H

C C

C

C O

H H

CH2ONO2

O

H

O

H

NO2 NO2

6

3 2

Figure 12. Triply nitrated NC repeat unit. 

products for each of the ingredients on the expanded list are then obtained from the 

decomposition-products database, and the net set of products computed according to the 

proportion of each ingredient. A trial value for the linear burning rate is then read in from the 

problem input file and converted to mass flux using the code-computed propellant density. The 

surface temperature for this mass flux is then computed from the pyrolysis law, and these values 

are passed to a modified version of the PREMIX (28) code, which is called as a subroutine to solve 

the gas-phase conservation equations. From this solution the heat feedback to the surface is 

computed and compared with the heat feedback required to transform the propellant ingredients at 

their initial temperature to the condensed-phase decomposition products at the trial surface 

temperature using the trial mass flux. New trial values are automatically provided by the code 

until the two heat fluxes are in satisfactory agreement. The criterion for convergence is actually 

expressed in terms of acceptable changes in successive guesses of the mass flux, which is the 

unknown sought. 

30
12.6 %N

mono di tri

Fractio

n of Repeat U

nits

0.0

0.1

0.2

0.3

0.4

0.5

0.6

Monte-Carlo Code

T&G 1996 (P7 L-334)

13.4 %N

Substitution Level

mono di tri

Fractio

n of Repe

at Units

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Monte-Carlo Code

T&G 1996 (P1 L-333)

13.15 %N

mono di tri

Fractio

n of Repeat Units

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Monte-Carlo Code

T&G 1996 (L10683Y)

T&G 1996 (F10692Y)

T&G 1996 (W10673Y)

T&G: Todd & Glasser

Virginia Tech

NMR Data

Figure 13. Monte Carlo model for distribution of repeat units among cellulose tri-, di-, and

mononitrates for an NC specimen of given percent nitrogen. Comparison is made of the 

model predictions with the NMR data of Todd and Glasser (35). 

31
Figure 14. Conceptual deconstruction of a propellant containing NC into subingredients and then into net

condensed-phase decomposition products entering the gas phase. 

Results for the burning rate of M10 and JA2 are shown in Figures 15 and 16 compared to the 

experimental data of Miller (22), Juhasz (37), and Atwood et al. (38) for M10 and Miller (22) and

Juhasz et al. (39) for JA2. An extensive analysis (29, 30) of the chemical-kinetic origin of the 

inflection in the M10 curve has been made. It was discovered that those reactions controlling the 

dark-zone length at low pressures (<10 MPa) have little influence over the burning rate at low 

pressures, but, to our surprise, a major effect on the burning rate at high pressures, where a dark 

zone does not even exist. Thus, studies of the dark-zone chemistry experimentally accessible at 

low pressures actually are probing those reactions with critical influence over the burning rate at 

high pressures. This finding provides new impetus to experimental dark-zone investigations, 

especially for nitramine propellants because their dark-zone chemistry is less well known than that 

of nitrate-ester propellants. How well CYCLOPS predicts the structure of the dark zone is shown 

in Figure 17, where the temperature and NO profiles are seen to compare very favorably with 

experiment. Table 5 gives a comparison of predicted and measured major species for a doublebase propellant similar to M9; again, the agreement is excellent. The experimental measurements

used for comparison are those of Heller and Gordon (40), Lengelle et al. (41), and Vanderhoff et 

al. (42), and they were renormalized for direct comparability by Vanderhoff et al. (42). 

32
Pressure ( MPa )

0.1 1 10 100 1000

Linear Burning Rate ( cm/s )

0.1

1

10

100

Experiment, M10: Miller (1985)

Experiment, M10: AFATL (1982)

Experiment, M10: Atwood, et al. (1988)

CYCLOPS Code

M10

Figure 15. Comparison of CYCLOPS-code calculations of burning rate of M10 

propellant with experimental data.

A final example is given of more recent work (43) on RDX in an energetic thermoplastic elastomer 

(ETPE) binder. The excellent comparison with experimental burning-rate data shown in Figure 18 

is tempered by the poor prediction of the thermal structure in Figure 19. This inconsistency might 

at first be thought to suggest that the good burning-rate prediction is simply fortuitous, but this is 

not necessarily the case. As we learned in the case of M10, the burning rate and the dark-zone 

length at low pressure may be controlled by different sets of reactions, so it is possible that the 

good burning-rate prediction and poor thermal structure prediction may be due to imperfections in 

the chemical mechanism describing the nitramine dark-zone chemistry. In all propellants with 

RDX, we assume that the RDX evaporates unchanged, just as I had assumed for neat RDX; 

however, binder ingredients are assumed to decompose in the condensed phase. 

33
Pressure (MPa)

0.1 1 10 100 1000

Linea

r Burning Rate (cm/s

)

0.01

0.1

1

10

100

Strand Burner: Miller (1993)

200cc Closed Bomb: Juhasz (1999)

CYCLOPS Code Prediction

JA2

Figure 16. Comparison of CYCLOPS-code calculations of burning rate of JA2

propellant with experimental data. 

The foregoing examples give evidence of the promise of the CYCLOPS code in providing both 

reasonable predictions of burning rate as well as details of the gas-phase flame structure and 

detailed chemistry. At present, the code is limited to certain families of propellant ingredients, but 

the number of families will likely increase as other systems are studied. However, perhaps more 

important is the growing probability that one might be able to predict pyrolysis laws from first 

principles considerations. Certainly, greater theoretical capabilities to predict the final 

condensed-phase decomposition products of propellant ingredients can be brought to bear. Thus, 

the CYCLOPS code may well provide the best computational vehicle for incorporating results of 

new theoretical approaches to the condensed-phase and surface processes. How this might be 

accomplished while maintaining mathematical tractability of the code is suggested in section 5. 

34
Distance from Surface in Gas (cm)

0.01 0.1 1

Ga

s Temperature (K)

500

1000

1500

2000

2500

3000

CYCLOPS Code Prediction

Vanderhoff, et al. (1992) 

 by NO Absorption

Vanderhoff, et al. (1992)

by OH Absorption

JA2 @ 1.6 MPa

Figure 17. Comparison of CYCLOPS-code predictions of dark-zone thermal

structure compared with the experimental data of Vanderhoff et al. (42) 

Table 5. Comparison of major species mole fractions in the dark zone of double-base propellant with various 

experimental measurements. 

Parameters 

Heller and 

Gordona Lengelle et al.b Vanderhoff et al.c CYCLOPS 

(Present Model) 

P (MPa) 1.6 0.9 1.7 1.7 

Dark-zone temperature (K) 1600 1500 1500 1543 

NO 0.24 0.21 0.24 0.25

CO 0.33 0.38 — 0.32

H2 0.08 0.08 — 0.08

N2 0.04 0.02 — 0.004

H2O 0.20 0.20 — 0.19

CO2 0.10 0.09 — 0.10

HCN 0.004 — — 0.004

CH4 0.008 0.026 — 0.009

C2H4 0.008 0.008 — 0.001

a

 Source: (40). b

 Source: (41). c

 Source: (42).

35
RDX/TPE Propellant @ 293K:

80% RDX, 20% BAMO/AMMO (50/50 Mix)

Pressure (MPa)

0.1 1 10 100

Linear Burning Rate (cm/s)

0.01

0.1

1

10

Data of Zenin (2000)

CYCLOPS: with NTRB pyrolysis law

CYCLOPS: with RDX/BA pyrolysis law

Figure 18. Comparison of CYCLOPS-code calculations of the burning rate of and

RDX/thermoplastic elastomer (TPE) propellant using two different pyrolysis 

laws with the experimental data of Zenin (32). 

Distance from Surface (mm)

0.0 0.5 1.0 1.5 2.0

Tem

perature (C)

0

500

1000

1500

2000

2500

Zenin Embedded-Thermocouple Data

CYCLOPS: NTRB pyrolysis law

RDX/BAMO-AMMO

(10 atm, 20 C)

Figure 19. Comparison of CYCLOPS-code predictions of dark-zone thermal structure

for an RDX/TPE propellant with the microthermocouple data of Zenin (32).

36
In summary, the CYCLOPS code is based upon a few simple propositions: 

• The surface regression may be adequately idealized as 1-D. This condition is probably met 

for most homogeneous propellants. It may also be approximately true for many composite 

propellants where the burning rate is not a strong function of oxidizer particle size. Surface 

melt layers, for example, may provide effective premixing of ingredients where the oxidizer

particle sizes are not too large. 

• The surface regression may be described by a pyrolysis law, equation 36. I have shown that 

this condition has been met by a large number of neat energetic materials and mixtures of 

typical propellant ingredients. 

• The overall products of condensed-phase decomposition may be estimated with sufficient 

accuracy. Chemical balance and results of thermal decomposition experiments are sources 

of guidance here, though only the former can be considered as an absolute constraint because 

of the high probability of secondary reactions in decomposition experiments. 

• The decomposition of one propellant ingredient does not affect the decomposition of other 

ingredients, i.e., the decomposition of a multi-ingredient propellant may be described as the 

superposition of the independent decompositions of each of its ingredients. This may well be 

a better approximation for some ingredients than for others. Short of describing the mixture 

decomposition theoretically on a molecular scale, there is little recourse to this assumption. 

However, it is often found that significant changes in detailed combustion mechanisms have 

a remarkably weak influence on the burning rate, that quantity being a highly integrated 

consequence of myriad underlying details. CYCLOPS exploits this relative insensitivity to 

mechanistic details. 

Finally, since the CYCLOPS code is the first code of its kind (though still undergoing 

development to improve robustness for general use), one can anticipate ways in which it might be 

put to practical use even at its current stage of development. First, the CYCLOPS code can be 

used to optimize ingredient proportions to achieve a target burning rate. Also, CYCLOPS could 

be used to explore subtle, previously unexplained effects of formulation on burning rate. For 

example, German-made JA2 burns about 20% faster than U.S.-made JA2 with identical ingredient 

proportions. However, the two differ by different specifications on the NC component. While the 

average percent nitrogen for the two propellants is the same, the two materials are made from

blends of two lots of NC with different percent nitrogen. Because our nitrate-state-distribution 

Monte Carlo code predicts that the surface products for these two propellants will be different, it 

will be worthwhile to see if CYCLOPS can predict the burning-rate difference. Further, 

CYCLOPS could be used to help set manufacturing tolerances in the specifications of ingredient 

proportions in military propellants. The CYCLOPS code could be used to determine how any 

given set of tolerances will map into variations of the burning rate and these burning-rate 

variations could be judged by using them as inputs to interior-ballistic codes. The burning rate at

low pressure has been shown to be sensitive to the heats of formation of propellant ingredients in 

37
certain cases (29, 30). To the extent that ingredient purity is reflected in the heat of formation 

(e.g., due to dinitroglycerine impurity in trinitroglycerine), the CYCLOPS code could be used to 

determine how ingredient purity affects the burning rate, thereby enabling rationalization of 

quantitative tolerances for these ingredients. Finally, the CYCLOPS code could determine the 

effect of chemical modifiers on the burning rate without mixing chemicals. Though not replacing 

the need to mix and test, using the code as a screening agent and test bed for ideas could improve 

efficiency and possibly suggest new compounds to test. 

5. Challenges and Opportunities

The semi-empirical aspects of the CYCLOPS code need be tolerated only until theoretical 

advances obviate their necessity. The next level of improvement may well come from MD 

descriptions of the condensed phase and surface phenomena. However, this more fundamental 

level of treatment will have its own set of approximations and limitations. In order to deal with 

condensed-phase reactions, for example, considerable progress will have to be made in the 

parameterization of, and experience with, reactive force fields. In my judgment, one will never 

want to treat the gas phase, with its dozens of species and hundreds of reactions, in this way. To 

do so would discard more than 50 years of hard-won kinetics research gains. One will therefore be

faced with merging the discrete and continuous descriptions into a single, tractable mathematical 

entity. I believe that this will be no easy task, as the MD approach is inherently time-dependent 

with (computationally intensive) explicit statistical averaging and the steady-state continuummechanics approach outlined earlier in this work is inherently time-independent with implicit 

statistical averaging. The most fruitful approach may be to develop simple, idealized continuum

submodels calibrated to stand-alone MD models of a limited set of phenomena. For example, a 

reactive-MD model of the condensed phase would provide a set of condensed-phase 

decomposition products needed by the CYCLOPS code. Similarly, an MD model of

multicomponent evaporation and/or pyrolysis could be used to determine a pyrolysis law for 

CYCLOPS. Results of the continuum submodels, depending on their scope, could then be easily 

incorporated into either a 3-phase model or a code such as CYCLOPS, retaining mathematical 

tractability and efficiency. Since, prior to the development of CYCLOPS, I had made a start on 

incorporating multicomponent evaporation into our 3-phase code, it will be used as an example 

here of the proposed approach. 

5.1 “Molecular” Continuum Model of Multicomponent Evaporation 

I consider that a molecule evaporating from a liquid surface (Figure 20) is at some distance d

above the surface. I further assume that the molecule follows a tubular path from deep within the 

liquid to points outside the liquid; the tube has radius a0. The conception here is that the radius of

the tube is of the order of the “radius” of the escaping molecule, and what is a tortured path 

geometry in reality is idealized as straightened to a cylindrical path. Since the escaping molecule 

38
Figure 20. Model for continuous-phase molecular forces experienced by a molecule evaporating from a liquid surface. 

is at a distance r from all the molecules in the differential volume element, dV = 2π a da dz , the 

potential energy of interaction is dΦ = φ(r)n dV, where n is the number density of liquid 

molecules (Figure 20). If I assume a Lennard-Jones (L-J) interaction potential, where ε is the 

well depth and σ is the collision diameter, i.e., 

( ) 













 



 



  − 

 



 =

12 6

4

r r

r

σ σ φ ε , (37)

and integrate over all the molecules in the liquid, one obtains

( ) [ ] ( ) [ ( )]































 − 













 Φ = g x d

a

f x d

a

d nLiq

3

0

9

0

2

4

1

10

1 8 σ σ

π εσ , (38) 

where 

( )

0 a

d

x d = − , (39)

( ) ( ) ( ) ( ) ( ) ( ) 768

105

arctan x

384

105

384 1 x

945

192 1 x

315

48 1 x

63

8 1 x

9

9

x f x 2 2 2 3 2 4 2

π

+ +

















+ +

+

+

+

+ = , (40) 

and 

( ) ( ) ( ) 4

arctan x

2

1

2 1 x

x

g x 2

π

+ +

+ = . (41)

39
Two important extensions can be made to this formulation. The first is to add the effects of

interaction between the escaping molecule and the gas-phase molecules. This modifies the 

potential to 

( ) ( ) [ ] ( ) [ ( )]































 − 













 = − g x d

4 a

1 f x d

10 a

1 d 8 n n

3

0

9

0

2

Liq Gas

σ σ Φ π εσ . (42) 

The other important extension is to the multicomponent case. This requires a sum over the 

interactions between the escaping molecule of species i with the other species j, of total different 

kinds, N. Thus, the multicomponent version of the interaction potential between an evaporating 

molecule of species i at a distance d from the surface and all other species both in the liquid and in 

the gas phase is 

( ) ( ) [ ] ( ) [ ( )]































 − 













 = − ∑=

g x d

4 a

1 f x d

10 a

1 d 8 n n

3

0

ij

9

0

ij

N

j 1

2

i Liq Gas ij ij

σ σ

Φ π ε σ , (43) 

with its ancillary definitions in equations 39–41. These equations may seem awkward but are 

easily and quickly evaluated in a computer code and represent, in a very general way, enormously 

complex phenomena. For example, they describe quantitatively how the heat of vaporization goes 

to zero at the critical point, where the liquid- and gas-phase densities become indistinguishable. 

A necessary test of the model is the accuracy with which the previous equations predict the heats 

of vaporization of pure substances. The heat of vaporization is obtained by taking the limit in the

previous equations as d 64 in equation 38. The result is 































 − 













 =

6

0

3

0

3

Liq

2

vap 32 a

7 1

a

H n

σ σ

π εσ . (44)

Application of this model is made to 61 data sets for polar and nonpolar molecules using the 

compilation of L-J parameters in Reid and Sherwood (44) and Reid et al. (45). Calculation results 

shown in Figure 21 are obtained by performing a non-linear least-squares fit of the heat-ofevaporation experimental data to equation 44 using Pe, defined as follows, as an adjustable 

parameter to mediate between the raw theoretical result Hcalc and the best estimate value Hest; 

Hest = PeHcalc . (45)

Assuming that 



 





σ

a0 is unity, the best-fit value of turns out to be 1.10. Its proximity to unity 

suggests that the physical basis of model assumptions are reasonable. The standard deviation of 

the error using the optimized value of P

Pe

e is about ±12%. This error is surprisingly small 

considering that the L-J potential is not generally as faithful as an exponential-6 potential for 

nonpolar molecules and is has even worse fidelity for polar molecules, which are abundantly 

40
Hvap (kcal/mole)

0 1 234 5 6 789 10

% Er

r

or in Estimate

d Hvap

-40

-30

-20

-10

0

10

20

30

40

Nonpolar Molecules

Polar Molecules

I

2

CHCl3

CCl

4

SO2 C2H2

N2

Ar

O O2 2

C2H2

C3H8

C2H6

N2O

CO2

H2S

HI

Br2

CS2

HCN

n-C4H10

Cl2 HBr

HCl

C2H4

Xe

CO2

C2H6

C2H4

CH4

Kr

CH4

Ar

F2

F2

N2

Ne

− σ

+ σ

NO

CO

H2O

NH3

n-C3H7OH

C2H5OH

CH3OH

C3H6

C2H5Cl

O2

CH2CHCH3

CH3Br

CH3COCH3

CH3COOCH3

CH2Cl2 CH3COOC2H5 C6H6

n-C6H14

CH3OCH3

CH3Cl

iso-C4H10

CH3CCH n-C5H12

C6H12

C2H5OC2H5

Figure 21. Accuracy of the simple heat-of-vaporization theory using 61 L-J 

parameter sets for both polar and nonpolar molecules. 

represented in the data set. In fact, even without the adjustment parameter Pe, the standard 

deviation of the prediction of Hvap is about 16%. This good performance might be improved upon 

by either using another, better interaction potential or by possibly modifying the model to 

recognize the discrete nature of close-neighbor molecular interactions. 

By itself, the previously mentioned model is not sufficient to treat the multicomponent evaporation 

problem. One needs also to describe the dynamics of the evaporative escape mechanism. 

Drawing upon kinetic theory to construct the flux of molecules leaving the liquid surface, Γout , 

one obtains 

e

RT

H

out Liq n e v

4

1 s

vap















 = −

Γ , (46)

where the quantity in parenthesis is the density of molecules with energies greater than Hvap, and 

e v is the average velocity of those molecules with enough energy to escape the heat-ofvaporization barrier. Assuming that the molecules in the surface are equilibrated at the surface 

temperature, Ts, it may be shown that the average escape velocity is 

41
( )

( )



 





+

+ =

−

−

e

v

e

2 v

e

2

e

erfc v

2

v e

1 v e

v 2

e 2

2

e 2

β π β β

β

β

β

, (47)

where 

RTs

W

2 β = , (48)

where W is the molecular weight of the escaping molecules and the minimum velocity for escape 

 is e v

W

H

v vap

e

2 = . (49)

The equilibrium vapor pressure is determined by equating the outward and inward fluxes at 

equilibrium (see arguments of section 3.3.1), i.e., 

W

RT

RT

p

s

e s out

π

8

4 Γ = . (50)

To test the predictions of this part of the model apart from imperfections in the predicted value of

the heat of vaporization, I use the experimental value of Hvap in the vapor-pressure formulas. With

no adjustment parameters at all, the standard deviation of the predicted from experimental values 

is 95%. Evidently, there are more serious shortcomings in the vapor-pressure model. I 

experimented with several empirical modifications of the model and got interesting results by 

using an adjustable parameter to scale the value of Hvap in computing the escape velocity of 

equation 49. This strategy results in a standard deviation of about 36%, a much improved 

accuracy but possibly not sufficient for use in the multicomponent evaporation code. The best-fit 

scaling factor has a value of 0.68, i.e., the calculation is significantly improved by assuming that 

only 68% of the full heat of vaporization must be overcome in order to escape the liquid surface. 

These results are illustrated in Figure 22. Noting the slight downward tendency of the error with 

increasing Hvap in the figure, I tried using a two-parameter fit to the fraction of Hvap used to 

compute ve. This improved the standard deviation slightly to 30%. The model may well be further 

improved and placed on a more sound theoretical basis by doing MD studies to help inform the 

assumptions. For example, perhaps decreasing the minimum escape velocity improves the 

idealized model because molecules tend to equilibrate, on average, at a value of potential energy 

somewhat above that in the bulk by means of collisions in the interfacial region closest to the bulk 

liquid. A great advantage of this bootstrapping partnership between discrete and continuum

descriptions is that, by using the same model potential in both, the physics of the evaporation 

42
Heat of Vaporization (kcal/mole)

0123456789 10

% E

rror in V

a

p

or Pres

s

ure a

t B

oilin

g Poin

t

-100

-80

-60

-40

-20

0

20

40

60

80

100

+σ

−σ

Nonpolar Molecules

Polar Molecules

Figure 22. Accuracy of the simple vapor-pressure model using 61 L-J parameter sets for both polar

and nonpolar molecules. 

process can be studied and built into the continuum model apart from the behavior of any real 

substance. Separate studies can then address the issue of the best potential model to use. It goes 

without saying that the continuum models described previously can be implemented with any 

potential-energy function while still preserving the relatively rapid computational qualities 

essential to incorporation into a multicomponent burning-rate model. 

One might argue that empirically based engineering correlations could do considerably better at 

predicting vapor pressures. In fact, I applied the Reid method (44, 45) based on the 

corresponding-states principle to the same database of molecules used previously and found that 

the standard deviation was only 4%. However, these methods are apparently much less successful 

when applied to molecule mixtures. This may be an important limitation for energetic-material 

combustion, where the surface is multicomponent even for the simplest case of ozone as explained 

in section 3.3.2. 

43
44

5.2 MD Simulations of the Condensed Phase 

Another fruitful area where the MD approach has played a role and enjoys bright future prospects 

is in the condensed phase. An example of past use is in determining the idealized mass densities 

(36) of putatively pure polymers of cellulose mono-, di-, and trinitrate. These values are needed in 

the semi-empirical burning-rate model previously discussed for nitrate-ester propellants. In 

addition, the method was also applied (36) to the unreacting solid form of JA2 propellant to check 

the accuracy of the JA2 mass density computation by additive molar volumes, the method 

employed by the CYCLOPS (29, 30) code. Figure 23 shows a particular relaxed configuration of 

JA2, believed to be the first published “anatomically correct” molecular view of a real propellant. 

The CYCLOPS-computed value of the density is 1.56 g/cm3

, which compares very well with the 

MD-computed value of 1.59 ± 0.02 g/cm3

, which, in turn, compares very well with experimental 

value of 1.57 ± 0.01 g/cm3

. 

NC2

NC3

NG

DEGDN

JA2 Propellant Simulation

NC2

NC3

NG

DEGDN

JA2 Propellant Simulation

Figure 23. One view of a relaxed configuration of JA2 propellant computed by a MD

simulation (36) consisting of 2 chains of 15 monomers representing NC, 16

molecules of diethylene glycol dinitrate (DEGDN), and 10 molecules of NG. 

Oxygen atoms are in red, nitrogen in blue, carbon in grey, and hydrogen in white. 

A number of the component molecules can be identified as indicated. This is

believed to be the first computed molecular representation of a real propellant 

formulation. 
Besides the purpose just described, such a MD simulation could be used to determine propellant 

heats of formation that properly account for molecular interactions between dissimilar ingredients. 

Currently, all such interactions are ignored for lack of computational capability and the expectation 

that heats of solution will generally be small, but we previously found (29, 30) that small changes 

in the heats of formation of propellant ingredients can have surprisingly large effects on the 

burning rate at low and intermediate pressures. This effect disappears at maximum gun pressures, 

and therefore will not affect thermodynamic equilibrium calculations, which presently must make 

the same assumption (that heats of mixing are negligible). The reasons for this unexpected 

sensitivity of the burning rate to ingredient heats of formation cannot be neatly isolated because it 

is a highly integrated effect, but it probably arises from a close competition between energy release 

by reactions and the locally dissipative processes of convection, conduction, and molecular 

diffusion. 

Finally, I will express my belief that, while a full MD description of the gas phase will never be 

able to compete with the continuum description in terms of calculational efficiency and accuracy 

(for reasons discussed at the beginning of section 5), neither will a continuum description 

ultimately be able to compete with a discrete molecular description of the condensed-phase 

processes. It may be a long time before this promise is fulfilled, but if a complete description of 

propellant combustion is to be realized, there appears to be little choice but to pursue the MD 

approach tenaciously. In my discussions with molecular dynamicists, I have often observed an 

irrepressible optimism that the method can be made to work for phase changes, for subtle transport 

effects such as thermal diffusion, and even for reactive processes. However, examples given are 

usually for systems that are idealized in the extreme. To be applicable to a practical propellant 

burning-rate model, these simulations will have to be made to work for very general systems, even 

particularly difficult cases involving very large numbers of particles. Perhaps it is best that the 

dynamicists not fully understand the difficulties that await them! 

6. Conclusions 

In this work, I have made a reasoned attempt to predict the most promising course of future 

research aimed at predicting the burning rates of solid propellants from their ingredients. The last 

15 years has seen the rise to dominance of models that treat the gas phase on the level of 

elementary reactions with full multicomponent transport. This explicit recognition of chemical 

specificity has been a necessary precursor to predictive capability. Energetic materials used as 

propellants are designed to produce gas pressure to accomplish work, and any model of burning 

rate must include descriptions of the condensed phase and the gasification mechanism. However, 

attempts to treat the condensed-phase and surface processes at the same level of rigor as the gas 

phase have been stymied by formidable difficulties in both experimental and theoretical 

approaches. It is possible that experiments may never be devised to provide the kind of accuracy 

45
and detail that we have come to expect in the gas phase; at the least, it will probably be a very long 

time in coming. Despite this sobering possibility, we have shown that some degree of 

predictability may be possible using the semi-empirical approach embodied in the CYCLOPS 

code. At present, this predictability appears to be limited to members of well-studied classes of

propellants, such as double-base and RDX/HMX/binder propellants. However, as continuing 

research encompasses ever wider classes of ingredients, generalizations for the semi-empirical 

aspects of the code may well enable wider applicability. For example, it has been suggested that a 

single pyrolysis law for most propellants might prove sufficient if the sensitivity of the burning 

rate to this law is low. 

In my opinion, the greatest hope for treating the condensed-phase and surface processes in a full 

3-phase, first-principles model lies with the developing field of reactive MD. Many obstacles, 

both known and unknown, will have to be overcome in treating condensed-phase reactions by this 

approach. It is not clear that reactive force fields with sufficient generality can be developed. It is 

not clear that methods for treating the many-body interactions can be developed with sufficient 

generality. However, there are ideas in the community about how to proceed; only time and effort

will prove their value. It is worth remembering that burning rate is a highly integrated 

macroscopic consequence of almost unfathomably numerous microscopic processes. Because of 

this, the burning rate has often proved to be extraordinarily insensitive to the underlying processes. 

Thus, one may be justifiably hopeful that even an imperfect description of the detailed processes 

may lead to predictions of macroscopic phenomena that will provide insights and guidance of a 

practical nature in the development of new and optimized propellant formulations. 

46
7. References 

1. Miller, M. S.; Vanderhoff, J. A. Burning Phenomena of Solid Propellants; ARL-TR-2551; 

U.S. Army Research Laboratory: Aberdeen Proving Ground, MD, July 2001. 

2. Adams, G. K. The Chemistry of Solid Propellant Combustion: Nitrate Ester or Double Base 

Systems. In Mechanics and Chemistry of Solid Propellants; Eringen, A.C., Liebowitz, H., 

Koh, S. L., Crowley, J. M., Eds.; Pergamon Press: Elmsford, NY, 1967. 

3. Miller, M. S.; Anderson, W. R. Energetic-Material Combustion Modeling With Elementary 

Gas-Phase Reactions: A Practical Approach. In Solid Propellant Combustion Chemistry, 

Combustion, and Motor Interior Ballistics; Yang, V., Brill, T. B., Ren, W. Z., Eds.; Progress 

in Astronautics and Aeronautics Series; American Institute of Aeronautics and Astronautics: 

New York, 2000; Vol. 185, Chapter 2.12, pp 501–531. 

4. Liau, Y.-C.; Yang, V. Analysis of RDX Monopropellant Combustion With Two-Phase 

Subsurface Reactions. Journal of Propulsion and Power 1995, 11, 729–739. 

5. Davidson, J. E.; Beckstead, M. W. Improvements to Steady-State Combustion Modeling of 

Cyclotrimethylenetrinitramine. Journal of Propulsion and Power 1997, 13, 375–383. 

6. Parr, R. G.; Crawford, B. L., Jr. A Physical Theory of Burning of Double-Base Rocket 

Propellants. I. Journal of Physical and Colloid Chemistry 1950, 54, 929–952. 

7. Rice, O. K.; Ginell, R. The Theory of the Burning of Double-Base Rocket Powders. Journal 

of Physical and Colloid Chemistry 1950, 54, 885–917. 

8. Streng, A. G. Combustion and Explosive Properties of Ozone. Explosivstoffe 1960, 10, 218–

225. 

9. Ben-Reuven, M.; Caveny, L. H. MAE Report 1455; Princeton University: Princeton, NJ, 

January 1980. 

10. Miller, M. S. Three-Phase Combustion Modeling: Frozen Ozone, A Prototype System. In 

Proceedings of the Materials Research Society Symposium: Decomposition, Combustion, 

and Detonation Chemistry of Energetic Materials; Brill, T. B., Russell, T. P., Tao, W. C., 

Wardle, R. B., Eds.; Materials Research Society: Pittsburgh, PA, 1996; pp 169–180. 

11. Sandri, R. On the Decomposition Flame of Liquid Ozone-Oxygen Mixtures in a Tube. 

Combustion and Flame 1958, 2, 348–352. 

12. Miller, M. S. In Search of an Idealized Model of Homogeneous Solid Propellant 

Combustion. Combustion and Flame 1982, 46, 51–73. 

47
13. Anderson, W. R.; Haga, S. W.; Nuzman, J. F.; Kotlar, A. J.; Anderson, R. J. U.S. Army

Research Laboratory, Aberdeen Proving Ground, MD. PREAD Computer Code: A 

Versatile, Portable FORTRAN Computer Code for Interpreting the Complex Chemical 

Kinetics and Transport Properties of a Premixed, Laminar, Steady-State Flame, unpublished 

work, developed 1990–2002. 

14. Anderson, R. J.; Nuzman, J. F.; Anderson, W. R.; Bitely, J. J. U.S. Army Research 

Laboratory, Aberdeen Proving Ground, MD. ChemPlot Computer Code: A Portable JAVA 

Computer Code for Rapid, On-line Visualization of Complex Chemical Kinetic Code 

Outputs, unpublished work, developed 1997–2002. 

15. Anderson, R. J.; Anderson, W. R. U.S. Army Research Laboratory, Aberdeen Proving

Ground, MD. ELEMAP: An Interactive, Portable, JAVA Computer Code for Rapid 

Visualization of Chemical Pathways Diagrams Related to Complex Chemical Kinetic Code 

Outputs, unpublished work, developed 1999–2002. 

16. Hougen, O. A.; Watson, K. M. Chemical Process Principles Charts, 1st ed.; John Wiley & 

Sons: New York, 1946. 

17. Fifer, R. A. Chemistry of Nitrate Ester and Nitramine Propellants. In Fundamentals of

Solid-Propellant Combustion; K. K. Kuo, M. Summerfield, Eds.; Progress in Astronautics 

and Aeronautics Series; American Institute of Aeronautics and Astronautics: Reston, VA, 

1984; Vol. 90. 

18. Melius, C. F. Thermochemical Modeling: II. Application to Ignition and Combustion of 

Energetic Materials. In Chemistry and Physics of Energetic Materials 1990; S. Bulusu, Ed.; 

NATO ASI Series, Kluwer Academic Publishers: Norwell, MA; Vol. 309; pp 51–78. 

19. Prasad, K.; Yetter, R. A.; Smooke, M. D. An Eigenvalue Method for Computing the 

Burning Rates of RDX Propellant. Combustion Science and Technology 1997, 124, 35–82. 

20. Zimmer-Galler, R. AIAA Journal 1968, 6, 2107–2110. 

21. Glaskova, A. P. The Effects of Catalysis on the Combustion of Explosives. Fizika Goreniya 

i Vzryva 1974, 10, 323–334. 

22. Miller, M. S. U.S. Army Research Laboratory, Aberdeen Proving Ground, MD. 

Unpublished data, 1985–1993. 

23. Zenin, A. HMX and RDX: Combustion Mechanism and Influence on Modern Double-Base 

Propellant Combustion. Journal of Propulsion and Power 1995, 11, 752–758. 

24. Ulas, A.; Lu, Y. C.; Kuo, K. K.; Freyman, T. In Proceedings of the 32nd JANNAF 

Combustion Meeting, Huntsville, AL, 23−27 October 1995; Chemical Propulsion 

Information Agency: Laurel, MD, 1995; CPIA Publication 631, Vol. I, pp 461–469. 

48
25. Atwood, A. I.; Boggs, T. L.; Curran, P. O.; Parr, T. P.; Hanson-Parr, D. M. Burning Rate of 

Solid Propellant Ingredients, Part 1: Pressure and Initial Temperature Effects. Journal of 

Propulsion and Power 1999, 15, 740–747. 

26. Homan, B. E.; Miller, M. S.; Vanderhoff, J. A. Absorption Diagnostics and Modeling 

Investigations of RDX Flame Structure. Combustion and Flame 2000, 120, 301–317. 

27. Yetter, R. A.; Dryer, F. L.; Allen, M. T.; Gatto, J. L. Development of Gas-Phase Reaction 

Mechanisms for Nitramine Combustion. Journal of Propulsion and Power 1995, 11, 683–

697. 

28. Kee, R. J.; Grcar, J. F.; Smooke, M. D.; Miller, J. A. A Fortran Program for Modeling 

Steady Laminar One-Dimensional Premixed Flames; SAND85-8240; Sandia National 

Laboratories: Albuquerque, NM, December 1985; reprinted March 1991.

29. Miller, M. S.; Anderson, W. R. CYCLOPS, A Breakthrough Code to Predict SolidPropellant Burning Rates; ARL-TR-2910; U.S. Army Research Laboratory: Aberdeen 

Proving Ground, MD, February 2003. 

30. Miller, M. S.; Anderson, W. R. A Burning-Rate Predictor for Multi-Ingredient Propellants: 

Nitrate-Ester Propellants. Journal of Propulsion and Power, submitted for publication, 2002. 

31. Zenin, A. A. Study of Combustion Mechanism of Nitramine-Polymer Mixtures; Final 

Technical Report; European Research Office of the U.S. Army: London, UK, November 

1998; Contract No. N68171-97-M-5771. 

32. Zenin, A. A. Study of Combustion Mechanism of Nitramine-Polymer Mixtures; Final 

Technical Report, European Research Office of the U.S. Army: London, UK, August 2000; 

Contract No. N68171-99-M-6238. 

33. Zenin, A. A. Study of Combustion Mechanism of New Polymer/Oxidizer Mixtures; Final 

Technical Report, European Research Office of the U.S. Army: London, UK, May 2002; 

Contract No. N68171-01-M-5482. 

34. Leider, H. R.; Seaton, D. L. Nitrate Ester Decomposition and Degradation of Molecular 

Weight in Nitrocellulose From Thermal Decomposition of PBX-9404 Below 100 °C; UCRL52776; Lawrence Livermore National Laboratory: Livermore, CA, May 1979. 

35. Todd, J.; Glasser, W. G. NMR Spectroscopy of Nitrocellulose Samples; In Final Report on 

Cellulose and Cellulose Nitrate Characterization; W. G. Glasser, A. G. Zink, Authors; 

Contract Report to AlliantTechSystems, Inc.: Radford, VA, February 1996; performed at 

Virginia Polytechnic Institute and State University, Blacksburg, VA, Section IV. 

36. Bunte, S. W.; Miller, M. S. Atomistic Simulations of the Physical Properties of Nitrate 

Esters; ARL-TR-2496; U.S. Army Research Laboratory: Aberdeen Proving Ground, MD, 

May 2001. 

49
37. Juhasz, A. A., Ed. Round Robin Results of the Closed Bomb and Strand Burner; CPIA 

Publication 361; Chemical Propulsion Information Agency: Laurel, MD, July 1982. 

38. Atwood, A. I.; Price, C. F.; Curran, P. O.; Zwierzchowshi, N. G. Burning Rate, Radiant 

Ignition, and Global Kinetics of a Nitrocellulose Propellant. In Proceedings of the 25th 

JANNAF Combustion Meeting, Huntsville, AL, 24−28 October 1988; Chemical Propulsion 

Information Agency: Laurel, MD, 1988; CPIA Publication 498, Vol. I, pp 69–81. 

39. Juhasz, A.; Bullock, C.; Homan, B.; Devynck, D. Micro Closed Bomb for Characterizing 

the Burning of Propellants at Gun Pressures. In Proceedings of 36th JANNAF Combustion 

Subcommittee Meeting, Cocoa Beach, FL, 18−21 October 1999; Chemical Propulsion 

Information Agency: Laurel, MD, 1999; CPIA Publication 691, Vol. I, pp 175–187. 

40. Heller, C. A.; Gordon, A. S. Structure of the Gas Phase Combustion Region of a Solid 

Double Base Propellant; Journal of Phys. Chem. 1955, 59, 773. 

41. Lengelle, G.; Bizot, A.; Duterque, J.; Trubert, J. F. Steady-State Burning of Homogeneous 

Propellants. In Fundamentals of Solid Propellant Combustion; Kuo, K., Summerfield, M., 

Eds.; Progress in Astronautics and Aeronautics Series; American Institute of Astronautics 

and Aeronautics: Reston, VA, 1984; Vol. 90. 

42. Vanderhoff, J. A.; Anderson, W. R.; Kotlar, A. J. Dark Zone Modeling of Solid Propellant 

Flames. In Proceedings of the 29th JANNAF Combustion Subcommittee Meeting, Hampton, 

VA, 19−23 October 1992; Chemical Propulsion Information Agency: Laurel, MD, 1992; 

CPIA Publication 593, Vol. II, p 225. 

43. Miller, M. S.; Anderson, W.R. Prediction of Advanced-Nitramine-Propellant Burning Rates 

With the CYCLOPS Code; ARL-MR-552; U.S. Army Research Laboratory: Aberdeen 

Proving Ground, MD, March 2003. 

44. Reid, R. C.; Sherwood, T. K. The Properties of Gases and Liquids, 2nd ed.; McGraw-Hill: 

New York, 1966. 

45. Reid, R. C.; Prausnitz, J. M.; Sherwood, T. K. The Properties of Gases and Liquids, 3rd ed.; 

McGraw-Hill: New York, 1977. 

50
Appendix. Propellant Ingredients and Formulations 

Propellant Ingredients:

ADN ammonium dinitramide 

AMMO poly 3-azidomethyl-3-methyl oxetane 

BAMO poly 3,3-(bis)azidomethyl oxetane 

CBIH copolymer of butadiene and isoprene with hydroxyl terminated groups 

CL20 2,4,6,8,10,12-hexanitrohexaazaisowurtzitane, rectangular crystal size 5 × 3 × 3 μm 

DEGDN diethylene glycol dinitrate 

GAP1U glycidyl azide-polyurethane copolymer 

GAP2 glycidyl azide polymer (molecular mass of 2000) 

HMX cyclotetramethylenetetranitramine 

PU polyurethane 

PUNE mixture of PU and NE 

NC nitrocellulose (cellulose nitrate)

NC1 cellulose mononitrate 

NC2 cellulose dinitrate 

NC3 cellulose trinitrate 

NE mixture of two nitroesters (dinitratdietileneglycole and dinitrattrietileneglicole) 

NG nitroglycerin (glycerin trinitrate) 

RDX cyclotrimethylenetrinitramine 

THF tetrahydrofuran 

Nominal Propellant Formulations Assumed by CYCLOPS Code:

JA2 60% NC (13.1 %N), 25% DEGDN, 15% NG 

M9 59.1% NC (13.25 %N), 40.9% NG 

M10 100% NC (13.15 %N) 

51
INTENTIONALLY LEFT BLANK.

52
2DEFENSE TECHNICAL ABERDEEN PROVING GROUND

 INFORMATION CENTER 

 DTIC OCA 2 DIR USARL 

8725 JOHN J KINGMAN RD AMSRL CI LP (BLDG 305) 

 STE 0944 AMSRL CI OK TP (BLDG 4600) 

FT BELVOIR VA 22060-6218 

 1 COMMANDING GENERAL 

US ARMY MATERIEL CMD 

 AMCRDA TF 

5001 EISENHOWER AVE 

ALEXANDRIA VA 22333-0001 

1 INST FOR ADVNCD TCHNLGY 

THE UNIV OF TEXAS AT AUSTIN 

3925 W BRAKER LN STE 400

AUSTIN TX 78759-5316 

 1 US MILITARY ACADEMY 

MATH SCI CTR EXCELLENCE 

 MADN MATH

 THAYER HALL 

WEST POINT NY 10996-1786

 1 DIRECTOR 

US ARMY RESEARCH LAB 

 AMSRL D 

DR D SMITH 

2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 1 DIRECTOR 

US ARMY RESEARCH LAB 

AMSRL CS IS R 

2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 3 DIRECTOR 

US ARMY RESEARCH LAB 

AMSRL CI OK TL 

2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 3 DIRECTOR 

US ARMY RESEARCH LAB 

AMSRL CS IS T 

2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

53
NO. OF

COPIES ORGANIZATION

ABERDEEN PROVING GROUND

 28 DIR USARL 

AMSRL WM BD

W R ANDERSON 

 R A BEYER 

A L BRANT 

 S W BUNTE 

C F CHABALOWSKI

L M CHANG

T P COFFEE 

 J COLBURN

P J CONROY 

R A FIFER 

B E FORCH 

B E HOMAN 

S L HOWARD

 P J KASTE 

A J KOTLAR 

 C LEVERITT 

 K L MCNESBY

 M MCQUAID

 M S MILLER 

A W MIZIOLEK 

 J B MORRIS 

 J A NEWBERRY 

 M J NUSCA 

 R A PESCE-RODRIGUEZ 

G P REEVES 

B M RICE 

R C SAUSA 

A W WILLIAMS 

54
NO. OF

COPIES ORGANIZATION

55

 1 M MILLER 

 1124 COWPENS AVE 

 TOWSON MD 21286