Abstract
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 [1โ3] 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 [6โ9] upon initiation by a disturbance. The simulation of undesirable nonlinear axial combustion instability symptoms in SRMs employing cylindrical and noncylindrical propellant grains [10โ13] 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 [10โ13], 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.
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]: where 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 (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: where is the equivalent heat input from the gas to the solid phase via 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]: Through substitution, one arrives at which conforms to
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 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] 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] The coefficients and 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] In finite difference format (in this example, for first-order accuracy in time), the temperature at a given internal location may be presented as For a first-order explicit scheme, note that the Fourier stability requirement stipulates that [21] 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 [10โ13]. The corresponding spatial step would typically be set to the Fourier stability requirement limit, namely
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 relative to spatial nodes deeper in the propellant, a given node's temperature may be updated via The boundary condition at the propellant surface (, ) to first-order accuracy may be applied through 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 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 with time when In the unusual case where was greater than , one would have leading , 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 () increment sizing below a certain threshold sizing, the surface boundary condition as relates to net heat input should be stipulated by
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, 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 is the universal gas constant, is the activation energy, and is the Arrhenius coefficient. In this example, the effect of local static pressure 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, is 1675?m/s, and 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]:
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 , such that As a result, for ?s-1, the true surface temperature remains unchanged (constant) relative to the nominal instantaneous Arrhenius value 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 (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 , that would also retain the physics of energy exchange as represented in a general form by (6).
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 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 , defined by is a pertinent parameter that can be charted with respect to driving frequency (when is unity, there is a flat response to the cyclic driving [16, 17]; in this case, 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 , 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 . The response function is defined in terms of mean and fluctuating components of static pressure and incoming mass flow from the propellant surface [15]: 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]
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 in the present model.
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 , rather than a nominal reference value, as discussed later). This presentation format conforms with that presented in [26, Figures 3 and 4].
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 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 (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 . 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.
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 ( 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].
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 was adjusted upward in moving from the lower operating pressure to the higher.
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.
Nomenclature| : | Arrhenius coefficient, m/s |
| : | De St. Robert coefficient, m/s-Pan |
| : | 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) |
| : | Local gas static pressure, Pa |
| : | 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 |