b t 

z r. 

3 

IMPROVED TWO-DIMENSIONAL-KINETICS 

COMPUTER PROGRAM 

. NAS8-35931 

Third Quarterly Pr gress Report 

1 July 1984 through 30 September 1984 

d3-16? 

, 

N&+# -P& /7/..7v 

(NASA-Ch--) IHPITCVED N87-11508 

TWO- DIfl Ebl SION AL-KI NETIC S CCMFE IE li PRUG B AM 

Cuarterly Progress Report, 1 Jul - 30 Sep. 

IS64 (Softuarf and Enqineerinq Associates, Unclas 

Inc,) 23 p CSCL 09B 63/61 43852 

Pr e pzr e d by : Sof%ware and Engineering Associates, Inc. 

1560 Brookhollow Drive, Suitz 203 

Santa Ana, California 32705 

1 
IMP ROVED T WO-DI MENS IONAL-KINETICS 

COMPUTER PROGRAM 

NAS8-35931 

Third Quarterly Progress Report 

1 July 1984 through 30 September 1984 

Prepared for: George C. Marshall Space Flight Center 

Marshall Space Flight Center 

Alabama 35812 

Prepared by: Software and Engineering Associates, Inc. 

1560 Brookhollow Drive, Suite 203 

Santa Ana, California 92705 

1 
IMPROVED TWO-DIMENSIONAL-KINETICS COMPUTER PROGRAM 

1 .O BACKGROUND 

Future Orbit Transfer Vehicles (OTV) presently under 

consideration need rocket engines delivering a high specific 

impulse. This high performance can be obtained with large 

area ratio thrust chambers using oxygen with either hydrogen 

or hydrocarbon fuels. In the projected nozzles the 

combustion products are expanded to low pressure and 

temperature levels at high Mach-number flow, a domain which 

nas not been experienced with exisiirig ~ucket engines. 

Some modifications have recently been incorporated in 

the Two-Dimensional-Kinetics (TDK) computer program, such 

as: Condensed phase flow simulation, treatment of shock 

waves induced by the wall curvature, and the coupling of the 

TDK code with the Boundary Layer Module (BLM) for an 

improved thrust chamber specific impulse prediction. 

Recent boundary layer calculations for large area ratio 

nozzles have revealed that the boundary layer becomes very 

thick in the region of high Mach-numbers. This result is 

quite different from rocket engines built up to now and 

requires an examination of the existing thrust . loss 

calculation method due to the viscous effects adjacent to 

the wall. Furthermore, the knowledge of the Knudsen-number 

is desirable which differentiates between continuum, slip 

and free molecular flow, and thereby identified when the 

analysis, presently based on Newtonian fluid flow, becomes 

questionable or invalid. 

2 
In order to obtain highly accurate performance 

predictions, the best supporting data for equilibrium and 

finite rate chemistry must be prepared. Furthermore, the 

large area ratio nozzles require many more calculation steps 

to cover the complete nozzle flow field. This increacses 

computation time significantly and favors an error 

accumulation which should be reduced as much as possible. 

Another point of importance connected with the OTV is the 

existence of ionized combustion products which permits 

vehicle tracing or interferes with the transfer of 

communication or command signals. 

m L - LII~ current effort of snocK wave rnoaeiing inaucea by 

the curvature of the wall should be expanded to permit the 

shock simulation caused by wall contour discontinuities. 

Also the Mach shock existence in the center of the flow 

field interacting with other shocks as well as the 

analytical treatment of multiple shock waves may 

significantly affect the interior nozzle flow and associated 

performance and need to be simulated. 

3 
2.0 OBJECTIVE -------- 

The objective of this effort is to increase the 

analytical capability of the existing TDKIBLM computer 

program for performance of OTV thrust chambers. Areas which 

need further examination and improvement are the thick 

boundary layer inviscid core flow interaction and the 

related thrust loss calculation. To warrant highly accurate 

results the provision of the best available program input 

data is mandatory as well as the use of sophisticated 

modeling techniques to reduce computation time with advanced 

error control criteria. The simulation of wall shocks, 

Mach-snocKs in adaition to S~~CKS induced by large concave 

wall curvature is essential since this interaction with each 

other produces flow field changes which in turn affect the 

nozzle performance. 

4 
The program plan for this work is presented in Figure A. 

Progress in each of the work tasks during the reporting period 

is presented in the sub-paragraphs 3.1 through 3.5. 

3.1 TDK-BLM INTERFACE 

Interaction between the inviscid core flow and the wall 

boundary layer is expressed through 1st order by the definition 

of displacement thickness, 6". We have rederived 6" and e for 

'cne nozzie waii boundary iayer situation for tne purpose of 

obtaining a more accurate representation, especially when the 

boundary layer is thick. In particular the radial coordinate, 

r, has not been factored. Integral expressions are obtained 

that are somewhat different than the classical expressions. For 

example, the expression fo$ 6* is found to be: 

r +r 

2 

* - = f6(l - pur ) dy , r = 1 e . ----- le --_-- 'a Jo - 

'euer1e 

This is an implicit expression, since conditions at e are 

displaced from the wall by a distance' 6;. The subscript 1 

refers to a displacement equal to the velocity thickness. It is 

found that the boundary layer thrust deficit does not use this 

integral, but a somewhat different integral, as follows: 

r +r 

2 

- * "(1 - pur ) dy , r = 1 w . e---- 1w ----- 'b a Jo - 

Peuer 1 w 

This is a explicit expression. The subscript w refers to 

conditions at the wall. 

Similar remarks apply to the momentum thickness, 0. 
------- ............. %>., ......... ........... - .. _-. --.-.-. c ..-.._." -+:I---- I 

a c 

z 

0 

--_.._.. ............. 

0 - flra ..... 

-.*.-.. ... ;..Q -.., -~.-- L 

I I , I' 

_C.L._..-.. 

... - ....... 1 

-r 

I y_. _...__ __ . -___ 1. ......... ~ ...... 

3t- t . - I 

.~ . 

A . r ... - ............... +--.- 

.I 

-_. 

Q 

E 

0 

W 

I1 '4 a 

Y Mt 

nao 

bGV l-l 

L__---- -_-__ .... ~ .... -.-. ............. .__I__-_ ....... J : 
Task 3.2b, ODK Error Control 

This task was completed during the reporting period. The 

numerical integration method was rederived, and is attached as 

Appendix A. 

A new error control formula was derived that is based on an 

estimate of the absolute error for the integration step. This 

is done by evaluating the next term (3rd order) in the series 

expansion. The absolute error is converted to a relative error, 

unless the variable is too small, in which case the absolute 

error is used. The error term is used to control step size 

doubling or halving as shown in the attached write-up for 

subroutine INT. 

It was found that accurate integration of the differential 

equations requires consistancy in the nozzle cross-sectional 

area and its derivatives. Values of A(x), dA/dx and dZA/dx2 are 

required for the supersonic flow region. The subroutine used 

for interpolation, SPLN, was modified so that y, y', and y" are 

found by cubic, parabolic and linear interpolation, 

respectively, such that y' is the integral y'l, etc. 

The method was tested using the TDK analysis for the SSME. 

Plots describing the supersonic wall contour and its derivatives 

are shown in Figure 1. Local error vs. x/y* and step size vs. 

x/y* are plotted in Figure 2 and 3. The maximum step size was 

input as .2. The run time for this case was 25 seconds on the 

VAX 780. Previous to implementing the new error control method, 

this case required 127 CPU seconds. The old error control 

forced a fixed step size of .005. Thus, the new method gives a 

factor of 5 improvements on run time for this case. 

The success of this task will prove very useful when task 

3.2a; TDK with ODK Tables, is done since the ODK computer time 

has been reduced substantially. 
6 

..e 

e 

0 

am 

d> 

N 

Nh 

0 

C 

iS 

k 

1, 

I? 

I 

! 

il '9 

OD 

I 

A 

X 

v) > / 

I I 1 

tl/ Hao' O 

02' I 08' 0 

X 

hl 

X 

3 

\ 

h hl 

-0 

00 

i 

I 1 I 1 I - 

0' 09 O'Ot 0'02 0% 00-'s 00' t 00'2 00'0 oo.2-0 

H AZO 
U 

h 

\ 

X 

cn 

7 

i 

ozoo~o aroo.0 9roo*d r!oo*o zroo*o croo~o eooo-i 90OC~O P000'0 2000~0 0000' 

tlBti83 ?ti301 

a 

? 

aa 

0 

(0 

0 

0 

* 0 

0 

e4 

? 

0 

0 

9 

0 

c 

p: 

9 

I' 

0 

v 
I 

021'0 0or.o 080'0 

l 

I 

090'0 oro*o OZO'O 000 091'0 OPi'd 

3ZIS d31S 

0 

N 

4 

N 

0 

.O 

N 

4 

0 

m 

s 

0 

? 

rn 

0 

.+ 

4 

0 

N 

3 

3 

3 

? 

3 

N 

? 

a 

? 

c 

I 
ORIGINAL P.’”., t22 

OF POOR QE-rn:.: h ‘1 

-------------- SUBROUTINE INT 

Provides control for the implicit integration procedure, 

determines the proper set of nonhomogeneous equations to solve, 

and, after each integration step, computes the next integration 

step size according to the following relations: 

I 

hn+2 = 2hn+ 

1 -h hn+2 2 n+ 

1 ----------------------- ki ,n+l - 2ki,n +k i,n-1 

3Ymax 

I ’’ma, 

6 

20 

---- < 

MAX 

>6 

MAX 

where = max [lyi,nl,10-51 ;i - 1,2,.. NSP+3 ’max 

On option, (JF=l) only the fluid dynamic variables are used 

in determining the next integration step size. 

If the step size is halved for the fourth step, the 

The correspondence between equation number and physical 

i-ntegration is restarted using one-half the original step size. 

property is: 

Property -------------- Equation Number -------- 

1 Velocity of Gas 

2 Density of Gas 

3 Temperature of Gas 

4 NSP+3 Gaseous species mass fraction 

(1 ... NSP) correspondence to 

(4 ... NSP + 3) 

When the flow is supersonic, continuity is used to control the 

integration step size to insure that: 

/ (PVA)N+l - (PVA)N I 

I ( pVA) N+ 1 

< C0NDEL i _-_-____-_-__-___ I 

where C0NDEL is an input relative criterion with a default value 
This task has been completed except for adding a default low 

temperature data file to the program. Section 6.1.1 of the TDK 

manual has been rewritten to describe the usage of this data, 

and is attached as Appendix B. Note that only Cp vs T is to be 

input. 

Example, plots obtained from the method are shown in Figures 

4,5, and 6. The properties are for the species H2. Figure 4 

shows the input table of Cp vs. T. Figures 5 and 6 show the 

integrated values derived from Figure 4 for S and h vs. T. 

't- 1 I I 8 1 

300 .o 300 80 400 =O 500 SO sbo .o 700 -0 300 0 T 

2 Figure 4: C vs T €or H 

P 
0 

0 

aD 

0 

0 

-t 

0 

r 

0 

0 * 

I 

0 

0 

co 

I 

0 

hl 

a 

1 1 1 1 1 1 

0.0 300 BO 400 a0 500 0 800 a0 700 10 800 I 0 

T 

2 Figure 5: S vs T for H 

1 1 1 1 1 1 

'ZOO IO 300 e0 400 .o 500.0 - 600.0 700 10 900.0 

2 Figure 6: H vs T for H 1 
Task 3.2f, Update JANNAF Data 

The JANNAF thermochemical Data has been updated by obtaining 

a complete new data file from NASA/LRC, and attaching it to the 

program. The letter of transmittal is attached. 

Task 3.2g, Update Kinetic Data 

The reaction rate study has been completed and the results 

were sent to the COR at NASA/MSFC in September. 

3.3 NOZZLE SHOCK WAVES 

The plan for doing Task 3.3c, multiple shock waves, has been 

approved by the COR. A small amount of work was done on this 

task during the quarter. 

3.4 VERIFICATION 

No work was done on this task during the reporting period. 

3.5 DOCUMENTATION 

Work was begun on this task near the end of the reporting 

period. 
All remaining work is to be completed during the next 

reporting period. 

5.0 CURRENT PROBLEM AREAS .................... 

The rate of spending on the contract is too low, about 20% 

per quarter rather than 25% per quarter. It is clear that a 

time extension will be necessary. A revised schedule is being 

prepared. 

6.0 COST STATUS --------- 

As of the end of the reporting period, the cost status of 

this contract is as follows: 

1) Total Cumulative Coats through 

30 September 1984 

2) Estimated Cost to Complete 

3) Estimated Percentage of Physical 

Completion of the Contract 

$63,300 

$55,700 

55% 

The costs on the contract to date correspond well with the 

estimated percentage of completion. 
Given a differential equation of the form dy = f (x,y), a 

family of implicit integration schemes is based on the following 

-- 

dx 

back-differencing differentiation formula. 1 

3 ' Yn+1 

h 

=v 

-x 

n+ 1 n 

=x 

etc. 

For example, using one term in the sum on the left gives 

+ hy;+l, which is the usual backward Euler method. Yn+l Yn 

In general, using the first k terms gives an implicit 

method of order k, with a local error approximately equal to 

the (k+l) t-erm. 

Thus, using the first two terms gives 

or 

with a local error approximated by 

--------- 

Ref. 1 Dahlquist, G., and Bjorck, A., Numerical Methods, 

Prentice-Hall, 1974. (see equation 8.3.12). 
The integration method used in ODK is a variation of the above 

second order method in that y'n+l is approximated by 

From (1) and (21, 

3 kn+l -fkn - =h(fn+ah+3k n n n+t ) 

2 2 

with 

kn VYn = Yn-Yn,l , 

'n = f (Xn,Yn). 

and - 

kn+ 1 'Yn+1 - Yn+l-YnP 

Solving for k n+t ' 

(4) 

The procedure can be extended easily for a system of n equations 

in n dependent variables. 

Given 

equation (1) becomes 

= fi(~n,~1,n9..,~n,n)+ afi(xn) (X~+~-XJ 

-- 

ax 
equation (4) becomes 

n 

k) 3 ki,n+l - ' ki,n (fi,n+ai,nh+ ' 'i,j,n j,n 

- 

j=1 i=l,. . ,n 

2 2 

and equation (5) becomes 

n 

(fi,n+ai,nh+ ' 'i,j,n k) j,n 

j#i 

j=l 

The maximum local error is defined by equation (21, 

1. - 

- ' (ki,n+l * ki,n+ki,n-l max 

l<i<n -- 3 
APPENDIX B 

6.1.1 THERMODYNAMIC DATA BELOW 30O0K. 

If the temperature at any point computed by ODE, ODK, or 

TDK is found to be below the Thermodynamic Data lower 

temperature limit, Tk, the polynomial curve fit data (see 

Section 6.1) will be extrapolated to obtain values for the 

thermodynamic data. The extrapolated data can be inaccurate. 

The LOW T CPHS data set can be used to input low temperature 

data for those species for which extrapolation of the JANNAF 

data is not accurate. Only Cp vs. T is input. Enthalpy and 

rr.C.-,n-... -...- -I-&- GI1cII upy ai = uubaiiied by integrati~ig 2p VY. T within the computer! 

program. 

The lower temperature limit, Tk, in the Thermodynamic Data 

supplied with the program is 300°K. Thermodynamic Data below 

the temperature, Ti, ma.y be input by data cards as described in 

Table 6-5. 

An example of this input is given in Table 6-6 which shows a 

card listing extending the Thermodynamic Data for an 02/H2 

propellant to 100OK. Data in Table 6-6 is taken directly from 

the JANAF tables (Reference 23), except for Argon which is taken 

from NASA SP-3001. 

------- 

Ref. 23 Stull, D.R., Prophet, H., et al., JANAF 

Thermochemical Tables, Second Edition, NSRDS-NBS 

37, National Standard Reference Data Series, 

National Bureau of Standards, June 1971. 
Table 6-5: Input Specifications for Low Temperature 

Thermodynamic Data. 

The first card is a directive card that identifies the start 

of the low temperature thermodynamic data input. It reads as 

follows, columns 1 through 10: 

LOW T CPHS 

The next card contains a species name consisting of 12 

characters, or less, left justified to column 1. An integer is 

placed in column 21 indicating the number of T, Cp values that 

follow. A value of 1,2, or 3 can be used. For example: 

Next, the thermodynamic data for the species named above is 

input. There must be one pair of T, Cp values per card. These 

cards are numbered consequently in column 45. They are read as 

2F10.OS20X,I5. For example, for H20: 

100. 

200. 

7.961 

7.969 

1 

2 

Species name cards and thermodynamic data cards for other 

species, if any, follow. 

A final directive card is used to identify the end of the 

low temperature thermodynamic data. It reads as follows, 

columns 1 through 14: 

END LOW T CPHS 

9 
TABLE 6-6. LOW TEMPERATURE Co ,H;,S; DATA FOR AN 02/H2 PROPELLANT 

PT 

LOU t CP”5 

AR 

1oo.o 4,9661 

200,O 0.96ai 

H 2 

lOO.0 4,968 

200.0 4,968 

Ht 2 

100.6 5.393 

200.0 6.518 

H20 2 

100.0 7.961 

200.0 7,969 

NZ 2 

100.0 7,074 

200.0 6.989 

0 2 

100.0 5,666 

200.0 5,434 

OH 2 

lO0.0 7.567 

02 2 

END LOU T CPHS 

200,o 7,309 

lOO*O 6,958 

200.0 6.961 

2 

1 

z 

1 

2 

1 

2 

1 

2 

1 

2 

1 

z 

1 

2 

1 

2