Mathematical models and computer program on nonlinear combustion instability in liquid propellant rocket engine

Burstein, S. Z.

Chinitz, W.

English

Gas dynamic model for nonlinear combustion instability analysis in liquid propellant rocket engine

NASA Technical Reports Server

Nonlinear combustion instability in liquid-propellant rocket motors

Fluid dynamic model and diffusion-controlled flame analysis data on nonlinear combustion instability in liquid propellant rocket engine

Nonlinear combustion instability in liquid propellant rocket motors

Nonlinear Combustion Instability in 
Liquid-Propellant Rocket Motors 
Samuel Z. Burstein and Wallace Chinitz 
Fourth Quarterly Report 
to 
Jet Propulsion Laboratory 
July 30, 1968 
JPL CONTRACT 
951946 
"This work was performed for the Jet Propulsion 
Laboratory, California Institute of Technology, 
as sponsored by the National Aeronautics and 
Space Administration under Contract NAS7-100." 
https://ntrs.nasa.gov/search.jsp?R=19680023323 2020-03-12T06:35:05+00:00Z
. 
"This report contains information prepared 
by Mathematical Applications Group, Inc. 
(MAGI) under JPL subcontract. Its content 
is not necessarily endorsed by the Jet 
Propulsion Laboratory, California Institute 
of Technology, or the National Aeronautics 
and Space Administration." 
, 
i 
TABLE OF CONTENTS 
INTRODUCTION .................................... 1 
I . Toroidal Motor . Nonlinear Differential 
and Difference Equations ................. 2 
I1 . Boundary Conditions ........................ 7 
A . Injector face (Z=O) .................... 7 
. B . Nozzle exit plane (Z=L) ............... 10 
I11 . Results of a Test Case .................... 12 
IV . Program COMB Status ....................... 17 
ii 
INTRODUCTION 
This report describes the activities of MAGI during the 
fourth quarter of the contract with the Jet Propulsion Laboratory. 
During this quarter the effort was directed towards completion of 
the toroidal motor. In addition, some refinements in the pancake 
motor have been completed and will be presented. 
the ability to represent streakline data graphically is now avail- 
able as an operator option in the running of program COMB. 
In particular, 
The authors would like to acknowledge the strong programming 
effort of Harold Schechter during the course of this research 
program. 
1 
I. Toroidal Motor - Nonlinear Differential and Difference 
Equations 
The differential equations describing the motion of a com- 
pressible fluid in the coordinates e-Z-t can be given by the 
vector form 
II aw 1 a~ aH * aLnA - -+ - - + - + € I  - -  at r a e  az az 
The source term, +, has four components corresponding to the rate 
of production of mass, momentum and energy per unit volume. 
vector W also is a four vector, the components of which are the 
mass, momentum in the 8 and Z directions and total energy all per 
unit volume. The vectors G and H represent the flux of these 
The 
quantities in the tangential and axial directions respectively. 
The term containing the logarithmic derivative corresponds to an 
approximate way of treating small radial variations of the toroidal 
geometry. It is analogous to the treatment of one-dimensional 
time dependent flows with an attempt to include area variations. 
This term can also be considered as a source term of flux propor- 
tional to derivative of the logarithmic variation of the area, 
A = A ( Z ) ;  the proportionality constant being the flux H* = H-pSi2 
1 for i = 2 
'i2 =[ 0 for i = 1 , 3 , 4  
This term allows us to treat, in a simple way, a toroidal 
subsonic-supersonic nozzle which has an area variation about r = 1 
(without loss of generality). This will be discussed at greater 
length in the section on boundary conditions. 
2 
, 
The exact form of W is 
while the components of the tangential flux and axial flux are 
given as 
H =  ) 
(E+p) u 
( 3 )  
The pressure p is a function of the density and specific 
internal energy through the equation of state 
p = P(p,e) ( 4 )  
The total energy of the system per unit volume is 
E = p(e++(u 2 2  +v ) ) .  
Hence, by using Equation 4 and the above definition of the total 
energy, this allows one to write the fluxes G = G(W) and H = H ( W )  
only. We may differentiate G and H with respect to W and write 
Equation 1 in the standard form for quasilinear partial differential 
equations, i.e. 
n-1 
We associate the matricies I, aG/aW and aH/aW with Ao, A1 and A2 
respectively. In this notation our variables become Xo = t, X1 = 0 
and X2 = Z while the inhomogeneous term is 
B = H* aLnA - Q 
az 
3 
Equation 5 will be used as a starting point for the construc- 
tion of a stable difference approximation of the subsonic inflow 
condition at the injector face. This will be discussed in the 
section on boundary conditions. 
We seek a solution to Equation 1 by the method of finite 
difference approximations. Equation 1 is a conservation law. The 
difference approximation to 1 will also be a conservative difference 
form; the accuracy of such an approximation will be at least 
second order. The reader should refer to NASA TR 32-1111 for an 
example of another application of the following method. 
On the periodic space D defined by t>O, 0s,;8.s2.rrI OsZiL, we 
introduce a uniformly spaced mesh or net which consists of the 
points 
Oi = iA0, Zj = jAZ, tn = nAt irjrn = 0,1,2,... 
The set of net points Dh is defined by 
Dh = {ei, zj, tnli=O,l, ...; j=0,1,...; n=O1l , . . . I  
The approximation to W(Oif Zj, tn) defined on the space D is given 
by V(Oil Zj, tn)= Vij defined on Dh. 
used is most easily written in two steps. 
n The approximation that is 
-n Let Vij be the first predicted value. It is given by 
4 
The final value is obtained from 
The system of equations (7) constitute the difference approxi- 
mation to Equation 1. One may verify that the System (7) is indeed 
a second order accurate approximation to Equation 1 by allowing 
6G = A1 6V and 6H = A26V. Here 6 corresponds to a spatial difference 
operator in 8 or Z and A1, A2 are to be taken as constant matrices. 
Then substituting 7a in 7b and using these definitions for the 
differences in G and H the result 
+ terms O(At*B) 
is obtained. 
about t = nAt in terms of first and second order space operators 
6 and 6 2 .  
order accurate. 
This expression is but the Taylor expansion for V 
This shows that the difference approximation is second 
5 
The periodicity condition in e on the solution W, i.e. 
w(e,z,t) = w(e+2n,z,t), 
reduces the problem of finding boundary conditions to just the 
inflow and outflow conditions at 2 = 0 and Z = L. This theory 
is presented in the next section.' 
. 
6 
I1 - Boundary Conditions 
c 
A.  Injector face ( 2  = 0) 
Only if we were to seek a solution in D where IZ(<m (rather 
than the finite interval O<Z<L) would boundary conditions not be 
required. Such a problem is called an initial value problem. 
However, since the range of Z is finite, boundary conditions must 
be prescribed at Z = 0 and Z = L. The natural conditions at the 
injector, Z = 0, is that the momenta m = pu and n = pv be 
prescribed 
m(e, 0, t) = mo 
n ( e ,  0, t) = no 
as well as the energy 
0 
Conditions 8a and 8b are used to simulate the boundary values 
that one would impose on a gas injected rocket motor. Only axial 
injection implies no = 0. Since we are only considering subsonic 
inflow, one variable cannot be prescribed on the boundary. This 
is most easily seen by writting Equation 5 in characteristic form. 
We rewrite Equation 5 as 
Wt + AlWe + A2Wz + B = O 
Let P be a matrix whose eigen vectors are that of Al. 
that W = PV and substitute in Equation 9 
Define V so 
(PV), + A1(PV)e + A2(PV)Z + B = 0 (10) 
7 
Carrying out the indicated differentiation, there results an 
equation for V 
Vt + P-'A2PVZ = -P'l p + P t V + A 2 P z V + A l  (PV)  e 1 
which is in the simple characteristic form. We obtained this 
result by premultiplication by P - I .  Now P - l A  P = A is diagonal 
with Bi the elements on the diagonal so that differentiation on 
the left hand side of Equation 11 is in the characteristic direction 
i.e., the total differential of V itself is the left member of 
Equation 11. 
The elements Bi are the eigenvalues of A2; u, u+c and u-c. 
The direction u-c is left running for u<c; 
istic runs to the boundary Z = 0 from the interior Z>O. The 
vector V has the components 
hence this character- 
p + v - p c  2 
P + PCU 
P - PCU 
p - v - p c  2 
Hence the data p - pcu appears as a total differentiated quantity 
in the direction u-c; it is given by 
( 1 3 )  (p- pCU)t + (u-c) (p- PCU) z = -v(P- PCU) 8 - PC 2 ve 
One simple approximation to Equation 13 is obtained by 
allowing 
8 0- 
then solving for the pressure by using forward differences for the 
time and 2 space differences ( 0  differences can be centered) 
Now that the pressure is known at the advanced time nAt on 
the injector face, Z = jAZ = 0, we solve for the density p from 
Equations 8a, 8 b and 4. The following result is obtained from 
the positive root of the quadratic equation in the density, 
An iterative procedure could now be used by redefining the barred 
quantities 
which is given above. The present status of the toroidal motor 
and F i n  terms of the latest estimate of the density, 
code does not do this, however it is to be included in a future 
version. 
all the data required to advance the solution to the next time 
Now that the density is known, as well as the pressure, 
9 
level is known. The next section considers the boundary condition 
at the plane of outflow, Z = L. 
B. Nozzle exit plane (Z = L) 
The simplest and perhaps best way of treating an outflow 
condition is to impose the same condition on all the character- 
istics: that they have positive slope (all Bi>O) when the flow 
velocity is u>O. This means that if at the boundary Z = L, the 
flow is supersonic, then all three characteristics with slope u, 
u*c will intersect the exit plane when originating from the 
interior of D. Then the boundary is completely specified by the 
interior flow field conditions. Conversely, information specified 
at the boundary because of the geometry of the characteristics 
can have an effect only on the field to the right, or outside of D. 
We impose this condition at the boundary by specifying a 
schedule of area variation in the axial direction, 6LnA/6Z. At 
some arbitrary point, say Z = / o  the chamber length, the area 
decreases simulating a converging nozzle section. At a section 
further downstream, R < Z t < L  the area then increases up to the point 
Z = L. 
section, of length L-Zt, produces a sufficiently strong supersonic 
outflow condition. The only conditions imposed on how one choses 
such an area variation is the restriction on an allowable flow 
rate for the steady state of the rocket motor. Hence, the ratio 
of the chamber area A(Z<I) to the throat area A(Zt) is determined 
by allowing the Mach number at the throat section to be unity for 
the desired rate of flow specified at the injector. The area at 
The throat of our motor corresponds to Zt and the diverginq 
10 
the exit plane of the motor, A(Z = L), is chosen to produce a 
Mach number of the flow sufficiently greater than one so that 
even under transient conditions the outflow will remain supersonic. 
Then the boundary condition specified at Z = L is 
w ( e ,  L, t) = w ( e ,  L-AZ, t). 
This corresponds to zero order extrapolation and is sufficient to 
produce smooth results. This completes the discussion of the 
boundary conditions. The next section presents some results of 
a calculation using these boundary conditions in the toroidal 
motor. 
11 
I11 - Results of a Test Case 
The previously described theory was coded and a test case 
was generated for the program. We wanted to see the convergence 
of the calculation to a solution from somewhat arbitrary initial 
data. The initial data prescribed is 
p = 1. 
Th 
oszsl!  
p = 1. 
u/fi = 0.4 
J v = o  
exit Mach number ras approximately two for the area 
variation chosen. The exact dependence of the density pressure 
and velocity was linear in Z, the endpoint values p(L), p(L) and 
u(L) being obtained from one-dimensional reversible flow tables 
at y = 1.2. Table one below shows the convergence of the solution 
at the exit plane of the combustor. 
Table 1 - Toroidal Motor Convergence 
Pressure Iteration cycle Number 
1.0 
0.57281 
0.52067 
0.51894 
0.52048 
0.52018 
0 
200 
400 
600 
800 
1 0 0 0  
12 
Figures 1, 2 ,  and 3 are CalComp plots of the density P ,  
2 2 4  pressure p, and velocity (u +v ) . In this test case v = 0 for 
all time so that the calculation is essentially one-dimensional 
as Figure 3 shows. In Figures 1 and 2 one observes a gradient 
in the density and pressure at the injector face. This gradient 
occurs through the low order accuracy of the treatment of the 
boundary condition at the injector face. As was stated in the 
previous section, an iterative procedure will be employed that 
will reduce the truncation error of the boundary calculation so 
as to improve on this inaccuracy. 
13 


. 
IV - Program COMB Status 
In this section we present some recent results of a test run 
using the pancake motor computer program. 
55 in the third quarterly report to the Jet Propulsion Laboratory, 
work was progressing on the ability of program COMB to produce 
instantaneous streakline pictures of particles moving in the r-9 
plane. This feature is now included in the computer program; a 
test case has been completed and the results are presented in 
Figures 4 through 15. 
As was stated on page 
Figures 1, 2 and 3 correspond to the initial data for this 
calculation, the density, the pressure and the velocity respec- 
tively. For this case the max (min) pressure is 590 psia (10 psia). 
Figure 7 is the velocity field after almost one cycle of rotation 
of the pressure wave. Figure 8 is a plot of a streakline obtained 
by integrating the motion of a point (starting at the midpoint of 
the first quadrant, 8 = 1 ~ / 4 ,  r=0.5) over the first 500 cycles of 
calculation. The point ends near its starting position moving in 
a clockwise fashion which is the same direction as the wave. 
Linear theory would have the point moving in closed elliptical 
paths. Figure 9 completes the motion of the particle through 1000 
cycles of calculation. Here, in Figure 9, the fluid motion can be 
considered to be 'steady'; at least the nonlinear aspects of the 
motion are well defined. The particle stays within the first 
quadrant of the motor. 
17 
Figure 7 shows the particle motion with the starting transient 
included. The velocity field at the instant t = 6.7596 is shown 
in Figure 10. 
Figure 11 is the nondimensional pressure field at 500 cycles 
and the dimensional pressure on the circumference of the motor at 
this same time is shown in Figure 12. 
the same corresponding quantities as Figures 11 and 12 but at 
1000 cycles, when the wave has completed almost two and one-half 
rotations. Finally, Figure 15 shows the density field at 1000 
cycles. The program is flexible in that it can print out any or 
every one of these graphed quantities at prescribed intervals of 
time. 
operator control via input data cards. 
Figures 13 and 14 also show 
The intervals and selection of graph combinations are under 
18 
c 

I- 1.13 - -63 -.38 --I3 13 a 38 * 63 . ae 1.13 
T=O 0000 0 CYChES 0 -00 REVS 
0 
* 2 3  
m GRFlPH NCI. 

. 
T= 3.1466 500 c p  ES .95 SEJS 

29 
3 I 


Nonlinear combustion instability in liquid propellant rocket motors

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server

NASA Technical Reports Server

NASA Technical Reports Server