About this Journal Submit a Manuscript Table of Contents
International Journal of Aerospace Engineering
Volume 2008 (2008), Article ID 826070, 10 pages
Research Article

Transient Burning Rate Model for Solid Rocket Motor Internal Ballistic Simulations

Department of Aerospace Engineering, Ryerson University, 350 Victoria St., Toronto, Ontario, Canada M5B 2K3

Received 8 January 2007; Revised 30 April 2007; Accepted 8 August 2007

Academic Editor: R. I. Sujith

Copyright © 2008 David R. Greatrix. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


A general numerical model based on the Zeldovich-Novozhilov solid-phase energy conservation result for unsteady solid-propellant burning is presented in this paper. Unlike past models, the integrated temperature distribution in the solid phase is utilized directly for estimating instantaneous burning rate (rather than the thermal gradient at the burning surface). The burning model is general in the sense that the model may be incorporated for various propellant burning-rate mechanisms. Given the availability of pressure-related experimental data in the open literature, varying static pressure is the principal mechanism of interest in this study. The example predicted results presented in this paper are to a substantial extent consistent with the corresponding experimental firing response data.

1. Introduction

An important aspect in the study of the internal ballistics of solid-propellant rocket motors (SRMs) is the ability to understand the behaviour of a given motor under transient conditions, that is, beyond what would be considered as quasisteady or quasiequilibrium conditions. Transient combustion and flow conditions arise for example during the ignition and chamber filling phase [13] prior to nominal quasisteady operation, during the propellant burnout and chamber emptying phase [4, 5] in the latter portion of a motor's firing, and on occasion when a motor experiences axial or transverse combustion instability symptoms [69] upon initiation by a disturbance. The simulation of undesirable nonlinear axial combustion instability symptoms in SRMs employing cylindrical and noncylindrical propellant grains [1013] has provided the motivation for the present study.

SRM internal ballistic simulation models incorporate algorithms for describing the internal flow and the mass input to the core flow from the burning surface of the solid propellant. More recent models may also incorporate the deflection of the surrounding structure, for example, propellant grain, casing, and heavyweight (e.g., steel) static test sleeve [1013], as reflected in the schematic SRM diagram of Figure 1. For numerical models, in the case of unsteady operation under transient flow conditions, one ideally would capture the dynamic characteristics of both the flow and combustion to a level of accuracy that would enable the prediction of inherent design limits for a given motor (e.g., stable or unstable operation at a given chamber pressure 𝑝𝑐 and initiating pressure disturbance of 𝑝𝑑). An example head-end pressure-time profile from [13] illustrating classical axial combustion instability symptoms, of a limit-magnitude travelling axial shock wave moving back-and-forth within the motor chamber superimposed on a base dc pressure shift (approaching 5 MPa) above the normal operating pressure of approximately 10 MPa, is given in Figure 2. While in some instances the assumption of a quasisteady (i.e., rapid kinetic rate) burning rate response of the propellant to local flow conditions (e.g., static pressure and core mass flux above the burning surface) may be adequate for predictive purposes, this may not be true in other cases where a time or frequency dependence in this combustion response may more directly influence the nonsteady internal ballistics of the motor.

Figure 1: Cylindrical-grain SRM schematic model setup for static test stand firing in a laboratory.
Figure 2: Predicted unsteady head-end pressure-time profile, for SRM undergoing axial combustion instability symptoms (𝐾𝑏=20000 s−1, Δ𝐻𝑠=150000 J/kg, Δ𝑝𝑑=2 atm).

The Zeldovich-Novozhilov (Z-N) phenomenological approach [14, 15], in its most general sense, was considered a good basis for the development of a numerical transient burning rate model that could function in an overall dynamic internal ballistic simulation environment. The Z-N energy conservation criterion in this context requires a numerical heat conduction solution with time in the solid phase (propellant beneath the burning surface), but conveniently, empirical or semiempirical steady-state burning rate information may be used in place of more complex dynamic flame-based reaction rate equations in tying in the gas phase above the propellant surface [15]. This is a distinct advantage for preliminary design and instability evaluation purposes, especially where quicker computational turnaround times are desirable. An approach of this kind can be easily adopted by motor developers, by fitting the model to observed response data for many kinds of propellants through a few parametric constants. In this paper, a working general numerical burning rate model is described. As a new development from previous Z-N models, the present model is general in the sense that the model may be incorporated for various propellant burning-rate driving mechanisms, for example, driven via static pressure, core flow velocity/mass flux, and normal acceleration, or some combination thereof. Past transient burning rate models have typically been constructed in terms of the specific driving mechanism of interest, for example, with equations derived explicitly as a function of pressure. In the present approach, the equations are derived as a function of the quasisteady burning rate, which in turn may be a function of one or several driving mechanisms as noted earlier. Note that the work outlined in this paper is a continuation of the initial burning-rate model development that has been reported in [16, 17] (the reader may find some of the early model developments and evaluations thereof of interest). Example results are presented in this paper in order to provide the reader with some background on the sensitivities of a number of pertinent parameters in the present study. Given the more ready availability of pressure-based burning experimental data in the open literature, comparisons are made to reported experimental pressure-based combustion response data for some composite and homogeneous solid propellants.

2. Numerical Model

The Z-N solid-phase energy conservation result may be presented in the following time-dependent temperature-based relationship [15]: 𝑇𝑖,e=𝑇𝑖1𝑟𝑏𝜕𝜕𝑡0Δ𝑇𝑑𝑥,(1) where 𝑇𝑖,e is the effective initial propellant temperature for instantaneous burning rate estimation, 𝑇𝑖 is the actual initial propellant temperature, and in this context, Δ𝑇=𝑇(𝑥,𝑡)𝑇𝑖 is the temperature distribution in moving from the propellant surface at 𝑥=0 (and 𝑇=𝑇𝑠) to that spatial location in the propellant where the temperature reaches 𝑇𝑖. Here, 𝑟𝑏 is the nominal instantaneous burning rate (later, this parameter is more specifically identified as the unconstrained instantaneous burning rate). For the purposes of the development outlined in this paper, a more direct and general equation is sought relating 𝑟𝑏 and the quasisteady burning rate 𝑟𝑏,𝑞𝑠 (value for burning rate as estimated from steady-state information for a given set of local flow conditions). This was considered potentially more advantageous than exploiting past Z-N correlations that utilized pressure-specific functions such as the pressure-dependent burning rate temperature sensitivity of the propellant (analytic advantage/“short-cut” for transient pressure studies [18]) with mixed success in predicting expected burning rates.

In a finite difference format, energy conservation in the solid phase over a given time increment may be represented by the following equation: 𝑞in𝜌𝑠𝐶𝑠Δ𝑡𝑟𝑏𝑇𝑠𝑇𝑖Δ𝐻𝑠𝐶𝑠=Δ𝑡Δ𝑇Δ𝑥𝑡+Δ𝑡Δ𝑇Δ𝑥𝑡,(2) where 𝑞in is the equivalent heat input from the gas to the solid phase via 𝑞in=𝜌𝑠𝐶𝑠𝑟𝑏,𝑞𝑠𝑇𝑠𝑇𝑖Δ𝐻𝑠𝐶𝑠.(3) One can note the inclusion of the net surface heat of reaction, Δ𝐻𝑠 (sign convention: positive when exothermic output of heat). The quasisteady burning rate may be ascertained as a function of such parameters as static pressure of the local flow, for example through de St. Robert's law [5]: 𝑟𝑏,𝑞𝑠(𝑝)=𝐶𝑝𝑛.(4) Through substitution, one arrives at 𝑟𝑏=𝑟𝑏,𝑞𝑠Δ𝑇Δ𝑥𝑡+Δ𝑡Δ𝑇Δ𝑥𝑡𝑇𝑠𝑇𝑖Δ𝐻𝑠/𝐶𝑠Δ𝑡(5) which conforms to 𝑟𝑏=𝑟𝑏,𝑞𝑠1𝑇𝑠𝑇𝑖Δ𝐻𝑠/𝐶𝑠𝜕𝜕𝑡0Δ𝑇𝑑𝑥.(6)

In (6), 𝑟𝑏 is the nominal (unconstrained) instantaneous burning rate, and its value at a given propellant grain location is solved at each time increment via numerical integration of the temperature distribution through the heat penetration zone of the solid phase. In this respect, the present numerical model differs from past numerical models, for example, as reported by Kooker and Nelson [19]. In those past cases, the thermal gradient at the propellant surface (𝜕𝑇/𝜕𝑥|𝑥=0) was explicitly tied to the true instantaneous burning rate 𝑟𝑏 by a specified function, and in turn, 𝑟𝑏 commonly tied to a variable surface temperature 𝑇𝑠 by an Arrhenius-type expression [18]. In following this established trend, earlier versions of numerical Z-N models did not follow through on using (6) directly, but switched to a burning rate temperature sensitivity correlation such as [18, 20] 𝜕𝑇|||𝜕𝑥𝑥=0=𝑟𝑏𝛼𝑠𝑇𝑠1𝜎𝑝𝑟ln𝑏𝑟𝑏,𝑞𝑠𝑇𝑖,(7) where 𝜎𝑝 is the pressure-based burning rate temperature sensitivity [20]. As reported by Nelson [18], the predicted 𝑟𝑏 augmentation using (7) was commonly lower than expected from experimental observation, for a number of transient burning applications, and relative to various flame-based expressions for thermal surface gradient as a function of 𝑟𝑏, of the general form [18] 𝜕𝑇|||𝜕𝑥𝑥=0=𝑔1𝑟𝑏+𝑔2𝑟𝑏.(8) The coefficients 𝑔1 and 𝑔2 may be found as a function of a number of fixed and variable parameters [18], depending on the given model being employed.

Returning to the present model, in the solid phase, the transient heat conduction is governed by [14, 21] 𝑘𝑠𝜕2𝑇𝜕𝑥2=𝜌𝑠𝐶𝑠𝜕𝑇𝜕𝑡.(9) In finite difference format (in this example, for first-order accuracy in time), the temperature at a given internal location may be presented as 𝑇𝑡+Δ𝑡,𝑥=𝑇𝑡,𝑥+𝛼𝑠Δ𝑡Δ𝑥2𝑇𝑡,𝑥Δ𝑥2𝑇𝑡,𝑥+𝑇𝑡,𝑥+Δ𝑥.(10) For a first-order explicit scheme, note that the Fourier stability requirement stipulates that [21] 𝛼𝑠Δ𝑡Δ𝑥212.(11) While a fourth-order Runge-Kutta scheme was adopted for the present work, the above criterion does prove useful as a guideline. Note that a time step of 1 × 10−7 seconds has been chosen as the reference value for this investigation, based on previous SRM simulation studies [1013]. The corresponding spatial step Δ𝑥 would typically be set to the Fourier stability requirement limit, namely Δ𝑥𝐹𝑜=(2𝛼𝑠Δ𝑡)1/2.

Allowing for the propellant's surface regression of 𝑟𝑏Δ𝑡 with each time step, where 𝑟𝑏 is the true instantaneous rate of regression, and maintaining the surface position as 𝑥=0 relative to spatial nodes deeper in the propellant, a given node's temperature may be updated via 𝑇𝑡,𝑥+𝑟𝑏Δ𝑡=𝑇𝑡,𝑥𝑇𝑡,𝑥Δ𝑥𝑇𝑡,𝑥𝑟Δ𝑥𝑏Δ𝑡.(12) The boundary condition at the propellant surface (𝑥=0, 𝑇=𝑇𝑠) to first-order accuracy may be applied through Δ𝑞e where 𝑇𝑠 represents the net heat input from the gas phase into the regressing solid (this parameter will be discussed in more detail later in the paper).

With respect to the burning surface temperature 𝑇𝑠, one has the option of treating it as constant, or allowing for its variation, depending on the phenomenological approach being taken for estimating the burning rate. In the past, it was not uncommon to encounter estimation models assuming a mean or constant 𝑇𝑠. More recently, usage of an Arrhenius relation for solid pyrolysis dictating a variable 𝑇𝑠 has become more prevalent [14, 18].

The numerical model, as described above and with a preliminary assumption of a constant 𝑟𝑏, is unstable (solution for transient 𝑟𝑏 being strongly divergent), using example propellant properties comparable to Propellant A (see Table 1; nonaluminized ammonium perchlorate/hydroxyl-terminated polybutadiene [AP/HTPB] composite), and over a range of different time and spatial increment sizes [16]. An additional equation limiting the transition of the instantaneous burning rate 𝑑𝑟𝑏𝑑𝑡=𝐾𝑏𝑟𝑏𝑟𝑏.(14) with time is required to physically constrain the model. From the author's previous background in general numerical modeling where lagging a parameter's value is a desired objective, a simple empirical means for applying this constraint is as follows: 𝐾𝑏 The rate limiting coefficient 𝑟𝑏 effectively damps or slows the change in value of the unconstrained burning rate 𝐾𝑏<1Δ𝑡.(15) with time when 𝐾𝑏 In the unusual case where 1/Δ𝑡 was greater than 𝑟𝑏, one would have 𝑟𝑏 leading Δ𝑟𝑏,𝑡+Δ𝑡=𝐾𝑏𝑟Δ𝑡𝑏,𝑞𝑠𝑟𝑏,𝑡+𝐾𝑏𝑇𝑠𝑇𝑖Δ𝑇Δ𝑥𝑡Δ𝑇Δ𝑥𝑡+Δ𝑡.(16), rather than lagging it. In the numerical scheme, an incremental change in burning rate over a given time step would as a result be Δ𝑡 In order to be consistent on input and output heat energy at the propellant surface, such that the converged solution is independent of time (Δ𝑥) and spatial (Δ𝑞e) increment sizing below a certain threshold sizing, the surface boundary condition as relates to net heat input Δ𝑞e=𝐾𝑏𝜌Δ𝑡𝑠𝐶𝑠𝑟𝑏,𝑞𝑠𝑟𝑏𝑇𝑠𝑇𝑖Δ𝐻𝑠𝐶𝑠.(17) should be stipulated by 𝐶𝑠

Table 1: Propellant characteristics.

The empirical coefficient 𝑟𝑏 will need to be set below a maximum permissible value for rendering a nondivergent solution for 𝐾𝑏, and adjusted even further downward in order to match up approximately with combustion response behavior for a typical solid propellant. A limitation of the present approach would of course be the nonavailability of experimental response data, say from T-burner experiments, for the specific propellant in question, to allow for this alignment of 𝑟𝑏 (as opposed to a direct theoretical derivation). An additional or complementary concern would be the strength of the assumption that (14) is in fact a practical means for describing the damped response, to a sufficient degree of accuracy for the purpose at hand, as compared to say, the usage of dynamic flame-based reaction rate equations that also intrinsically limits the movement of the propellant burning rate (e.g., usage of (7) or (8)). Comparing the model's results to a number of different experimental results, for different propellants, would help establish the degree of confidence that might be warranted. Some comparisons are presented later in this paper, that do lend support to the present approach. The structure of composite propellants (as opposed to homogeneous [double-base] propellants) at the local microscale is physically more heterogeneous as a solid mixture of oxidizer crystals and polymer binder than the underlying assumption of solid homogeneity implied by the Z-N model. Using bulk-average propellant properties and adjusting the rate limiting coefficient and net surface heat release to an appropriate setting may nevertheless produce a reasonable and pragmatic predictive capability for motor design and evaluation at the macroscale for both classes of propellants.

In the present approach, the surface thermal gradient is free to find its own value at a given instant, via the numerical scheme for the regressing solid phase. This contrasts with past approaches that dictated an analytical function tying the surface thermal gradient to surface regression rate. One can argue then that the use of (14) or some comparable damping function, while empirical, parallels the approach taken by past researchers in enforcing a stipulated surface thermal gradient (a form of (7) or (8)). Both approaches act to constrain the exchange of energy through the burning surface interface, allow for some variability (through one or more coefficient settings) in better comparing to a given set of experimental data, and prevent so-called burning-rate “runaway” (unstable divergence of 𝑇𝑠 with time [19]).

As part of the model development studies, a variable propellant surface temperature 𝑇𝑠 was given some consideration. In practice, 𝑟𝑏=𝐴𝑠𝑝𝑛𝑠𝐸expas𝑇𝑠,(18) tends to increase with increasing burning rate (although in relative terms, the numerical value change is typically small), as may be expressed as a function of burning rate via the following Arrhenius relation for solid pyrolysis [22]: where 𝐸as is the universal gas constant, 𝐴𝑠 is the activation energy, and 𝑝 is the Arrhenius coefficient. In this example, the effect of local static pressure 𝑛𝑠0 is also included, that is, 𝑛𝑠, as occasionally done in the literature [14, 22]. As demonstrated by the numerically predicted temperature curves within the condensed phase (Propellant A of Table 1, where exponent 𝐴𝑠 is 0, 𝐸as is 1675 m/s, and 𝑇=𝑇𝑖+𝑇𝑠𝑇𝑖𝑟exp𝑏||𝑥||𝛼𝑠.(19) is 100 × 106 J/kg) for differing steady burning rates (Figure 3), the heat penetration zone into the solid propellant becomes substantially thinner as burning rate increases. The profiles in Figure 3 would be expected to conform to the following relationship for a homogeneous solid [14]: 𝐾𝑏

Figure 3: Steady-state temperature distribution in regressing solid propellant.

Analogous to constraining the burning rate to some degree via 𝑟𝑏, one in turn may wish to damp the response of the surface temperature to the change in 𝑑𝑇𝑠𝑑𝑡=𝐾𝑡𝑇𝑠𝑇𝑠.(20), such that 𝐾𝑡=0 As a result, for 𝑇𝑠 s−1, the true surface temperature 𝑇𝑠 remains unchanged (constant) relative to the nominal instantaneous Arrhenius value 𝐾𝑡=1/Δ𝑡 established via (18) for a given burning rate, ranging up to the case when 𝐾𝑡, where there would be no lag between the two values.

As reported in [17], the results for when 𝐾𝑡 is of appreciable value do not coincide with what one would expect from the literature. The appearance of a low-frequency and a high-frequency peak in the frequency response graphs, or an initial peak dropping down to a prolonged plateau in magnitude at higher values for 𝑀, is not commonly observed experimentally (see Figure 4 for an example predicted result for limit magnitude 𝐾𝑡, defined by (21), as a function of driving frequency, at differing values for 𝑇𝑠). The more variable surface temperature in essence (within the predictive model) acts to suppress the degree of augmentation of the principal response peak (as it moves the peak to a higher frequency) for a given cyclic driving mechanism, and introduces a secondary resonant response peak at a low frequency. From a modeling standpoint, one can appreciate that a moving 𝐾𝑡 would tend to counteract the incoming/outgoing energy (to/from the solid) and affect the numerically predicted near-surface temperature and energy transfer. This undoubtedly is playing a role in producing the odd results noted in [17] as 𝑇𝑠 is increased. Given the better overall comparisons to experimental profiles, to date, with an assumed constant 𝐾𝑡=0 (i.e., 𝑇𝑠 s−1) in the predictive, phenomenological model, the numerical results reported in this paper will be generated in a similar fashion. In doing so, it is understood that this does not rule out the possibility of a future model development that would allow for a variable 𝑟𝑏,𝑜=0.01, that would also retain the physics of energy exchange as represented in a general form by (6).

Figure 4: Frequency response of Propellant A, 𝐾𝑏=170000 m/s, Δ𝐻𝑠=0.0 s−1, 𝑟𝑏,𝑞𝑠 J/kg.

3. Results and Discussion

Examining propellants with characteristics as provided in Table 1, one means for comparing and potentially aligning the numerical model to actual firing data is to examine the frequency response to cyclic ±0.001 input, at different values for the burn rate limiting coefficient. For example, application of a 𝑟𝑏,𝑜 m/s sinusoidal cycle on a reference base burning rate 𝐾𝑏 of 0.01 m/s at different driving frequencies on the reference propellant, with a given value for 𝑀, produces a set of limit-amplitude cycle results. The nondimensional limit magnitude 𝑀=𝑟𝑏,peak𝑟𝑏,𝑜𝑟𝑏,𝑞𝑠,peak𝑟𝑏,𝑜,(21), defined by 𝑀 is a pertinent parameter that can be charted with respect to driving frequency (when 𝑟𝑏,peak is unity, there is a flat response to the cyclic driving [16, 17]; in this case, 𝑟𝑏,𝑞𝑠,peak would equal 𝐾𝑏 (0.011 m/s)). See Figure 5 for an example predicted result for Propellant A, for differing values of Ω. The limit magnitude profile may also be described in terms of the dimensionless frequency 𝛼Ω=2𝜋𝑠𝑟2𝑏,𝑜𝑓.(22), defined by [15], 𝑝 To allow for further comparison with frequency response data available in the literature, the simulation results can be presented in a form that relates to a specific flow parameter, in this case static pressure 𝑅𝑝. The real part of the pressure-based response function Ω is typically presented as a function of 𝑅𝑝=𝑚/𝑚𝑝/𝑝.(23). The response function is defined in terms of mean and fluctuating components of static pressure and incoming mass flow from the propellant surface [15]: 𝑅Re𝑝=𝑟𝑏,peak𝑟𝑏,𝑜𝑟1𝑏,𝑞𝑠,peak𝑟𝑏,𝑜1/𝑛11.(24) In the context of this study's application, the real part of the pressure-based response function in this example application can be estimated as [16] 𝑟𝑏,𝑜=0.01

Figure 5: Frequency response of Propellant A, Δ𝐻𝑠=0.0 m/s, 𝑟𝑏,𝑜 J/kg.

Referring to Figure 6, one can compare the pressure-based response for an AP/HTPB propellant (B in Table 1; when not available, some propellant characteristics listed in Table 1 may be approximated from the literature). The base burn rate 𝐴 of this propellant at the test pressure is assumed to be 1.13 cm/s, from available information. While there is some degree of variability in the experimental data [23] commonly encountered with T-burner results [24], in this particular case there is an appreciable level of agreement between the predicted curve of Figure 6 and the test data points. In [23], the authors found that setting the coefficient 𝐵 to 8.0 and coefficient 𝐾𝑏 to 0.69 in the Denison-Baum A-B pressure-coupled response model [22] produced a relatively good fit to the available test data points (the Denison-Baum predicted curve is quite similar to the present model's curve). In Cohen's study of the application of the A-B response model to homogeneous propellants [25], he notes that an increase in value for coefficient A results in a higher peak response magnitude, at a higher dimensionless frequency. This corresponds to the effect of an increase in the value of 𝐵 in the present model, as reflected by Figure 5. Similarly, he notes that a decrease in the value of coefficient 𝐵 results in a higher peak response magnitude, with a relatively small change in dimensionless frequency (at lower values for Δ𝐻𝑠). This corresponds to the effect of an increase in the positive value of 𝐾𝑏=29000 in the present model.

Figure 6: Frequency response of AP/HTPB Propellant B, theory (Δ𝐻𝑠=75000 s−1, 𝑟𝑏,𝑜 J/kg) and experiment [23], in terms of real part of pressure-based function (theory & experiment, pressure at 1200 psig).

Experimental combustion response data for an AP/PBAN (polybutadiene acrylic-acid/acrylonitrile binder) composite propellant designated as A-13 in the literature [22, 24, 26, 27] is also available for comparison. As displayed in [24] (Figure 3 of that paper), test data points produced from various institutions' experiments for that propellant range considerably for any given frequency. For Figure 7 of the present paper, two institutions' data points are chosen (they more or less reflect upper and lower bounds for the collected data; UTC refers to United Technologies Corp., and CIT refers to California Institute of Technology) for comparison to the present model's predicted results. Given the range of the experimental data, it was decided to also produce example upper and lower predictive model curves, rather than a single median curve. The assumed characteristics of this propellant are listed under C in Table 1. The base burn rate Ω of this propellant at the test pressure is assumed to be 0.5 cm/s, an estimate from available information. The log-linear format of the graph corresponds to that used in [24]. Given the range of variability of the experimental data, the agreement with the predictive model in terms of qualitative trends and quantitative magnitudes for some portions of the two predictive curves is encouraging. Figure 8 provides an alternative linear-linear format for the above results, where the pressure-based response is in terms of the dimensionless frequency Ω (output 𝛼𝑠 based on the propellant's actual thermal diffusivity 𝐾𝑏=6700, rather than a nominal reference value, as discussed later). This presentation format conforms with that presented in [26, Figures 3 and 4].

Figure 7: Frequency response of AP/PBAN Propellant C, theory (Δ𝐻𝑠=115000 s−1, 𝐾𝑏=18000 J/kg; Δ𝐻𝑠=68000 s−1, 𝐾𝑏=6700 J/kg) and experiment [24], in terms of real part of pressure-based function (theory & experiment, pressure at 300 psig).
Figure 8: Frequency response of AP/PBAN Propellant C, theory (Δ𝐻𝑠=115000 s−1, 𝐾𝑏=18000 J/kg; Δ𝐻𝑠=68000 s−1, Ω J/kg) and experiment [24, 26], in terms of real part of pressure-based function (theory & experiment, pressure at 300 psig).

In [22, Figure 4], the experimental (T-burner) pressure-based response for the A-13 propellant is displayed for differing pressures (and therefore, differing base burn rates). For the present model, the base burning rates were estimated from available information, as follows: 0.32 cm/s at 100 psig (0.79 MPa abs.), 0.565 cm/s at 400 psig (2.86 MPa abs.), and 0.766 cm/s at 800 psig (5.62 MPa abs.). As noted in [28], and as done as well in [29] for propellants with differing properties, it has been commonplace in the past to present output data in terms of 𝛼𝑠,ref as a function of a reference thermal diffusivity value Ω, namely 3 × 10−4 in2/s or 1.9355 × 10−7m2/s. This reference value differs substantially with the assumed actual value for the A-13 propellant (C in Table 1, i.e., 8.28 × 10−8 m2/s). While it was not indicated explicitly in [22] which value for thermal diffusivity was used for the graph's output Ω=𝑓(𝛼𝑠,ref) (see (22)), a better comparison between the present model's prediction and the test data is attained when, for output only, the reference value for thermal diffusivity noted above is used (of course, for the actual numerical calculations in the solid phase that produce the displayed response results, the actual thermal diffusivity would be used, not the reference value). Conforming to the log-linear format of [22], Figure 9 presents the experimental data and the corresponding predicted curves for the three pressures using 𝐾𝑏=18000. Accounting for the inherent variability of T-burner data, one can see some positive correlation between the predictive curves and the experimental data points. Unfortunately, on the experimental side, no data was collected in these cases at lower frequencies, so the comparisons here for Figure 9 must be acknowledged as being incomplete. However, for the same propellant at 300 psig (2.17 MPa abs.), experimental response levels as shown in Figure 8 can in fact be demonstrated to be quite high (approaching a value of 6 at a lower frequency). A similar propellant to A-13, designated A-35, at 200 psig, shows one experimental data point exceeding a value of 5 at a lower frequency [22]. The predicted curves of Figure 9 clearly show the effect of a lower-base burning rate on augmenting the propellant's response, with the experimental data consistent with this trend at the higher tested frequencies.

Figure 9: Frequency response of AP/PBAN Propellant C, theory (Δ𝐻𝑠=68000 s−1, 𝜇 J/kg) and experiment [22], in terms of real part of pressure-based function, at differing pressures.

In [24, Figure 5], the authors report experimental response data for AP/HTPB propellants, displaying profiles that are less typically observed in the literature. The response at lower frequencies rises as one goes up in frequency, and then tends to plateau as frequency further increases, and in some cases, starts dropping off to some degree at even higher frequencies. In Figure 10, one has experimental data points shown for one such propellant (84% concentration of 5-Δ𝐻𝑠m-diameter AP crystals). The corresponding predicted curve is based on the assumed or estimated propellant characteristics listed under Propellant D in Table 1. A base burn rate of 0.5 cm/s was assumed from the available information. In order to obtain the reasonably good agreement with the experimental data points, a strong endothermic surface heat of reaction is being applied (𝐾𝑏=140000 at −1 × 106 J/kg). This contrasts with the smaller exothermic values applied in the earlier cases. Logically, one would expect an endothermic process to produce such a deadening in the propellant's response behavior. Note further that the presence of reactive or nonreactive particles in the flow can in practice through particle damping (i.e., primarily aerodynamic drag losses reducing the limit-amplitude travelling pressure wave strength) produce reductions in a propellant's apparent combustion response as measured by a T-burner. This possibly may produce in appearance what would be considered an inordinately (in chemical kinetic terms) strong endothermic process [24, 28].

Figure 10: Frequency response of AP/HTPB Propellant D, theory (Δ𝐻𝑠=1000000 s−1, 𝑟𝑏 J/kg) and experiment [24], in terms of real part of pressure-based function (theory and experiment, pressure at 200 psig).

In [30, Figure 3], the authors report experimental pressure response data for a homogeneous (double-base) propellant, designated as JPN in the literature (comprised of approximately 52% nitrocellulose, 43% nitroglycerine, plus additional additives). In Figure 11, one has experimental data points for two operating conditions (800 psig (5.62 MPa abs.), base 𝑟𝑏 of 1.3 cm/s; 1600 psig (11.1 MPa abs.), base 𝐾𝑏 of 2.1 cm/s). The corresponding predicted curves are based on the assumed or estimated propellant characteristics listed under Propellant E in Table 1 and a value for Δ𝐻𝑠 of 39 000 s−1. In order to provide a better comparison between theory and experiment, the value for 𝐾𝑏=39000 was adjusted upward in moving from the lower operating pressure to the higher.

Figure 11: Frequency response of double-base Propellant E, theory (Δ𝐻𝑠=45000 s−1, Δ𝐻𝑠=90000 J/kg [800 psig], 𝑅𝑝 J/kg [1600 psig]) and experiment [30], in terms of real part of pressure-based function.

4. Concluding Remarks

A general numerical model exploiting the Z-N energy conservation approach has been presented. While past numerical models for transient burning rate have tended to use a surface thermal gradient boundary condition for problem resolution, the present model employs the integrated temperature distribution in the solid propellant directly for instantaneous regression rate calculations. The introduction of a burn rate limiting coefficient was necessitated initially by numerical model stability considerations, but in turn this, in conjunction with adjusting the net surface heat of reaction value, allows one to potentially line up the model's response behavior to that observed experimentally for a given solid propellant. The example results presented here (within a certain range of burn rate limiting coefficient and net surface heat release values) are to a substantial extent consistent with corresponding experimental firing response data. They clearly confirm the effect of lower-base burning rate in augmenting the propellant's response to a given driving mechanism.

While static pressure was chosen as the local driving mechanism for example unsteady burning results, the model could just as easily be set up for other flow mechanisms such as core mass flux and radiation. Other mechanisms, such as local normal acceleration resulting from motor spinning or structural vibration, could be modeled in this general scheme. One may need to establish whether the rate limiting coefficient selected through observation of one set of experimental tests (say for 𝐾𝑏) would in turn be consistent for other mechanisms. In addition, it is acknowledged that for some cases, one might encounter some variability in Δ𝐻𝑠 and 𝐴𝑠 in moving from one operational condition (e.g., a given chamber pressure and base burning rate) to another for a given propellant. Having said this, the relative simplicity of the present numerical burning rate model makes it a potentially useful candidate for transient internal ballistic studies of solid-propellant rocket motors, for example, for studies requiring shorter simulation computational turnaround times for the prediction of undesirable axial combustion instability symptoms in a given motor design.

𝐶:Arrhenius coefficient, m/s
𝐶𝑠:De St. Robert coefficient, m/s-Pan
𝐸as:Specific heat (solid phase), J/kg-K
𝑓Activation energy, J/kg
𝑓𝑟:Frequency, Hz
Δ𝐻𝑠:Resonant frequency, Hz
𝐾𝑏:Net surface heat of reaction, J/kg
𝐾𝑡:Burn rate limiting coefficient, s−1
𝑘𝑠:Surface temperature damping coefficient, s−1
𝑀:Thermal conductivity (solid phase), W/m-K
𝑚:Limit magnitude (cyclic input)
𝑛:Mass flow, kg/s
𝑝:Exponent (de St. Robert's law)
𝑞in:Local gas static pressure, Pa
Δ𝑞e::Equivalent heat input, W/m2
Net heat input, W/m2
𝑅𝑝:Universal gas constant, J/kg-K
𝑟𝑏:Pressure-based response function
𝑟𝑏,𝑜:Instantaneous burning rate, m/s
𝑟𝑏,𝑞𝑠::Reference burning rate, m/s
𝑟𝑏Quasisteady burning rate, m/s
𝑇𝑖:Unconstrained burning rate, m/s
𝑇𝑠:Initial temperature (solid phase), K
Δ𝑡:Burning surface temperature, K
Δ𝑥:Time increment, s
Δ𝑥𝐹𝑜:Spatial increment, m
𝛼𝑠:Fourier limit spatial increment, m
𝜌𝑠:Thermal diffusivity (solid phase), m2/s
Ω:Density (solid phase), kg/m3
𝐿:Dimensionless frequency


  1. E. W. Price, “L instability,” in Nonsteady Burning and Combustion Stability of Solid Propellants, L. De Luca, E. W. Price, and M. Summerfield, Eds., vol. 143 of Progress in Astronautics and Aeronautics, pp. 325–362, American Institute of Aeronautics and Astronautics, Washington, DC, USA, 1992.
  2. Y. M. Timnat, Advanced Chemical Rocket Propulsion, Academic Press, New York, NY, USA, 1987.
  3. D. R. Greatrix, J. J. Gottlieb, and T. Constantinou, “Numerical model for pellet-dispersion igniter systems,” Journal of Propulsion and Power, vol. 4, no. 5, pp. 412–420, 1988.
  4. H. Krier, J. S. T'ien, W. A. Sirignano, and M. Summerfield, “Nonsteady burning phenomena of solid propellants: theory and experiments,” AIAA Journal, vol. 6, no. 2, pp. 278–285, 1968.
  5. G. P. Sutton and O. Biblarz, Rocket Propulsion Elements, John Wiley & Sons, New York, NY, USA, 7th edition, 2001.
  6. A. L. Karnesky and S. E. Colucci, “Recent occurrences of combustion instability in solid rocket motors—an overview,” Journal of Spacecraft and Rockets, vol. 12, no. 1, pp. 33–38, 1975.
  7. M. Barrère, “Introduction to nonsteady burning and combustion,” in Nonsteady Burning and Combustion Stability of Solid Propellants, L. De Luca, E. W. Price, and M. Summerfield, Eds., vol. 143 of Progress in Astronautics and Aeronautics, pp. 17–58, American Institute of Aeronautics and Astronautics, Washington, DC, USA, 1992.
  8. L. Strand, “Laboratory test methods for combustion-stability properties of solid propellants,” in Nonsteady Burning and Combustion Stability of Solid Propellants, L. De Luca, E. W. Price, and M. Summerfield, Eds., vol. 143 of Progress in Astronautics and Aeronautics, pp. 689–718, American Institute of Aeronautics and Astronautics, Washington, DC, USA, 1992.
  9. F. S. Blomshield, “Pulsed motor firings,” in Solid Propellant Chemistry, Combustion, and Motor Interior Ballistics, V. Yang, T. B. Brill, and W.-Z. Ren, Eds., vol. 185 of Progress in Astronautics and Aeronautics, pp. 921–958, American Institute of Aeronautics and Astronautics, Washington, DC, USA, 2000.
  10. S. Loncaric, D. R. Greatrix, and Z. Fawaz, “Star-grain rocket motor—nonsteady internal ballistics,” Aerospace Science and Technology, vol. 8, no. 1, pp. 47–55, 2004. View at Publisher · View at Google Scholar
  11. D. R. Greatrix, “Nonsteady interior ballistics of cylindrical-grain solid rocket motors,” in Proceedings of the 2nd International Conference on Computational Ballistics, pp. 281–289, WIT Press, Cordoba, Spain, May 2005.
  12. J. Montesano, D. R. Greatrix, K. Behdinan, and Z. Fawaz, “Structural oscillation considerations for solid rocket internal ballistics modeling,” in Proceedings of the 41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Tucson, Ariz, USA, July 2005, AIAA 2005-4161.
  13. D. R. Greatrix, “Predicted nonsteady internal ballistics of cylindrical-grain motor,” in Proceedings of the 42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Sacramento, Calif, USA, July 2006, AIAA 2006-4427.
  14. B. V. Novozhilov, “Theory of nonsteady burning and combustion instability of solid propellants by the Zeldovich-Novozhilov method,” in Nonsteady Burning and Combustion Stability of Solid Propellants, L. De Luca, E. W. Price, and M. Summerfield, Eds., vol. 143 of Progress in Astronautics and Aeronautics, pp. 601–642, American Institute of Aeronautics and Astronautics, Washington, DC, USA, 1992.
  15. M. Q. Brewster, “Solid propellant combustion response: quasi-steady (QSHOD) theory development and validation,” in Solid Propellant Chemistry, Combustion, and Motor Interior Ballistics, V. Yang, T. B. Brill, and W.-Z. Ren, Eds., vol. 185 of Progress in Astronautics and Aeronautics, pp. 607–638, American Institute of Aeronautics and Astronautics, Washington, DC, USA, 2000.
  16. D. R. Greatrix, “Numerical transient burning rate model for solid rocket motor simulations,” in Proceedings of the 40th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Fort Lauderdale, Fla, USA, July 2004, AIAA 2004-3726.
  17. D. R. Greatrix, “Numerical modeling considerations for transient burning of solid propellants,” in Proceedings of the 41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Tucson, Ariz, USA, July 2005, AIAA 2005-4005.
  18. C. W. Nelson, “Response of three types of transient combustion models to gun pressurization,” Combustion and Flame, vol. 32, pp. 317–319, 1978. View at Publisher · View at Google Scholar
  19. D. E. Kooker and C. W. Nelson, “Numerical solution of solid propellant transient combustion,” Journal of Heat Transfer, vol. 101, no. 2, pp. 359–364, 1979.
  20. K. K. Kuo, S. Kumar, and B. Zhang, “Transient burning characteristics of JA2 propellant using experimentally determined Zel'dovich map,” in Proceedings of the 39th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Huntsville, Ala, USA, July 2003, AIAA 2003-5270.
  21. F. Kreith and W. Z. Black, Basic Heat Transfer, Harper & Row, New York, NY, USA, 1980.
  22. F. E. C. Culick, “A review of calculations for unsteady burning of a solid propellant,” AIAA Journal, vol. 6, no. 12, pp. 2241–2255, 1968.
  23. J. D. Baum, J. N. Levine, J. S. B. Chew, and R. L. Lovine, “Pulsed instability in rocket motors: a comparison between predictions and experiments,” in Proceedings of the 22nd AIAA Aerospace Sciences Meeting, Reno, Nev, USA, January 1984, AIAA 84-0290.
  24. M. W. Beckstead, K. V. Meredith, and F. S. Blomshield, “Examples of unsteady combustion in non-metalized propellants,” in Proceedings of the 36th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Huntsville, Ala, USA, July 2000, AIAA 2000-3696.
  25. N. S. Cohen, “Combustion response functions of homogeneous propellants,” in Proceedings of the 21st AIAA/ASME/SAE/ASEE Joint Propulsion Conference, Monterey, Calif, USA, July 1985, AIAA 85-1114.
  26. M. Shusser, F. E. C. Culick, and N. S. Cohen, “Combustion response of ammonium perchlorate composite propellants,” Journal of Propulsion and Power, vol. 18, no. 5, pp. 1093–1100, 2002.
  27. M. D. Horton, “Use of the one-dimensional T-burner to study oscillatory combustion,” AIAA Journal, vol. 2, no. 6, pp. 1112–1118, 1964.
  28. F. S. Blomshield, R. A. Stalnaker, and M. W. Beckstead, “Combustion instability additive investigation,” in Proceedings of the 35th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Los Angeles, Calif, USA, June 1999, AIAA 99-2226.
  29. R. S. Brown and R. J. Muzzy, “Linear and nonlinear pressure coupled combustion instability of solid propellants,” AIAA Journal, vol. 8, no. 8, pp. 1492–1500, 1970.
  30. R. W. Hart, R. A. Farrell, and R. H. Cantrell, “Theoretical study of a solid propellant having a heterogeneous surface reaction I—acoustic response, low and intermediate frequencies,” Combustion and Flame, vol. 10, pp. 367–380, 1966. View at Publisher · View at Google Scholar