Multisized Inert Particle Loading for Solid Rocket Axial Combustion Instability Suppression

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.


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 lowmagnitude 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][2][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][5][6][7].On occasion, researchers have employed a numerical modeling approach, to work towards a more comprehensive quantitative understanding of the physics involved (the nu-merical 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][11][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 International Journal of Aerospace Engineering 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 statictest 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 d opt 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 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.

Numerical Model
A simplified schematic diagram of the physical system of an SRM, that is placed on a static test stand, is provided in Figure 1.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 statictest 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.

Equations Relevant to Two-Phase Flow in SRM.
One defines b 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] where ).The gas phase void fraction α g is defined by For a lower particle loading, one can assume that α g is close to unity in value (V p is the volume occupied by the particles in an elemental total volume of V ).The particle-gas mass flux ratio β g is stipulated as One can define an average two-phase density in a given volume as where α p is the particle-loading fraction in the flow, and N is the number of particles of average mass m p in the elemental volume above, then one can show the correlation between particle density ρ p and gas density ρ When one has two or more sets of particles of differing sizes (comprised of the same material; thus, the same solid specific heat C m ), one can use the following correlations: so that the nonequilibrium sound speed can be estimated via Under nominal flow conditions, u/b would be unity at the nozzle throat.The effect of particle mass loading ΔM p 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 The new effective solid propellant overall density becomes By substitution, one can show that

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 rate r b 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 Here, the total specific energy of the gas is defined for an ideal gas as E = p/[(γ−1)ρ] + u 2 /2.The corresponding equations of motion for an ith inert (nonburning) particle set within the axial flow may be found from International Journal of Aerospace Engineering Here, the total specific energy of a local grouping of particles from an ith set is given by E pi = C m T pi + u 2 pi /2, where T pi 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 D i , and the heat transfer from the core flow to a particle from the ith set is defined by Q i .In the case of drag between the gas and a representative spherical particle at a given axial location, one notes that where C d 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: where the Nusselt number Nu can be found as a function of Prandtl and relative Reynolds number for a sphere of mean diameter d mi .One will need to solve (13) for each of N set particle sets as part of the calculation process, where for dual or triple particle set loading, N set has a value of 2 or 3 in the present study.
Longitudinal acceleration a l 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, x positive moving right on structure as per Figure 1; acceleration of local surrounding structure rightward is designated positive a l ) 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.

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 finiteelement numerical solution, rather than via theoretical approximations [24].

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]: where r b,qs is the quasisteady burning rate (value for burning rate as estimated from steady-state information for a given set of local flow conditions), T i is the initial propellant 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 at this juncture the inclusion of a net surface heat release term, ΔH s , 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, Δy Fo = (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), r * b 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 r b may be found as a function of r * b through the empirical rate-limiting equation [25] The rate-limiting coefficient K b effectively damps the unconstrained burning rate r * b when for a finite time increment Δ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 r b with time) [26].As discussed in [25], the use of K b 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 rate r b,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: The pressure-based burning component may be found through de St. Robert's law The flow-based erosive burning component (negative and positive) is established through the following expression [27]: 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 r o , while at higher flow speeds, the positive erosive burning component r e , established from a convective heat feedback premise [27], should dominate: For the above case, where the base burning rate r o is a function of the other mechanisms (pressure and acceleration), one finds the velocity-based component of burn rate from (21) via r u = r b − r o .At higher flow speeds, r u becomes equivalent to r e .The effect of normal acceleration a n resulting from radial propellant/casing/sleeve vibration may be determined via [28] 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 a l are stipulated through the accelerative mass flux 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 a n relative to the burning propellant surface is present) [28].One should also note that a n 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 r o is a function of the other flow mechanisms (pressure and core flow), one finds the acceleration-based component of burning rate from (23) via r a = r b − r o .
With respect to the 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 for estimating the burning rate [25].While in the past a number of estimation models might have used a constant value for T s, more recently the usage of a variable T s has become prevalent.However, based on good comparisons in general to experimental data as reported in [25], a constant T s was employed in the present Z-Nbased phenomenological numerical combustion model, for the present investigation.

Results and Discussion
The characteristics of the reference motor for this study are listed in Table 1.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 M l , 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 (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 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 Figure 3 for head-end pressure p c 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 Δp d of 2 atm, the sustained compression wave front arriving about every 1 ms, oscillating at the fundamental frequency f 1L 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 (Δp w , 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.One can refer to Figure 4 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 [Δp w ], as compared to 1.42 MPa for the 0% loading case noted earlier [Δp w,peak ], giving a nondimensional attenuation M a , defined by International Journal of Aerospace Engineering 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.
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 M a 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 M a 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 M a Figure 10: Nondimensional attenuation as function of particle diameter and loading of a single particle set, reference motor, acceleration active.is around 0.52, while, for 40 μm particles at a 5% overall loading, the value for M a is around 0.55.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 Figure 8, allowing for the effect of vibrationinduced 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. 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 M a 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 M a is less than 0.15, using 10 μm particles at a 1.5% overall loading produces a value for M a is less than 0.2, and for 60 μm particles at a 1.5% overall loading, the value for M a 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.
One adjustment in particle size produces a significant change in the result of Figure 9, 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 M a 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.
Along the lines of Figures 7 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 (d m1 ), 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.
An alternative format for illustrating trends associated with using two sets of particles is provided in Figures 15, 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.

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 compositepropellant 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 vibrationinduced 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.

Nomenclature
A: local core cross-sectional area, m 2 a: gas sound speed, m/s a l : longitudinal (or lateral) acceleration, m/s 2 a n : normal acceleration, m/s 2 b: nonequilibrium sound speed of two-phase mixture C: de St. Robert coefficient, m/s-Pa n C m : particle specific heat, J/kg-K C p : gas specific heat, J/kg-K C s : specific heat, solid phase, J/kg-K D i : drag of gas on a particle from ith particle set, N d: local core hydraulic diameter, m d mi : mean particle diameter for ith particle set, m E: local total specific energy of gas in core flow, J/kg E pi : local total specific energy, ith particle set in flow, J/kg f : frequency, Hz, or Darcy-Weisbach friction factor G a : accelerative mass flux, kg/m 2 -s h: convectiveheattransfercoefficient, W/m 2 -K ΔH s : net surface heat of reaction, J/kg K b : burn rate limiting coefficient, s −1 k: gas thermal conductivity, W/m-K k s : thermal conductivity, solid phase, W/m-K M a : magnitude of attenuation M : limit magnitude, cyclic input m pi : mean mass of a particle from ith particle set, kg N i : number of particles from the ith set in a given volume N set : total number of particle sets n: exponent, de St. Robert's law p: local gas static pressure, Pa Δp d : initial pulse disturbance step pressure, Pa Q i : heat transfer from gas to a particle from ith particle set, W R: specific gas constant, J/kg-K r b : instantaneous burning rate, m/s r b,o : reference burning rate, m/s r b,qs : quasisteady burning rate, m/s r * b : unconstrained 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 pi : temperature of particle from ith set, K T s : burning surface temperature, K Δt: time increment, s u: core axial gas velocity, m/s u pi : core axial particle velocity for particle from ith set, m/s v f : nominal flamefront velocity, m/s V : elemental volume, m 3 x: distance from head end, m Δx: spatial increment in axial direction, m y: radial distance from burning surface, m Δy: spatial increment in radial direction, solid phase, m Δy Fo : Fourier limit spatial increment, m α g : gas phase void fraction α p : total particle mass fraction of overall core flow α pi : particle mass fraction of ith particle set in overall core flow α s : thermal diffusivity, solid phase, m 2 /s β g : particle-gas mass flux ratio γ: gas ratio of specific heats δ o : reference combustion zone thickness, m δ r : resultant combustion zone thickness, m k: vibration-based wall dilatation term (1/A • δA/δt), s −1 μ: absolute gas viscosity, kg/m-s ρ: gas density, kg/m 3 ρ p : total density of particles in core flow, kg/m 3 ρ pi : density, ith particle set in core flow, kg/m 3 ρ m : solid density of particle, kg/m 3 ρ s : solid density of propellant, kg/m 3 ρ so : reference solid density of propellant (no particle loading), kg/m 3 φ: acceleration orientation angle, rad φ d : longitudinal/lateral-acceleration-based displacement orientation angle, rad.

Figure 7 :Figure 8 :
Figure 7: Nondimensional attenuation as function of particle diameter and loading of a single particle set, reference motor, acceleration nullified.

Figure 15 :
Figure15: Nondimensional attenuation as function of particle diameter (d m2 : d m1 ), for a 1 : 1 loading distribution for two sets of particles at total loading of 3%, reference motor, acceleration nullified.