# Multisized Inert Particle Loading for Solid Rocket Axial Combustion Instability Suppression

Authors: David R. Greatrix
Journal: International Journal of Aerospace Engineering
 (2012)

---

## Abstract

In the present investigation, various factors and trends, related to the usage of two or more sets of inert particles comprised of the same material (nominally aluminum) but at different diameters for the suppression of axial shock wave development, are numerically predicted for a composite-propellant cylindrical-grain solid rocket motor. The limit pressure wave magnitudes at a later reference time in a given pulsed firing simulation run are collected for a series of runs at different particle sizes and loading distributions  and mapped onto corresponding attenuation trend charts. The inert particles’ presence in the central core flow is demonstrated to be an effective means of instability symptom suppression, in correlating with past experimental successes in the usage of particles. However, the predicted results of this study suggest that one needs to be careful when selecting more than one size of particle for a given motor application.


---

## Body

## 1. Introduction

Over the last number of decades, a multitude of research efforts have been directed towards understanding the physical mechanisms, or at least the surrounding factors, behind the appearance of symptoms typically associated with nonlinear axial combustion instability in solid-propellant rocket motors (SRMs). The principal symptoms are the presence within the motor chamber of stronger finite-amplitude traveling axial pressure waves that may be shock fronted, commonly (although not always) accompanied by some degree of base chamber pressure rise (dc shift). Note that low-magnitude pressure waves due to vortex shedding from segmented/gapped components in the motor chamber are not included (here) in this more traditional category of nonlinear axial instability. Studies of nonlinear axial combustion instability have ranged from numerous experimental test firing series on the one hand [1–3], and linear/nonlinear acoustic theory modeling on the other (largely, the analysis producing frequency-based standing wave solutions for a given chamber geometry, but without some useful quantitative information) [4–7]. On occasion, researchers have employed a numerical modeling approach, to work towards a more comprehensive quantitative understanding of the physics involved (the numerical model producing a traveling wave solution to a limit wave amplitude and corresponding small or larger dc shift, typically a time-based result evolving from an initial pulse disturbance introduced into the chamber flow) [8, 9]. Available computational power and associated result turnaround times commonly forced some simplifications in the given numerical model.The motivation for the experimental, analytical, and numerical studies noted above was and is of course to bring this better understanding to bear in more precisely suppressing, if not eliminating, these axial instability symptoms. For example, it has been long known that inert (nonreactive) or reactive particles in the internal core flow can help to suppress axial combustion instability symptoms [10–12]. As pointed out by Blomshield [13], in his wide-ranging review of a number of cases of different motors experiencing combustion instability over the years, it is not always clear as to the quantity of particle loading (and corresponding particle size) that is needed to adequately suppress the given symptoms, if that is the suppression technique being exploited. An additional potential complication, pointed out by Waesche [14], is that it is not always clear that the effectiveness of particle or additive loading is due entirely to particle/drag effects within the central core flow, or in fact in part or in whole due to altered combustion response of the propellant, as a result of the presence of the particle/additive at the burning surface before entering the core flow region. Waesche suggests that this effect may be more readily observed for reactive particles, rather than inert ones, considering the heat transfer effects in the solid phase [14]. Given this background from past experimental observations, it would be advantageous to have a predictive numerical simulation model that would help establish the particle loading/sizing requirement for a given SRM, in this case for using inert particles of differing sizes.An effective numerical model combines the effects of the unsteady one- or two-phase flow, the transient combustion process, and the structural dynamics of the surrounding propellant/casing structure. A case study reported by Blomshield [13], where the changing of a heavyweight static-test motor casing to a flightweight casing structure led to the appearance of combustion instability symptoms, provides one motivation for inclusion of structural effects in the numerical model. Experimental observations by Dotson and Sako [15] on in-flight fluid-structure interaction effects lend further weight in this regard.In the present investigation, an updated numerical model incorporating the above attributes is used in the prediction of the unsteady instability-related behavior in a cylindrical-grain motor  and allows for an evaluation of the corresponding effectiveness of using two or more sets of inert spherical particles (same material (nominally aluminum), differing diameters) in suppressing instability symptoms. While aluminum as a common solid propellant fuel addition is reactive (noninert) in practice (and its burning and other behavior at and away from the propellant surface in the central flow may have a significant influence on the given SRM’s combustion stability), the properties of aluminum are assumed for the inert particles in this study so as to allow for comparison to the results of future studies where the aluminum particles are modelled as reactive. In practice, one can note that inert particles composed of such materials as aluminum oxide (which forms from the combustion of aluminum and oxygen) or zirconium carbide do see usage for combustion stabilization purposes. The present study is a followon to the study reported in [16] (where the use of a single set of inert spherical particles is examined). In practice, one might see the use of two or more particle sizes in a given motor. This is sometimes done to target two or more different pressure wave frequencies that have been identified as problematic (e.g., one longitudinal and two transverse if using three particle sets). The Dobbins-Temkin correlation [17] indicates that the best particle diameter dopt for suppression is a function of the inverse square root of the target frequency; that is, likely a smaller diameter particle is more effective at a higher transverse frequency, versus a lower axial frequency, everything else being equal(1)dopt=9μπfρm.In the present paper, the focus for presented results will largely be on those cases where the transient burning response of the propellant is the primary mechanism for sustaining appreciable traveling pressure waves in the combustion chamber. A few additional results will illustrate the effect of normal acceleration (through radial vibration) as a complementary mechanism acting on the transient combustion process.


## 2. Numerical Model

A simplified schematic diagram of the physical system of an SRM, that is placed on a static test stand, is provided in Figure1. In this case, the cylindrical-grain motor is free to vibrate radially without any external constraint (i.e., only constrained as indicated by the thick steel static-test sleeve surrounding the aluminum flightweight motor casing), while axial motion is constrained to a large degree by a thrust-measuring load cell at the lefthand boundary. Under normal (nominal) quasiequilibrium operating conditions, the internal gas flow (or gas-particle flow, if two-phase) moves smoothly from the burning propellant surface into the central core flow, heading downstream to eventually pass through and beyond the exhaust nozzle.Figure 1
Schematic diagram of sleeved cylindrical-grain SRM, showing referencex-direction.


### 2.1. Equations Relevant to Two-Phase Flow in SRM

One definesb as the nonequilibrium sound speed of a 2-phase mixture, which for a single set of monodisperse particles within the gas can be estimated via [18](2)b=a[1+βg(Cm/Cp)(Tp/T)X+αgβg(γ-1)(up/u)(1-(up/u))]1/2,
where 𝒳 denotes αg2(1+βg(up/u))(1+γβg(Cm/Cp)(Tp/T)). The gas phase void fraction αg is defined by
(3)αg=V-VpV=1-VpV.
For a lower particle loading, one can assume that αg is close to unity in value (Vp is the volume occupied by the particles in an elemental total volume of V). The particle-gas mass flux ratio βg is stipulated as
(4)βg=ρpupρuαg.
One can define an average two-phase density in a given volume as(5)ρ2ph=mV≈N⋅mp+ρVV=ρpV+ρVV=αpρ2ph+ρ,
where αp is the particle-loading fraction in the flow, and N is the number of particles of average mass mp in the elemental volume above, then one can show the correlation between particle density ρp and gas density ρ(6)ρp=N⋅mpV=αpρ2ph=αp1-αp⋅ρ.When one has two or more sets of particles of differing sizes (comprised of the same material; thus, the same solid specific heatCm), one can use the following correlations:(7)ρp=∑i=1NsetNimp,iV=∑i=1Nsetρp,i,αp=∑i=1Nsetαp,i,βg=ṁpṁg=∑i=1Nsetρp,iup,iρuαg,βgCmCpTpT=CmCpT⋅∑i=1Nsetρp,iup,iTp,iρuαg=Σ1,βgupu=∑i=1Nsetρp,iup,i2ρu2αg=Σ2,βgup2u2=∑i=1Nsetρp,iup,i3ρu3αg=Σ3,
so that the nonequilibrium sound speed can be estimated via(8)b=a[1+Σ1αg2(1+Σ2)(1+γΣ1)+αg(γ-1)(Σ2-Σ3)]1/2.
Under nominal flow conditions, u/b would be unity at the nozzle throat.The effect of particle mass loadingΔMp into a solid propellant, that originally was of solid density ρs,o, on the loading mass fraction into a solid volume V is given by(9)αp,s=ΔMpΔMs,o+ΔMp=ρmVmρs,o(V-Vm)+ρmVm=1(ρs,o/ρm)(V/Vm-1)+1.
The new effective solid propellant overall density becomes
(10)ρs,new=ΔMp+ΔMs,oV=ρs,o(V-Vm)+ρmVmV=ρs,o(1-VmV)+ρmVmV.
By substitution, one can show that(11)ρs,new=ρs,o1-αp,s+αp,s(ρs,o/ρm).


### 2.2. Equations of Motion-Governing Two-Phase Flow

The equations of motion describing the nonsteady core flow within the SRM must be solved in conjunction with the local pyrolysis raterb of the solid propellant, and the surrounding structure’s instantaneous geometric deformation. As pertains to the present study of a small motor having a larger length-to-diameter ratio, the quasi-one-dimensional hydrodynamic conservation equations for the axial gas flow are given below(12)∂ρ∂t+∂(ρu)∂x=-1A∂A∂xρu+(1-αp)ρs4rbd-(4rbd+κ)ρ,∂(ρu)∂t+∂∂x(ρu2+p)=-1A∂A∂xρu2-(4rbd+κ)ρu-ρal-∑i=1NsetρpimpiDi,∂(ρE)∂t+∂∂x(ρuE+up)=-1A∂A∂x(ρuE+up)-(4rbd+κ)ρE+(1-αp)ρs4rbd(CpTf+vf22)-ρual-∑i=1Nsetρpimpi(upiDi+Qi).
Here, the total specific energy of the gas is defined for an ideal gas as E=p/[(γ-1)ρ]+u2/2. The corresponding equations of motion for an ith inert (nonburning) particle set within the axial flow may be found from(13)∂ρpi∂t+∂(ρpiupi)∂x=-1A∂A∂xρpiupi+αpiρs4rbd-(4rbd+κ)ρpi,∂(ρpiupi)∂t+∂(ρpiupi2)∂x=-1A∂A∂xρpiupi2-(4rbd+κ)ρpiupi-ρpial+ρpimpiDi,∂(ρpiEpi)∂t+∂(ρpiupiEpi)∂x=-1A∂A∂x(ρpiupiEpi)-(4rbd+κ)ρpiEpi+αpiρs4rbd(CmTf+vf22)-ρpiupial+ρpimpi×(upiDi+Qi).
Here, the total specific energy of a local grouping of particles from an ith set is given by Epi=CmTpi+upi2/2, where Tpi is the mean temperature of that group. As outlined in [19], the viscous interaction between the gas and a particle from the ith particle set is represented by the drag force Di, and the heat transfer from the core flow to a particle from the ith set is defined by Qi. In the case of drag between the gas and a representative spherical particle at a given axial location, one notes that(14)Di=πdmi28Cdρ(u-upi)|u-upi|,
where Cd is the drag coefficient for a sphere in a steady flow with low-flow turbulence (determined as function of relative Reynolds number, relative flow Mach number, and temperature difference between the particle and the gas). In the case of heat transfer from the core flow to a representative particle at a given axial location, the following applies:(15)Qi=πdmik⋅Nu⋅(T-Tpi),
where the Nusselt number Nu can be found as a function of Prandtl and relative Reynolds number for a sphere of mean diameter dmi. One will need to solve   (13)  for each of Nset  particle sets as part of the calculation process, where for dual or triple particle set loading, Nset has a value of 2 or 3 in the present study.Longitudinal accelerational appears in the gas and particle momentum and energy equations as a body force contribution within a fixed Eulerian reference (fixing of x= 0 to motor head end, xpositive moving right on structure as per Figure 1; acceleration of local surrounding structure rightward is designated positive al)  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, explicit, and finite-volume random-choice method (RCM) approach [19, 20]. The RCM solver employs a Riemann-solution technique noted for low artificial dispersion with time of wave activity in tubes, and so forth. The equations of motion of the gas and particles will be solved over a given time step Δt (on the order of 1 × 10−7 s for the present study, given the motor solution node allocation in the axial direction from head end to nozzle exit plane), in sequence with additional equations for structural motion and propellant burning rate as described below.


### 2.3. Equations for Structural Motion

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, and material properties) [13, 15, 21]. The level of sophistication required for modeling the motor structure (propellant, casing, static-test sleeve, and nozzle) and applicable boundary conditions (load cell on static test stand) can vary, depending on the particular application and motor design. Loncaric et al. [22] and Montesano et al. [23] employed a finite-element approach towards the structural modeling of the given motor configuration. In the present study, a cylindrical-grain configuration allows for a simpler finite-difference approach via thick-wall theory, as reported in [20, 24]. The radial deformation dynamics of the propellant/casing/sleeve are modeled by a series of independent ring elements along the length of the motor. Axial motion along the length of the structure is modeled via beam theory, and bounded by the spring/damper load cell at the motor’s head end. Viscous damping is applied in the radial and axial directions. Reference structural properties are assumed for an ammonium-perchlorate/hydroxyl-terminated polybutadiene (AP/HTPB) composite propellant surrounded by an aluminum casing and steel sleeve. For greater accuracy, some properties like the propellant/casing/sleeve assembly’s natural radial frequency may be predetermined via a finite-element numerical solution, rather than via theoretical approximations [24].


### 2.4. Equations for Propellant Burning Rate

With respect to transient, frequency-dependent burning rate modeling, the Z-N (Zeldovich-Novozhilov) solid-phase energy conservation approach used in the present simulation program may be represented by the following time-dependent temperature-based relationship [25]:(16)rb*=rb,qs-1(Ts-Ti-ΔHs/Cs)∂∂t∫-∞0ΔTdy,
where rb,qs is the quasisteady burning rate (value for burning rate as estimated from steady-state information for a given set of local flow conditions), Ti is the initial propellant temperature, and in this context, ΔT=T(y,t)-Ti is the temperature distribution in moving from the burning propellant surface at y=0 (and T=Ts) to that spatial location in the propellant where the temperature reaches Ti. One may note at this juncture the inclusion of a net surface heat release term, ΔHs, in the calculations. 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, ΔyFo=(2αsΔt)1/2, which is a function of the chosen time increment Δt [25]. The time increment itself must be coordinated between the flow and structural model solution systems [23].In   (16), rb* 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. The actual instantaneous burning rate rb may be found as a function of rb* through the empirical rate-limiting equation [25](17)drbdt=Kb(rb*-rb).The rate-limiting coefficient Kb effectively damps the unconstrained burning rate rb* when for a finite time increment Δt(18)Kb<1Δt.
In the present approach, the surface-thermal gradient is free to find its own value at a given instant. One can argue that the use of  (17) or some comparable damping function, while empirical, parallels the approach taken by past researchers in using a stipulated surface-thermal gradient; both approaches act to constrain the exchange of energy through the burning surface interface, allow for some variability in better comparing to a given set of experimental data, and prevent so-called burning-rate “runaway” (unstable divergence of rb with time) [26]. As discussed in [25], the use of Kb at a set value does allow for a converged solution that is independent of the increment size for Δt and Δy, as long as one respects the Fourier stability requirement noted earlier.The quasisteady burning raterb,qs may be ascertained as a function of various parameters; in this study, as a function of local static pressure p, core flow velocity u (erosive burning component), and normal/lateral/longitudinal acceleration, such that:(19)rb,qs=rp+ru+ra.
The pressure-based burning component may be found through de St. Robert’s law(20)rp=Cpn.
The flow-based erosive burning component (negative and positive) is established through the following expression [27]:(21)rb=rbro|δr⋅ro+re,
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 burn rate ro, while at higher flow speeds, the positive erosive burning component re, established from a convective heat feedback premise [27], should dominate:(22)re=h(Tf-Ts)ρs[Cs(Ts-Ti)-ΔHs].
For the above case, where the base burning rate ro is a function of the other mechanisms (pressure and acceleration), one finds the velocity-based component of burn rate from   (21) via ru=rb-ro. At higher flow speeds, ru becomes equivalent to re. The effect of normal acceleration an resulting from radial propellant/casing/sleeve vibration may be determined via [28] (23)rb=[Cp(Tf-Ts)Cs(Ts-Ti)-ΔHs](rb+Ga/ρs)exp⁡[Cpδo(ρsrb+Ga)/k]-1,
where the compressive effect of normal acceleration and the dissipative effect of steady or oscillatory longitudinal (or lateral, if say for a star grain configuration) acceleration al are stipulated through the accelerative mass flux Ga(24)Ga={anprbδoRTfrorb}ϕ=0∘cos⁡2ϕd.
Note that the longitudinal/lateral-acceleration-based displacement orientation angle ϕd is greater than the nominal acceleration vector orientation angle (ϕ; zero when only normal acceleration an relative to the burning propellant surface is present) [28]. One should also note that an is negative when acting to compress the combustion zone, and treated as zero when directed away from the zone. For the above case, where the base burning rate ro is a function of the other flow mechanisms (pressure and core flow), one finds the acceleration-based component of burning rate from   (23) via ra=rb-ro.With respect to the burning surface temperatureTs, 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 [25]. While in the past a number of estimation models might have used a constant value for Ts, more recently the usage of a variable Ts has become prevalent. However, based on good comparisons in general to experimental data as reported in [25], a constant Ts was employed in the present Z-N-based phenomenological numerical combustion model, for the present investigation.


## 2.1. Equations Relevant to Two-Phase Flow in SRM

One definesb as the nonequilibrium sound speed of a 2-phase mixture, which for a single set of monodisperse particles within the gas can be estimated via [18](2)b=a[1+βg(Cm/Cp)(Tp/T)X+αgβg(γ-1)(up/u)(1-(up/u))]1/2,
where 𝒳 denotes αg2(1+βg(up/u))(1+γβg(Cm/Cp)(Tp/T)). The gas phase void fraction αg is defined by
(3)αg=V-VpV=1-VpV.
For a lower particle loading, one can assume that αg is close to unity in value (Vp is the volume occupied by the particles in an elemental total volume of V). The particle-gas mass flux ratio βg is stipulated as
(4)βg=ρpupρuαg.
One can define an average two-phase density in a given volume as(5)ρ2ph=mV≈N⋅mp+ρVV=ρpV+ρVV=αpρ2ph+ρ,
where αp is the particle-loading fraction in the flow, and N is the number of particles of average mass mp in the elemental volume above, then one can show the correlation between particle density ρp and gas density ρ(6)ρp=N⋅mpV=αpρ2ph=αp1-αp⋅ρ.When one has two or more sets of particles of differing sizes (comprised of the same material; thus, the same solid specific heatCm), one can use the following correlations:(7)ρp=∑i=1NsetNimp,iV=∑i=1Nsetρp,i,αp=∑i=1Nsetαp,i,βg=ṁpṁg=∑i=1Nsetρp,iup,iρuαg,βgCmCpTpT=CmCpT⋅∑i=1Nsetρp,iup,iTp,iρuαg=Σ1,βgupu=∑i=1Nsetρp,iup,i2ρu2αg=Σ2,βgup2u2=∑i=1Nsetρp,iup,i3ρu3αg=Σ3,
so that the nonequilibrium sound speed can be estimated via(8)b=a[1+Σ1αg2(1+Σ2)(1+γΣ1)+αg(γ-1)(Σ2-Σ3)]1/2.
Under nominal flow conditions, u/b would be unity at the nozzle throat.The effect of particle mass loadingΔMp into a solid propellant, that originally was of solid density ρs,o, on the loading mass fraction into a solid volume V is given by(9)αp,s=ΔMpΔMs,o+ΔMp=ρmVmρs,o(V-Vm)+ρmVm=1(ρs,o/ρm)(V/Vm-1)+1.
The new effective solid propellant overall density becomes
(10)ρs,new=ΔMp+ΔMs,oV=ρs,o(V-Vm)+ρmVmV=ρs,o(1-VmV)+ρmVmV.
By substitution, one can show that(11)ρs,new=ρs,o1-αp,s+αp,s(ρs,o/ρm).


## 2.2. Equations of Motion-Governing Two-Phase Flow

The equations of motion describing the nonsteady core flow within the SRM must be solved in conjunction with the local pyrolysis raterb of the solid propellant, and the surrounding structure’s instantaneous geometric deformation. As pertains to the present study of a small motor having a larger length-to-diameter ratio, the quasi-one-dimensional hydrodynamic conservation equations for the axial gas flow are given below(12)∂ρ∂t+∂(ρu)∂x=-1A∂A∂xρu+(1-αp)ρs4rbd-(4rbd+κ)ρ,∂(ρu)∂t+∂∂x(ρu2+p)=-1A∂A∂xρu2-(4rbd+κ)ρu-ρal-∑i=1NsetρpimpiDi,∂(ρE)∂t+∂∂x(ρuE+up)=-1A∂A∂x(ρuE+up)-(4rbd+κ)ρE+(1-αp)ρs4rbd(CpTf+vf22)-ρual-∑i=1Nsetρpimpi(upiDi+Qi).
Here, the total specific energy of the gas is defined for an ideal gas as E=p/[(γ-1)ρ]+u2/2. The corresponding equations of motion for an ith inert (nonburning) particle set within the axial flow may be found from(13)∂ρpi∂t+∂(ρpiupi)∂x=-1A∂A∂xρpiupi+αpiρs4rbd-(4rbd+κ)ρpi,∂(ρpiupi)∂t+∂(ρpiupi2)∂x=-1A∂A∂xρpiupi2-(4rbd+κ)ρpiupi-ρpial+ρpimpiDi,∂(ρpiEpi)∂t+∂(ρpiupiEpi)∂x=-1A∂A∂x(ρpiupiEpi)-(4rbd+κ)ρpiEpi+αpiρs4rbd(CmTf+vf22)-ρpiupial+ρpimpi×(upiDi+Qi).
Here, the total specific energy of a local grouping of particles from an ith set is given by Epi=CmTpi+upi2/2, where Tpi is the mean temperature of that group. As outlined in [19], the viscous interaction between the gas and a particle from the ith particle set is represented by the drag force Di, and the heat transfer from the core flow to a particle from the ith set is defined by Qi. In the case of drag between the gas and a representative spherical particle at a given axial location, one notes that(14)Di=πdmi28Cdρ(u-upi)|u-upi|,
where Cd is the drag coefficient for a sphere in a steady flow with low-flow turbulence (determined as function of relative Reynolds number, relative flow Mach number, and temperature difference between the particle and the gas). In the case of heat transfer from the core flow to a representative particle at a given axial location, the following applies:(15)Qi=πdmik⋅Nu⋅(T-Tpi),
where the Nusselt number Nu can be found as a function of Prandtl and relative Reynolds number for a sphere of mean diameter dmi. One will need to solve   (13)  for each of Nset  particle sets as part of the calculation process, where for dual or triple particle set loading, Nset has a value of 2 or 3 in the present study.Longitudinal accelerational appears in the gas and particle momentum and energy equations as a body force contribution within a fixed Eulerian reference (fixing of x= 0 to motor head end, xpositive moving right on structure as per Figure 1; acceleration of local surrounding structure rightward is designated positive al)  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, explicit, and finite-volume random-choice method (RCM) approach [19, 20]. The RCM solver employs a Riemann-solution technique noted for low artificial dispersion with time of wave activity in tubes, and so forth. The equations of motion of the gas and particles will be solved over a given time step Δt (on the order of 1 × 10−7 s for the present study, given the motor solution node allocation in the axial direction from head end to nozzle exit plane), in sequence with additional equations for structural motion and propellant burning rate as described below.


## 2.3. Equations for Structural Motion

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, and material properties) [13, 15, 21]. The level of sophistication required for modeling the motor structure (propellant, casing, static-test sleeve, and nozzle) and applicable boundary conditions (load cell on static test stand) can vary, depending on the particular application and motor design. Loncaric et al. [22] and Montesano et al. [23] employed a finite-element approach towards the structural modeling of the given motor configuration. In the present study, a cylindrical-grain configuration allows for a simpler finite-difference approach via thick-wall theory, as reported in [20, 24]. The radial deformation dynamics of the propellant/casing/sleeve are modeled by a series of independent ring elements along the length of the motor. Axial motion along the length of the structure is modeled via beam theory, and bounded by the spring/damper load cell at the motor’s head end. Viscous damping is applied in the radial and axial directions. Reference structural properties are assumed for an ammonium-perchlorate/hydroxyl-terminated polybutadiene (AP/HTPB) composite propellant surrounded by an aluminum casing and steel sleeve. For greater accuracy, some properties like the propellant/casing/sleeve assembly’s natural radial frequency may be predetermined via a finite-element numerical solution, rather than via theoretical approximations [24].


## 2.4. Equations for Propellant Burning Rate

With respect to transient, frequency-dependent burning rate modeling, the Z-N (Zeldovich-Novozhilov) solid-phase energy conservation approach used in the present simulation program may be represented by the following time-dependent temperature-based relationship [25]:(16)rb*=rb,qs-1(Ts-Ti-ΔHs/Cs)∂∂t∫-∞0ΔTdy,
where rb,qs is the quasisteady burning rate (value for burning rate as estimated from steady-state information for a given set of local flow conditions), Ti is the initial propellant temperature, and in this context, ΔT=T(y,t)-Ti is the temperature distribution in moving from the burning propellant surface at y=0 (and T=Ts) to that spatial location in the propellant where the temperature reaches Ti. One may note at this juncture the inclusion of a net surface heat release term, ΔHs, in the calculations. 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, ΔyFo=(2αsΔt)1/2, which is a function of the chosen time increment Δt [25]. The time increment itself must be coordinated between the flow and structural model solution systems [23].In   (16), rb* 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. The actual instantaneous burning rate rb may be found as a function of rb* through the empirical rate-limiting equation [25](17)drbdt=Kb(rb*-rb).The rate-limiting coefficient Kb effectively damps the unconstrained burning rate rb* when for a finite time increment Δt(18)Kb<1Δt.
In the present approach, the surface-thermal gradient is free to find its own value at a given instant. One can argue that the use of  (17) or some comparable damping function, while empirical, parallels the approach taken by past researchers in using a stipulated surface-thermal gradient; both approaches act to constrain the exchange of energy through the burning surface interface, allow for some variability in better comparing to a given set of experimental data, and prevent so-called burning-rate “runaway” (unstable divergence of rb with time) [26]. As discussed in [25], the use of Kb at a set value does allow for a converged solution that is independent of the increment size for Δt and Δy, as long as one respects the Fourier stability requirement noted earlier.The quasisteady burning raterb,qs may be ascertained as a function of various parameters; in this study, as a function of local static pressure p, core flow velocity u (erosive burning component), and normal/lateral/longitudinal acceleration, such that:(19)rb,qs=rp+ru+ra.
The pressure-based burning component may be found through de St. Robert’s law(20)rp=Cpn.
The flow-based erosive burning component (negative and positive) is established through the following expression [27]:(21)rb=rbro|δr⋅ro+re,
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 burn rate ro, while at higher flow speeds, the positive erosive burning component re, established from a convective heat feedback premise [27], should dominate:(22)re=h(Tf-Ts)ρs[Cs(Ts-Ti)-ΔHs].
For the above case, where the base burning rate ro is a function of the other mechanisms (pressure and acceleration), one finds the velocity-based component of burn rate from   (21) via ru=rb-ro. At higher flow speeds, ru becomes equivalent to re. The effect of normal acceleration an resulting from radial propellant/casing/sleeve vibration may be determined via [28] (23)rb=[Cp(Tf-Ts)Cs(Ts-Ti)-ΔHs](rb+Ga/ρs)exp⁡[Cpδo(ρsrb+Ga)/k]-1,
where the compressive effect of normal acceleration and the dissipative effect of steady or oscillatory longitudinal (or lateral, if say for a star grain configuration) acceleration al are stipulated through the accelerative mass flux Ga(24)Ga={anprbδoRTfrorb}ϕ=0∘cos⁡2ϕd.
Note that the longitudinal/lateral-acceleration-based displacement orientation angle ϕd is greater than the nominal acceleration vector orientation angle (ϕ; zero when only normal acceleration an relative to the burning propellant surface is present) [28]. One should also note that an is negative when acting to compress the combustion zone, and treated as zero when directed away from the zone. For the above case, where the base burning rate ro is a function of the other flow mechanisms (pressure and core flow), one finds the acceleration-based component of burning rate from   (23) via ra=rb-ro.With respect to the burning surface temperatureTs, 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 [25]. While in the past a number of estimation models might have used a constant value for Ts, more recently the usage of a variable Ts has become prevalent. However, based on good comparisons in general to experimental data as reported in [25], a constant Ts was employed in the present Z-N-based phenomenological numerical combustion model, for the present investigation.


## 3. Results and Discussion

The characteristics of the reference motor for this study are listed in Table1. The motor, based in large measure on a similar experimental motor [20] 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. 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 2 (positive value, exothermic heat release). The general response is given in terms of the nondimensional limit magnitude Ml, defined by(25)Ml=rb,peak-rb,orb,qs,peak-rb,o,
where the reference burning rate rb,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 fr is set via the value of Kb (20000 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 fr is in fact relatively close to the fundamental longitudinal acoustic frequency f1L 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.Table 1
Reference motor characteristics.
ParameterValuePropellant grain length,Lp52 cmInitial port diameter,di3.6 cmNozzle throat diameter,dt1.6 cmGrain/nozzle-conv. length ratio,Lp/Lc16 : 1Propellant specific heat,Cs1500 J/kg-KPropellant density,ρs1730 kg/m3Propellant thermal conductivity,ks0.4 W/m-KPropellant thermal diffusivity,αs1.54 × 10-7m2/sPropellant flame temperature,Tf3000 KPropellant surface temperature,Ts1000 KPropellant initial temperature,Ti294 KGas specific heat,Cp1920 J/kg-KSpecific gas constant,R320 J/kg-KGas thermal conductivity,k0.2 W/m-KGas absolute viscosity,μ8.07 × 10-5kg/m-sGas specific heat ratio,γ1.2De St. Robert exponent,n0.35De St. Robert coefficient,C0.05 cm/s-(kPa)nParticle solid density,ρm2700 kg/m3Particle specific heat, J/kg-K900 J/kg-KParticle mass fraction,αp0%Propellant elastic modulus,EA45 MPaPropellant Poisson’s ratio,νA0.497Casing inner wall radius,rm3.24 cmCasing wall thickness,hB0.127 cmCasing material density,ρB2700 kg/m3Casing elastic modulus,EB80 GPaCasing material Poisson’s ratio,νB0.33Sleeve wall thickness,hC0.47 cmSleeve material density,ρC7850 kg/m3Sleeve elastic modulus,EC200 GPaSleeve material Poisson’s ratio,νC0.30Casing/prop. rad. damping ratio,ξR0.35Casing/prop. long. damping ratio,ξL0.10Figure 2
Frequency response of  reference propellant (rb,o=1.27 cm/s, Kb=20000 s−1,  differing ΔHs ) in terms of nondimensional limit magnitude.An initial pulsed-firing simulation run was completed as a starting reference for this study, in which no particles are present or any other suppression technique being applied. In Figure3 for head-end pressure pc as a function of time, one can see that, at some point, the principal compression wave reaches its quasiequilibrium strength from an initial disturbance pressure Δpd of 2 atm, the sustained compression wave front arriving about every 1 ms, oscillating at the fundamental frequency f1L of 1 kHz. The base pressure is not appreciably elevated over the nominal operating chamber pressure. The effect of normal acceleration on the burning process (related to the radial vibration of the motor propellant/casing; see [16, 20, 23]) has been nullified for this simulation (in order to isolate frequency-dependent Z-N combustion response as the predominant instability symptom driver), a factor in reducing the development of a dc shift. One can note that the limit pressure wave magnitude (Δpw, peak to trough) is decreasing gradually with time after first reaching its quasiequilibrium level, as the cylindrical grain burns back and the base pressure rises.Figure 3
Predicted head-end pressure-time profile, reference motor (Kb=20000 s−1, ΔHs=150000 J/kg,  Δpd=2 atm, αp=0%), acceleration nullified, no particles.One can refer to Figure4 for the pressure-time profile for the same motor, but now with 5% particle loading (by mass) of inert spherical aluminum particles having a mean 10-μm diameter. Of course, in practice, the aluminum particles would in fact be reactive (burning, if sufficient reactants like oxygen or chlorine are present in the surrounding gas), and as a result, in general continually decreasing in diameter with particle surface regression as they move aft towards the nozzle. There is also the possibility of particle agglomeration, or the coming together of two or more particles. Given the scope of the present investigation, calculations for particle regression or agglomeration were not to be done; one can consider the mean inert aluminum particle diameter for a given set as a reference size, providing results which may prove useful as a guideline when one does move to inclusion of particle burning in the computational model. Observing the results of Figure 4, suppression of axial wave development after an initial 2-atm pulse is near-complete (limit magnitude of the sustained pressure wave, at 0.26 s, is about 0.045 MPa [Δpw], as compared to 1.42 MPa for the 0% loading case noted earlier [Δpw,peak], giving a nondimensional attenuation Ma, defined by(26)Ma=Δpw,peak-ΔpwΔpw,peak,
a value of 0.97, noting that a value of unity is complete suppression). Historically, suppression of high-frequency tangential and radial pressure waves in SRMs by the use of particles in the range of 1 to 3% loading by mass has been in general largely successful. In the case of axial pressure waves, the effectiveness of particles from 1% to over 20% loading in suppressing wave development has been less consistent, relative to the previously mentioned transverse cases. In the case of Figure 4, remembering that acceleration as a factor has been nullified in the combustion process, a loading of 5% at 10 μm does appear to effectively suppress axial wave development in this particular motor, at this point in its firing.Figure 4
Predicted head-end pressure wave profile, reference motor (Kb=20000 s−1, ΔHs=150000 J/kg, Δpd=2 atm), acceleration nullified, 5% 10-μm Al particle loading by mass, single particle set.In considering an example of two particle sets being used (5 and 10μm diameters for particles of the same material), in Figure 5 one sees the result of an evenly split 1 : 1 distribution of the two sets at an overall loading of 5%. The limit magnitude of the sustained pressure wave, at 0.26 s, is about 0.47 MPa, as compared to about ten times less for the uniform 10-μm-loading case noted for Figure 4, giving a nondimensional attenuation Ma value of 0.67 as compared to 0.97. In considering an example of three particle sets being used (5, 10, and 40 μm diameters for particles of the same material), in Figure 6 one can observe the result of an evenly split 1 : 1 : 1 distribution of the three sets at an overall loading of 5%. The limit magnitude of the sustained pressure wave, at 0.26 s, is about 0.57 MPa, giving an Ma value of 0.6, or a bit less than the previous two-particle example (0.67). Referring to Figure 7 [16], one can observe that, in the limit of just using 5 μm particles at a 5% overall loading, the value for Ma is around 0.52, while, for 40 μm particles at a 5% overall loading, the value for Ma is around 0.55.Figure 5
Predicted head-end pressure wave profile, reference motor (Kb=20000 s−1, ΔHs=150000 J/kg, Δpd=2 atm, αp=5%, dm1=10 μm, dm2=5 μm, αp1=2.5%, αp2=2.5%), two particle sets,      acceleration nullified.Figure 6
Predicted head-end pressure wave profile, reference motor (Kb=20000 s−1, ΔHs=150000 J/kg, Δpd=2 atm, αp=5%, dm1=10 μm,  dm2=5 μm, dm3=40 μm, αp1=1.67%, αp2=1.67%, αp3=1.66% ), three particle sets, acceleration nullified.Figure 7
Nondimensional attenuation as function of particle diameter and     loading of a single particle set, reference motor,  acceleration nullified.Let’s consider the case when vibration-induced acceleration is active as a mechanism working in conjunction with the transient response of the burning solid propellant. Referring to Figure8, allowing for the effect of vibration-induced acceleration on combustion, one has a much more active motor in the absence of particles in the flow, with a substantially bigger limit magnitude (around 10.2 MPa at 0.26 s) for the now shock-fronted axial pressure wave moving back and forth within the chamber than that seen in Figure 3, and a much more visible dc shift is present. In considering an example of three particle sets being used (5, 10, and 60 μm diameters for particles of the same material), in Figure 9 one can observe the result of an evenly split 1 : 1 : 1 distribution of the three sets at an overall loading of 1.5%. The limit magnitude of the sustained pressure wave, at 0.26 s, is about 8.97 MPa, giving an Ma value of 0.121. Referring to the single-set chart of Figure 10 [16], one can observe that in the limit of just using 5 μm particles at a 1.5% overall loading, the value for Ma is less than 0.15, using 10 μm particles at a 1.5% overall loading produces a value for Ma is less than 0.2, and for 60 μm particles at a 1.5% overall loading, the value for Ma is less than 0.2. In qualitative terms, one can state that the low-suppression attenuation observed in Figure 9 is reasonably consistent with what one would expect from the aggregate of the three individual diameter results of Figure 10  although, quantitatively, at least in this case, the 3-set attenuation seen in Figure 9 is a bit worse (lower) than any of the three in isolation.Figure 8
Predicted head-end pressure-time profile, reference motor (Kb=20000 s−1, ΔHs=150000 J/kg, Δpd=2 atm, αp=0%), acceleration active, no particles.Figure 9
Predicted head-end pressure wave profile, reference motor (Kb=20000 s−1, ΔHs=150000 J/kg, Δpd=2 atm, αp=1.5%, dm1=10 μm,  dm2=5 μm, dm3=60 μm, αp1=0.5%, αp2=0.5%, αp3 = 0.5% ), three particle sets, acceleration active.Figure 10
Nondimensional attenuation as function of particle diameter and   loading of a single particle set, reference motor,  acceleration active.One adjustment in particle size produces a significant change in the result of Figure9, as evidenced by Figure 11. Using 40 μm particles in place of the 60-μm particles at a 0.5% loading as part of the overall 1.5% loading, the limit pressure wave magnitude is significantly decreased, down to 1.22 MPa from the earlier value of 8.97, or producing an Ma of around 0.88. Referring to the single-set trend chart of Figure 10, this indicates that the system has decided to shift to the high-suppression domain by this one adjustment, moving from the low-suppression domain that the system preferred in Figure 9.Figure 11
Predicted head-end pressure wave profile, reference motor (Kb=20000 s−1, ΔHs=150000 J/kg, Δpd=2 atm, αp=1.5%, dm1=10 μm,  dm2=5 μm, dm3=40 μm, αp1=0.5%, αp2=0.5%, αp3=0.5% ), three particle sets, acceleration active.Along the lines of Figures7 and 10, one can produce charts showing trends as relates to using two sets of particles of the same material, but differing diameters. As the baseline case, Figure 12 shows an attenuation chart at three different overall particle mass loadings (αp) for various loading distributions (αp2:αp1) of 10-μm particles. Acceleration effects on combustion are nullified for these examples. The loading ratio in these examples is referenced to the 10 μm size (dm1), since as illustrated in Figure 7, it is nominally around the best size for attenuation of pressure-wave development in the single-size case. Because the particles of the two sets are the same size, there would not be any change in attenuation, as reflected by the resulting curves being horizontal lines from left to right. Figure 13 by contrast shows downward sloping lines when using a distribution of 5 and 10 μm particles at different overall loadings. Each of the 3 curves would continue rightward to their respective asymptotic limit, whereby only 5 μm particles are ultimately present (see Figure 7 for limit values of attenuation at the three loadings of the 5-μm size). A chart such as that shown by Figure 13 would also be instructive towards the case of using reactive particles that would reduce in size (due to burnback) during their lifetime in the motor chamber. The level of attenuation may not be as effective as one might expect from the baseline nominal starting value for the particle diameter. Figure 14 shows downward sloping lines when using a distribution of 20 and 10 μm particles at different overall loadings, but of a less severe nature (relative to that seen in Figure 13), given the similar effectiveness of the two particle sizes below an overall loading of 5% as per Figure 7.Figure 12
Nondimensional attenuation as function of particle diameter and loading distribution[αp2:αp1] for two sets of particles (10 μm [dm1] and 10 μm [dm2] particles), reference motor,  acceleration nullified.Figure 13
Nondimensional attenuation as function of particle diameter and loading distribution[αp2:αp1]  for two sets of particles (10 μm [dm1]  and 5 μm [dm2] particles), reference motor,  acceleration  nullified.Figure 14
Nondimensional attenuation as function of particle diameter and loading distribution[αp2:αp1] for two sets of particles (10 μm [dm1] and 20 μm [dm2] particles), reference motor, acceleration nullified.An alternative format for illustrating trends associated with using two sets of particles is provided in Figures15, 16, and 17, for overall particle mass loading percentages of 3, 4, and 5%. Here, the three graphs are restricted to evaluating an even 1 : 1 split in the distribution of the two sets, but with the particle diameters varied relative to each other. A more dramatic change in attenuation effectiveness as a function of varying the particle diameters is seen in the curves of Figure 17, for the 5% loading case, which again (referring to the previous paragraph), tends to correlate with the trends associated with Figure 7.Figure 15
Nondimensional attenuation as function of particle diameter(dm2:dm1), for a 1 : 1 loading distribution for two sets of particles at total loading of 3%, reference motor, acceleration nullified.Figure 16
Nondimensional attenuation as function of particle diameter(dm2:dm1), for a 1 : 1 loading distribution for two sets of particles at total loading of 4%, reference motor, acceleration nullified.Figure 17
Nondimensional attenuation as function of particle diameter(dm2:dm1), for a 1 : 1 loading distribution for two sets of particles at total loading of 5%, reference motor, acceleration nullified.


## 4. Concluding Remarks

A numerical evaluation of the use of two or three sets of different-sized nonburning particles within the flow as a means for suppressing axial pressure wave development has been completed for a reference cylindrical-grain composite-propellant motor, in cases where the transient burning response of the propellant is the primary mechanism for driving the instability symptoms, and in cases where both the transient burning response in conjunction with vibration-induced acceleration is playing a role. The ability of the particles to suppress axial wave development is evident, at relatively low loading percentages, results that are consistent with experimental experience. This is clearly reflected by the respective attenuation maps. If this or a comparable numerical model proves to be suitably accurate, such maps could prove a useful tool for motor designers evaluating their own motor configurations for instability behavior.The adverse effect of loading a second (or a third) particle set that has a size that is less effective in suppressing pressure wave motion is also made evident in the present results. This study also gives some indication of what might be expected with reactive particles, where due to particle size reduction with time (under burning) the suppression effectiveness may not be quite what one would have ideally expected. In a similar vein, the present results are an indicator of the potential adverse effects of the agglomeration of particles in producing particle sizes bigger than what one would ideally find desirable for suppressing wave activity. These and other issues, some of them potentially quite complex, remain to be explored in regards to the use and modelling of reactive particles for suppression of combustion instability in SRMs.



---

*Source: 173129-2012-04-04.xml*