Modeling of nonlinear longitudinal instability in solid rocket motors

Acta Astronautica Vol. 13, No. 6/7, pp. 339-348, 1986

0094-5767/86 $3.00+ 0.00 © 1986 Pergamon Journals Ltd

Printed in Great Britain. All rights reserved

MODELING OF NONLINEAR LONGITUDINAL INSTABILITY IN SOLID ROCKET MOTORSt JOSEPH D. BAUM~ and JAY N. LEVINE Air Force Rocket Propulsion Laboratory, Edwards AFB, CA 93523-5000, U.S.A. (Received 28 January 1986)

A~traet--A comprehensive model of nonlinear longitudinal combustion instability in solid rocket motors has been developed. The two primary elements of this stability analysis are a finite difference solution of the two phase flow in the combustion chamber and a coupled solution of the nonlinear transient propellant burning rate. A new combination finite difference scheme gives the analysis the ability to treat the type of multiple travelling shock wave instabilities that are frequently observed in reduced smoke tactical solid rocket motors. Models for predicting the behavior of both gas ejection and solid ejecta pulses were developed and incorporated into the analysis. Extensive comparisons between model predictions and experimental data from pulsed solid rocket motor firings have been carried out. The nonlinear instability analysis was found to be capable of predicting the complete range of nonlinear behavior observed in actual motor firing data. Good agreement between measured and predicted initial pulse amplitude, pulse evolution, limit cycle amplitude and mean pressure shift was obtained. This investigation has also provided new insight into the nature of nonlinear pulse triggered instability and the factors which influence its occurrence and severity. This new instability analysis should significantly enhance our capability to design tactical solid rocket motors that are free from troublesome and expensive nonlinear combusion instability problems.

1. INTRODUCTION Tactical solid rocket motors are frequently subject to a combustion instability problem at some point in the design cycle. When instability is encountered it can take one of several forms, e.g. linear or nonlinear, longitudinal or tangential. Over the last twenty years, considerable resources have been expanded to understand, predict, control and eliminate combustion instability in solid rocket motors. Most of this effort has been devoted to linear instability problems, and as a result, such problems can now be treated in a rational, cost effective manner. Comparatively little work has been accomplished towards the understanding and resolution of nonlinear combusion instability problems. Thus, when nonlinear instabilities are encountered, the solution is too often an expensive cut and try process. Linear instabilities are characterized by small amplitude, sinusoidal oscillations that originate from the amplification of infinitesimal random disturbances in the motor chamber. The growth or decay of random small amplitude pressure oscillations is determined by a delicate balance between the driving and damping of oscillatory energy. The processes which dissipate acoustic energy include convection and radiation of acoustic energy through the nozzle (nozzle damping), the viscous and thermal losses produced by interaction of the condensed phase combustion products tPaper presented at the 36th Congress of the International Federation, Stockholm, Sweden, 7-12 October 1985. :~Currently with the Laboratory for Computational Physics, Naval Research Laboratory, Washington, D.C. 20375, U.S.A. Astronautical

with the combustion gases (particle dampings), the nonlinear viscoelastic characteristics of the grain case (structural damping), viscous losses in the gas phase and at inert surfaces, gas phase vibrational relaxation effects, radiation of energy through the motor case and inelastic acceleration of combustion products leaving the propellant surface (flow turning). The primary source of oscillatory energy considered to date has been the response of the propellant combustion zone to acoustic pressure and acoustic velocity oscillations (termed pressure and velocity coupling, respectively). Other sources of energy include the combustion of incompletely reacted products in the chamber, distributed combustion of metal droplets in the gas phase, and conversion of mean flow energy to acoustic energy due to vortex shedding. The determination (both by theoretical predictions and experimental measurements) of the various energy sources and sinks in a combustion chamber is critical to our ability to a priori predict the stability characteristics of proposed solid rocket designs. Nonlinear instabilities are usually characterized by large amplitude oscillations having steep-fronted, shock-like waveforms. Nonlinear axial mode instability in solid propellant rocket motors is initiated by random finite amplitude evens such as the expulsion of an igniter or insulation fragment through the nozzle. When instability is initiated in this m a n n e r in a motor that is otherwise linearly stable (i.e. stable to infinitesimal disturbances) it is said to be a "triggered" instability. Nonlinear behavior can also become important when "linear" instabilities reach an amplitude sufficient to cause a transition in the nature of the oscillations. The existence of triggered in-

339

340

JOSEPH D. BAUMand JAY N. LEVINE

stabilities is a direct result of the fact that all of the aforementioned acoustic energy gain or loss mechanisms are nonlinear, i.e. amplitude dependent to some degree. These same nonlinearities also ensure that a nonlinear instability will not grow without limit, but rather will eventually reach a limit cycle amplitude at which the net gains and losses are balanced. Nonlinear axial mode instabilities usually result in pressure oscillations that propagate as steep-fronted waves which are actually weak shock waves. The acoustic pressure and velocity oscillations are frequently accompanied by an increase in mean chamber pressure (usually referred to as a DC shift) and increased mean propellant burn rate. These oscillations and pressure shifts can alter the thrust time characteristics of the motor, affect the operation of sensitive electronic equipment, and, in extreme cases, result in motor case failure. Certain trends and characteristics of nonlinear instability have been documented. However, attempts to form generally applicable conclusions have been stymied by the number, complexity of. and mutual interactions between the governing physical phenomena. The objectives of this research work were to further our understanding of the physical mechanisms which control the initiation and severity of nonlinear combusion instability and the development of a comprehensive analytical model capable of predicting the relative nonlinear stability of proposed solid rocket motor designs. In flight, nonlinear instability is triggered by random events. To facilitate the timely conduct of research and motor development programs, artificial, but controllable, laboratory pulses are used to study triggered instabilities. Thus, in order to enhance the value of the analysis to solid rocket motor designers, an additional objective of this work was the development of models to predict the performance of laboratory pulsers, and the motor response to pulsing. 2. APPROACH

Efforts to understand and model nonlinear instability date back to the 1960s[1-3]. The most recent work has been divided between so-called "exact" and "approximate" mathematical approaches. The "exact" methods of Levine and Culick[4] and Kooker and Zinn[5] numerically solve the nonlinear partial differential equations governing both the mean and time dependent flow in the combustion chamber, as well as the combustion response of the solid propellant. The "approximate" methods of Culick[6] and Poweil et al.[7], utilize expansion techniques to reduce the problem to solution of sets or ordinary differential equations. All of these analyses have deficiencies which preclude them from treating the multiple shock, triggered type of instabilities that tNomenclature is given in Appendix at end of paper.

occur in modern variable cross-section area reducedsmoke motors. In this investigation, the most advanced model, (Levine and Culick[4]) was modified and extended to overcome its deficiencies. The new developments include incorporation of: an advanced shock capturing finite difference scheme; the ability to treat velocity coupling effects; and models for predicting the behavior of laboratory pulse units. The two primary elements of the nonlinear stability analysis are a finite difference solution of the two phase flow in the combustion chamber, and a coupled solution of the nonlinear transient propellant burning rate. The developed models and the solution techniques are described in the following sections. 3. MODEL DEVELOPMENT 3.1 Equations o f motion The following assumptions were employed in formulating the two phase equations of motion: the gas in the motor is ideal, nonreacting and inviscid (except for particle-gas interactions) with constant specific heats, constant Prandtl number, etc.: the flow is quasi-one dimensional; since the cycle to cycle area variation is small, the time derivative of the local cross-sectional area can be neglected: the metal particles are completely oxidized at the burning surface, the resulting inert particles are spheres, divided into a number of discrete size groups; the weight fraction of each size group is invariant with location; the particles have uniform internal temperatures and occupy negligible volume; the different particle size groups do not interact; all particles of a given size, at a specific location, have the same velocity and temperature (this implies that the entering particles instantaneously exchange heat and momentum to acquire the local particle velocity and temperature); the gas and particles from the burning propellant enter normal to the burning surface, at the local instantaneous flame temperature of the propellant; and the effects of propellant surface inclination and nonisentropic flame temperature variations are accounted for. Subject to these assumptions, the equations of motion in conservative form aret Continuio Gas phase c~p + I 8(puA) -- (~) & A 8x

(l)

Particles ) (pp,) + 1 O(pp,upA A 8x 8t

-

fl~

(2)

Momentum Gas phase 8(pu) 8t

1 8 PdA A ~ [(e + p u 2 ) A l = A d x + u ' °

+ ~ [G + A~)(u.,--.~,)l

(3)

Instability in solid rocket motors Particles

1 a(ppiu2A) ~(pp,up,) + ~ ' = -Fp. + fl, o~up, (4) ~t A ~x Energy Gas phase

341

The term in parenthesis is the Nusselt number correlation suggested by Carlson[9]. To facilitate the shock capturing finite difference solution of the gas phase equations, the conservative form of the equations was employed. The particle equations were, however, formulated and solved in nonconservative form. 3.2. Boundary conditions

u2

1 8(uPA) N

+o~(cpTf+ ~IU2, ) + Z {Qp, + UpFp, i=l

+/~:[cm(:r

~,,)

'

~

~ )]}

(5)

Particles 0 ~,(cmrp ,

1 0 ,2 .___ , + ~uu,,)] + AOx

XLOpUp.(c.Tp. I 2 A ]=-Qp-up:F,, , , ,+~Up,) +~:o(c,.T,. , + ~up,) 12

(6)

The equations of motion are supplemented by the following equations State

P = pRT

Steadystateburnrate

-r = ~ ( ~ "~"(1 +CkU) (7) \ Pref /

The quantity co is related to the burn rate by the following relation 09 =

p~rl A(1 + fl)cos 0

(9)

where I is the local perimeter of the propellant grain. Reynolds number

Rei = p [up, - u [Di,

(10)

where p, the gas viscosity, is given by

"-=(5) °" #o

(ll)

\To]

The quantity Fp, that represents the drag force exerted by the ith particle group on the gas, per unit volume, is defined as / 18pp.\ Fp., = t ~ ) l a ( U p , _

/'I Re~/3"~ - u ) t +--g---)

(12)

The drag coefficient, ca, is defined as follows cd=~e \

+

(13)

This definition was obtained from a correlation of the "standard" drag coefficient vs Reynolds number data for a single sphere in steady flow[8]. The heat transfer term Qp, that defined as Qp, :

f

6cp "~

x (re, - T) (2 -4- 0,459Re°55Pr °33)

04)

Only a single boundary condition is required when the head end of the motor is inert, i.e. u = 0, up~ = O. If a motor has propellant burning at the fore-end, both the gas velocity and temperature at the edge of the flame zone must be specified. It is assumed that the particles and the gas are in thermal and dynamic equilibrium, thus Up, = u and Tp, = T. The particle density Ppi, is given by the particle to gas weight of the propellant. The most accurate treatment of the nozzle e n d - integration of the equations through the throat into the supersonic portion of the nozzle--does not require the specification of any downstream boundary conditions. For the sake of computational efficiency, however, an option to terminate the solution at the nozzle entrance was also developed. This option utilizes the concept of an equivalent perfect gas and the quasi-steady short nozzle approximation to define a Mach number boundary condition at the nozzle entrance. For linear waves, it has been shown[10] that this is equivalent to specifying a real nozzle admittance equal to (7 - 1)M/2. 3.3. Nonlinear transient burn rate model Due to existing deficiencies in the understanding of velocity coupling mechanisms, the present approach was to develop a complete nonlinear pressure coupled response model and simplified (ad hoc) velocity coupled models. Past investigations of nonlinear burning rate response have been carried out utilizing both expansion procedures[11] and numerical schemes[12]. The expansion schemes are too cumbersome and are difficult to extend to higher orders. Thus a numerical approach was selected for use in this investigation. The model utilized represents a nonlinear extension of the Denison and Baum model[13]. The assumptions employed in the current model are: the gas and solid phase are treated as homogeneous; only variations in the direction normal to the burning surface are accounted for; the solid is converted to gas at an infinitesimally thin surface; subsurface reactions are not considered; the gas phase responds in a quasi-steady manner; the flame is anchored at the surface; heat release is uniformly distributed from the surface to the edge of the flame; the gas reaction rate is a function of pressure only; the presence of metal in the propellant is not directly accounted for; gas and propellant specific heats and thermal conductivities are taken to be constant, but different; and the surface pyrolysis rate is given by an

342

JOSEPH D. Baum and JAY N. LEVINE

Arrhenius law modified to account for pressure coupled surface reactions. With these approximations the calculation of the transient burning rate reduces to the solution of the following nonlinear heat conduction equation, written in a coordinate system fixed to the burning surface FsT

~T

~2T

--+r~-=~ gt cx

15)

~?_r:

The boundary condition at the back wall is T~...... , = T,

16)

The boundary condition at the burn surface that incorporates pressure coupling only (velocity coupling effects will be treated later in this paper) is obtained by performing a mass and energy balance for an infinitesimally thin control volume at the gas/solid interface. The mass balance yields the relation p~r = pgUg while the energy balance can be shown[4] to yield the following boundary condition:

\

=

p,a

+ T -- T~ - 0-2

+

[0~. + ( c , - cp)(T~- T,.)]

The nonlinear heat conduction equation subject to the boundary conditions at the backwall and at the burning surface is solved using the Crank-Nicolson finite difference scheme[14]. This scheme has temporal and spatial second order accuracy. Several other schemes, including the fourth order (spatial) accuracy Keller-Box scheme[15] and compact differencing scheme[16] were tried. Nevertheless for the desired accuracy (_+ 1%), the Crank-Nicolson scheme requires less computational time. 3.4. Shock capturing.finite diflerence scheme The finite difference scheme of Rubin and Burstein[17] was utilized originally[4] to integrate the two phase equations of motion. However, as shown in Fig. la, when this method was utilized to solve a traveling shock wave type of instability (i.e. propagation of traveling shock waves in the combustion chamber), the solutions obtained demonstrated that post-shock oscillations developed in time into a number of discrete "'humps.'" It has been demon-

(o)

0.30 ,

0.20

g co

O.lO

P,

--'

o.oo E z

-O.tO --~

-0.20

I 0

I

10

1

I

J

3O

40

5O

1 40

1 50

]

Nondimensional time

0.20

(b)

o.Io

i

0.00

i

-O.lO -0.20 -0.30 0

t 10

I zo

i 30 Nondimensional 'time

7 60

Fig. 1. Time evolution of normalized pressure oscillations at an end of the chamber. (a) Rubin-Burstein, (b) L W + H + ACM.

Instability in solid rocket motors strated[18] that these wiggles result from erroneous energy transfer between modes (i.e. dispersive errors that cause pressure signals to travel at the wrong speed) and from dissipative errors that result in excessive damping of energy contained in the high frequency modes. In order to remedy this problem the capabilities of several recently developed shock capturing finite difference schemes were evaluated[19] by utilizing them to solve several test problems. The best results were obtained with a combination of the Lax-Wendroff scheme[20], the Hybrid scheme[21], and the Artificial Compression scheme[22] (termed here: L W + H + ACM). This combination method involves two steps: in the first step the second order Lax-Wendroff scheme is hybridized with a nonoscillatory first order accurate method to allow a monotonic (i.e. nonoscillatory) transition across discontinuities; in the second step, an artificial compression correction is applied to sharpen transitions of discontinuities, since the hybridized first order accurate method is too dissipative. A switch value based upon flow gradients (density gradients were used herein) is utilized so that the artificial compression and the first order monotonic schemes are activated only in the immediate vicinity of admissible discontinuities. This combined method preserves the second order truncation error of the Lax-Wendroff scheme (in the smooth parts of the solution) and yet, as shown in Fig. lb, yields sharp and oscillation-free shock transitions. The L W + H + A C M finite difference integration scheme is utilized only at interior mesh points. Numerical stability and uniqueness problems are avoided by using the method of characteristics to obtain solutions at the boundary points and at interior geometrical and combustion discontinuities (e.g. partial length grains, uninhibited grain ends, etc.). 3.5. Velocity coupling model The results of a J A N N A F workshop on velocity coupling[23] indicated that existing velocity coupling models have significant deficiencies and are not based on an understanding of the fundamental physical mechanisms. The objective of this research work was not to develop a new velocity coupling model but

343

to demonstrate the capability of the nonlinear combustion instability model to predict triggering phenomena. Thus it was decided to treat velocity coupling by utilizing several different ad hoc functions of velocity to augment either the heat transfer to the propellant surface or the transient burn rate itself. The heat transfer augmentation approach yielded significant nonlinear effects (DC pressure shifts and modulated limit cycle amplitude) only at unrealistically high values of response function (Rye= 60, where Rvc is the velocity coupled response function), as shown in[24]. In contrast, it was shown[24] that the burn rate augmentation model yields strong nonlinear effects with velocity coupled response functions as low as 3.2. The relative ineffectiveness of the heat transfer augmentation model was shown to be a result of the characteristics of the response function implied by the Denison-Baum type model utilized. It was also concluded that a realistic velocity coupling model has to provide strong driving at the higher harmonics in order to predict such phenomena as triggering, mean pressure shift, modulating limit cycle, etc. The burn rate augmentation model is implemented by simply modifying the instantaneous mass burning rate (including pressure coupled effects) to account for the presence of acoustic velocity fluctuations as follows Wpc = Woe[1 + RwF(u)]

where W is the instantaneous propellant mass burning rate (W = I,P + W') and W~ is the instantaneous mass burn rate computed from the pressure coupled model only. Several forms for F(u) were tested. Results presented herein were obtained by utilizing F ( u ) = l u ' [ . In the limit, as u' approaches zero, equation (18) reduces to the linear velocity coupling model used in[4]. 3.6. Pulser modeling Several types of pulser units have been employed in previous nonlinear instability studies. These pulsers can be categorized as either gas injection pulsers, e.g. pyrotechnic and piston pulsers, or solid ejecta devices. A schematic of the pyro pulser is shown in Fig. 2. A squib initiator ignites a pyrotechnic charge (Red Dot double-base powder) in the combustion chamber. The gaseous combustion prod-

Shield

o Motorod

:;ooo0,er--\ ~.~.~.

--Pyro 0u,serhoosin "~ ~ .° ~

Orifice plate--/' / . . . . Burst diophragm_ /

(18)

Pyrotechniccharge Holex-1,196squib

~ ~ "--Pressure transducer \~ - \~" Polyurethonefoam restrainer ~ Votumespacer

Fig. 2. Schematic of the pyro pulser unit.

344

JOSEPH D. BAUMand JAY N. LEVINE

leeve

7/

~

~lnitiotor squib

~

~ / ~

0 ring seals Piston Teflon stopper

~.-Pre,ssure transducer ~ ' - " Burst diapnrogm "-- Styrofoam retainer Pyrotechnic charge

Fig. 3. Schematic of the piston pulser unit.

calculated using these pulser performance models were utilized as boundary conditions for the chamber flow problem, which was solved by modifying the nonlinear combustion instability program[25]. A combined experimental and analytical study resulted in the development of two methods for predicting initial nozzle ejecta induced pulse amplitudes[26]: a simple model based on linear wave propagation theory and the assumption of quasi-steady nozzle behavior; and a numerical model, which utilizes the quasi-steady assumption to provide a nozzle entrance boundary condition to the nonlinear instability program. 4.

ucts increase the pressure in the chamber to the rupture pressure of the burst diaphragm. As the diaphragm bursts, the combusion products plus a fraction of the remaining unburned pyrotechnic charge expand into the pulser barrel and are then vented into the chamber. In order to allow for more variation in the range of pulse characteristics available for motor test and research purposes, another type of pulser called the piston pulser was designed. A schematic of the piston pulser is shown in Fig. 3. After rupture of the burst diaphragm, the combustion gases and the unburned fraction of the charge enter the breech volume. As the pressure increases, the piston accelerates into the bore volume. The gases in the bore are compressed by the piston and are vented in the motor. A schematic of an ejecta pulse test apparatus is shown in Fig. 4. At predetermined times during the firing, the pulser ejects a sphere into the chamber which is carried out through the nozzle by the flow of combustion gases. To simulate the range of material densities found in combustion chambers, spheres made from RTV, rubber, nylon, teflon and steel are utilized. Ballistic models of the pyrotechnic and piston pulsers were developed utilizing a simple lumped volume treatment. The mass and energy flow rates

RESULTS

A large number of nonlinear instability solutions have been obtained using the model outlined above. In most cases these solutions have been compared to experimental data obtained in both laboratory scale devices and full scale motors. These solutions have provided insight into the nature of nonlinear instability, and have served to demonstrate the validity of the modeling. Highlights of the results are summarized in this section. A more complete description of the results obtained to date may be found in[27],

4.1. Limiting amplitude studies Figure 5a shows the pressure time history at the head end of a cylindrically perforated (C-P) motor (chamber pressure of 18.03 MPa/2616psi) initiated by perturbing the steady state with a fundamental mode disturbance with an amplitude of 40% of the mean pressure. The wave reached a limit cycle amplitude of 21.73% of the mean pressure (peak to peak). The calculation was performed for a propellant without particles and a linear pressure coupled response function of 3.2. Other solutions were obtained with initial disturbance amplitudes of 60%, 20%, 10% and 2% of the mean pressure. All of these solutions reached the same limit cycle (i.e. same limiting amplitude and waveform). The 2% initial disturbance solution is shown in Fig. 5b. Several other series of Transducer locationcm . in ) 45

Fwd closure

I

"

(48.0) 30.5 (

'I I

q

2o)

{~ JK 4 /

T]'

( 610 ) ""1 K5 /

/ K6

ir

/ tO pulser K~

:

Fig. 4. Schematic of the ejecta pulse test apparatus.

o

Instability in solid rocket motors

345 b)

0.30.

0.45]

(a I

O.lO]

._~ 0.20.

g

~ o.,io.

~ o.o5

OlOO

o.oo,

¢_ -o2o

~_ -040.

-0.30

-0,45-

-0.40

-0.20

0.0

flOiO

200

300

400

50'.0

600

O0

200

40.0

600

800

1000

Nondimensionol time

Nondimensionol time

Fig. 5. Time evolution of normalized pressure oscillations at the head-end of the motor, first harmonic initial disturbance. (a) P~ = 0.4 ?, (b) P6 = 0.02 _P.

calculations were performed with other motor/ propellant combinations and with varying sizes and amounts of particles. The same limit cycle is obtained for a defined size/amount of particles, independent of the amplitude or characteristics of the initial disturbance. Thus it has been demonstrated that pressure oscillations in unstable rocket motors (with metallized as well as unmetallized solid propellants) reach a limit cycle that is independent of the characteristics of the initiating disturbance. Extensive experimental results obtained in a recently completed Air Force Rocket Propulsion Laboratory program with Aerojet Tactical Systems[28] support the analytical conclusion regarding limit cycle uniqueness. The only exception being that, in some cases, when a motor was pulsed after it reached limit cycle oscillations[14] a second traveling shock was developed and was supported along with the original wave. 4.2. Geometry effects In order to meet difficult mission requirements, tactical motors often have relatively complicated grain geometries that contain one or more discontinuous jumps in cross-sectional area. In the presence (0 )

0

of traveling wave type instabilities such discontinuities result in internal reflections that greatly complicate the observed waveforms. No other existing analysis has the capability to treat such actual motor geometries. The ability of the present analysis to treat such geometries is demonstrated below. The first example shown in Fig. 6 is for a small partial length grain laboratory motor. The motor was pulsed during the firing. The very sharp shocked initial waveform decays rapidly to an almost sinusoidal waveform indicating that at this time the motor is stable. Figure 6 shows a comparison of the measured and predicted pressure perturbations at the fore-end of the motor. Excellent agreement is demonstrated between the measured and predicted initial pulse amplitude and waveform, the decay rate of the pulse, and the temporal evolution of the waveform. In particular, notice the ability of the analysis to reproduce the generation of a second shock due to partial reflection from the area discontinuity at the end of the grain. The second example is for a motor with a so-called finocyl geometry that created a sizable grain area discontinuity approximately 2/3 of the distance from the fore-end to the nozzle. As shown in Fig. 7, this motor was pulse triggered into

46 psi per division

0.01

002

003

0.04

Time, s

o.o8

"

(b) I

,

p, oo43 I [IJJAII

L,

,AAA

VVvvvvvvvvvvv

0.00

o

1"o

Eo

~o

Nondimensional time

Fig. 6. (a), (b) Comparison of measured and predicted pressure perturbations at the fore-end, test PCC3, first pulse, 6 = 10.3 MPa (1500 psi).

346

JOSEPH D. BAUM and JAY N. LEV1NE (a )

-~0.

P

96 psi Per division

(bi d92 psi Predicted

24 ps, per dwiston

per

division ~2 r~s ~er divsion

0-

0

!

I

i

OOdO

l

0020

'1

ol3

O. 0

i

' "1

0040 O0

Time,s

Fig. 7. (a), (b) Comparison of measured and predicted pulse induced pressure oscillations at the fore-end, test PSMA-! 1, pulse 2. nonlinear instability. Reflection at the grain area discontinuity immediately caused a second wave to appear. This wave was amplified until the resulting instability contained a dual shock wave system. The numerical results faithfully reproduced this complex wave development. 4.3. Triggering Before this work was initiated, there was much speculation about the cause of triggered instability, but very little real information. No analysis was capable of predicting triggering, and certainly not with all of the attendant complex wave behavior that is frequently observed in motors. Results obtained with the present model have helped to provide new insight into the physical mechansms that produce triggering. As shown in Fig. 8, it was found that nonlinear gas dynamic effects, the effects that produce the steepening of finite amplitude waveforms into shock waves, result in the convection of energy from the fundamental acoustic mode of the chamber to higher harmonics. This is typically a wave dampening process, since the higher harmonics are usually considerably more stable than the fundamental mode. Nevertheless, it was found that in the presence of sufficiently high harmonic amplitudes, nonlinear combustion processes associated with the effect of acoustic velocity-propellant combustion interactions

I Fundomentol mode frequency Nonlineor fluid dynomic ~'oce~s

.,

Nonlinear combustion processes I

Higher harmonics

Fig. 8. Schematic of triggering process.

I

I

0 o2

I

|

O0 ~. Time~s

I

I

o(36

~

' -"~ © 08

Fig. 9. Comparison of measured and predicted pulse induced pressure oscillations at the head-end, test 7. second pulse. can overcome this damping by providing more energy to the high harmonics and by transferring energy from the higher harmonics back to the fundamental mode. This two element triggering scenario is vividly depicted in Fig. 9. This figure shows the measured and predicted head-end pressure oscillations in a cylindrically perforated solid rocket motor that was pulsed by ejecting a steel ball out of the nozzle. Being heavy, the steel ball moved through the nozzle relatively slowly, producing a pulse that is almost sinusoidal at the fundamental mode frequency. The gasdynamically produced transformation to a steep fronted wave takes about nine wave cycles, during which time the amplitude continuously decays. Only when the wave has steepened and contains sufficient energy at the high harmonics, does it begin to grow, as a result of nonlinear combustion feedback to the fundamental mode. The ability of the analysis to correctly reproduce this complex behavior is particularly gratifying. Knowledge of the multimodel character of the triggering process is also of significant practical importance. The energy that sustains an instability is derived from the interaction of the oscillations with the propellant combustion process. This interaction is highly frequency dependent. Prior to this work efforts to stabilize a motor exhibiting fundamental mode instability centered around reducing the driving of the propellant at the fundamental mode frequency. While this often provided the solution to linear instability problems, more often than not, it did not result in the elimination of nonlinear triggered instability. The present results point out the need to consider the propellant response for at least the first several harmonics. In support of this conclusion it should be mentioned that a nonlinear instability problem in a reduced smoke motor was solved by changing the oxidizer from a bimodal to a trimodal

Instability in solid rocket motors •

L

.

~

.

i

i

8 E

000

005

040

0.45

020

Time,s ~2o~

347

strate that the model consistently predicted amplitudes within 10% of the measured pulse amplitudes, on an a priori basis. Thus the present analysis provides the needed ability to design pulsers that will yield the pulse amplitudes and waveforms needed to verify the nonlinear stability of motors during preproduction testing.

(b) 5.

4

oo

•

.

.

•

.

.

4o0

2oo

300

400

5oo

6oo

Nondimensional time

Fig. 10. (a), (b) Pyro-pulser. Comparison of measured and predicted pulse induced pressure perturbations at the headend of the chamber. blend. The measured combustion driving of the two propellants was almost identical at the fundamental mode frequency, however, the trimodal blend had significantly lower driving at the first several harmonic frequencies. 4.4. Pulser modeling As previously mentioned, one of the objectives of the present work was to produce models of all the types of pulsers currently used to study nonlinear instability. The objective was to develop a capability to predict the pulse amplitude, pulse waveform and motor response to pyrotechnic, piston and ejecta pulsers. Results cited in the previous two sections demonstrate that the goal was achieved for piston (Fig. 6) and ejecta pulsers (Figs 7 and 9). Figure 10 shows a comparison of predicted and measured headend pressure oscillations for a stable pyrotechnic pulse. Here again excellent agreement was achieved. The present model was used to predict the pulse amplitude for a series of laboratory scale pulsed motor firings[29] involving both pyrotechnic and ejecta pulsers. The results shown in Table 1 demonTable 1. Comparison of measured and predicted pulse amplitudes (in percentage of mean pressure) at the head end closure Head End Test Number 1 first pulse 2 first pulse 3 first pulse 4 first pulse 5 first pulse 5 second pulse 6 first pulse 6 second pulse 7 first pulse 7 second pulse 8 first pulse 9 first pulse

Measured

Predicted

3.4 3.6 1.7 1.6 3.9 9.6 4.4 9.7 5.9 9.2 6.6 19.9

3.5 3.7 1.8 1.6 3.8 10.8 4.7 10.2 6.1 10.3 7.8 22.7

C O N C L U S I O N S

Nonlinear combustion instability has been a recurrent problem during the development of solid rocket motors. The current trend towards lower burn rate, long L/D motors serves to further exacerbate this problem. In response to the need for tools to help motor designers cost effectively prevent or eliminate such problems, a new nonlinear combustion instability model has been developed. This model describes the two phase flow field in the chamber and its transient interactions with the burning propellant surface. The nonlinear instability analysis was found to be capable of predicting the complete range of nonlinear behavior observed in motors. Good agreement between measured and predicted initial pulse amplitudes, waVeform evolution, growth and decay rates and instability induced mean pressure shifts was obtained. An advanced shock capturing finite difference scheme based on a combination of the Lax-Wendroff, Hybrid and Artificial Compression schemes was developed. This scheme provided the capability to treat the multiple shock, triggered type of instability that occurs in modern variable crosssectional area reduced smoke motors. An ad hoc velocity coupling model was developed, which, when incorporated into the analysis, provided, for the first time, the capability to predict triggering and the associated shifts in mean operating pressure. Pulser performance models were developed for all of the gas ejection and solid ejecta pulsers currently used to study or test for nonlinear instability. These pulser models provide an ability to design pulsers that produce the pulse amplitudes and waveforms desired to verify the nonlinear stability of motors during motor development testing. In addition to the aforementioned practical significance of the present work, it has also provided insight into the nature of nonlinear pulse triggered instability, and the factors which influence its occurrence and severity. In doing so, the analysis has identified new approaches for solving or preventing future nonlinear instability problems. REFERENCES

1. E. W. Price, Axial mode intermediate frequency combustion instability in solid propellant rocket motors. AIAA Paper 64~146, January (1964)• 2. W. G. Brownlee, Nonlinear axial combustion instability in solid propellant motors. AIAA J. 2, 205-284 (1964). 3. G. A. Marxman and E. E. Wooldrige, Finite-amplitude

348

4.

5.

6.

7.

8.

9. 10. 1l. 12. 13. 14.

15. 16. 17. 18.

19.

20. 21. 22.

23.

24.

25.

26. 27

28.

JOSEPH D. BAUM and JAY N. LEVINE axial instability in solid rocket combustion. Twelfth Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, Pa, pp. 115 127 (1969). J. N. Levine and F. E. C. Culick, Nonlinear analysis of solid rocket combustion instability. Air Force Rocket Propulsion Laboratory, Edwards AFB, Calif. AFRPLTR-74-45, October (1974). D. E. Kooker and B. T. Zinn, Numerical investigation of nonlinear axial instabilities in solid rocket motors. BRL-CR-141, March (1974). F. E. C. Culick, Nonlinear behavior of acoustic waves in combustion chambers. Proceedings of the 10th JANN A F Combustion Meeting, Vol. I, CPIA Publication 243, pp. 417436 (1973). E A. Powell, M. S. Padmanabahn and B. T. Zinn, Approximate nonlinear analysis of solid rocket motors and T-burners: Air Force Rocket Propulsion Laboratory, Edwards AFB, Calif., AFRPL-TR-77-48 (1977). G. Rudinger, Effective drag coefficient for gas-particle flow in shock tubes. Project Squid Technical Report CAL-97-PU, March (1969). D. J. Carlson and R. F. Hoglund, AIAA J. 2, 1980, (1964). B. T. Zinn, J. Sound Vibration 22, 93 (1972). R. H. Brown and R. J. Muzzy, A I A A J. 8, 1492 (1970). H. Krier et al.,: AIAA J. 6, 278 (1968). M. R. Denison and E. Baum, A R S J. 31, 112 (1968). J. Crank and P, Nicholson, A practical method for numerical integration of solutions of partial differential equations of heat-conduction type. Proc. Cambridge Philos. Soc. 43, 50 (1947). H. B. Keller, Annual Review qflquid Mechanics (Milton Van Dyke, Ed.) Annual Reviews Inc. (1978). Y. Adam, J. o f Comp. Phys. 24, 10 0977). E. L. Rubin and S. Z. Burstein, J. o[Comp. Phys. 2, 178 (1967). J. D. Baum and J. N. Levine, A critical study of numerical methods for the solution of nonlinear hyperbolic equations for resonance systems. J. Comput. Phys. 58, 1 (1985). J. D. Baum and J. N. Levine, Numerical techniques for solving nonlinear instability problems in solid rocket motors. A1AA J. 20, 955-961 (1982), P. D. Lax and B. Wendroff, Commun. Pure Appl. Math. 13, 217 (1960). A. Harten and G. Zwas, J. Comput. Phys. 9, 569 (1972). A. Harten, The artificial compression method for computation of shock and contact discontinuities: III, self adjusting hybrid schemes. AFOSR TR-77-0659, March (1977). M. W. Beckstead, Report on the workshop on velocity coupling. Proceedings of the 17th J A N N A F Combustion Meeting, Langley, Va. CPIA Publication 324, November, pp. 195 210 (1980). J. N. Levine and J. D. Baum, A numerical study of nonlinear instability phenomena in solid rocket motors. AIAA J. 20, 557 564 (1982). L D. Baum, J. N. Levine and R. L. Lovine, Pulsing techniques for solid propellant motors; modeling and cold flow testing. J. Spacecraj? Rockets 20(2), 150- 157 (1983). R. L. Lovine, J. D. Baum and J. N. Levine, Ejecta pulsing of subscale solid rocket motors. AIAA ,1. 23, 416-423 (1985). .I.D. Baum and J. N. Levine, Modeling of nonlinear combustion instability in solid propellant rocket motors. AFRPL-TR-83-058, Air Force Rocket Propulsion Laboratory, Edwards AFB, Calif. February (i984). R. L. Lovine, Nonlinear stability of tactial motors: Vol. V: Propellant screening methodology. AFRPLTR-84017, Sept. 1985.

29. J. D. Baum, J. N. Levine and R. L. Lovine, Pulse triggered instability in solid rocket motors. AIAA `1. 22, 1413-1419 (1984). APPENDIX NOMENCLATURE A a B, c ? Q c, c., c; c~ D,

cross section area gas only. sound speed constant in Arrhenius pyrolyss law. equation (17) ratio of particle to gas specific heats, c.,/cp constant in steady state burn rate, equation (8) particle drag coefficient erosive burning constant, equation (8) specific heat of metal oxide specific heat of gas at constant pressure specific heat of solid propellant diameter of the ith particle group activation energy of surface reactions Fp, particle- gas interaction force per unit volume for the ith particle group, equation (12) k thermal conductivity ! perimeter of the grain M,, Mach number at nozzle entrance m mass flux from burning surface N number of particle groups Nu Nusselt number n pressure exponent in steady state burn rate law n,. exponent of pressure dependence on surface reaction rate P pressure Pr~f reference pressure in steady state burn rate law Pr Prandtl number Qt heat release, per unit mass, in flame Qp, Particle--gas heat transfer rate per unit volume, for the ith particle group Q., net heat of reaction for processes at burning surface ~ net heat of reaction at T R gas constant R0 universal gas constant R< Reynolds number for the ith particle group, based on particle diameter, and particle gas relative velocity r linear burn rate T temperature 1]~ reference temperature in viscosity law T~ backwall temperature of the propellant t time u velocity u, velocity of the combustion products as they enter the main flow u~ axial component of u, u, threshold velocity x axial distance fi total particle to gas weight flow ratio. Eft, fi, particle to gas weight flow ratio of ith particle group gas only isentropic exponent O angle with respect to axial direction h~ thermal diffusivity iI viscosity p density f~ nondimensional frequency, ~.~h/g2 ~,J mass burning rate, per unit length, per unit cross sectional area % power in viscosity vs temperature law Subscripts f flame g gas P ith particle group .~ at the burning surface Superscripts ()' denotes a perturbation (-3 denotes a steady state value