Mathematical modeling of liquid-rocket combustion instability

Math/ Compur. ModellinK, Vol. 1 I, pp. 57&574, Printed in Great Britain

MATHEMATICAL

MODELING

1988

0895-7177/88 $3.00 + 0.00 Pergamon Press plc

OF LIQUID-ROCKET

M.A. Farvin and John Peddieson

COMBUSTION

INSTABILITY

Jr

Department of Mechanical Engineering Tennessee Technological University Cookeville, Tennessee 38505 USA A mathematical model of liquid-rocket combustion instability Abstract. is developed which is of intermediate complexity. The application of the theory is illustrated by three examples involving velocity-sensitive combustion instability. Keywords. Combustion velocity sensitivity.

instability,

Galerkin

Nonlinear combustion instability in liquid-fuel rocket motors is a highly complicated process. Mathematical modeling of this phenomenon can be carried out on several levels of complexity. Thus a variety of mathematical models are possible, each having its own uses. For final design calculations a highly accurate model is required. For preliminary design calculations or the study of idealized cases to gain insight into the physical processes at work to create instability, on the other hand, a less accurate model may be of great use provided that it has a relatively simple mathematical structure.

The mathematical balances

In the present work a mathematical model of intermediate physical complexity is developed. It results in a single second-order partial-differential equation (governing a velocity potential) which lends itself to numerical solution by standard techniques for hyperbolic equations and to approximate analytical solution in important special cases. This relatively simple mathematical structure is a result of definite choices concerning the level of physical complexity to be employed. This interaction between the physical details accounted for and the desired properties of the mathematical model is discussed in detail.

atP+?-

methods,

model

is

based

on

the

mass

(P;) = w

atpL*?.(pxGk)

= wQ = -w,

the linear momentum P(at;t;.G)t$p/v a;,+;,.G;,

=

the thermal

(1)

balances

= 0

0,

(2)

energy

balances

P(atT+G.?T)-(Y-l)p(at(lnP)+v.?(lnp))=O (3)

atTl+ce.ifT, = 0,

Three aoolications of the model are oresented. "These involve both closed-form 'and numerical solutions and illustrate the ability of the model to predict the basic features of velocity-sensitive combustion instability.

and the equation

of state

p = pT.

(4)

In the above equations t denotps time, 3 the gradient operator, P density, v velocity,w rate of mass addition per unit time due to burning, In writing p pressure, and T temperature. (3) all interphase transfer terms (2) and This is reasonable have been neglected. because interphase mass transfer is normally for instance, of primary importance (see, Powell 1970) and also necessary to achieve the type of model described in the introduction.

EQUATIONS

A two-phase flow exists in a liquid-fuel rocket-motor combustion chamber. The two phases are a cloud of fuel drops (liquid phase) and a gas consisting of oxidizers and products of reaction (aas ohase). For the purposes of mathematical modeling, both phases will be treated as continua. It will be assumed that the gas phase is a non viscous non heat conducting, calorically perfect gas with gas constant R and constant ratio of specific heats y. The initial quiescent state of the has associated with it a uniform P,. a uniform temperature T,, and a pressure pr = Rp,T,. These are used

perturbation

erence quantities tonondimensionalize densities, temperatures, and pressures respectivelyin the equations that follow. Similarly the speed of sound vr = (yRT,-)4, a typical length L, associated with the chamber geometry, the time t, = L,/v,, and the quantity wr = prvr/Lr are used to nondimensionalize velocities, lengths, times, and burning rates respectively. In the equations that follow all quantities are dimensionless, unsubscripted dependent variables are associated with the gas phase, and dependent variables having the subscript a are associated with the liquid phase.

INTRODUCTION

GOVERNING

method,

By combining (3a) resulting equation, relations

chamber density uniform as ref-

p = TY/(Y-l)

and the

, p = Tl/(y-l)

(4) and solving the usual isentropic-flow

(5)

By combining these equations can be derived. with (2a) and taking the curl of the result

570

Proc. 6th Int. Co@ on Mathematical

571

ModrIling

it can be shown that if the gas flow is initially irrotational it will remain so. Assuming this to be the case allows one to define a velocity potential $ such that ; = %$

(6)

Substituting (6) into (2a), integrating with respect to the spatial variables, and using (5) yields T = J-I

= p(Y-I)/v

=

I-(v-I)(at++$.G/e)

Substituting (6) and (7) the differential equation

into

(1)

(7) produces

w = 0

(8)

The burning rate w is typically regarded as a function of the gas pressure and the relative velocity between the phases. Equation (Zb) indicates that the liquid-phase velocity is determined by the conditions at the injector. This together with (6) and (7) indicates that w = w(e)

(9)

Equation (9) shows that (8) is a single partial differential equation to be solved for the gas-phase velocity potential. A formulation of this kind would not have been possible without the assumptions mentioned previously. An approximate only quadratic as

form of (7) and (8) containing nonlinearities can be written

T = l-(y-l)(at~+&$~~~/2) p = l-v(at~+~~.~~/2)+Y(at8)'/2 P = 1-(at~+~~.~~/2)+(u-2)(at~)*/2

(10)

and

+(I-(y-e)at+)

w = 0

where E may be interpreted as a measure of the maximum amplitude of the initial disturbance and will' later be used as an ordering In (12) superposed bars indicate parameter. steady-state quantities and superposed carrots denote perturbations. Substituting the steady forms of (12) into the steady forms of (6), (lo), and (11) yields

+Fj,.?(2at,+it+.?,/2) +(1-(~-l)(at~+C+.?+,2))(Y-2)/(Y-1)

T=?(z)+Ei(x,y,z,t)

(11)

Equations (10) and (11) have the virtue of retaining the most important nonlinear effects in (7) and (8) while greatly simplifying the form of the equations. They are, therefore, well suited to the solution of model problems.

In the equations discussed so far no decomposition of variables into steady-state and fluctuating components has been carried out and, in fact, none is necessary. These formulations can be employed to investigate instabilities which occur before a steady state is established. It is now desired, however, to obtain a set of equations governing the behavior of initiallv small disturbances superimposed on a state of established steady flow and combustion. For simplicity the steady flow will be assumed to be one dimen;;onal (with the fl ow direction denoted by . Toward this end let

; = l-(EV)Z/2 T = l-(v-l)(EV)Z/2

(13)

Equations (13) constitute the formulation of the steady-state problem and hold for It should be pointed out all values of E. that rather than specify w as a function of equivalent 5 and solve (13a) for 5; it is This to simply specify W as a function of z. will be-the approach employed in the present paper. For simplicity the perturbation problem will This together be formulated only for ~<
= 0

p = -,(at~+,(vat~+(~~.~~-(at~)*)/2)) P = -(at~+,(vat~+(~~.~~-(2_y)(at~)~)/2)) 'r = -(v_l)(at~+E(va,~+~~.~~/2)) ; = 3,

(14)

where it is understood that all dependent variables in (14) now represent perturbations. To complete the formulation of the stability problem a relationship between w and Q is A simple burning-rate law suggested required. herein is w = w(-npy(at~-jatm,)+nv(~~.~~-j~~~

.h,))

(15)

512

Proc. 6th Int. Conf. on Mathematical

where &(X,Y*z,t)

Modelling

in a manner equivalent Powell (1970). =o (x.Y*Z.t-T)

to

that

employed

by

(16)

constants. and n,,, nv, T and j are material Equation (15) (to the order of accuracy inherent in (14)) superimposes a pressuresensitive response of the type used by Powell (1970) and many previous investigators (terms multiplied by np) and a velocity-sensitive response of the type employed by Wong et. If j al. (1979)(terms multiplied by n,). = 0 w depends only on the current value of If i = 1 w deoends 0 (instantaneous resoonse). on the history of' $ (history-dependent‘ reThe sponse) in the following crude way. response exhibits perfect memorv back to t-T'and no memory beyond that. It-is believed model consisting of that the mathematical (14) and (15) is a useful one for a varietv of stability .studies. The governing equation found by combininq (14a) and (15) is a sinqle differential second-order hyperbolic partial equation containing quadratic nonlinearities. As such it can be solved by a variety of standard analytical and numerical techniques. For the sake of brevitv onlv the case of instacombustion purely velocity-sensitive bility will be considered in the remainder To obtain this special of the present paper. case one lets np = 0 and nv = n which reduces (15) to (17) No real system is, of course, purely velocity sensitive. One benefit of the mathematical model employed herein is that it allows the properties of such an idealized system to be determined and then compared with the behavior of real systems. In this way the systems associated with properties of real velocity sensitivity can be determined. Typical investigations focusing of previous on velocity-sensitive combustion instability are the reports by Prien and Guentert (1962), Burstein and Schechter (1970), Wong et. al. (1979), and McDonald (1979). APPLICATIONS

consider As a first application standing The waves in a cylindrical chamber (h=l). hand-wall boundary condition a&(I,e,t)

= 0

is satisfied

(19)

by the Fourier-Bessel

@(r,e,t)=fo(t)Jo(solr)+fl(t)J,(sllr)cos(e) +f2(t)J2(s21r)cos(2e)+...

(20)

where sab is the bth root of equation Ja(s)=O. Substituti ng a three-term truncation of (20) into (18) and aoolvina the Galerkin orthoaonalization for instaice, procedure " (see, Finlayson and Striven, 1966) leads to a system of three coupled nonlinear ordinary differential equations (not shown for the sake of brevity) governing the modal amplitudes fo, Figures 1 and 2 show typical and f,. fl) results computed using the initial conditions

f,(O)=O,

f1(0)=1,

f,(O)=0

f,(O)=O,

f,(O)=O,

i2(0)=0

(21)

where a superposed dot denotes differentiation with respect to t. The figures show the maximum absolute value of each of the modal amplitudes (which occur once per cycle) plotted versus time. Figure 1 represents a stable case and Figure 2 represents an unstable case. It can be seen that the modal amplitudes exhibit very strong growth in the unstable case. As a second application consider standing waves in a thin annular chamber (~<
r = 1

(22)

(18) to

a2~-a~~+E(w(at~+n((a~~)'-j(a~~~)')) +(T-I)a&at~+2ao~a&~)

Consider an annular combustion chamber with inner radius l-A and outer radius 1 (imolvina that the reference length L, is the 'oitet radius in this case) and uniformly distributed This steady-state combustion (w = constant). system is best described by polar coordinates (r,e,z) with the z axis being the axis of symmetry and the plane z = 0 being coincident with the injector face. For transverse vectors disturbances (perturbation velocity perpendicular to the direction of the steady flow) $I = $(r,e,t) and the combination of (14a) and (15) yields

series

(23)

= 0

In addition it will be assumed that the combustion response is instantaneous (j=O) and (last that the nonlinear convection terms two terms in (23)) are negligible. This reduces (23) to the approximate form (24) Truncating

the standing-wave

solution

0 = f~(t)cos(e)+f,(t)cos(2e)+~~~

(25)

it into at two terms, substituting (24), and using the Galerkin method results in two coupled nonlinear ordinary differential Solving these equations governing fl and fp. equations approximately by the two-variable perturbation method (see, for instance, Nayfeh, 1973) for the initial conditions (which happen to lead to a particularly simple solution)

(18) should be pointed out that the present formulation does not include nozzle damping. This could easily be included, if needed,

f1(0)=1,

fz(O)=O

f,(O)=O,

f,(o)=-l/(zp

It

(26)

Proc. 6th Int. Conf. on Mathematical

Modelling

yields

REFERENCES

f,=exp(-&/2)sec(n(l-exp(-&/2))/2

3'2)cos(t)

f2=-exp(-Ei&2)tan(n(l -exp(-&/2))/2

3/2 )sin(2t)/23'2

The secant and tangent both become when their arguments reach n/2. the corresponding arguments start and achieve maximum values of n/23/2.

infinite In (27) at zero

Thus for n<2% the disturbances eventually decay to zero while for n>2% they grow to infinity in a finite time. This aooroximate solution illustrates in closed fo;m the same type of behavior exhibited by the numerical results presented in Figures 1 and 2. In addition it can be seen from (27) that even in unstable cases the disturbance initially decays. This same behavior is present in Figure 2, although it is difficult to see to the scale of the drawing. Finally, consider an approximate traveling-wave solution to (24). Using a procedure similar to that discussed by Crighton (1979), (24) can be reduced (for the case of traveling waves only) to a differential eauation qoverning v=a&_ The solution of this equation subject to the initial condition

v(e,O)=v,(e) the work

of Mickens

Mickens, R.E. (1985). Exact Difference Schemes for the Non-Linear Unidirectional Wave Equation. J. Sound Vibr., 100, 452-455. McDonald

G.H.,

J. Peddieson, and M.B. Ventrice Stability Analysis of a Liquid NASA Fuel Annular Combustion Chamber. Contractor Report CR-159734, Tennessee Technological University.

(1979).

Perturbation Mayfeh, A.H. (1973). John Wiley and Sons, New York.

Methods.

(29)

Wong, (30)

(28) becomes v=sin(e-t)exp(-E&/2)/(1-nsin(e-t)(l -exp(-&/2)))

Finlayson, B.A. and L.E. Striven (1966). The Method of Weighted Residuals-A Review. Appl. Mech. Rev., 19,735-744.

(1970). Nonlinear Combustion Powell, E.A. Instability in Liquid Propellant Rocket Ph.D. Engines. Dissertation, Georgia Institute of Technology.

case of

v,(e)=sin(e)

(1979). Model Equations of Crighton, D.G. Ann. Rev. Nonlinear Acoustics. Fluid Mech ., 11, 11-33.

(1985),

v=v,(e-t)exp(-E&2)/(1-nv,[a-t)(l -exp(-E&2)))

and H.S. Schechter (1970). Burstein S.Z. A Study of High Frequency Nonlinear Combustion Instability in Baffled Annular Liquid Propellant Rocket Motors. Report, Mathematical Applications Group, Inc.

(1962). Priem, R.J. and D.C. Guentert ComInstability Limits bustion Determined by a Nonlinear Theory and a One-DimenNASA Technical Note TN sional Model. D-1409, National Aeronautics and Space Administration.

(28)

can be shown,using to be

For the special

513

(31)

It can be seen that the disturbance dies out for nl. Also, unstable disturbances initially decay. The pattern of initial decay followed by catastrophic growth exhibited by the solutions discussed above has been observed in many practical combustion cases of instability. The foregoing analysis shows that this is a property of a purely velocity-sensitive combustion This response. suggests that velocity sensitivity should not be neglected when attempting to model such behavior. CONCLUSION In this paper a mathematical model of combustion instability in liquid-fuel rocket motors has been developed. The model is one of intermediate complexity. The interaction between the desired form of the governing equations and the physical assumptions employed was discussed. Three applications of the theory to velocity-sensitive combustioninstability problems were given.

and M.B. Ventrice K.-W., J. Peddieson, Analysis of Combustion Insta(1979). bility in Liquid Fuel Rocket Motors. NASA Contractor Report CR-159733, Tennessee Technological University.

514

Proc. 6th Int. Conf. on Mathematical

Modeliing

P

.r

c,