Scale Effects on Solid Rocket Combustion Instability Behaviour

The ability to understand and predict the expected internal behaviour of a given solid-propellant rocket motor under transient conditions is important. Research towards predicting and quantifying undesirable transient axial combustion instability symptoms necessitates a comprehensive numerical model for internal ballistic simulation under dynamic flow and combustion conditions. A numerical model incorporating pertinent elements, such as a representative transient, frequency-dependent combustion response to pressure wave activity above the burning propellant surface, is applied to the investigation of scale effects (motor size, i.e., grain length and internal port diameter) on influencing instability-related behaviour in a cylindrical-grain motor. The results of this investigation reveal that the motor’s size has a significant influence on transient pressure wave magnitude and structure, and on the appearance and magnitude of an associated base pressure rise.


Nomenclature:
A = local core cross-sectional area (m 2 ) a  = longitudinal (or lateral) acceleration (m/s 2 ) D = drag force of gas on particle (N) d = local core hydraulic diameter (m) d p,i = initial grain port diameter (m) E = local total specific energy of gas in core (J/kg) f = frequency (Hz), or Darcy-Weisbach friction factor f r = resonant combustion response frequency (Hz) f 1L = fundamental resonant axial acoustic chamber wave frequency (Hz) h = convective heat transfer coefficient (W/m 2 •K)

OPEN ACCESS
H s = net surface heat release (J/kg) p  = propellant grain length (m) M  = limit magnitude, cyclic input m p = mass of particle (kg) p = local gas static pressure (Pa) p d = magnitude of initial pressure disturbance (Pa) Q = heat transfer from gas to particle (W) r b = instantaneous burning rate (m/s) r b,o = reference burning rate (m/s) r o = base burning rate (m/s) T f = flame temperature, gas phase (K) T i = initial temperature, solid phase (K) T s = burning surface temperature (K) t = time (s) u = core axial gas velocity (m/s) u p = axial particle velocity (m/s) v f = nominal flamefront velocity (m/s) x = axial distance from head end (m) y = distance from propellant surface into the solid phase (m)  = propellant surface roughness (m) κ = vibration-based wall dilatation term, 1/A  A/t (s −1 ) ρ = gas density (kg/m 3 ) ρ p = particle density in flow (kg/m 3 )

Introduction
With respect to the internal ballistics of solid-propellant rocket motors (SRMs), one ideally should be able to understand the behaviour of a given motor under transient conditions, i.e., beyond what would be considered quasi-steady or quasi-equilibrium conditions.Over the last number of decades, there has been a number of research efforts directed towards a better understanding of the physical mechanisms, or at least the surrounding factors, behind the appearance of transient symptoms associated with nonlinear axial combustion instability in SRMs; see [1,2].The motivation for these studies was and is of course to bring this better understanding to bear in more precisely suppressing, if not eliminating, these symptoms.The primary symptom for traditional nonlinear axial combustion instability is a sustained cyclic axial pressure wave presence of significant amplitude in the combustion chamber; a second symptom is the occasional appearance of a base pressure rise in the combustion chamber that accompanies the axial or transverse pressure wave motion that is present, referred to as a dc shift.A comprehensive numerical model for simulation of nonlinear dynamic flow and combustion conditions is essential in the quest for the ability to predict and quantify axial combustion instability symptoms in SRMs.An effective model combines the effects of the unsteady flow, the transient combustion process, and the structural dynamics of the surrounding propellant/casing structure.On the numerical prediction side, as various component models evolve, say for transient burning rate or structural vibration, their incorporation into an overall internal ballistics simulation program allows for new motor firing simulations to take place, which in turn allows for updated comparisons to available experimental firing data, and a better understanding of the influence of various factors.Examples of pulsed motor firing simulations completed previously, at interim stages of earlier simulation model development, may be found in [3][4][5].
In the present investigation, an updated numerical model incorporating the above attributes is used in the prediction of the unsteady instability-related behaviour in cylindrical-grain motors of differing properties and scales.A more recent version of a transient burning rate model ( [6]; based on a general Zeldovich-Novozhilov [Z-N] approach) is employed for the present study.Additionally, the combustion module for the simulation program employs an updated erosive burning model ( [7]; allows for prediction of negative erosive burning at low flow speeds, in addition to positive augmentation at higher flow speeds).Example results are presented in this paper in order to provide the reader with some background on the sensitivities of motor size, as relates to cylindrical-grain SRMs.

Numerical Model
In the present investigation, the SRM is assumed to be positioned on a static test stand, as schematically represented in Figure 1.In this case, the cylindrical-grain motor is free to vibrate radially without any external constraint (i.e., only constrained in some cases by the thick steel static-test sleeve surrounding the aluminum flightweight motor casing), while axial motion is constrained to a large degree by the thrust-measuring load cell (represented in the predictive model as a spring/damper system) at the lefthand boundary of the motor's head end.Under normal (nominal) quasi-equilibrium operating conditions, the internal gas flow (or gas-particle flow, if two-phase) moves smoothly from the burning propellant surface into the axial flow, from left-to-right as per Figure 1 through and beyond the exhaust nozzle.The equations of motion describing the nonsteady core flow within the SRM must be solved in conjunction with the local surface regression rate r b of the solid propellant, and the surrounding structure's instantaneous geometric deformation.As pertains to the present study of a motor having a larger length-to-diameter ratio, the quasi-one-dimensional hydrodynamic conservation equations for the axial gas flow are given below: (2) Energies 2011, 4 Here, the total specific energy of the gas is defined for an ideal gas as E = p/[(−1)ρ] + u 2 /2.The viscous interaction between the gas and particle phases is represented by the drag force D, and the heat transfer from the core flow to the particles is defined by Q (in this format, assuming one inert [unchanging] particle size; [8]).Here, longitudinal acceleration a  appears in the gas momentum and energy equations as a body force contribution within a fixed Eulerian reference (fixing of x = 0 to motor head end, x positive moving right on structure; acceleration of local surrounding structure rightward is designated positive a  ), and may vary both spatially along the length of the motor and with time.The effects of such factors as turbulence can be included through one or more additional equations that employ the information from the bulk flow properties arising from the solution of the above one-dimensional equations of motion.The principal differential equations themselves can be solved via a higher-order random-choice method (RCM) approach [3][4][5]9].Given that the current investigation will not be incorporating particles within the flow (particle loading mass fraction a p ≈ 0%), for brevity, the particle-phase differential equations of motion are not presented here; the reader is referred to [8].
Structural vibration can play a significant role in nonsteady SRM internal ballistic behavior, as evidenced by observed changes in combustion instability symptoms as allied to changes in the structure surrounding the internal flow (e.g., propellant grain configuration, wall thickness, material properties).The level of sophistication required for modeling the motor structure (propellant, casing, static-test sleeve, nozzle) and applicable boundary conditions (load cell on static test stand) can vary, depending on the particular application and motor design.Montesano et al [4] employed a finiteelement approach towards the structural modeling of the given motor configuration.In the present study, a cylindrical-grain configuration allows for a simpler approach via thick-wall theory, as reported in [3].The radial deformation dynamics of the propellant/casing/sleeve are modelled by a series of independent ring elements along the length of the motor.With that representation, the following ordinary differential equation applies for the movement of a radial position r within the assembly: where r o is the reference radius under no chamber pressure, and  is the local instantaneous radial displacement under static loading.Equation 4 may be solved using a backward difference scheme.
Allowing for viscoelastic behavior,  R is the radial damping ratio.The fundamental natural radial frequency,  nR , of the sleeve/casing/propellant assembly may be determined using a commercial finite element analysis software program, if the theoretical approximation is not considered adequate.Axial motion along the length of the structure is modelled via beam theory, and bounded by the spring/damper load cell at the motor's head end.Viscous damping is applied in both the radial and axial directions.Reference structural properties are assumed for an ammonium-perchlorate/hydroxylterminated polybutadiene (AP/HTPB) composite propellant surrounded by an aluminum casing and steel sleeve.
With respect to transient, frequency-dependent burning rate modeling that allows one to solve for the local r b value for the burning propellant surface at a given time, a generalized Z-N solid-phase energy conservation approach is used in the present simulation program.It may be represented by the following time-dependent temperature-based relationship [6]: where r b,qs is the quasi-steady burning rate (value for burning rate as estimated from steady-state information for a given set of local flow conditions above the propellant surface), T i is the initial propellant temperature (e.g., room or outside air temperature), and in this context, ∆T = T(y,t) − T i is the temperature distribution in moving from the burning propellant surface at y = 0 (and T = T s ) to that spatial location in the propellant where the temperature reaches T i .One may note the inclusion of a net surface heat release term, ΔH s , in the calculations, to allow for a better comparison to transient-response experimental data.This term is commonly included as an add-on correction to the solid-phase energy contribution in an energy-conservation model involving solid propellants (e.g., see Equation 12below for erosive burning).On a pragmatic level, adjusting the value of H s positively upward (more exothermic contribution) increases the burning response magnitude of the propellant at a given frequency, and visa versa for a more endothermic contribution, which acts to deaden the response [6].For the numerical solution of the integral of Equation 5, the transient heat conduction in the solid phase can be solved by an appropriate finite-difference scheme.One needs to take care in setting the solid-phase spatial increment y, to be in accordance with the Fourier stability limit y Fo , which is a function of the chosen time increment t [6].The time increment itself must be coordinated between the flow and structural model solution systems [6].In Equation 5, r b * is the nominal (unconstrained) instantaneous burning rate.The actual instantaneous burning rate r b may be found as a function of r b * through the rate limiting equation provided below [6], The rate limiting coefficient K b effectively damps ("slows down") the unconstrained burning rate r b * to allow for a better comparison to experimental transient response data, and to prevent artificial burning-rate divergence ("runaway"), when for a finite time increment t: A higher value for K b , moving closer to but still less than 1/t, tends to shift the peak combustion response frequency f r to a higher value.Adjusting the value for K b and H s permits one to better compare to a given solid propellant's experimentally observed transient response behaviour, for example as determined from a T-burner [1].An example result, with a comparison to corresponding experimental data for the pressure-based combustion response of a conventional composite solid propellant, is provided in Figure 2. The value for f r is around 1150 Hz in that example.H s = +75,000 J/kg) and experiment [10], for chamber pressure at 1200 psig (8.27 MPa).
Displayed response is the real part of the pressure-based combustion response function [6].
The quasi-steady burning rate may be ascertained as a function of various parameters; in this study, as a function of local static pressure p and core flow velocity u (erosive burning component), such that: While some solid propellants can also exhibit a strong sensitivity to normal acceleration as relates to burning [5], this component will be left out of the present study, for brevity.The pressure-based burning component may be found through de St. Robert's empirical law: The coefficient C will tend to vary somewhat with the propellant's initial temperature T i , getting lower as T i gets lower in value.The flow-based erosive burning component (negative and positive) is established through the following expression [7]: where at lower flow speeds, the negative component resulting from a stretched combustion zone thickness (δ r > δ o ) may cause an appreciable drop in the base burning rate r o , which in this study is the pressure-dependent burning component.The stetching of the flame zone at low speed may be viewed as being the result of a laminar-like sliding process of the local axial flow in the boundary layer acting to extend and curve the path of a representative particle moving up from the burning surface towards the flame front, such that the effective reactive length is increased.From [7], in the low-speed regime, the following expression may be applied, such that the core flow Darcy-Weisbach friction factor f is below the limit value f lim at which point negative erosive burning is no longer in effect: The parameter K δ is a shear layer coefficient, whose value set at 2600 m −1 , along with a value for f lim of 2.5 × 10 −4 , produced a good comparison to experimental data for various propellants and motors [7].At higher flow speeds, the positive erosive burning component r e , which can be established from a phenomenological convective heat feedback premise [7], should dominate: For practical modelling purposes, the value for net surface heat release term H s as specifically employed in Equation 12 is usually small in comparing well to available experimental quasi-steady erosive burning data [7], and thus is commonly set to zero when calculating r e (as is done in this study).The comparable H s term in Equation 5, for transient combustion response, is treated separately, and may have a significant positive or negative value to better compare to transient experimental data.One can find the effective convective heat transfer coefficient h of Equation 12 as a function of the zero-transpiration (zero blowing) value h* and overall propellant burning rate r b [7]: In turn, the value for h* may be found as a function of the zero-transpiration Darcy-Weisbach friction factor f*, via Reyonolds analogy [7]: In the realm of fully-developed turbulent flow, Colebrook's [11] semi-empirical expression for f* has proven accurate over a wide range of Reynolds numbers: Through Equation 15, one can see the correlation of propellant surface roughness height  and grain port diameter d p in influencing the magnitude of flow-dependent erosive burning.With composite propellants using ammonium perchlorate crystals, roughness height can commonly be correlated to the AP crystal size.
With respect to the solid propellant's burning surface temperature T s , one has the option of treating it as constant, or allowing for its variation, depending on the phenomenological approach being taken, and solution closure requirements, for estimating the local burning rate r b [6].In practice, T s 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 [12]: where  is the universal gas constant, E as is the activation energy, and A s is the Arrhenius coefficient.
In this example, the effect of local static pressure p is also included in the Arrhenius equation, i.e., n s ≠ 0, as occasionally done in the literature [12,13].Based on good comparisons in general to experimental data as reported in [6], a constant T s was employed in the Z-N phenomenological combustion model, for the present investigation.

Results and Discussion
The characteristics of the baseline reference motor for this study are listed in Table 1.In the case where the given propellant characteristic was not directly measurable or available, such as mean propellant burning surface temperature T S , assigned values are estimates, for the most part taken from the open or confidential literature, for comparable propellants.Some of the structural properties listed in Table 1, such as damping values, are estimations arising from previous studies involving comparisons to experimental vibration data.The motor, based in large measure on a similar experimental motor [3], is a smaller cylindrical-grain design with an aluminum casing and static-test steel sleeve, with a relatively large length-to-diameter ratio.The motor at the time of pulsing has a moderate port-to-throat area ratio, with a considerable propellant web thickness remaining.As noted earlier, the focus of the present study will be on cases involving single-phase flow (no substantial inert or reactive particle loading).The predicted frequency response for the AP/HTPB propellant at three different settings for the net surface heat release value may be viewed in Figure 3 (positive value for ΔH s , exothermic heat release).The general response is given in terms of the nondimensional limit magnitude M  , defined by where the reference burning rate r b,o in this case is the motor's approximate mean burn rate at the point of pulsing (1.27 cm/s).The propellant's resonant frequency f r is set via the value of K b (20,000 s −1 ) to be on the order of 1 kHz (a value within the range of what might be expected for this type of composite propellant at that base burning rate).This value for f r is in fact relatively close to the fundamental longitudinal acoustic frequency f 1L of the combustion chamber, in providing examples later in this paper that are close to the worst-case scenario for susceptibility to axial combustion instability symptoms.Figure 4 illustrates the effect of changing the base burning rate, whereby a lower base rate increases the propellant's transient response, and lowers the resonant frequency thereof.The selected base burning rate values used for Figure 3 are relevant to results presented later in this paper.Pertinent key factors can be varied in a series of test simulations, and evaluated for any identifiable trend as relates to the prediction of axial combustion instability symptoms in actual motors.Figure 5 illustrates a segment of a simulated reference motor firing's head-end pressure-time profile, showing a sustained axial wave system after the initial 2-atm pulse disturbance is introduced.In this case (and for subsequent runs, in order to narrow the scope of the present paper), the effect of acceleration [5] on the combustion process has been nullified.The sustained pressure wave is not shock-fronted [5], and one can note the apparent absence of an appreciable dc shift (rise in base chamber pressure).Moving to a somewhat higher base pressure-dependent burning rate (approx.1.5 cm/s vs. 1.27 cm/s) by increasing de St. Robert coefficient C (where r b,o = Cp n ) , one can clearly see the result on the pressure-wave profile as provided in Figure 6, for the case of acceleration being nullfied.After the intial disturbance, the axial pressure wave essentially dies out.At least two factors are at play in reducing the wave strength: (1) the combustion response resonant frequency f r has been shifted higher, away from the acoustic chamber frequency f 1L , and (2) the combustion response magnitude has been lowered at f 1L , due to the higher base burning rate, as reflected by Figure 4. Moving to a lower base pressure-dependent burning rate (approx.0.93 cm/s vs. 1.27 cm/s) by lowering de St. Robert coefficient C, one can observe the reverse from that seen earlier.The result is shown in Figure 7.The limit pressure-wave magnitude later into the firing is now approaching 5 MPa, as compared to about 1.4 MPa for the case of Figure 5.The 1L (fundamental longitudinal acoustic resonance) axial pressure wave becomes shock-fronted as observed in Figure 7(b), and there is now also a substantial dc shift, even though no acceleration effects are being included.Thus, although the lower base rate will shift the combustion response resonant frequency to a value lower than f 1L , the combustion response magnitude is heightened sufficiently at f 1L to produce a substantial increase in limit pressure wave amplitude Δp w , relative to the earlier cases.In addition, because erosive burning becomes more pronounced at lower base burning rates, it is likely also contributing to this increase in wave strength and in turn the dc shift, noting that there will be a substantial induced flow behind the traveling axial pressure wave's compression front.With respect to erosive burning, the propellant surface roughness is set at 400m for this case, a relatively high value that will certainly encourage a higher erosive burning component [7].By comparison, looking at the above case but with a lower propellant surface roughness of 100 m, one can observe in Figure 8 a lower limit pressure wave magnitude (1.5 MPa), and significantly diminished dc shift.However, unlike the case of Figure 5, the pressure wave does retain a shock front (albeit a relatively weak one), as may be seen in Figure 8  At this juncture, returning to the reference propellant surface roughness setting, one can now look at the effects of scaling the motor up in size, beginning with doubling the length of the reference motor's propellant grain.In order to retain the same reference base chamber pressure of about 9 MPa, one can lower the de St. Robert coefficient C down to 0.02 cm/s•(kPa) n .The resulting mean pressure dependent burning rate of about 0.5 cm/s is relatively low as a result, bringing the motor well into the domain of high combustion instability (c.i.) sensitivity, even with acceleration effects nullified.In order to reduce the resulting c.i. response somewhat as shown by Figure 9, the setting for the net surface heat of reaction H s has been lowered to 30,000 J/kg (from 150,000).The resonant transient combustion response frequency of 500 Hz corresponds closely to the acoustic 1L frequency for the 1-m long motor (peak response M  of 7.675), to further exacerbate the result (high, and growing, limit shock-fronted pressure wave magnitude, and high corresponding dc shift).Like the result of Figure 7, the limit pressure wave magnitude grows with time, as the port diameter increases.One factor likely contributing to this trend is the increasing port-to-throat area ratio producing a stronger reflected shock wave, even though the erosive burning component induced by the gas velocity behind the 1L shock wave front may be less with the higher port diameter.In other cases, one might in fact observe the reverse trend, due (at least in part) to the dominant effect of the lowered erosive burning with the increasing port diameter.The result of further doubling the length of the above motor, to approximately 2 m, but also increasing the starting port diameter and nozzle throat diameter to retain the same approximate base chamber pressure of 9 MPa, may be seen in Figure 10.The limit 1L pressure wave magnitude and dc shift are approximately one-quarter that seen earlier in the simulated firing for the smaller motor of Figure 9, with the quasi-equilibrium limit wave magnitude not growing nearly as much with time as seen in Figure 9.Given that the acoustic f 1L for the present motor is now approximately 250 Hz (as compared to resonant combustion f r being 500 Hz, at a base burning rate of approximately 0.5 cm/s, as before), one thus sees the effect of a reduced combustion response M  of 4.47 at 250 Hz.A further increase in motor size, as illustrated in Figure 11, shows a further reduction in the magnitude of the c.i. symptoms.The limit 1L pressure wave magnitude and corresponding dc shift are now about half the previous motor's, given the reduced combustion response M  of 2.77 at 125 Hz.While not shown, it can be noted that the breadth of the 1L wave (from shock front down the expansion tail to a return to the base pressure level, relative to the full cycle period) is thinner here (about one-third of the cycle period), as compared to earlier cases like Figure 7(b), where the breadth is on the order of half the cycle period.H s = 30,000 J/kg; p d = 2 atm;  = 400 μm), acceleration nullified.H s = 30,000 J/kg; p d = 2 atm;  = 400 μm), acceleration nullified.
A further increase in motor size, to a propellant grain length of 10 m and a starting port diameter of 70 cm, produces the interesting pressure-time profile as presented in Figure 12.The dc shift is substantially lower, about one-third that of the previous motor of Figure 11.The limit 1L pressure wave magnitude is about one-half that of Figure 11, but as revealed by Figure 12(b), the dominant 1L wave front is accompanied by a transitioning set of one, and for periods of time, two accompanying superimposed secondary waves that appear to be moving through the system at a different speed than the 1L wave front.The appearance of a multiple wave system has been observed in the past in experimental test firings [2].While in some cases a 2L wave system can be attributed to the fact that resonant combustion f r is close to f 2L [5], in this case it is more likely that the wave system is trying to find more energy at the higher harmonic frequencies where M  values are higher.The breadth of the dominant 1L pressure wave is quite thin, about one-fifth of the cycle period.The reduced combustion response, at the chamber fundamental axial acoustic frequency of 52 Hz, is 1.81.Moving to a larger motor having a propellant grain length of 20 m and a starting port diameter of 1.4 m, the dc shift appears to be almost absent in Figure 13.At the fundamental axial frequency of 26 Hz for the chamber, the 1L pressure wave magnitude is also substantially smaller than the previous case, for a combustion response of 1.41 at that frequency.As revealed by Figure 13(b), there is still some evidence of secondary waves transitioning through the system, although these waves are not nearly as significant relative to the 1L wave front, compared to the previous case of Figure 12(b).The breadth of the dominant 1L wave is on the order of one-eighth of a cycle period.Finally, at a grain length of 40 m and a starting cylindrical port diameter of 2.8 m, the motor appears to be inherently stable with respect to nonlinear axial instability symptom development.The result for Figure 14, with acceleration effects on combustion nullified, shows the case of a stronger initiating pressure disturbance of 24 atm (vs. 2 atm), with the motor still demonstrating its resistance to c.i. symptom development.The combustion response at an axial 1L acoustic chamber frequency of 13 Hz is 1.17.With acceleration effects included, with a proportionally sized motor casing wall (without a static test steel sleeve) surrounding the propellant grain, the result is similar.This is not surprising, given radial acceleration levels are quite low with such a large motor structure.In order to demonstrate the appearance of some c.i. symptom behavior at the present motor size, the net surface heat of reaction of the propellant was increased to 100,000 J/kg (from 30,000), with the result presented in Figure 15.One observes the presence of a dc shift that is a significant fraction of the limit pressure wave magnitude.As revealed by Figure 15(b), a two-wave (2L) system (with secondary superimposed wave activity) is present (which would explain the relatively larger dc shift, with the enhanced induced gas velocities of two dominant waves, plus additional secondary waves, producing a higher erosive burning component, vs. one dominant wave).The combustion response M  is 1.43 for this propellant, at 26 Hz, which for a 2L wave system is sufficient for sustained existence, while at 13 Hz, there would have been insufficient energy ( M  < 1.2) for a 1L system.The breadth of each dominant 2L wave is on the order of less than one-tenth of a 1L cycle period.

Conclusions
The implications of such factors as motor size on nonlinear axial combustion instability symptom development have been demonstrated by the example numerical simulation results presented in this study of cylindrical-grain solid rocket motors.The results presented in this paper confirm the importance of motor scale and other factors in augmenting or mitigating instability symptoms, as these factors work in conjunction with the influence of the particular propellant's transient frequency-dependent combustion response.For example, one clear trend is indicated as the fundamental acoustic chamber axial frequency f 1L becomes progressively smaller than the solid propellant's peak resonant combustion response frequency f r as one moves to larger motor lengths, namely the sustained pressure wave magnitude becomes less, as does the primary axial wave's effective width.At smaller pressure wave magnitudes, the example results suggest that one may also observe a one-wave system evolve into a two-or three-wave system as the quasi-equilibrium state of choice.Undoubtedly, further work remains to be done, both experimentally and computationally, in establishing a more complete understanding of the various mechanisms involved in driving instability symptoms in these motors, in both the axial and transverse directions.

Figure 1 .
Figure 1.Reference SRM model setup, for static test firing in laboratory.

Figure 3 .
Figure 3. Frequency response of reference propellant (K b = 20,000 s −1 ; differing values for ΔH s ; r b,o = 1.27 cm/s) in terms of nondimensional limit magnitude.

Figure 4 .
Figure 4. Frequency response of reference propellant (K b = 20,000 s −1 ; H s = 150,000 J/kg; differing values for r b,o ) in terms of nondimensional limit magnitude.