A digital computer program for calculating nonequilibrium expansions flows of dissociated oxygen and ionized argon around a corner

A DIGITAL COMPUTER PROGRAM FOR CALCULATING NONEQUILIBRIUM EXPANSION FLOWS OF DISSOCIATED OXYGEN AND IONIZED ARGON

AROUND A CORNER

1

8 SEPt 19R

by

..

A DIGITAL COMPUTER PROGRAM FOR CALCULATING NONEQUILIBRIUM EXPANSION FLOWS OF DISSOCIATEDOXYGEN ANDIONIZEDARGON

AROUND A CORNER

AUGUST, 1963

by

J. Galipeau* and A. Takano

UTIAS TECHNICAL NOTE NO. 69 AFOSR 64-0075

..

..

ACKNOWLEDGEMENT

We wish to thank Dr. G. N. Patterson for his encouragement and interest in the present work.

The advice and assistance received from Dr. 1. I. Glass in preparing the present note are acknowledged with thanks .

The Institute of Computer Science. University of Toronto. was generous in providing us with the use of the IBM-7090 Computer and we are pleased to express our appreciation.

The work was supported by the U. S. Air Force Office of

Scientific Research. the Canadian National Research Council and the Defence Research Board.

..

'

.

SUMMARY

A digital computer program coded in the Fortran II language suitable for an IBM-7090 computer is described that permits calculations of nonequilibrium expansion flows of dissociated oxygen and ionized argon around a corner. The procedure of calculation is based on the numerical-graphical method of characteristics. Details of the derivation of the flow equati'ons and their applications to particular problems are given in separate reports that are cited in the list of references.

An overall flow chart of the program is included to indicate clearly the sequence of logical and arithmetic tasks which the program per-forms. A listing of the program is also given to reveal fully its sequence of commands and operations. The program is divided into five phases. The first two phases designate the input data, which contain the statement of the free-stream values and the prescribed wall deflection angle as weU as the control of calculation and printing. In the subsequent phases calculations are performed to derive the flow quantities. There is rarely any need to print all the results computed and therefore only the flow quantities required to adequately indicate the current flow properties are printed at appropriately selected points on chosen characteristic lines. The skipping or printing is controlled by the input data.

TABLE OF CONTENTS

Page

1. INTROC UTION 1

2. INPUT DATA 2

3. CALC ULATION IN PHASE 3 8

4. CALC ULATION IN PHASE 4 11

5. CALC ULATION IN PHASE 5 16

6. DISCUSSION AND CONCLUSIONS 17

REFERENCES 20

TABLES FIGURES

1. INTRODUCTION

Theoretical studies were made of nonequilibrium expansion flows of dissociated oxygen and ionized argon around a corner in Refs. 1 and 2. respectively. The D;,lethod of solution was based on the

numerical-graphical method of characteristics and the practical procedure of calculation was described in detail in Ref. 1.

The characteristic mesh was constructed by starting from points on the frozen wave head and the corner. The distance of two neigh -bouring points on the frozen wave head and the frozen Mach number interval prescribed at the corner determine the size of characteristic mesh in the entire flow region. A finite mesh size introduces an error in the values of the flow quantities obtained. As pointed out in Ref. 1, the accuracy of the results can be improved by applying an iterative calculation at each charact-eristic mesh point. However. in the present calculations, by using a

characteristic mesh of small size. sufficient accuracy for practical purposes was obtained without repeating the above iterative computation.

The system of fundamental equations was reduced to the system of characteristic finite difference equations, which were solved numerically subject to the given free-stream conditions along the frozen wave head and the boundary conditions along the wall surface. At each characteristic mesh, by using the known values of the flow quantities at two initial data poin.ts, the above difference equations were solved simultaneously to provide the values of the flow quantities at a new point. At the very corner. the flow is frozen and of the Prandtl-Meyer type and therefore the flow quantities corresponding to the prescribed frozen Mach number can be obtained analytically. However. these numerical values are calculated to provide the boundary values at the corner for the calculation of the nonequilibrium expansion flow between the frozen wave head and the wave tail. The numerical computations were run on the IBM-7090 Computer at the Instituté of Computer Science. University of Toronto.

In the present note. a digital computer program is described that permits calculations of nonequilibrium expansion flows of dissociated oxygen and ionized argon around a corner. An overall flow chart of the pro-gram is included (Table lIl) to indicate clearly the logical and arithmetic tasks which the program perforrns. A listing of the program is also given. (Table IV) to reveal its sequence of commands and operations. The program is stated in the Fortran II languagesuitable for an IBM-7090 computer.

The calculation scheme is divided into five phases. The first two phases (PHASES 1 and 2) are assigned to the statement of controlling the calculations and printing of the results as weU as to the input of the given free-stream values and waU angle. In PHASE 3, the frozen Prandtl-Meyer expansion flow at the very corner is computed to provide the boundary values for the calculations along characteristic lines originating from the corner

(betweenthe frozen wave head and the expansion wave tail. PHASE 4). The flow quantities behind the expansion wave tail. including those along the wall surface. are calculated in PHASE 5. where the boundary condition states that the flow deflection angle must be equal to the prescribed wall angle along the waU surface.

The computations are of course made at every point of the characteristic mesh constructed in the entire flow region. However. for practical purposes the results at denser mesh points in the vicinities of the frozen wave head, the wall surface (including the corner) and the expansion wave tail are required to adequately indicate the flow properties . Therefore, 1t is not necessary to select for printing the results at all mesh poLnts and. in order to save printing time. the results of the flow quantities are printed at appropriately selected points on chosen characteristic linea. Th.e skipping or choosing of points and characteristic lines where the output data are to be printed is controlled by the sta.tement of the input data (PHASES 1 and 2).

In describing the program a different notation from that

employed in Refs. 1 and 2 is used for convenience. For example. dimension-less flow quantities without a prime are always considered in the present note. The cbaracteristic lines originating from the corner and the points on the wall surface are numbered as CO. Cl .. C2 •...• CW and CW1. CW2, ... J

re-spectively; here CO is the frozen wave he ad and CW is the expansion wave ta11 (see Fig. 1).

As mentioned in Ref. 2. in the case of ionizing argon, only a few equations differ from the corresponding relationa for dissociating oxygen and the overall procedure of calculation is similar for every case. There~ fore. the flow chart of the program prepared can be applied to computations of both dissociating oxygen and ionizing argon. Thus, the program contains the statement of calculating formulas for every case and the selection of the gas to be computed can only be made by theinput data.

In Tables land II) the input data are listed and their detailed explanation is given in Section 2. The necessary equations in the ca,lculations in PHASES 3. 4 and 5 are presented in Sections 3. 4 and 5. respectively. AJJ. overall flow chart of calculation and a liBting of the program are shown in Tables ill and IV, respectively.

2. INPUT DATA

In PHASES 1 and 2. the given free-stream values of the flow quantities (with subscript 0) and the prescribed wall angle Qw (deg. ) are

transmitted to the program and the statement of controlling the printing of the results are made. Theyare called the "INPUT DATA"., which are Usted in Table I. The list in Table I is made up of 12 lines and contains 36 diatinct quantities. The items on. a single llne go on on.e DATA CARD. Therefore, the statement of the INPUT DATA consists of 12 INPUT DATA CARDSR con

-'

.

taining 36 quantities. When going .. from one case to the next, the difference

may be only in a single DATA CARD. For simplicity and ease of preparing and checking the INPUT DATA, each one of 36 items are written down on a DATA SHEET, from which the DATA CARDS are subsequently punched. An example of the DATA SHEET for oxygen is shown in Table li. where the

free-stream values of Case 1 with gw = -150 (see Table Ili in Ref. 1) are

listed.

Symbols and a detailed explanation of the INPUT DATA are

given as. follows (see Tables land li):

Case 11 KCYCLE IDGAS go w MAXIT

An INTEGER to identify the various cases prepared for the

IBM-7090 computation (see Table UI in Ref. 1 and Table II in

Ref. 2).

An INTEGER to serve as a PRINTING FREQUENCY, which

controls the SKIPPING or PRINTING of ANSWERS in PHASE 3,

where CORNER VALUES are calculated. Therefore, if

KCYCLE

=

10, only the ANSWERS associated with Mo,

Mo

+

10 ~Mc , Mo

+

20 AMc, ... will be printed. There is

rarely need to print all ANSWERS; instead, there is need for

flexibility of usage and economy of output or printing.

An INTEGER to identify the gas considered; if !DGAS

=

1, the

program performs calculations for oxygen and, if !DGAS::: 2, for argon.

The INITIAL Mach number, i. e., the free-stream frozen

Mach number (= Mfl, see Refs. 1 and 2).

The Mach number INCREMENT at the corner, i. e., the

com-putation interval of step size in PHASE 3. As shown in Table

UI in Ref. land Table U in Ref. 2, áMc {=AMf> was taken

as 0.01 for all cases.

The WALL ANGLE given in degree. When it is reached. the

calculations in PHASE 3 are ended.

A TOLERANCE which will determine the precision of the

com-puted value Mcw. i. e .• the Mach number corresponding to g~

(see (3.18»). The program uses a Newton iteration io find

Mcw (see Appendix D. 2 in Ref. 1). Thus, if _AMTOL :: 0.00001,

Mcw will be correct to 4 decimal places. This value is used for all cases (see Table li).

An INTEGER defining the maximum numher of Newton Iterations

allowable. but still providing with a computed Mcw as precise

a safety measure to prevent the computer from ever being caught in a loop when attempting to satisfy LlMTOL' The value of MAXlT =,50 is used, although an actual maximum ever needed for hMTOL = 0.00001 is probably less than 10. The program would print an appropriate comment if MAXIT should ever be exceeded and the AMTOL should not be satisfied. Po' T o. Po. 0<0'

T'

0

"" C

t l r

INITIAL VAL UES used in the problem. It should be noted that . the free-stream pressure Po. temperature T o and frozen Mach number Mo are the necessary and sufficient initial values. The free-stream values of the other flow quantities can be calcu-lated from equations listed in Section 3 (also see Refs. 1 and 2). Nevertheless. for convenience additional free-stream values

( Po.

0(0 and rol are listed in the INPUT DATA. Here, 0( 0

denotes the degree of dissociation or ionization. For

dissociating oxygen. ro (= rfl) is the frozen isentropic index given by (2.46) in Ref. 1 and for ionizing argon. ro = 5./3 (see (2.55) in Ref. 2). In the program the dimensionless flow quantities are always considered without a prime.

Constants occurring in equations (3.12) and (4.43). but ex-pressed as symbols in the program to allow to change them if needed. For oxygen

~ =

,Jvl

~D = O. 038272 (see (3.26) in Ref. 1) C =

î

n

= 0.097818 (see (3.31) in Ref. 1)

However. for argon these data are not needed (see (3.12') and (4.43') ) and therefore set ~ = 0, C = O.

Constants occurring in equation (3. 14'). but expressed as symbols in the program to allow to,them if needed. For argon

CT1

=

Texcl #1 = O. 73284

CT2 = (~l - T_exc

)/1J

_I = 0.26716

(see (A. 19) in Ref. 2). However, for oxygen these are not needed and therefore set CT 1 = _{0, C T 2}

=

O.

The MESH SIZE used along the CHARACTERISTIC CURVE CO (see Fig. 1), when starting the calculations in PHASE 4

( Ar = ~r'f' see Table lil in Ref. land Table II in Ref. 1).

nAr

NOTALK

The TOT AL NUMBER OF MESHES of SIZE .ör along the

CHARACTERISTlC CURVE CO (n Ar =

nr,

see Table IIl in

Ref. land Table II in Ref. 2). Thus, if Ar

=

0.01 x 10 12 and

n.o r = 200, then the length of the CHARACTERISTIC CURVE

CO.is n.ör . .6 r

=

2.0 x 10 12 (for example, see Fig. 1). .

An INTEGER which should be always set to zero. It was used

in order to limit or re strict the program during checkout. DlST ANCES (radial) selected to control the SKIPPING or PRINTING of ANSWERS along a given CHARACTERISTIC

CURVE (Cl, C2, ... CW and CW1, CW2', . .. , in PHASES

4 and 5, respectively) as shown below, in order to save PRINTING time.

NCYCB1, NCYCB2, NCYCB3

INTEGERS used for the same purpose and in conjunction with rBl and rB2. The coordinates x and y (dimensionless) of any mesh point on a characteristic line are used to find the distance

r

=

-I

x 2

+

y2 of the mesh point from the beginning or first

point on the curve. The program evaluates x, y, and.r. For

example, suppose that

rBl

=

O. 10 x 10 12 rB2

=

1. 0 x 10 12

NCYCB1 = 1 NCYCB2 = 10 NCYCB3

=

20

Then on a given characteristic curve which has been selected for printing, all mesh points such th at

will have their ANSWERS SKIPPED or PRINTED according .

to the PRINT FREQUENCY NCYCB1; hence in the given example. every mesh point satisfying

r ~ O. 10 x 10 12

will have its associated ANSWERS printed. Thus, in the range

o

< r ~ O. 10 x 10 12

no mesh point would be skipped. But in the range

only every 10th mesh point encountered on the curve would have its associated ANSWERS printed; and in the range

i. e. ,

r

>

1. 0 x 1012

only every 20th mesh point would do. Ordinarily. rBI and rB2 are selected such that

and

NCYCB 1

<

NCYCB2

<

NCYCB3

since denser ANSWERS will be needed near the corner and the wall surface. where physical changes are more abrupt and important.

FRACTIONS selected to control the SKIP PING or CHOOSING of CHARACTERISTIC CURVES in PHASE 4 as shown below.

NFCYCI. NFCYC2. NFCYC3

INTEGERS used for the same purpose and in conjunction with fMc 1 and fM

c2· Just as all the mest points on a given curve need not their associated ANSWERS PRINTED. it is not

necessary to select all CHARACTERISTIC CURVES for PRINT-ING ANSWERSo _ In PHASE 4. THE CHARACTERISTIC CURVES CO. Cl, C2, . . .• CW have a first point (at the corner) that is associated with Mo. Mo

+

AMC' Mo

+ 26Mc

, ...• Mcw' re-spectively. To designate only some of these CURVES for PRINTING. the: _{fractions fMcl and fM c 2 are chosen such that}

and the progra_{m finds two Mach numbers Mc land M c 2 such that}

and hence

Mcl = Mo

+

fMcl (Mcw - Mo) Mc 2

=

Mo

+

_{fM c 2 (Mcw - Mo>}

In practice. the characteristic curve CO is not printed since the flow quantities along it are constant and equal to the free-stream values. Curve Cl, which is the first significant curve.

•

is always selected for printing. But the curves C2, C3, C4.

have their first point (at the corner) examined for its

associated Mach number. If it is less than or equal to Mc 1.

the curves are SKIPPED or CHOSEN according to the

FREQUENCY NFCYCl. lf it _{is less than or equal to Mc 2.}

the curves are SKIPPED or CHOSEN according to NFCYC2.

lf greater than Mc2' then according to NFCYC3. Thus, one

mechanism SKIPS or CHOOSES CURVES and when CHOSEN,

another mechanism SKIPS or PRINTS POINTS along the

chosen CURVE. By PRINTED POINTS it is meant that the

ANSWERS such as P. T,

P,

etc. associated with these points

are printed. The FINAL CURVE of PHASE 4, called CW, is associated with 9~ and Mcw and is always printed by the pro

-gram. NCURVE1' NCURVE2

INTEGERS selected to control the SKIPPING or PRINTING of entire CHARACTERISTIC CURVES in PHASE 5 as below.

NCYC 1, NCYC2, NCYC3

INTEGERS used for the same purpose and in conjunction with

NCURVE1 and NCURVE2. In PHASE 5, the first characteristic

curve is CWl. The program always prints it. The subsequent

CURVES CW2, CW3, CW4, . .. win be SKIPPED or CHOSEN for printing. In PHASE 5, starting with CW2, the program

keeps count of the CURVES which have been calculated. Thus.

if it had just completed the entire computation of the CURVE

CW15, it would have tallies 14, i. e. , set TALLY = 14.

Therefore, TALLY is compared to the INTEGER NCURV1.

lf the TALLY is less than or equa.l to NCURV1, the curves in

PHASE 5 are SKIPPED or CHOSEN according to the FRE-QUENCY NCYC1, L e. , only every NCYC1-th curve is

CHOSEN for printing. When

NCURVl

<

TALLY ~NCURV2

only every NCYC2-th curve in PHASE 5 is CHOSEN .. When

TALLY finally exceeds NCURV2, only every NCYC3-th

Curve will be PRINTED. Once a CURVE in PHASE 5 has been

CHOSEN by the above mechanism, POINTS along it will be SKIPPED or PRINTED by means of the mechanism in PHASE 4 mentioned before. Thus, along a given CURVE, mesh points

will still be SKIPPED or PRINTED according to a FREQUENCY which depends on the distance of the current point from the

ICW An INTEGER which will arbitrarily cause the first ICW curves, following CW1, to be prînted. Thus, if leW = 3, the CURVES eW2, CW3 and eW4 will be PRINTED, no matter what other

specifications may state. Once ICW has been satisfied, the other specifications resume control. If ICW = 0, the speci-fications in PHASE 5 alone have complete control of printing.

3. CALC ULA TION IN PHASE 3

The flow field and the characteristic lines ( ~ and ~ families) are shown in Fig. 1. For illustrative purposes it is chosen that n Är = 8 and therefore, as described in Section 2, the length of the eHARACTERISTIC CURVE CO (frozen wave head) becomes 8 . ..6 r. The number of CHARACT-ERISTIC CURVE (C 1, C2, ... ) between the frozen wave head CO and the expansion wave tail CW can be determined by designating .àMc for given values of Mo and 9w . The frozen Prandtl-Meyer expansion flow at the corner is calculated in PHASE 3. The flow quantities obtained (CORNER VALUES) provide the boundary values for the calculation (in PHASE 4) along the characteristic curves originating from the corner. The flow quantities behind the expansion wave tail ew (along eW1, CW2, ... ) is calculated in PHASE 4 subject to the boundary condition along the wall surface that 9 = 9w . The practical method of these calculations was detailed in Ref. 1 for the case of dissociating oxygen and also shown in Ref. 2 to be applied for the case of ionizing argon. In the present and subsequent sections, the

method of calculations performed by the program is outlined. The equations required for calculations are also listed. These equations can be applied for both oxygen and argon. However, for ionizing argon the equations with a prime. musfbe uséd instead ofîhe" corresponding eq.uations without a prime, which can be applied only for dissociating oxygen.

In PHASE 3, by designating the froz,en Mach number incre-ment AMc(= O. Ol, throughout the present calculations, see Table II), the frozen Mach number Mc can be prescribed at the corner for a given free-stream Mach number (Mc = Mo

+

i AMC' i

=

1, 2, ... ). Through the equa-tions listed below, the values of the flow quantities corresponding to the prescribed frozen Mach number can be calculated. Wh en the wall angle 9w is given, the frozen Mach number Mcw (= Mfw in Refs. 1 and 2) correspond-ing to the final expansion angle (9 = 9w) can be found by uscorrespond-ing a Newton iteration (see (3. 18) ).

EQUATIONS IN PHASE 3

(3. 1)

where

r

= 7

+

3

cX

o o 5

+

d-_o for oxygen for argon

~M

2 - 1 - D arc tan e . D 1

}l'e

= arc tan

,JM_e 2 - 1 Óe 2 = r e

+

Sc 1

+

0.5 ( r 0 - 1) Mo 2 1

+

O. 5 ( Po - 1) Me 2

fo

T'o

-

1 1 T' 0 - 1

c

T 1.5 (1 e - e

--VIT

e) e -1/T e T_o for oxygen (3.2) (3.3) (3.4) (3. 5) (3. 6) (3. 7) (3.8) (3. 9) (3. lO) (3. 11) (3. 12)

1 Kc = _{T 2.5 e- 1 / Tc} c

ei

_{ec =} 1 "'Kc Pc + 1 _ 2

(o(e~

1 - 0(0 Wc -

ie

(1 +0(0) 1 - O(ec

(CT! )

(e

T2)

= ~+2 exp ~ x .r/o(o

[0/

ec2 Be

=

1. 5 + 1 Tc Tc

=

2.5 + 1 Tc Ac = 3. 5

+

1. 5

ei

0 = 2. 5 (1

+

~) 1 - 0(0 1 - olec for argon (3. 12') (3. 13) - 0<'02] for oxygen (3. 14) - 0(02 ] for argon (3. 14') for oxygen (3.15) for argon (3.15') for oxygen (3.16) for argon (3. 16') (3. 17)

FIND Mcw CORRESPONDING TO 9~J THE GlVEN WALL ANGLE. SUCH THAT

(3. 18)

USlNG A NEWTON lTERATlON UNTlL Mcw SATlSFIES THE GlVEN TOLER-ANCEAMTOL· MeI

=

Mo + fMcl (Mcw - Mo) M c 2 = _{Mo + f Mc2 (Mcw - Mo)} (3. 19) (3.20)

.

4. CALCULATION IN PHASE 4

The nonequUibrium flow between the frozen wave head CO and

the expanaion wave taU CW is computed in PHASE 4. By the INPUT DATA of

.Ar and nAr' the points on the frozen wave head CO, from which characteristic

network begins to be constructed, can be determined and therefore the

coordinates of these starting points are given by

Xc = Yc = 0 (corner points) (4. 1)

and

xoi = i .Ar . cos

ro

(4.2)

Yoi

=

i . .6r . sin~ö (4.3) where i = 1, 2, 3, . . . , n,Ar 1 (see (3. 4) )

ro

= arc tan

From these points, the characteristic linea of ~ -family are drawn with

angle of dl = - .)A- o' The characteristic lines of the other family ( "l -family)

are drawn from the corner points (see (4. 1) ) with inclination angles Sc2

given by (3.5) (Cl, C2, ... , CW~ see :Fig. 1). At the first intersection of

Cl and ~ -characteristics drawn from the first point on the frozen wave

head CO, the flow quantities calculated by using the met.hod described sub

-sequently (for details, see Ref. 1) with the known values at the corner

(CORNER VALUES, corresponding to Mo

+

.6Mc) and on the wave head

(free-stream values). The calculations proceed along the first CHARACTERISTIC

CURVE Cl and then continue to the subsequent CURVES C2, C3, .. . , CW.

It should be noted that, in the case of dissociation oxygen.

the partially frozen sound speed (see (4. 34) and (4.35) ) specifies the

characteristic directions in the non--equilibrium flow, except at the very

corner and along the wave head. where the frozen sound speed (see (3. 9) )

does (Ref. 1). Therefore, only along the wave he ad CO and at the corner

the slopes

8

_{1 and}

á

_{2 of} characteristic Hnes are determined by using (3. 4)

and (3.5), respectively. At general points in the flow field, except the wave

he ad and the very corner, the slopes of characteristic lines must be obtained

by using (4.4) and (4.5) together with (4. 34) to (4.38). However. in the case

of ionizing argon (Ref. 2), the characteristic directions are always specified

by the frozen sound speed (see (3.9) or (4.34') and (4.36) ) and therefore SI

and ~2 can be calculated by using (4.4) and (4. 5) together with (4.34') to

A method of calculation for a general characteristic mesh

(see Fig. land, for details, Fig. 2) was detailed in Ref. 1. By using the

known values at PI and P2 (INITIAL DATA points, see Fig. 2) the system of the characteristic finite clifference equations are solved simultaneously to

give the flow quantities at a.NEW point P3. In the following equations, the

subscript 1 and 2 are referred to the flow quantities (KNOWN or FOUND

INITIAL DATA) at PI and P2' respectively, and the subscript 3 to the flow

quantities obtained at P3 (ANSWERS). From PI and P2' the characteristic

lines of ~ and

"t

families. respectively, are drawn to give the intersection

P3' The slopes ~1 and ~ 2 of these characteristic lines and the coordinates

x3 and Y3 of P3 can be calculated by using

5

1 = - ~

+

9 1

<5

₂ = )12

+

92 xltan _{ó l - X2tanó 2}

+

Y2 - Yl tan ó1 - tan 62 (4.4) (4,5) (4.6)

Y3

=

Y2 tanó 1 - YftanJ2

+

(xl -x2)tan 61tanÓ2

tan

dl -

tan ~2 (4.7)

respectively. For convenience of controlling the PRINTING OF ANSWERS

(see INPUT DATA in Section2). the polar coordinates R3 and 9_{yx3 of the}

point P3 are obtained by using

9~x3

= arc

tan

Y3 - Yc

x3 - Xc

(4.8) (4. 9)

The characteristic distances

.6'

and a'l( (see Fig. 2) are obtained from the

relations

Ll' ::

~(X3

- X2)2

+

(Y3 - yÛ

2

.à1l

=

1(X3 - X1)2

+

(Y3 - Y1)2

(4. 10) (4. 11) and therefore the values of P3 and 93 at the NEW POINT P 3 can be caiculated from the following equations:

P3 = GlPl

+

G2P2

+

92 - 91 - F1À

'1 -

_F2A

l

G 1 +G2

(4. 12)

In order to obtain the value of 0( at P3 by using the mass pro-duetion rate equation (4.24). the streamline through P3 is drawn with the flow

deflection angle 93 given by (4. 13) and the interseetion P 4 of this streamline and PlP2 are found:

X.4 = xlY2 - x2Yl

+

(xl - x2) (x3 tan 93 - Y3)

.(xl - x2) tan 93

+

Y2 - Yl

(4. 14)

(4. 15)

(4. 16)

When xl - x2' the following equations must be used instead of the above equations: L1X

=

Y4 - Y2 Yl - Y2 (4. 14') (4.15') (4.16')

The flow quantities at P 4 ean be ealeulated by using the fol.lowing linear

inte rpola tions:

0(4= ci₂

+

ÄX (0(1 - 0( 2) (4.17) w4 = w2

+

AX (wl - w2) (4. 18) u4 = u2

+

.l1X (uI

-

u2) (4. 19)

f 4=

f

2 + ~x ~

f

1

- .r

2) (4.20) P4 = P2 + 6x (PI

-

P2) (4. 21) T 4 = T2 + _{4X (Tl - T 2 )} (4. 22)

The distanee along the streamline is given by

AS

=

,J(X3 - X4)2

+

(Y3 - Y4)2 (4.23)

and therefore the value of 0( at P₃ ean be obtained by using the following mass production rate equation:

(4. 24) Consequently, the values of the flow quantities at P3 (except P3, 93 andO< 3) are calculated from the following equations:

L

z4 = for oxygen only

T 4 2 z4 z4 e

P4= for oxygen only

(eZ4 _ 1)2

A4 = 3.5 + 1. 5 o{4 + (1 - d._{4 )} ~4 for oxygen

= 2.5 (1 + 0(4) for argon Z4 ) - - _{T 4} eZ4 - 1 for oxygen for argon

j>

3

=

P3 T 3(1 +<X3}

,,9

z3

=

T 3 for oxygen only

2 z3 ?3

=

z3 e

(eZ3 _ 1}2

for oxygen only

r

3 = 7.0 +3.00(3 +2.0 (1-o(3}~3 5. 0

+

0( 3

+

2. 0 (1 - 0(3) ~ 3 = 5 3 for oxygen for argon (4. 25) (4. 26) (4.27) (4. 27') (4.28) (4.28') (4. 29) (4.30) (4.31) (4.32) (4.33) (4.34) (4.34' )

ARG

=

r

₃ (1

+

rJ._{3 ) T3} a3

=

"'ARG M3=

~

a3 )'<3 = arc tan 1 A3

=

3.5

+

1. 5

oe

3

+

(1 - 0(3) ~ 3 for oxygen

=

2.5 (1

+

()(3) for argon

=

1.5+ for oxygen =2.5+.!. T3 for argon B = B3 T 3 T 3 3

c

T31. 5 (1 _ e-;;'/T3) e-1/T3 1 for oxygen

=

for argon for oxygen (4.35) (4.36) (4. 37) ,(4. 38) (4.39) (4.39') (4.40) (4. 40') (4.41) (4.42) (4. 43) (4.43') (4. 44) (4. 45)

=

(CT:!

+

2)

exp

(CT~2)

X .P3 2 0( 3 rO<e32 1 - 0(3

L

1 -

~e3

(4.45') (4.46) 5. CALCULATION IN PHASE 5

The flow field behind the expansion wave tail CW (along CW1,

CW2, . .. , see Fig. 1) is caleulated in PHASE 5, where the ealculating pro

-cedure and the equations given in Section 4 are uSllally applied, exeept when

the calculations go to charaeteristic mesh on the wall surfaee, as shown in

Fig. 3. The eomputation for such a charaeteristic rnesh ean be performed

as a special case of the program used in PHASE 4, where P2 and P3 are on

the wall surface (WALL POINTS) and P 4 is identical _{with P 2.} In th is section.

the method of caleulation is outlined for obtaining the flow quantities at the

WALL POINT P 3 (for details, see Ref. 1). The numbering of equations

listed in th is section indicates the corresponding equation used in PHASE 4.

In other words, for the calculations of the WALL POINT the equations listed

in th is section must be used in the program instead of those applied in

PHASE 4 with the respective numberings. (Not all of the relevant equations

will be listed in this section. )

In PHASE 4, the very first points on the CHARACTERISTIC

CURVES CO, Cl, C2, ... CW are just at the corner (CORNER POINT) and

therefore their eoordinates are given by (4. 1). However, the CHARACTERISTIC

CURVES CW1, CW2, ... in PHASE 5 originate from the points on the wall

surfaee (WALL POINTS, see Fig. 1) and therefore

Xc' Yc

=

x, y eoordinates of the VERY FIRST POINT

(on the wal! surface) of the eurrent CURVE (5. 1) In the eharacteristic mesh just on the wall surface, the alVEN POINT P2

and the NEW POINT P3 are on the wall surface (WALL POINTS, see Fig. 3).

Therefore, in this case the value of S2 at P2 must be given by

(5. 5)

me.asured from the VERY FIRST POINTS (WALL POINTS) of the current CURVES, so that for P3 on the wa.l! surface

(5.8)

9yx3 :: 0 (5. 9)

Since the point P2 is also on the wall surface(see Fig. 3). then the following

equationa must be used:

(5. 14) (5. 15) Ax= 0 (5.16) (5.23 ) (j. 3 =

ot

₂

+

w 2 AS (5. 24) u2

Particularly, when the NEW POINT P 3 is the VERY FIRST

POINT of CW1 (FIRST WALL POINT). the point P'? is just at the corner

(x2 :: Y2

=

0) and therefore the flow quantities at P2 have the CORNER

VALUES corresponding to Mcw (see PHASE 3). Since they are the frozen

values, the values of A2 (= A4) and B2 (= B4) must be given by

A2 ::: A4 :: 3.5

+

1. 5 0(0 (5. 27) (5. 28)

for oxygen only, instead of (4.27) and (4.28) (see (3. 16) and (3. 15) for the

CORNER VALUES). For argon, equations (4.27') and (4.28') can be

applied in the pres.ent calculations in PRASE 5 as weU.

6. DISCUSSION AND CONCLUSIONS

. A digital computer program coded in the Fortran II language

slJitable for an IBM·-7090 was described that permitted calculations of non

-equilibrium expansion flows of dissociating oxygen and ionizing argon around

a corner. A detailed derivations of equations. listed in the present note was

given in·Refs. 1 and 2, and particularly. the method of calculation was

de-tailed in· Ref. 1. The program was conveniently divided into five PHASES;

in the first two PHASES the INPUT DATA were designated and in the subse ..

quent PHASES calculations. of the frozen flow at the corner, the nonequilibrium

corner expansion flow (along Cl, C2 • . . . • CW) and the flow behind the

expansion wave taU (along CW1. CW2. ) were performed.

The INPUT DATA (PHASES 1 and 2) contG!,in the prescribed free-stream values and wall angl.e as weU as the statement of controlling the ca1culations and the PRINTING of ANSWERS. An overall flow chart of the program, which .indicates the logical and arithmetic tasks performed by the program, is shown in Table III, where the numbers written on each state-ment are referred to the addresses in the program. A listing of the above program is shown in Table IV, which reveals the sequence of commands and operations. For ionizing argon, only by modifying some equations used for oxygen, calculations can be performed by using the flow chart of the program shown in Tables UI and N. The selection of the gas to be treated in the pro-gram is made by the INPUT DAT A only (designation of IDGAS, see Tables I and II).

It is expected that the intersection of CHARACTERISTIC CURVES belonging to the same family (It{ -family) would take place in the present calculations to form a recombination shock wave (see Ref. 1). If

such an intersection would appear, the present method of calculation cannot be applied in the current characteristic mesh as well as beyond it, because the NEW POINT P 3 should go upstream of the KNOWN POINT PI

(x3

<..

xl' Y3) Y1). As a consequence, if the calculations would proceed be -yond the characteristic mesh where the above intersection appears. then physica11y unreasonable values of flow quantities wi11 result. Therefore, for practical purposes. by considering the flow conditions imposed by the present free-stream conditions listed in Table IU in Ref. 1 and Table U in

. Ref. 2, it is stated in the program that the following 10 rules must be satisfied:

o

<

p

<

10 X 10- 6

o

<.

T

<

O. 1.

a>

0 M> 1 w

<0

current CURVE would STOP and PROCEED to the subsequent CURVE. In

this case the program would PRINT an appropriate comment.

The time required for calculations on an IBM-7090 depends

of the size of characteristic mesh ( ~r and LlMc), the length of the frozen

wave head (= _{n-o r .,Ar) and the wan angle (9w) chosen} in each case of the

calculation. In the computation with n br

=

200, .Mi_c

=

0.01 and 9_w = -150,

it took about 10 minutes on an average to calculate one case. The time re-quired for PRINTING can vary with SKIP PING or CHOOSING of POINTS or CHARACTERISTIC CURVES PRINTED. Usually, in order to save the

PRINTING time, denser POINTS were PRINTED only in the region close to

the frozen wave head, the wave tail and the corner, where the variations of

the flow quantities were significant and important.

In conclusion, the present program was found to be very use-ful for calculating nonequilibrium expansion flows of both dissociating oxygen

and ionizing argon. By changing the statement in the INPUT DATA CARDS,

the program permitted computations for various free-stream values and wan angles, as wen as for various characteristic mesh sizes.

1. Glass, I. 1. Takano, A.

2. Glass, I. I. Takano, A.

REFERENCES

Nonequilibrium Expansion Flow of Dissoc~ted Oxygen Around a Corner, UTIA Report No. 91, June 1963.

N onequllibrium Expansion Flow of Ionized Argon

TABLE 1

List of the Input Data

Case #

KCYCLE

IDGAS

Mo

.6.

Me

gO _w

_MTOL

MAXIT

ro

_. (j. 0

T o

Po

fo

~

C

_CTl

_CT2

6r

_n.ór

NOTALK

rBl

rB2

NCYCBl

NCYCB2

NCYCB3

fMel

fM e2

NFCYCl

NFCYC2

NFCYC3

NCURVl

NCURV2

NCYCl

NCYC2

NCYC3

1 2.5785 1. 46143 0.038272

o.

01E 12 O. 10E 12 1 0.25 10 20 10

o

10 0.01 0.19964 O. 097818 200 1. OE 12 50 0.50 20 60 50 TABLE II Input Data Sheet 1 -15.0 0.062636

o

o

50 50 50

OXYGEN

1

0.00001 50 2 0.20597E-6 2.7411E-6 3

o

4 5 6 7 8 9 10 11 12

'rABLE lil

SL T-I FLOW eH/IR T: FROi<1 p/tGt:-S: _ _ _ I L<J/II> '''''P,"'~ PIITA I FflCl1( PAGt!-Z-:- - .-d._e Ac. ~.j re Ac.

r;.

$I\V.e ~5 A I~-Il/~f~f) RJ!c~,?J) Tm:. cVI?~J!A/r

Con.NI1f'1. VhL<.h!$: Me f"c f'"'! Be 0; f.,. f'"

'"'e

Ge:. '->le

7; J~2 J~:

.lL'é { TI(CJ 5 ;

PR/A/TIN& ïl7'"'L~:;

1

.', /5-If,'tfli!l> flt::coR I>

~ cof\'t!-3 BLeeK

TABLE III - 1

IA/crU"115,! 11-,,, il~ AII"'e

2+ _{TABLE III - 2}

1. ~

tic -ra

1Jy IV""' r.", • T"'(II TIO 111 I~'ND;

I'RIN T c"IOf",~lIIrl Me cUINr!~"'''' PING 7lJ

tr

_e - (l"

....

t

~AL<:UL.llrl! ANI SA"~ r",~ 1~'-W.t?D

R"~"

"'»<-,

C"~~5J i1E~n~ Ccl{fU~~hAl/)lt4I(/, r" Me ~ MclN' !'tv' C "/fA ,'I! MII(!e 31 ( I

us".,;; 1JwUR. ~ mp" ... r,ll' II"'I!.J

-!--z,t!C.ifL/'s'11 CoRI! C:-,4PAC"1 T'V

';:.~~::('"'

, . .

R"~

::~~~~~z:/;:wLh;:sr:1

FIA//) : Me, /i( M~2.1

'53 45/ JJa m/fl!< < j33

f

'1éS WC 5Î/~V T. Cc/{IfS ' . '58 IPlfIlS/!-.J r-",'D(.fj)

I

7ll!GIIV PlfIlS~·+ Ji~,,;;:.i~·-C'-'f/''''I.!

I

'S~r _?C'" ; ~'(. =01 cu rl'" iJ;' '~~

e.

Je,;.

I

.J

I r-.Il i...; ' .. :J./ :J", ' " (11A.k '

~

[ ~/lL~~'L"T(! /J,LL

/17//,:e Cufl.vt!s C/.J C2 _.I < ~-' .' , _.I ('W iJr p//IIS/ó'-f-

I

Till/S: F'RnM PAGe-3:_

-

_- ,

r '~

..3ao~ OR .3~o+

IJ ~'NG- TTII! 1(111"' .... / ~I I, ;:) ., ~ "Pc" V/IL..}~S oi-lP () 1C _:1

r/lv D. .I, _{4'1. '):" '/,1 R3}

0;"

_{A ~ IJ",}

-!-U",A/G n//! f(1v°(iW/V' • p,"

e

" 1'2" VIlt.UI!S C;/~

G

-+

r-

0< A - ft. (" I

r/lv j} " .f~ (;" ~~1:

_:t+

. .

~: 0(+ .u,-".. A+- ('~ 11- ~

AS '-'<i Z _..- (lot- Á+- J3+ _.!:::t. T; ~J z.~ (3"

r;

/fIN. dJ /113 1'-" A~ ];, _,l G_~ I(j o(e _b.,r'

--

- ----1 A

--~

"i/I,r:~/ 'IC, ~)C ... ) ,4< rell,NATtvIF. t!C;cH lfaA/,S AII./é u~r:p Tc po,'p

")(1- !1 ... L\X

rNt: VAI?f /)1I,,~ S L:",'PI:!.l. Lllv'tEiJ (A/ I?t!/)

Mil! r S/lTISI't' TIft3 10 I?Ul-r! S ar- P/cGt!

2-1-T"Hû ":i;

---~---

_.--I

IMI't~"'h"""L. 'Y FrrI!V'OI/!>

PO,,,,r p'"

AAltJ o~ eJ~ql!Alr

Ir",:! L TY PIJ rw r

"J

Tl/!! tJi.//tf lS.é~ CF Pil'" $

eN 14 Af V ct/R.VI! rllAr r-OLL'/I!$

I\I,~L 'ZIt! P~'Po'lTI"",'I\ r~LY ,-éss PII"~~ - 1- /EII/J/!J) '8éG/A/ PtVIISIf -5 '''' R3 ~ .,.k~, Ir- ""11, < R. ~ ft~2. Ir' f? ~ ./r, 1'ABLE 1lI " 3 AL.w/'Y$ pf?INr ~1Yii~1---IIof 1'111.' T I"o,AI r 's .', ,,~T pr- - NCYC t" IIcycB2 A/CYC IJ3 A4'S"'~<? s : ;<J 'J3 1<,." .. , G

cr.-5~I.~CT eN'-1 F-VI!~y pp Pdll.lT' F;

F.,? MIN'T/AlG- IT~ A N'!JwteJl5 :

8;~:o ft;) G,

"'<e" éVJ':'f'l

Ir- Alo r ... i,fIt~/'D., PIIOJrl!D

Pilt"'/.' Trf~'Àk~W~RS..c'r- j rilt! 1./15 r .~.In ~ 41t. P4IAr, Rt:e;IIII.DL~S" .F PF 4;)

"",

1".: (9" J o<~ ""3 r-.. .3

.r. ..

... 11. LIA/I!$ o ro

~/!L~ T" nu! ",~.IC r cu~1/1!

AMONe;- Cl.; C2" C 3, ' • <, cW'

3002-':.0 R -300+

-WI"!'" tV~ I/ICLI'o ",JAIS Tè .TIf~ 10 ,?ULI!S oF /"'Á& - 2~

PlllIsff-S WILL /l1lC/e t11₄A:. cuRvéS CW2.

J 11AC/INt; 1é't!51'~<: r,vt!t...Y

t?1",~ J 1'YI1J,4c. -1.1 ~.k - 2) . , . r-1"C... p/)llI/iS CAl r/l(~,.,. TiI~ I/ell.'( L~S;- ,,"1 t?1 ti::. Ct.l~tI~ lIAS 'lJI.JI I PO/N";-•

.<IA IJ SI Nc:r ~ 17 FINO: W/l~N WHif",' ~

+

FIN P; T1ff3. ')( J:

.

p~ = p~

.. ~ .300a O/~ .JO/C

KIII(JCv',V 'i J~ - --W/4Lt.. f't:IAJ, k

..

~ 'P' 2 fllr:KN(J",,'AI"r:" ~ "PL " V/k,j'"" OF F o( ,<..<r""""' (' T _f!. _ (,I~ _")(~ _ ~~ _ " : : . 0( 't- ... ';..

lil/feil! ~i - WAL/.. POO.' r 0 P cvl«l/'t! c W I

.'. /!Lrr=lllIIl.rtVr= ':GUltrtOA/S 1J.'~t! IJSCO IC PlAID: A+ ~ 13+

Fofl. :

ÁLS.

fI1t..1S r SI1TISFY Tllr. IC RL/LES o r= PAGr= 2

+-r/lU!>

:-~l

- -

P'lIvr ,4A1S.wl6rt S OF

I,..lI1rEl1,1'>T7:!LY P/1.f5//OU5

".'101""

PS.)

AND oF CultltaA/r

FAUL.TY Po'kr PIJ

Til/! I/tJ~~ IF p</II,/rs

ON I/NV cuRve ""f}T FoL .. QW$

IVIL.L tJF- (",H Po~ TI.'" .. r!!L. Y LI!S5 15 Ge T4 PA~·I: ~ ______ -L ____ ~ Ir- R3 ~ ./(11, Ir- . .kt), < f?3 E "'+:a~ Ir 1?!> > Al' TABLE III - 5 ,4LINIlYS PIl/NT 1---1--1 Fm,si p",,"r's ,4Ai'SW"'RS: 'X') ';j' R, .•• ,.·G

r

.', ;!IéT 'Pr" NcyCBI 11 " " I, #CycBZ AlCVC63

~I!L~cr elolL't ~veq.v PF do. fclA/ï p~

Pcr,z. PRINrlA/G IT'S AIVSwrER5; 11', 'j~ _~3. e;)(3 ft;,

h

T; ('3 G, <><'e,3

13o/1!1t Y

Ir-"/Cï ALRrEl4DY fa,,,, ""~J).J PRIA/T Tl/a /III/.s.vr:~ S 01« Til/!. VEl'ly LilS r POlIoI""

REG-4RDL.rESS ol" pr-a._J .'"'3 Il. M, f'.; F. .) .f." 0:. I-IA/I!S t&" 0(3 Ge TO PAGe -4-; ,3008 'O~' 3010

'}-_____ +SI!,I!C i rl/I! /./f!:;t .. T

CuRvC-AMcJAlG e",,'I, ~'W2J cw3 " .

TABLEIV

*

*

C 10 1 2 3 C 15 1 2 3 16 1 C 17 C 18 19 1 2 3 4 22 1 C C 23 C 24 26 C 1 2 25 1 2 3 4 5 1 2 3 4 30 1 2 3 DR. I. GLASS FQ07+2 GLTI 10 4 ID CARDS

COMMON D, E, KFAZE5, KWALL

DIMENSION CORE3(45161, CORE4(681SI, CORE8(68181

NEW PAGE, PRINT HEADINGS WRITE OUTPUT TAPf 6, 10

TABLE IV - 1

FORMAT ( lHl, 39HNON-EQUILIBRIUM PRANDTL-MEYER EXPA~SIONI

lHO, 35HOF DISSOCIATING PURE DIATOMIC GASESI

lHO, 35HPROBLEM DR. GLASS, DR. TAKANO UTIAI

lHO, 23HPROGRAM J.GALIPEAU NCS I

LOAD I"JPlJT DATA

READ INPUT TAPE 5,16, '(ASE, KCYCLE, IDGAS,

VM0, DMC, THAWDG, DMTOL, MAXIT,

GAMO, ALFO, TO, SMPO, RHOQ,

V, C, CTI, CT2

FORMAT ( 3110, 42X I 2FIO.6, FIO.4, FIO.6, 110, 22X I

3FIO.6, 2EIO.5, 22X I 4FIO.6, 32X I

IF (IDGAS-21 17, 18.18

SEE IF GAS IS OXYGEN •••

GA.S = 6HOXYGEN

GO TO 19

••• OR ARGON

GAS = 6H ARGON

CONTINUE

READ INPUT TAPF 5, 22, DR, NDR, NOTALK,

RBl, RB2, NCYCBl, NCYCB2, NCYCB3,

FMC1, FMC2, NFCYCl, NFCYC2, NFCYC3, NCURVl, NCURV2, NCYCl, NCYC2, NCYC3,

ICW

FORMAT ( EIO.4, 2110, 42X 12EIO.4, ·52X I 3110, 42X I

2FIO.6, 52XI 3II0, 42XI 2110, 52XI 311n, 42XI 110, 62X I

IF NDR LESS OR

=

400, CORES ONLY.

OPTIMIZ~, HENCE DECIDE IF CORES ONLY •••

TF (NDR-4001 23,2~, 24

KORES = 1

GO TO 26

••• OR SLOWER BINARY TAPES 3,4,8

KORES

=

0

CONTINUE

NEW PAGE, PRINT CASE IDENTIFICATION

WRITE OUTPUT TAPE 6, 25, KASE, KCYCLE, IDGAS, GAS, VMO, DMC, THAWDG, DMTOL, MAXIT,

GAMO, ALFO, TO, S~PO, RHOO

FORMAT ( lHI, 6X,4HCAS~, 4X,6HKCYCLE, 5X,5HIDGAS, 7X,3HGAS I

IX, 3110., 4X,A6 I

lHO, 8X,2HMD, 7X,3HDMC, 4X,6HTHAWDG, 5X,5HDMTOL, 5X,5HMAXITI

IX, 2FIO.6, FIO.4, FIO.6, 1101

1H0, 6X,4YGAMO, 6X,4HALFO, 8X,2HT~, 13X,2Hpn, 11X,4HRHOOI

IX, 3FIO.6, 2E15.6 I

WRITE OUTPUT TAPE 6, 30, DR, NDR, NOTALK, RBl, RB2,

NCYCB1, NCYCB2, NCYCB3,

FMCl, FMC2, NFCYC1, NFCYC2, NFCYC3,

, NCURVl, NCURV2, NCYC1, NCYC2, NCYC3,

ICW

FORMAT ( lHO, 13X,2HDR, 7X,3HNDR, 4X,6HNOTALK I

IX, EI5.5, 2110 I

lHO, I?X,3HRBI, I2X,3HRR2 I IX, 2E15.5 I

C C C C C C C C F TABLE IV - 2

6 1HO, 4X,6HNCURV1, 4X,6HNCURV21 IX, 2110 I

7 1HO, 5X,5HNCYC1, 5X,5HNCYC2, 5X,5HNCYC31 IX, 31101

8 1HO, 7X,3HICW / IX, 110 ) 40 42 50 52 55 170 175 180 190 200 1 2 ~OO 1 305 306 1 2 3 308 IF (KOR~S) 50,5D,40

WRITE OUTPUT TAPE 6, 42

FO~MAT (lHO, 25HPROGRM~ USING COPES ONLY. )

L

= 1

GO TO 55

WRITE OUTPUT TAPE 6,· S2

FORMAT ( lHO, 3~HPROGRAM USING BINARY TAPES 3.4,8. l

REWIND '3 INITIALIZE MCWOK

=

0 NREC

=

0 KFAZE5 = 0 NCURV = 0 ICWPR = 0 KARTAL

=

0 KWALL

=

0 KWALLI = 0 KPRINT

=

KCYCLE-1 KPTITL = 15 1 = 0 VMC

=

VMO GAMC = GA~O ALFC = ALFO

VMC, GAMC MUST BOTH RF GR~ATER THAN l.O.

OTHERWISF STOD CALCS, PRINt COMMENT h~~ ~O ~EXT CASC. IF (VMC-l.0) 170, 170. 180

WRITF OUTPUT TAPE 6. 175, VMC, GA;\K

FORMAT ( lHO, 8X,2HMC. 5X,2HOR, 5x, 6X.4HGAMC,

5X, 42HNOT GRFATER THAN 1.0, HFNC~ ~o NEXT CASE. /

IX, FIO.5. 12X, FIO.I) l

GO TO 15

IF (GAMC-1.0) 170. 170. 190 THACMC = THETA ( VMC. GAMC ) D = SQRTF ( (GAMC+l.Oll(GAMC-1.01 THAW = THAWDG*3.141593/180.0 E = THACMO + THAW

D. E ARE TO AF PASSED ON AS IMPLICIT ARGUMENTS FOR tHF SUSRTN FXFIX

BY MEANS OF COMMON

THACMC = THFTA ( VMC, GA~C ) THAC = THACMC - THACMO

THACDG

=

THAC*lBO.0/'3.141593 IF (THA,C-THA.W) 300. 3l"), 320

AS SOON AS THAC IS LESS THAN THAW, APPLY

NEWTON ITERATION TO GET FINAL MACH NLJ'v1BFR ~1CW. FXF1X

CALL NITER ( VMCP, VMC, VMCW, DMTOL, MAXIT, KEY. ITERS. FXF1X )

IF (KEY) 305, 305, 308

WRITE OUTPUT TAPE 6, 306, VMCP, VMC, KFY. ITERS FORMAT ( 1HO, 33HNEWTON ITFRATION DID NOT CONVERGE/

IX, 33HHENCE ~CW NOT FOUND, DO NEXT CASE/

IX, 4X,6HMC-DMC, 8X,2HMC, 7X,3HKEY, 5X,5HITERSI IX. 2F10.6, 2110 )

GO TO 15 VMC

=

VMCW

THACDG

=

THAWDG MCWOK

=

1

GO TO 320

TABLE IV - 3

C TF THAC = THAW, CURRENT MACH NUMBER = FJNAL VALUE MCW

C

c

C C C ":310 M(\.<JOK = 1 '320 334 335 GO TO 320

THAC GREATER THAN THAW, CONTINUE CALCS

ANG~C = ANGM CVMCl

ANGMCD

=

ANGMC*180.0/3.141593 DELTC2 = ANG~C + THAC

DELC2D = DELTC2*180.0/3.141593 TOPI = 1.0 + O.5*CGAMC-l.Ol*VMO*VMO DENI

=

1.0 + O.5*(GAMC-l.Ol*VMC*VMC TCTO = TOPI/DENI TC

=

T(TO*TO ~XPI = GA~(/(GAMC-l.n) SMPCP0

=

(TCT01**EX DI $MPC = SMDCPO*SMPO EXP2 = 1.0/(~AMC-l.0) ROC~OO = (TCT01**~XP2 RHOC = ROCQOO*RHOO

SMAC = SMA C GAMC, ALFC, TC 1 UC = SMAC*VMC

GC

=

G ( VMC, GAMC, SMPC 1 VKC = VK ( TC, V, C, IDGAS ALFEC = ALFE ( VKC, SMPC 1

WC = W ( RHOC, ALFC, ALFEC, IDGAS, TC, CTl, CT2 )

SEE TF GAS IS OXYGF~ OR ARGON IF (IDGAS-21 335, ":334,334 RCTC = 2.5 + 1.0/TC AC

=

2.5*(1.0 + ALFC) GO TO 338 CONTINUE BCTC

=

1.5 + 1.0/TC AC

=

3.5 + 1.5*ALFC 338 CONTINUE 340

FCV

=

FV C WC, UC. VMC. BCTC, AC, ALFC CORES ONLY. OR BINARY TAPE 3

IF CKORES) 350,350, 340

TO STORE 15 CORNER VALUES, INTO CORE3 CORE3CL) • VMC CORE3CL+1)

=

ANGMC CORE3CL+ZJ

=

ANGMCD CORE3CL+3) z THAC CORE3CL+4)

=

THACDG CORE3CL+5) • TC CORE3CL+6) • DELTC2 CORE3CL+1)

=

DELC2D CORE3CL+8)

=

SMPC CORE3(L+9)

=

RHOC CORE3CL+10)

=

UC CORE3(L+11) c GC CORE3CL+12) • ALFEC CORE3CL+13) c WC CORE3CL+14)

=

FCV

PREPARE FOR NEXT RECORD

L = L + 15 GO TO 360

1 C 360 400 410 415 1 2 420 430 1 435 1 C 450 C C 500 C C C C 600 610 1 C C C 650 620 C C C 651 653 654 1 2 3 4 658 C C 660 TABLE IV - 4 WRITE TAPE 3,VMC,ANGMC,ANGMCD,THAC,THACDG,TC,DELTC2,DELC2D,SMPC,

RHOC,UC,GC,ALFEC,WC,FCV

INCREASE BY 1 THE NUMBER OF RECORDS SAVED 50 FAR NREC

=

NREC + 1

KPRINT

=

KPRINT + 1

IF (KPRINT-KCYCLEI 450, 400, 400 IF (KPTITL-15) 420. 410, 410 WRITE OUTPUT TAPE 6. 415

FORMAT ( 1HO, 6X,2HMC, 2X,6HTHACDG, 8X,2HTC, 2X,6HDELC2D,

11X,2HPC, 9X,4HRHOC, 8X,2HUC, llX,2HGC. 3X,5HALFfC, 11X,2HWC, llX,2HFC )

KPTITL

=

0 GO TO 430

KPTITL

=

KPTITL + 1

WRITE OUTPUT TAPE 6, 435, VMC,THACDG,TC,DELC2D,SMPC, RHOC,UC,GC.ALFEC,WC,FCV

FORMAT ( 1HO, F8.5, F8.3, FI0.5, F8.3, 2E13.5, FIO.S, E13.5, F8.5, 2E13.5 I

KPRINT :a

°

HAS LAST VALUE OF MACH NUMBER BEEN DONE IF (MCWOK) 500, 500, 600

IF NOT, SAVE CURRENT MC, THAC FOR POSSIBLE INTERPOLATION AND INCREASE MC BY DMC

VMCP

=

VMC THACP

=

THAC VMC

=

VMC + DMC AND GO ON WITH CALCS GO TO 200

BUT IF ALL MACH NUMBERS MC HAVE BEEN DONE FROM MO TO MCW, IN STEPS OF DMC,

SEE IF LAST LINE NEEDS TO BE PRINTED IF (KPRINT) 650, 650, 610

WRITE OUTPUT TAPE 6, 435, VMC,THACDG,TC,DELC2D,SMPC, RHOC,UC,GC,ALFEC,WC,FCV

PHASE 3 COMPLETED.

KNOWING MCW, FIND AND PRINT MCI. MC2 FOR SKIPPING OR PRINTING PHASE 4 CURVES VMC1

=

VMO + FMCl*(VMCW - VMOI

VMC2

=

VMO + FMC2*(VMCW - VMO)

WRITE OUTPUT TAPE 6, 620, VMCI, VMC2

FORMAT ( lHO, 4X,4HVMCI, 4X,4HVMC2/ IX, 2F8.5 I

WHEN CORES ONLY ARE USED, SEE IF NREC AT END OF PHASE 3 IS TOO LARGE FOR CORE3(45161. IF SO REPEAT PHASE 3 USING BINARY TAPES 3,4,8 AND THEN DO PHASES 4,5.

IF (KORESI 658,658, 651 IF (NREC-300) 658, 653,653 WRITE OUTPUT TAPE 6, 654

FORMAT ( IH1, 32HCAPACITY OF CORE3f45161 EXCEEDED / IX, 43HWITH NREC

=

OR GREATER THAN 300 IN PHASE 3, / IX, 38HHENCE CANNOT CONTINUE WITH PHASES 4,5. /

IX, 44HTHUS REPEAT PHASE 3 USING BINARY TAPES 3,4,8 / lX, 23HAND THEN DO PHASES 4,5. I

GO TO 24 CONTINUE

BEGIN PHASE 4

CORES ONLY, OR BINARY TAPE 3 IF (KORESI 670,670, 660

L

=

1

ANGMlD = CORE3IL+21 THAI = CORE3IL+31 THAlD = CORE3IL+41 Tl = CORE3 I L+",l DELl = CORE3IL+61 DELlD = CORE3IL+71 SMPI = CORE3IL+8) RHOI

₌

CORE3IL+91 UI

₌

CORE3IL+lO) Gl

₌

CORf3 I L+ll) ALFEI

₌

CORE3IL+121 Wl = CORE3IL+13) Fl

₌

CORE3IL+14) LI NK32 = 1 GO TO 3200 670 REWIND 3 REWIND 4

READ TAPE 3, VMl, ANGMl, ANGMI D, THAI, THAlD, 1 Tl, DELlt DELlD, SMPl, RHOl, UI, Gl, ALFEI, Wl,

READ TAPE 3, VM2, ANGM2, ANGM2D, THA2, THA2D, 1 T2, DEL2, DEL2D, SMP2, RH02, U2, G2, ALFE2, W2, C PRINT CORNER VALUES

3201 MONITR = 1 LI NKC

=

1 GO TO 8000 8001 CONTINUE

C FIND X,y COORDINATES ALONG INITIAL C SET CORNER X,Y

=

0 IN PHASE 4

XC

=

0

YC

=

0

XOI

=

DR*COSFIANGMl) VOl

=

DR*SINFIANGMl) C SET XOI

=

YOl = 0

XO I

=

0

YOI = 0

( ESTABLISH CURRENT XOI, YOI XOI

=

XOI + XOl

YOI

=

YOI + VOl Xl

=

XOI

Yl

=

YOI

DELTI

=

-ANGMI + THAI X2

=

0

Y2

=

0

DELT2

=

ANGM2 + THA2 C FIND X3, Y3, R3, THAYX3

\oJAVE FRONT

CALL X3Y3 I Xl, Yl, X2, Y2, DELTI, DELT2, XC, YC, 1 X3, Y3, R3, THAYX3 )

C FIND DXI

=

P2-P3. DETA

=

PI-P3

CALL DXIETA I Xl. Vl. X2, Y2, X3. Y3.

1 DXI, DE TA )

C FIND SMP3, THA3. THA3D

TABLE IV - 5

Fl F2

CALL P3THA3 I Gl, G2. SMPI, SMP2, THAI, THA2, FI, F2, DXI, DETA, 1 SMP3, THA3, THA3D I

( FIND X4, Y4. DX

CALL X4Y4DX I Xl, YI. X2. Y2. X3. Y3, THA3, I X4, Y 4, DX )

ALFI

=

ALFC ALF2

=

ALFC

TABLE IV - 6 1 ALF4, W4, U4 )

( FIND RH04, SMP4, T4

(ALL RHOPT4 ( RHOl, RH02, SMP1, SMP2, Tl, T2, DX, 1 RH04, SMP4, T4 ,

C FIND DS, ALF3

(ALL DSALF3 ( X3, Y3' X4, Y4, ALF4, W4, U4, 1 DS, ALF3 )

( FIND A4, B4

(ALL A4B4 ( ALF4, T4, V, IDGAS, 1 A4, B4 ) ( FIND U3, T3, RH03 ( ( ( ( ( 5001 680 ( 1 1 690 I' 2 ( ( 3401 ( ( 1001 ( 2001 ( 5002 700 34()2

(ALL UTRH03 ( U4, T4, RH04, SMP4, ALF4, A4, 84, S~P3, ALF3, U3, T3, RH03 ,

FIND GAM3, SMA3, VM3, ANGM3, ANGM3D, A3, 83/T3, 83

(ALL GAMUAB ( V, T3, ALF3, U3, KEYSMA, KEYVM, IDGAS, GAMO, GAM3, SMA3, VM3, ANGM3, ANGM3D, A3, 83T3, 83 )

ALSO FIND DELT2D

DELT2D = ANGM3D + THA3D FIND G3, VK3, ALFE3, W3, F3 G3

=

G ( VM3, GAM3, SMP3 , VK3

=

VK ( T3, V, (, IDGAS ) ALFE3 = ALFE ( VK3, SMP3 ,

W3 = W ( RH03, ALF3, ALFE3, IDGAS, T3, CTl, CT2 F3

=

FV ( W3, U3, VM3, 83T3, A3, ALF3 )

PRINT TITLES, ANSWERS OF FIRST POINT P3 ON (1 KPRINT

=

N(Y(B1 - 1

KPTITL

=

15 LINKCY

=

1 GO TO 5000

(ORES ONLY, OR BINARY TAPE 4 IF (KORES) 690,690, 680

J = 1

LINK34 = 1 GO TO 3450

CURRENT P3 vALUES NOW KNOWN, SAVE THESE 17 VALUES ON 9INARY TAPF 4 WRITE TAP~ 4, X3, Y3, ANGM3, ANGM3D, THA3, THA3D,

SMP3, ALF3, U3, T3, RH03, SMA3, VM3, G3, ALFE3, W3, F3

HENCE FIRST POINT P3 ON Cl HAS BEEN FOUND. FIND THE REMAINIG NDR-1 POINTS P3 ON Cl DO 30 0 1 J = 2, N D R

ESTABLISH (URRENT XOI, YOI XOI = XOI + X01

YOI = YOI + Y01 Xl

=

XOI

Yl = YOI

DELT1

=

-ANGM1 + THAI

Ta GET NEXT P3 VALUES, SET NEW P2 VALUES

=

LAST P3 VALUES LI NK

=

1

GO TO 1000 GO TO 2000

PRINT ANSWERS OF NEXT POINT P3 ON Cl LINKCY

=

2

GO TO 5000

CORES ONLY, OR BINARY TAPE 4 IF (KORES) 710,710, 700 LI NK34

=

2

GO Ta 3400 GO TO 3001

710 3001 C C C C C 15000 15010 C C C C C 2468 120 4101 130 4001 C C 140 150 1 C 3202 C C 9001 8642 8002 C 2002 C 19Q9 100nO NTAPE

=

4 GO TO 3000 CONTINUE RESET I

=

0 1=0

HENCE ALL NOR POINTS P3 ON Cl HAVE BEEN FOUNO. CHARACT. CO. Cl NOW OONE

HENCE SET NKARAC

=

2 NKARAC

=

2

FOR CHECKOUT. DO ONLY A FEW CHARACT.

=

NOT ALK tF (NOTALK) 15010. 15010. 15000

tF (NKARAC - NOTALK) 15010. 3600. 3600 CONTINUE

DO THE CHARACT. C2. C4. C6. C8 • ••• IN SUCCESSION READING FROM TAPE 4. WRITING ON TAPE 8

CORES ONLV. OR BINARY TAPE 4 IF (KORES) 130.130. 120 J

=

1 L1NK41

=

1 GO TO 4150 GO TO 4001 REWIND 4 REWINO 8

FIND FIRST POINT P3 ON C2 THUS.

tNSERT FIRST PI VALUE ON Cl (FROM TAPE 4) LINK

=

1

NTAPE

=

4 GO TO 4000

DELT1

=

-ANGM1 + THAI

CORES ONLV. OR BINARY TAPE 3 IF (KORES) 150.150. 740 L1NK32

=

2

TABLE IV - 7

GO TO 3200

INSERT FtRST P2 VALUE ON C2 (FROM TAPE 3. MO+2.0MC VALUES) READ TAPE 3. VM2. ANGM2. ANGM20. THA2. THA20.

T2. OEL2. DEL2D. SMP2. RH02. U2. G2. ALFE2. W2. F2 PRINT OR SKIP CZ. C4. C6 ••••

LINKAR

=

1 GO TO 9000

t F (KEVPRT) 800Z. 8002. 8642 CONTI NUE

PRINT CORNER VALUES MONITR

=

NKARAC LINKC = Z GO TO 8000 CONTI NUE ALFZ

=

ALFC X2

=

0 Y2

=

0

DELTZ

=

ANGM2 + THA2

FINO FIRST POINT P3 ON C2 LINK

=

Z

GO TO 2000

PRINT OR SKIP C2. C4. C6 •••• IF (KEVPRT) 5003. 5003. 1999

PRINT TITLES. ANSWERS OF FIRST POINT P3 ON C2

IF (R3 - RBl! 10000. 10000. 10010

10020 10030 10040 ( 5003 160 3801 ( ( ( 110 3002 ( ( 180 4102 190 4002 ( ( 1002 2003 ( 4501 ( 2999 ( 5004 800 3802 ( 810 3003 ( ( 4601 ( 4101 ( ( ( KPRINT

=

N(Y(B2 - 1 GO TO 10040 KPRINT

=

N(Y(B3 - 1 KPTITL

=

15 LINK(Y

=

3 GO TO 5000

(ORES ONLY, OR BINARY TAPE 8 IF (KORESI 710,170, 760 LlNK38

=

1

GO TO 3850 GO TO 3002

SAVE 11 VALUES AT P3 (ONTU TAPE 81 NTAPE = 8

GO TO 3000

HEN(E FIRST POINT P3 ON (2 HAS BEEN FOUND.

FIND THE REMAINING NDR-1 POINTS P3 ON (2 DO 3003 I

=

2, NOR

(ORES ONLY, OR BINARY TAPE 4 IF (KORESI 7909790, 780

LI NK41

=

2 GO TO 4100 GO TO 4002 LINK = 2

INSERT NEXT PI VALUES ON (1 (FROM TAPe 4) NTAPE = 4

GO TO 4000

DELT1 = -ANGM1 + THA1

ON C2 SET NEW P2 = OLO P3 VALUES LINK = 2

GO TO 1000

FINO THE NEW POINT P3 ON (2 LINK = '3

GO TO 2000

IF (MOREl 4601,4601, 4501 PRINT OR SKIP (2, C4, (6, •••

IF (KEYPRT) 5004,5004, 2999

PRINT ANSWERS OF NEXT POINT P3 ON C2 LINKCY = 4

GO TO 5000

(ORES ONLY, OR BINARY TAPE 8 IF (KORESI 810,810, 800 LlNK38 = 2

GO TO 3850 GO TO 3003

SAVE 11 VALUES AT P3 (ONTO TAPE 81 NTAPE

=

8

GO TO 3000 (ONTINUE NORMAL EXIT GO TO 4701

ERROR EXIT, REDUCE NOR TO I-I NOR = NOR1

RESET I

=

0 I

=

0

HEN(E ALL NOR POINTS P3 ON C2 HAVE BEEN FOUNO. SEE IF ANY MORE CHARACT. TO 00 IN PHASE 4 NKARAC

=

NKARAC + 1

FOR CHE(KOUT, DO ONLY A FEW CHARACT.

=

NOTALK IF (NOTALK) 16010, 16010, 16000

16010 C 25000 C C C C C C C C 3':;79 820 8101 830 C C 4003 C 840 C 850 C 3203 9002 9753 C 8003 C C 2004 C 30QQ 20000 1 CONTINUE IF (NKARAC - NREC) 3579, 25000, 25000

IF NOT, 5TART PHASE 5 = CWI, CW2, CW3, •••

KFAZE5

=

5

ITAPE

=

8

JTAPE

=

4

FOR CORES ONLY, SET IJ

=

2

TO READ FROM CORE8, WRITE INTO CORE4

IJ

=

2

WJTH CWl

=

ODD CHARACT.

CW WAS EVEN, READ FROM 4, WROTE ON 8 GO TO 26000

IF SO, DO THE NEXT CHARAC. AMÓNG C3, C5, C7, C9, DO THE CHARAC. C3. C5, (7, C9, ••• IN SUCCESSION

READING FROM TAPE 8, WRITING ON TAPE 4

CORES ONLY, OR BINARY TAPE 8

IF (KORESI 830,830, 820 J

=

1 LINK81

=

1 GO TO 8150 GO TO 4003 REWIND 8 REWIND 4

FIND FIRST POINT P3 ON C3 THUS.

INSERT FIRST P1 VALUE ON C2 (FROM TAPE 8)

LI NK

=

3

GO TO 4000

DELTI

=

-ANGMI + THAI

CORES.ONLY, OR BINARY TAPE 3

IF (KORESI 850,850, 840

LlNK32

=

3

GO TO 3200

TABLE IV - 9

•••

INSERT FIRST P2 VALUE ON C3 (FROM TAPE 3, MO+3.DMC VALUESI READ TAPE 3, VM2, ANGM2, ANGM2D, THA2. THA2D,

T2, DEL2, DEL2D, SMP2, RH02, U2,G2, ALFE2, W2, F2

PRINT OR SKI~ C3. C5, C7, •••

LI NKAR

=

2

GO TO 9000

IF (KEYPRT) 8003. 8003, 9753 CONTINUE

PRINT CORNER VALUES

MONITR

=

NKARAC LI NKC

=

3 GO TO 8000 CONTINUE ALF2

=

ALfC X2

=

Û Y2

=

0

DELT2 = ANGM2 + THA2

FIND FIRST P3 POINr ON C3 LI NK = 4

GO TO 2000

PRINT OR SKIP C3, C5. C7, •••

IF (KEYPRT) 50U?, 5005, 3999

PRINT TITLES, AN$WERS OF FIRST POINT P3 ON C3

IF (R3 - RBl) 20000, 20000, 20010

KPRINT = NCYCBl - 1

20030 20040 C 5005 860 3403 C 810 C C 3004 C C 880 8102 890 4004 C C 1003 2005 C 4502 C 4999 C 5006 910 3404 C 920 3005 C C 4602 C 4702 C C C 17000 GO TO 20040 KPRINT

=

NCYCB3 - 1 KPTITL

=

15 LINKCY

=

5 GO TO 5000

CORES ONLY. OR BINARY TAPE 4

IF (KORES) 810,810, 860

LINK34

=

3

GO TO 3450 GO TO 3004

SAVE 11 VALUES AT P3 ONTO TAPE 4

NTAPE

=

4

GO TO 3000

HENCE FIRST POINT P3 ON C3 HAS BEEN FOUND. FIND THE REMAINIG NDR-I POINTS P3 ON C3

DO 3005 I

=

2, NDR

CORES ONLY, OR BINARY TAPE 8

IF (KORESI 890,890, 880

LINK8I

=

2

GO TO 8100 GO TO 4004

LINK

=

4

INSERT THE NEXT PI VALUE ON C2 (FROM TAPE 8 )

NTAPE

=

8

GO TO 4000

DELT1

=

-ANGM1 + THAI

ON C3 SET NEW P2

=

OLD P~

LI NK

=

3

GO TO 1000

FIND NEW POINT P3 ON C3

LI NK

=

5

GO TO 2000

IF (MORE) 4602,4602, 4502 PRINT OR SKIP C3, C5, C7, •••

IF (KEYPRTl 5U06,5006, 4999

PRINT ANSWERS OF NEXT POINT P3 ON C3

LINKCY

=

6

GO TO 5000

CORES ONLY, OR BINARY TAPE 4

IF (KORES) 920,920, 910

LINK34

=

4

GO TO 3450 GO TO 3005

SAVE 11 VALUES AT P3 ONTO TAPE 4

NTAPE

=

4

GO TO 3000 CONTINUE NORMAL EX IT GO TO 4702

ERROR EXIT, REDUCE NDR TO 1-1

NDR

=

NDR1

RESET I

=

0

I

=

0

HENCE ALL NDR POINTS P3 ON C3 HAVE BEEN FOUND.

SEE IF ANY MORE CHARACT. TO DO IN PHASE 4

NKARAC = NKARAC + 1

FOR CHECKOUT, DO ONLY A FEW CHARACT. = NOTALK

IF (NOTALK) 17010, 17010, 17000

IF (NKARAC - NOTALK) 17010, 3600, 3600