NASA Contractor Report 3369 

Hyperbolic/Parabolic Developmen 

for the GIM-STAR Code 

L. W. Spradley, J. F. Stalnaker, 

and A. W. Ratliff 

CONTRACTS NASl-15783 and NASH-15795 

DECEMBER 1980 

NASA 

CR 

3369 

c. 1 
TECH LIBRARY KAFB, NM 

NASA Contractor Report 3369 

Hyperbolic/Parabolic Development 

for the GIM-STAR Code 

L. W. Spradley, J. F. Stalnaker, 

and A. W. Ratliff 

Lockheed Missiles G Space Company, Inc. 

Huntsville, Alabama 

Prepared for 

Langley Research Center 

under Contracts NASl-15783 and NASl-15795 

NASA 

National Aeronautics 

and Space Administration 

Scientific and Technical 

Information Branch 

1980 
. 

: - 

,. ‘. 
CONTENTS 

Section Page 

1 INTRODUCTION AND SUMMARY ................ 1 

2 CALCULATION OF TWO-DIMENSIONAL INLET 

FLOWS WITH SPILLAGE .................... 7 

2.1 Introduction ......................... 7 

2.2 Development of the Solution ................ 9 

2.3 Results and Discussion ................... 11 

3 INVESTIGATION OF LINEARIZED BLOCK IMPLICIT 

METHODS FOR THE GIM CODE ................ 35 

3.1 Introduction ......................... 35 

3.2 One-Dimensional Unsteady Development ......... 36 

3.3 One-Dimensional Linearized Block Implicit Results ... 45 

4 DEVELOPMENT OF GIM/STAR SPATIAL MARCHING 

ALGORITHMS .......................... 51 

4.1 Introduction ......................... 51 

4.2 Classical Parabolic Approaches .............. 52 

4.3 The Quasi-Parabolic Idea ................. 58 

4.4 Results of Computation ................... 66 

5 REFERENCES .......................... 75 

A ppendixe s 

A Use of the GIM SE-2 Code ................... 78 

B Three-Dimensional LB1 Schemes for the Navier-Stokes 

Equations. ............................ 104 

C Multilinear Interpolants for GIM Marching Methods ...... 136 

D Vectorized Block Tridiagonal Equation Solver for the 

GIM/STAR Code. . . . . . . . . . . . . . . . . . . . . . . . . 152 

. . . 111 
1. INTRODUCTION AND SUMMARY 

The General Interpolants Method (GIM) code was developed to analyze 

complex flow fields which defy solution by simple methods. The code uses 

numerical difference techniques to solve the full three-dimensional timeaveraged elliptic Navier-Stokes equations in arbitrary geometric domains. 

The equations are cast in strong Conservation law form and written in an 

orthogonal Cartesian coordinate system. Included are a continuity equation 

for global mass conservation, three components of momentum conservation, 

total energy conservation and an equation for conservation of individual species 

of a binary gas. Pressure is related to the conservation variables through 

the ideal gas law for a binary mixture. A generalized geometry package is 

used to model the flow domain, generate the numerical grid of discrete points 

and to compute the local transformation metrics. Computation is done in 

physical space by explicit finite-difference operators. The GIM approach 

essentially combines the finite element geometric point of departure with 

finite difference explicit computation analogs. This provides a capability 

which takes advantage of the geometric flexibility of an element description 

and the superior computation speed of difference representations. 

The numerical analogs of the differential equations are derived by 

representing each flow variable with general interpolation functions. The 

point of departure then requires that a weighted integral of interpolants be 

zero over the flow domain. By choosing the weight functions to be the interpolants themselves, the GIM formulation produces identically the classical 

implicit finite element discrete equations. These forms are not used in 

the GIM code due to their fully implicit nature and inherent inefficiencies. 

Rather, the weight functions are chosen to be orthogonal to the interpolant 

functions which produces explicit finite difference type discrete analogs. 

By appropriate choice of constants in the weight functions, the GIM becomes analogous to such finite difference schemes as centered, backward, 
forward, windward and multi-step predictor-corrector schemes such as the 

MacCormack method. The GIM analogs, however, are automatically produced 

for arbitrary geometric flow domains and,:h,ence is a-more general point of 

departure and provides greater flexibility in choosing difference schemes. 

I . . ’ _ 

A motivation for developing this code on these. principles was to pro- I 

vide an analytical ,tool,tihich is more user ,orientedL than the basic research 

tools which exist.. A. fully production-line code to solve the complex NavierStokes equations .does not exist today. In dev.eloping the GIM code, an atte.mpt 

was made to bridg.e the.gap somewhat between the pure research codes and 

the ,ultimate production tool. The code was originally develpped for a CDC 

7600 computer system. It has been subsequently converted to vector FORTRAN 

for the CDC STAR-100 system at NASA Langley Research Center. Reference 1 

provides documentation for the GIM/STAR code designated version SE-l (STARElliptic No. 1). This version of the code has been used to compute a number 

of complex flow fields including nozzle flows for both subsonic and superson.ic 

regimes, and.two and three-dimensional Scramjet exhaust fLow simulations 

(Ref. 2). 

The current contract work involves utilization and extensions of the 

GIM/STAR code. The objectives of the study are to: 

0 Compute flow fields in supersonic inlet configurations 

using the SE- 1 code 

l Upgrade the technical capability of the SE-1 code 

l Develop a hyperbolic and a parabolic version of the GIM/ 

STAR code to supplement the elliptic capability.’ 

This report presents the progress to date on the development and application 

of the GIM/STAR code. 

Section 2 presents the res’ults of an application of the code to a two: 

dimensional super sonic inlet. The calculation was started upstream of the 

compression surface which turns at 25 deg to the-horizontal. ’ The Mach = 5 

2 
freestream flow generates a bow shock off the leading edge of the ramp. The 

calculation involved t&o primary considerations’; (1) determine the &no&t ’ 

of flow captured and the amount spilled into the freestream and (2) c-ompute 

the inlet flow field and predict the shock wave/boundary layer interaction. 

: 

The problem was run in ttio parts with the GIM code on the STAR machine. 

The ingested flow was determjned,(invisc,idly) first and found to be 66% of 

the incoming stream. This agrees well with the numbers for which the simulated inlet was designed. The flowfield distribution at ‘the nozzle’entrance 

was then used to drive the internal flow allowing the performance parameters 

to be determined; The flow angularity produces a shock wave off the cowl 

lip which propagates into the nozzle. The ultimate interaction of this shock 

and the laminar boundary layer on the upper propulsion surface were computed. All shock waves were determined using the “capture” mode of calculation. Section 2 shows the computed solution for the spillage part of the flow 

and for the internal nozzle portion. The separation of the boundary layer due 

to the adverse pressure gradient is clear from the velocity and pressure 

contour plots. Radial distributions of the steady state flow field are given 

and a “time” history of the shock/boundary layer interaction calculation is 

also shown. 

Section 3 of this report describes an investigation of linearized block 

implicit (LBI) finite difference schemes for the GIM code. The current 

explicit MacCormack schemes are relatively efficient for flows with inviscid boundary conditions. In anticipation of other requirements to compute three-dimensional viscous flows, the, necessity of eliminating the explicit 

stability limit becomes apparent. However, the extreme inefficiencies inherent in “fully I1 implicit methods, due to the large band-width matrices, 

make them unrealistic for large three-dimensional viscous flow problems. 

The most promising concept is the linearized block implicit (LBI), or approximate factorization, schemes. These methods retain the Conservation Law 

equation form while *I splitting” the spatial dependence in the manner of the 

AD1 schemes. The resulting matrix bandwidth is- once again small-(usually 3) 

3 
?and is practical to use. The study of LB1 schemes in this work was concentrated on: 

l Stability requirements of the block tridiagonal scheme 

of Beam-Warming (Ref. 3) 

l Accuracy of the LB1 scheme itself and more precisely, 

the accuracy and speed of linear equation solvers for 

vet tor machines 

l Shock wave resolution of LB1 schemes used in a capture 

mode and artificial damping requirements 

l Techniques to vectorize LB1 schemes for use on the 

STAR- 100 machine. 

‘The study was carried out with a one-dimensional code that uses the BeamWarming formulation. 

Results of the LB1 investigation are discussed in detail in Section 3. 

The stability of the scheme was found to be strongly coupled to the accuracy 

of the linear equation solver used and to the artificial damping added to the 

explicit side of the scheme. The “unconditional” stability indicated by the 

theory could not be achieved numerically using centered differences. Schemes 

based on one-sided windward differences did prove to be unconditionally stable. 

The LB1 scheme was shown to be as good as the explicit MacCormack for resolution of shock waves. The overall conclusion of this part of the study is that 

LB1 schemes appear to be very promising for three-dimensional viscous flows 

but they are not as outstanding as the literature indicates. 

The third part of this study reported here is the development of hyperbolic and parabolic methods to supplement the elliptic code. Section 4 describes 

the details of the work on the GIM maching algorithms and the current status 

of the code. The basic idea of the GIM code marching scheme is to combine 

the classical parabolized Navier-Stokes methods with a “quasi-time” relaxation. The term “quasi-parabolic” (QP) will be used to refer to this algorithm 

although the scheme applies equally well to hyperbolic, supersonic inviscid 

flows. The QP algorithm is contrasted to a fully elliptic method in that down- 

‘stream effects cannot be felt upstream and that a full flow domain need not be 

4 
stored for the QP scheme. The QP algorithm is also contrasted to classical 

parabolic methods in that mixed subsonic/supersonic flows do not produce a 

multiple “decode” root and that real-wall no-slip boundaries can be treated 

with the QP algorithm. The equations are the classical parabolized NavierStokes but with a psuedo-time derivative added back to them. The solution is 

known at upstream data planes 1,2,. . . . N-l and the solution is sought at plane 

N with no influence from plane N+l. Time relaxation is used to solve for plane 

N from only the (converged) solution at upstream planes. Backward differences (second order) are used, of course, in the quasi-marching coordinate, 

As the algorithm is formulated, either explicit or linear block implicit time 

relaxation can be incorporated. 

The resulting algorithm then requires much less computer storage than 

a GIM elliptic flow field calculation and does not have the “singularities11 inherent in classical parabolic marching algorithms. The QP scheme has been 

coded and partially checked out on the STAR system. At the time of this 

writing, the GEOMETRY,MATRIX and INTEG modules of the SP-1 GIM code 

(STAR Parabolic, Version 1) have been run successfully for several sample’ 

cases. 

Some details of the current contract work are appended. The most current version of the GIM elliptic code (SE-2) is discussed in Appendix A authored 

by L. W. Spradley. Differences in SE-l and SE-2 are described and reasons for 

the changes explained. New INPUT data sheets for SE-2 are given to replace 

the ones in the “Blue Book” (Ref. 1). This basic guide should be used in conjunction with the Blue Book for inputting the GIM code on STAR. Appendix B 

by Jiirgen Thoenes, contains a derivation and list of the complex linearization 

Jacobians for three-dimensional LB1 schemes. The GIM-Marching code (SP- 1) 

requires a special set of weight/shape functions. These are derived in Appendix 

C, which is authored by John F. Stalnaker. The final item to be covered here is 

a description of the vectorized linear algebraic equation solvers which were 

developed on the STAR system for use with the LB1 schemes. The mathematical 

development and performance of several techniques, both direct and iterative 

are shown in Appendix D, authored by S. J. Robertson. 

5 

I 

2. CALCULATION OF TW-&DIMENSIONAL INLET 

FLOWS WITH SPILLAGE 

2.1 INTRODUCTION 

Figure 2-:l shows the model two’-dimensional supersonic inlet for 

which the flow field was computed using the elliptic GIM/ST,AR code. 

The compression ‘surface makes a sharp 25 deg turn at x= 0. It turns 

50 deg through a circular arc centered about x= 5 into the 25 deg expansion surface. The expression surface and the lower cowl from the 

nozzle. The freestream flight conditions are also shown in Fig. 2- 1. 

All flor variables are made dimensionless with the freestream quantities. 

For inlets with fixed geometry it is important to know the amount of 

flow captured by the inlet.and the amount that spills into the freestream. Thus, 

special emphasis was placed on calculating the mass flow rate at the inlet 

throat (x = 5). The model inlet was designed inviscidly to capture 66.6% of 

the incident flow. 

It is felt that a brief history of the development and an outline of the 

pitfalls incurred in obtaining the final solution would be of benefit to future 

users of the GIM/STAR code. This discussion appears in Section 2.2. A 

complete analysis of the final solution is given in Section 2.3. These discussions are divided into two parts: (1) the external flow field below the 

compression surface and including the freestream flow which spills below 

the cowl, and (2) the internal (nozzle) flow field. 

Use of trade names or names of manufacturers in 

this report does not constitute an official endorsement 

of such products of manufacturers 

istration. 

_~- ___ -- .ym either expressed or 1 

implied, by the National Aeronautics and Space Admin- 

-- 
0.67 J 

Freestream 

y = 1.4 

P= 1 

p = 1.4 

M = 5.0 

t 

1 

.l 

I I I I I I I I I I 

1 2 8 9 10 3 4 5 6 7 

x (axial) 

Fig. 2-l - The Model Two-Dimensional Inlet 
2.2 DEVELOPMENT OF THE SOLUTION 

2.2.1 External Flow Field 

To limit the problem size, the computational grid was constructed 

originally with the input boundary lying along the 35.7 degree bow shock line., 

This resulted in computational difficulties with the grid points along a horizontal line from the leading edge of the cowl to the shock line. In order to 

wrap the grid around the cowl, a discontinuity in the shapes of the elements 

arose along this line. Relatively uniform rectangular elements were mated 

to severely skewed elements. It is believed that the computational problems 

arose from the finite difference analogs generated along this line of nodes. 

These improper influences were caused by either sharp discontinuities in 

the transformation metrics or inadvertent extrapolation in the transformations. 

As a result, the post-shock grid was abandoned and it was decided that 

the geometry should be constructed to allow the bow shock to be captured. The 

analogs for this grid were thoroughly examined using a coarsely spaced version 

of the final computational mesh (shown in Fig. 2-2). A “double-valued” node (i.e., 

two nodes at the same spatial location) was used to allow the proper splitting of 

the flow at the cowl lip, Due to the small shock angle at the bow which would not 

permit a sufficient number of nodes between the shock and the surface, the first 

eight nodes along the compression surface were held fixed at the inviscid postshock conditions. This eliminated numerical disturbances which were otherwise 

generated at the bow and propagated downstream leading to instabilities. 

2.2.2 Internal Flow Field 

As originally modeled, the upper body of the inlet had a sharp 50 deg 

expansion at x = 5. In the initial inviscid analysis of the nozzle it was found 

that the flow overexpanded around this turn leading to pressure undershoot 

and instability. The sharp turn was rounded to alleviate this overexpansion. 

However, subsequent analysis revealed the problem to be excessive damping on the continuity equation. This resulted in an artificial dissipation 

9 
.Fig.2-2 - Computational Grid fsr External Flow in the Two-Dimensional Inlet : 
of mass away from the wall. Reduction of this damping allowed successful 

computation of the expansion; however, the rounded surface remained. The 

primary difficulty with the nozzle calculation was with the inviscid treatment 

of the expansion surface. This was first indicated by the failure of the inviscid SEAGULL code (Ref.4) to converge in the nozzle. Imposing a viscous 

boundary layer on the upper wall allowed the GIM code to develop a strong 

shock-boundary layer interaction which made evident the fallacy in the inviscid treatment. 

In arriving at the final solution it has become increasingly clear that 

solutions with the GIM code are strongly dependent on two factors: (1) the 

structure of the computational mesh, and (2) proper modeling of the physics 

of the problem. 

2.3 RESULTS AND DISCUSSION 

2.3.1 External Flow Field 

The 3557 node computational grid for the external flow field is shown 

in Fig. 2-2. The solid boundaries were treated inviscidly. The USERIP option 

in the GIM/STAR code was used to initialize the flow field in order to lay in the 

bow shock as closely as possible to the inviscid 35.7 degree line. The solution 

converged to steady state in 900 iterations. The integrated mass flow rate 

indicated that the inlet captured 66.5% of the incident mass flow which compared almost exactly to the theoretical value. Figures 2-3 through 2-5 show 

the velocity vectors, pressure and Mach number contours for the complete 

flow field. Figure 2-6 shows a comparison of the mass flow rate (m=pu) 

across the inlet plane as calculated by the GIM/STAR code to that calculated 

by the inviscid SEAGULL code (Ref. 4) for a similar inlet with the same imposed capture rate. The agreement is excellent with the only apparent differences resulting from the different treatment of shocks in the two codes. For 

computational economy the deteched shock effects at the cowl lip were not 

treated here, Rather, after 100 iterations the values of the flow variables at 

the lip node were held fixed at attached post-shock conditions determined 

11 
Fig. 2-3 - Velocity Vectors for the Two-Dimensional Inlet (Maximum Velocity = 5.0). 
Fig.2-4 - Pressure Contours for the Two-Dimensional Inlet. 
ID Mach 

1 0.21 

2 0.64 

3 1.00 

4 1.49 

I? 

1.91 

2.34 

7 2.76 

8 3.19 

9 3.61 

10 4.04. 

Fig. 2-5 - Mach Number Contours for the Two-Dimensional Inlet. 
-2, 

-2 

* 

d 

0 .d 4-l .d 

0” 

n -2 

?i 

.mo 

a 

k 

z 

-2 

-2 

.6- 

.8- 

. 

l GIM 

0 SEAGULL 

I 

10 

I I 

20 30 

Mass Flow per Unit Area, Pu 

I 

40 

Fig.2-6 - Mass Flow per Unit Area vs Vertical Position Across the Inlet 

Plane - (St a t ion 5) Comparison with SEAGULL Solution (Ref. 4). 
from conditions immediately upstream of the cowl lip, Further, when the 

lip shock was allowed to detach, the blunt body effects did not extend beyond 

one grid width from the lip. 

2.3.2 Internal Flow Field 

As noted in Section 2.2.2, it was necessary to impose a viscous boundary 

layer on the upper wall of the nozzle to obtain the solution. The boundary 

layer profile at the throat was estimated by a quadratic laminar profile with 

a Reynolds number of 1 x 104. The inlet flow variables were input by linear 

interpolation of the external flow results and the remaining nodes exterior 

to the boundary layer were initialized via the USERIP option by an isentropiclike area expansion along streamwise rows of nodes. It was found that allowing the code to develop the shock in the nozzle was preferable to estimating 

the shock position as was done in the external flow calculation. Figure 2-7 

shows the 3000 node computational mesh for the nozzle. The solution converged to steady state in 1200 iterations. The final velocity vectors, pressure and Mach contours are shown in Figs. 2-3 through 2-5. Figures 2-8 

through 2-20 show the time development of the solution from iteration 0 

through 1200. Figures 2-21 through 2-24 show variations of the Mach number 

and pressure in the nozzle compared with the available SEAGULL results. It 

is apparent from these last figures that the boundary layer is artifically too 

thick (a result of the choice of Reynolds number). However, the result obtained provides considerable insight into the physics of the problem as well 

as the reasons behind the failure of the inviscid analysis. 

16 
Fig. 2-7 - Computational Grid for Internal Flow in the 

Two-Dimensional Inlet. 

17 
c- - 

* -e --- c) -cI 

CcI 

Fig.2-8 - Nozzle Velocity Vectors (ITmax = 4.91; No Iterations). 

18 
-- 

c 

c 

- 

-- 

Fig. 2-9 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V 

max 

= 4.92; Iteration 100). 

19 
, 

Fig. 2-10 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V - 4.98; Iteration 200). 

max - 

20 
Fig. 2-11 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V 

max 

= 5.00; Iteration 300). 

21 
Fig.2-12 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V 

max 

= 4.84; Iteration 400). 

22 
Fig.2-13 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V max 7 4.88; Iteration 500). 

23 
Fig.2-14 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V 

max 

= 4.82; Iteration 600). 

24 
Fig.2-15 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V 

max 

= 4.87; Iteration 700). 

25 
11.-- CCL- - - 

----4--o-- - - 

- - 

Fig. 2- 16 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V 

max 

= 4.90; Iteration 800). 

26 
Fig.2-17 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V 

max 

= 4.91; Iteration 900). 

27 
- 

- - 

Fig.Z-18 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V - 4.91; Iteration 1000). max - 

28 
Fig.2-19 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V max 

= 4.91; Iteration 1100). 

.* 

29 
Fig.2-20 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V max 

= 4.91; Iteration 1200)- 

30 
. . . . 

GIM 

SEAGULL 

No SEAGULL Results 

T Upper Wall Available. Failure to 

x = 6.0 Converge at xz 6.4 . 

0 1 2 3 4 0 1 2 3 4 0 1 2 

Mach Number Mach Number Mach Number 

Fig.221 - Mach Number Variations Across the Nozzle. 
Upper Wall 

- /x=5.4 

!LLu 

I 

I 

GIM 

0 l l SEAGULL 

I 

20 

Pressure, P/P 

co 

Fig. 2-22 - Pressure Variations Across the Nozzle. 
xx 6.0 

GIM 

l l l SEAGULL 

-2. 

No SEAGULL Results 

Converge at x Z 6.4 

J 

-3. I 

0 10 20 

Pressure, P/P 

co 

10 

Pressure, P/P 03 

Fig.2-23 - Pressure Variations Across the Nozzle. w 

w 
- GIM 

0 a l SEAGULL 

l.O- 0 SEAGULL fails to 

Converge at this 

0 

0.5- 

0.0, I 

. 

9.. 1. . .~ooop ..I. l 

l 

, , 

5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7. 

Axial Distance from Throat 

0 

F&.2-24 - Pressure Variation on Upper Nozzle Wall. 
‘3. INVESTIGATION OF LINEARIZED BLOCK 

IMPLICIT METHODS FOR THE GIM CODE 

3.1 INTRODUCTION 

Numerical solution of the unsteady Navier-Stokes equations by explicit 

finite difference techniques has a number of disadvantages. The most serious 

one, from a practical engineering viewpoint, is the small time steps which are 

usually required to maintain stability. Computation of boundary layer flows at 

high Reynolds number requires fine grids near solid boundaries, hence very small 

time steps and long computer run times. One apparent cure for these difficulties 

is the use of implicit methods some of which are unconditionally stable for any 

size time step. These schemes are not without problems of their own in terms 

of their practical use. Among the major difficulties are the following: 

1. Implicit finite differences, in general, lead to systems of nonlinear 

algebraic equations when applied to the Navier-Stokes equations. 

These must either be solved directly or linearized in some manner. 

2. Direct linearization, via classical ADI processes, will destroy the 

Conservation Law Form of the Navier-Stokes equations and hence 

shock capture algorithms cannot be used. 

3. Multi-dimensional implicit methods lead to very large systems of 

simultaneous algebraic equations. Even for linear systems, the 

efficient solution is not practical due to large size of the matrix 

coefficients. 

4. Fully implicit methods cannot be programmed for efficient use on 

advanced vectorized machines such as the STAR, ILLIAC, or NASF. 

Numerical treatment of the steady state parabolic form of the NavierStokes equations face many of the same difficulties as the elliptic form. The 

spatial marching step size is constrained by the small grid required to resolve 

boundary layers normal to a solid wall. Marching downstream great distances 

can result in impractically long run times. Implicit finite differences have 

the potential to eliminate the difficulties mentioned above. 

35 
3.2 ONE-DIMENSIONAL UNSTEADY DEVELOPMENT 

The first item to be developed is the formulation of an implicit scheme 

which results in a linear algebraic system yet retains the conservation law 

form of the Navier-Stokes equations. This idea can be explored by considering the equations in one space variable, x, and the time coordinate, t. 

t I 

4tlx 

4x 

Direct linearization is usually done by “lagging” certain of the nonlinear contributions by one time step. This destroys the conservation nature of the 

Navier -Stokes equations . 

The case considered here is an elliptic boundary value problem in space 

and an initial value problem in time. The equations considered are: 

Governing Equations 

(1) 

7 = 

36 
where 

P = mass density 

P= pressure 

p = viscosity parameter 

(2 ps + A) 

t = time coordinate 

v = flow velocity 

6 = total energy 

k = thermal conductivity 

x = space coordinate 

P = (y- 1) p [&- v2/2] ideal gas law 

General Finite Difference Form 

-- 

This analysis will use the “delta” form of the flow variables 

4Un = un+l - un 

AEn = En+’ - En (2) 

fG n+l = 7 - Tn 

where n is the time step index. All data are assumed known at n= 0. Solving 

for AUn then allows the data at level n to be advanced to level n + 1: 

U n+l = un + 4un 

The class of finite-difference schemes considered can be written as follows: 

AUn = +$ s (AU”) + at 1 (U”) + -& 4U”-’ + 0 li-E at 

(0-f -E ) At2 t At3 1 (3) 

The parameters 0, E are used to generate a specific type of scheme. For 

8 = 0, the scheme is fully explicit; 8 > 0 gives an implicit method. If E = 0, 

the schenle requires two data sets of storage at time levels n, n+ 1. If 

E > 0, then three levels are required to be stored, n - 1, II, n+ 1. 

37 
IIf@= l/2 + E, the scheme is second order accurate; and is first order 

otherwise. In this work, we are primarily concerned with 

cl = 1 E = l/2 

which is a second order, implicit, three level scheme. 

Development of the Scheme 

The differential equation (1) is substituted into the general scheme (3) 

to get the following form: 

(-A.? t AT”) t +$ 1 [ & (-En t 7n) t & .4Un-’ 1 (4) 

By approximating the spatial derivatives, a/ax, by finite-differences, we get 

a set of nonlinear algebraic equations. The incremental flow variables, AEn, 

AT’~ are nonlinear functions (1) of the independent vector AUn. 

For our implicit scheme, this would require a simultaneous nonlinear 

algebraic equation solver. The best known methods are iterative ones which 

require long computer runs. 

For this work, we will perform a linearization as follows to obtain a 

sot of linear algebraic equations and use matrix methods for their solution. 

The main idea here is to linearize the algebraic equations, but retain the 

fully conservative nature. 

Expanding E, 7 in a Taylor series, we get 

E ntl 

- U”) t O(At’) 

38 
or 

and 

or 

d-.+- O(At2) 

= An AUn + .0(At2) 

n+l 

T = 7”. + ($ AU” + ($) AU: t O(At2) (5) 

.A? = Pn AUn t Rn AU; t O(At2) 

where 

U X = au/ax 

The expression for A7 can be rewritten in a more convenient form by expanding the x-derivative to get 

A? = (P-Rx) AUn t & (R AU)” t O(At2) 

where 

R X = aR/ax 

This form produces a linear system of equations with the same formal accuracy (At 2 ) as the nonlinear set. It does however, require evaluation of the 

Jacobians 

(6) 

and 

aR/ax 

39 
Putting the Taylor series (5) into the scheme (4) gives the following 

expression: 

-An Aunt (P-Rx) AUn t +-(RAU)n 

I 

t e &(-En+ ~~)t -& AUn-' 

The last two terms on the right hand side of Eq. (7) are all explicit at time 

levels n, n- 1. Denote this by Dn, and write Eq. (7) as follows: 

A$+ 2 (A-PtRx)nAUn-$$RAU)n 

1 

= Dn 

For convenience, let 

h = s B=A-p+R X 

and write Eq. (8) as follows 

Aunt h a2 &(BnAUn) -- 

ax2 

(RnAUn) 1 = Dn 

(7) 

(8) 

(9) 

To see that the form Eq. (9) may be useful, we will now write it for node point 

l’it’ in space and use second order centered finite differences 

af 

I 

fitl - fi-l 

axi= 2Ax t O( x2) 

(10) 

a2f= fitl - 2 fi t fi 1 

ax2 I i 

&2 - + O(h2) 

40 
With these difference expressions, Eq. (9) can be written as follows: 

1 

)I 

= D; (11) 

Combining coefficients of each AU:, AUrtl, AU:-1 terms gives 

h n 

-- &2 Ri,l n Aui-l = Dr 

(12) 

(where I is the 3x 3 identity matrix) 

Boundary values i = 1, and i=K must be treated separately due to the 

centered differences. For now we will let i= 2, 3, . . . k - 1 and worry about 

boundary conditions later. 

The coefficients of the AU terms are 3 x 3 matrices which couple the 

three governing equations at each node point. There is an equation (12) for 

eachnode i=2, 3, . . . k-l. 

To readily see the character of this system of linear algebraic equations, let 

Lfl = 

-h n h n 

ZG Bi-l 

-- 

Ax2 Ri-l 

21~ n My = It- Ax2 Ri (13) 

NE1 = h II h n -- 2Ax Bitl Ax2 Ritl 

41 
The finear algebraic system then has the form ~. 2 . 

:: 

or in matrix notation: 

. 

M2 

L2 

0 

. 

. 

. 

. 

0 

N3 0 0 . . . 0 

M3 N4 0 0 0 

L3 M4 N5 O . 

. 

. 

. . . 

. . . 

- n 

D2 

D; 

Dqn 

. 

. 

G-2 

n 

DK-l 

(14) 

The system (14) will be termed “block tridiagonal.” The individual 

matrices are full 3 x 3 arrays but they are arranged in a tridiagonal manner 

in the full matrix. The block arrangement occurs due to the linearization 

scheme used. This effectively couples the three differential equations at 

each node point. The boundary values for i = 1, and i = K have not been 

treated. This is an additional development item. 

The advantages of a system like Eq. (14) are: 

1. The Conservation Law form has been retained. 

2. Block tridiagonal systems are not much more costly to solve than 

pure tridiagonal systems. 

3. Operations like (14) can be vectorized for use on STAR-like machines. 

Formulation of this scheme requires the analytical evaluation of the Jacobian 

matrices A, P, R. A brief look at these operations now follows. 

42 
Calculation of the Matrices 

The final matrices needed are L, M, N in Eq. (13). These are made up 

of combinations of B, R matrices from Eq. (5). 

Bn = An - P” + R; 

n 

n 

The matrices A, R will have relatively simple elements (as we shall see), 

but the P, Rx matrices will be quite complex. For now we will assume that 

the viscous coefficients are constants; hence we will see that 

P-R = 0 X 

(See Beam-Warming paper, Ref. 3). 

We then need to analytically evaluate A, R, where 

aEi a7i 

A.. = au 33.. = av 

‘J j ‘J X. 

J 

The algebra for these operations is straightforward and is not included here. 

The final results are: 

43 
A = 

1 

R=; -PV 

: 

0 0 

I-1 

1 -(p+)v2-$- w-&v V V V 

0 

Y-1 

YV 

0- 

0 

k 

cv 

where Cv is the (constant) specific heat at constant volume, y is the ratio 

of specific heats and k is the thermal conductivity. 

Summary of Computational Procedure 

1. Set initial data at t = 0 for all nodes i = 1, 2, . . . K. 

2. Form the vectors U, E, 7. 

3. Compute D by explicit differences (Eq. (8)). 

4. Evaluate A, R matrices from Eq. (1.6). 

5. Form the L, M, N matrices from A and R (Eq. (13)). 

6. Modify for boundary values. 

7. Call TRIDAG in the GIM code logic to solve the block 

tirdiagonal system for AUn. 

8. Advance solution vector to (n+l) 

U x-ii-1 

= Un + AUn 

9. Repeat the process to step 2 for a specified number of 

steps or until IAU”l < 6 for convergence to steady state. 

44 
3.3 ONE-DIMENSIONAL LINEARIZED BLOCK IMPLICIT RESULTS 

The procedure outlined in Sections 3.1 and 3.2 was subsequently coded 

and checked out. The equations were modified slightly to handle the problem 

of an expanding duct, quasi one-dimensional, by the inclusion of the area terms. 

This permits the computation of flows other than just the trivial case of constant property flow through a constant area duct. Three cases were considered in order to check out and prove the method. Consider Fig. 3-l where 

the simplest case is when the inflow conditions are fixed at the upstream end 

of the duct. For. completely subsonic flow, elementary considerations indicate 

that the outflow at the downstream end of the duct has a unique solution. Consider for the moment that the flow is controlled entirely by the inflow conditions 

and the out-flow conditions are permitted to develop freely. Of course, it is 

known that physically one could change the back pressure at the downstream 

end and this would affect conditions at the upstream end. But, computationally, 

we specify the inflow conditions and therefore all the flow properties are 

uniquely determined. The same reasoning applies to the case where the flow 

is completely supersonic. In this case there exists the choking effect which 

means that when the back pressure is lowered below the limiting value no 

upstream effect is felt. If however, the back pressure is raised, the situation develops where a normal shock moves into the duct with its positioning 

depending upon the back pressure. Thus for fixed inflow conditions for the 

supersonic case an unique solution depends upon the outflow pressure. Since 

so much is known analytically about this quasi one-dimensional case it was 

deemed a reasonable test with which to evaluate the linearized block implicit 

(LBI) scheme. 

In order that the LB1 scheme could be applied to all three cases, i.e., 

including the strong shock case, pseudo viscous effects were included in the 

original coding in terms of numerical diffusion cancellation (NDC) terms. The 

first case computed was for fully supersonic flow through the expanding duct. 

The LB1 scheme worked well reproducing the analytical results within about 

2’7” over the length of the duct. The case was initially run at a Courant number of one. Subsequent runs were made at larger Courant numbers up to 3 

45 
-* M>l 

Normal 

Fig. 3- 1 - One Dimens ional Duct Configuration 

46 
with good results. Increasing the CFL multiplier further caused rapid deter-ioration of the solution and ultimate destruction of the case (it blowsL up). 

Theoretically the implicit solution should work for very large Courant 

numbers,. Mathematically this is, of c.ourse, true but it ignores the. physics 

of the situation.. -To verify-that the solution was correct the complete der.ivation-was double checked, the coding was rechecked and.nothing was found 

wrong. At first it was thought that the non-dominance of the main diagonal 

might be causing matrix ill-conditioning. The super- and sub-diagonals are 

both proportional to the step size while’the main diagonal remains constant 

(at least for the inviscid.equations). A natural conclusion might then be. 

drawn that, as the step size is increased, non-dominance could occur such 

that the solution of the block matrices loses accuracy thus destroying the 

solution. 

To test out this theory some numerical experiments were carried out. 

First, an unblocked scalar matrix with three diagonals was used. A known 

solution was fed into the matrix reduction scheme and the non-dominance 

factors between the main and other diagonals were increased gradually. The 

case was run on the PDP-11, single precision arithmetic and inaccuracies did 

show up in the sixth place for even a 2 to 1 non-dominance ratio. At 10 6 to 1, 

inaccuracies occur in the first and second places and at lo9 to 1, order of magnitude inaccuracies were produced. Using double precision arithmetic’on the 

PDP-11 or running the case on a CDC 7600 produced no inaccuracies whatsoever. 

Thus it is concluded that scalar matrices manipulated on high precision computing 

equipment have no accuracy problems associated with diagonal non-diminance. 

The same type of numericai experiments were then conducted with the 

block matrices. .. -The CDC ,760O was used in order to eliminate.any inaccuracy due to less precise computing equipment. Non-dominance ratios on the 

order of lo6 to 1 were necessary to generate errors in the fifth and sixth 

place. Since’the suspected non-dominance caused by increasing-the step 

size would only be of order 10, -it is concluded that the reason the case would 

47 
not run at large Courant numbers is due to the problem physics and is not 

related to the mathematics of diagonal non-dominance. 

Subsequent consultations with NASA-Ames personnel (Robert Warming 

and Richard Beam) indicated that they saw no accuracy problems related to 

non-dominant diagonals and they believe the problem with using large Courant 

numbers is due to physically unrealistic propagation of pressure signals which 

then cause oscillatory behavior and eventually a negative pressure. Two 

different solvers were used to eliminate the possibility of an error in the 

coding. The two solvers, one from Lockheed-Huntsville and one from Ames 

Research Center, produced identical results. A fourth order damping term 

was appended to the RHS to help alleviate some of the oscillatory behavior. 

Ames indicated that in all their calculations with centered differences, 

fourth order damping was used. A fourth derivative term was therefore 

approximated and added to the RHS of the equations. The numerical diffusion 

cancellation terms were then dropped, except for the cases with shock. Use 

of the damping term eliminated some of the spatial oscillation but is highly 

dependent upon the value of an arbitrary coefficient which can vary between 

0 and 2. If too small a value is used the parameters oscillate, if too large 

a value is used the solution is overdamped and becomes linearized. A compromise value used throughout this study was 0.1 which worked quite well for 

most of the cases analyzed. 

Another idea that was investigated involved the use of a MacCormack 

operator to compute the RHS. It is well known that the two-step MacCormack 

operator gives second order accuracy and is very stable. This scheme worked 

quite well and eliminated the necessity of including the fourth order damping 

term. 

To this point all the calculations were done using three point centered 

difference approximations to the derivatives. This results in the basic block 

tridiagonal scheme. One can also formulate the equation set based upon a 

48 
backward difference approximation which then results in a block bidiagonal 

scheme. This approach worked very well and, as is well known, has excellent 

stability characteristics. Its major limitation is that it is only first order 

accurate and generally is applicable only to supersonic flows. This scheme 

is inherently stable for any step size, and several cases were run at Courant 

numbers of 1000. 

Instead of using a pure centered scheme which has stability problems 

or a backward difference approximation that is only first order accurate, a 

combination of weighted differences was evaluated. Sever al combinations of 

weight factors were investigated, such as 2/3 centered plus l/3 backward, 

and generally it was found that a slight increase in the Courant number could 

be obtained over that required for the pure centered scheme. Accuracy remained about the same as the centered differencing scheme. 

As Fig. 3-2 shows, the solution technique previously discussed produces 

very reasonable results including the location of the normal shock in the 

diverging duct. 

49 
3.0 

2.0 

C 

Dissipation 

plus Selective NDC 

Damping 

\ 

NDC Damping 

ytic. Sdlution 

1.0 2.0 

Area Ratio 

3.0 

Fig. 3-2 - Quasi-One Dimensional Results. 

50 
4. DEVELOPMENT OF GIM/STAR SPATIAL 

MARCHING ALGORITHMS 

4.1 INTRODUCTION 

The GIM/STAR SE-l code treats the full elliptic flow field us,ing explicit 

finite difference methods. This technique is applicable to a large range of 

fluid dynamics problems and has been successful in computing a number of 

these. The current code can be an “overkill” for some problems of interest 

in that a full elliptic treatment is not necessarily required. A parabolized, 

spatial marching algorithm could provide accurate flow fields for these situations and would be considerably more economical. 

The elliptic code is constrained by two items which restrict its use on 

large three-dimensional viscous flows: 

1. 

2. The large amount of data storage needed for threedimensional viscous flows causes large “page faulting” 

on the STAR machine. Any finite difference method, 

explicit or implicit, still requires the large data base. 

A GIM/STAR code with a parabolic spatial marching 

algorithm would not attack as many kinds of problems 

as the elliptic version but would allow large threedimensional viscous flows to be treated with no page 

faults on STAR. 

The time step in explicit schemes is restricted by the 

CFL and viscous stability limit. This is usually controlled 

by the small grid sizes normal to no-slip boundaries. 

If inviscid, free slip boundaries can be used, i.e., ignore 

the boundary layer, then the severity of this constraint 

decreases. The implicit, linearized block methods 

described in Section 3 provide a possible remedy for 

the time step difficulty in the elliptic code. 

The intent of this research is to provide both an elliptic, time-dependent 

GIM/STAR code and a hyperbolic/parabolic spatial marching version. It 

51 
is not too difficult to conceive of the future codes which could contain switching logic to automatically change from elliptic to parabolic etc., depending 

on the phys its of the flow, -however this will not be attempted here. The 

section will review the available classical parabolized methods and their 

problems, and then present an idea believed to be new for computing quasiparabolic flows. 

4.2 CLASSICAL PARABOLIC APPROACHES 

Most of the literature on spatially marching schemes, hyperbolic or 

parabolic, treat equations which have been transformed to a Cartesian computation grid which is uniform. The space marching can then be done in 

much the same manner as time marching. This is a good approach if a 

single transformation exists for the full flow domain. The GIM code strategy 

has been to compute in the physical domain whereby completely arbitrary geometries can be treated. This approach presents a problem in developing a space 

marching algorithm, i.e., the fact that the geometry changes in the marching 

coordinate direction. This is akin to a GIM unsteady time marching scheme 

whereby the geometry is allowed to change with time. If we are to keep the 

GIM strategy of arbitrary geometries, then a space marching scheme must 

be developed which will account for the geometric variations in the streamwise direction, i.e., non-uniform computational domain. 

The recent work of Roberts and Forester (Ref.5) use a boundary-fitted 

computational mesh in a parabolic c.ode for ducts of arbitrary cross section. 

Their algorithm for solving the equations appears to be a refinement of the 

classical method of Patankar and Spalding (Ref. 6). Rubin and Lin (Ref. 7) 

presented a nonlinear, iterative finite difference method for three-dimensional 

viscous flows. A parabolic method using a block implicit type scheme was 

given by Hirsh (Ref. 8). The solutions were restricted to supersonic flow 

(shear layers) of the free mixing type. Lubard and Helliwell (Ref. 9) calculated flows on cone at angle of attack using a parabolized method. This paper 

discussed some of the inherent difficulties with singularities, ambiguities and 

52 
departure solutions which arise in parabolized algorithm. The paper discusses explicit and implicit schemes for parabolic marching flows. 

Lin and Rubin (Ref. 10) presented a method using psuedo-time relaxation 

with a space-centered implicit differencing technique. They discuss many of 

the problems inherent in “pure” parabolic marching and show how time relaxation can eliminate departure solutions. The GIM technique, although developed 

independently of Lin and Rubin, also employs time relaxation but with an explicit, one-sided, predictor-corrector scheme and arbitrary three-dimensional 

geometries. The second order backward-forward, backward-backward explicit 

scheme of the GIM code is also a unique approach to parabolic marching solutions. 

For problems in which viscous terms can be neglected entirely and the 

main flow direction remains supersonic , we would like the capability in the GIM 

code to resort to a simple hyperbolic algorithm. The classical methods presented in the literature for parabolic and hyperbolic flows are drastically different because of the treatment of the pressure terms in the marching direction. 

As long as the flow is inviscid and supersonic, the axial pressure terms can be 

treated exactly. However, for subsonic flow, for example, several problems 

arise in applying a hyperbolic algorithm to the parabolic equations. This is, 

however, the approach that would be most general, if the “parabolic pressure” 

problem can be treated, 

Certain assumptions must be made in using a spatial marching technique. 

l There must exist a dominant flow direction in which to march. 

There can be no flow back upstream, i.e., no recirculation in the 

streamwise direction. 

l Stress terms are not allowed to act on the cross planes: i.e., there 

can be no second order terms (diffusion, viscosity) in the marching 

coordinate. 

l The downstream pressure field must not be allowed to propagate 

upstream. 

There are a number of strong implications in these assumptions. A supersonic, inviscid flow satisfies them all. A supersonic viscous flow will conform to the assumptions if the viscous terms are dropped in the marching 

coordinate direction. 

53 
Consider now the problem of spatial marching in a subsonic viscous 

flow. The first two of the assumptions can be met by. simply,not allowing any 

flow reversal problems to be attempted and dropping all streamwise diffusion 

terms. The downstream pressure field can still feed back through a subsonic 

stream. One obvious approach is to drop the streamwise pressure gradient 

term. This would satisfy the third assumption, but it appears a serious 

matter to simply drop this important term. 

Another approach commonly used is to provide a separate, explicit equation for the pressure and use windward, one-sided differences. The most exact 

way is to compute the conserved flux parameters and then l’decodel’ for the 

velocity, density and energy and compute the pressure from a state relation. 

The ideal gas law, a set of equilibrium thermodynamic relations or Boussinesq 

equations, is used to couple the state variables. Each of these approaches contains inherent difficulties which render their general use questionable. The 

following is a summary of some of these problem areas with classical parabolic 

s theme s. 

Zero Axial Pressure Gradient 

This does not cause any significant numerical problems in computing a 

flow field. It does however create a major problem in that the computed answers 

are probably wrong for most flow fields. A mixed supersonic/subsonic flow, for 

example, with a shock wave crossing the flow field cannot be computed at all because of the large axial (and radial) gradients. Some researchers still proceed 

to use this approach and try to justify it. 

Exact Pressure Treatment 

The rigorous way to compute the parabolized equations is to include the 

pressure in the conservation variable state vector for the momentum equations. 

A state equation can then be used to “decode” for the pressure. The advantages 

of the approach are that: (1) fully conservative differencing can be used; (2) shock 

capture algorithms are applicable; and (3) an auxiliary differential equation for 

the pressure is not needed. However, there are major problems with the “exact” 

treatment of the parabolic pressure. 

54 
l One-sided upstream differences must be used 

l The “decodel’ is ambiguous at Mach = 1 since two 

roots appear for the velocity (or pressure) 

. 

l Real viscous, no- slip walls cannot be tre&ted since 

/’ the,deco’de is singular: . 

1 l ’ Flows w?th a quiescent part, such as jets exhausting 

._ 

into an ambient, motionless atmosphere cannot be 

treated because of the singularity for zero velocity. 

Consider the two-dimensional parabolic system 

where x is the marching coordinate, y the cross plane (or radial) coordinate, 

E is the state vector of conservation variables and F is a nonlinear (viscous 

plus convective terms) function of E. 

A typical state vector E, for the parabolized Navier -Stokes equations is 

E = 

P= (Y - 1)p 

[ 

($ - u2 f 1 

Here, u is the axial velocity, v is the cross plane coordinate velocity, p is 

density, 6 is the total energy per unit volume and P is the pressure. Suppose 

a calculated value for the E vector exist at a plane X = Xi. It is now required 

to “decode” for the primitive variables. The following is one decode procedure that can be used in computer codes. 

55 
(1) v = E3/El 

(3) P = El/u 

E4 

(4) b? = yE + 

1 

9 (u2 t v2) 

(5) p = (y=l)P E [ - “2t v’] 

Two problems are immediately obvious: 

The radical in the u velocity decode causes an ambiguity. 

It can be easily shown that the correct decode is to take 

the t sign for u supersonic and the - sign if u is subsonic . In mixed flows, the sonic nature of a grid point 

is not known a priori. This Mach= 1 ambiguity prohibits 

a general parabolic marcher from being developed using 

the classical notions. See Section 5 of Appendix B. 

The axial velocity, u, cannot be zero or the decode is 

singular. The axial component must be zero, however, 

if a real wall is to be put into the problem. All classical 

parabolic codes simply use some wall functions or resort 

to inviscid slip conditions to avoid the singularity. 

A third difficulty, which is not so obvious, is that attempts to use implicit 

methods to march the solution downstream often fail. The reason is that 

the boundary conditions are not treated exactly, and these errors build up as 

the streamwise coordinate is traversed. Often, the explicit differencing of 

points near the wall is used as a patchwork way of circumventing the boundary 

condition difficulty. 

56 
Explicit Treatment of Pressure ’ 

This is the most widely used. of the parabolic procedures and its origination is usually attributed to D. B. Spalding. The idea is to -provide an explicit 

differential equation for the pressure field in addition to the basic conservation 

laws. This is usually a Poisson-type relation obtained from combining continuity and momentum equations. Satisfaction of local mass conservation is 

generally the criteria used for convergence of the elliptic Poisson equation. 

In general, a state vector will have the following appearance: 

where the E5 component now represents the differential equation solution for 

pressure from whatever means. 

Now note the difference in the “decode” from the exact treatment case: 

P = E5 

E2 - Es 

u = 

E1 

V = E3/E1 

2 

P = 

E1 

E2 - E5 

57 
The Mach = 1 ambiguity is no longer present as the radical does not appear. 

Thus mixed subsonic/supersonic flows .can be computed -without a’ priori knowledge, of the Mach number. Note, however, tha.t the. decode still contains the 

axial velocity in the. denominator.. Real solid walls cannot.enter if viscous 

boundary conditions are used. 

I 

This “explicit pressure 11 treatment requires solution, at each plane, of 

a Poisson-type equation. Thus, an iteration between planes is required before 

moving on down to the next plane. Even with its inherent bad points, this approach remains the most successful and widely used parabolic algorithm. 

4.3 THE QUASI-PARABOLIC IDEA 

The results of the initial investigation of a parabolic/hyperbolic GIM 

code led to the conclusion that there’.just is not a good approach being used 

today that fits the GIM code strategy. Three basic requirements were placed 

on a GIM/parabolic algorithm: 

l The geometric treatment must be applicable to arbitrary 

shapes. 

l The same basic algorithm should be applied to both hyperboiic and parabolic flows and be capable of eventual coupling 

with an automated algorithm for switching back and forth to 

the elliptic solver. 

l The algorithm should be readily vectorizable to realize the 

speed gain from using the STAR computer. 

In terms of a “classical” parabolized spatial marching algorithm, several 

geometric approaches were investigated. 

The first approach considered would generate the geometry plane by 

plane as the solution evolves, assembling the elements locally at each step. 

This would of course mean that the GEOMETRY module would be called at 

each integration step, thus coupling the geometry and the flow. An advantage 

58 
of this approach is that only the amount of geometry needed would be computed and stored at any cycle of the calculation. This would reduce the 

computer storage and reduce the large input/output problems. However, 

this approach would also require considerable reprogramming to make the 

GEOMETRY module of GIM a subprogram to the INTEGRATION module. 

A second approach appears to be a treatment of the geometry uncoupled 

from the flow. This means that all geometry, transformations and element 

assembly would be done before the flow field is integrated. The matrix data 

would be read from a stored file for each cross plane as needed. The advantages of this approach are the geometry is computed only once for a 

given configuration, the geometry module can be separate (as it is now), 

from the integration module and the grid could be inspected for desirable 

character prior to computing an expensive flow field. Disadvantages are 

that the basic character of the flow must be analyzed a priori to place grid 

planes in desirable locations, and data must be read from files at each cross 

plane which could effect the thru-put time on the computer. 

A third possiblity is to switch to computing in a transformed computational space. This makes the marching algorithm straightforward but 

forfeits one of the major advantages of GIM - completely arbitrary 

geometries. 

Approach 2 was selected as the best compromise and also provides the 

ultimate capability of elliptic-parabolic switching discussed earlier. 

The classic algorithms for treating the parabolic pressure field 

were deemed unsatisfactory. The following idea evolved from this 

research. The approach is termed “Quasi-Parabolic” and arose from 

the requirement of eliminating the ambiguities and singularities of 

existing methods. 

59 
The basic idea is to combin,e the classical parabolic marching approach 

with a II qua& time” relaxation. The parabolic-march procedure greatly reduces the amount of computer storage compared to a fully elliptic field. The 

time relaxation form of the equations eliminates the “decode” ambiguity associated with the parabolic pressure problem and allows velocity boundary conditions at solid walls to be treated. The equations used in the QP method are 

the time-averaged full Navier-Stokes, but with all second order terms dropped 

in a quasi-marching coordinate. Another way to view the QP equations is to 

take the parabolized Navier -Stokes and add back t’psuedo time” derivatives. 

The QP solution procedure, as any parabolic marcher, thus allows no downstream diffusion effects or pressure wave feedback through a subsonic flow. 

The solution is assumed known at upstream data planes, 1, 2, . . .N-1, and the 

solution is sought at plane N with no knowledge of plane N t 1. “Psuedo Time” 

- 

relaxation, is used to obtain the solution at plane N in terms of the (converged) 

solution at a number of upstream data planes. Backward differences of some 

type, (second order) must be used to prohibit downstream feedback. So the 

QP algorithm is not a classical space marching scheme, and is also not a 

time-dependent elliptic method. It is somewhat of a hybrid technique which 

combines the better features of two approaches and eliminates the bad ones. 

The GIM/STAR elliptic code will converge a case in 500 to 1000 steps 

if the initial guess is chosen reasonably close to the answer. Also, GIM/STAR 

is relatively cheap to run, if the problem size is small enough to fit into memory 

and not require large page faults. The QP algorithm relieves both of these 

difficulties to some extent. By storing only a small number of data planes 

(and not the entire elliptic field) the large page fault problem is gone. The 

QP marching procedure can also assign a reasonable guess to the Nth data 

plane since it knows the upstream converged solution, i.e., guess it is equal 

to the N-lSt plane or extrapolated in some way. This should allow the time 

relaxation to converge very rapidly. If an implicit time-relaxer is used 

with the QP algorithm (with steps many times the CFL), the relaxation should 

go even faster. 

60 
The QP method allows an exact treatment of the Parabolic pressure 

field. No ambiguity exist in the QP decode at Mach = 1 (since it is “quasitime” dependent) and no-slip walls can be treated exactly, i.e., boundary 

layers . The QP algorithm eliminates many of the bad features of pure 

parabolic methods. 

One obvious disadvantage of the QP approach is the planewise iteration 

(time relaxation) which must be done. This can be time consuming on the 

machine, and a good criterion for convergence must be used to avoid error 

propagation downstream. Spalding’s method suffers from this same planewise iteration to correct the pressure as well; other linearization schemes 

such as Roberts and Forester (Ref. 5) which use the conservative equations 

also suffer. Planewise iteration is not uncommon in most parabolic methods, 

thus the QP scheme is no better or worse in this respect. A linearized block 

implicit scheme, as discussed in Section 3, appears to be very attractive for 

performing the quasi-time relaxation. 

Figure 4- 1 shows the QP form of the three-dimensional Navier-Stokes 

equations in Cartesian coordinates. Note that these are the classical parabolized form plus a psuedo-time derivative. Included are global mass conservation, three components of momentum conservation, total energy and 

an equation for conservation of individual species in a binary mixture. 

Figure 4-2 is a typical computation molecule for a QP type marching. 

Assume that all flow variables are known at planes 1, 2.. .K and the solution 

is sought at plane K t 1. If backward differences in x are used, (first order, 

second order, etc.), then the scheme of Fig. 4-3 allows no downstream 

feedback, and allows plane K t 1 to be uniquely determined from upstream 

information, i.e., quasi-parabolically. 

Now consider the ultimate, not immediate, implications of such an 

algorithm. A flow field could be marched out quasi-parabolically from an 

initial data plane 1 to plane K, where K is set a priori ‘by the users. An imbedded elliptic region is encountered between planes K and K tM. The 

number of planes that the QP algorithm can treat on any given sweep is not 

61 
av+E+~+*-o at ax ay a2 

P 

PU 

u= pv i1 PW 

Pt 

PC 

E= 

F= 

bu 

pu2 + P 

PU v 

. 1 PU w 

(pE + Pb 

PUC 

- VT -wr -q 

YY YZ 

-Pw 

pWU-T 

z 

PWV- T 

Y= 

pwZtP-Tez 

(p&t P)W-UTZ - VTYZ-WT - q 

zz z 

PwC - RZ 

Y 

t- X 

Z 

(X is QP coordinate) 

Fig.4-1 - Three-Dimensional Quasi-Parabolic Navier -Stokes System in Cartesian 

Conservation Law Form. 
I 

I 

t 

Dominant 

Flow Direction 

---., 

X 

z 

I 

t 

(k- l)St 

Plane 

k th 

Plane 

Plane 

Fig. 4-2 - Symbolic GIM Code Computation Molecule. 
Equation 

Predictor 

A w-1 U. L,j,k 

= un * 

i,j,k - Ei,j,k-l ] t $e [E;,j,k-l - E’;,j,k-2] 

n 

j,k - Fi,j,k 1 [ At Gn 

- x i,j+l,k - Gn 

Lj,k 1 

Corrector 

ntl 

u. LUn 

I 

3 At 

L,j,k = 2 i,j,k 

+cntl wzz 

i,j,k 

-ET 

L,j,k-l] ’ t% [E’;,j,k-l - E?,k,k-2] 

- 6” 

i-l,j,k I [ 

At &n 

- z i,j,k 

Fig. 4-3 - The Quasi-Parabolic Scheme Finite Difference Equations. 
restricted to one. Simply specify single plane marching up to plane K, 

switch to the elliptic operator on the next M planes, and then parabolically march from the (K+M)th plane to the final th N plane. The switching 

can then be done II automatically, 1’ but the user must still determine the location to perform the switching. Eventually, perhaps, an algorithm could be 

written to detect the onset of a separation bubble, flow reversal, or other 

elliptic phenomena. This is not being considered at this time, but only the 

fact that the capability is within the framework of the QP algorithm. 

Advantages of the QP Algorithm Outlined 

There is no special treatment required of the parabolic 

pressure field. It is handled exactly except for the usual 

assumption of no downstream feedback. 

No ambiguity exists in the decode procedure at Mach= 1. 

Thus mixed flows can be treated with no a priori knowledge 

of the relative velocity magnitudes. 

Solid wall boundaries can be handled in the QP method 

with no-slip values. Regular parabolic procedures 

must avoid these type boundaries. 

Inclusion of more than one upstream plane will allow 

second order accuracy to be maintained inthe quasimarching coordinate. 

The QP scheme can accommodate either explicit or 

implicit It time” relaxation finite-differences. 

Within the basic framework of the QP scheme, an 

elliptic region could be treated before, during or 

after a marching integration simply by including 

k-data planes (instead of 2) during the relaxation. 

The QP algorithm requires very little addition storage over a classical parabolic method; and requires 

many times less storage than a fully elliptic treatment. Thus on STAR, the GIM/QP code could 

march out very large flow fields with no large page 

faults. 

By dropping the cross-plane viscous terms, the QP 

procedure becomes Quasi-hyperbolic with no further 

coding changes. Thus one algorithm can accommodate 

either parabolic or hyperbolic flows. 

65 
4.4 RESULTS OF COMPUTATION 

The QP c,ode has been essentially completed and a number of test. I 

cases exercised. Three test problems are shown in this section for illustration of the Quasi-Parabolic code. These cases are: 

1. Flow in a three-dimensional duct with an expansionrecompression and interaction of two shock waves. 

2. Flow over a 10 degree planar wedge. 

3. Two-dimensional viscous flow resulting from interaction 

of a nozzle exhaust with a supersonic freestream. 

Other cases are currently in progress, including a boundary layer calculation, 

containing subsonic and supersonic flow. 

The first problem shown is depicted in Fig. 4-4. The 1 x 1 square nozzle 

expands via a trignometric variation to 10 units and then has a constant 2 x 2 

cross section. The supersonic Mach number expands up to the 10 unit plane 

then, due to recompression, shock sheets form at,the top and outer side wall. 

The two shocks intersect as depicted in the figure. The QP code was used 

essentially in a hyperbolic mode with free slip inviscid solid walls and contained approximately 24,000 grid points. The purpose of the case is to determine the ability of the QP marcher to handle rapid expansions and strong 

compressions (shock capturing). 

The solution for Mach number at the lower wall corner is shown on 

Fig. 4-4. A comparison is attempted here with a forward-marching hyperbolic code of the classical variety (a MacCormack code). The GIM/QP solution shows a strong shock wave while the other marcher would not solve for 

the large gradients at all. Figure 4-5 shows additional profiles for this case. 

The pressure ratio (local to inlet) is shown for both the upper and lower wall 

corners. Comparison is made to a published solution (Ref. 11). Excellent 

agreement is seen for the smooth upper wall profiles and for the expansion 

portion of the lower wall corner. At the axial location where the shock intersection occurs, the two solutions differ considerably. The GIM/QP code, using 

a first order finite difference scheme agreed very well with the ATL results. 

66 
(2.y) = 1.5 - 0.5 .a+ t f$, x 5 10 

(2.T) = 2. x >lO 

Configuration 

GIM QP 

0 

M 

01 I I I I I I I I I I 

0 2 4 6 8 10 12 14 16 18 20 

X 

Mach Number vs Axial Distance 

on “Lower Wall” Corner 

Fig. 4-4 - Three-Dimensional Square Duct GIM Hyperbolic Computation. 
P 

PO 

Upper Wall Corner 

o 2 4 6 8 10 12 14 16 18 20 

X 

1.5 r- GIM QP 

Lower Wall Corner 

2 4 6 8 10 12 14 16 18 20 

Fig. 4-5 - Three-Dimensional Square Duct GIM Hyperbolic Computation 

(Axial Pressure Distribution). 

68 
However, the second order QP algorithm produced the curve shown in Fig. 4-5, 

i.e., a larger pressure rise. As a check on the accuracy of the QP shock wave, 

15,000 grid points were placed between 16 and 20 units. Very similar results 

were obtained as with the coarser mesh (11 x 11 x 81). It is thus felt that the 

GIM QP code is calculating the correct pressure rise across the shock. 

In order to test the shock-capture capabilities of the QP finite difference scheme, an oblique shock on a lo-degree two-dimensional wedge w&s 

computed. Two example cases were run with incident Mach numbers of 1.8 

and 2.4. The same 60 x 5 1 node grid was used for both calculations. Each 

case required about 26 seconds to converge. The results are shown in Fig. 

4-6 as the pressure ratio through the shock as a function of vertical position ” 

and pressure rise from the NACA 1135 shock tables. The shock was characteristically smeared over five grid points. The excellent agreement indicates a good shock-capturing capability with the QP second order backward 

difference scheme. 

Case three consists of a parabolic, viscous flow in the configuration of 

Fig. 4-7. A nozzle with high pressure exhausts into a lower pressure, hypersonic freestream flow. This case was solved with the GIM elliptic code with 

940 nodes and reported in Ref. 2. The QP algorithm gives virtually identical 

results as given by the full Navier-Stokes code. The grid used and the steady 

state Mach and pressure contours are shown on Fig. 4-7. Comparison of this 

solution with the reported values of Ref. 2 and with the inviscid SEAGULL code 

of Ref. 12 are shown in Figs. 4-8 and 4-9. The SEAGULL is an inviscid, slipline, shock fitting, forward marching code. Figure 4-8 shows vertical pressure distributions at three axial stations in the shear region, and Fig. 4-9 

gives the corresponding Mach number plots. As seen by the comparisons, 

the GIM marching algorithm does indeed work as expected and gives quantitatively the same answers as the other codes. Application to a boundary layer 

problem is currently under way. 

69 
5.0 

4.0 

a, 

c” 3.0 id u 

: 

;;i 0 .d &I 

fii 2.0 

:, 

1.0 

0.0 

m GIM 

- - - Shock Tables 

0 

Moo= 1.8 

8 = 44O 

1.0 2.0 

x= 3.519 

I 1 I 

1.0 2.0 

Fig. 4-6a - Wedge Shock Case to Verify Capture Technique in GIM/QP Code 

(Mm = 1.8, 9 = 44 deg). 
4.0 

3.0 

2.0 

1.0 

0.0 

- GIM 

--- Shock Tables 

M = 2.4 

em = 33O 

x= 3.519 

x= 1.944 

i 

I 

I -- - 

-- -1 

I 

‘I 

- I 

1.0 2.0 1.0 2.0 

P/P, 

I I 

1.0 2.0 

Fig.4-6b - Wedge Shock Case to Verify Capture Technique in GIM/QP Code 

(M, = 2.4, 8 = 33 deg). 
.L 

of-2 le Nozzle m 

30 k 1 ” 

‘1’ 

,A0 

r I I I 

2 3 4 5 6 

3 Freestream Axial Distance 

I Y I M I P I 

Nozzle 1.27 1.657 1924 

Freestream 1.27 5.0 106 

Configuration 

Mach Number 

Fig. 4-7 - Quasi-Parab o 1 ic Code Applied to Shear Flow. 
x = 3.4 x = 4.3 x = 5.6 

3.0 

2.5 

2.0 

* 

& 1.5 

: 

2 1.0 

.d i .5 

‘;;; 

u .d a 

b 

0 

> 

-0.5 

-1.0 

-1.5 

-2.0 

- 

- I I I I I I 

GIM NS 

SEAGULL 

GIM QP 

4 

I I I I I I I 

0 .1 .2 .3 

Pressure, P/p0 

Fig. 4-8 - Quasi-Parabolic Shear Flow Computation (Vertical Pressure Distributions 

at Three Axial Stations). 
3.0 

2.5 

2.0 

* 

1.5 

. 

: 1.0 

2 d 

: 0.5 

‘;;; 0.0 

u .d c) 

2 -0.5 

> 

-1.0 

-1.5 

-2.0 

r 

GIM NS 

x= 5.6 

2 3 4 5 3 4 5 3 ‘4 5 

Mach Number 

Fig. 4-9 - Quasi-Parabolic Shear Flow Computation (Vertical Mach Number Distributions 

at Three Axial Stations). 
5. REFERENCES 

1. Spradley, L. W ., and M. L. Pearson, “GIM Code User’s Manual for the 

STAR-100 Computer,” NASA CR-3157, 1979. 

2. Spradley, L. W., P. G. Anderson and M. L. Pearson, “Computation of 

Three-Dimensional Nozzle-Exhaust Flow Fields with the GIM Code,” 

NASA CR-3042, 1978. 

3. Beam, R. M., and R. F. Warming, “An Implicit Factored Scheme for the 

Compressible Navier-Stokes Equations, I’ AIAA Paper 77-645, Third 

Computational Fluid Dynamics Conference, Albuquerque, N. M., June 

1977. 

4. Salas, M.D., “Shock Fitting Method for Complicated Two-Dimensional 

Supersonic Flows,” AIAA J., Vol. 14, No. 5, May 1976. 

5. Roberts, D. W., and C. K. Forester, “Parabolic Procedure for Flows in 

Ducts with Arbitrary Cross Section,” AIAA J.,Vol. 17, No. 1, January 1979. 

6. Patankar, S. V., and D. B. Spalding, “A Calculation Procedure for Heat, 

Mass and Momentum Transfer in Three-Dimensional Parabolic Flow,” 

Int. J. Heat Mass Trans., Vol. 15, 1972, pp. 1787-1806. 

7. Rubin, S. G., and T. C. Lin, “A Numerical Method for Three-Dimensional 

Viscous Flow; Application to the Hypersonic Leading Edge,” J. Camp. phys., 

Vol. 9, 1972, pp. 339-364. 

8. Hirsh, R.S., “Calculation of Supersonic Three-Dimensional Free-Mixing 

Flows Using the Parabolic-Elliptic Navier-Stokes Equations,” Aerodynamic Analyses Requiring Advanced Computers, Part I, NASA SP-347, 

1975, pp. 543-565. 

9. Lubard, S. C., and W. S. Helliwell, “Calculation of the Flow on a Cone at 

High Angle of Attack,” AIAA J., Vol. 12, No. 7, 1974. 

10. Lin, T. C., and S. G. Rubin, “A Numerical Model for Supersonic Viscous 

Flow over a Slender Reentry Vehicle, II AIAA Paper No. 79-0205, New 

Orleans, January 1979. 

11. Dash, S. M., and P. D. Del Guidice, “Numerical Methods for Calculation 

of Three-Dimensional Nozzle Exhaust Flow Fields,” Aerodynamic Analyses 

Requiring Advanced Computers, NASA ~~-347, Part I, Mar. 1975, pp. 659-679 . 

12. Salas, M.D., “Shock Fitting Method for Complicated Two-Dimensional 

Supersonic Flows,~~ AIAA J., Vol. 14, No. 5, May 1976. 

75 

Appendix A 

USE OF THE GIM SE-2 CODE 

bY 

L. W. Spradley 

77 

Appendix A 

A.1 THE GIM/STAR SE-2 CODE 

The Blue Book (Ref. 1) describes the version of the GIM/STAR code 

designated SE-l (STAR Elliptic Version 1). This reference manual contains 

input guides and user information for the code. Since the publication of this 

Blue Book there have been a number of changes to the code which have not 

been documented. Some of these changes were necessary to allow large 

problems to be run with a minimal number of large page faults while others 

were made to reduce the possibility of wasting computer time generating 

MATRIX analogs on a bad grid. 

The input changes are not extensive but the user should use this 

Appendix in conjunction with the Blue Book when running a GIM/STAR 

problem. The following subsections describe the changes for the program 

modules and file usage. 

A.2 GEOM MODULE 

\ Module 1 of the GIM SE-1 deck was titled GEOMAT as it contained both 

the geometry and grid generation and the matrix coefficient assembly. The 

SE-2 version has the two operations broken out into separate modules. The 

first module of SE-2 is titled GEOM, as it now only performs the geometric 

description and grid generation (see Fig. A- 1). 

The user should be aware of the differences ,in this module between 

versions 1 and 2: 

l Input cards 16 and 17 (in the Blue Book) are no longer used - just 

omit cards 16 and 17. 

79 
Card Type 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

Parameter List/Format 

HEADER(I), I = 1,72 

(12A6) 

NZONES, IDIM, ISTEP, IMATRX, IMATE 

(515) 

IWRITE, LWRITE, NWRITE 

(315) 

KC(I), I= 1,6 

(6A5) 

NSECTS 

(15) 

MAPE(I), I = 1, 12 

(1215) 

MAPS(I), I = 1,6 

(615) 

(IBWL(I), I = 1,6), ITRAIN 

(715) 

(NNOD(I), I = 1, 3), (ISTRCH(I), I = 1, 3) 

(615) 

DIVPI(I), I = 1, 3 

(3E10.4) 

[AETA(J,I), I = 1, NNOD(J)], J = 1, IDIM 

(8E 10.4) 

[(Ac(I,K,J), I = 1,8), J = 1,4 or 121, K = 1,5 

(8E10.4) 

[AS(I, J), I = 1,8], J = 1,6 

(8E 10.4) 

(PT(I, J), I = 1, 5), J = 1,4 or 12 

(8E 10.4) 

Fig.A-1 - Input Guide for GEOM Module (SE- 2 Code). 

80 
l File 17 is not output from the geometry module, rather File 18 is 

now to be saved. This new File 18 is to be subsequently input to 

the new MATRIX module. 

l Card Type 4 has been changed. The values of a are no longer 

input, but rather a set of flags to retrieve the correct a’s are 

now used. The parameter KC .is set to alphabetic characters 

F, B, etc., for forward or backward differences, KC(l) is for 

x step 1, KC(2) for y step 1, etc., through KC(6) which is z 

step 2. The format is 6A5. 

A. 3 MATRIX MODULE 

The new MATRIX module of the SE-2 code now performs the analog 

coefficient calculation and file creation. This module should be executed 

following GEOM and before INTEG. The File 18 which was output from GEOM 

is now input to the MATRIX module. File 17 needed by INTEG is to be saved 

from MATRIX. Figure A-2 gives the storage requirements for MATRIX. 

The input to the MATRIX module consists of the Cards 16 and 17 which 

were omitted from the GEOM module, plus one new card. Each of these cards 

is now described: 

Card Parameters 

1 NDX, NDY, NDZ, ISNOPT 

2 KC(I), I= 1,6 

3 Nl, IC, NT 

Card Type 1 Format (415) 

Same as Card Type 16 (GEOM SE-l) p. 4-27 

NDX nodal decrement in the nl-coordinate system 

NDY nodal decrement in the n2-coordinate system 

NDZ nodal decrement in the n3-coordinate system 

ISNOPT special node treatment flag 

If ISNOPT = 1 the MATRIX module will calculate the 

number of special node terms placed 

on File 17 for input to INTEG 

Format 

(415 1 

(6A5) 

(315) 

81 
Matrix 2D - 

3D - 

/ACOM/ 

/PCOM/ 

(Q3MAP/ 

/IRFBC/ 

/JCFBC/ 

/PAFBC/ 

50*NN 196cNN 

4*NN 8 *NN 

24Hl-N + 18*%SPEC 48*NN t 18agNSPEC 

t 65542 t 65542 

6 *NSPEC 6::NSPEC 

6 :%NS PE C 6*NSPEC 

6*NSPEC 6*NSPEC 

where 

NN = total number of nodes 

NSPEC = number of special node terms allowed 

for in DIMENSION statements (DYNMAT 

input) 

The common block sizes may be calculated for each problem 

size to determine the ideal grouping on the LOAD card. If in 

doubt assign each block to a new large page boundary as below. 

Do not use GRLPALL, but use 

GRLP=*ACOM, GRLP=+PCOM, GRLP=W3MAP, GRLP=IRFBC, 

etc. 

Fig.A-2 - Module Common Block Sizes. 

82 
If ISNOPT = 0 the entire array of special node terms 

will be placed on File 17. The size of 

the array is determined by DYNMAT 

input (NSPEC parameter). 

Card Type 2 Format (6A5,) Analog Choice Card 

This card consists of a sequence of six characters (F or B) 

identifying the difference direction (forward or backward) 

for X, Y, and Z Step 1 and X, Y, Z Step 2, respectively. 

Examples : 

FFFBBB 

forward, forward, forward, backward, backward, backward 

for three-dimensional problems and 

F B - B F - 

forward, backward, backward, forward for two-dimensions. 

This is a new card for version SE-2 and is identical to GEOM card type 4. 

Card Type 3 Format (315) 

Nodal analog print control card. 

Nl first node of a print sequence 

IC print increment 

NT total number of nodes to print for this sequence. 

(See page 4-29 SE-1 manual for complete description.) 

Any number of cards of this type may be input. 

Place a -1 in Columns 4 and 5 on last card to terminate. 

Dynamic Dimension for MATRIX 

The new MATRIX module has its own dynamic dimensioning sequence. 

The same deck DYNMAT is used for MATRIX and GEOM, but the input of a 

third parameter is optional in MATRIX. 

83 
The DYNMAT deck is to be executed before running the MATRIX 

deck. 

The input consists of one card: 

Format 

(315) 

Card Parameters 

1 NN, IDIM, NSPEC 

The definition of these input variables are the following: 

Card 1 Format (315) 

NN number of nodes 

IDIM dimensionality (2 or 3) 

NSPEC number of special node terms to allow for in 

DIMENSION statements. If left blank, or zero, 

the arrays will be dimensioned to NN. The 

actual number of special node terms will be 

calculated and printed out in MATRIX. This 

value is then input to INTEG. (Not used by 

GEOM.) 

A.4 INTEG MODULE (SE-2) 

This module has remained virtually intact from a user standpoint. 

Three additional options have been added sine-e the Blue Book was issued. 

These are: 

l Capability to compute a CFL time step automatically over 

a multi-zoned grid. 

l Treatment of downstream subsonic boundary conditions using 

a mass balance condition. This option was added under another 

NASA contract and is documented here for completeness. 

l Input of a set of flags denoting the finite-difference direction. 

This aids in a more complete set of difference options and allows 

for full vectorization of all schemes. 

84 
Figure A-3 is a summary Input Guide for the INTEG SE-Z module. 

Note that the new input cards are designated Za, Zb, Zc, 3a and 6a. All except 

6a are optional and existing data decks will still work as they did for version 

SE-l. Card 6a must now be input in version SE-2. Card 2 has two additional 

inputs, IDS, IBOUND which control the optional input of 2a, 2b, and 2c. Zero 

values for the parameters on Card 2 signifies omission of the remaining Cards 

Type 2a, 2b and 2c. Figure A-3 is a description of the available options and 

each parameter that is to be input. 

Figure A-5 describes each parameter that is input on the optional Cards 

2a, 2b and 2c. An example of the use of the subsonic boundary condition option 

is shown in the sample grid of Fig. A-6. 

The time step calculation option is controlled via the value of KZONES 

read as the last data on Card 3. If this is omitted (=O), then one zone is 

as sumed. If KZONES > 0, then this signals the code that a multiple zone 

problem is being run. In this case, the value of KZONES should be equal to 

the number of zones used in the geometry module. If KZONES= 0, then Card 

Types 3a are not used, but any value of KZONES > 0, requires the input of 

come Cards 3a. The number of Cards 3a to be input is equal to KZONES-1. 

The time step information for zone number 1 is input on Card 3 itself. 

Figure A-7 describes the input of this time step information. 

Card 6a is simply the KC values used in GEOM and MATRIX, i.e., 

, FFF BBB 

in format 6A5. This card must agree with the previous module’s usage. 

One additional Fig. A-8 is included in this subsection. This chart 

shows formulas for determining the COMMON block sizes for the INTEG 

module. These values are needed for large problems to set up the LOAD 

card as described in the Blue Book (Ref. 1). 

85 
Card Type 

1 

2 

2a 

2b 

2c 

3 

3a 

4 

5 

6 

6a 

7 

8 

9 

10 

11 

12 

Parameter List/Format 

ICASE, IITITLE(I), I= 1, 78) 

(12, 78Al) 

IDIM, METHOD, ITMAX, IPRNT, ITSAVE, ISTART, 

IOTYPE, IUNITS, ITSTRT, IVISC, IDIST, ISPEC, 

IDS, IBOUND 

(1415) 

INFOUT, IJUMPO, JJUMPO, NIOUT, NJOUT, 

ICALC, AMFLW 

(615, E1O.O) 

INFINL, IJUMPI, JJUMPI, NIIN, NJIN, ICALC, 

OUTMFL 

(615, E1O.O) 

INFINL, IJUMPI, JJUMPI, NIIN, NJIN, INFOUT, 

IJUMPO, JJUMPO, NIOUT. NJOUT, ICALC 

(1115) 

NN, NNX, NDX, NNY, NDY, NNZ, NDZ, NPM, 

KZONES 

(915) 

KST, KNX. KDX, KNY, KDY, KNZ, KDZ 

(715) 

DTIME, DTFAC, INCDT 

(2E10.0, 15) 

REALMU, REALK, GAMSl, GAMSZ, WMl. WM2, 

DK, RK 

(8ElO.O) 

EMU, ELAM, ERHO, ESPEC 

(4ElO.O) 

KC(I), I= 1, 6 

(6A5) 

NNPM(I), NCPM(1). (NNCPM(1, J). J = 1, 5), 

ANGPM(1): I = 1, NPM 

(715, E1O.O) 

(NCT(I, J, K), PXPM(1, J, K), PYPM(1, J, K), 

K = 1.4); J = 1, NCPM(1); I = 1, NPM 

(15, 2ElO.O) 

RHOZ. PZ, ASTAR, NINC, A, B 

(3El0.0, 15. 2ElO.O) 

NJ, INC, NTOT, ITAN, ITYPE 

(515) 

RI, UI, VI, WI, PI, CSI 

(6ElO.O) 

Nl, IC, NT 

(315) 

Fig.A-3 - Input Guide for INTEG Module (SE-2 Code). 

86 
Card Cal. 

Type 2 l-5 

6- lo 

11-15 

16-20 

21-25 

26-30 

31-35 

36-40 

41-45 

46-50 

51-55 

56-60 

61-65 

Format 

15 

Variable 

ID IM 

METHOD 

ITMAX 

IPRNT 

ITSAVE 

ISTART 

IOTYPE 

IUNITS 

ITSTRT 

IVISC 

IDIST 

ISPEC 

IDS 

65-70 I5 IBOUND 

Description 

See Blue Book 

V’ 

Boundary Condition Flag 

= 0, one-sided differences 

(supersonic) 

= 1, mass balance technique 

(subsonic) 

Note: If IDS. Eq. 1 IBOUND 

should be set to either - 1, 0, or 1. 

If IDS. Eq. 0, IBOUND is left blank 

= -1, input inlet mass flow and 

calculate exit mass flow 

= 0, input exit mass flow and 

calculate inlet mass flow 

= 1, calculate both inlet and exit 

mass flow 

Note: If IDS. Eq. 1 card types 2a, 2b and 2c must follow type 2 card. The use 

of types 2a, 2b and 2c depends on the value of IBOUND. 

If IBOUND = - 1, use Type 2a 

If IBOUND = 0, use Type 2b 

If IBOUND = 1, use Type 2c 

Fig.A-4 - Definition of Parameters for Card 2. 

87 
Card 

Type 2a 

Type 2b 

Type 2c 

Cal. 

l-5 

6-10 

Format 

I5 

I5 

11-15 I5 

16-20 I5 

21-25 15 

26-30 I5 

31-40 E1O.O 

l-5 15 

6-10 15 

11-15 I5 

16-20 I5 

21-25 15 

26-30 I5 

31-40 E1O.O 

l-5 I5 

6-10 I5 

11-15 I5 

16-20 I5 

21-25 15 

26-30 I5 

31-35 I5 

36-40 I5 

41-45 15 

46-50 I5 

51-55 15 

Variable 

INFOUT 

IJUMPO 

J JUMP0 

NIOUT 

NJOUT 

ICALC 

AMFLW 

INFINL 

IJUMPI 

J JUMP1 

NIIN 

NJIN 

ICALC 

OUTMFL 

INFINL 

IJUMPI 

JJUMPI 

NIIN 

NJIM 

INFOUT 

IJUMPO 

JJUMP~ 

NIOUT 

NJOUT 

ICA LC 

Description 

Starting node on exit plane 

Nodal increment in ith direction on exit 

plane 

.th Nodal increment in J direction on exit 

plane 

Number of elements in i th direction on 

exit plane 

Number of elements in j th direction on 

exit plane 

Velocity update flag 

= 1, update inlet velocities 

= 2, update exit velocities 

Inlet mass flow rate (input by user) 

Starting node on inlet plane 

Nodal increment in the i th direction on 

inlet plane 

Nodal increment in the jth direction on 

inlet plane 

Number of elements in i th direction on 

inlet plane 

Number of elements in jth direction on 

inlet plane 

Velocity update flag (see Card Type 2a) 

Exit mass flow rate (input by user) 

Starting node on inlet plane 

Nodal increment in the i th direction on 

inlet plane 

Nodal increment in the jth direction on 

inlet plane 

Number of elements in i th direction on 

inlet plane 

Number of elements in jth direction on 

inlet plane 

Starting node on exit plane 

Nodal increment in the i th direction on 

exit plane 

.NNpilnncerement in the jth direction on 

Number of elements in i th direction on 

exit plane 

Number of elements in jth direction on 

exit plane 

Velocity update flag (see Card Type 2a) 

Fig.A-5 - Description of Input Parameters for Optional Card Types 2a, 2b, 2c (Subsonic 

Boundary Conditions). 

88 
6 18 30 42 54 66 

8 20 32 44 56 68 

10 22 34 46 58 70 

12 24 36 48 60 72 

14 26 28 50 62 74 

J=l J=2 J=3 J=4 J=5 

I= 1 

I=2 

I=3 

I=4 

Example: 

IBOUND = 0 INFINL = 6, IJUMPI = 2, JJUMPI = 12, NIIN = 4, N JIN = 5 

Fig. A-6 - Example of Subsonic Boundary Condition Usage. 

89 
Card 

3 

col. 

5 

10 

15 

20 

25 

30 

35 

40 

45 

Format 

15 

I5 

I5 

I5 

I5 

I5 

I5 

15 

I5 

Parameter 

NN - 

NNX 

NDX 

NNY 

NDY 

NNZ 

NDZ 

NPM _ 

KZONES 

Description 

See Blue Book 

The number of zones that was used to construct the grid. This is used to allow 

a CFL time step to be computed over multiple zones. Set to 1 for a single zone 

problem. 

Card 

3a 

Cal. 

5 

10 

15 

20 

25 

30 

35 

Format 

15 

I5 

I5 

I5 

I5 

I5 

I5 KDZ 

Parameter 

KST 

KNX 

KDX 

KNY 

KDY 

KNZ 

Description 

Starting node number 

of this zone. 

Number of nodes in q 

direction for this zon2 

Nodal decrement in r] 

direction for this zonk 

Number of nodes in V2 

direction for this zone 

Nodal decrement in ?j 

direction for this zong 

Number of nodes in q3 

direction for this zone. 

Set to 1 for 2-D flow. 

Nodal decrement in ‘7, 

direction for this zone. 

Set to 1 for 2-D flow. 

Note: Input Card Type 3a for each multiple zone to be used in computing a CFL 

time step. The number of Cards 3a is equal to KZONES- 1, where KZONES 

is input on Card 3. 

Fig.A-7 - Description of Parameters for Optional Card Type 3a Input. 

90 
Common Block 

Names 

/PRIM/ 

/TAU/ 

/TMvEC/ 

/VPROP/ 

/vBu~/ 

/BOUND/ 

/EBUF/ 

/xBUF/ 

/STEP/ 

/Ax~YM/ 

/Q3MAP/ 

Axisymmetric 

1 Gas 

5*NN+l 

9=*NN+3 

2sNN 

2 ‘:cN N 

8 =XNN 

5*NBtl 

8 *NN+4 

7 :kNN t4 

3*NNtlO 

8*NN 

2 Gases 

6:cNNtl 

9*NNt3 

2 ::NN 

4CNN 

10 *NN 

5*NBtl 

10 *NNt5 

7*NN+4 

3 *NNt 10 

9 :::NN 

24*NN+ 18cNSPEC 

t 6t COMP 

2-D 

1 Gas 2 Gases 

5*NNtl 6*NNtl 

9*NNt3 9*NNt3 

2 :%NN 2 :kNN 

2*NN 4 :<NN 

8 a%NN lO*NN 

5>gNBtl 5::NB+l 

8 4cNN+4 10 +NNt5 

7 *NNt4 7 ::<NN +4 

3ZcNNtlO 3::NNtlO 

8 9 

24*NNt18kNSPEC 

+ 6+ COMP 

3-D 

1 Gas 2 Gases 

6 :kNN 7 G;:NN 

12 x:NN 1 2 :cNN 

2 ::NN 2*NN 

2 ::NN 4 :RNN 

10 :kNN 12*NN 

5+NBtl 5 :cNB+l 

15 :::NN 18 *NN 

lO*NN+l 10:kNNtl 

4 *NNt9 4 ::NN t9 

9 10 

48:sNN+18+NSPEC 

+6+COMP 

NN = total number of nodes. 

NB = number of boundary nodes. 

NSPEC = number of special nodes. 

COMP = amount of storage need to complete a large page. 

Fig.A-8 - INTEG Module Common Block Sizes for SE-2 Code. 

91 
A.5 GIM SE-2 FILE DESCRIPTIONS 

Following is a brief description of the files used in the STAR SE-2 

system. In all but very small problems setups, (a few hundred nodes), a 

REQUEST card must be used for each file. The form of the REQUEST 

card is as follows: 

REQUEST (FILEXX/NSPGS, T = P) 

where 

NSPGS = 

Formulas for calculating NSPGS are now given for each file. 

the number of small pages of disk space 

allocated to the file 

1 small page = 5 12 words 

FILE 16 GEOM 

Work file used in GEOM only (Binary) 

NSPGS = 15+NN/512 NN = no. of nodes 

FILE 17 MATRIX/INTEG (Binary) 

Nodal analog file created in MATRIX and used in INTEG. 

2D - 

NLPGS = (16*NN t 18’:NSP t 6)/65536 

rounded up to next whole integer 

NSPGS = NLPGS*l28 t 1 

3D - 

NLPGS = (48::NN t 18*NSP t 6)/65536 

rounded up to next whole integer 

NSPGS = NLPGS’::128 t 1 

where 

92 

NN = total number of nodes 

NSP = number of special node terms. 
FILE 18 GEOM/MATRIX (Binary) 

File containing matrix assembly data. Created in GEOM and used in 

GA TRIX. 

2D - 

NSPGS = 5O%NN/512 + 1 

3D - 

NSPGS = 196*NN/512 + 1 

FILE 20 GEOM/INTEG/GIMTEK (Formatted) 

Nodal geometry file created in GEOM and used in INTEG and GIMTEK 

2D - 

NSPGS = 14+NN/5 12 

NSPGS = 2O+NN/5 12 

FILE 22 INTEG/GIMPLT (Formatted) 

Flowfield solution file created in INTEG and used both as a restart 

file and in GIMTEK. 

2D 

NSPGS = 1 WNN/5 12 per record 

3D .- 

NSPGS = 14*NN/5 12 per record 

Multiply by the number of iteration increments saved. 

93 
Controllee File Sizes 

The size of the controller file is specified on the LOAD ‘card in small 

pages. Formulas are given below for calculating the size required for a 

given problem. 

GEOM 

NLPGS = 

50*NN/65536 + 3 2D 

196::NN/65536 t 3 3D 

rounded up to next whole integer 

NSPGS = NLPGS*l28 This is the value that goes 

on the LOAD card. 

MATRIX 

No single formula exists to calculate the controllee file size for the 

MATRIX module. The procedure is to calculate the number of large pages 

(65536 words each) required for each GRLP parameter on the LOAD card, 

add these up, add 2 for other storage and multiply the result by 128. 

INTEG 

The same rule applies to the INTEG controllee file size as to the 

MATRIX module. Use the common block sizes to compute the number of 

large pages, add them up, add a couple and then multiply by 128 to get the 

controllee file size number. 

A.6 PLOT MODULE (GIMTEK) 

The GIM SE-2 code plotting module is now titled GIMTEK. This reflects the modifications which were made to the GIMPLT SE-l module in order 

to use the Tektronix 4014 for graphic output. The user need not be aware of 

the internal program changes that were made. The input data are identical to 

the SE-l version. Three items of significance to the user are now described: 

94 
The CM field length specified on the job card is calculated by 

CMlO = K-Xl0 + 2300010 

Notes: 1. 

2. 

3. 

This parameter must be set in the program and the array “A” 

dimensioned to this value. 

CM must be converted to octal for specification on job card, 

and RFL card 

Example 

NN = 2000 

KMAXIO = max 

1 

504~l~<2~oo~ = max /;:I;;/ 

= 22000 

CMlO = 22000 + 23000 = 4500010 

use CM = 130000 (1300008 = 45056 10'4500010) 

Fig. A-9 - GIMTEK Core Requirements. 

95 
l The formulas on page 6-25 of the Blue Book for computing core 

sizes for the plot module are no longer valid. Figure A-9 gives 

the revised formulas and an example calculation. 

l The plot save command was changed on the software system. The 

new save name is 

SAVPVF. 

l The routine that we use for obtaining GIMTEK plots from the 

Tektronix 4014 is 

PLIST. 

This allows enough options to select only those plots needed and 

also allows an unlimited time to examine a plot before proceeding. 

The input data for GIMTEK is ‘the same as described in the Blue Book. 

Figure A- 10 is a summary of the required input data presented here for 

comple tene s s . The user is referred to the Blue Book for a definition of the 

parameters. 

A.7 EXAMPLE RUNSTREAMS FOR THE SE-2 CODE 

The following pages show example runstreams that have been used for 

the SE-2 code on the STAR-100 machine. These should aid the new user in 

setting up a deck for GIM SE-2: 

Fig.A-11 - GEOM Module Only 

Fig.A-12 - MATRIX Module Only 

Fig.A-13 - GEOM/MATRIX Combination Run 

Fig. A-14 - INTEG Run Only 

Fig.A-15 - GIMTEK Run 

Note: The “blanks” which show up on the card listings are 7-8-q punch cards, 

i.e., end of record. 

The files for GIM SE-2 are cataloged under user number 838700C as GEOMB, 

MATRIXB, INTEGB and GIMTEKB. 

96 
Card Type 

1 

2 

3 

Specs. 

S-l 

s-2 

Grid 

G-l 

G-2 

VVEC 

V-l 

v-2 

I-l 

Contours 

C-l 

c-2 

I-l 

L-l 

Parameter List/Format 

ITITLE( l), ITITLE(2) 

(2A40) 

NX, ITERAD, ITRBLK, KDIM. ISP 

(515) 

GAMMA, FACTOR, RK, PO, TO, RHO0 

(6ElO.O) 

NPLT, STITLE, IVIEW, ISYM, ITHETl, IAXISl, 

ITHETZ, LAXISZ, IXTABL, VFAC 

(15, 5x, A20, 815, E1O.O) 

NTYPE, JO, IJUMP, JJUMP, NI, NJ, IPRNT 

(715) 

‘GRID’, IOPT, ICSCLE, NSPECS, (ISPEC(I), I = 1, 

NSPECS) 

(A4, IX, 15, 25X, 215, 715) 

(ISPEC(I), I = 8, NSPECS) if NSPECS > 7) 

(45X, 715) 

‘VVEC’, IOPT, NITER, ICSCLE, NSPECS, (ISPEC(I), 

I= 1, NSPECS) 

(A4, IX, 215, 20X, 215, 715) 

(ISPEC(I), I=8, NSPECS) (if NSPECS > 7) 

(45X, 715) 

(ITER(I), I= 1, NITER) 

(1615) 

ITYPE, IOPT, NITER, NC, ITABLE, INCR, ICSCLE, 

NSPECS, ISPEC(I), I= 1, NSPECS) 

(A4, lx, 515, 5X, 215, 715) 

(ISPEC(I), I= 1, NSPECS) (if NSPECS > 7) 

(45X, 715) 

(ITER(I), I= 1, NITER) 

(1615) 

(CVAL(I), I= 1, NC) 

(8ElO.O) 

Fig. A- 10 - GIMTEK Summary Input Guide. 

97 
GEuDCT~CM60000rT103. 

USEK*838700C. 

CHA~GEe101857*LKC. 

GETtULDPL=GEOM/UN=838i’OOC) 

GET ( DYNI~IAT=DY~I~~AT/UI~=~~&~OOC) 

UbDATE(F,C=TAPi%) 

DYWiAT. 

TOSTAI~~INPUTITAPE~) 

*ID NODS 

254 1 3 

STuKE 838700 400SDs TESTDECK t3 

STKSIDE~TlOO. 

FOKTKANt I=TAPE3r~=GEOl~i~r~=Ls) 

REUUEST(FILE16/75rT=P) 

RE~UEST(FILEl&/~74,T=P) 

REOUEST(FILE2U/lOO,T=P) 

LOA~(GEOMBICN=GEO,~G~,~~~~~GHLPALL= 1 

GEbVlGO. 

TOAS(Z=838700C~FILE18=~I~FIL~2O) 

*** GEOMETKY DATA *it*** 

Fig.A-11 - GEOP.4 Runstream. 

98 
MATRIX*CM6000urTlO0. 

USER 983t1700Ce 

CHARGE* 101857rLHC. 

GET t OLDPL=MAT~IX/UN=~~~~OOC) 

GET(DYNI~IAT=DYN~IAT/~~=~~~~OOC) 

UPDATE(F,C=TAPEB) 

COPYStiiF(TAPE8rOUTPUT) 

REWIND(TAPE8) 

DYNl”lAT. 

ATTACH(FILE~~=FILE~&B) 

TO~IAK(INPUTITAP~~,~ILE~~=~I//,U) 

*ID NONE 

254 1 3 1‘714 

STcJKE &38700 400SDs TtSTDECK B 

STHbIDE*TlOO. 

FOkiKAN(I=TAP~3rB=b’iATRB~O=LB) 

HEQUEST(FILE17/385,T=P) 

LOAU(~ATR~rCN=~~lATRGO,l~20 

. GRL!-‘=*AcOM, GRLP=-%Q~;~IA~ 9 GRLP=*t-,C;dw ,*k I RF BC .X JCFBC , XPAI- 0~ ) 

,vrATRGO. 

TOA~(L=~387OOC~FlLEl7=SI) 

121 11 1 1 

F F F d tl id 

1 1 20 

-1 

Fig.A-12 - MATRIX Runstream. 

99 
GEulilAT. CF16OOOb. T 100. 

USEKr83ti700C. 

CHAKGEI 101857qLRC. 

GET ( ULDPL=GE~;I~I/UN=~~~~OOC 1 

GET ( DYN~‘~IAT=DYI~IYA~-/UI~=~ jt17OOC) 

uPDATt(F,C=TAPt8) 

DYNl+lAT. 

COPYCF(TAPE3rGE( VIC) 

KEW 1 ND ( (;EOltJlC) 

KETURN(OLDPL) 

RETURN ( TAPE3 ) 

RETUKN ( T-APE8 ) 

GET( OLL)~L=;~AT~<IX/UIJ=~~~~OOC 1 

UPDATE(F,C=TAPtB) 

DYNMAT. 

COPY(3F(TAPE3rI”lATC;) 

KEW I iLCI ( I’wTC ) 

KETUkN(OLDPL) 

RETUKN(TAPE3) 

RETuKN(TAPE8) 

TOST AK ( INPUT 9 GEUI~IC 9 I/IA iC: ) 

254 1 3 

* ID I\ruNE 

254 1 3 1714 

STukE 838700 400SD5 TtSTDECI< B 

STKbIDk*TlOO. 

RE~JEST(FILE~~/~~,T=P) 

kEuutST(FILE17/38S,T=P) 

REU~EST(FILEl~/Y74,T=P) 

REUUESTCFILE2U/lOO,T=P) 

FOHTRAN( I=GEOlilC.Lj=GEOtqt)rC:=LB) 

LOAD ( GEOMUI CN=GEOMGu 9 1408. GRLPALL= ) 

GECmGU. 

- ***it GEOiQtTKY DATA -2 * -2 -2 * .* Q 

**** MATKIX DATA st*****+ 

Fig.A-13 - GEOM/MATRIX Runstream. 

100 
INTEGA*CM600Ud*T200. 

USEkr012839C. 

CHAtiGEv 1021101LKC. 

ATTACH(FILE17=FILE17A) 

ATTACH(FIL~~O=FILEZ~A) 

GET(OLDi=L=INTEG/UN=838i’OOC) 

GET(DYNDII~I=DYNDIIWUI~=~~&~~OC) 

UPDATE(F,C=TAPE8) 

DYND 114. 

COPYCF ( TAPERS I NTEGX ) 

REwINDi INTEGX) 

HETURNt OLDPL 1 

KETUKN(TAPE3) 

HETUKN(TAPEB) 

TO~TA~~(I~~PUT*II~T~GX~~IL~~~=UI//U) 

1225 1 0 175Y 

Fig.A-14 - INTEG Runs tream. 

101 
GIMTEK*CM120OW*T43~. 

USER*012839C. 

CHARGE* 1021101LRC. 

GET{ OLDPL=GIMTEK/UN=4Y2425C) 

UPDATE(F) 

FTN( I=COMPILE*L=O) 

ATTACH(TAPE2O=FILE23A) 

ATTACH(TAPE22=FILE22E) 

RFL(120000) 

ATTACH(LI~FTER~LKC(~OS~/C/I\~=LIBHARY) 

LDSET(LIB=LIBt.TEK/LRCGDSt-.~~ESET=NGINF) 

LGO. 

SAVE(SAVPVF=SCRJEl-) 

*ID KOQCHG 

2x. I G I I”IPLT .744 

ISET=I~ET+~ 

*o GIMPLT .Y 

COMILlON A( 14840 1 

*cl) GIMPLT. 15 

I<blAX= 14tj40 

**it** G I MTEK DATA ?-k 9 * x * * 4 

STOP 

Fig.A-15 - GIMTEK Runstream. 

102 
Appendix B 

THREE-DIMENSIONAL LB1 SCHEMES 

FOR THE NAVIER-STOKES 

EQUATIONS 

bY 

Jiirgen Thoenes 

Appendix B 

B. 1 INTRODUCTION 

Algorithms are developed for the solution of the three-dimensional 

compressible Navier-Stokes equations in conservation form. This work 

represents an extension of the methodologies outlined by Beam and Warming 

(Ref. B- 1) and Spradley (unpublished information) and familiarity with the 

cited literature is assumed. A time-dependent algorithm for the unsteady 

equations is developed first and then a spatial marching scheme for the 

three-dimensional parabolized steady equations is obtained. Algorithms 

for one- or two-dimensional problems are easily obtained by simply deleting appropriate terms from the equations. 

105 
B.2 THREE-DIMENSIONAL UNSTEADY ALGORITHM 

The three-dimensional compressible conservation equations can be 

written in conservative form 

=~[vlluJ’ux) + v 12 (U* Uy) + v13 (Us u,,] 

+3v21 (Us Ux) + vz2 (Us Uy) + vz3 (U, UJ 

t&P31 ( wJx) + 32 w* Uy) + V3,wJz)] 

where U is the vector of conserved variables and E, F, G and Vij are flux 

vectors. 

A generalized single-step temporal finite difference scheme for advancing the solution of Eq. (B.l), is the following. 

AUn = 

BAt 8Aun + At 8Un + L Aun-l 

iqat mat; 1+6 

(B-1) 

(B.2) 

where U n 

= U(nAt) and AUn = Untl - Un. (Terms of order At2 and At3 have 

been neglected, for simplicity. ) 

.106 
8U If Eq. (B.l) is solved for at and inserted in Eq. (B.Z), the resulting expression for AUn is 

(-AEn t AVTl t AVr2 t 4,) 

+ ,$& (-AEn f AV& + AVZ2 + AVg3) 

+ 2 (-AGn t AV;, t AVi2 t A$,) 1 

At 

+ 1+5 a(-, $ Vll t VI2 t V13)n 

ax 

+ ay ?- (-F + V21 t V22 + V23)n 

+% (-G + V31 + V32 + V33) 

n 

t & At?-' 

where AE n = En+’ _ En, etc. 

(B.3) 

n 

Note that AEn, AFn, AGn and AVij are nonlinear functions of the conserved variables U. A linear equation with the same temporal accuracy as 

Eq. (B-3) can be obtained by expanding AE”, AFn, AGn and AV: in a Taylor 

series, thus 

E ntl 

or 

AEn = An AUn (B.4a) 

107 
Similarly 

Bn AUn 

Cn AUn 

(B.4b) 

(B.4c) 

where A, B, and C are the linearization Jacobians. 

Strictly speaking, the same procedure should be applied to the viscous 

terms. However, as pointed out in Ref. B.l, treating the spatial cross-derivative 

terms, i.e., AVij (if j), in this manner would lead to considerable difficulties 

in constructing an efficient spatially factored algorithm. Therefore, spatial 

cross-derivative terms will be evaluated explicitly (without loss of accuracy, 

Ref. l), i.e., 

A?. 

‘J 

= AV;-l (if 3 (B.5) 

while the linearization is applied to the AVkk (k = 1, 2, 3). Remembering that 

Avkk = f(U,Ux ), 

k 

AVEk = ($)n AUn t (27 AUEk 

= Pkk AUn t RLk AUzk 

Application of the product differentiation rule shows that 

R;k AUn Xk = (Rkk Au);k - & Rkk m AUn 

(B.6) 

(B-7) 

108 
and therefore from Eq. (B.6) and (B.7) 

AVLk = Pkk Aunt (Rkk AU)n - Rkk x AUn 

Xk ’ k 

n 

= tpkk - Rkk,xk)n AUn + (Rkk AWxk (B.8) 

where Pkk and Rkk are the linearization Jacobians as defined in Eq. (B.6). 

Evaluation of these Jacobians will show that for constant transport coefficients 

-LR 

Pkk 8xk kk = ’ 

and thus 

AVtk = & k (Rkk Auf 

(B-9) 

(B.lO) 

If the approximations outlined above are introduced into Eq. (B.3), we obtain 

AUn = 

@At 8 

iqax 

1 [ 

-An Aunt & (Rll AU? 

I 

8 n n tay-B AU t 

[ 

5 (Rz2 AWn 

I 

t& 

[ 

-Cn AUnt&(R33 AU)" 

+ z &W,, + Av13) 

n-l 

+ 6 (AV21 + AV23) n-l 

t&(AV,, t AV32)n-1 

I 

(Continued) 

109 
At 

+ G-g (-E t Vll t VI2 t V13)n 

+ aY 

2 (-F t V21 + V22 + V23)n 

t &t-G + V31 + V32 + V33 )“3 

t & AUn-’ (B.ll) 

Thus, for constant transport coefficients only the Rkk linearization Jacobians 

are needed in addition to the A, B and C Jacobians. 

Expanding and rearranging Eq. (B.ll), we obtain 

Rtlll + 

a2 & Bn-- 

aY2 

a2 tgcn-- az2 RY3 * AUn 

=+$-$& (9, t AV,,)“-’ t & (Av21 + Av23) 

n-l 

t& (AV,, t AV,,)“-‘1 

At 

+1+5 C ax a (-E t Vll t V12 t V13)n 

+ w a (-F t V21 t V22 + V23)n 

6 n 

l-G + V31 + V32 + V33) 

3 

+ + Au”” (B.12) 

110 
Note the special notation used in writing the left hand side (LHS) of Eq. (B.12) 

which really must be considered an ‘toperator,lt operating on AU 

n 

, and which 

is of the form 

LHS(12) = Itatbtc (B.13) 

This can be written in a spatially factored form 

LHS(12) = I 

(I t a) (I t b) (I t c)’ 1 

* AUn 

= (Itatbtc) 

1 

t ab t ac t bc t abc * AUn I (B-14) 

if we note that ab, ac, bc, and abc all are at least an order of magnitude (in At) 

smaller than a, b, c. Thus, without loss of accuracy, 

LHS(12) = 

= 

+ 8At 8 An a2 n 

---(ax 

-- 

1+5 ax2 R1l) 1 * 

iLBT.-- 

( ay 

a2 

ax2 

R;2) * 1 

Ri3) * AUn (B.15) 

Following Beam and Warming (Ref. B.-l), in practice Eq. (B.15) is implemented 

by the sequence 

a2 -- 

ax2 R;l 

)I 

* AU** = RHS(12) 

a2 -- 

w2 

RT2 

)l 

* AU* = AU** 

(B.16a) 

(B. 16b.) 

where RHS(12) means the right hand side of Eq. (B.12). 

111 
* AUn = AU* 

U ntl 

= Un t AUn 

The remainder of the analysis follows that of Ref. B-l. 

(B.16c) 

(B.16d) 

112 
B.3 SUMMARY OF EQUATIONS FOR UNSTEADY ALGORITHM 

The vector of conserved variables, U, and the flux vectors of Eq. (B.l) 

are 

E = r Pu 

2 Pu +P 

F = 

G = 

PUV 

PUW 

GE + P) 1 = 

-Pv 

~ Puv 

2 

Pv +P 

Pm 

v(PE + P) 

-Pw 

Puw 

Pm 

Pw2+ P 

w(PE + P) 

= 

I = 

m 

m2/p+ P 

mn/P 

mq/p 

W/P) (r + P) 

n 

mn/P 

- I 

n2/p+ P 

W/P 

@/PI (r + P) 

- 

9 

m% 

W/P 

s2/p+ P 

(S/P) (r + P) 

113 
where the Pressure is given by the equation of state 

2 2 2 

P - u + “2 + w 

= (y-l) r - 

( 

n-l2 + II2 + q2 

2P ) 

The viscous flux vectors are 

vll = 

114 

0 

(2P + N ux 

pvx 

pwx 

(2P + A) uux t pwx t pwx + kTx 

0 

AV 

Y 

L I p”Y 

0 

AuvY + pvu Y 

t 

52 = 

v13 = 

v21 = 

- 

0 

AW Z 

0 

PUZ 

Auwz + pwu Z 

- 1 
p"Y 

v22 = 

(2/J+ A) vy 

pwY 

puuy t (2p t A) WY t pww yt kT Y 

'23 = 

V 32 = 

0 

PUZ 

V33 = PV, 

tzc-l+ N wz 

puu, t pw, t (2~ t A) wwz t kTZ 1 

where u = au/ax, etc. 

115 

x 
In order to write the viscous flux vectors in terms of the conserved 

variables, the temperature gradients must be expressed in terms of the conserved variables. It is easily shown that 

( 

am ap 

1 ( 

an 8~ --n - m P,5-mz -n Pa* ag 

1 

( 

as g -9 P,5-qac 

)I 

Using this equation it can be shown that 

- mP,) 

PP 

vll = 

-2 (PS, - 9P,) 

PP-3 I (z+$- $$ In (Pmx - mPx) 

I 

i (1 -j$ II nCPnx - nPx) + 9 (P9, - 

+ &Z PW X - W,) 

1 

W,) 3 

r0 1 

AP -2 W - nPy) 

V i-we2 WY - 

12= o 

WY) 

1 @3 m(Pn Y - nP,) + pP -3 n(Pm Y - mP,) i 

116 
v13 = 

v21 = 

v22 = 

0 

&f2 (PS, - ¶Pz) 

0 

E-LP -2 m 

-3 (p Z 

- mP,) 

w m(Pqz - SP,) + I-LP-3 9(PmZ - W,) 1 

0 

PP-2 (Pnx - nPx) 

Mm2 CPmx - mP,) 

0 

lJP 

-3 -3 

4Pnx - nPx) + AP n(Pmx - mp,) I 

I 1 

0 

IJP -2 (Pm - mp,) 

(2P+ 3 P 

-i (Pn, - W,) 

IJP -2 (PS, - spy) 

/Apm3 ((2 t h - Fr) n(pn 

1 

- nP,) 

+ (l-& [mlPmyY- mPy) + q(Pqy - 4Py)] 

- 

+ &i P (Pry - rp,) I 

v23 = AP -2 (P9 - 9P,) 

i-lP-2 (P” 

c3 

- np,) 

n(Pqz - ¶P,) + luP-3 WnZ - nPz) 

117 
v31 = 1 

0 

CLP -2 (P51x - qp,) 

0 

Ap:; (Pm, - mpx) 

IJP m(P9, - wx) + w-3s CPmx - mp,) L I 

‘32 = 

- 

0 

0 

pP-2 (PS, - WY) 

Mm2 (pn, - nP,) 

cLP-3 Wq y - 9p,) + c3 Wny - nP,) 1 

- 

0 

PP -2 (Pmz - mP,) 

I-LP -2 (Pn, - np,) 

V 33 = (2p+W P-2 (PS, - w,) 

I-lP -3 (2 + $ - g$, S(PS, - W,) 

1 

t (1 - g$’ m(Pm, II - mp,) + WnZ - W,) 1 

+ gr P(Pr Z - W,) 

I A 

118 
The linearization Jacobian6 A, B and C are: 

+ -; (PE t p)] 

0 

B+= -uv 

q(u2+v2tw2)-v2 

c = g= 

0 

-uw 

-VW 

J$J (u2tv2tw2)-w2 

wyE-+Et p,] I 

1 

- (Y- 3)u 

w 

@Et PI-(y-W2 

0 

- w- l)v 

U 

0 

- (y-1)uv 

0 

- o-l)w 

0 

U 

- (y-1)uw 

0 

t-Y-1 

0 

0 

1 

I F 

0 1 0 0 

V U 0 0 

- t-r-l)u - (y-3) v - (v-l)w (y-1) 

0 W V 0 

- (y-1)uv ;(PE* p) ty-l)v2 

- - 

(y-l)vw YV I 

0 0 1 0 - 

W 0 U 0 

0 W V 0 

- (y-l)u - (y-l)v - (Y’3)W (r-1) 

- (y-1)uw - (y-l)vw ;(PE+ p)-(y-l)w2 YW 
w 

The linearization Jacobians Rll, R22, R33 are: 

8 

r 

0 0 0 

(z+;)u - (2+$) 0 

avll 

R1l = aux = - 

f 

V 

W 

0 

0 

-1 

0 

(1 -&, (u2 + v2 + w2) - (2t;-&pl - (1 -&)v 

t gr E + (1 t;) u2 

r 

0 0 0 

I u -1 0 

av22 R22 = r = - 

(2$)v 

Y W 

(1 -&) (u2 + v2 t w2) 

t gr Et (lt;)v’ 

0 - (2+;) 

0 0 

-(l-&u - (2t; -+v 

U 

av33 V 

R33 = m = - 

2 (2$)w 

(1 -&r) (u2t v2 t w2) 

t gr Et (1t+)w2 

0 

0 

0 

-1 

- (1 -$p 

0 

0 

0 

-1 

- (1 -p+jw 

0 0 

0 0 

-1 0 

0 -(2$) 

-(l-&)v - (2+; -$gw 

0 - 

0 

0 

0 

-2 

Pr 

0 

0 

0 

0 

-2 

Pr 
- 

B.4 THREE-DIMENSIONAL STEADY PARABOLIZED ALGORITHM 

A flow model which uses a spatial forward marching procedure in the 

principal direction of flow to obtain a solution cannot tolerate the upstream 

propagation of any flow phenomena. Such a model is obtained by deleting 

those viscous terms from the governing equations which contain gradients 

in the marching direction, and the resulting set of equations is termed 

“parabolized.” 

The three-dimensional, steady state compressible conservation equations for parabolic flow (in the x-direction) are 

+: + E + +?$ = a ay ay II V22(E,Ey) + V23(E, Ez) 1 

+ V33(Es Ez)] (B.17) 

where E is the vector of conserved variables, and F, G and V.. are flux 

l&l 

vectors. 

In complete analogy to the treatment of the unsteady problem, a generalized single-step spatial finite difference scheme for advancing the solution of 

Eq. (B. 17) can be written as 

AEn = !?h&&,&.~:+~ a-1 itg ax it5 ax l-+5 (B.18) 

where E n 

= E(nAx,y, z) and AEn = E ntl 

- En (terms of order Ax2 and Ax 3 

have been negelected for simplicity). 

Solving Eq. (B.17) for aE/ax and substituting into Eq. (lB.18) yields 
AEn = +!$ J$ (-AFt AV22t Av23)n 

ts (-AGt AV32t AV33)n 

I 

& 

+ m ay [ 

a (-F t V22 t V23)n 

t$ (-GtV32t V33)n 1 t & AE~-~ (B.19) 

Again, AFn, AGn and AV: are nonlinear functions of the conserved variables 

E. A linear equation with the same spatial accuracy as Eq. (B..19) can be obtained by expanding AFn, AGn and AVLk in a Taylor series while treating 

AVFj (i f j) explicitly. Accordingly, 

An AEn (B.ZOa) 

Bn AEn (B.20b) 

where A and B are linearization Jacobians. 

While assuming that 

AVf: 

1J 

= AV;-’ 0 f j) 

we can write 

(B.21) 

= Pkk 

n 

AEn + R”kk AExk 

122 

(B.22) 
Since 

R!kq a AEn = & (REk AE") - & RFk AEn 

we can rewrite 

AVFk = (p;;tk -&R" 

k kk 

) AEn+ &(RFk AEn) (B.23) 

k 

where Pkk and Rkk are the linearization Jacobians defined in Eq. (B.22). 

Using the viscous flux vectors V22 and V33, it can be shown again that, 

for constant transport coefficients, 

p:k - axk aRzk = 0' 

which allows us to simplify Eq. (23) to 

(B.24) 

Av:k = axk -?- (REk AEn) (B.25) 

Substituting Eqs. 

AEn = 

t 

t 

t 

eh a 

iqay I [ 

-AnAEntF a (Rz2 AEn) 

I 

a 

+az [ 

-BnAEnt & (Rt3 AEn) 

1) 

+$f 

[ 

& (Av23)n-1 t 2 (AV32)n-1 I 

g [ ++ V22+ VZ31n 

t&(-G t V32 t V33)n I & AEn-’ 1+5 (B.26) 

123 
From Eq. (B.26) it is concluded that in addition to the A and B Jacobians we 

only need Sk ( k- - 2 , 3 ). if constant transport coefficients are assumed. 

Expanding and rearranging Eq. (B.26) to combine all terms containing 

AEn, we obtain 

a2 -- 

aY2 

A$;1 t & AV;,’ 

I 

(-F t V22 t V23)n t $ (-G t V32 t V33)n I 

In analogous fashion to the unsteady formulation and without loss of 

accuracy, we can write 

t AE~ (B.28) 

where LHS(27) means the left hand side of Eq. (B.27), which in practice is 

implemented by the sequence 

-- a2 

w2 

‘; A$ = RHS(27) 

* AEn = A$ 

(B.27) 

(B.29a) 

(B.29b) 

and, finally 

E ntl 

= En t AEn (B.29c) 

where RHS(27) means the right hand side of Eq. (B.27). 

124 
B.5 SUMMARY OF EQUATIONS FOR STEADY STATE ALGORITHM 

For steady flow, the vector of conserved variables, E, and the flux 

vectors of Eq. (B.17) are 

E = 

F = 

G = 

‘PV 

-Pu 

Pu2 + P 

Puv 

Puw 

u(PE + P) 

- I I 

1 r E3/u 

Puv 

E2 + E;/(EIu) - EMU 

E3 E4/(Eiu) 

E3 E5/(E p) 

E4/u 

E4 

E3 E4/(E p) 

E2 + E;/(E14 - Elu 

E4 E5/(E p) 

where u = f(E1, E2, E3, E4, Es) as obtained by decoding the E vector for the 

primitive variables. It should be noted that although the flux vectors for the 

steady case are the same as for the unsteady case, their form in terms of the 

125 
conserved variables differs from that in the unsteady case. It is assumed 

that decoding the E vector was accomplished by obtaining u as a function of 

the E vector components, i.e., 

and 

v = E/El 

w = E4’E 1 

P = E/U 

p = E2 - El u 

E = (Es - W/E1 

It can be shown that in the decode procedure the (t) and the (-) sign apply to 

supersonic and to subsonic flow, respectively. As long as the flow is strictly 

supersonic or subsonic, choice of the sign should be no problem. It is obviously a problem for transonic flows, boundary layers and supersonic flows 

with imbedded subsonic pockets. 

In terms of the conserved variables, the heat conduction term becomes 

Using this equation the viscous flux vectors V.. 

1J 

in Eq. (B.17) can be written 

as shown on the following page. 

126 
 

aE 

I-1 &(&$ i=l’ ” 

v22 = 

Ml+g(*)~ 

i=l 

+ lJ.t2+ z-&)E;~E~ El% h 

( 

+ Y(1 -&)Ei3E4 

+ ~ik E-2 aE5 

- _ 

Pr 1 1 ay 

aEl 

-E3ay 

> 

aE,\ I 

and 

v23 = 

127 
v32 

and 

0 

0 

k-i2 

A E;” 

I-1Ei3 

tA 

( aE4 E- a*l 

1 ay - E4ay 

> 

( 

aE3 E- aEl 

1 aY - E3F 

. > 

( aE4 BE1 

E3 Elay -E4ay 

> 

-3 

( 

aE3 aEl 

El E4 Elay -E3ay 

> 

1 

128 
The linearization Jacobians A and B are obtained by differentiating 

the flux vectors F and G with respect to the conserved variables. The 

result is 

A = A!. -Ayj (&) 

‘J j 

B = B!. -Bij(&) 

‘J j 

where (au/aE.) is a column vector obtained by differentation. 

J 

i 

1 1 1 1 

0 0 0 0 

pz 

2 2 2 2 2 2 2 2 A!‘. = u tv u tv u tv u tv 

LJ V V V V 

W W W W 

(Et;) (E + ;I @‘+;I (Et;) 

129 
.and 

B!. = + 

U 

L 

B!‘. 

LJ 

= pz 

Ic 

( 

0 

0 

-VW 

-(u2 + w2) 

-w(E t ;) 

0 ci 1 0’ 

0 0 U 0 

0 w V 0 

U 0 2w 0 

0 0 (E+pE) W 

1 1 1 1 1 

0 0 0 0 0 

V V V V V 

2 2 tw 2 2 2 2 2 2 2 2 1 u tw u +w u tw u tw 

W W W W W 

E + ;I W;) (E + ;I (Et;) (E + ;I 

! 

The components of the (au/aE) column vector are given by 

2!L= 

aE2 

au 

- = (+)y 2 aE3 - 1 0 1 

130 
au = (+)y 2 aE4 - 1 0 1 

where 

The linearization Jacobians Rkk (k= 2, 3) are obtained by differentiating 

the viscous flux vectors Vkk with respect to the spatial derivatives of the 

conserved variables Ei (i= 1,5). The result is: 

R22 = Ri2 t Q(pE) 

R33 = Ri3 t Q(pE) 

0 

0 

1 J!L 

R22 = pu 

-(2+$x7 

-w 

-+v2 

-(l -& (v2 t w2) 

- gr (E + ;I 

0 0 

0 0 

0 v+3 

0 0 

0 (2+ ; - &v 

0 

0 

0 

1 

(1 - + 

0 

1 

Pr 
0 

0 

-V 

(z++v 

1+$v2 

-(l -& (v2 + w2) 

-& (E+$ 

0. 0 

0 0 

0 1 

0 0 

0 (1 -+p 

1 

0 0 0 

1 1 1 

0 0 0 

0 0 0 

(l-j+ (1 -&)u (l-+ 

- 

0 

0 

0 

v+3 

(2+$ -$-)w 

0 0 

1 1 

0 0 

I 0 0 

Cl-& (1 +u 

i0 

0 

0 

0 

1 

Pr 

132 
B.6 REFERENCES 

B-l. Beam, Richard M., and B. F. Warming, “An Implicit Factored Scheme for 

the Compressible Navier-Stokes Equations ,‘I AIAA Paper No. 77-645 (1977). 

133 

Appendix C 

MULTILINEAR INTERPOLANTS FOR 

GIM MARCHING METHODS 

bY 

John F. Stalnaker 

Appendix C 

C. 1 INTRODUCTION 

The following is a study to determine the nature of the multilinear 

weight functions which will generate certain implicit, spatial marching 

finite difference schemes within the GIM framework. The derivations are 

performed using rectangular two- and three-dimensional grids for simplicity of understanding. In these grids the local and Cartesian coordinates 

coincide; however, it should be realized that this is not always, perhaps 

seldom, the case and that the finite difference scheme generated by the 

GIM code is in terms of the local coordinates. 

With the above caveat aside let us proceed by setting in one place 

the notation to be used herein: 

sa 

wP 

Dtk) 

“P 

f, g, h 

E, F, G 

E,@,f$ 

a, b 

AJ 

shape function for element point a 

weight function for element point p 

element difference operators for T)~ direction 

unique components of D (1) 

, Dt2) , and Dc3), 

respectively 

vectors of conserved variables 

Beam-Warming marching parameters 

finite difference parameters in the n2 and n3 

directions, respectively 

determinant of the Jacobian of transformation. 

137 
Subscripts 

Lj,k assembled grid point indices for 171, q2, and q3, 

respectively. 

= 1 ,..., 9 2-D 

= 1 ,..., 27 3-D 

asi3 element point indices = 1,. . . ,4 2-D; = 1,. . .8, 3-D 

In all the following q 1 (and x in the case of the rectilinear grid) is 

assumed to be the marching direction. 

C.2 TWO-DIMENSIONAL BILINEAR WEIGHT FUNCTIONS 

For the two-dimensional case the shape functions are assumed to be 

the same bilinear shape functions now in the GIM code (Ref. C-l): 

1 

(1 - T)l) (1 - 7-4 

sa = 

r)l (1 - 7-Q 

n1 n2 

(1 - rll) Q 

(C-1) 

The weight functions are assumed to have the form 

where the 

individual 

The two-dimensional element difference operators are 

(C-2) 

cpi are to be determined. The nodal numbering system for the 

elements and the nine node box are shown in Fig. C-l. 

(C-3) 

138 
Y 

t- X 

(a) 

i-l, jtl i, jtl i+l, j+l 

T v l 

i-l, j D ..i, j 0 i+l, j 

t 

AY 

i-l, j-1, A .i+l, j-l 

I- nr-4’j-l 

(b) 

Fig. C-l - Two-Dimensional Nodal Numbering System 

(a) Element System; (b) Nine-Node Rectangle 

139 
(C.3) 

(Conclfd) 

Noting that for the rectilinear element 

ax -= 

aQ 1 

Ax, -@ = Ay, = = 

an2 

-2.y = 0, 

aq2 h, A, = Ax Ay 

Substituting Eqs. (C.l) and (C.2) into (C.3) results in the following element 

difference operators 

k& 3 D(l) = - D(l) = 12 ‘a, [6 C 

P-2 Pl PO - 3cp1 - 2 cp2 + cp3 1 

& 3 D;:’ = - D;‘d = & [6 CPo - 3 CPl - 4 CP2 t 2 cp3] 

(C-4) 

?& E D12) = 

P4 - D;“1’ = & L.6 cpo - 2 cp1 - 3 732 + Cp31 

The differential equation 

3Et” = 0 

ay 

is modeled with a general spatial marching scheme (after Ref. C-2) 

140 

& [( 1 +E) E. Li-1, j - (‘+26) Eij + E Ei-1 j] 

, 

(C.5) 

+ ay ’ [a Fi+l,j+l + (1-2a) “i+l,j + (a-1) Fi+l,j-lI 

+ w [a Fi jt 1 t (1-2a) Fi j 

, , 

+ (a - 1) Fi j 1] 

, - 
-& [a Fi 1 j+l + (1 -2a) Fi-1 j + (a-1) Fi-1 j-11 - , , , 

(C.5) 

(Conclld) 

=o 

Examination of Eq. (C.5) reveals the significance of the difference parameters: 

a cross-plane (I) ) difference parameter 2 

= 1 for forward cross-plane differences 

= 0 for backward cross-plane differences 

= f for centered cross-plane differences 

E marching (I) ) difference parameter 1 

= 0 for forward differences 

= -1 for backward differences 

= -$ for centered differences 

8 weighting parameter for plane (itl) 

> 0 .th cross-plane differences are taken in the L 

and (i+l) st plane 

= 0 cross plane differences are not taken in the 

(itl) st plane 

4 weighting parameter for plane (i-l) 

.th 

> 0 cross-plane differences are taken in the L 

and (i-l) st plane 

= 0 cross-plane differences are not taken in the 

(i-l) st plane 

Note that for the present explicit differences in the elliptic GIM code 0 = 4 = 0. 

Assuming a form of the weight functions similar to that presently in 

the code, i.e., 

f!Pwp = AJ O (dpl - YIl) 52 - tl,) (Cha) 

141 
results in 

dPl = + [v gp1 - gp2)/(gp1 - gp)l (C.6b) 

and 

fPl + fpz = gp1+gp2 

Assembling the elements as in Ref. C-l in terms f’s and g’s, equating 

these to the difference coefficients in Eq. (C.5) and substituting into Eqs. (c.6) 

:: 

yields weight functions of the following form: 

“‘=irJ 1 “l (P - tl,) (2/3 - rl2) 

w2 = 

a2 

Tp2 - 171) (2/3 - rl2) 

w3 = 

“3 

nJ (P2 - ‘71) u/3 - rl$ 

w4 = 

“4 

aJ @ 1 - tll) (l/3 - T-t,) 

where 

“1 = 36a(l +E- 20) 

“2 = 36a(E- 24) 

(C.7a) 

(C.7b) 

Note that al s and p’s are substituted here for the d 

Pi 

in Eqs. (c.6) to show the 

similarity of these to the current GIM weight functions. 

142 
“3 = 36 (a-l) (e i 24) 

“4 = 36 (a-l) (lte-28) 

(C.7b) 

(Conclld) 

and 

1 2(1+e) - 38 

Pl = 3 1te-28 [ I 

(C.7c) 

1 e-34 

p2 = 3 e-2t#J [ I 

These apply except in the case where 

E = 2e- 1; ‘Go 

(e.g., Crank-Nicholson or e = 0, 8 = l/2). In these cases the weight functions 

do not maintain the same symmetry as the present weight functions. Values 

of the a’s and p’s are available from the author. 

C.3 THREE-DIMENSIONAL TRILINEAR WEIGHT FUNCTIONS 

The procedure for the three-dimensional case is much the same as the 

two -dimens ional problem. We assume the same trilinear shape functions now 

in the GIM code: 

sa = 

(l-r,,) (1-r,2) (1’r)3) 

171 (1-7?2) (1-r13) 

‘11 772 (l-r/3) 

b-ll) 772 (h3) 

h--q (1-r12) T)3 

t71 (l’r12) 773 

‘11 q2 t73 

b-tl) r)2 t73 

(C-8) 

143 
The weight functions are assumed to have the form: 

A, - W 

P 

=‘~po’~plr)l~cp2~2’c~~3~~+c~4~~f12 

(C.9a) 

or 

dO wp = ++ (d Pl - rll) (dp2 - rl2) (dp3 - 03) (C.9b) 

The relations between the Cl s and the dl s is obvious. The numbering system 

for the rectilinear element and the 27 node box are given in Fig. C-2. 

The full expressions for the element difference operators are given in 

Ref. C-3. For the rectilinear box, where 

A, = Ax Ay AZ, * = kc, e = Ay, $ = AZ 

ax. 1 

arl 1 

0 i f j, 

2 3 ‘arlj = 

These operators become 

D(l) 

pa = AY Az]1d,,/Ldq2/ld”3 wp 2 

Dt2) = P” Ax Az /Id?$‘dq20$dq3 WP 2 

Dt3) 

Pa 

= Ax Ay of&l ;‘drj2/:q3 Wp 2 

0 0 0 

(C. 10) 

Substituting Eqs. (C.8) and (C.qa) into Eq. (C-10) results in four unique difference operators for each direction 

144 
7)2 

A, 

Plane (i- 1) , 

‘13 

ql L Plane (itl) 

- \ 

Fig. C-2 - Three-Dimens ional Nodal Numbering Sys tern 

(a) Element System; (b) 27-Node Rectilinear BOX 

145 
It can be seen from the above that 

fPl + fp2 + fp3 + fp4 = gpl + gp2 + gp3 + gp4 = hpl + hp2 + hp3 t hp4 (C. 12a) 

and 

fPl + fp2 = gpl + gp: $3 + $4 = gp + gp4 

fpl + fp3 = hpl + hp2; fp2 + fp4 = hp3 + hp4 

gPl + gp3 = hp1 + hp4; ppz + gp4 = hp2 + hp3 

Now, Eqs. (C.qb), (C.ll) and (C. 12) can be combined to yield 

(C. 12b) 

fal 

fP2 

=!@A. %L=%L; &,G 

fp4 ’ 732 gP4 hi= hP3 

(C.13a) 

dPo 

= 216 

tfp1+fp2)-tfp3ffp4) Pg1+hp2)- th 3th 4) 

fpl+fp2+fp3+fp4 I[ hp1+hp2+hp3fh@4 13 I[ (gp1+ gp3) - (gp2+ gp4

1 

dpl = 3 L 2(g 

[ p1733) - (gp 2 + gp4) I[ / kp 1-b gp3) - (gp2 + gp4) 1 

dp2 = 3 L 2(h 

I QPh@2) - (hp3+hp4) / I[ (hpl + hp2’ - (hp3+ hp4) 1 (C.13b) 

dp3 = 3 1 2(f 

I p1+ fpz) - (fp3+ fp4) IL / (fpl + fp2’ - vp3+ fp4) 1 

The differential equation 

is modeled by 

146 
(1 i-e) Et;‘: - (lt2e)Ef kteE;-;] 

I , 

. 

+ gy La =;;;, k 

it1 t (I-2a)Fi kt (a-l)Fttl 

, 

] 

J-1, k 

t w [a F! 

. 

J+l, k t (1 - 2a) Fi , 

k t (a - 1) Fi- 1 k] 

, 

A -- 

. 

- Ay [a Fit:, k t (1 - 2a) Fir; t (a-l) FiB1 ] j-l, k (C. 14) 

it1 

’ & Lb =j, k+l t (1 - 2b) Gitl t (b - 1) Gitl 

j,k j, k-l ] 

t w [bGf 

J, k+l t (l-2b)Gf kt (b-l)Gt 

J, J, k-l 

] 

-AZ [ 

de bGi-’ 

j, ktl t (1 - 2b) Gi-’ 

j,k 

t (b - 1) Gi-’ j,k-11 = ’ 

where the difference parameters have the same significance as before and b 

is equivalent to a for differences in n 3’ 

Assembling the elements and equating coefficients does not lead to 

expressions for the weight functions in as straightforward a manner as in 

the two -dimensional problem. In order to resolve several ambiguities, the 

following considerations along with Eqs. (C.13a) were used: 

1. The weight functions should reduce to the form presently in the 

code for 8 = I$ = 0. 

2. The first four weight functions should readily reduce to the twodimensional case. That is, the internal symmetries of the twodimensional element difference operators should carry over to 

three dimens ions. 

3. The weight functions derived here should be applicable to the 

elliptic solver with 8 = 4 = 0. Thus, the boundary terms must 

be differenced consistantly. 

147 
From this, the following weight functions result: 

w1 

al 

= aT;(Pl -Q (5- v,) (3 2 - 773) 

“2 2 2 

w2 = Ajp2Vl) (T--r,,) (7-113) 

“3 1 2 

w3 = *;(P2-q1) t?j-T,& (T--rl3) 

“4 1 2 

w4 = *y ml- u1) t?j - 112) (3 - rl3) 

a5 2 1 

w5 = *py171) (j372) (3’r,3) 

a6 

w6 = hi+-74 (5- 7-/2) ($ - r,,) 

a7 1 1 

w7 = *y(P2-9) (~972) (3’1,3) 

“8 1 1 

w8 = ,--;@I-ql) (T-r,,) (T-773) 

where 

a1 = 

(a-b+ b2) (1 t E - 20); a5 = b(b- 1) (1 tc-28) 

a2 = (a-btb') (c-24) ; a6 = b(b - 1) (E - 2#) 

a3 = (a-2btb’) (e-24) ; a7 = (1 -b)2 (E- 24) 

“4 = 

(a - 2b t b2) (1 t E - 2@); a8 = (1-b)2 (H-E-2& 

and Pl and p2 are defined in Eq. (C.7c). 

148 
C.4 REFERENCES 

C-l. Spradley, L. W., P. G. Anderson, and M. L. Pearson, llComputation of 

Three-Dimensional Nozzle-Exhaust Flow Fields with the GIM Code,” 

NASA CR-3042 (1978). 

C-2. Beam, R. M., and R. F. Warming, ItAn Implicit Factored Scheme for 

the Compressible Navier-Stokes Equations; II, the Numerical ODE 

Connection,l’ AIAA Paper No. 79-1446 (1979). 

C-3. Spradley, L. W., and M. L. Pearson, “GIM Code User’s Manual for the 

STAR-100 Computer,” NASA CR-3157 (1979). 

149 

Appendix D 

VECTORIZED BLOCK TRIDIAGONAL EQUATION 

SOLVER FOR THE GIM/STAR CODE 

by 

S. J. Robertson 
, 
Appendix D 

An attempt was made to develop a vectorized algorithm for solving 

large systems of tinear equations of the form: 

L i-pi-1 t MiUit NitI Ui+l= Di (D-1) 

where the L, M and N elements are 3x3, 4x4 or 5x5 matrix blocks, and 

the U and D elements are three-, four- or five-component column vectors. 

The subscript i in Eq, (D.l) corresponds to nodal points in a computational 

grid, and the three, four or five dimcnsionality of the matrix and vector 

elements depend on whether the system of equations are for a one-, two- or 

three-dimensional flow field problem (see Section 3). The system of linear 

equations represented by Eq. (D.l) forms a block tridiagonal system. 

A solution algorithm was sought that would make use of the parallel 

processing capability of the STAR-100 vector computer. The Gauss-Seidel 

relaxation technique, based on an iterated solution of 

Ur+’ = (1 -w)U; - wMfl (Li-lU;-l t NitlUc+l - Di) (D-2) 

where w is an over-relaxation factor, is the only technique which we could 

find that permits a straightforward use of vectorized computer programming. 

-1 The inverse matrix Mi in Eq. (D.2) is evaluated for all Mi prior to entering 

the iteration loop. Since each matrix block is dimensioned only up to 5 x 5, 

the inverse can be evaluated by direct algebraic manipulation or by a Gauss 

153 
elimination technique. For the time being, we have coded only the algebraic 

inversion, since vector programming can be used in this method. 

Separate subroutines were programmed for 3 x 3, 4 x4 and 5 x 5 block 

tridiagonal Gauss-Seidel equation solvers. These are listed in Tables D-l, 

D-2 and D-3 as subroutines EQSOL3, EQSOL4 and EQSOL5, respectively. 

The argument list in these subroutines is (U, L, M, N, D, NODES, W, EPS, 

MAXI). The vector U is the solution vector which enters the subroutine as 

an initial or trial solution and returns as the updated or final solution. The 

matrices L, M and N and the vector D enter the subroutine as constants. 

The scalar NODES is the number of nodal points, W is the over-relaxation 

parameter, EPS is the error tolerance in the convergence test and MAXI is 

the maximum allowable iterations. The U and D vectors are doubly dimensioned, and the matrices L, M and N are triply dimensioned. The first subscript of both vectors and matrices corresponds to the nodal point index. The 

second subscript of U and D corresponds to the vector components, and the 

second and third subscripts of L, M and N corresponds to the matrix elements. 

As of this writing, these subroutines have not been evaluated, except for 

some very simple test cases. They have not been applied to realistic fluid 

dynamics problems where their usefulness can be determined. 

154 
Table D - 1 

LIST OF SUBROUTINE EQSOL3 

50 CONTINUk 

NC4 1 =NODES- 1 

NM2=NODES-2 

I TEK=O 

10 Cc)NTINUt 

1T~k!=1TE&+l 

DO 100 I=193 

UP(lrI)=(l .-;41*ut 19 I) 

UP(l*I)=UP(19I)-v~*l~II(l9I~l)*(N(l*l*l 

3. +~(lrlr3)*~(2v3)) 

1 

(Continued) 

155 
Table D- 1 (Concluded) 

UP(~*I)=JP(~‘I)-~~*I”II(~~I~~) ;;-(I\( 1939 1 )*b(2* l )+I\( 1*3*2)*U(2*2) 

s +N(lr3*3)*U(2,3)) 

up{29 1ow12)=( 1.--d)*~(2~ Isw2) 

UP ( 2 9 I LbhiY2 ) =UP ( 2 9 I blUidi2 ) - ‘iv++14 I ( 2. I 9 l robi;‘“12 ) * (L ( 2 9 1 9 l aN1b.12 ) 3&u ( l . l aai\lls12 ) 

3 +L ( 2 9 1 9 2biup42 ) *cl ( 1 q 251~12 ) +L ( 2 9 1 9 3wh2) *u ( l 9 3ww2) ) 

UP{ 29 IbNiY2) =clP ( 2 3 I~I‘JI$~)-Lv*vI’*II (29 I l lbiLM2) 7% ( I\; ( 2 * 1 * 15biWl2 ) *u ( 3 9 1 PIJIb ) 

51 +N(2r 1 r21~M2)*J(3*2~1~1~12)+~(2. 1 r3Ll‘\ll~i2)*C;(3r3~1\11~12) ) 

ut= ( 2 9 I W41’42 ) =d? ( 2 9 1 btW12 ) - W-~IY I ( 2 * I 9 25~NV12 ) T* (L ( 2 9 2 1 l brub’l2 ) %u ( l , l SNtll2 ) 

3 +~(2*.292bN~2)*J( l r281\1r~i2)+L(2r2~3~iu11~2)*u( l r3biw~2) ) 

uP~2~13[LI~12~=~rl(~~1bN,~12)-~*l~:1 (29 I ~2e~1~12)*(~(2r2ribN~~~2)-~~(3~ isw12) 

3~ +N ( 2 9 2 * 2bNivl2 ) ++J ( 3 9 2bi’U42 ) +N ( 2 l 2 9 3LIuiw2 ) *U ( 3 * 3Liwl2 ) ) 

UP(~~I~~~I~I~)=U~-‘(~~IL~I\~I~~~)-~J*I~II(~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

s +L ( 2 9 3 9 2bNp12 ) *-cl ( l 9 28hb12 ) +L ( 2 * 3 t 3brw2 ).-!+u ( l 9 3ardvi2 1 ) 

uP~2rIOlul42)~U~J(2,I~~l~2)- Lb*14 I ( 2 9 I 9 3LiYi’*i2) * ( C-4 (2 * 3 1 1 JJI\~+I~ ) *U ( 3 9 1 BNh2 ) 

0 +N(2*3*2~NF12)*~(3~231~~.12)+N(2*3.3~Nl~12)*U(3r33N;M2)) 

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

d *D ( 2 9 2sNM2 ) +14 1 i 2 9 I 9 3ibNr42 1 +D ( 2 I 3bNM2 ) ) 

UP(N~DESrI)=(l.-;~)*U(r~~DE~~I) 

UP(NODE~~I)~=~~~(I~~D~~~I)-~*~~~II~~~DE~~I~~)~~L~NUDE~*~~~~~U~N~I~*~) 

5 +L~NOD~S~~,~)*U~NM~~~)+L~NOD~~~~*~~~U~NI~I~*~) ) 

UP(NO~E~~I)~~P(~~~DE~~I)-~~~~~I~I~~D~~*I*~)~~L~~~D~~*~~~)~U~~~~~~~*~) 

B +L(NOD~~~~,~)*U(N~~I~~~)+L(N~DE~*~~~)*U(~I~I~~~)) 

UP ( NODES 9 I ) =ilr’ ( tudDkS * I 1 -bti*lh 1 ( htiDk 39 - 1*3)7~(L(lNu~cb93r l )ii-U(1\1,~11* 1) 

5 +L(NOD~~t~,~)~u~~J~~Ilr2)+L(~uD~~~3~3)*u(l\il~i1~3~) 

UP(I~~~L;E~~I)~=~~~(.LUDE~~I)+~L*(I~I~(NODE~*I~~)*D(~\~~~UE~*~)+~~I(N~~~~~~*I~~) 

S *D(NODE~~~)+IJII (,JODEs*I *3)*U(NODEb*3)) 

lci0 CcjNTINUE 

IJcj 350 1=1*3 

DIF(~~I~NCJ~=S)= 3P(lrI~N~DES)-U(l*IBNODES) 

350 CurdT I NLt 

wis = 0 . 

DEL=O. 

I.20 360 I=113 

Di) 360 J= 1 rr\130Ea 

DcL=utL+~lF(JII)~~~IF(jri) 

R,u’lS=kMS+UP ( J , 1 ) jkUt-‘( J 9 I 1 

3~0 CLINT I Nut 

i)EL=LUtiT CljEL I 

RlbS = SOtiT ( KdS ) 

TEST=DEL/kkiS 

00 400 I=193 

U(l.IBN~D~~)=JP(lrI~N~D~~) 

400 CUNT I NUt 

IF( TtST.Lt.EPs.Ul~.ITEH.Gk.viAXI )ktl-Ukik 

Gc, TO li; 

END 

156 
Table D-Z 

LJSTING OF SUBROUTINE EQSOJA 

(Continued) 
Table D-2 (Continued) 

(Continued) 

158 
Table D-2 (Continued) 

0 +N(2r2r3)*U(2,3)+N(2~2*4)*24)) 

UP ( 1, 1 ) =UP ( 1 9 1) -Ll*vlI ( 1 l 1*3)*(1~(2*3rl)*U(2r1)+~(2~3*2)~U(2~2) 

3 +N(2r3~3)*U(2,3)+N(2*3*4)*24)) 

UP(lrI)=UP(191) -,~*MI(lrI*4)*(~(2*4*1)*u(2*1)+1~(2*4*2)~-U(2~2) 

tl; +N(2r4*3)*~(2,3)+N(2*4*4)*U(2*4)) 

UP~l*I)=UP~~,1)+~~~~l~lI~~~I~l)~~~l*l)+l~l~l*l~2)~~~l~2)+ 

5 1~lI~1r1*3~*~~1,3)+l”~1~1~1~4~*~~1~4~~ 

UPC29 Iswl2)=( I.-d)*U(2r ltbN1Y2) 

UP ( 2 9 I SI\,.*‘l; ) =cllJ ( 2 , 18N1‘/12) -‘W’tw I ( 2 l 1 9 1 ilWl2 ) -2 ( L. ( 1 9 1 * 1 blwl2 ) +lJ ( 1 * 1 %l\t1’12 ) 

Ic, +L(lrlr 2%1ufi,2 ) ++d ( 1 9 ~LI’v,~‘I~ 1 +L ( 1 9 1 9 jLILe*i2) *b ( 1 * 3~1wl2) +i ( 1 3 1 9 4dh1v12 ) * 

0 u ( 1 , 48rwl2 ) ) 

UP ( 2 9 1 L61\1l”I2 ) CUP ( 2 9 I IplUM ) -‘ti*ih I ( 2 9 I . 1 amI* ) ‘SC (h ( 3 9 1 l 1 aN1v12 ) -zU ( 3 l 1 ~i’d12 2 

s +~(3* 1~23Np12) 3ttJ(3r2tLrw~2)+lu(39 1 r3Bi\llv12)*C,(3r35EulY2)+i’u(3r 1 r4aiw~2)* 

56 U(3r4tbNM2) ) 

~P(~~~BNI~~I~)~U~((~~IN~~I~)--VJ-~~~YII~~~I~~~N~I~~*~L~~~~~~~N~~~~)-~U~~~~~~M~~ 

5 +LC 1r2*2aiu~l2)~~cl1~2~~~i~~2)+L(1~2r3LL;L;v~2)*u(1r3bPuivi2)+~(1r2r4rr~l~~2)~ 

‘b u( 1 r4k+h~2) ) 

cJP( 2 * I rihd2 ) =dP ( 2, I ILi‘wl2) -LV*lJl I (2 * 1 r3druM2)-~(L( 1939 lbNtd2)*U( 1 * 1ON1vl2) 

3 +L( 1 r3*2all\lly2)i%u( I r2b1dly2)+L( 1 *3*3~~1~I~l2)*ll( 1 r3biWl;l)+L( 113rqri,N,vi;?)* 

Lb cl ( 1 ,4BlUlY2 1 ) 

UP(2r I81\1,*12)=c1P( 2,I~i\lM2)-W~hI (29 I *3bl\lw’l2)*(lU(3~~r 1b1\1<12)*U(3r 18NA2) 

‘b +N ( 3 9 3 9 21oNf.12 ) ?%U (3 9 2sb.1~12 ) +Iu ( 3 9 3 9 3:m1*~2 ) *u ( 3 9 3Lolwi2 ) +,A ( 3 13 9 4dw12 ) * 

5 U(3r48NVl2)) 

UP ( 2, I b1\radl2 ) =dp ( 2, I LNlv’l2) --\ni*tm I (2 9 I * 4-l’Jl”12 1 * (L ( 1 9 4. 1 aIJ142) +lJ ( 1 * 18Nk12 1 

5 +L( 1 r4r2kbx~2) -*d L I *2bluM2)+L( 194 r3aNivl2) *u ( 1 9 3brd2) +L ( 1 9 4 9 4r6;~~2 I* 

3 U( 1 9 4Sl’Wlil ) ) 

ut-‘( 2, 1 ~OI\,VQ ) =uP ( 2, I aNlv12 ) -W*I”I I ( 2 9 I * 4bI+12) 9t(1\(3*4r 1lul\llJl~)~u(3* ldJhlr’l2) 

b +N(3,4~~~l\r~~1~)*U~3~2bl~lvl2)+N(3*4’~~~~~”lil) UC, ( 3.3.dw,;i ) +l~ ( 3 9 4 9 4blulv12 j * 

5 U ( 3 9 481\11”12 ) ) 

U~(2rI~,\1,.l2)~J~((,l1,U;Vl2)+b~*(I~II(~*l’1~iUI”l2)*LI(2 . 1 ~r~l~~2) +IVI I ( 2 9 I 9 2olw12 ) 

b sD ( 2,2a,dpl2 ) +,.l 1 ( 2, 1 , 3b1ql*12 ) its ( 2 9 3ardvl2 ) +I’\ I ( 2 l I 9 4hlU142) *lJ ( 2 q 4~&1~(2 ) ) 

UP(l~ODE~,I)=(1.-H)~U(NODEjrI) 

uP(NODEsrl)=u~(~\ri)~E~:,rI) -w-!+& 1 (i\rODEb 9 I 9 1 ) * (L ( Nh 1 9 1, 1 )sU(NMl * 1 ) 

3 +L(Nt’?1 * 1,2)‘~L1(Ni~l1r2)+L(l\ilvll , 193)*u(N1b1 *3)+L(iuM1 l 1*4) 

B *U(NM1 94) ) 

UP(I\IODE~~ I )=Up(juilDEs* I I -~w++l”,I (NUDEbr I r2)*‘(L(lWll r2r 1 )+u(Ivlfil 9 1) 

2.3 +L(NMl , ~,2)+cl~rw11 r2)+L(wd1 r2r3)++U(Nwl r3)+L(NFi1 r2*4) 

B -;cu ( ,\I\1 1 9 4 ) ) 

uc-‘(luolJ~s*l )=dP(19dUlibr I) -wU:[v, 1 ( l\ltiDtb 9 1 93) -X- ( I- ( hl”l1 939 1 )*U(iwll 9 1 ) 

B +L ( P4I.l 1 ,3,~]~u(I‘~~.~1.2)+L(Nf~ll ,3*3)*UCiNlvi1*3)+L(1\11~11 l 3*4) 

3 9UtNMl 94) I 

~C-‘(r\l~u~d, 1 )=,~P(,\1~)3tbq I)-hf*~vlI (I’dU;)t;b* i94)*(L(Iulv11 9 4 1 1 ) *u ( I\ll6 1 * 1 1 

k +L ( i\iW 1 l 4,2) *~Uiw1112)+L(tw11 r4*3)*lJ(NiG1 r3)+L(IUl%1 *4*4) 

(Continued) 

159 
Table D-2 (Concluded) 

3 *U(Nlbil 94) ) 

UP~NODES~I)~U~~,~ODES~I)+~~*~~~I~NODE~*I~~)~D~NODES*~)+MI~NODES*I*~ 

s SD(NUDES~~)+,~II (&UljEbr 1 r3)*D(i\UDEbr3)+lvll (ILODEbw I r4)*D(NODEbr4)‘) 

100 CONTINUE 

DU 350 I=194 

DIF(l*IsNODES)= UP(lrIBNGDES)-U(l*IBNODES) 

350 Cv\T I NUE 

kMS= o . 

DEL=0 l 

Di, 360 I=194 

I30 360 J=l.NL)tiES 

DEL=DEL+C)IF(J,I)~DIF(JII) 

&I~IS=HI~L+UP (J, 1 ) i’tcll=( J, I ) 

360 C3NT I NUE 

D~L=SQ~T.(DEL) 

HIV’ISZSG~HT ( RIv’iS ) 

TEST=DEL/klvlS 

WRITE(6*450)Up,3,DIF 

430 F3f?r~iAT(YE10.3) 

WHITE(~~SOO) 1TtiRqT~STI~EL9H~S 

500 FOHI”IAT( 15*3E10.3) 

DC) 400 I=134 

U(lrlBN~DES)=UP(lrIbN~DE~) 

4uO CCjNTINUt 

IF(TEST.LE.EP~.~l~.ITE~.G~.i~AXl )HETUHN 

GO TO 10 

El\lD 

160 
Table D-3 

LISTING OF SUBROUTINE EQSOL5 

SUBI~OUTINE E~~~L~(U~LIM*NIDINOD~~*W*E~~*MAXI) 

DIMENSIUN U(9,5),L(9rSr5)*M(9*~~5)*1~1(9*5*5)* 

3 N(gr5r5)*D(gv5) ,DET(9)*UP(9*5) l DIF(9*5) 

IdEAL Lrl”l~l”lI ,,\ 

CWll~ON/Xi’rl I /M 1 

l”lI( ~~~~~~NCIDE~)=(I”I(~Y~Y~~~\~UDE~)*I~I(~*~*~~I~~U~~)--IVI( lr3*2%db0kb)* 

b @l( lr2e35NODEb) )*(M( l r4r4rt;hC~DEb)*ld( lr5r5bhODkb)--M( 1*5*48NODEb)* 

8 M( 1.4rSBNODES) 1 

Ml(lrlrlsNODE5)=1~1(l+l*l%NODES) -(M( 1 r2r2WdDk5)*M( 1 r4r38NODEb) 

3 -M( 1 rL+r25iN()DES) *rd( 1 r2r3WdODE5) )*(ivi( l r3*4dYODkb))+Id( l r5r5aNLDEL) 

3 -iY( lr5r48N30ES)*Yl( l r3rSBNODE5)) 

MI ( 19 1 Y ~SINODE~) =,vil ( 1 Y 1 Y 15NUUtb)+(M( 1 r2*2bN0Dta)*M( 195r38NOUEs) 

3 -ly ( 1 9 5 t 2BiJOL)tb ) *VI ( l 9 2 9 3aNUUt: 5, )*(PI( 1 l 3*4bNu~tb)*l”1( 1 r4r5bNuUk5) 

3 -M( lr4*4~N3DE~)*M(lr3~5~i~UD~S)) 

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

b -&I ( 1 ,L+ ( 2aNOL)ib ) WY ( 193r3%NuDEs ) ) * ( PI ( 1 q 2 9 4ducIIkb) *lY ( 1 l 5 l 5aNWEb 1 

YB -MC lr5q48NODES)*M( lr2r58NODES)) 

1”II(lrlrl~N~DE~)=1~~11(l*l~l~N~D~~) - ( Pi ( 1 * 3 9 2blLUDtb ) *A( 1 .5 l 3bNUDEb I 

Lb -,q ( 1 15 Y 2bNODtb ) +,“I ( 1 9 3.3aNbUEs ) ) * (PI ( 1 .2 9 4bNUUtb I *PI ( 1 .4 Y 5bNUDEL ) 

3 -M( 1 r4r4sNODES) +M( 1 r2158NODES)) 

MI(~*~~~BN~DE~)=I*~I(~~~~~~NCJDE~)+(I~I(~*~*~~~~D~~~*I~~(~~~~~BIUUDE~) 

3 -M(~~~*~~NOD~S)~M(~~~*~~N~D~~))*(I~(~~~*~~NOD~~)*M(~~~*~~NODE~) 

3 -lvl ( 1 r3 9 43NOlJES ) *PI ( 1 * 2 * 5rt;NODEb ) ) 

N I ( 1 9 2 t 1 !&NUDES ) = ( 14 ( 1 Y 2 t 1 &NODE5 ) *I% ( 1 9 3 Y 38hUJtb 1 -iv1 ( l 9 3 9 l b;NODkS ) * 

b PI( 1 r2.3Yb1’4UD~b) ) *(l-1( 1 r4.4bNiJDLb)*l”l( 1 l 5*5~i’k~L~)-Pl( 195r4!-f+JUDE2)* 

3 M( 1 .4,58NODES)) 

~~lI(lr2~1BNOD~~)=,~I(l~lr2SNODE~~-(Ivl(lr2*1~l\ri)Dt~~*1~1(1*4r3~N~DEs~ 

b -M( 1 r4r l%tnrOatS)*~‘d( lr2.3bNOtiti5)) *(PJI( 1*3*4%NODks)*M( 1 r5*58NOUE5) 

3 -M( 1 r5-43NODES) *M( 1 r3r53NODES) ) 

I~~I(~~~~~ONUD~~)=,~~(~~~~~BNODES~+(IVI(~~~~~~NCID~~~~I~(~~~~~~~;I~ODE~) 

3 --;v ( 1 Y 5 Y 1 BNoL)L~ ) *IV ( 1 Y 2* 3bNODEb ) ) * (I”1 ( l * 3 9 4W’dOUEb) *Pi ( 1 l 4 * 5aNOL)Eb ) 

3 -M(lr4r48NoDES)*M(lr3~53NODES~~ 

MI( ~~~~~GNL)DE~)=I’II(~.~~~~NODES)+(I”~(I~~*~.~NC)DE~)*IY(~*~~~~~~~DE~) 

3 -M( 1949 1SNOD~h;-) *M( 1 r3r3SNODEb) )*(l”l( 1 l 2*4b\r0Dkb)*l”i( 1+5*58NODEa) 

3 -M(l.5r48NODES)*M(lr2r58NODE5)) 

MI(l~2~l3NOD~S)=MI(l~l~2%NODE~~- (M(lr3*l~NODE~)*ld(l*5*33NCIDE~) 

3 -M( lr5rl~tiNO~~~)*Pl( 1~3r3bNOUE5) )*(14( 1*2’46NUDEb)*M( lr4r5djhUDE~) 

0 -MC 194949 WDES)*M( 1 r2rSsNODES)) 

l”l1(1*2*1dd’A~DE3) =piI (1 Y 1 r28N0DEs)+(l~i( 1 r4* 1bNC)Dta)*w( 1 *5*33NUDE>) 

3 -rq( 1959 lbluOfJts) -liPit 1 r4r3.&NODEb) )*(n( 1 *2*4bhOOiJES)*i~l( 1 r3r52d’JULjEb) 

3 -M(I~~~~~NODES)*I~(~.~*~~N~DES)) 

~lI(lr3rlBNi)DES)=(i~l(lr2rl5NC;DE5)+~Y(lr3r2~i~WUE~)-~~I(lr3~lBNUDE5)* 

(Continued) 

161 
Table D- 3 (Continued) 

16 f’l( 1*2rZbNODES) )*(M( 1 *~*~~NOUES)*IY( 1 r5r5bhODEb)-lq( 1 r594$NoDkb)jt 

8 +‘I( 1 r4953N0DES) ) 

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

3 ~~~(l~4*l~i~O~~S~~~v~(ir2r28NOD~S~~~(~~(lr3~4~l~OD~S)~~l(l~5~5~NOD~S) 

5 --M ( 1 Y 5 9 4aluOL)ES ) +fi\‘i ( 1 Y 3,5ZblUODES ) ) 

~1(l~3*1~N~D~~)=,~I(l*l~3~NODES)+(M(l.2~1~NOD~5)*,~(l,~,2~N~DES) 

3 -M( 1*5*13lWDES)*i~( 1*2*2BNODES) )*(M( Ir3r4&NOUES)++M( lr4r5h~ODES) 

5 -M(lr4~48NO~~S~*M(l~3~5~NOD~~)) 

1~11 ( 1 l 3* l%NODES) =tdI (1 9 1 r3bl\lUUEs)+(W1( 1939 l&i\UDt,,)*~til( 1 ,4,2~,~uD~~) 

3 --M(l*4~1~~~~~S)Qivl(1*3*2tN~DE5))*t~~i(1r2~4~i\l~Ut~)*~vl(1,~y5~N~~t~) 

3 -M( 195948N00tS) *l”‘i( 1 r2r5bNODES) ) 

~~lI(1r3rl~NOD~~)=t~lI(l~l~38NODES;-(~~i(l~3~l~bf~UDt~)*1~l(l~5~2~~~OD~~) 

3 -Iv( 1 *5* lhl\ODt3) *l”l( 1 r392bN’JDEb) )*(lil( 1 r2~4%PdO~E~)Stlb$( 1 *4,5bNUDE=) 

B -M( 1 *4r4%~NOjsJtS) *M( 1 r2r5sNODES)) 

rviI(l~3*1~i\l~D~~)=,vlI(1*1r3bNODES)+(~l(1r4rl3lu~gts)9,~(1r5r23Nrj~E;r,) 

0 -~~(l*5’l~NO~t~)~~~~(l~4*2~NUD~~))*(l~(l~2~4~NUDt~)~~~~(~,3~5~NOD~S) 

3 -,y( lr3r4~NoDES)*M(lr2r59NODES)) 

~YI( lr4rl~~u~~~)=(~~(1’2’ 18NbDES)*ld( lr3r2~bluOUt~)-lrl( 1*3*lslNuDE~)* 

P ,“I( 1 r2r2bNU~ks) )*(lh( 1 r4~32iWDtb)*l”l( 1 *5*5=NUDt~)--ivl( 1 r5*3sNQDtb)* 

5 I’I( 1 .4.5rbiu3DES) 1 

~~I(lr4rl~~UDE~)=i~I(l*I*4~N~D~S) - ( I./I ( 1 .z. 1 BIuU~EL ) *I’d ( 1 I 4 * 2BNODEL 1 

0 -~(l,4r1~N~i)t;5)~~~(lr2r2~N~~E5))*(lv’l(1*3~~~N~G~5)-~~vl(1*~~5~~~UD~~) 

3 -M( lr5r3sNODES)*M( lr3r5%NODES) 1 

,v,I(l,4~loNUD~~)=~~~I(lrlr4~NOD~S)+(1~l(1r2~l~~UD~~)*~~~(1*5*2~~~OD~S~ 

x2 -~~(l~5rl~~~O~~S)~i~l(lr2~2~NUD~~))*(l~1(l*3~3~NUDt~)*~~~(l*4*5~NODE~) 

Lb -iq( 1*4r3&NOIJES) *M( 1 r3r5sNODES) ) 

~~I(l,4~l~~~D~~)=~“iI(lrl,4~NODES)+(1”1(l*3*l~NCjDtS~*~~~(l*4*2~NODES~ 

3 -M ( 1 ,4 Y 1 %;ioDEb ) *IW ( 1 Y 39 2biuUDEb ) ) * ( lvl ( 1 9 2 * 3~lNOUt~) *PI ( 1 9 5 l 5ai’JCDES ) 

3 -M(~,~~~%NODES)*M(~~~*~~I~ODES)) 

1‘41 ( 1949 laFuUDEa>) =,iiI ( 1 * 1 rLlbi\JbDtb) -(Iv~( 1 ,3rlbhtiDca)*l’h( 1 r5*2%NtitiEb) 

22 -l\1(l,5,l~,~O~~~))t~~~(l,3r2~NODE~)~*(l~(l*2*~~l\luut5)*(.M(1~4*5~NODES~ 

5 -i~( 1 r4r39NODtS) sM( I *2*5bNODES) I 

~I(lr4rl~,~~D~S)=~~~I(l,l~4~N~~~5)+(l~l(l*4*l~NUD~5)*~(l*5*2~~~D~S~ 

3 ~-I~(lr5rl~~~O~~S)+~~(l~4~2~~~~D~S))*(F/~(1~2*3~i~0Dt;5~+~~(1~3*5~NUDES~ 

3 -M(lr3r3~NO~tS)*M(lr2*5~NOD~5)) 

~~I(lr~rlfi;NOD~S)~((IV~(l~2rl~N~D~S)*lVl(l~3*2~NUD~S)--l~(l*3~1~NOD~S)* 

c& !q( 1 ,~,~~N~I-JE~) )+(kl( 1 ,4r3blVODtL)*W( 1 •5*4~NUDt~~)-i’~~( 1 *5*3sNbDE5)* 

3 M( 1 r4.43NODES) ) 

162 
Table D-3 (Continued) 

l”II(l*5*l~NUD~~)=;~~I(l*l~5~N~DES)-(~~l(lr3*l~~~D~~)*1~(1*5*2~N~D~S) 

3 -M(lr5rlBNOD~~)~l~(lr3*23N~DE~))*(1~(1*2*3~N~~~S)*M(l~4*43N~D~S) 

3 -M(lr4*3%NODES)*M(le2r48NODES)) 

MI(l*5*l~N0DES)=~~I(1*1*5~NODES)+(M(l*4*l~N~DE~)*i~(l*5*2~NODES) 

3 -M( 1~5*15IiNODE~)*M( 1*4*23NODES))*(M( 1*2*3~NODES)*M( 1*3*4BNODES) 

3 -M( lr3*38NODES)*M(l*2*43NODES)) 

I~l1(f*1*2~N0~E~)=~(rvl(1*1*231\100E~)*I~(1*3r~3N~DE~)-l~(1*~*2~N~DEb)* 

3 M( 19 lr3aNO~Eb) )*(M( l l 4*4%NODES)*M( 1 l ~*~~NUDES)--PI( l l 5*48NCJDEb)+f 

3 M( 1 r4*5sND~Es) ) 

MI(~*~*~~NODE~)=I”II(~*~*~~NODES~- (Ivl( 1 l 1 l 2CbNODES)*M( 1*4*3gNODES) 

a ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

3 -rL)( 1,5*4SNODES) *M( 1 e3r5sNODES) 1 

i~~I(lrl~2~f;CJOD~S)=~~I~lr2~l~NUDE~)+~I~i~1*l~2~NODt~)*~~~lr5r33NL)DE5) 

3 -r~(l*5*2sNO~t;~)*1~l(l*l*3~N~DES))*(1~l(1*3*4~lU0DtT5)*l~l(1*4*5~NuDES) 

3 -,v( 1 r4*43NODliS)*Mf I l 3*58NODES) ) 

MI ( I l I r2bNODES)=,dI ( 1 r2r I~NODES)+(M( I l 3*2~l\t0Dtb)*r~l( I l 4*33NODES) 

5 -,y( lr4*23NODth)*M( 1-3, 3bNUDEb) )*(lQ( 1.1 l 4sNODES)*M( 1 r5r5tbNLJDES) 

3 -M(l*5*4!6NOD~S)*M(l*l*53NODES)) 

MI(l*l*2%d’;)D~~)=,SI(l*2*l~Ni)DE~)- (@I( I l 3*2bNODta)+1Jl( I l 5*3kbNtiiJEb) 

3 -M(1*5*2~~~O~cS)~lVi(l*3*3~N~DES))~(l~(l*l*4~i~~~~~)*~(1*4*5~N~DE~) 

3 -bq( 1 r4*43NoDES) *M( l* lr5sNODES) ) 

~~lI(l*l*2~NOD~~)=~~II(1*2*1~NODES)+(i~l(1*4*2~~~Dt;S)*~~l(1*5*3~NODES) 

Lb -M ( 1 * 5 l 2aNoDt;.s ) *PI ( I * 4 * 3aNUDEs ) ) * ( PI ( I . I * 4br\lOUtb) *MI ( 1 93 l 53NUDES ) 

3 -M(I*3*4~~0DES)~M(l*I*5~NL)DELi~) 

~I(I*~*~~NC~~~S)=(I~(I*I*I~NOD~S)*I~I(~*~*~~N~UE~~-~~(I*~*IBNUD~S~* 

B IYI(~*~*~~NUDE~))*(~~I(~*~*~~NODES)*CM(~*~*~~~UC)DES)-~*~(~*~*~~~\~C)DE~,* 

3 l*l( 1+4,53NO~Es) ) 

MI(1*2*2oNOD~S)=,~I(1*2*2~NODES) -(MC 1.1 l ~SNODES)*IM( 1 *4*38NODES) 

rs -M(1,4*l~N~~E~)*;~~(l*l*3~NOD~S))*(M(l*3*4~N~D~S)*M(l*5*5~NDDES) 

3 -~~(1*5*4%NO~ES)*M(l*3*5~NUDES)) 

(\,I 1 ( 1 , 2 ,2ib~~~~ ) =l4 I ( 1 ,2, ~AN~DE~) + ( 1’4 ( i 9 1 9 A bNGDtS 1 *c,d ( 1 15 9 3WJLiDES 1 

kb -M ( 1 15 9 1 hi\Ol)tb ) *I”1 ( 1 l 1 r3bNO(~Es) )*C(PI( 1*3*4bNuDts)*l~l( 1 *4*5bNUDES) 

3 -M( l*4*48NODES)++M(l*3*5~NODE.S)) 

iv1 ( 1 ,~,~&NODE~)=,~II (1*2*2bNODES)+(t”l( 1*3* lsNUDtS)*M( 1 l 4*3sNoDES) 

3 -M( 1 ,,q, lk!sNOC)t~) *I’dit 1 r3,38NODEb) )*(l;l( I * I l 4~l’J0DEb)*lti( 1 *5*581\0DEb) 

3 -~q(l*5*4~NO~ES)W”l(1*l*53NODES)) 

MI(~~~~~~NODE~)=I~II(~*~~~~N~D~~) -(kc 1 ,3rlsNODtb)*1~l( 1 *5*3BNCJDES) 

5 rM ( 1 ,5r 1 SNoDES ) *IV ( 1 l .3* 3aNODES) 1 * (P’l ( 1 l 194SJODEb)*M( 1 l 4*58NODES) 

3 -M(lr4r43NODES)W’~(lrIr58NODES)) 

IV11 ( 1 l 2*2!bNODES)=idI ( 1*2*28NODES)+(M( l l 4* ~B~ODE~)*I”I( 1 l 5*38NODES) 

3 -M(~,~,~~NOD~S)*~~I(~*~~~~N~DE~))*(I~I(~ l 1 rL+oNODEb)*b’l( 1 r3*5br‘duDl%b) 

3 -M(lr3r43NO~E~)*M(1~1~5~NODES)) 

~(I(l,3*2rf;Ni)D~3)~(l\/1(l*l.l~NDDES)*lYI(1*3*2~~ODES~-~~(1*3*1~~~~D~~)* 

3 I‘,i( 1, 1 •~~I\~CJDE~) )S-c(M( 1 r4943NOD~b)*lW( 1 r5*5~l\LJDk~)-l~1( 1 *5*43N~JUih)* 

3 IVI( 1 r4r53NGDES) ) 

MI(I*~*~~NODES)=;~I(~*~*~~N~DES) -(M( l*l*lBNUDES)+M( 1*4*23NODES) 

(Continued) 

163 
Table D-3 (Cont,hued) 

3 -~~(~~~*~BNO~~S)*I~(~~~*~~NOD~S~,~~(M(~*~*~~NODES~*~~~(~*~~~~NODES~ 

3 -M( lr5r43NODES) *M( 1 r3r58NODES) 1 

MI(l*3,28NOD~~)=,~I(1*2*3~NODE~)+(M(l*l*l~N~DEs~~l”~(l*5~2~~~DE~~ 

3 -~(1~5119\10~~~)*M(1*1*281\10D~S))*(I’v’l(1*3*48N00~~~*1~1(1r4r5~bl\l~DE~) 

3 -M( ] r4*42=l\lODrS) *MC 1 r3r59NODES) ) 

lilI(l*3r28NOD~~)=1~1(lr2*3~NDDES~+(1~(1*3*l4;NODES~*M(1*4*28N~DES~ 

3 -MC ~r4rl3NO~~S)~ivl~~r~r~3NODES))*~il~lrlr43~ODES~*l~~1~5.53NODES~ 

3 -~(lr5*43NoDES)*M(lrl*53~~DDES)) 

MI ( 1 l 3*23NODEa) =,\I1 ( 1*2*3~NOOE~)- (M( 1 l 3* laluODEb)*~~l( 1*5*231\ltiDEb) 

3 -~(1*5~l~NO~~~)*tiv~(lr3*25NODE~~)~~M(l*l*4~NCUE~)*M(l*4*5~N~DE~) 

3 -M(l~4*48NO~~~)~~l(lrl*53NDDES~~ 

141 ( 1 r3r2h~OD&)=,qI (1*2*3bNUDEb)+(lQl( 1 r4* ld~NODtb)*~~l( 1 l 5r231\10DES) 

5 -M( 1.5. I8NODES)*b’l( 1*4*28NODES) )*(M( 1 l 1*43NODtb)*M( 1*3*58NODES) 

3 -M~l~3*4~NOD~~)*~l~1~l~5~NODES~~ 

MI(I*4*2~NODE~)=~M(1*l*l~NODES~~M(l*3*2~NODES~-M~l*3*l~NODES~~ 

5 M( 191 r25frNO~E~) )*(M( 1 *4*39hCDEb)Wl( 1 *5*5BhODEb)-ld( 1 l 5*33NODES)* 

3 M(lr4rS~NODE.S)) 

MI(lr4.2%+lGDE~)=l~I(l*2*4~NODES~ -(M(l*l*lYaNODE~)*M(l*4*23NODE~) 

3 -~(l*4*l3NO~~~)*i~(l~l*2~NODES~~*~M(l~3*3~NODES~*M~l*5~5~N~DES~ 

3 -M( lr5*33NoDES)*M( Ir3*53NODES) 1 

~~~I(lr4r2~~~D~S)~~~II(l~2*4~NOD~S~+(l~~l*l*l~NU~~~~*M(l*5*2~N~DES~ 

3 -~(l*5*l’~NO~~~)*M(l*l*2~N~D~~)~*(l~~(l~3*3~N~~t~~*l~(l*4.5~N~DE~~ 

5 -M( I*4*3~[~ODES)~M(Ir3,5~~~DES)) 

IMI ( 1 ,4,23NOOgES)=i~I ( 1 l 2*4SNODES)+(M( 1 r3* laNtiDts)*M( 1*4*28NODES) 

3 -M( 1 ,4,1bNCD~5)W;( 1 ,3,2”1\1ti(JEb) )*(M( 1 * 1 l Ljbl”\IDDtb)*I’l( 1 *5*5bNODEL) 

3 -~~l~5*33l~oDES)*i~i~lrl*5~NODES~~ 

MI ( 1 *4*28NODES)=,vlI (1*2*4~NUDES) -(Ml 1 r3r lrt;luODE~.,*lvl( 1 *5*23NODES) 

5 -~(Ir5*1~NO~~~)~M(l*~*2~~~~ES~~*(~~(l*l*3~N~D~S~*M(l*4*5~~~~~S~ 

3 -M(lr4*38NO~~S)*M(l*l~5~NODES~) 

l.lI( 1 *4*2bNODES)=,II (1 *2r48NODEb)+(M( 1 *4*lrLNC’Dts)++M( 1 *5*2sNODES) 

0 -lvj( 1 r5* 13N()DEb) ++1’4( 1 r.qr23NtiDEb) ) -2 (PI ( 1 , 1 9 ~jbC\ODES) *VI ( 1 ,3,58i”&~Eb ) 

3 -M(lr3*3~i\100~S)91\‘(11115~NODES)) 

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

5 Ivl( 1, 1 r23LNODE~) ) *(I”I( 1 *4*3bhODES)*Vl( 1 *5*4bNOUEb)-1~l( 195*38NUDEb)* 

3 I”I( 1 r4r43NODES) ) 

MI ( 1 *5*2aNODEa)=r~lI ( 1 l 2*5bNODEb) -(lA( 1 * 1 * ltL~‘~UDts)*t~l( 1 *4*23NCJDEb) 

3 -~~(lr4*l~NO~ES)+l~i(1~l*28NODES~)-~(M(1*3r3~iVOUtS~*i\ll(1*5*4SN~DE~~ 

3 -M ( 1 * 5*3si~OO~S ) ii-l”1 ( 1 *3* 4WLODES 1 1 

~I(l,5*2~NOD~~)=~~I(l*2*5~~OD~S~+(M~l*l*l5NODE~~~~~~(l*5*23NODES~ 

3 -M(I,5,I81\roDES)~M(I*l,28NC;DES))*(iivl(l*3*3~N~DES)~~~(l*4*4lbN~DES) 

5 -~(1,4*38h00ES)“M(l*3~4~NODES~) 

IV,I ( 1 ,5,2&NdDEb) =,:I1 ( 1 r2r~~d\ICIDEb)+(lVl( 1 r3* lbNbDLb)~h( 1 *4*2~~~UDE~) 

3 -i ( 1 + 4 1 1 aj\(~Dtb ) *I’I L 1 13 + 2aNdDEs ) ) * ( l*i ( 1 * 1 * 3brLl~DEb )*14 ( 1 9 5 9 4-~;l\iuDES ) 

3 -~(lr5*3~~O~~~)*l~l(1*1*4~~~DE~~~ 

141 ( 1 .5*2~1\10DEa) =,dl ( 1 •2rEiWJLUE~) -(pit 1939 1al\ruDks)*w1( 1 r5r23NtibLa) 

(Continued) 

164 
Table D-3 (Continued) 

16 -M(lr5~~3NO~tS)~M(1~3~23N~D~S))*~~~~l~l~3~NODES)*M~l.4~4~N~D~S) . 

3 -M(lr4*3~NOD~S)*M(1~l~43NODES~) 

MI( 1~5~2&NODEb)=i~I (lrzr5bbNoDES)+(1VI( 1~4*181\10D1Lb)*l-rl( 1*5*2351N~DESl 

ii -M( lr5rIBNODiS)*M( 1 r4r2%NODEb) )*(Pl( lrlr3~NUUE5)*t.‘l( 1 r3r4SNUDEs) 

3 -M(~~~~~~NOOES)~I~(~,~,~NODES)) 

lylI(1rl~3~NOD~S)=(M~lr1r~~NODES)*l~i~lr;?r381U~U~~,-~~l~lr~r~43lVODES)* 

B M( 1* lr34NODES) )*(M( 1 r4r43NODES)*M( 1*5r53NCJDES)-M( 1 l 5*431d0DES)* 

4; M(lr4rTBNODES)) 

MI(~~~~~~N~)D~~)=I~~I~~~~~~BNODES~ -(M( 1 t 1 r23NODEa)*M( 1 r4r38NODES) 

8 -M( 1 r4r2!bikODEL) *Pl( 19 1’3kbb~ODES) )*c(l”l( 1*2r4bNUDtS)++1~1( 195*53NUbE5) 

56 -M(1~5r43~00ES)*M(1.2.5~NODES~) 

MI ( 19 1 r3rt;NU~E~i)=j41 ( 1.39 l%NODE5)+(M( 111 r2bNUDkS)*iV( 1 r5933NUDESl 

B -M(~~5~23NO~~~)+M(lrl~~~NUDES))~(l~~~*2*4~N~Dt~~*l~~~~4*5~~UDES~ 

ii -M(1~4*4bNODES)*M(lr2rSBNODES)) 

I”11 ( 1 * 1 *38NUDEb) =;glI (1 r3r lbNu~Es)+(M( 1 r2r25NUlJEb)*M( 1 r4r3aNuDEb) 

32 -1’4 ( 1 9 4 9 23NDDtb ) *WI ( 1 92 9 3bNODES ) ) Ji ( IY ( 1 l 1 t 431u0DEL) *IV ( 1 l 5 * 5sNUDEs 1 

3 -M(lr5r43l~ODtS)~M(l,lr58NOOES)) 

MI(~~~~~~NODE~)=III~~~~~~~N~DES) -(14( 1 r2r2bNuDES)*M( 1 r5r3d.+JC)D~a) 

2b -M( 1 r5r28NODts) *M( 1 r2r33hODEs) )*(lvl( 19 1 r4%NUDLS)*M( 1 *4*5d~\dDEb) 

3s -M(lr4r43NODES)“M(l~I~53NODES)) 

MI(lrlr35NODES) =A1 (lr3r ~ZJNODES)+(I~I( 1 r4r23NODES)*M( 1 r5r33NUDES) 

a ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

ft; -M( 1r2*4bNO~~~)**r”l( lr1.50l\1UDkS) 1 

I~II(~~~~~~NDDE~)=(I~I(~~~~~~NUDES;*I*I(~~~~~L~~\~UUE~)-IY(~~~~~~N~DES)~ 

3 M( 19 1 r3oNUDEa) ) *(MC 1 .4r4bNODtS)+ld( 1 r5*5LbNuLjEb)-1~1( 1 r5*43NODEb)* 

s M(lr4rS*NODEb)) 

MI ( 1 c2r3bNO~E~)=rlI ( 1 r3r2%~NUDEb) - ( M ( 1 9 1 9 1 bNcjDts ) *PI ( 1 9 4 l 33Nl)DEa 1 

3 -M(lr4rl~NOD~s)*M(lrl~3~NUDt~))*(~~(l*2~4~NuDtS)*~1(1*~*5~1~~DE~) 

5 -M( 115943 ~ODES)*~~I(~~~~~~NUDES,) 

b11( 1 r2*33NJDEs) =,qI ( 1 r3.2bNUDES)+(b’l( 1 l 19 ldWJDE~)*Pl( 1 r5r33NUDEb) 

3 -M(~~~~~BNODES)~~I(~~~*~~N~DES))~(I~(~~~*~~~~D~S)*~I(~*~*~~N~DE~) 

3 -M(lr4*43NODES)*M(lr2r58NOOES)) 

MI( 1*2r3SNUDEa) =IyI (lr3*28NODES)+(M( 1 329 ~~NC)DES)*IV’~( 1 .4.33NODEs) 

5 -M( 1 r4*1bi\1ODt;b) +*@I( 1 -2+3bNUDES) )*(+‘I( 1 l 1 r4?bl\rObtS)*tvl( 1 r5rgANODEs) 

B -&I ( 1 9 5 l 43lNODE.5 ) *+‘i ( 1 9 1 9 53NCJDE5 1 1 

I”11 ( 1 r2*Z.$bNUD~a) =A1 ( 1 r3~2bN3~ES) -(lrl( lr2r laNClDtS)*M( 1 r5’33NUDES) 

3 -rvl( 1959 l~l\ro~tb) *wl( lr2r3bNbDES) )*c(M( 1 t 1 r4bNuUt5)*M( 1 l 4,5kbNODEs) 

4; -M( lr4r48NODES)*M( lrlr53NODES) 1 

MI (lr2r33NODEh)=i~I (lr3r23NODES)+(M( 1*4*13NUDES)*M( 1 r5*33NODEb) 

32 -M(lr5~lBNDD~S)*M(lr4~33NODES))*(~(l*l*4~~ODE~)*~Vl(l*2*5~N~DES) 

si -M(lr2r43NODES)*M(lrl~53NODES~~ 

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

3 ~~(~~~~~~NODE~))*(M(~.~~~~NODES)*~~(~*~*~~N~DES)-I~(~*~~~~NODE~)~ 

4; C4( 1 rbr58NO~ES)) 

165 
Table D- 3 (Continued) 

I”l’I(~*3*~8N~D~~)=,~~I(~*3*3~N~~i~,+~FI~1*1*1~~N~D~S))t-l”l~~*5*24;N~DES) 

3 -IV ( 1 .5* l%NODES ) *PI ( 1 * 1 * 2bNC)DEb ) ) *(lvl( 1 r2r.+bl\bDE~j*1’~1( 1*4*5aNCJDEs) 

8 -M ( 1 .4 r43NOL)ES ) *Ivl’( 1 ,2,5!GNOUES I 1 

Ml ( 1*3*38NODEb) =pi I ( 1 * 3*3b1\OUEb I+ ( 1% ( 1 l 2 * 1 r;NODtb ) *M ( 1 l 4 * 2bLcjDEj ) 

3 -~(~*4*~3N~~t~)~M(~*2*2~t\1~DE5~)*(lIvl~1*~*4~~~~D~S;-~l~l~1*5*5~l~ODES) 

3 -M~~*~*~~NOD~S~~I~(I*]*~~N~DES~) 

MI(l*3*3~l~OD~~)=,~lI(1*3*3~N~DES)- (lW( 1 r2r 1aNUDEb)++‘I( 1 *5*23NtiDES) 

s -~(~*~*l~N~D~~)~1~(~*2*2~l~i)DE~))*(1~~1*~*4bi\li)DtS)*~~(1*4*5~l~ljDES) 

s -M( 1*4*4b1\10DES) *l”l( 19 1 r5bNCjDES) 1 

PII ( 1*3*d%NODEd)=~~iI (1*3*3bNUDES)+(i~l( 1 l 4* 1 bNtiDEb)++l’vl( 1 rfjr2SNdDEs) 

!b -M( 1 r5r 1b1uOi)ta) +M( 1 *4*2uNU0Eb) )*(M( 1 * 1 *4~~\ILJUt~)~l’~1( 1 *2*5~l’JU0Eb) 

3 -~~(~*2*48~NO~ES)*M(l*l*5ssNODE~)) 

I~/~((*~*~~~~~D~~)=(~I(~*~*~~N~DES)~J~I(~*~*~~~U~DES)--~~(~*~*~~N~DES)* 

% I’I( 19 lr2bN3DE~) )9(l%( 1 *4*331NODk~)*1~1( 1 l 5*5bh0Dtb)-b’\( 1 l 5*33N~DE>1+’ 

3 +‘I( 1 r4r55NODES) ) 

MI( 1 r4r38NODEb)=1~11 (1 r3*4~NbDES)- (itic 19 19 l~ruuoEs)*~~l( 1 *4*231\ODLb) 

0 -~~(~*4*~~N0~t~)*M(l*~*2~N~Dt~~)~(~l(~*2*3~l~~D~~)~~~(l*5*5~l~~~t~~ 

3 --14( 1 r5*3%NoDES) ++M( 1*2*5BNODES) ) 

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

s -M(~*5*~3NODES)~~~i(~*~*2~~~OOES))*(M(l*2*3~~~~~ES)~~~(l*4*5~N~DES) 

s -M(~r4*3~NODE~)*M(l*2*53NODES)) 

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

B -M( 1 l 4* lki\r(j~i;a) *IY( 1 r2?2bI\lti3Eb) )*(Vi( 19 1 l ~zdlUUtb)*l~l( 1 r5r5oNUDEs) 

3 -M ( 1.5 *3bNoDES ) *lVl ( 1 * 1~ 53NODES 1 1 

vi1 ( 1 r4r3sNODEa) =,d]I (1 l 3v43NODES)- (I%( 1 r2rlsNODE5)*1’4( 1 r5*231\IUDEb) 

b -M( 195, lL~\()L)t;b) ++~y’lI 1 *2*2aNOOEb) )*(l”l( 1 l 1 l 3&NUUEb)*PI( 1 *4*5Bt\ti~JLb) 

3 -i~(l~4r38~00tS)~M(lrl~~~N~D~~)) 

MI ( 1,4*3SNOD~s)=:~I ( ~*~*~SJN~CES)+(FI( 1 l 4* lBNODlL~)*M( 1 r5*23NODEb) 

32 -rq( 115r 1bNoDt.i) +lvl( 1*4*2bNUDEb) )*(IY( 1 * 1 *Lj~lk0UEb)~M( 1 r2r50NGCEs) 

If; -M(lr2.~ri;NOi)tS)*l*l(l,Ir53NODES)) 

PiI( ~*~*C~~NOD~~)=~,~(~*~*~~NODES)*I~~(~*~*~~N~DE~)-~~~(~*~*~BNODES~* 

0 lV(( 1, lr2;blAO~Ei~) )*(fW( 1 l 4*3bNODEa)*i~i( 1 *5*4Gl\ODEb)-iVI( 1 *5*39NC~DE5)9 

3 P’I( 1 r4r43NODES)) 

~1I(l*5*3~NUDE~)=r~I(l*3*5~N~~tS)-(l~~~l*l*l bNLDts)*ivi( 1 *4*2~bNtiDES) 

3 -M(1*4*~8N~1)~~)*l~i(~*~*2~N~D~~))*~i~l(1*2*~~iLCilit~)~1~(~*~*4~~~UD~~) 

3 -M( lr5r3%NoDtS)*M( lr2r4bNUDES) 1 

1Y~(~r~r~~N~~~~)=,~~I(~r3r~~N~~ES)+(l~~(lrl*l~~~~~~)*l~l(l~5*~~N~Lj~~) 

3 -p.q( lr5r loNoDES) *141( 1* 1*2sNODES) )*(iY( 1 r2r3hNUDkS)*M( 1~4r43NC~~Jtb) 

3 -M ( 1 14 +33NoDES ) -YclVl( 1 l 2 * 4sNODES 1 ) . 

IV\1 ( 1 l 5*3bNtiDE>) =i+lI ( 1 *,~*5~1\13DES)+(l’d( 1 *2* ldNODES)*i~l( 1 *4*2bNODEa) 

3 -,&q( 1 l 4* lkb,lqoDts) *14( 1 *2*2~l~UDE~) )Jk(l~l( 1 * 1 *3WJUDtb)*l”I( 1 r594bNLi)Es) 

3 -1~(~*5*391\1oDES)*~l(1.~*4~NODtS)) 

PiI ( 1*5*35N(JDES)=plI (1 l 3*5=bN3DES)- (id( 1 *2* lLbiAGDEb)~I~l( 1 *5*23NODES) 

s -~~(~*~*~~NOD~~)~;~(~*~*~~N~~DE~))*(I~I(I*~*~’~NODES)~IV~(~*~*~~~‘J~)DES) 

3 -MC i*4*381~oDES)*Fl( 19 1*4LNCIDEb)) 

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

(Continued) 

166 
(Continued) 

167 

I Ill1 
Table D-3 (Continued) 

3 -M( 1 ,5q 1~1~~0~~s) ++igl( 1 * 1 r2SNODES) )*(1’4( 1 *2*4aNUDES)*M( 19395”“DES) 

3 ‘M ( 1 ,3 r4LbNOD~S ) *b’l( 1 * 2r 55luODkS) ) 

,%I ( 1 ,3,4&i\u~E~) =piI (1, 4,3~i\lU~~~)+(l”l( 1 *29 lbl~UDtL~)*i~( 1 r3*2°‘\iu~~s1 

8 -M ( 1 13 9 1a1uODtb ) *$I ( 1 l 2 * 2LNUL)E:=, ) ) * (Fi ( 1 l 1 * 4W\lUl~kb) ++I”1 ( 1 * 5 * 5%1’JUDEs 1 

0 -~1(lr5r4~NO~E~;*M(l*l~5~i~~D~S)) 

IV’II( 1*3r4~i~UD~S)=,JlI(1*4*3~i\l~Dt~) - ( IY ( 1 .2 1 1 >ruuDtb ) +I;1 ( 1 * 5 * 2iNi)iJES ) 

B -M( 1 *5*1&iqO/)tb) +Pl( 1 *2*2bNUUkb) )+ (M( 1 Q 1 l 4uuOC)tS)++crJl( 1 r3r5bNL~OLb) 

3 -M( lr3r4bNODES)*M( 1, IrSbNODES) 1 

MI ( 1,3r4SNODE~)=;41 (1*4*3bNODEb)+(WI( 1 *3* lSNODt.b)++l( 1 l 5*25NODES) 

s -14 ( 1 15 9 1 !bNol)t5 ) *,$I ( 1 13, z&NUDES 1 ) * (1’4 ( 1 l 1 * 4biulJLJtL) +lvi ( 1 9 2 * 5bNODEb 1 

3 -I~(~*~*~~NO~E~)*M(~*~*~~~~OD~S)) 

,“~~(~*~~~~I~OD~~)=~~;‘I(~*~*~BN.ODES~*~~LI(~~~*~~N~DE~~-~~(~*~~~BNODES~~ 

B &‘l( 1, 1 r23NODEi) ) *(M( 1 *.3*3:1:,ODEb)Wv;( 1 *~~~L+GI’JUDE.~)-I~I( 1 l 593LdUDEb)s 

3 4 ( 1 1 3 * ~%INODES ) ) 

1’41 ( 1 *~~L+!~JNODE.ZJ)=~‘~~I ‘i ~~~+t.+SNUDtd)--(lh( 1 * 1 * fbr\rUDts)s~#l( 1 *~*~SB~JUUEL) 

5 -M ( 1 93. 1 Lbi’JC>(‘Jtb ) *Iv1 ( 1 * 1 * 2&brYb~~C~) ) * (14 ( 1 9 2 * 3kI’JObtb) *I*‘1 ( 1 q 5 * SbPdUbES ) 

I -M( 1r5*3SNODE’~)‘tM( lr2r~SNODEs)) 

1~11 ( 1*4*4~bNODEb)=1~\1 (1 ,4*4bN0DE;5)+(M( 1 t 1 q lbNODtb)+W( 1 l 5*25NUDES) 

Ya -M( 1 *5* ~~IWJ~L>) -kh( 1. 1 *2+NuDEs) )+(+‘I( 1 *2*3LAUUtL)+lV( 1 r395bkUDEb) 

B -M( 1*3r35NoDES) *i%( 1*2*5&NODkS) ) 

~~~(~r4r4~NOD~~)=~~~~(~r4,4~~O~~S~+(Ir1(~*2*~~[\iUDt~)*1vi(~*3*28l\;ODE~~ 

Lb -iq( 1 r39 lbI‘JODE3) iklW( 1 *~*~LNULIE~) )*(ld( 1 l 1 ~~LIuLLEL)~M( 1 .5rfjLNtiDEL) 

53 -,P~(~*~*~~NODES)~M(~*~*~~NODES)) 

MI( 1*4.4!&NUDEa)=AI (1*4*4kJ’JODES)-(iY( 1*2* lBt’dJDE~)*ivl( lr5r25NODEb) 

Lb -M ( 1 ,5* 1 b1‘40D~> ) *PI ( 1 * Z* 2blUiiDEb ) ) ++(I.I( 1, 1 ,~~~‘~OIJE~)*I’JI( 1 r3r~BNUDE5) 

3 -rY( lr3*38NO~ES)~M(l*l*5~NODE5~) 

i”lI ( 1 *494slqi)DEd)=,v\I (1 9494~l\U~t~)+(l”l( 1 *39 ltNObts)+,‘~‘l( 1 *5*2kbrd(JI)Ea) 

B -~(~r5*~~1uO~~~)~crvl(~r3~2~N~D~5~)*(M(l*~*3~bl’JC)LjtC~~~*l(1*2*5~NOD~~) 

B -M( 1 r2*331NOC)tj)+hd( 1 l 1 r5bNUDkS) 1 

~~I(l*5r4~NOD~~)~(i~l(l*l*l~1~~D~S~*l~l(l~2*2~IVUU~S~‘l~1(1*~*1~N~~~S~* 

S 1’4( 1 * 1 r261\10~E~) ) *((Iv’I( 1 * 3 * 3sNODts I ++I”1 i I v 3 l 4zbiubDts ) -PI ( 1 l 5 9 ~SI\~UDES ) x 

~5 Jvi( 1 r3*4sNODES) ) 1 

Jvi1 ( 1 r5r4bNbDES) =I+II (1*495W\1UDEb)-(r’l( 1 * 1, lSl’~bDtL)*~vl( 1 ,3,2sNODZs) 

5~ -M( 1 r3r l~iuoDES)*YI( 19 1*2~NODkS) )s(M( 1 *2*3h~ur\l~j~j-~)G+JVl( 1 ,5,4$&JuDES) 

B -M(l*5*3~NOD~~)~lM(lr2~4~NOD~S)) 

I”i 1 ( 1 15 14SN0DES ) =,.i 1 ( 1 , 4, ‘SdNo~E5 ) + ( IL’1 (1 , 1 , 1 b[\lC~Dk>) *ii;; ( 1 ,5,2$i\ibaES ) 

I, -Pl( I*59 l~~~OOt~)~~~dl( I*1 r2~PJUDtbl )*(Pl( 1*2*3b1uubka)+t1q( 1 *l*4>iqLUtzz) 

b -M( 1*3*3bNOUt~)++M( 1*2*4kb~ODkbl ) 

MI ( 1 *5*48iNODEb)=I+iI (1 *4*5bN0DES)+(v1( 192* l!+l\rOD~~);;~~l( 1 ,3*2$~bct,~) 

lo -b+l( l*3*l~~~O~t~)~lvl(1*2t2~iU~DES)~*(~i(l*l*~~~~~~~S;)*~~~( 1,5,4&NO~~S) 

5 --P-1( 1*5*32~NODE5) *I”(( 131 r4LNODES) ) 

1~11 ( 1 *5r4siuODEs)=.~lI ( 1 *4*5bNODEb)- (rvl( 1929 l%brbcj[jkL)~,4( 1 r5*25NGDEj) 

0 -M( 1 *5’1a1dODts) ~;-IvI( 1 r2-2aNuDtb) )*(h( 1. 1 *~~N~~~~)~I~~I( 1 ,3,4al\~~k.j 

5 -M(l*3*3~NO~~s)*l~(l*l~4~NOD~~)) 

(Continued) 

168 
Table D-3 (Continiied) 

(Continued) 

169 
Table D-3 (Continued) 

(Continued) 
Table D-3 (Continued) 

b -M( 1 ,4r l~til\oDES).+h’~( 1 r3r2bNLUtb) )*(~“~( 1 , 1 9 3aJi\i)i>Eb) .%lVi ( 1 ,2,4ti;hCIL)E5 1 

8 -M(lr2~3~NO~ES)*i~l~l~lr4BNOOE5)~ 

DET(~~N~~E~)=,~~~(~,~,~SI\JOCE~)*MI(~*~*~~N~DE~) 

3 +~,(l,~,28hODES)*MI(l~~~l~t~~~~~) 

3 +pl ( 1 , 1 ,331\i0DES 1 -KM I ( l 93 9 l ?AODES ) 

5 4 1 +,q( ~,I,~SNODES)*MI (1 9 lr4sNCJDES) 

5 +I~~(~,~.~~~~OD~S)*MI(~*~*~~N~DE~) 

NMl=NODES-1 

NM2=NODES-2 

DO 50 1=1*5 

DO 50 J=lq5 

1’4 I ( 1 9 I 9 JdidODEs ) =I., I ( l * I ~JYNO~E~)/~~T(~BI~~U~~) 

50 C;UNTINU~ 

I TEH=O 

10 CONTINUE 

ITEK=ITEG+l 

Da 100 I=195 

UP(lrI)=(1.-d)*U(lrI) 

U~(l~I)=UP(l,I)-~~*l~iI(l rI~1~“~~~~2rlrl)*U~2~l~+l~~2~l*2)~u~2*2) 

3 +N(2*1*3)~U(2,3)+1~(2*l~4)*U(2~4)+~(2*1*5)~~(2~5)) 

UP(l*I)=UP(1,I)-“~“l’lI(l 9 I *2)9(1\(2r2* 1 )+b(;lc 1 )+I\“(zr2r2)*;:u(292) 

3 +N(212r3)*U(L,3)+1~(2~2*4)*~(2*4)+N(2*2~5)~U(2*5)) 

UP(lrI)=UP(191) -d++PII ( 1. I *3)“(N( 2939 1 )*u(29 1 )+ru(~*3~2))-u(2~~) 

s +N(2*3~3)*U(2,3)+h( 2r3r4)~U(2r4)+i‘~(2*3*r5)3U(2r5)) 

uPclrI)=UP(1,I)-.r~l’~lI( 1 •I*4)-~(,\(2*4rl)*u( 29 I )+rq(294*2)*u(292) 

3 +N(2r4*3)~U(2,3)+~(2~4~4)~U(2*4)+1~(2*4*~)-~~(~~~)) 

UPC l* 1 )=JP( 19 I )-;J*vlI ( 1 * I *5)*(N(2*5* 1 )-1tU(;lr 1 )+id(295~2)9tb(2r2) 

3 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

UPC 1) I )=UP( 1, I )+“J*(v!I (1 *1*1~~~0~1r1~+M1~1*1*2~*~~1~2)+ 

9; 1~1I(l~Ir3)gD(l,3)+,~1I(l l 114)*0(1*4)+MI(l*I15)*~(1*5)) 

UP ( 2 9 18wl2 ) = ( 1. -d 1 xu ( 2 9 I BIUM2 1 

UP ( 2 9 15l;rul”l;! ) =uP ( 2, I aJiWVl2) --vL+l~lI (2 I I 9 l*[ul”l2)* (L( 19 1. l’uldJl2) XU( 1 * lLiUl’l2) 

3 +L ( 1 9 1 l 2bN1v12 ) *u ( 1 t 25bI~iv)2 ) +L ( 1 3 1 9 3$b\r1~12) *Y ( l 9 3eikM2) +L ( l , l q4~iL~q2 ) it 

3 U ( 1 9 4rLhd12 1 +L ( 1 , J, l 5Yd1’12) *cI ( 1 .5Ll\v12) 1 

UP ( 2. I BN&‘l2 ) =LJP ( 2 , I biwl2 ) -u*c,‘v, I ( 2 9 I 9 1 aw’I2) i[- (I\ ( 3 9 1 1 l~lIul”l~) ‘-cl (3 9 1 !bi\11”12 ) 

b +N ( 3 q 1 * 2bN;\;12 ) *u ( 3 12alui’vi2 ) +N ( 3 9 I 9 38Nl~12 ) *b ( 3 9 31bd~i2 ) +h ( 3 9 l 9 40lLw2 ) * 

B U(3*4~N1~12)+~(3,lr5~N1~,2)*U(3*~~Ni~l2)) 

UP(2rIBN,~2)=UP(Z,I~~,~2)- W*MI (29 I r2bd\ih2)*(L( lr2r 14~N1’42)-*lJU 1 q lBiLi~l2) 

3 +L ( 1 9 2 9 25Nf42 ) jtd ( 1 9 2shM2 ) +L ( 1 9 2 9 3d~r;l2 ) *U ( 1 9 3dwi2) +L ( 1 9 2 9 45rqivl2 ) * 

3 U ( 1 .4bNIY2 ) +L ( 1 , 2 l 5sbN1Y12 ) *U ( l . S%i’F\.1”12 ) ) 

UP ( 2 9 I %NM2 ) =UJp ( 2 , I biWI2 ) - JJ*.‘~ I ( 2 9 I .251w12 ) s- (rd(3929 lkrw2)iiu(3r ilN1~2) 

3 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

3 U(3r4~NM2)+N(3,2*58Nlvl2)*U(3~5~N1~2)) 

UP ( 2 * I BiWl2 ) =up ( 2 9 I bNM2 ) -h*NI (2* I *3NuPl2)++:-(L( 193. lAku1~12)-~UJ( 1. 1alWl2) 

(Continued) 

171 
Table D-3 (Continued) 

3, +L( l l 3*2&l’J,&q2)*U( I r218luM2)+L( I *3*3~Nb’i2)*U( l *3tlL1~2)+L( l *3r4*lkl”l2)* 

B U(lr4s~~~12,+L(1,3*5~NM2)*U(l*5~NiV12)) 

uP(2r I3N1~12,=bP( 2,1SNiY2)-W*clYI (29 I l 3~Nl”l2)*(N(3*3~lbNl~2)~U(3~ 13NM2) 

& ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

3 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

UP ( 2 9 I8wd2 ) =JP ( 2 3 I Ln’1vl2 ) --id*l’vi I ( 2 l I 9 4aNivl2) * (L ( I 9 4 * 1 SNlri2 ) *U ( 1 9 1 LNM2 ) 

3 +L( 1 qL+r23NM2)*U( 1 *2Sl\K2)+L( 1 r4r3bNli12)*U( 1 *39V12)+L( I r4r4*huld2)* 

B U( I *4lil’wi2)+L( 1 v.$r5SNM2)*U( 1 r5bKlb2) ) 

UP(29 I~iwl2)=d?( 29ISi\i+2)-bv*ClI (29 I l 4bNlh2)*(h(3*4r1uNh%2)*U(3~ 18Nb12) 

B ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

3 U ( 3 9 45hvi2 ) +iV ( 3 9 4 9 5aNp12 1 *U ( 3 9 5~Nl~‘12 ) ) 

UP(~~I~N,V~~)=JP(;Z,I~I~I~I~)--*~*~~II (29 I r5ww12) -;;-(L( I-59 1wh1~12)x-u( 19 isrw2) 

lu +L ( 1 9 5* 23lUlvi2 ) *u ( 1 9 2~l\rl~l2 ) +L ( 1 v 5 9 3m\rt~ ) *U ( I 9 31Liw12 ) +L ( I 9 5 t 4u1utv12 ) * 

3 U ( 1 ,481\11’d2 ) +L ( 1 , 5. 5aN1’412 ) *U ( 1 9 5LNwl2 1 1 

uP(2r 18111d2) =clp’( 29 IaN~vl2)-~+1”11 (29 I *5%1\ll”l2) *(lu(3rtjr ltl\rlvi2)~3(3* lrNr*12) 

3 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

3 U(3~4~l’wl2)+l\1( ~,~*E~~~NI~;~)*U(J*S~NM~)J 

UP ( 2 q I 81~1b2) =tiP ( 2 9 I w162) +h* ( I*‘( I (2 9 I . 18w12 ) KD ( 2 9 1 w.bl2 I +li/l I ( 2 9 I 92~hv~2 I 

e *D ( 2 * 2~~vh2 ) +,>I I i 2 t I 9 3tialv12 ) *;-D ( 2 9 3ww1 2 ) i-b’1 I ( 2 9 I 9 48NM2) i;-D ( 2 * 4ai\bi2 ) 

3 +,&II (2~1*53Nt~2)*D(2*55NM21 1 

UP(NODEbr I )=( 1 .-d)*U(NODEL* I) 

U~(NODESr1)~U~~,\r~~DE~.I)-W*i~I(l~~GE~*ir1);;-(L(Nli.~1 l 19 1 )>FU(NMl 9 l ) 

3 +L ( NM l , 1 ,2)~tillI\li”ilr2)+L(N1~ll , 1 9 3) *U ( NtylI .3 1 +L ( NIM l l 1*4) 

3 +~(NMlr4)+L(1ui~l 9 l *5)*U(l\ll~l1*5) ) 

uP(l\rbDts*I )=U~‘,~~~)DE~,I)--~~*~~~~I(IUJUE~*I~~)*(L(I~I’~~ *;?r 1 )+~U(lu8~llr 1) 

3 +L ( NM 1 ~~,~)~uU(IUIYI~~)+L(N~‘~II l 2~3)*U(~t~iI l 3)+L(NMI l 2*4) 

rb ‘-‘U(KM1*4)+L(l\?l1 r2*5)*U(NFll 95)) 

UP(I~ODES~I)=~~(,~JDES~ I)-w*l~lI (I\IL)DEa* Ir3)*(L(i\ilvlI 939 1 )++U(NiJlI 9 1 ) 

E +L(hlvil 9 3, 2 ) -it U ( ILl*i I 9 2 ) +L ( NM 1 , 3 l 3 ) *b ( iw l t 3 ) +L ( ~I’I l 9394) 

3 “u(NA1*4)+L(~~i~?l 93rES)+U(NMl 95)) 

uP(NoIJEbr 1 )=Up(l~dUEbr I)-Ud*dI (ILoDEs* I r4)*(L(hlvli *4r 1 )*U(Nlv:l. 1) 

3 +L( NM1 9 4. 2 ) *U ( Ni’% I + 2 ) +L ( NM I r4*3)*U(N,viI r3)+L(NMI *4*4) 

3 *U( NM1 94 )+L (IviM r4r5Y*U(NMl *5)) 

UP(NODEbr I )=LJP(;~~DESI I)--w*MI (h0DEbr I r5)*(L(Nlr11 959 1 )-*U(tuidl 9 1) 

3 +L ( IdNl 1 9 5, 2 ) *b ( i\llb; I * 2 ) +L ( NI’~ 1 9 5 9 3 ) +u ( hl”l1 * 3 ) +L (i\IdI 1 *5*4) 

E *U ( NH 1 9 4 1 +L ( 1 didi 1 r59tj)*u(rkb91 95) ) 

UP( NODES, I ) =IJp ( ,~dDLbr I ) +bJ* (VI 1 (l\jODEb* I * 1 )*u(NuLJtsr 1 )+[vlI (NcjDLSic I r2) 

5 *‘D ( ,NCjDc.~ * 2 ) +I’I 1 t ,uUlJCb * I 13 ) *lJ ( ~‘JUDEJ * 3)+lviI (iwtitsr 194)-xD(NODEsr4) 

% +jvlI (NiJtjEb* 1 v5)+:-O(l\rODESr5)) 

lU0 CONTINUt 

00 350 1=1*5 

DIF( 1,101\100~S)= dP( 1, IsNODES)-U( 19 IaNODES) 

(Continued) 

172 
Table D-3 (Concluded) 

350 CONTINUE 

w4s=o. 

DEL=00 

DO 360 I=lrEj 

Di, 360 J=l*NOi)Es 

UEL=DEL+UIF(J,I)*~IF(JII) 

RI~IS=RMS+UP(J~I)*UP(JII) 

360 CcJNTINUE 

DEL=SQHT(DtL) 

kMS=SUKT ( KMS ) 

TEST=DEL/KMS 

WKITE(6*4E l)UP,dvDIF 

450 F~Rl”~AT(YE13.3) 

WRITE(6*5UO) ITERITE~TIDELIKIMS 

500 FORI”IAT( 19q3~10.3) 

DO 400 i=lrg 

U(1~IBN~DE~)=UP(lrl~i~~D~~) 

400 CuNTINUt’ 

IF(TtST-Lte~P~.ti~2. ITEK.GE.MAXI )K~TuKN 

Gu TU 1U 

END 

173 
1. Report No. 2. Government Accerrion No. 

I 

3. Recipient’s Catalog No. 

a-.-. -- 

4. Title and Subtitle 5. Report Date 

December 19 80 

HYPERBOLIC/PARABOLIC DEVELOPMENT 

FOR THE GIM-STAR CODE 

7. Author(s) 

6. Performing Organization Code 

__ __ 

8. Performing Organization Report NO. 

LMSC-HREC TR D697882 

L. W. Spradley, J. F. Stalnaker, A. W. Ratliff 

9. Performing Organization Name and Address 

10. Work Unit No. 

Lockheed Missiles & Space Company, Inc. 

Huntsville Research & Engineering Center 

11. Contract or Grant No. 

4800 Bradford Blvd. NASl-15783, NASl-15795 

_---.-c?_ 13. Type of Report and Period Covered 

2. Sponsoring Agency Name and Address Contractor Report 

National Aeronautics and Space Administration 14. Sponsoring Agency Coda 

Washington, DC 20546 

5. Supplementary Notes 

Langley Technical Monitors: James L. Hunt and J. Phillip Drummond. 

Progress Report 

Appendix A by L. W. Spradley, Lockheed Missiles & Space Company, Inc. 

Appendix B by J;‘rgen Thoenes, Lockheed Missiles & Space Company, Inc. 

Appendix C by John F. Stalnaker, Lockheed Missiles & Space Company, Inc 

Appendix D by S. J. Robertson, Lockheed Missiles & Space Company, Inc. 

6. Abstract 

A study was conducted to compute flow fields in supersonic inlet configurations. The 

problem was run with the elliptic GIM code on the STAR computer. Spillage flow 

under the lower cowl was calculated to be 33% of the incoming stream. The shock/ 

boundary layer interaction on the upper propulsive surface was computed including 

separation. All shocks produced by the flow system were captured. An investigatio! 

of linearized block implicit (LBI) schemes was also conducted to determine their 

application to the GIM code. Pure explicit methods have stability limitations and 

fully implicit schemes are inherently inefficient; however, LB1 schemes show 

promise as an effective compromise. This study also included the development of a 

quasi-parabolic version of the GIM code. The basic idea was to make use of classical parabolized Navier-Stokes methods combined with quasi-time relaxation. This 

scheme is referred to as quasi-parabolic although it applies equally well to hyperbolic supersonic inviscid flows. Second order windward differences are used in the 

marching coordinate and either explicit or linear block implicit time relaxation can 

be incorporated. 

17. Key Words (Suggested by Author(s)) 

Inlet flows 

Shock/bound ary layer interaction 

Linearized block implicit schemes 

Parabolized Navier -Stokes 

18. Distribution Statement 

Unclassified - Unlimited 

Subject Category 34 

19. Security Classif. (of this report1 20. Security Classif. (of this page1 21. No. of Pages 22. Price 

Unclassified Unclassified 175 A08 

For sale by the National Technical Information Service, Springfield. Virginra 22161 

NASA-Lang1 ey, 1980