NASA Contractor Report 175037

Viscous Compressible Flow Direct and

Inverse Computation With Illustrations

(HASA-CB-175037) - VISCOUS COHEBESSIBLE FLOH

DIfiECT AND INVERSE COBPUTATIOB AND

ILLOS1RATIOHS Final Report (Clemson Oniv.)

102 p HC A06/HF A01 CSCL 21E

N86-20391

Unclas

Tah-Teh Yang and Francois Ntone

Clemson University

Clemson, South Carolina

January 1986

Prepared for

Lewis Research Center

Under Grant NAG 3-538

NASA

National Aeronautics and

Space Administration
TABLE OF CONTENTS

PAGE

List of Nomenclature iii

Preface 1

Summary 2

1. Introduction 3

2. Analysis 5

3. Results and Discussion 29

i|. Conclusions 56

5. Recommendations for Future Work 56

6. References 57

Appendix A - Body - Fitted Coordinate System 60

Appendix B - Matrices Used in the Numerical Algorithm 6^

Appendix C - Solution of Block Tridiagonal Systems 69

Appendix D - Boundary Conditions Extrapolation 72

Appendix E - Computer Program 75
List of Nomenclature

A matrix defined on page 13 or variable in turbulence model in page

19.

A matrix defined on page 15.

B matrix defined on page 13.

§ matrix defined on page 15.

Ccp constant in turbulence model on page 20.

^kleb constant in turbulence model on page 20.

Cp specific heat at constant pressure.

Cv specific heat at constant volume.

CW|< constant in turbulence model in page 20.

C2 constant defined on page 7.

D matrix defined in Appendix B.

E matrix defined in Appendix B.

_E matrix defined in Appendix C.

F matrix defined in Equation (1).

FI matrix defined in Equation (1).

F matrix defined in Equation (2).

FI matrix defined in Equation (2).

F matrix defined in Equation (4).

_F vector defined in Appendix C.

Fkleb Function in turbulence model on page 20.

Fwake Function in turbulence model on page 20.

F(y) Function defined on page 20.

G matrix defined in Equation (1)

G-) matrix defined on page 6.

G matrix defined in Equation (2).

(G~1 matrix defined on page 9.

iii

° PRECEDING PAGE BLANK NOT FILMED
G matrix defined in Equation CO.

G-j matrix defined on page 11.

G

1 matrix defined on page 11.

G matrix defined on page 12.

H matrix defined on page 6 or the height between two plates in

Couette flow.

H matrix defined in Equation (2).

H' matrix defined on page (2).

I identity matrix

J Jacobian defined on page 9.

K variable defined in Equation (12).

L matrix defined in Equation (12).

M Mach number

M matrix defined in Appendix B.

MI matrix defined in Appendix B.

N matrix defined in Equation (12).

Pr Prandtl number

Pr^ turbulent Prandtl number

Q matrix defined in Equation (12).

R matrix defined in Appendix B.

Re Reynold number

s spacing defined in Appendix A.

S matrix defined in Equation CO.

T temperature

To stagnation temperature

V velocity, along rpline defined on page 9.

W matrix defined in Appendix B.

iv
J( vector defined in Appendix C.

Y matrix defined in Appendix B.

a coefficient in equation (AH) of Appendix A.

a^ variable defined on page 65.

&2 variable defined on page 65.

a^ variable defined on page 65.

an variable defined on page 65.

b coefficient in Equation (AH) of Appendix A.

c coefficient in Equation (AH) of Appendix A or speed of sound in

page 35.

d constant in Equation (AH) of Appendix A.

e internal energy

f function in Equation (A1) of Appendix A.

f] coefficient in Equation (AH) of Appendix A.

?2 coefficient in Equation (AH) of Appendix A.

g function in Equation (A1) of Appendix A.

g-| coefficient in Equation (AH) of Appendix A.

g2 coefficient in Equation (AH) of Appendix A.

h-| variable defined in Equation (A12) of Appendix A.

\\2 variable defined in Equation (A12) of Appendix A.

k thermal conductivity

k^j constant in turbulence model on page 19.

]<2 constant in turbulence model on page 19. ®

1 the mixing length in turbulence model in page 19.

o order of magnitude in difference equation.

p pressure

po stagnation pressure

Pi

guessed pressure value in Equation (15).
P2 guessed pressure value in Equation (15).

' q matrix defined in Equation (1).

o

q matrix defined in Equation (2).

q matrix defined in Equation

r radial coordinate

r0 radius of circular pipe

rw radial distance at wall Fig. 17,19

t time

t

1 time in transformed system

u velocity component in x direction.

uo inlet velocity for sample calculation.

u average velocity on page 36.

U(jif variable defined on page 20.

v velocity component in r direction.

w total velocity defined on page 35.

x axial coordinate

r Riemann invariant defined on page 35.

A forward difference operator in Equation (6).

V backward difference operator in Equation (6) or gradient of a

function defined in Appendix A.

ic positive constant defined in page 27.

a variable defined in Appendix A.

a-| variable defined on page 65.

a 2 variable defined on page 65.

013 variable defined on page 65.

•

oij variable defined on page 65.

8 variable defined in Appendix A

vi
Y variable defined in Appendix A or specific heat ratio.

Y variable in turbulence model on page 19.

6 Boundary layer thickness

E! dissipation coefficient in Equation (1H)

E2 dissipation coefficient in Equation (TO

C finite difference scheme selector on page 12.

n . coordinate in transformed plane.

6 finite difference scheme selector on page 12.

X constant defined on page 51.

X

I constant on page 51.

y dynamic viscosity

yo constant in turbulence model on page 19.

MJ- turbulent viscosity

\> kinematic viscosity

E; coordinate in transformed plane

p density

p time - average density

Op convergence criterion for wall pressure residual on page

Trr shear stress defined on page 6.

TW shear stress at the wall

Txr shear stress defined on page 6.

TXX shear stress defined on page 6.

TOO shear stress defined on page 6.

<(> dependent variable in turbulent flow.

ij> mass average quantity on page 7.

<J> time average dependent variable in turbulent flow.

<|>' fluctuating quantity in turbulent flow.

vii
fluctuating quantity defined on page 7-

vorticity in turbulent model on page 20.

superscript

d design condition

n the n-th time step

subscript

exit

i

inner

j

outer

r

ref

t

t'

w

exit condition

the i-th point in difference scheme

inner region

the j-th point in difference scheme

outer region

derivative with respect to r

reference

derivative with respect to t

derivative with respect to t1

wall

derivative with respect to x

derivative with respect to n

derivative with respect to £

free stream condition

viii
PREFACE

This work was conducted under a grant NAG 3~338 from NASA Lewis Research

Center. Mr. Ron Steinke is the current project monitor. Mr. Richard Roelke

and Dr. Louis Povinelli were monitors when the project was started. Mr.

Harold E. Rohlik, the chief of the turbine branch, provided us with this

research opportunity. During the course of this research, we had discussion

with Rod Chima, James Schmidt and Leo Donovan, of NASA Lewis Research Center.

Their comments were invaluable. Professor William Joe Rae at State

University of New York at Buffalo and Professor Gino Moretti at Brooklyn

Polytechnic Institute also provided suggestions and encouragements when the

progress was slow. Additional funds from the South Carolina Energy R & D

Center and computer time provided by Clemson University were significant.

All this support made it possible for us to conclude the project.
SUMMARY

An algorithm for laminar and turbulent viscous compressible

two-dimensional flows is presented. For the application of precise boundary

conditions over an arbitrary body surface, a body-fitted coordinate system is

used in the physical plane. A thin<~layer approximation of the Navier-Stokes

equations is introduced to keep the viscous terms relatively simple. The

flow field computation is performed in the transformed plane. A factorized,

implicit scheme is used to facilitate the computation. Sample calculations,

for Couette flow, developing pipe flow, isolated airfoil, 2-D compressor

cascade flow, and segmental compressor blade design are presented. To a

certain extent, the effective use of the direct solver depends on the user's

skill in setting up the gridwork, the time step size and the choice of the

artificial viscosity. The design feature of the algorithm, an iterative

scheme to correct an assumed geometry for a specified surface pressure

distribution, works well for subsonic flows. A more elaborate correction

scheme is required in treating transonic flows where local shock waves may be

involved.
VISCOUS COMPRESSIBLE FLOW DIRECT AND INVERSE COMPUTATION WITH ILLUSTRATIONS

1. Introduction

1.1 Background

In the late 1960's, we extended Stanitz's (1) inverse solution method

for planar potential flows to axisymmetric flows. Curved-wall diffusers were

designed by using the extended method, and they are referred to as Griffith

diffusers (2,3). In Griffith diffusers, the measured wall pressure

distributions correlate very well with the prescribed distributions. In the

early 1980's, we extended the design procedure to include shear flows to

design curved-wall diffusers for short ejectors (M). Even though record high

thrust augmentation ratios were observed, the diffuser pressure distributions

no longer correlated as well as they did in the sixties. The major

difference is that in the Griffith diffusers, there is always a distinct

potential core flow while in the ejector diffusers, the entire flow field is

often dominated by the viscous forces. It is a natural progression that we

began to take on the development of an inverse solution method for flows

without neglecting viscosity. The momentum equation considered here-in is

the Navier-Stokes equation with the thin-layer approximation. The energy

equation is also included as one of the governing equations to account for

the compressibility of the fluid. In the Reynolds stress terms, both

Baldwin-Lomax (5) and Cebeci's (6) turbulence models were considered. The

Baldwin-Lomax model was used in the cascade computations and Cebeci's model

was used in the pipe flow computations. The NASA Lewis Research Center

initially was interested in the possibility of applying this inverse

procedure for gas turbine blade redesign. As the algorithm and computer code

were being developed, its interest was shifted to the applicability of the

code to compressor blade design. Therefore, a compressor blade of known
surface pressure distribution, reported by Schmidt et al (7), was used to

verify the algorithm and the code.

1.2 Outline of the Computational Procedure

The computation begins with an assumed geometry. The flow field

computation, including pressure distributions over the boundaries, is

performed in a transformed plane. The transformation is from a body-fitted

coordinate system. Thompson's (8) general coordinate transformation is used

in conjunction with Sorenson's (9) method to provide the orthogonality of the

grid work at the solid boundaries. An implicit scheme for solving the

compressible Navier-Stokes equation was developed, by the utilization of

Beam and Warming's (10) scheme. The direct solver used in this report

differs fromSteger's (11 ) in that the present procedure allows solutions of

axisymmetric flows and planar flows while Steger's is for planar flows only.

The method of obtaining the desired geometry of either an axisymmetric flow

passage or a blade cascade is achieved by the Secant method (12), and the

virtual velocity method (13). These are progressive, iterative procedures.

Thompkins and Tong (13) used an iterative method in correcting geometries,

for inviscid flows. Modifications and interfacing of the above

techniques were required to yield this comparatively complex procedure. In

this report, we consider planar flows as sub-cases of axisymmetric flows. In

the sub-cases, appropriate terms are deleted from the axisymmetric flow

equations. Flags are used in our computing code to direct the computation to

either axisymmetric flow or planar flow. The listing of the computer code is

given in Ntone's dissertation (14 ) and will be deposited at Cosmic. Requests

for copies of the code can be directed to Cosmic, Computer Software

Management and Information Center, Suite 112, Barrow Hall, Athens, GA 30602.
2. Analysis

In this chapter, the form of the Navier-Stokes equations selected for

computation will be presented. The numerical algorithm selected for their

solution will then be derived. A brief discussion on boundary conditions

will be given here and additional details will be provided in the chapter of

sample calculations. A brief presentation of Cebeci's and Baldwin-Lomax

turbulence models for Reynolds stress calculation will be given. For

improving computational stability at high Reynolds numbers, an artificial

viscosity was introduced. It will also be discussed in this chapter.

2.1 Governing Equations

For a viscous, compressible axisymmetric flow, the Navier-Stokes

equations are written in the vectorial form similar to Peyret and Viviand

(15), as

3F . 3G 1/-3F, . 3G,'

H (1)

where

P

pu

pv

e

pu

pu2 + p

puv

(e+p)u

pv

puv

pv2 +

(e+P)v

G1 = r

0

T

xx

T

xr

_£ $7 -I <""»•> ]*««» '

" 0

T

xr

T

rr

Pr 3r p 2 " +V * UTxr *

k VTxr

"•

vtrr_
H =

lxx

"o

0

p -

_0

"I*1

1

Re

 Tee

3iS _ 2 y 3(rv)

^ ~ 3 r 3r

lxr

3v irr

— C_3jj v_i

~>V V O« »i J

lee

2

T

3u

The above equations are written in the physical plane with cylindrical

coordinates (x, r, t) for axisymmetric flows, where x denotes the axial

coordinate, r the radial coordinate, and t the time. If the multiplication

1 3 factor r is set equal to one, the — terms go to zero, and the -5— terms are

a

replaced by -~— then the sygtem of equations becomes useful for planar flows,

dy

In the above, p is the density of the fluid, u and v the velocity components

in x and r direction, respectively, e the sum of thermal energy CVT and the

kinetic energy ^ (y2 * y2)» P tne pressure, Y the ratio of specific heats and

y the viscosity. The equations are non-dimensionalized using reference

values of length, velocity, density and viscosity. The Prandtl number Pr is

defined by

Pr =

urefcp

kref

where Cp is the specific heat at constant pressure, and kref is a reference

thermal conductivity. The Reynolds number, Re, is defined by
where L, is the reference length. Pressure, p, is obtained from the equation

of state

p = (Y-1) [e - 1/2 p(u2 + v2)]

For laminar flow, the dependence of u on temperature can be accounted for by

using Sutherland's relation:

3/2 C T 2 + ref

"

 Tref C2 * T

where T denotes temperature, and

C2 = 198.6 °R

For turbulent flows, it is customary to express each dependent variable as

the sum of a mean and a fluctuating quantity, such as <j> = <)> + <}>' , and then

time-average the Navier-Stokes equations. For compressible flows, this

procedure leads to the presence of second and third order moments of the

fluctuating variables due to the density fluctuations. To avoid this, the

concept of mass averaging (16-18) is introduced. For example, if $

represents the instantaneous value of a dependent variable, then the

following decomposition is used:

<j> = <|) + <j>"

where $ is a mass averaged quantity, defined as $ = p<|>/p, <}>" is a fluctuation, and the <)> denotes the time averaged quantity. The <t>" is then

related to <J>', which is the customary fluctuation quantity, by

For Navier-Stokes equations, the mass averaged decomposition is applied to u

and v while the customary time averaging is applied to p, e and p. The time
averaged turbulent Navier- Stokes equations are essentially the same as those

for laminar flows, except that one must add to each laminar shear stress and

heat flux term its corresponding turbulent contribution resulting from

fluctuations. Thus, for turbulent flows, a turbulent viscosity y^ ^3 added

to the molecular viscosity u» and the coefficient — - is replaced by

Y (~ + — ), where Pr is a turbulent Prandtl number. However, proper

r » i iV L

If

closure of the equations is not achieved until a method of calculating the

values of ut and Pr^ is introduced. This subject will be dealt with in the

section on "Turbulence Model."

In many problems, one wishes to have boundary conditions satisfied

exactly on an arbitrarily shaped body surface. Therefore, the need for a

coordinate transformation from the cartesian coordinate system (x, r, t) to a

more general curvilinear system (£, n, t1) arises. With the coordinate

transformation

5 = 5(x, r, t)

n = n(x, r, t)

t

f

 = t ,

the axisymmetric Navier-Stokes equations can be rewritten as:

jiq 8? 30" 1 r3F, 3Gn A TT

3f 31 "5^ " Re

 l 35 3nJ

where

q = q/J

F = 5t Q + 5X F/J + Cp G/J
G = n«.q + n F/J + n G/J

t X P

VVJ

H = H/J

J

 - 1/(YV

If one introduces the velocities, V, along n

= constant line, then F and G become:

'pU

constant line and, U, along

F

 = J

G

 = J

where

puU + E p A

pvU + €r P

(e + p) U - 5fcp

PV

puv + nx p

pvM l v + n P r K

_(e + p)v - nfcp

t x r

In the above transformed equations, the presence of the matrices £x,

nx, nr and the Jacobian, J, implies that the coordinate transformation is

known. The details of how the transformation was obtained are given in

Appendix A. For a stationary boundary, ^ and nt

 sfQ zero. For a moving

boundary £t and nt must be specified.

The viscous part, (RHS of Equation (2)), contains a large number of

terms. These items can be greatly reduced if one uses the "thin layer"

9
assumption, according to which £-derivatives in the viscous part are

neglected as small compared to the n~derivative. Such an approximation is

different from that of the usual boundary layer equation assumption, in the

sense that it allows a pressure gradient across the streamlines even inside

the boundary layer. Degani and Steger (19) recently presented a comparison

of calculations using the thin layer assumption and the full Navier-Stokes

equations to validate the approximation. When the thin layer approximation

is incorporated into the governing equations, Equation (2) becomes:

la. + il + 3G . -1

1 9C ~3n Re

(3)

where

and

°, -7

M

0

(

 2 2)u + 1 (u + ) _ 1 1

(

 2 2)y + J_ , + ) _ 2 v

u r/ 2X 2 , 2x 2 „, s T 3 u.. + x -fr C(u ) n ••- (v ) n + 2(uv) n n J ~ 4 y— (un + vn ) o n x i r rixr2r x r

0

0

1 r4 v 2 , ,,

P ~ -5- t-r y — - -=• y(u n +v n,,)] Re3r 3 n x n ^

0

If we let the multiplication factor r be one in Equation (3), 1/r terms

be zero, then H will vanish, and Equation (4) the governing equation for

planar flows will be

10
la + 3t'

A A

JL i§.

Re 3n

(1)

where

A A A

q = q/r , F = F/r , G = G/r

and

V

v + 1

n 3

•?

I(u2

+v

2

):

n ? ? ? ? | (nx

2

 * np

2

) (u2

 + v2

)n

nr + 2(uv)n nx nr]

The axiayrametric equations can also be written in a form very similar to

the planar flow governing equations, and is shown as follows:

A

It'

A A

A

3G

_

Re 3n

(5)

where

A

G.

H' = ± (±

1 1 A

- G')

pv

puv

2

pv

(e+p)v

11
and

'1 - J

-?)

nx)

1

u v nr - - u v

2 v

nx - j y F

In this report examples are presented to show the solutions for

Equations (3) and (4). Equation (5) is suited for more complicated

geometries than equation 3. However, the computing time will be longer when

Equation (5) is used.

2.2 Numerical Schemes

In this section the derivation of the finite difference equations used

in the solution of Equation (3) is presented. Following a similar procedure,

the corresponding finite difference equations for Equations (4) and (5) were

also derived but are not included in this report.

The algorithm selected here is the implicit factored scheme of Beam and

Warming. As a starting point, the time differencing formula is written in

Fade's form as,

(6),

where A and V are forward and backward difference operators, and superscript

"n" represents the nth time step. This is demonstrated by:

.—n —n+1 —n

Aq = q - q

n—n —n —n-1 Vq = q - q

12
and t1

 = n(At') with positive integer of n. Recognizing Aqn~^ = Vqn,

Equation (6) can be rewritten as

Aqn

M

+ 0 (8 - C - -5) (At'2) + 0 (At'3) . (7)

Using Equation (3) we can show that Equation (7) takes the following form:

5!

= Aq" - ^ Aq""1 + (6-C~) 0 (At'2) * 0 (At'3). (8)

2.2.1 Linearization

Equation (8) is non-linear because AFn, AGn, AGn, and AHn

contain terms depending on qn+1 . However, those terms can be linearized,

as shown below for AFn:

F" = r V + r " (F(q)/J)" + L n (G(q)/J)n

U X *

then

AFn . U.

 n+1

 qn+1 + e"

+1

 (F(q)/J)n+1 * £ "+1 (G(q)/J)n+1

I/ X I

- ktq" + 5X" (F(q)/J)n + £r

n (G(q)/J)n} .

From the definitions of F and G, it can be shown that:

F(q)/J = F(q/J) = F(q)

13
G(q)/J = G(q/J) - G(q),

and

3 F(q) _ 3F(q)

3 G(q) 3G(q) . _

= —o,' ' = D.

3q

3q

One may also show from matrix multiplication, that

F(q~) = ACJ and G(q) = B^.

The Jacobian matrices, A and B, are given in Appendix B. Expanding F(q~)

in a Taylor series yields:

pF(q/ c-c~ ) i x A = F(q) + Aq

3q

and since

— n — n+1 -n

Aq = q - q rt,A.0(At

then,

_,—.n+1 _,—.n .n .—n+1 -n. n,tl.2.

F(q) = F(q) + A (q - q ) + 0(At )

0(At2)

similarly,

= B" q"+1 + 0(At2)

lit
then

AFn+1 - [Cn+1 I * Cn +V * Bn] tf * ( ?

n + 1

 - Sn )qn

A similar expression for AGn can also be obtained, if £ is replaced by

Equation (8) can be now expressed as follows:

9At>

 *

~ 1

At' 3Fn tjF 1

1TC

 [I5 aTT Re

(8 - C - -1) 0 (At'2) + 0(At'3) , (9)

where

n-*-1T n+1,n n+10n

. I +n A-t- p B X P

F

n

 = Cd-e)ct

n

 - ee "+1] 71 * [d-e)£

 n

 - ec

 n+1] Fn

(q)

t * t / X X

9 S ] Gn(q)

and

~n r/« «\ n „ n+1, —n r/. ... n - n+1 n _n,—.

G = [(i-e)n.. + en,. ] q + [d-e)n + 9n ] F (q)

t i f X X

C(i-e)nr

n

 * e nr

n+1] cn

(q") .

15
In the present research, the three-point backward implicit scheme is

used, where 9 = 1 and t, = V2. In the direct solver, the boundaries are fixed,

Therefore, the metrics do not change with respect to time. It can then be

shown that

£

n

 T?n „».,* ?n 7^n

F = F and G = G ,

and Fn and Gn no longer depend on (q)""1"

1. One can approximate AG^ in

the following form:

. "T^ri /1 TT "TT \ n . ~~TMn

AG1 = (-j M + MI ) Aq ,

and, similarly,

AH" = (D + E^Aq*1 .

These are linear because M, Mlf D and E are expressed at time n, and

are therefore known, just as in the case of AFn and AGn. The matrices

M!, M, D and E are listed in Appendix B. Then the linearized algorithm

appears as follows:

2.2.2 Factorization

A major gain in computational efficiency may result if Equation (9) is

written in the following factored form:

R.H.S. of equation (10)

16
and

— n+1 — n .— n

q = q + Aq .

Equation (11), the factored equation, differs from Equation (10) only by

terms which do not affect the stated accuracy. Using Equation (11), the

two-dimensional operator of Equation (10) was reduced to two one-dimensional

operators .

For spatial discretization, the central difference scheme is used, and

the terms (-£} and -5— (far*-} are represented by the following approximations

1,J

and

2A£

o(An2),

When AC

 and An are taken to be 1, Equation (11) becomes:

where

1 e(At') -n

2 1+c 1-1,

Qi = I

1 9(At') -n

2 1+c 1-1,J

and

17
2

1 il

The Kjj represents the right hand side of Equation (10), and

are 4x4 matrices. The solution of Equation (12) for Aq^ j requires the

inversion of a block tridiagonal matrix. This procedure is shown in Appendix

C. The solution of Aq's, of Equation (11), which is referred to as the

rrsweep, is obtained from the following:

where

I _ 6At'

 ( _rn __L w

n

 + _L

L

j-1 1H

 ( i,J-1 Re Wi,J~1 Re

M

 6At>

 fRn - — L Y" - — NMn

" 1+c l.J + 1 Re I.J + 1 Re i.J + 1

and

2

 1 J 1 Jmax - 1 •

The matrices W, Y, R are given in Appendix B. For each time step, the

solution of the previous time step is used to calculate the R.H.S. of

Equation (10) and the L.H.S. of Equations (12) and (13). The solutions of

Equations (12), (13) and (11a) yield the flow field at interior points of

the new time step. At the boundaries, either the flow variables are

specified, or they are calculated from the known values at the interior

points. To initiate the calculation, an initial set of values must be

assumed. The three-point backward scheme uses three time levels; the third

level is calculated from the previous two levels. Therefore, one could not

start the calculation with just an initial solution at n=0. The trapezoidal

18
rule is employed for the first time iteration, and the three-point backward

scheme is used thereafter.

2.3- Boundary Conditions

The importance of treating the boundary conditions correctly has long

been pointed out by authors such as Moretti (20). The 1981 Symposium on

Numerical Boundary Condition Procedures (21) was a recognition of the

importance of this aspect in numerical computations. Essentially, boundary

conditions affect both the well-posedness of initial-boundary value problems

and the stability of their numerical solution (22).

At in-flow either the velocity components and the density (or

temperature) or the stagnation temperature and pressure and the flow angle

were specified. At solid boundaries, both u and v are zero and either a

constant surface temperature or an adiabatic wall condition is specified.

Fluid properties which are not specified at boundaries are obtained from

solutions at interior points by extrapolation. At the outflow boundary,

where p is specified, p, u and v are obtained by linear extrapolation.

Details on the extrapolation are provided in Appendix D.

2.4 Turbulence Model

As stated earlier, a complete closure of the governing equations

require that a turbulence model be provided. In this research, two

turbulence models, namely Cebeci's algebraic model and the variation of it

introduced by Baldwin and Lomax (5), were examined.

Cebeci's model (6) is a two-layer model in which the calculation of the

eddy viscosity, p^, depends on whether one considers points in the inner

region or the outer region. In the inner region, the expression for p^ is

based on Prandtl's mixing length theory:

19
'* \ I? I •

where y is the normal distance from the wall and I the mixing length. They

are obtained from the Van Driest expression

with

A+ . 26 v (.IH)"*2 .

Pw

Here, v is the kinematic viscosity, and TW and pw are the shear stress and

density at the wall, respectively.

In the outer region, pt *3 calculated from:

KWer

 = %c

klebp

where

/ (ue-u)dy

o

and

c

kleb

In the above, ue is the velocity at the edge of the boundary layer, 6 is the

boundary layer thickness, Y is the Klebanoff intermittency factor, and k2 =

0.0168.

The inner region extends from the wall to the crossover point yc, where

(Dinner " Wo-

 For

 V > yc» values of (utWter

 are

 applicable. One of the

difficulties in applying Cebeci's model lies in the accurate determination of

6. Baldwin and Lomax (25) pointed out that significant errors would result

20
from the inaccurate evaluation of 6 and modified Cebeci's model to avoid this

calculation.

In the Baldwin and Lomax model, the eddy viscosity for the inner region

is obtained from:

(uthnner - P *2 I » I •

where

„ jta j*y_

3r 3x .

For the outer region, the following formula is used

(ut)outer =

k2 ccp P Fwake Fkleb

where

CCp = 1.6

y F or

= { raax max } use the smaller of the two.

wake Cwk ymax U dif / Fmax

Ckleb ' 0-3 Cwk = 0.25

Fmax is the maximun of F(y) defined as

F(y) .-l|M|l .

and ymax ia the y-location where F(y ) maximum occurs, u^f is given by

udif =

The success of this model relies on the existence of a well-defined peak

for the function F(y) and not having to calculate the boundary layer

21
thickness, 6. In the current cascade computation, the Baldwin-Lornax model

was used in the computer program. No difficulty was encountered.

2.5 Artificial Viscosity

The algorithm described in section 2.2 works well for low Reynolds

number calculations, but often requires additional damping terms for

computational stability at high Reynolds number. The form of the damping

terms used here is the same as Steger's (11). When these damping terms are

added, Equation (11) becomes:

- r J c, V A 2 1+£ n n r'

3F^ . 30^

L Re

In Equation (HO, a fourth order dissipation term has been added to the right

hand side of the equation, while second order terms are added to both

implicit factors. Operators such as V£ and AS are backward and forward

operators. In a central difference approximation,

and

V

2* '

22
are used for second and fourth order derivatives, respectively. The

coefficients, e, and ez, are dissipation coefficients.

Artificial viscosity terms of this nature were studied by Desideri et

al. (22). Their study was based on simplified forms of the governing

equations and certain types of boundary conditions. Nevertheless, their

limits for the ratio of e2 to e, are used as a guideline for the stability

behavior of the present damped algorithm.

2.6 Inverse Design

The inverse design consists of computing a wall geometry that satisfies

a specified pressure distribution. One may initiate the inverse design with

an educated guess-geometry and a reliable procedure for wall modification.

Each newly selected wall geometry should yield a wall pressure closer to the

target pressure than the previous ones. In this research, two methods have

been examined. One is an adaptation of the so-called secant method, often

used in the numerical solution of non-linear equations. The other is based

on the virtual wall velocity idea previously used by Thompkins and Tong (13).

2.6.1 Secant Method

The secant method is illustrated in the following sketch. In this

figure, it is desired to find the zero of a function of <f>(x).

% *z *,

Starting from two initial guesses (x-|,<j>i) and (X2i4>2)

 on the curve, a

straight line is drawn through those two points, intersecting the x-axis at

23
xg. Then another straight line is drawn through (x2t<(>2)

intersecting the x-axis at xjj. If this process is continued, the

intersection of each new straight line with the x-axis should get closer and

closer to the correct value of x for which <J> vanishes. If the two initial

guesses are sufficiently good and the function <j> sufficiently well-behaved,

convergence of the method is quite rapid.

For the extension of the method to inverse design problems, the

following formula is used:

(rv

In Equation (15), all quantities are evaluated at the wall. Subscripts "1"

and "2" denote the two guesses used to calculate the next wall geometry and

subscript "d" the design pressure. At each wall correction, Equation (15)

must be used at every wall point in the region of interest. Clearly, this

means smooth wall geometries will be obtained only if wall pressures are

smooth. Since this may not necessarily be so, a wall geometry smoothing

method must be introduced. In this case, all calculated wall geometries were

smoothed by a least squares method, using a low degree polynomial.

In order to apply the above wall correction scheme to a design problem

using the thin layer Navier-Stokes equations, it is possible to perform a

complete direct calculation on each of the two geometries used to obtain the

next one and repeat the process until a satisfactory geometry is obtained.

However, it seems computationally more efficient to treat the whole problem

as a single direct computation in which wall corrections will be made after a

specified number of iterations (time steps) of the direct solver. One may

start the calculations using one guess for the geometry. The direct

calculations are stopped when the maximum absolute difference between the

21*
current wall pressure and the wall pressure at the previous time step (wall

pressure residual) becomes smaller than a specified value op. At this point,

the second guess-geometry is used, and direct calculations proceed with this

second geometry until the wall pressure residual falls below Op. Then using

the previous two geometries and their corresponding pressures, a new or the

third geometry is calculated using equation (15). This geometry and one of

the two guess-geometries whichever has a pressure closer to the design value

are used to obtain the 4th geometry. The process repeats with 3rd and *Jth

geometries as the first and second guess geometries. For accuracy, a small

value of Op should be specified. However, since Op determines the number of

direct solver iterations between two wall corrections, the Op value may be

increased to speed up the computation. For a given Op, the number of direct

solver iterations decreases as one approaches the target geometry.

Therefore, to maintain accuracy, one may decrease Op as the target geometry

is approached.

The form of Equation (15) makes it suitable for the design of flow

passages such as nozzles, diffusers, etc. where the inlet radius is specified

and remains fixed throughout the calculations. However, at the exit, where

pressure is usually specified, the term P2~Pi vanishes, which makes the

formula inapplicable. This inconsistency affects points in the neighborhood

of the exit, which must be obtained by extrapolation. Here, the

extrapolation was done so that the flow can be made nearly parallel at the

exit by maintaining the exit wall radius equal to the wall radius of the

adjacent point. This, of course, imposes a limitation on the kinds of exit

flows to be considered, but there does not seem to be a general rule of

extrapolation for the exit wall nodal points.

25
Due to the presence of ?2 ~ PI in the denominator in Equation (15), it

is clear that numerical problems will arise if the difference between

calculated and target pressures becomes very small. Therefore, calculations

must be stopped when that difference becomes smaller than a pre-determined

value. Also, initial guesses must be selected such that their corresponding

pressure distribution curves do not cross each other.

2.6.2 Virtual Wall Velocity Method

Referring to the flux vector, G, at the right hand side of Equation

(10) and imposing the no-slip condition at the wall, it can be shown that:

r o

w

V

V

. o

at all times for the direct calculations (for a fixed wall). The subscript,

w, refers to the location at the wall. However, it is possible to imagine a

problem where the wall would be moving towards a steady state position

corresponding to the target (design) pressure. At that steady state

position, Gw would be:

0

V

V

0

where superscript "d" denotes the design conditions. Assuming that the wall

moves with a finite velocity towards its steady state position, u, v, V and

nt are no longer zero at the wall. They can be calculated by requiring Gw

for the moving wall to be equal to G , which leads to:

w

26
/

r

J

\

pv

puv + n P

X

pW + npp

(e+p)V - nfcp

\ /

\ I n

1 - 1 -

\

J

/ " \

0

n p

X

n Pr

0 w

Equation (16) represents a system of four equations. The first equation

is automatically satisfied when the steady state is reached. If

geometric changes in time are disregarded, the last equation yields:

nt - (£ + i) v,

and the second and third equations yield:

\

u = v —

One may then obtain v from the second equation as:

2

n. 1/2

2 A . 2 p

Now from the equation for nt and the definition of V, it can be shown that

V = - | (unx + vnr) .

By applying the rules of partial differentiation, one can show that:

\ " "Vf " W -

Assuming that the wall only moves in the r-direction, then

n

t

 =

 ~

 nrr

t' '

and the magnitude of the wall correction is given by:

2

r\ ~ - d - V2

Ar |= At- (| + 1) ( ( -^ + 1)

n

r

(17)
If the calculations are initiated with an assumed geometry, Equation (17) can

be used to correct the wall at each time step. However, since this would

require the calculation of a new grid at each time step, it is useful to

introduce Op as described in 2.6.1 to set the number of time steps between

two wall corrections. In practice, it may be necessary to multiply the right

hand side of Equation (17) by a constant X. The value of X should be

selected small enough to keep the calculation stable. A value of 0.05 made

the computations very slow and a value of 0.2 was satisfactory. Also, the

computed r's at the wall are smoothed out by adding to the newly computed r's

the expression:

where j=1 denotes the wall, and K a positive constant. In the calculations

performed here, K = 0.2 was used.

28
3. Results and Discussion

In this chapter, some of the cases used to test the algorithm of the

numerical solution described in the analysis chapter are presented. Since

the inverse design feature of the computational program relies entirely on

the direct solution, it is important to assess the capabilities and

weaknesses of the direct solver. Therefore, a great part of the chapter will

be devoted to the results of the direct solver. Then illustrations will be

provided for the inverse design procedure. These illustrative cases have

either known analytical or numerical solutions or experimental data. All

calculations were performed on an IBM 3081 .

3-1 Direct Solution Results

3.1.1 Laminar Flows

3.1.1.1. Flow Formation in Couette Motion

The first case is the two-dimensional flow between two parallel plates

at low Reynolds number. One of the plates is held fixed and the other one is

suddenly set into motion with a constant velocity, uo, at a Reynolds number

—**- of 6.2. Such a flow was considered by Beam and Warming (23), when they

extended their algorithm, originally derived for hyperbolic systems in

conservation law form, to the Navier-Stokes equations. However, their

governing equations were the full Navier Stokes equations in physical

coordinates, while a coordinate transformation and a thin-layer approximation

are used here.

29
j * j „ X / / X // //////// / / ////////// / x

1 - 1 I. I u - a

J ~ / / s s / / / / / 7 / s / f/s/'SS S //{/ . ~

' M.

Figure 1 : Flow formation in Couette motion

The geometry of the problem is shown in Figure 1 . The no-slip condition is

applied at the walls (j = i , and j=jn) along with an adiabatic wall condition.

At the ^-boundaries (1=1 and i=im), periodic boundary conditions are imposed

as follows:

Essentially, the periodic boundary conditions are used here to maintain the

non- dependence of the flow field on the x-direction (or ^-direction) . The

solution of Equation (4) for this case is compared to the exact solution (24)

in Figure 2. The selected grid was uniform and had 6 points in the

{[-direction and 11 points in the n-direction. AS expected, the solution

showed no variations in the x-direction, and only one x-station needed to be

considered. Due to the use of a relatively large time step, only 40

iterations were required for convergence. Because of a large time step, the

computed results of 5 and 10 iterations, labeled as 5n and 10n in Figure 2,

are in poor agreement with the exact transient solution, while the exact

steady state solution and the computed result labeled 40n are in good

agreement. This is an indication that the steady state solution is

independent of the time step and of the initial conditions. Better

30
1.0

0.8

~ 0.6

>

O

o

_J

LU

Exact (Schlichting)

C. U. Numerical

40n

0.4-

0.2-

0.0 0.2 0.4 0.6 0.8 1.0

DISTANCE FROM STATIONARY WALL (y/h)

Figure 2. Flow formation in Couette motion.

4 C. U. Numerical, constant wall temp

o CD. Numerical, adiabatic wall

— Friedman, Gillis, and Liron,

Numerical

2.2

1.8

0.04 0.08 OJ2 OJ6 020

x/(r0-Rero)

Figure 3. Centerline velocity for the entry region of a pipe at Re(j=40.

31
agreement with the exact transient solution is obtained with smaller time

steps, but more iterations are required.

3.1.1.2. Entry region of a circular pipe at Re^ = 40

The second case is the entry region of a circular pipe at low Reynolds

number, with a uniform velocity uo of 100 ft/sec prescribed at the inlet.

Friedman et al. (25) used a numerical solution of the incompressible

Navier-Stokes equations to show how the Reynolds number affects the flow. In

Figure 3, the present solution of the centerline velocity is compared to

theirs for a Reynolds number of 40. The present solution (Eq. 3) was

obtained with a 21x15 grid and converged in MOO iterations from an impulsive

start from uniform conditions. Both a constant temperature wall and an

adiabatic wall were considered. It is seen that the constant temperature

condition yields a better agreement with Friedman et al.'s solution.

3.1.1.3. Entry region of a circular pipe at Re^ = 3000

The third case is the laminar flow in the entrance of a pipe at Re,j =

3000. The uniform velocity prescribed at the inlet was again 100 ft/sec.

Because the entry region at this Reynolds number is quite long (above 300

radii), only part of it was considered in the present calculations. A 41x18

grid was used, and the solution converged in 1200 iterations from an

impulsive start at uniform conditions. The solution is shown in Figure 4,

along with Nikuradse's experimental curves (26) and the analytical results of

Campbell and Slattery (27). Also shown is the only point from Reshotko's

experimental data, as presented in Reference 27, that falls within the range

considered here. Figure 4 shows that the agreement of the present solution

with the solution by Campbell and Slattery (27) is generally good, except at

the centerline. Nikuradse's data are also presented for comparison. The

validity of these data for comparison in this problem, however, has been

32
1.8

— Nikuradse (experimental)

— Campbell and Slatery (Analytical)

& Reshotko (experimental)

o C.U. Numerical (Red = 3000)

0.0 0.01 0.03 OX>4 0.05

x/(r-Rer)

Figure 4. Entry region of a circular pipe. Laminar flow velocity distributions.

o 8 CL Q.

— Schiller (analytical)

o C. U. Numerical (Red = 3000)

0.01 0.02 0.03 0.04 0.05

x/(r0-Rero)

Figure 5. Entry region of a circular pipe. Laminar flow pressure distribution.

33
questioned by Campbell and Slattery, Shapiro (see discussion in Reference 27)

and Sparrow et al. (28). The data point from Reshotko's experiment seems in

good agreement with the present solution near the centerline. Figure 5 shows

a comparison of the computed non-dimensional present pressure distribution

versus the non-dimensional axial distance divided by the Reynolds number with

the analytical results of Schiller reported by Tietjens (26), are in very

good agreement with experimental data. It should be noted that the present

results were obtained using Equation (3) and a rather coarse grid with the

spacing in the x-direction taken as one diameter. For pipe flow problems,

the solution of Equation (3) at high Reynolds numbers seems to allow such

large spacings. However, making Ax much larger than one diameter would

result in a computationally unstable situation.

3.1.1.1. Isolated airfoil

As a test of the algorithm for more complicated geometries, the laminar

flow over an isolated NACA 0012 airfoil was examined. The calculations were

performed at a Reynolds number of 10,000 with respect to the chord length and

a free stream Mach number of 0.2. The grid around the airfoil was obtained

using Sorenson's code (9) which allows specification of a uniform rpspacing

and near-orthogonality at the wall. For the computations, a 79x31 grid of

the C-type, as shown in Figure 6, was used. The distance between the wall

and the nearest £-line was approximately 0.0011, and an elliptic shape outer

boundary was specified at 7 chords away from the airfoil. A zero-degree

incidence angle (measured from the x-axis) was selected, which means symmetry

of the flow field existed about the x-axis. For a geometry such as shown in

Figure 6, the elliptic outer boundary is the inflow boundary, while the

vertical line at the right is the outflow boundary. At the outflow boundary,

the pressure was specified, while p, u and v were extrapolated. The
Ill// / I I I I I I I /

\\\\ \ \ \ \ \

\\\ \ \ V ~

\ \ \ \ \

-3.50 -2.50 -1.50 -0.50 0.50

X

1.50 2.50 3.50

Figure 6. Body-fitted coordinates, NACA 0012 airfoil.

35
extrapolation was done by taking 8/8n as zero, where n is the direction

normal to the outflow boundary. This involves the inversion of three scalar

tridiagonal matrices to obtain p, u and v at the boundaries. At the inflow

boundary, two methods of specifying the boundary conditions were considered.

The first consisted in specifying p, u and v and extrapolating for P by

3P

taking -r— = 0. The second consisted of specifying the stagnation pressure

ori

Po, and temperature, To, and the flow inlet angle, then the left running

Riemann invariant, given by:

2c

where c is the speed of sound and ID the total velocity (u> = /u2+v2), was

linearly extrapolated from the interior to the boundary (see Reference 29)

For adiabatic flows, u at the inlet can be obtained from:

-

(Y-1) r~ + (it(Y+1)CpTo -

where all quantities are evaluated at the inlet. Then u and v are obtained

from w using trigonometric relations and the flow inlet angle. Also, knowing

w and To, the inlet static temperature and speed of sound are obtained. The

inlet Mach number can then be calculated, and the inlet pressure is obtained

from an isentropic relation between P0 and the Mach number. This method is

valid for inviscid flows and only approximately correct for viscous flows.

Therefore, the inlet boundary should be taken far enough from the airfoil.

For the isolated airfoil problem, both methods of treatment of the inflow

boundary were found adequate.

At the wall, the no-slip condition (u=v=0) was specified, along with an

adiabatic wall condition. The wall pressure was obtained from the interior

pressure, which requires the inversion of a scalar tridiagonal matrix. At

36
the common boundary in the wake region, all flow variables were extrapolated

and the averages of the extrapolates from above and below were used. The

extrapolation procedure, again, assumed no gradients in the normal direction.

The results of this calculation are shown in Figure 7, where the

computed pressure distribution is compared to Mehta1s incompressible solution

as presented in reference 11. The present solution, obtained after 500

iterations from an impulsive start from free stream conditions, agrees well

with Mehta1s solution.

3.1.2 Turbulent flows

3.1.2.1. Flow in a circular pipe

In order to assess the validity of the turbulence models presented in

Chapter 2, a turbulent flow in a pipe was considered. Calculations were

performed in a pipe section 5 radii long. A logarithmic velocity profile was

prescribed at the inlet to simulate fully-developed conditions, and the

objective of the test was to verify that such a profile would be maintained 5

radii downstream. This test was performed at a Reynolds number, —, of

110,000, with an average velocity , u, of 100 ft/sec. Both the Cebeci

and the Baldwin-Lomax models were examined, and the velocity profiles

obtained at the pipe exit using both models are compared to the logarithmic

profile in Figure 8. A better agreement with the logarithmic profile is

obtained with Cebeci's model. It was also observed from the computed data

that the Baldwin-Lomax model tended to underpredict the wall shear. For that

reason, Cebeci's model was used in pipe-flow computations.

As a more severe test for turbulent-flow calculations, the flow in the

entrance region of a pipe was considered. The calculations were performed

for a Reynolds number of 388,000 and an inlet velocity of 200 ft/sec using a

Ar

31x3^ grid. The grid spacing near the wall was — = 0.0005, with r

r0 o

37
Mehta (imcompressible, analytical)

o C.U. Numerical (Compressible),

after 500 iterations.

l.O

Figure 7. Pressure coefficient distribution, NACA 0012 airfoil, zero

degree incidence, M=0.2 Re =10,000, laminar.

o

0.0 02 0.4 06 0.8 10

Logarithmic profile

° Cebeci's model

Baldwin-Lomax Model

wall

Figure 8. Fully developed flow in a circular pipe at

38
being the pipe radius. This grid was stretched exponentially towards the

centerline, in order to reduce computational time. In the x-direction, a

uniform increment of one diameter was used. Due to the small grid near the

wall, the time step for this problem had to be taken small enough to allow

stable calculations from an impulsive start at uniform conditions. However,

to accelerate the computations, the time step was increased by a factor of 20

after 100 iterations. The velocity results after 600 iterations are shown in

Figure 9, along with the experimental data of Barbin and Jones (29), obtained

at the same Reynolds number. The agreement is generally good, although the

numerical solution (of Equation 3) indicates a slightly larger acceleration

near the centerline. In Figure 10, the computed pressure distribution is

compared to the experimental one. The experiment was performed on a pipe

length of M5 diameters, while the calculations were performed on a pipe

length of 30 diameters. For comparison, the reference pressure was

arbitrarily taken as the exit pressure of the computations. In this case,

the numerical solution deviates only slightly from the experimental data.

This could partially be attributed to the imperfections of the turbulence

model, and partly attributed to the difference in pipe length between

computational and experimental cases.

3.1.2.2. Compressor cascade

The last case of the direct solution is the flow in a compressor cascade.

The cascade blade geometry used in the calculations was the result of an

analytical and experimental investigation by Schmidt et al. (7). Sorenson's

code (9) was used to obtain the finite difference grid for the cascade. The

grid is shown in Figure 11. It is a 79x31 grid of the C type, in which the

vertical line on the extreme right represents the outflow boundary, and the

vertical line at the extreme left, the inflow boundary. The outflow and

39
— Barbin and Jones, experimental (Red = 388,000)

°1

G } C.U. Numerical (Re, = 388,000)

oj

i £ r/r0 = 0.499

a—s—B r/T0- 0.749

€> r/rc

s

 0.936

30.0

Figure 9. Entry of a circular pipe. Turbulent flow velocity distributions.

1.012

1.008

1.004

— Barbin and Jones (Experimental Red = 388,000)

o C.U. Numerical (Re =388,000)

d

10.0 20.0 30.0

X/d

Figure 10. Entry region of a circular pipe. Turbulent flow pressure

distributions.
-2.00

Figure 11

2.00 3'. 00 14.00

Body-fitted coordinates, compressor blade (ref. 7)
inflow boundaries are connected by two periodic boundaries maintained at a

vertical distance equal to the cascade spacing from each other. It should be

noted that the periodic boundaries are joined to the inlet boundary by

circular arcs. For a clean separation between inlet and periodic boundaries,

it seems desirable to keep the radius of those circular arcs as small as

possible. However, it was found that trying to use a very small radius could

make the grid generation calculations unstable.

The spacing maintained between the inner-boundary and the nearest £-line

was approximately 0.00056. The blunt trailing edge of the airfoil was

changed to a sharp trailing edge in order to allow a smooth stretching of the

grid in the wake region (^-direction), without the need for additional points.

This modification was not expected to significantly affect the solution away

from the trailing edge.

The boundary conditions for the cascade problem were treated as

described in the case of the isolated airfoil, except for the periodic

boundaries. At the periodic boundaries, extrapolation from the interior

•i

(using -7T- = 0) was used for all variables. In order to maintain periodicity,

on

the average values of the extrapolators from above and below were imposed at

both boundaries.

The calculations were performed at a Reynolds number of 170,000, which

is roughly a tenth of the Reynolds number of the experimental investigation

of Reference 7. It was assumed that the solution at the lower Reynolds

number would still represent the essential features of the solution at the

higher Reynolds number. Attempting a high Reynolds number calculation would

have required the use of a finer grid near the wall, and, therefore, a much

larger number of grid points. This, combined with the smaller time steps
associated with finer grids, would have led to a much more expensive

solution.

Convergence for the cascade problem was found to be quite slow, with

large oscillations of the dependent variables in time. In particular, when

p, u and v were specified at the inlet and the pressure extrapolated, large

pressure fluctuations were observed at the inlet. Since those fluctuations

did not seem to die out fast enough, this treatment of the inlet was

abandoned in favor of the specification of the stagnation pressure and

temperature and the flow angle (35.3°). Although oscillations were present

in this inlet treatment as well, they were not as large, and they died out

faster. However, in this case, the inlet Mach number cannot be specified but

is determined from the specified inlet stagnation quantities and outflow

static pressure. For given inlet stagnation quantities, the exact inlet Mach

number will be obtained if the exact downstream pressure is specified. Since

the exact values were not available at the required locations, the

experimental inlet Mach number (0.75) could not be obtained without going

through a trial-and-error procedure.

The results of the calculations shown in Figure 12 are for inlet Mach

numbers of 0.734 and 0.77 measured at the inlet midpoint. Referring to

Figure 12, comparison of the two computed Mach number curves suggests that

for a given upstream stagnation conditions, a small change in downstream

static pressure may result in a small change in inlet Mach number which in

turn may have a significant effect on the solution. Therefore, one may

expect that the non-uniform inlet Mach numberswill affect the solution.

Since a uniform inlet Mach number only occurs when the inlet boundary is

taken far enough from the airfoil, it is reasonable to expect that doing so

might improve the results. Figure 13 shows a linear interpolated result for
Upper Surface

izo C.U. numerical Mtall, -0.77

C.U. numerical Ml(l(tt

s

0734

0.0 1.0

Figure 12. Surface Mach number distributions.

Upper Surface

1.2

1.0

« 0.8

0.6

0.4

0.2

— CU Numerical, linear interpolated

results from M«.734 and M».77

Nasa experimental data ref. 7.

Lower Surface

00 0.2 0.4 06 0.8

x/c.

to

Figure 13. A comparison of surface Mach number distributions between

interpolated numerical results and experimental data.
Figure 12 for inlet Mach number of 0.75. It is compared to the experimental

data. The wall Mach numbers plotted in these figures were obtained from an

isentropic relation between the inlet stagnation pressure and the wall static

pressure. In general, Figure 13 shows that a good agreement exists between

the computed and the experimental values with an inlet Mach number of 0.75.

However, in all cases the computed peak Mach number seems too high.

Over the pressure side, which is the lower surface, the agreement

between computed and experimental Mach number is very good for the first 60$

of the chord. A maximum deviation of 10$ occurs around 80$ chord location.

Over the suction side, the upper surface, deviations occur in the 15 to >40%

chord region. The experimental peak seems to occur at 25$ chord while the

computed peak occurs at 30$ chord. The maximum deviation, however, only

amounts to 7$.

The flow compression takes place rapidly from 30$ chord onward over the

upper and lower surfaces. Some accleration occurs over the lower surface

toward the trailing edge. This flow behavior is shown both computationally

and experimentally. Because of the rapid compression, it is suspected that

there might be a possibility of flow separation over the upper surface toward

the trailing edge. The velocity field plot, shown as Figure 1i», indicates

that there is a small flow reversal and then a reattachment over that region.

Experience gained with the cascade calculations shows that care must be

exercised in obtaining an adequate grid for the computations. The

distribution of points along the airfoil must be such that the geometry is

adequately represented.

Convergence for cascade direct computations may require 3000 to 5000

iterations from an impulsive start. Acceleration of the computations

requires the use of large time steps. The non-dimensional time step At1
-2.00 3.00 4.00

Figure 14. Velocity field between compression blades. Enlargement over

the upper surface near the trailing edge showing flow reversal.
considered here varied from 0.001 to 0.02. However, large time steps require

more smoothing for stability. Since excessive smoothing may lead to less

accuracy, a careful balance must be kept in the choice of the time step and

the smoothing coefficient (referred as EI, in section 2.5). Also, subsonic

computations require less smoothing than computations involving transonic

effects. It should be noted that upwind differencing has not been

incorporated in the code. Therefore, the code is not quite ready for

transonic computations. This may explain some of the discrepancies observed

in Figure 13 where the wall Mach number exceeds 1. On the IBM 3081 and for a

79x34 grid, roughly U.2 seconds of CPU time are required for each iteration.

3.2 Inverse Solution Examples

The illustrations presented for the direct solution were designed to

show that good results can be obtained for a variety of cases. For the

design problem, it is assumed that the direct solver is good, and the success

of the inverse methods presented earlier is judged by their ability to

produce a wall geometry that yields a pressure distribution close to a target

distribution.

3.2.1 Secant Method

An axisymmetric nozzle with a wavy wall geometry as shown in Figure 15

was selected for sample calculations. The corresponding pressure

distribution, as calculated by the direct solver, is shown in Figure 16.

This pressure distribution, along with the two initial geometry guesses shown

in Figure 15, were input to the design calculations. Figures 15 and 16 show

the convergence histories for the wall geometry and the wall pressure,

respectively. It is seen from those figures that a good design can be

obtained with a small number of wall corrections. In this case, the

calculations were started with uniform flow conditions prescribed on the

1*7
•S 100

o

o

o

o

'•5

o

tr

"5

'55

I

'•5

o

Design

1st wall correction

2nd wall correction

3rd wall correction

4th wall correction

00 10 20

Non-dimensional Axial location x/r0

Figure 15. Wall geometry at various interation levels. Secant method

applied to a nozzle at Re<j =24.

40-

£

Q.

i

Q.

•»

c

Cll

o

1

30

0)

'o

"S

8 20

s

3o>

0)

at 1 0

— Design

• 1st guess (iteration 130)

e 2nd guess (iteration 199)

* 1st wall correction (iteration 291)

o 2nd wall correction (iteration)

o 4th wall correction (iteration)

00 10 20

Non-dimensional Axial location x/r0

Figure 16. Wall pressure at various iteration levels. Secant method

applied to a nozzle at re<j=24.
first guess-geometry. The value of Op was kept constant throughout the

calculations, but was selected small enough to allow accurate results. The

W9 direct solver iterations required for the problem are not excessive,

given that the direct solution for a single geometry would require about 200

iterations.

3.2.2 Virtual Wall Velocity Method

3.2.2.1. Flow nozzle

The nature of the virtual wall velocity method restricts its use to

problems in which the region to be corrected lies between two fixed points.

For example, its use for a nozzle design problem may assume that the

locations of the inlet and exit radii are fixed, and the interior part is to

be redesigned.

This method was applied to a nozzle design problem as shown in Figure 17.

In that figure, the target geometry is shown as a solid line, and corresponds

to the pressure distribution shown as a solid line in Figure 18. The

pressure distribution was input to the design calculations as a target

pressure, and the geometry shown as a broken line in Figure 17 was used as an

initial guess. The best geometry, represented by the circles, was obtained

after 15 wall corrections and 670 iterations of the direct solver, which

shows that this method is somewhat slower than the previous method. The

corresponding pressure distribution, shown by the circles in Figure 18, is

quite close to the target pressure.

In order to check if results could be obtained with an initial geometry

taken farther from the correct geometry than in the above case, the same

problem was restarted with an initial geometry as shown in Figure 19. The

best computed geometry, obtained after 22 wall corrections (700 direct solver

iterations), is also shown in Figure 19. Here as before, the corresponding
— Target geometry

— Initial geometry

o Final geometry

0.0 1.0 2.0

Non-dimensional Axial location x/r0

Figure 17. Wall geometry at various iteration levels. Virtual velocity

inverse method applied to a nozzle at Re(j=24.

Of

w

0.

113 40-

«

"S

o

o

o>

h.

3

W

I 0-

— Prescribed

o Computed

Figure 18.

0.0 1.0 2.0

Non-dimensional Axial location x/r0

Wall pressure distribution prescribed versus computed. Virtual

velocity inverse method applied to a nozzle at Re<j=24.
J>

v»

>

'•5

o

cc

o

o

V

o

Target geometry

Initial geometry

Final geometry

00 10 20

Non - dimensional Axial Location x/r«

Figure 19. Wall geometry at various iteration levels. Virtual

velocity inverse method applied to a nozzle at Red=24.

Unmodified blade

Target

Computed (10th wall correction)

o.ej ««>" o Computed (26th wall correction)

2

a> °^

£t

o

o

0.4

0.3-

0.0 0.2 0.4 0.6 0.8 I.O

Non-dimensional Axial Location x/cx

Figure 20. Wall mach number distribution at various iteration levels.

Virtual velocity inverse method applied to a portion of a

compressor cascade. Inlet mach number 0.44.

51
pressure was quite close to the target pressure. In the above cases, A = 0.2

was used and the sign of the R.H.S. of equation 17 was taken as positive if

p

d was greater than p and negative otherwise. More illustration may be found

in reference 30.

3.2.2.2. Redesign of Compressor Cascade Blade

Since the relationship between pressure and wall geometry may become

quite complicated in a cascade problem, the selection of two initial

guess-geometries may not be as easy in this case as it is in nozzle or

diffuser flows. For that reason, the virtual wall velocity method was

selected for use in cascade computations. In order to test the method, the

compressor cascade previously described for direct computations was again

considered, with the inlet Mach number reduced to O.M so that transonic

behavior could be avoided. The wall Mach number from a direct computation on

the unmodified cascade is shown in Figure 20. For the design calculations,

it was decided that the upper surface Mach number "bump", as seen in Figure

20, be removed. The wall correction method was applied only in that portion

of the blade, and the rest of the geometry was kept the same as for the

unmodified blade. Such a modification only serves the purpose of

illustrating the inverse design, and it is understood that in most cases, one

may need to increase the area inside the Mach number curve rather than

decrease it.

For the cascade computations, A was not kept constant, but was

calculated as

d y2

x . v (JL-1J.)

where A1 is specified. The desired effect of this modification was to

emphasize the difference between small and large values of | Ar | calculated

from Equation (17). Particularly, geometry corrections near the junctions

52
between the specified geometry region and the specified pressure region were

to be made smaller, while corrections inside the specified pressure area

(where | P^ - P | is larger) were to be made comparatively larger.

The design calculations were initiated using the unmodified cascade and

its converged solution. After each wall correction, the new geometry was

manually input to Sorenson's Grape code to obtain a new grid. It should be

noted that making the process fully automatic requires coupling of the grid

program with the Navier-Stokes code). After the 6th wall correction, X' was

increased from 0.5 to 0.9 to accelerate the process. The final Mach number

distribution shown in Figure 20 was obtained after 26 wall corrections (1,800

iterations). It is quite close to the target distribution, and could have

been made even closer by allowing for more wall corrections. The

corresponding geometry modification is shown in Figure 21.

Examination of the wall Mach number history after each wall correction

reveals that the wall Mach number in the region being modified behaves in a

wave-like manner. A wavy wall Mach number distribution such as the one shown

in Figure 20 (10th wall correction) is typical. The waves are generated at

the right side of the region being modified and move to the left, where they

tend to vanish. The process will be stable if the amplitude of each new wave

is smaller than the previous one.

The speed of the calculations is controlled by both At' and X'. Since

At' is usually limited by stability requirements for the direct solver,

acceleration of the inverse process is obtained by varying X'. At this

point, the upper limit for X' has not been established as yet.

At the junctions between the region of specified geometry and the region

of specified pressure, little control can be exercised over the pressure

distribution, as shown by the Mach number spike on the left side of the

53
K

u

0)

"5

i 0.3

k_

o

o

o

I

Unmodified blade shape

Final blade shape (26th wall correction)

0.2-]

c

^o

"55

0.1-

•-• 3% of a chord length

O.I 0.2 0.3 0.4 0.5 0.6

Non-dimensional x- coordinate x/cx

Figure 21. A portion of blade shapes before and after modification.
corrected region. Similar behavior has been reported in Reference 31.

Removal of such a spike can be done by independent smoothing of the geometry

or by extending the region of specified pressure.

In the cascade computations, the sign of the right hand side of Equation

(17) was taken as positive for Pd < P and negative for Pd > P. The method

works well for subsonic cases. However, for transonic calculations where

small changes in geometry may lead to substantial changes in flow behavior,

particularly when shocks are present, a more elaborate way of determining the

sign of the wall correction may be required.

55
4. Conclusions

The following objectives of this research were achieved:

1. An implicit factored algorithm was applied to the two-dimensional

(axisymmetric and planar) thin layer Navier-Stokes equations. Good results

were obtained for a variety of cases.

2. Two methods for using the direct solver as an inverse design tool have

been examined. The method of virtual velocity was used for compressor

re-design calculations. In the re-design illustration, 26 wall corrections

or 1*800 time-steps were required to yield the final geometry. The rate of

convergence is considered quite acceptable.

5. Recommendations for Future Work

It is believed that if one gains more experience in using the direct

solver one may obtain more accurate results with less computation time.

Also, improvements can be made for more efficiency and better accuracy. For

example, adaptation of the algorithm to vectorized computers may lead to

substantial savings in computer time. Introduction of upwind differencing

would lead to a better handling of shock waves in transonic flows. For flows

with supersonic inflow or outflow boundaries, conditions compatible with the

supersonic nature of the flow will have to be imposed. Also, higher level

turbulence models can be used, and extensions to the three-dimensional case

are desirable. ' .

56
6. References

1. Stanitz, J. D., "Design of Two-Dimensional Channel with Prescribed

Velocity Distributions Along the Channel Walls," NACA Report 1115, 1953.

2. Yang, Tah-teh, Hudson, William G., and Nelson, Carl D., "Design and

Experimental Performance of Short Curved Wall Diffusers with Axial

Symmetry Utilizing Slot Suction," NASA CR-2209, March 1973.

3. Yang, Tah-teh and Nelson, Carl D., "Griffith Diffusers," Journal of

Fluids Engineering, Vol. 101, No. 1, December 1979.

1. Yang, Tah-teh, Ntone, Francois, Jiang, Tong, and Pitts, D. R., "An

Investigation of High Performance, Short Thrust Augmenting Ejectors,"

Journal of Fluids Engineering, Vol. 107, No. 1, March 1985.

5. Baldwin, B. S., and Lomax, H., "Thin Layer Approximation and Algebraic

Model for Separated Turbulent Flows," AIAA Paper 78-257 presented at

AIAA 16th Aerospace Sciences Meeting at Huntsville, Alabama, January

1978.

6. Cebeci, T., Smith, A. M. 0., and Mosinskis, G., "Calculation of

Compressible Adiabatic Turbulent Boundary Layers," AIAA Journal, Vol. 8,

No. 11, November 1974.

7. Schmidt, James F., Gelder, Thomas F. and Donovan, Leo F., "Redesign and

Cascade Tests of a Supercritical Controlled Diffusion Stator

Blade-Section," NASA TM 83635, 1984.

8. Thompson, J. F., Thames, F. C., and Martin, C. W., "Automatic Numerical

Generation of Body-Fitted Curvilinear Coordinate System for Field

Containing Any Number of Arbitrary Two-Dimensional Bodies," Journal of

Computational Physics, Vol. 15, No. 3, July 197*1.

9. Sorenson, R. L., "A Computer Program to Generate Two-Dimensional Grids

About Airfoil and Other Shapes by the Use of Poisson's Equation," NASA

TM-81198, 1980.

10. Beam, R. M. and Warming, R. F., "An Implicit Finite-Difference Algorithm

for Hyperbolic Systems in Conservation Law Form," Journal of

Computational Physics, Vol. 22, No. 1, Sept. 1976.

11. Steger, J. L., "Implicit Finite-Difference Simulation of Flow About

Arbitrary Two-Dimensional Geometries," AIAA Journal, Vol. 16, No. 7,

July 1978.

12. Forsythe, G. E., Malcolm, M. A. and Moler, C. B., Computer Methods for

Mathematical Computation, Prentice-Hall Inc., 1977.

13. Thompkins, W. T., Jr. and Tong, S. S., "Inverse or Design Calculations

for Nonpotential Flow in Turbomachinery Blade Passages," Journal of

Engineering for Power, Vol. 104, No. 2, April 1982.

57
14. Ntone, Francois, "Direct and Inverse Calculations with Viscous

Compressible Flow Theory," Ph.D. Dissertation, Clemson University,

1985.

15. Peyret, R., and Viviand.H., "Computation of Viscous Compressible Flow

Based on the Navier-Stokes Equations," AGARD-AG-212, 1975.

16. Farre, A., "Equations des Gaz Turbulents Compress!bles," Journal de

Mecanique, Vol. 4, No. 3, Sept. 1965.

17. Rubesin, M. W. and Rose, W. C., "The Turbulent Mean-Flow, Reynolds

Stress, and Heat Flux Equations in Mass-Averaged Dependent Variables,"

NASA TM X-62248, March 1973.

18. Rubesin, M. W., "Numerical Turbulence Modeling," AGARD Lecture Series

No. 86, April 1977.

19. Degani, D. and Steger, J. L., "Comparison Between Navier-Stokes and

Thin--Layer Computations for Separated Supersonic Flow," AIAA Journal,

Vol. 21, No. 11, November 1983.

20. Moretti, G., "Importance of Boundary Conditions in the Numerical

Treatment of Hyperbolic Equations," The Physics of Fluids Supplement II,

1969.

21. "Numerical Boundary Condition Procedures," NASA Conference Publication

2201, A symposium held at AMES Research Center, Moffett Field, Oct.

19-20, 1981.

22. Desideri, J. A., Steger, J. L., and Taunehill, J. C., "On Improving the

Iterative Convergence Properties of an Implicit ApproximateFactorization Finite Difference Algorithm," NASA TM-78495, June 1978.

23. Beam, R. M. and Warming, R. F., "An Implicit Factored Scheme for the

Compressible Navier-Stokes Equations," AIAA Journal, Vol. 16, No. 4,

April 1978.

24. Schlichting, H., Boundary Layer Theory, 6th Edition, McGraw-Hill, 1968,

pp. 83-84.

25. Friedman, M., Gillis, J. and Lison, N., "Laminar Flow in a Pipe at Low

and Moderate Reynolds Numbers," Applied Scientific Research, Vol. 19,

No. 6, Nov. 1968.

26. Tietjens, 0. G., Applied Hydro- and Aeromechanics, First Edition,

McGraw-Hill, New York, 1934.

27. Campbell, W. D. and Slattery, J. C., "Flow in the Entrance of a Tube,"

Journal of Basic Engineering, Vol. 85, No. 1, March 1963.

28. Sparrow, E. M., Lin, S. H. and Lundgren, T. S., "Flow Development in the

Hydrodynamic Entrance Region of Tubes and Ducts," The Physics of Fluids,

Vol. 7, No. 3, March 1964.
29. Bar bin, A. R., and Jones, J. B., "Turbulent Flow in the Inlet of a

Smooth Pipe," Journal of Basic Engineering, Vol. 85, No. 1, March 1963,

30. Proceedings of International Conference on Inverse Design Concepts in

Engineering Sciences, Oct. 17-18, 198*1, pp. 387~396.

•

31. Greitzer, E. M., Thompkins, W. T., Sr., McCure, J. E., Epstein, A. H.,

Tan, C. S., and Hawthorne, W. R., "Current Problems in Turbomachinery

Fluid Dynamics," AFOSR-TR-82-0027, November 1981.

59
APPENDIX A

Body-Fitted Coordinate System

In order to have boundary conditions satisfied exactly over an arbitrary

shaped body, the Navier-Stokes equations are written in a transformed

coordinate system such that the boundaries of the flow field correspond to

grid lines of the finite difference network. In this Appendix, a method of

obtaining such a transformation is explained.

Following Thompson et al. (8), one may define a coordinate

transformation from a physical plane (x, r) to a computational plane (£, n)

so that the computational and physical coordinates are related by the Poisson

equation:

^

 +*

rr

\x

 + V

In the above, f and g are forcing functions to be specified. In practical

problems, it is often desired to define a computational grid work in a

rectangular, uniform, transformed plane. Equations (A1 ) are therefore

rewritten so that £ and n become the independent variables:

a

 \t, ~

 2 e% + Y xnn +

 -

f

\

 + g V - °

Y rnn + 2(frC + g

 V ' ° (A2)

where

_x + r_r , J = = S n ~ 5 n C ri £ n xr r - x r_ *x r r 'x

60
and

£s = J r , £ = - Jx ,

x n r n

nx = ~ J r^ , np = Jx^ (A3)

Equations (A2) are a set of elliptic partial differential equations. If

\

a set of values for x and r are specified at the boundaries of the region of

interest, these equations can be solved numerically by using, for example, an

SOR (successive over relaxation) or an LSOR (line successive over relaxation)

technique .

As mentioned earlier, f and g must be specified. If one selects f = g =

0, equations (A1 ) reduce to a set of Laplace equations, and the solution of

Equations (A2) yields a highly uniform distribution of grid lines as viewed

in the physical plane. If a non-uniform grid is desired, for example when

one needs a fine grid near a wall, which gets coarser as one moves away from

the wall, then one must find a method of specifying f and g. Steger and

Sorenson (9, 11) using the exponential form of f and g suggested by Thompson,

derived a method of specifying the forcing functions so that orthogonality

\

and a specified spacing can be maintained near the boundaries. According to

this method, f and g have the following form:

f(C.n)

 = f U) e"an + f.(0 e"

gU,n) = gO e"bn + g(O ema x , (AH)

where flt flt glf gt are positive constants which must be calculated so that

four geometric constraints are satisfied, and a, b, c and d are positive

constants to be specified.

61
The third and fourth constraints are the same as the first two, except

that they are applied to the boundary j=jmax (outer boundary). They lead to

two equations similar to (10) for x I. . and r I. .

nlj=Jmax nb=jmaxNow considering the boundary j = 1 (n=0), it is clear, since r\max is

large, that Equations (A4 ) reduce to

fU.O) =

g(C,0) = g

Assuming that Equations (A2) apply at the boundary j=1 , one may obtain

expressions for fj(S) and fa(C) by combining (A2) and (A1 1 ) :

(A12)

where

Yxnn)]j=l

In the above equations, x.^ and rn are obtained from Equations (A10)..

A similar procedure applied to the outer boundary allows one to obtain

expressions for f2(£) and g2(£). Then for specified values of a, b, c, d,

f(5,n) and g(C,n) are known functions which can be used in the solution of

Equations (A2). For details on the solution procedure, the rea*der is

referred to Reference 9.

In this research, a computer code was written to solve Equation (A2),

both with f=g=0 and with f and g specified as explained above. The code was

62
tested on simple cases. However, for complicated viscous grids such as for

airfoil cascade calculation Sorenson's Grape code is used (Ref. 9).

63
APPENDIX B

Matrices used In the Numerical Algorithm

In this

sr-Stokes

derived

ces invc

A = -

B = -

From the

A =

appendix, the matrices used in the algorithm for the thin layer

i equations are presented. Equations (12) and (13) of the text

for the numerical solution of Equation (1) of the text. All the

>lved are listed as follows:

_Q I Zl I 0 | Q "

-^ u

2

 + ^- v

2

 (Y-3)u (Y-1)v 1-Y

. : L i

r 1 1

uv -v I -u 0

_ „ i _ . I

p -j ^

Yeu + (1-y) u (l|2 + v2) - Ye + Y~1 (3^2 +y2) (y-l )llV -YU

P P 2

f Q 1_l _ 0 I| -1 I 0 1

 t

uv i -v -u i 0 | |

3~2Y2 2 \ 1-Y 2 rv\— \ ii )nu ( vv T T J jW v I i-YT

|

YPM ? ? YP Y— 1 ? ?

— + (1-Y) V (iT+ v^) (Y-1) uv -

 Y^- ^-L

(3v^+ u^) -Yv

P p 2 |

definition of A and B we can show that

^»

0

 I Cx I 5

—— -_^-*-_ »^.XM.~.«.«......^^.»L^T . J. _

r ° ^...«. « - - -«. -

3— Y ? 1 — Y ? I ii .1 n + lv \ r (-Y •^^ll ^r CY^I ^ u r ( ^\

±y^ 2 2 x x x

"

5

i

v

 !!:v ,

 HrU

 |

-r,

— c uv— E ( u + v ) £. v— E (Y— 1 )u lE u— E (Y— 3)v — E (1— Y)

x r 2 2 x r I x r r

_ ____ _ — __~— ______-i _ ^ __ i._

j t~~ " t

^.Yeu . .2 2v, -Ye Y-1 2 2 I .

^r r |_ ; •—— i i i ^ M| I^UVU^ Y^II^ H +Vv ) ^ JT C, r rL * +^ VTI V jUI * +vV ^l ^ Jl'

 r L f IYT 1 1^ll ^UvV Y^ , T Vr t n, U

x p xp 2 J x-

 x

1 ~1 YPV ? ? —YP Y—1 ? ?

r r - + f1 Y^irfl l +ir |T r ^Y 1^nr r r 4— Y "Jir -Hi 1 HI

r p r r p 2 1

*^

v)
B is obtained from A by replacing £ by n and E, by n .

2 2 Let Ol

 = u (nx - np) 2

fe

a,r

a, r a,r

f— !— 1 + fJLl

J 1J+1 J

a.r a.r

^^T^1 i + ^~T~ J 1 » J J

a0r

l>J

 - (4-).,,., * a r

•i.J + 1

a^r a^r

a- is obtained by replacing the subscript 1 of a to 3

3i,J i,j

a is obtained by replacing the subscript 1 of a to 3

J1.J + 1 i.J + 1

a. is obtained by replacing the subscript 1 of a to 3

^i.J-1 i.J-1

Similarly a^ , a^ and aj. may be obtained.

I,j i.J + 1 i.J-1

65
then

M =

r __ Q ji

.3/JUx _ 3,Ju>

-°3 Trnfy-°2 3^>

... 2 2

3.-J,U +V . Gs,

Q

4 3n r1- p ; 2 ;J

P

2 2

_l(^H_)_a

 3fjv

 )

~

a

2 Tn(r^~)

L-

 c

 - J

«i /<->

-i^)

L 0 l O

3

^ ^ i ^ a

2 9^(?p-) 1 °

^3^ ] °

fa -a ) 3f

JV) a 3f J (01 )

3

 ar9n(

rP

} ai| 3n(

rp}

_i(Ju.) 1

0 i 0 1 0

1 _ _ _ 1

| T

_'_ J i>J , !'J _ I

T 1 (a v + a?u). , -a 1 -a

._ !__:!_ -i-Jiii |.Ji:i

, 2 2 ST t / ^ T ,

Aj ^*'J

 (ai

»

U

-

a

1

u)

l>:j

 (Va

3

v)i,j

••-(a u

2

+a v

2

) ~

(a2v)i,J ~(a

2

u)

1 3 i.J ' 1>J

+2(a2uv)

0

0

I

0

~\i

66
0

-(a3Y + a2u)lj + 1

(

 2+v2

-^)

-(a

'^T +a 2

3v )

-2(a2uv)

r .», , , - . i . i

c.

**i j-1

 p *

-(a u +a_v2)

^ i.J-1

-2(a?uv)

i.J-1

! 0 (.

~

a

!i

2

i,j + 1

<V>1'Jt1

I ° 1

. "" '

|.

0

2

1J + 1'

-a3i>j + i

^/W^j-i

2

 i.J+1

i ° - J 1 _ 1

—

A J-1

1 J~ 1

+ (a u)

i.J-1

I 0

0

L

l o

I

r —

4

i.j + 1

0

0

0

Vj-i
Matrices D, M , E and E. . 1 are given by

i i ,j+ i 1 1 j i

0

0

0

0

0

0

0

0

Y-1

0

H. •= - 1 r

2 v

2 v

3 ReJ,

0

0

_0_

0

2

 nx

2

 nx

0

0

E. . . is obtained by changing . ... of E. to . . ,.

i ,J~i

 1t J*' itJ"1

"' iiJ~ i

0

0

0

0

68
APPENDIX C

Solution of Block Tridiagonal Systems

The use of the Beam and Warming implicit factored algorithm (23) in the

solution of the thin layer Navier-Stokes equations leads to the inversion of

tridiagonal systems. This appendix describes how the solution of such

systems was obtained in this research.

Consider the following system:

Q1 N1

L>2 Q2 N2

L3 Q3 N3

"

 X1 ~

X

2

X

3

x

n-1

_

 xn_

~ F 1

F

2

F

3

F _

_

F

n .

(CD

L

n-1 Qn-1

 Nn-1

as may result from Equations (12) or (13) of the main text. Here, L^, Q^ ,

and MI are 'Jxl matrices, and X^ and Fj are 4 element column vectors.

Equation (C1) is rewritten as:

EX = £ , (C1a)

where "=" denotes a matrix and "-" a vector. A Lower-Upper decomposition of

•

E yields:

69
L

- Q -

 N n-1 n--1

I

i2 I

• •

U1 N1

U2

u3 N3

• •

(C2)

un-1 Nn_.|

where u^ and lj are ^Ix1! matrices, and I is the ^1x4 unity matrix. Carrying

out the multiplication on the right hand side of Equation (C2) and

identifying terms, one obtains:

UT = Qx (C3a)

NJ--,

(C3b)

(C3O

TO
with i = 2, 3, ... n

Now replacing J^ in Equation (C1a) by its Lower-Upper decomposition,

carrying out the multiplication and identifying terms, one obtains:

g1 = fi (Ola)

and g-L = fj - S-i gi--) , i = 2, 3, ... n (CUb)

Xn = un-1 gn (C5a)

Xi = Ui"1 (gi - NiXi+1) , i = n-1 ... 2, 1 (C5b)

Therefore, knowing u^ and lj from Equations (C3) and gj from Equations

(C4), one may calculate the elements of X from Equations (C5). It should be

noted that all required matrix inversions apply to 4x1* matrices and can be

easily performed.

71
APPENDIX D

Boundary Conditions Extrapolation

The following sketch shows the boundary conditions specified for pipe

flows.

at J = Jmax u = v = 0 -I1

 =

P

u

specified

J ~ Jmax

p specified

J = 1

i=1 atj = 1 v = 0 -5T = ° TT ="

For the explicit part, the R.H.S. of Equation (10), at the inflow, the

pressure p is obtained from linear extrapolation:

, . Asi. As-i

where Asi = S2~s-|, A&2 = 33-82 and s is the length along a E, = constant line.

The energy terra, e, is obtained from the knowledge of p, u, v and p. At the

outflow, u, v and p are obtained from an expression similar to the above

extrapolation formula. At the wall, j = jmax» the following backward second

order approximation is used to obtain T and p at wall from interior points.

* »Jmax *»Jmax-1 * »Jmax~2 = •• - i in . -..

• »Jmax

2An

where, -^- = 0 and -r— = 0 for an orthogonal grid,

or) or)

At the centerline, or line of symmetry, v = un = Tn = 0 are specified. Again

= 0. The following second order approximation is used to obtain u, T and
2

2An

0(An

For the implicit part of the algorithm, Equation (13) shows that at j = Jmax~1»

the knowledge of N >m q^ .ffl or a typical term such as §" , Aq^ ,m is required.

At the wall, u = v = 0 and the matrix multiplication of (B Aq ). . i , jm

yields the following:

(Y-1)

For pn = 0 one may approximate Ae^fjm as follows,

= a Ae" + b Ae?

then

(B Aq ). . may be rewritten as

73
ooo o

0 0 0 (nx)i>Jra_1(Y-1)a

ooo o

In such a manner, N ^ lg

Jm

 flqi,Jm

 l3

OOO O

L° o o o

.—n

as a linear combination of

. and Aq1

-1 i,jm-2 .' . _,

 which are the values at the interior

points.
APPENDIX E

Computer Program

The computer program developed for the solution of the two dimensional

(axisymmetric and planar) thin layer Navier-Stokes equations is presented in

this appendix. The structure of the code is discussed, and a description of

the code input variables is given. Output from the code is also described,

and suggestions concerning how to run the code are given.

Code Description

The main routine of the code is very short, and merely allows to select

between Equations (3), (4) or (5) of the text. It contains a parameter

called IAXISY. For IAHSY = 0, Equation (4) is used. For IAXISY = 1,

Equation (3) is used, and for IAXISY = 2, Equation (5) is used. For these

three values of IAXISY, calls are made to Subroutines AXISYO, AXISY1 and

AXISY2, respectively. These three subroutines are similar to each other in

their internal structure: they all have an initialization part, or solution

part, a boundary condition part and a wall correction part.

In the initialization part all constants and variables needed for the

calculations are assigned values. A detailed description of input data is

given below. During the initialization, a call is made to Subroutine TRANSF

for the calculation of the finite difference grid. Also, during the

initialization, initial values are given to the flow variables p, u, v, e, P

(in the code, RH02 (I,J), U2(I,J) , V2(I,J), E2(I,J), P2(I,J)) in the entire

flow field considered.

In the solution part, the right hand side of Equation (1) of the text

(the K^ j's) is formed, with the explicit artificial viscosity terms included.

The right hand side is stored in arrays HR(I,J) , HM(I.J), HN(I,J) and HE(I,J).

Then the £-sweep is carried out: matrices L^, Q^, and N^ (represented by

75
C(I, K1, K2), A(I, K1, K2) and B(I, K1, K2), where K1 and K2 have values from

1 to *0 are constructed, with the implicit artificial viscosity terms

included, and the system of equations represented by Equation (12) of the

text is solved. During the solution of Equation (12), a call is made to

Subroutine TRIDIA which performs block tridiagonal inversions. Subroutine

TRIDIA itself makes calls to Subroutine INVERT, which performs the inversion

of H x 1 matrices. The solution of Equation (12) yields the Aqfj 's,

which are stored in HR(I,J), HM(I,J), HN(I,J) and HE(I,J) for use in the

^-sweep. During the rrsweep, matrices Lj, Qj and Nj are constructed with

implicit artificial viscosity terms included, and a call to Subroutine TRIDIA

leads to the solution of Equation (13) of the text for the Aq^ i 's.

In case the flow is turbulent, a call is made, at the beginning of the

solution part, to Subroutine TURB or SUBROUTINE CEBECI, for the Baldwin-Lomax

or the Cebeci turbulence model.

In the boundary condition part, the values of the flow variables are

obtained at the boundaries according to the guidelines given in the text.

Only the explicit treatment is done in this part, since the implicit

treatment is involved in the formation of matrices L^, Q^, N^ and Lj, Qj, Nj.

In cases where scalar tridiagonal inversions are to be done (as in airfoil

problems as explained in the text), Subroutine STRID is called.

The wall correction part is only used when an inverse design is

desired, and only if the wall pressure residuals become smaller than rp.

When the secant method is used, a call is made to Subroutine SMOOTH, for the

least squares smoothing of the computed geometries.

Subroutine TRANSF is the grid generation subroutine. In this

subroutine, the coordinates of the boundary points are input and the

coordinates of the interior modal points (X(I,J), R(I,J)) are obtained in two

76
ways: by interpolation from boundary points or by solution of elliptic

partial differential equations (Appendix A). In addition, the grid may be

read from a separate file (i.e. when the grid has been calculated from a

different code and stored in a file). It is also possible to obtain a

coarse, but smooth grid by solving the equations of Appendix A and refine it,

say near a solid wall, by a spline interpolation procedure.

Input variables

The following is a description of the input variables, their meanings,

and how they should appear on data cards. A complete description is first

given for the planar case. Then changes related to the axisymmetric case are

shown.

For the planar case, data are input as follows (M, L and all variables

beginning with I or J are right justified integers):

Card 1

Column

1-10

Card 2

11-20

21-30

31-^0

Variable

IAXISY

Column

1-10

Variable

STEP

THETA

M

Description

This flag is zero for planar flow. It is one

for axisymmetric flow if Eq. (3) is used, and 2

if Eq. (5) is used. Eq. (3) should be used only

for straight channels. Eq. (4) can be used for

more complex shapes, but has not been tested for

Re>100000.

Description

Time step for the numerical algorithm. Input in

the form X.XXXE+XX, right justified.

Numerical scheme selector. For the

three-point-backward scheme, Z = 0.5.

Numerical scheme selector. For the

three-point-backward scheme, THETA = 1.0.

Maximum number of grid points in the i-direction

(^-direction). M should be less than 80.

77
41-50

Card 3

Column

1-10

11-20

21-30

31-40

41-50

51-60

61-70

Card 4

Column

1-10

11-20

21-30

Variable

UMAX

IFREQ

CRIT1

CRIT2

CRIT3

CRIT4

CRIT5

Variable

PO

TO

AMACH

31-40 RAD

Maximum number of grid points in the j-direction

(rrdirect!on). L should be less than 45.

Description

Maximum number of iterations (time steps)

allowed.

Frequency with which intermediary data are

output. If no intermediary output is desired,

make IFREQ greater than ITMAX. If output of

the dependent variables in the entire flow field

is desired every 10 iterations, IFREQ = 10.

Value below which the density residual must fall

for convergence.

Value below which the residual of pu must fall

for convergence.

Value below which the residual of pv must fall

for convergence.

Value below which the residual of e must fall

for convergence.

Value below which the wall pressure residual

must fall before a wall correction is made in

design calculations. If no inverse design is

desired, make CRIT5 smaller than CRIT1, CRIT2,

CRIT3, CRIT4. CRIT1 to CRIT5 are input in the

form X.XXE-XX and right justified.

Description

Inlet stagnation pressure (Psia)

Inlet stagnation temperature (°R).

Inlet Mach number. For IFLAT=0, (see below),

AMACH and the stagnation quantities are used to

calculate the inlet velocity UINF, which is kept

constant and is used as a reference velocity.

The same thing applies when IFLAT=1 and IPOUT=0.

For IFLAT=1 and IPOUT=1, this Mach number is

discarded and recalculated from PO and POUT.

Characteristic length (i.e. chord length for

airfoil, inlet radius for nozzle, diffuser,
11-50

51-60

61-70

Card 5

Column

1-10

11-20

21-30

31-40

11-50

51-60

61-70

Card 6

Column

1-10

ALPHA

IV IS

IFLOW

Variable

IREAD

ITURB

IFLAT

INP

IPRT

ISYM

ICASC

Variable

IB

etc...) in feet. RAD is input in the form

X.XXXE^XX.

Artificial viscosity coefficient for the

explicit part (c-|).

IVIS=1 for viscous flow. IVIS=0 indicates an

inviscid flow. However, adjustments have not

been made for inviscid calculations.

Use zero if only grid computations are needed,

and one if both grid and flow computations are

needed.

Description

Use zero if initial flow data are internally

prescribed (usually uniform conditions) and one

if initial flow data are read from a file.

Use zero for laminar calculations and one for

turbulent calculations.

Use one for airfoil calculations, and zero for

flow between parallel plates, symmetric nozzles

or diffusers, etc.. (where the centerline is a

boundary).

For IFLAT=0, use INP=0 if inlet conditions are

internally prescribed (usually uniform

conditions) and INP=1 if inlet velocity

distribution is manually input. Not used if

IFLAT = 1.

Use IPRT=1 if initial flow data throughout the

flow field are to be output, and IPRT=0

otherwise.

Use ISYM=1 for symmetric airfoils at zero

incidence angle and ISYM=0 otherwise. Not used

for IFLAT=0.

Use one for cascades, and zero for isolated

airfoils. Not used for IFLAT=0.

Description

Number of points along each periodic boundary in

cascades. For isolated airfoils, IB=1.

79
11-20

21-30

31-^0

Card 7

Column

1-10

11-20

21-30

31-^0

141-50

Card 8

Column

1-10

11-20

21-30

1C

IG

IGP

Variable

ITEMP

IPOTO

THIN

THOUT

XTE

Variable

ITWALL

TWALL

POUT

IB + number of points along the inlet boundary

+1. For isolated airfoils, 1C = M.

Number of points along the wake boundary

(between and including exit and trailing edge).

IG + Number of points along the airfoil surface

(counting the trailing edge only once). IB, 1C,

IG and IGP are disregarded for IFLAT = 0.

Description

Use zero if temperature is not maintained fixed

at the inlet, and one if it is (in which case

velocity is also kept fixed at the inlet and

pressure is extrapolated). Not used for

IFLAT=0.

Use zero if stagnation pressure and stagnation

temperature are not kept fixed at the inlet, and

one if they are (in which case flow angle is

prescribed and Riemann invariant is

extrapolated). To maintain p, u and v fixed at

the inlet, use ITEMP=IPOTO=0. For IFLAT=0, use

ITEMP=IPOTO=0.

Prescribed angle between the inlet flow

direction and the x-axis (degrees). Not used

for IFLAT=0.

Rough estimate of the exit flow angle (degrees).

Used in internal calculation of initial velocity

field. Not used for IFLAT=0.

The x-coordinate of the trailing edge,

non-dimensionalized with chord length. Not used

for IFLAT=0.

Description

Use zero for an adiabatic wall and one if wall

temperature is prescribed. Not used for

IFLAT=0.

Prescribed wall temperature for ITWALL=1 (°R).

Prescribed exit pressure (Psia). Not used if

IFLAT=1 and IPOUT=0.

80
IPOUT

Card 3

Column

1-10

Variable

INVERS

Use IPOUT=1 if the exit pressure is input as

POUT. Use IPOUT=0 if the exit pressure must be

internally specified as the free stream

pressure.

Description

Use one if inverse calculation is desired, and

zero otherwise.

11-20 NEWAL For IFLAT=0, use NEWAL=0 for secant method and

NEWAL=1 for virtual wall velocity method (for

IFLAT=1 only virtual wall velocity method is

used).

21-30 IDES1 Value of i corresponding to the beginning of the

controlled region in inverse design. For

IFLAT=0 use IDES1=1.

31-40 IDES2

41-50

51-60

FACT

RLOW

Value of i corresponding to the end of the

controlled region in inverse design. For

IFLAT=0, use IDES2=M. For airfoil calculations,

the controlled region must be on the upper

surface, unless modifications to the code are

made. For IFLAT=0 (and for the axisymmetric

case), only the upper boundary (J=L) can be

modified.

Multiplication factor used in smoothing for

virtual wall velocity method (0.2 has been

used).

Minimum allowable value of wall distance from

centerline in wall correction calculations with

secant method. Used in smoothing, to discard

erratic values.

61-70 RHIGH Maximum allowable value of wall distance from

centerline (see RLOW).

Card 10

Column

1-10

11-20

Variable

IM

IM

Description

Same as M for grid calculations where ICLUS1=0

or when ISOREN=1 (see below).

Same as L for grid calculations where ICALC=0 or

ICLUST=0 or when ISOREN=1 (see below).

21-30 ITMAX Maximum number of iterations for grid

calculations.

81
31-40 ICALC

41-50

Card 11

Column

1-10

W

Variable

ISOREN

Card 12 (skip if ISOREN=1)

Column Variable

1-10 ICLUST

11-20

21-30

JMC

JMLL

If ICALC=0, the grid is obtained by

interpolation from specified boundary values.

If ICALC=1, a PDE solution (APPENDIX A) gives

the grid.

Relaxation coefficient in SOR solution of grid

PDE's (0<W<2.0).

Description

If ISOREN=1, the grid is read from a file (e.g.

a file obtained from Sorenson's program). If

ISOREN=0, the grid is internally calculated.

Description

For simple channels (diffusers, nozzles, etc...)

the grid can be obtained as follows: calculate

a coarse but smooth grid (ICALC=1) and then

refine the grid by adding grid lines, say near

the wall. If ICLUST=0, no grid refinement is

needed. If ICLUST=1, the coarse grid is

calculated using the input values of X(1,J) and

R(1,J) (inlet coordinates). Then for the finer

grid, a new distribution of the X(1,J)'s and

R(1,J)'s must be given (the XG(J) and RG(J)

values). With the new inlet distribution, the

maximum number of points in the J-direction

changes from JM to JMC. The new distribution of

inlet points does not have to cover the entire

inlet. For example, one may need to keep points

from J=1 to J=JMLL as they are, and input a ne.w

distribution for J=JMLL+1 to J=JMC (see Figure

E-1). The XG's and RG's are used to read in

these JMC-JMLL values, and the grid in the rest

of the domain is obtained by spline

interpolation from the old grid.

Maximum number of points in J-direction for

refined grid when ICLUST=1 is used. For

ICLUST=1, JMC is equal to L.

Number of inlet grid points (starting with J=1)

which are left unchanged when a refined inlet

distribution of points (ICLUST=1) is prescribed.

If the refined distribution covers the whole

inlet, make JMLL=1.

82
• Coarse distribution (X( I,J), R( I,J))

— Finer distribution

J=JM

or J=JMCffr

J=JMLL+I --•>

J = JMLL"

J= I

XG(J )

RG(J ) 1

J=JMLL+ I to J = JMC

Figure E-1 : Grid refinement option (ICLUST=1)

Card 13 (Skip if ISOREN=1)

Column Variable

1-10 ICLUS1

11-20

21-30

Card 1U

Column

1-10

11-20

21-30

31-^0

IMC

IMLL

Description

Similar to ICLUST, but applies to I-direction.

Associated with IMC, IMLL, XGI(I), RGI(I). The

ICLUS1 feature is not used for inverse design

computations where grid computations are

performed many times during the same run, after

each wall correction.

Similar to JMC, but applies to I-direction. For

ICLUS1=1, IMC is equal to M.

Similar to JMLL, but applies to I-direction.

(Skip if ISOREN=1)

Variable Description

P1

PX1

P2T

PX2

Initial value of f-| (Eq. AU). A value of zero

is usually used.

Coefficient a in Eq. (AU).

Initial value of

is usually used.

Coefficient c in Eq. (AH).

A value of

83
Card 15 (Skip if ISOREN=1)

Column Variable Description

1-10 Q1 Initial value of g-( in Eq. (A4). A value of

zero is usually used.

11-20 QX1 Coefficient b in Eq. (AH) .

21-30 Q2 Initial value of g2

 in

 Eq. (AM). A value of

zero is usually used.

31-40 QX2 Coefficient d in Eq. (AH) .

Card 16 (Skip if ISOREN=1)

Column Variable Description

1-10 WP Relaxation parameter in solution of grid

equations (0<WP_<2.0).

11-20 WQ Similar to WP

21-30 PLIM Limitation factor to help stability of grid

solution (PLIM=1.0 has been used).

31-40 QLIM Similar to PLIM.

Card 16 (Skip if ISOREN=1)

Column Variable Description

1-10 WR Similar to WP

11-20 WS Similar to WP.

21-30 RLIM Similar to PLIM

31-40 SLIM Similar to PLIM.

Next input

Values of X(I,1), the x values of grid points along the boundary J=1. IM

values should be input, 5 per card, in 13 columns fields. Data are input in

the form X.XXXXXXE+XX. All values are non-dimensionalized by characteristic

length and are right justified. Skip if ISOREN=1.
Values of R(I,1), the r values of grid points along the boundary J=1.

IM values must be input. Skip if ISOREN=1.

Value of X(I,JM), the x values of grid points along the boundary J=JM.

IM values must be input. Skip if ISOREN=1.

Values of R(I,JM), the r values of grid points along the boundary J=JM.

IM values must be input. Skip if ISOREN=1.

Values of X(1,J), the x values of grid points along the 1=1 boundary.

JM-2 values are needed (J=2 to J=JM-1). Skip if ISOREN=1.

Values of R(1,J) the r values of grid points along' the 1=1 boundary.

JM-2 values are needed (J=2 to J=JM-1).

Values of X(IM,J), the x values of grid points along the I=IM boundary.

JM-2 values are needed (J=2 to J=JM-1). Skip if ISOREN=1.

Values of R(IM,J), the r values of grid points along the I=IM boundary.

JM-2 values are needed (J=2 to J=JM-1). These values are used for the

initialization but are not kept during the calculations. Skip if

ISOREN=1.

Values of XG(J), the x values of nodal points along the inlet boundary

for refined grid. JMC-JMLL values are needed. Skip if ICLUST=0 or if

ISOREN=1, or ICALC=0.

Values of RG(J), the r values of nodal points along the inlet boundary

for refined grid. JMC-JMLL values are needed. Skip if ICLUST=0 or if

ISOREN=1, or ICALC=0.

Values of XGI(I), the x values of nodal points along the J=1 boundary

for refined grid. IMC-IMLL values are needed. Skip if ICLUS1=0 or

ISOREN=1 or ICALC=0.

85
Values of RGI(I), the r values of nodal points along the J=1 boundary

for refined grid. IMC-IMLL values are needed. Skip if ICLUS1 = 0 or

ISOREN=1 or ICALC=0.

Values of PDES(I), the target pressure distribution, in Psia, for

I=IDES1 to I=IDES2 (from 1=1 to I=M for axisymmetric case). Skip if no

design calculations are desired.

Values of U2(1,J), the inlet velocity distribution, for J=1 to J=L, in

ft/sec. Skip if INP=0.

Values of the r coordinate at the wall (R(I,JM)) for second guessed

geometry if secant method is used. Data are input for 1=1 to I=M.

Skip if secant method is not used.

In the above, all boundary data are input 5 per card, in the form

x.xxxxxE+xx, and are right justified. Geometry data are

non-dimensional.

For the axisymmetric case, data are input in the same way, except

for the following cards:

Card 3

Column Variable Description

1-10 PINF Exit static pressure, used as initial pressure

everywhere (Psia)

11-20 TINF Initial temperature, used as reference

temperature (in °R).

21-30 UINF Inlet velocity, used as reference velocity (in

ft/sec).

31-^0 RAD See planar case.

11-50 ALPHA See planar case.

51-60 IVIS See planar case

61-70 IFLO Same as IFLOW in planar case.

86
Card

Column

1-10

11-20

21-30

31-40

Card 5

Column

1-10

11-20

21-30

Card 6

Variable

ITURB

INP

IREAD

IPRT

Variable

1TEMP

TO

ICENT

Variable

Description

See planar case.

See planar case.

See planar case

See planar case.

Description

Use ITEMP=0 if density is kept fixed at the

inlet, and ITEMP=1 if temperature is kept

constant. In both cases the inlet velocity is

kept fixed and pressure is obtained by

extrapolation.

This stagnation temperature and UINF are used to

obtain the static temperature which is

maintained fixed at the inlet.

Use ICENT=1 if the centerline is a boundary of

the flow field (J=1) and ICENT=0 if there is a

solid wall at J=1 (annular type). For ICENT=1,

only laminar calculations can be performed,

unless subroutine CEBECI is modified to account

for two walls.

Description

See planar case

See planar case

See planar case

See planar case

See planar case

INVERS

NEWAL

FACT

RLOW

RHIGH

For the axisymmetric case, skip cards 7 to 9 of the planar case.

Output Description

The output of the code first shows all the input parameters in the

order in which they are input. Reference quantities such as velocity,

87
density, temperature, viscosity, etc which are calculated by the code,

are also printed.

Once the input parameters have been printed, the coordinates of the

nodal points in the finite difference network are output: the X(I,J)'s are

printed first, as M sets of data, each set corresponding to values of J from

J=1 to J=JM. Then the R(I,J)'s are printed in the same manner. If ICALC=0,

these values are the ones obtained by interpolation. If ICALC=1, these are

the initial data for the grid computations. For ICALC=1 , the residuals DELX

and DELR for X and R are then printed for the first and the last iterations

of the grid computations. Then the final distributions of the X(I,J)'s and

the R(I,J)'s are printed.

If wall design calculations are to be done, the non-dimensional design

pressure distribution (along the controlled portion of the solid wall to be

modified) is printed next. If no design calculations are to be performed,

the above is skipped and the initial velocity, density, pressure and energy

fields are printed next if IPRT=1. The form of the output is the same as the

form of the X(I,J) and R(I,J) output. The printed values are

non-dimensional, and the title before each set of data indicates how the

non-dimensional values can be used to obtain the dimensional values.

The residuals DEL1P, DEL2P, DEL3P, DELUP and DEL5P are printed next.

Residuals DEL1P to DEL4P are defined as the maximum differences between the

•

values of p, pu, pv or e, respectively, between the current and the previous

iteration in time (ITIME). DEL5P is defined the same way, but only applies

to the wall pressure. Also printed are the locations (values of I and J

such as IHR and JHR for p, IHM and JHM for pu, IHN and JHN for pv, IHE and

JHE for e) where the maximum difference occurs. The residuals are printed

for each iteration, and their overall trend should be decreasing.
If wall corrections are performed, the convergence history given by

the above residuals will be interrupted whenever DEL5P becomes smaller than

CRIT5. At this point, the wall pressure in the controlled region

(corresponding to the current geometry) is printed. A new wall geometry is

obtained, and the corresponding r-coordinates are printed. Also printed is

DEL6P, the maximum difference between the current pressure and the target

pressure. DEL6P is printed after every wall correction. However, it

usually does not give an adequate indication on whether or not the current

pressure is approaching the target pressure, since the wall correction scheme

may not perform adequately at some points (such as the end points of the

controlled region). A plot of the pressure distribution in the controlled

region after each wall correction would give a better indication on the

success of the wall correction scheme. For airfoil calculations, the

calculations are stopped after each wall correction because a new grid must

be externally obtained.

When DEL1P, DEL2P, DEL3P and DEL1P become smaller than CRIT1, CRIT2,

CRIT3 and CRITi), respectively (or after the number of iterations ITIME

becomes larger than ITMAX), the calculations are stopped, and the final

density, velocity, pressure and energy fields are output.

In addition to the output described above some data are automatically

stored by the code on files. At the end of the grid calculations, the grid

is stored in a data set with the reference number 13. At the end of the flow

calculations, the entire flow field values of u, v, p, e, pu, pv, P are

stored in a data set with the reference number 12. This way, the grid and

flow data can be read from data sets 13 and 12 and used, for example, in a

plotting routine written by the user.

89
Additional Considerations

The present code has been run in FORTRAN IV on an IBM 3081 computer.

Double precision is used in the calculations. Adequate job control language

statements must be provided in order to run the code.

As mentioned in the input data description section, a grid file is

read by the code if ISOREN=1 and a flow field data file is read if IREAD=1.

These files correspond to data set reference numbers 8 and 11, respectively.

Adequate job control language statements must be incorporated to account for

these input data sets. Also, in cases where the data sets are not available

(which may be the case when IREAD=0 or ISOREN=0), if may still be necessary

to create them (even if they contain no useful information) in order to make

the execution possible.

As mentioned in the output description section, some data are stored

in data sets 12 and 13. Adequate job control language statements must be

incorporated to account for these output data sets. One interesting way of

using data sets 11 and 12 would be, for example, to start the computations

with IREAD=0. The final values of u, v, p, e, pu, pv for the calculation are

automatically stored in a data set named DATA1 for example, corresponding to

reference number 12. If more iterations are needed for convergence, one may

restart the calculations using IREAD=1 and, through the appropriate job

control language statements, associate data set DATA1 to reference number 11

and use a new data set name, i.e. DATA2, in association with reference number

12. Then for the new run, the initial flow field is read from 11 and the new

results are stored in 12 (DATA2). This allows for a more efficient use of

the code in very long calculations.

In order to read flow data stored by the code in a data set, a

statement such as the following must be used:

90
DO SN1 J=1,L

DO SN1 1=1,M

SN1 READ(11,SN2), U2(I,J), V2(I,J), RH02(I,J), E2(I,J), M2(I,J),

1 N2(I,J), P2(I,J)

SN2 FORMAT(5X, 7E15.6)

In the above, SN1 and SN2 are statement numbers. The same format must be

used in writing data into a file to be used by the code if IREAD=1. In order

to read grid data stored by the code in a data set, a statement such as the

following must be used:

DO SN1 J=1,L

DO SN1 1=1,M

SN1 READ(8,SN2), X(I,J), R(I,J)

SN2 FORMAT(2E14.6)

Sample Input and Output

The following is a simple example of a set of input data, and a

computer plot of the velocity field (see Figure E-2). Notice that the plot

was obtained by a small, separate code. Since plotting software varies with

facilities, no plotting routine was incorporated in the code.

91
0.500E+00 0.5 1.0 17

450 1001 0.10E-OU 0.10E-04

14.70 5UO.O 500.0 0.335E-02

000 1

1 542.0 1

0 0 0.0 0.0

17 11 100 1

0

1 19 1

0 0 1

0.0 0.3 0.0 0.7

0.0 0.3 0.0 0.7

0.3 0.3 1.0 1.0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0.3 0.

.OOOOOE+00

.62500E+00

.12500E+01

. 18750E+01

.OOOOOE+00

.OOOOOE+00

.OOOOOE+00

.OOOOOE+00

.OOOOOE+00

.62500E+00

.12500E+01

.18750E+01

.10000E+01

.10342E+01

.10999E+01

.11472E+01

.OOOOOE+00

.OOOOOE+00

.10000E+00

.60000E+00

.20000E+01

.20000E+01

.10000E+00

.60000E+00

.OOOOOE+00

.OOOOOE+00

.OOOOOE+00

.OOOOOE+00

.15473E+00

.72896E+00

. 92463 E+00

.98875E+00

3

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

1.0

. 12500E+00

.75000E+00

.13750E+01

.20000E+01

.OOOOOE+00

.OOOOOE+00

.OOOOOE+00

.OOOOOE+00

.12500E+00

.75000E+00

.13750E+01

.20000E+01

.10000E+01

.10473E+01

.11130E+01

.11472E+01

.OOOOOE+00

.OOOOOE+00

.20000E+00

.70000E+00

.20000E+01

.20000E+01

.20000E+00

.70000E+00

.OOQOOE+00

.OOOOOE+00

.OOOOOE+00

.OOOOOE+00

.30946E+00

.78717E+00

.94371 E+00

.99500E+00

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

1.0

.25000E+00

.87500E+00

.15000E+01

.OOOOOE+00

.OOOOOE+00

.OOOOOE+00

.25000E+00

.87500E+00

.15000E+01

.10000E+01

.10604E+01

.11261E+01

.OOOOOE+00

.OOOOOE+00

.30000E+00

.80000E+00

.20000E+01

.20000E+01

.30000E+00

.80000E+00

.OOOOOE+00

.OOOOOE+00

.OOOOOE+00

.OOOOOE+00

.U5157E+00

.83374E+00

.95896E+00

.10000E+01

0

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

19

. 10E-04 0

2.0

0.0

1.55

37500E+00

10000E+01

16250E+01

OOOOOE+00

OOOOOE+00

OOOOOE+00

37500E+00

10000E+01

16250E+01

10079E+01

10736E+01

11393E+01

OOOOOE+00

OOOOOE+00

40000E+00

88000E+00

20000E+01

20000E+01

40000E+00

88000E+00

OOOOOE+00

OOOOOE+00

OOOOOE+00

56525E+00

87099E+00

97117E+00

. 10E-04 0.10E-05

1 1

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

0.

50000E+00

11250E+01

17500E+01

OOOOOE+00

OOOOOE+00

OOOOOE+00

50000E+00

11250E+01

17500E+01

10210E+01

10867E+01

11472E+01

OOOOOE+00

OOOOOE+00

50000E+00

95000E+00

20000E+01

20000E+01

50000E+00

95000E+00

OOOOOE+00

OOOOOE+00

OOOOOE+00

65620E+00

90079E+00

98094E+00

92
o

o

0)

Q>

r-(

0)

o

CJ

k 1 1 «M •

O O

0

0) O

(— 1 O

\ T t ' 1

I 1 I

1 j j J i (

tH

•0 fl)

<u •*->

•O "S

3 e

r— i *M

c; x

c o

if 1

. -Hi.

1

i A i ;t / t / t y\

 a

' c a

nJ (u

£< 05

jj «r-4

A A Ji

A A A i\ A A 1 L. P

nt

V

c

(0 i.

,

3

 £

i

*M **

<U !

•o c

Ittt

A A A i / I A A L, 03

53 -o

C rM

-4 O

B

3

13

>>

r-» 0)

(MA /

1

JftA

1

11

Wih

1

A A A /

A A A /»

A A A A

i /

A

A

i 1

t

n

i ;

A

t

\ I

A

<£

1

•^

x •

(0

£.

o

t S0)

c

5*

"»~l

o

o

r^l

0)

>

• •

CJ

t

M

3

bO

o

n

v

«->

^.

4H

o

rt

s

<u

c

3

0}

""'

93
1 Report No

NASA CR-175037

2 Government Accession No 3 Recipient's Catalog No

4 Title and Subtitle

Viscous Compressible Flow Direct and Inverse

Computation with Illustrations

5 Report Date

January 1986

6 Performing Organization Code

7 Author(s)

Tan-Teh Yang and Francois Ntone

8 Performing Organization Report No

None

10 Work Unit No

9 Performing Organization Name and Address

Clemson University

Mechanical Engineering Department

Clemson, South Carolina

11 Contract or Grant No

NAG 3-538

12 Sponsoring Agency Name and Address

National Aeronautics and Space Administration

Washington, D.C. 20546

13 Type of Report and Period Covered

Contractor Report

14 Sponsoring Agency Code

505-62-21

15 Supplementary Notes

Final report. Project Manager, R.J. Stelnke, Internal Fluid Mechanics Division,

NASA Lewis Research Center, Cleveland, Ohio 44135.

16 Abstract

An algorithm for laminar and turbulent viscous compressible two-dimensional flows

1s presented. For the application of precise boundary conditions over an arbitrary body surface, a body-fitted coordinate system 1s used 1n the physical

plane. A thin-layer approximation of the Navler-Stokes equations 1s Introduced

to keep the viscous terms relatively simple. The flow field computation 1s performed 1n the transformed plane. A factorlzed, Implicit scheme 1s used to facilitate the computation. Sample calculations, for Couette flow, developing pipe

flow, an Isolated airflow, two-dimensional compressor cascade flow, and segmental

compressor blade design are presented. To a certain extent, the effective use

of the direct solver depends on the user's skill 1n setting up the grldwork, the

time step size and the choice of the artificial viscosity. The design feature

of the algorithm, an Iterative scheme to correct geometry for a specified surface

pressure distribution, works well for subsonic flows. A more elaborate correction scheme 1s required 1n treating transonic flows where local shock waves may

be Involved.

17 Key Words (Suggested by Authorfs))

Viscous compressible flows;

Computational fluid dynamics; Inverse

solutions; Compressor cascade

19 Security Classlf (of this report)

Unclassified

18 Distribution Statement

Unclassified - unlimited

STAR Category 07

20 Security Classlf (of this page)

Unclassified

21 No of pages

101

22 Price'

A06

For sale by the National Technical Information Service, Springfield, Virginia 22161
ivlational Aeronautics and

Space Administration

Lewis Research Center

Cleveland Ohio 44135

Official Business

Penalty for Private Use $300

SECOND CLASS MAIL

ADDRESS CORRECTION REQUESTED

Postage and Fees Paid

National Aeronautics and

Space Administration

NASA-451

NASA