General Disclaimer

One or more of the Following Statements may affect this Document

This document has been reproduced from the best copy furnished by the 

organizational source. It is being released in the interest of making available as 

much information as possible.

This document may contain data, which exceeds the sheet parameters. It was 

furnished in this condition by the organizational source and is the best copy 

available.

This document may contain tone-on-tone or color graphs, charts and/or pictures, 

which have been reproduced in black and white.

This document is paginated as submitted by the original source.

Portions of this document are not fully legible due to the historical nature of some 

of the material. However, it is the best reproduction available from the original 

submission.

Produced by the NASA Center for Aerospace Information (CASI)

https://ntrs.nasa.gov/search.jsp?R=19840018935 2018-11-29T21:18:16+00:00Z
f

.z

r y'

A

F

a pl

"" k

NASA CONTRACTOR REPORT 166565 - P1-

IV I,lo

Solution of 'Three-Dimensional Time-Dependent Viscous Flows

Part 3: Application to Turbulent and Unsteady Flows

(14ASA-CR- 166565-.Pt-3) SOLUTION OF

3-DIMENSIONAL TINE-DEPENDENT VISCOUS FLOWS.

PART 3: APPLICATION TO TURBULENT AND

UNSTEADY 21OWS (Scientific Research

A sso ciltes, Inc.) 93 p HC A 05/MF A01

N84-27003

Unclas

G3/34 13670

Bernard C. Weinberg

Henry McDonald

CONTRACT NAS2- 10016

September 1982

4

__ - _

_

- . 1.1

 ....1.
NASA CONTRACTOR REPORT 166565

Solution of Three-Dimensional Time--Dependent Viscous Flows

Part 3: Application to Turbulent and Unsteady Flows

Bernard C. Weinberg

Henry McDonald

Scientific Research Associates, Inc.

Glastonbury, Connecticut

Prepared for

Ames Research Center

under Contract NAS2-1001{,

NASA

National Aeronautics and

Space Administration

Ames Research Center

Moffett Field, California 94,035

i
r

^...^^,.......^,_. _ s ".C. s . e4k.., r. ..ryi.;.F. iri' OpL- .s^e► ^'. ^ .. , ..: -^' ,y. ...

TABLE OF CONTENTS

Pape

SUMMARY. . . . . . . . . . . . . . . . . . . . . . . . . . 1

INTRODUCTION. . . . . . . . . . . . . . . . . . . . . . . . . . 2

LIST OF SYMBOLS . . . . . . . . . . . . . . . . • . • . • . . • . . . 4

ANALYSIS. . . . . . . . . . . . . . . . . . . . . . . . 7

Background . . . . . . . . . . . . . . . . . . . . . 7

Coordinate System . . . . . . . . . . . . . . . . . . . . . 11

Governing Equations . . . . . . . . . . . . . . . . . . . . . . 13

Continuity Equation . . . . . . . . . . . . . . . . . . . . 14

Momentum Equations . . . . . . . . . . . . . . . . . . . . 15

Diffusion Terms . . . . . . . . . . . . . . . . . . . 17

Convective Terms . . . . . . . . . . . . . . . . . . . . . . 18

Energy Equation . . . . . . . • . . . • . • • • • • . . 19

Energy of State . . . . . . . . . . . . . . . . A . . . . . 20

Linearizations . . . . . . . . . . . . . 20

The Turbulence Model . . . . . . . . . . . . . . . . . 21

Spatial Difference Approximations . . . . . . . . . . . 23

OR Operator Notation . . . . . . . . . . . . . . . . 23

Application of Coupled Nonlinear Parabolic Equations . . 27 

Linearized Block Implicit Scheme . . . . . . . . . . . . . . . 28

Implementation of the LBI Scheme Employing the OR

Operator Technique . . . . . . . . . . . . . . . . . . . . . 34

Boundary Conditions and Initial Conditions . . . . . . . . . . 36

The Computer Code . . . . . . . . . . . . . . . . . . . . . . 37

Numerical Results . . . . . . . . . . . . . . . . . . . . 40

Two-Dimensional Steady Turbulent Flow , , , , , , , , , . . . 40

Weighardt Flat Plate . 4G

Newman Airfoil Case . . . . . . . . . . . . . . . . . . 44

Two-Dimensional Unsteady Flow . . . . . . . . . . . . . . . 46

CONCLUSIONS. . . . . . . . . . . . . . . . . . . . . . . . . . 51

ACKNOWLEDGEMENT . . . . . . . . . . . . . . . . . . . . . . . . . . 52
TABLE OF CONTENTS (Cont.)

Page

APPENDICES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

Appendix A _. Linearization Technique . . . . . . . . . . . . . n 53

Appendix B - Geometric Properties of Airfoils . . . . . . . . . . 58

Appendix C - Gengralized Operator Compact Implicit Schemes . . . 70

REFERENCES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

TABLES. . . . . . . . . . . . . . . . . . . . . . . . . 78

FIGURES . . , . . . . . . . . . . . . . . . . . 83
I

1

SUMMARY

The flow over a helicopter rotor is an important example of three-dimensional

time-dependent viscous flow. The boundary layers that develop on the rotor blades

play a significant role in that they set loss levels and control retreating blade

stall. Consequently, there is considerable interest in developing a numerical scheme

for solving the time-dependent, three-dimensional compressible viscous flow equations

in order to predict such Flow fields, and which can be used as an aid in the design of

helicopter rotors.

In order to further investigate the numerical procedure, we have exercised the

computer code that was developed under a previous phase o the current research

program to solve an approximate form of the three-dimensional unsteady Navier-Stokes

equations empluying a Linearized Block Implicit technique in conjunction with a QR

operator scheme. Results of calculations are presented for several two-dimensional

boundary layer flows including steady turbulerh and unsteady laminar cases.

A comparison of fourth order and second order solutions indicate that increased

accuracy can be obtained without any significant increase in cost (run time). The

results of the computations also indicate that the computer code can be applied to

more complex flows such as those encountered on rotating airfoils. Finally, th

geometry of a symmetric NACA four digit airfoil is considered and the appropriate

geometrical properties are computed.

h

M
INTRODUCTION

The behavior of boundary layers on wings and bodies has long been of interest

to aerodynamicists. In both steady and unsteady flows, the boundary layers are known

to govern a major portion of the lasses and to significantly influence the vehicle

lift and moment coefficients. When the flow is steady, boundary layer prediction

schemes based on numerical solution to the governing partial differential equations

of motion have reached a high level of sophistication and predictive accuracy, even

in three space dimensions. In unsteady flows, such as are commonly encountered in

rotary winged aircraft, ;ome progress has been made in two space dimensions but

little to date has appeared on unsteady three-dimensional boundary layers.

Two particular problems arise with time-dependent three-dimensional boundary

layers relative to the steady case. -The first of these is the rather obvious one of

time integration with its added requirements of transient accuracy coupled with an

increase in the computational labor. The second of these is the so-called negative

cross flow problem, which to some extent has troubled the steady boundary layer

prediction schemes. Kendall, et al (Ref. 1) discuss the negative cross flow prob]em

for steady three-dimensional boundary layers in a very illuminating fashion. This

particular problem arises when the spanwise component of velocity changes sign and

will be discussed in detail subsequently. Because of the interest by external aerbdynamicists in swept wing boundary layers where the negative cross flow problem

(in this case flow from tip to root) is not usually encountered, the negative cross

flow problem has not received a great deal of attention to date. However, in transient

flows, particularly those encountered on rotor blades in forward flight, negative cross

flows are frequently encountered. For instance, the advancing rotor blade has cross

flows of one sign during the first ninety degrees of rotation and these can change sign

over part of the blade during the second ninety degrees.

Thus to be of practical value, time-dependent three-dimensional boundary layer

prediction schemes require high computational efficiency and transient accuracy

coupled to the ability to treat arbitrary cross flow profiles.

In this report we describe the development of a computer code for the efficient

solution of three-dimensional time-dependent viscous flows on fixed and rotary

aircraft. The Linearized Block Implicit (LBI) technique of Briley and McDonald (Ref. 2)

in coordination with a tridiagonal QR operator scheme (Ref. 3) is employed to solve

the reduced turbulent Navier-Stokes equations which are derived for nonorthogonal
i

A

coordinates in generalized tensor form. The rationale for the choice of this

approach is discussed in detail in Ref. 3 and k.

The basic assumptions made in the derivation of these equations are that, the

pressure does not vary normal to the shear layer and that in the energy equation

the square of the normal velocity is neglected with respect to the other velocity

components (To - constant). The latter assumption is included only for computational

simplification purposes and is not essential in the analysis. For turbulent flows,

a two-layer mixing length model is employed and its formulation in generalized

tensor notation is given. A novel method is employed for solving the continuity

equation in conjunction with the reduced Navier-Stokes equations. The continuity

equation is split by employing the Douglas -Gunn procedure to obtain a consistent

approximation to the full equation which is then solved as an integral, Results of

several boundary layer calculoti. ons are presented and comparisons with experimental

data and other reported computations are made.

3

w

r

,.^^ti ^.,w...^ ,ems _^._.^._.______..____.__.^._..s;^•^
LIST OF SYMBOLS

a, c, d coefficients of differential equation 

^

aid coefficient of differential operator

B Btnornal vector

Cf skin friction coefficient 

nonlinear three:-dimensional differential operator 

ei covariant basis vector 

f^, f(x) right-hand side of differential equation

ij

gij , g components of metric tensor

H shape factor

h step size

hl , h2 scale factors

J JacobiAA 

K curvature 

Kg geodesic curvature

Kn normal curvature 

R mixing length

L(t&) linear differential operator

,441 linear three-dimensional differential operator

Li product operatori i

Mn difference between n and n

.x 4.

N normal vector

P pressure

Q tridiagonal difference operator

'.

_

q , q c , q+ components of Q

q velocity vector

R tridiagonal difference operator 

`

4
LIST OF SYMBOLS (Continued)

r', rc , r* components of R

f Re Reynolds number

Rci

cell Reynolds number

i S source term

r,

` t time

r

T static temperature

To stagnation temperature

'f

ui contravariant velocity component

U, v, w physical velocity component

xi coordinate direction in i direction

` GREEK SYMBOLS

h^

factor for centering time step

tL y ratio of specific heats

rid

Christoffel symbol

' d boundary layer thickness

d * displacement thickness

di Kronecker delta function

A velocity divergence

At time step

Ax step size in i direction

6 momentum thickness

ejk strain tensor

T1 boundary .layer coordinate

a viscous stability parameter

' u molecular viscosity

P density

5
LIST OF SYMBOLS (Continued)

CT 

ik stress tensor

dissipation function

vector of unknowns in LBI scheme

source term in LBI scheme.

T wall shear

w

W frequency of oscillation

SUBSCRIPTS

e edge

p physical

00 free stream

SUPERSCRIPTS

n time level

*) ** 1 *** intermediate time level in ADI procedure

6

--7

r
i

It

x

ANALYSIS

Background

Three-dimensional boundary Layers occur on the wings and fuselages of both

fixed and rotary wing aircraft. In both types of vehicles, the boundary layers are

important in setting loss levels and determining useful operating ranges. As is

well known, boundary layers are sensitive to pressure gradients. In tune-dependent

flow the temporal acceleration terms appear in the momentum equation in a form very

similar to the conventional imposed pressure gradient and so for qualitative evaluation purposes can be regarded as "pseudo" or "auxiliary" pressure gradients. Viewed

in this manner, the temporal acceleration terms can be seen to influence quantities

of practical importance such as skin friction, displacement thickness and the onset

of separation. At the range of frequencies typically encountered in rotary wing

aircraft aerodynamic problems, it is clear, for instance, from the extensive review

of McCroskey (lief. 5), that very significant transient boundary layer effects can

be observed.

In examining the flow problems of prantical interest such as loss levels or the

onset of separation it is evident that all three space dimensions must be considered.

In conventional aircraft, the sweep effect is of interest and inherently threedimensional. In rotary wing aircraft, in forward flight clearly very substantial

transient changes occur in what might be termed the local sweep angle. However,

generally speaking, the boundary layers remain thin unless catastrophic flow

separation occurs or the flow at the wing or rotor tip is considered.

Conventional boundary layer integration schemes have developed by forward

marching the streamwise velocity u in the streamwise x direction and simultaneously

marching out along the span in the z positive direction. In general, the spanwise

marching scheme does not normally encounter negative w, i.e., spanwise inflow. This

is very fortunate because it is difficult, indeed it could be argued impossible, to

structure a physically satisfactory unconditionally stable noniterative scheme which

permits forward marching in the spanwise direction with a negative w cross flow.

At least intuitively the problem of negative cross flow implies information being

transferred upstream against the spanwise marching direction. Conventional stability

analyses confirm the inability to forward march into regions of significant negative w.

From experience with attempts to march the two-dimensional boundary layer equation

into a region of separated flow and its obvious relationship to the negative cross

7
flow problem, it is not surprising that spanwise marching into a negative crous

flow regi n 16 not accomplished witli p ut special treatment. Recently, conventional

boundary layer devoJ,opers have been turning to performing an implicit spanwise

construction to remove the restriction of only positive cross flows (Kendall, et al,

Ref,l). Lin and Rubin (Ref. 6) in their predictor-corrector boundary region

solutions for flow over a yawed cone at moderate incidence have also shown that

retaining diffusion in the spanwise direction not only eliminates the problems

associated with negative cross flows, but improves upon the solutions obtained by

standard three-dimensional boundary layer techniques. Furthermole, boundary

conditions applied at the tip can influence the flow inboard, if required by the

physics of the flow.

For these reasons, the implicit spanwise construction has been a feature of the

three-dimensional duct flow analysis of Briley (Ref. 7) and McDonald and Briley

(Ref. 8). As a consequence of these observations and the need to remove the negative

cross flow restriction, a spanwise implicit formulation seems mandatory for

rotary wing applications and at least desirable for fix ,^d wing applications,

expecially as it can be had for a very modest increase in code computational labor.

Based on the experience in Reft. 7 and 8, the spanwise implicit sweep would only

result in a moderate increase in computational effort relative to the explicit

spanwise marching approach. The extension of the conventional three-dimesional

boundary layer equations to allow spanwise diffusion is easily accomplished, and in

view of the improved physical representation which thus follows, It is recommended

and has been implemented in this effort.

As a matter of course, it has been assumed that normal to the wall an implicit

formulation would be structured. In recent years for boundary layer type problems

there has been little dispute as to the efficiency gains to be had from an implicit

formulation normal to the wall (Ref. 9). However, in the streamwise direction for

steady 2-D flow, the equations are normally forward marched and the implicit

stability obtained entirely from being implicit in the normal to the wall direction.

In time-dependent flows a similar structure is to be had so that at each time level

one streamwise (explicit) forward marching sweep could be made with two implicit

sweeps in the spanwise and normal directions to give the desired unconditional

stability. As mentioned earlier, the explicit sweep would probably require less

computational effort by about 20% than an implicit streamwise sweep and of course

less storage. However, since the solution is being time marched the opportunity to

take a streamwise implicit sweep at roughly the same cost as the explicit sweep

8
t

sloes arise. If one does perform a streamwise oxplicit sweep, then the linearization

of nonlinear terms is performed about the known serial marching level. If an

implicit #treamwise structure is adopted, then full time linearization can be

utilized, That is the linearization of the nonlinear terms is performed about the

known time level.. As is pointed out in Ref. 8 1 it is easier to obtain a consistent

high order accurate spatial-temporal linearization by marching in time than in space

(in time the marching derivatives have the form u whereas in space they are nonlinear and have: the form u iui ). Further, by structuring implicitly in the space

marching direction, small regions of axial reverse flow would be permitted. As a

result of the combined benefits of consistent hkgh order linearizations and the

inclusion of small separated zones, a streamwise implicit structure is advocated

and has been employed in this effort,

Since a full 3-1) spatial ;integration is carried out at each time step of a

transient calculation, spatial accuracy plays a very important role in the overall,

efficiency of the numerical method. By choosing a spatial scheme of sufficient

accuracy, one can expect to reduce the total number of grid points to an acceptable

level. However, for the chosen scheme one must account not only for the spatial

discretization errors but those intraduc;ed us a result of the linearization procedure.

In order to get VI-A mos t out of a given spatial difference formula, the errors

incurred in reprr, „ '^-•,I1g nonlinear terms by linear combinations of terms should be

less than -;kr ecvn.1 to the spatial discretization errors. If the linearization introduces a greater error than the spatial differencing, then either a coarser spatial

mesh could be used, or iteration, or some form of linearization improvement is called

for. Iteration across a time step is not recommended since this only reduces the

linearization error and computationally costs as much as a complete time step.

Cutting back the time step would be prefereable to iterating to preserve the linearization error at some acceptable level, since cutting back on the time step would improve

both the transient error and the linearization error. This point is clearly demonstrated

in Ref. 3. To obtain a linearization, which introduces errors of at most the same

as the spatial difference formulae, a Taylor series expansion about the known time

level can be performed. This process clearly demands a formal block, i.e., coupled,

treatment of the system of equations. For instance, in the streamwise momentum

equation a typical term is linearized;

(uW)no . Un+1Wn + u

n wno _ un wn + O( Ate)

9
and clearly one cannot lag w

n*1 at the old time level n without introducing a firstorder time error in order to get an uncoupled system, i.e,, wn*1 not appearing in

the streamwise momentum equation, ':tius, formal linearization and consideration of the

resulting errors indicate the coupled system a- yht to be treated from the accuracy

point of view. This is further reinforced when it is realized that block, i.e.,

coupled, systems are not computationally expensive (in a relative sense).

Additionally, a second type of approximation arises unconnected with linearization but arising from basic coupling terms in the original equations and if indeed

some terra!, in an equation are time lagged in order to uncouple the equation system

and these terms are of equal importance to the terms retained, then again an iterative

updating is called for in order to achieve stability, accuracy and consistency.

This could be termed ad hoc equation uncoupling, Blottner (lief. 9) has shown that

many iterations around the ad hoc uncoupled set (> 10) are sometimes required in

order to achieve an overall solution accuracy commensurate with the local difareuce

molecule accuracy.

The linearization technique discussed above is described in Ref. 8, together

with its application to block coupled splitting schemes. Schemes of this general

type are here referred to as "split linearized block implicit" or split LBI schemes,

and are reviewed in detail by Briley and McDonald (Ref. 2).

As a general observation, care is required to obtain acceptable transient

accuracy for long time integration with conventional finite difference schemes.

A Crank-Nicolson centered time- implicit scheme for instance, although second order in

time, shows quite a disp ersion problem (relative to other schemes) on the simple

pure convection problem. However, the problem of transient accuracy is significantly

reduced in the typical boundary layer problem since the time dependency is continuously input through boundary conditions and the concern is with relatively

"short' time integrations. The computational problem is more of what the Phase lag

of the wall shear is, relative to the prescribed free stream disturbance, than

concern over the convection velocity oY a wave in a shear after a long propagation

time. The interest is in forcrd oscillations with a minimum scale of the boundary

layer thickness over a few cycles of the motion, just enough to obtain repetition

cyclically. It is, therefore, expected: that a significant dispersion problem will

not arise with a conventional implicit scheme.

The governing equations that are considered here are the Navier-Stokes ei;uations,

continuity, energy and the equation of state which are written in generalized tensor

form for a body oriented coordinate system (boundary layer coordinates). In

10
accordance with the boundary layer assumptions, the normal momentum equation is

eliminated and the pressure is specified throughout the viscous layer in its stead.

For the energy equation, constant stagnation temperature T o is assumed. This

assumption is a good approximation for the flow fields considered, and is thus

included here only for purpoies of simplification. In the an.,zlysis that follows,

the full energy equation could equally well have been used. Employing the equation

of state which relates the pressure p to the velocity components u and w by an

algebraic equation, the problem can be reduced to one involving only the three

velocity components, u, w and v and three equations, the streamwise and spanwise

momentum equations and the continuity equation. Hence, we consider a block-three

system rather than a block-four system which leads to a significant reduction in

computer time. If the full energy equation were to be considered, a block-four

system would. result due to the inclusion of the temperature as an additional unknown.

Coordinate System

Since the goal of this effort.is to solve for the flow over airfoils an understanding of the type of geometries to be considered is essential to guide the

choice of the coordinate system and the structure of the computer code. The coordinate system is not only dependent upon the geometry of the airfoil, but also upon the

approximations that are made to the governing Navier-Stokes equations. -As in boundary

layer theory, we also assume that in the approximate form of the Navier-Stokes

equations the pressure is constant normal. to the shear layer. Inherent in the

assumptions is that the shear layer is thin. As pointed out by Howarth (Ref. 10)

the boundary layer assumptions lead to the conditions that one coordinate direction

must be normal to the body surface (this being straight) while the other coordinate

directions must lie on the body surface. These conditions uncouple the metric data

on the surface from that in the normal direction. Hence, the metric data for the

surface coordinates are functions of the surface coordinates alone, while the metric

data for the normal. coordinate direction are functions of that coordinate alone.

The choice of the surface coordinates is rather arbitrary and is based on considerations such as the ease of construction or the grid distribution on the wing surface.

' In the numerical solution of the flow over an airfoil there are many advantages

to be gained by the judicious choice of coordinates. The most obvious advantage is

that the physical boundaries of a flow region can be represented by coordinate

surfaces. '"his removes the need for fractional cells in general; hence, the complications and loss of accuracy associated with a boundary interpolation are removed.

11
Another advantage is that a uniform numerical m,-t',od can be used. The solution can

then be performed with a fixed number of cells in any given direction and with a

uniform meFb spacing. For a rectangular plan form wing a. Cartesian coordinate

system would be adequate. However, for more general plan forms such art a ewept wing

a nonorthogonal grid which conforms with the boundaries is preferred s-Lnce it covers

the entire airfoil while a Cartesian grid would not.

Another consideration is the selection of a coordinate grid distribution;

the major objective being the resolution of large solution gradients. The approach

taken here is to construct coordinate transformations that contain distributions for

physical mesh points. In this context, tha uniform mesh of computational space is

simply mapped into a suitably distributed mesh in physical space. When the transformation contains the mesh point distribution, there is no need to construc the

apparatus for the discrete approximation of derivatives on a non,i iform mesh. This

results in a savings in both computer logic and storage.

Hence, in this work, a coordinate system is chosen that conforms with the

boundaries of the physical domain, i.e., the wing surface which in general will be

nonorthogonal. In addition, provisions are made for analytical grid transformations

(Ref. 11) in each coordinate direction, in order- to suitably distribute grid points

in regions of large gradients.

In view of the type of geometries to be considered and the assumptions made to

obtain the approximate form of the Navier-Stokes equations, a specialized nonorthogonal

coordinat- system is advocated where the metric tensor which has four independent

components is given by

911 912

0

9 1) gig 922 0

0 0 9 33

The subscripts 1 and 2 refer to the directions on the surface of the body while

subscript 3 refers to the direction normal to the body. Since the metric data in

the coordinate directions on the airfoil surface are not functions of the normal

direction, the metric data in a i - 2 surface above the body are evaluated on the

body surface (Ref. 10). Furthermore, due to the use of nonorthogonal coordinates

12

-• -+,F .. ..8, „ ^.:..; .mod: ®a9-,......^ F - ^:r - " - --'"iY:cLs
it becomes advantageous to derive the equations in general nonorthogonal coordinates

employing generalized tensors. Details of the generalized tensor notation can be

found in Refs. 12, 13, and 4.

An important feature of the analysis to follow is that the governing equations

which are derived, under the prescribed assumptions, are invariant for any coordinate

system or any grid transformation (although, of course, the physical approximations

are coordinate dependent). The grid transformations are absorbed into the geometrical

coefficients, leaving the equations unaltered in form. This feature allows for the

construction of the geometric data to be contained in one subroutine with the definition of the metric data and their derivatives as input.

Other forms of the governing equations for nonorthogonal coordinates which do

not rely on generalized tensor notation were considered, i.e., such as the steady

three-dimensional boundary layer equations as given by Cebeci (Ref.14 ). Since in

that set of equations the only viscous terms retained are these in the normal direction,

any nonorthogonal effects are introduced solely through the convective terms. This

is a direct consequence of the assumptions which decouple the normal coordinate

from the surface coordinates. In those equations, the geometric terms that explicitly

appear are the metrics, h l and hZ , the geodesic curvatures kgl and k9Z , the angle

between the surface coordinate lines and an additional nonorthogonal curvature term.

In contrast to Cebeci's equations, the set of equations considered here, the

approximate form of the Navier-Stokes equations, allow for all the viscous terms to

be retained. A description of these terms without generalized tensor notation would

be cumbersome. Hence, the generalized tensor approach was chosen.

In Appendix B the geometric properties of a surface in three dimensions are

discussed and where appropriate, the generalized tensor equivalents are given.

In addition, a symmetric (uncambered) NACA four-digit airfoil is considered and the

pertinent geometric coefficients are presented.

Governing Equations

In view of the ultimate goal of this program., to solve an approximate form

of the unsteady three-dimensional Navier-Stokes equations on airfoil shapes, the

governing equations are derived in general nonorthogonal coordinates and are

given in generalized tensor notation.

In the following derivation the governing equations are nondimensionalized as

follows, xi with respect to the characteristic length L, the velocity with respect to

Um , density, pressure and temperature with respect to p . , p .U. and U. /cp respectively

and time with respect to L/UCO . Viscosity is nondimensionalized with respect to u...

13
Continuity Equation

Consider the continuity equation wirtten in vector form so that it is independent of coordinate system, i.e.,

dp

at 

+ v pq = n (^>

where p is the density and q is the velocity vector, which can be written

in a covariant basis as

q - ui e i (2)

where u is the i-th contravariant velocity component and ei is the covariant basis

vector in the x i direction. The velocity vector could have been expressed in a number of

different forms, each with its own attribute. Here for the moment the velocity vector

is expressed in a covariant basis, for simplicity. Subsequently, the velocity vector

will be transformed into its physical components for the numerical solution of the

governing equations. The reason for this is that the contravariant basis exhibits

variation in its components if the coordinates are such that the metric varies.

For boundary layer flows, the physical. velocity components are roughly aligned with

the coordinates and exhibit no variation with the metric per €e. As such, it is

felt that the actual computations are better performed on the pnysical components.

The divergence of a vector is given by

P pq = puk l k = (puk k

 + pug rk = J ( Jpuk ), k

(3)

where puk I k is the covariant derivative, pu k , k is the partial derivative in the xk

direction, J is the Jacc'aian and rik is the Christoffel symbol.

In Equation 3 two forms of the divergence are presented, one involving the

Christoffel symbol or curvature term directly and the other the Jacobian. The

former is perhaps more restrictive since it requires additional smoothness of the

geometrical quantities. However, both forms are used depending on which is more

convenient for the given application. For the momentum equations the form involving

the Christoffel symbols is employed in the evaluation of the explicit (lagged) diffusion terms while for the continuity equation the form involving the Jacobian is

used. Thus, the continuity equation can be expressed as

dp + I (Jpuk ),k = O

(4)

at 

14

4
ORIGINAL PAGE ig

OF POOR QUALITY

Momentum Equations

The momentum equations in vector form can be written as i

J

^q PP aq ^q .p) q V , a

.fit at+

(5)

where a is the stress tensor.

,'

4

In generalized tensor notation, Equation 5 becomes

4

P^ + Ukui lk I e i =

a ikl kei dt i

(6)

The stress tensor is defined as

srY

aik _ -

elk

3 R 

Q1Sik+

lP+ i

(7)

6.;

where U is the viscosity, p is the pressure, A is the velocity divergence and eik

is the strain tcn$or, and the Reynolds number, Re, is defined as p,,U.L/}a.. The

strain tensor expressed in terms of velocity components is

.t

Elk = Ui Im mk + uklm Sim (8)

where gmk and gmi are the components of the metric tensor. Employing the fact that

dkj = gka and substituting the definition of the strain into the stress tensor,

we obtain

a,ik _ + 2 μ Eik+ μ u i I gmk+ '4 u k ( gmi

\p 3 Re J Re m Re m (9)

Substituting Equation 9 into the momentum equation and employing the relationship

A

k

j

g I k = 0 (10)

r

t

r

we obtain for the i-th momentum equation in the e direction

K P r du i + u k u il 1 gik(p + 2 μ 1^

at k ` 3 Re J k (ll)

mk

+ 9 [W ui lm] Ik+ 9m^^ R uk lm^ Ik

15
In Ref, 3 it was pointed out that the QR operator scheme requires that derivatives

in any direction operate on only one variable. In the momentum equation this requirement prevents the implicit treatment of certain diffusion terms that arise due to the

curvature of the body. Although these terms are often treated explicitly anyway the

use of standard finite difference techniques 'instead of the QR operators would give

one the opportunity to treat these terms implicitly, if so desired, This and the

use of the quasi-linear form of the governing equations are the only major limitations that arise in the QR operator treatment of the approximate form of the

Navier-Stokes equations. In the usual boundary layer approximations, these explicitly

treated terms would not appear in the equations since they are of order 0 (Re-1/2 ) or

smaller, and should, therefore, be of little consequence. In principle, the quasilinear and, for instance, the full conservative forr- of the differential equations,

are equivalent. In discrete form, various formulations of the governing equations

exhibit different properties (Ref. 15). In the present problem, no distinct disadvantage appears to arise from the required use of a quasi-linear form of the

governing equations.

The requirement that derivatives in any direction operate only on one variable

would be more restrictive in the treatment of the pressure gradient term in the full

Navier-Stokes equations. The linearization of this term introduces derivatives of

all the velocity components in a given direction. According to the limitations of

{ the QR operator scheme described above, some of these terms must be treated explicitly.

Since an explicit treatment of these terms could reduce the stability bound of the

calculation scheme, an alternate procedure should be considered. This would involve

the addition of an auxillary equation relating the pressure gradient term to the

i derivatives of the velocity components and would increase the block size of the

system. An assessment of the efficiency of such a procedure has not been carried

out and further work in this area would appear to be warranted.

In the discussion that follows, the governing equations are first split into

an explicit part and an implicit part in accordance with the QR operator requirements.

Thereafter, the resulting equations are cast into "standard form", so that the

equations can be appropriately linearized and treated with the LBI technique.

Since mixed partial derivatives are commonly treated explicitly in orthogonal

coordinate systems, a similar procedure will be used for generalized nonorthogonal

coordinates where the explicit treatment is extended to include mixed second covariant

derivatives. All other second covariant derivatives are retained as implicit.

Although such a procedure would automatically treat more terms explicitly than one does

for orthogonal coordinates, it simplifies the bookkeeping requirements in the

16
OF POOR QUALIV

construction of the computer code, and is thus adopted here. Conversely, retaining

the covariant derivative as a unit (whether implicit or explicit) could prevent

instabilities that may arise due to time splitting, This can occur when two

portions of one term should cancel identically, but do not due to their being

w

split between two sweeps.

Aifusion Terms

Consider the term

(μuml1 ) ! k (12)

If j = k the term is retained as implicit, and if j i k then it is treated explicitly.

We will consider the case j 0 k first. Upon expanding the explicit pain of the 

diffusion term, it becomes 

(μui ll ) Ik - C/L(ui

^ +u n riJn,k 

+., [ U m^ , +unr'mnlj`^mn

(13)

u i , m + un

rimnI 

r k

Note that the first term on the right-hand side of the equation is in conservation form.

Although the implicit equations are treated in quasi-linear form, for the purpose of

evaluating the explicit terms the most convenient representation is used. The

implicit terms, with j = k become (note there is no sum on j)

U

!

II u^ ) + (( ( u n ), r i. + un , • ri•n 

μui,nr'n + S^ ^ n ,)u + Ti. (/I ! J)1 °μ^ L N J J n μ 1 J 1J μ 1 J J (14)

where

Si 1n1 = r' m r mJ mn i - r^ r m 11 + 

r'i 

Jn,J Jn 

(15)

and

T^JJ - 2μ^ um, ! rm1 + um'J r

ml + M ,^ umrmJ+ un r, 1

(16)

m = i + 1, n = 5 - 2i and there is no sum on m and n

Since T^j involves velocity compo •.;ients and derivatives in directions gather than the 

i-th direction, the term is also treated explicitly.

Hence the total diffusion terms for the i-th momentum equation is given in

17
WGINAL PAGE: IS

quasi—linear form as Of POOR QUALITY

^[ C^ 9jj/i)u,Sjj +[Cjgjj(jL,j +2μrjE1) -kE C^gkkμI'kklul,j

1° 1 l J (17)

+ 8I [ 

Z C kgkkμ,k 

rkl ^ ul + [ 21 C k gkk/1 ikjk] u, J

+7 rCfgkkTkk + Z gnkCμuiln]lk+ Z gnl [μ uk lnJlk + ^ I ^ s lo g [i Lukl n=1 n=1 S^ J

k=1 t k

kOn n$1

where

1 0 0

8^ = 0 1 0

(18)

0 0 1

and 2 1 1

C = 1 2 1 (19)

12

and repeated indicies do not indicate summation.

Convective Terms

The convective term for the i-th equation can be written as

Puj u 1 I

j = P u j [ u

I ,

j + u m rI j J (20)

which becomes upon expansion

i 3 j( 2

j m i j3 n j i

P u ^ u Ij = jE { Pu u ,j+^ P u u mj +^1pu u nj} (21)

where the last term is nonzero only when j = 3. The full momentum equation is

obtained by substituting Equations (17) and (21) into Equation (11) and treating

the pressure gradient and velocity divergence as explicit terms. Since the pressure

is specified and impressed upon the viscous layer, its specification replaces the

normal momentum equation. Thus, the streamwise and spanwise momentum equations are

the only two retained.

18
ORIGINAL PAGE fJ

The momentum equation now becomes OF POOR QUALITY

aP E{IC19 tL]U'i

I 

jj+ LCi 9^μ^j+ 2 rjii) C+g kkμrkk 1ui'j

J I W

r j 3 k kk i i 3 k kk 1 j +si[z C i9 ^pk rkl^u + ^^ 1 C I g f Skjklu

Pujui ,j+^E P u j u m rmj + SbElFunu^ 

nij J

+kE l C 1 9kkTkk + Zl 9nk ^ μui In] Ik+ZlgniI/IU'InIIk

k#n nOl

+(I- 8k)g/rukliJk

Energy Equation

The energy equation employed here states that the stagnation temperature is

constant throughout

To = T + .2 q2 (23)

The general;Lzed tensor notation q 2 is given by

q 2 - u u - i j 9 j 

i!

where u and uj are the contravariant velocity components. Incorporating the assumptions made concerning the coordinate system employed, i.e.

= 0

g 13 - 923

we obtain

CI = (u i ) 2 911 + 2u1 U291? + ( u 2 ) 2 922 + (u3 ) 2 933

Neglecting the term involving (u 3 ) 2 with respect to the other terms, and defining

physical velocity components, i.e.

U P = u 1 h 1 WP = u2h2

19

(22)
0

we obtain

9

To T 2 (^pz + Wp + ht ' h 2 upw P (24)

This is the form of the energy equation used.

Equation of State

The equation of state assumes a perfect gas and is given by

P = Y pT (25) 

Linearizations 

The fallowing analyses assume a set of linear partial differential equations.

However, the continuity equation and the convective part of the momentum equation 

are nonlinear, containing terms that involve the product of density and velocity

components. in order t:o overcome this difficulty we employ the linearization procedure

(described in Ref. g and reviewed in Appendix A) to linearize the aforementioned terms

by Taylor series expansion about the known time level solution.

The density is first eliminated by employing , the equations of state and energy, 

and thereafter the resulting terms are linearized: These terms are of the following

form 

('04/8), +/3 (pn gn )an +/3+ (pnen)

vpn

+(3

n nen

P 

Tn

(u^)n[ +

ht 

h (u 2 ) n 

T. 

I (ut) n +^ (26)

2

+ Pn^nen 

(^ u2)n + 91 2 (ut)n] (,u2)n+Q

T

h, h'2 J

n nen n+f3

+ P 

Tn L 2(T n - T°) 

+ ^r-

P 

— 

2pn,,nen

where all velocity components are the contravariant ones, and 0 is always a velocity 

component, (ul , u 2 , u3 ) while ^ can be either a velocity component or a derivative of

a velocity component. In the case when only one velocity term appears, i.e., pug,

we set 

*n = n+$ = 1.

20

.^ .rte..
w

x:

It is important to note that in the preceeding equations the contravariant

velocity components are used. However, as noted in Ref, lei, it appears advantageous

to solve for the physical velocity components. Therefore, when the governing

equations are subsequently cast into a form amenable to the application of the LBI

scheme, they are transformed so that the physical velocity components appear.

The Turbulence Model

We treat the set of three-dimensional ensemble-averaged turbulent reduced

Navier-Stokes equations. Ensemble averaging permits the appearance of low frequency

(relative to the turbulence) time-dependent "mean" flaw. It is, therefore, necessary

to specify a turbulence model suitable for this problem.

The approach taken in the present effort assumes an isotropic turbulent viscosity

UT , relating the Reynolds' stress tensor to mean flow gradients.

Reynolds stress = r, ey = Re elk 3 Sj^ Q1 (27)

Using Favre averaging (Ref. 16), the governing equations then are identical to the

laminar equations with velocity and density being taken as mean variables and viscosity being taken as the sum of the molecular viscosity, p, and the turbulent

viscosity, UT•

At this point we require additional closure assumptions for the Reynolds

stresses, i,e., the evaluation of p T . There are a variety of approaches available,

from the simpler mixing length models to the more complicated one and two-equation

models. Since the intention here is to verify the code's performance in wall bounded

cases, we chose the mixing length models which have worked well in the past for simila

flow environments (Ref.17 ). The extension to more complex models could be undertaken

at a later time if warrn,,&Led. At that time, the LBI procedure that is used for the

solution of the momentum equation could be applied to the k and a equations.

Employing the Prandtl mixing length concept, the turbulent viscosity is

given as

^L

T c P l 2C*

(28)

where R is the mixing length, and e* is the mean strain, which for the general case

is related to the dissipation function.

e* _

(29)

21
ORIGINAL PAGE is

OF POOR QUALITY

where ^ is the dissipation function, in Cartesian tensor notation ( is given as

-- 2elj a l

CP = 2 e Ij e lj (30)

where

,

i

z 1 I'X' j + 

axl

Sul

In the case of two dimensions, -P does not reduce directly to jy^-j (where x3 is normal

to the walk.) without the additional assumption that the other components are small,

In nonorthogonal coordinates, with generalized tensor notation employed, an

analagous e%pression can be obtained for the dissipation function, i.e.

Cl

) = e l) UI (31)

j

where

e1j = g lm O j ( m + 

g jmU 1 IM (32)

After some algebra, 1V reduces to

2 e jj e jj (33)

where

I

e lj = g kl uj fm + gjm U IIm

For computation purposes, D can be expressed somewhat simpler, i.e.

CD = ( g IM U j 

1m + 

gjmU1 (m J g lk ukjj (34)

On expanding, and recalling that gm3 r 0 m 0 3, 0 can be expressed as a summation

^= E Uj) ^kl ^^Ukl(LZ9 9jnui^^+933^^g Ui I Uk) +gj u3l U3I

k j k j k j j I m 1k m k I lk 3 3 33 m j m

j (35)

+ g g33U31 

U31 3

22
ORIGINAL PACE 1

;

F OF POOR QUALITY

As in the Cartesian tensor formulation, 4, does not automatically reduce to the

dominant shear for two-dimensional boundary layers. Hence, provisions were made in

°k the computer code that on option retains only the dominant component of the strain,

14 The mixing length formulation is based on the two-layer model

inner layer _ k y

(36)

outer layer 10 ° jC0

a

where y is the normal distance from the wall, "b is the van Driest damping term

.2i = I — exp (-y+/A+)

and the ion Karman constant k .4 and the van Driest damping coefficient A+ R 26.0.

the dimensionless distance y+ is defined as

y+Z P

and the fxiation velocity u is taken as

u r = (Tt

In the outer layer ROD is given by

lam= ^,s

where d is the local boundary layer thickness defined as .995 u and X _ ,09.

This form of the mixing length has a discontinuity in its first derivative

at the matching point between k and R,o . An alternate form which avoids the discontinuity and varies smoothly between the two layers has been proposed by

McDonald (Ref.18 ), and is given by

1,,)tonh

^y (37)

of

In tre turbulent calculations to be presented we e ploy Eq. (37).

Spatial Difference Approximations

QR Operator Notation

In this section, implicit tridiagonal finite difference approximations to the

first and second derivatives and to the spatial differential operator are considered.

The QR operator procedure for generating a variety of spatial discretizations is also

introduced. As special cases, standard second-order finite differences, firstorder upwind differences, fourth-order operator compact implicit (OCI), fourthorder generalized OCI and exponential type methods are obtained. Since all these

schemes are of the same form (cf. below), a single subroutine which defines the

difference weights is all that is required to identify the method, while leaving

23
OF POOR QUALITY

the basic structure of the program unaltered. Subuequently, the results of numerical

experiments employing some of these schemes will be presented. The rationale for the

use of the QR approach in the present problem is discussed in dat* p ii in Ref. 3.

The QR formulation allows for ADI methods and permits the treatment of systems of

coupled equations, i.e,, LBI methods. Although variable mesh schemes can be employed

within the QR framework, it is believed preferable to use analytic transformations to

obtain a uniform computational mesh, hence attention is restricted to uniform menh

formulations.

The general concepts and notation will be introduced for two-paint boundary

Avalue problems and then the methodology will be extended to more general linear and

nonlinear parabolic partial differential equations in one dimension. The application

, fof the QR operator method to multidimensional problems is discussed in the section

pertaining to the LBI scheme.

Consider the two-point boundary value problem

L(u) - a(x)u XX + b(x)u x + c(x)u - f(x) (38)

with boundary values u(0) and u(1) prescribed, Derivative boundary conditions,

although not discussed here, can easily be incorporated into the framework of the

QR operator notation. Let the domain be discretized so that xj = (j-1)h, j = 1, 2,

J + 1, and U • u(xj ), F N ux (xj ), S ^ uxx (xj ) and h = 11J is the mesh

width. The numbering convention was chosen here to be compatible with FORTRAN coding.

Without loss in generality for a(x) # 0, Eq. (38) can be divided by a(x) so that

we may treat instead the following equation

L(u) - UXX + b(x)u x +c(x)u - f(x) (39)

where

b(x) - b(x.)/a(x), c(x) - c(x)/a(x) and f(x) - f(x) /a(x)

Substituting the finite difference approximations to the first and second

derivatives

Do

2h Uj U1+2hUj_I Fj " ux( x j ) + 0( h 2 ) (40)

Oho_ 

Uj

Uj-t `n2 +Uj+i . `'

^, 

R 

u xx (x j ) + 0( h2) (41)

into Eq. (39) and rearranging, we obtain

L(u) ^' S j + bj Fj I - bh J U 1-t +^cj h2 ^U j + _hT + 2 Uj+t a f j (42)

24
ORIGINAL PAGE IS

OF POOR QUP,l.1"('Y

or

^I

C 1 u j

_ i + Ch2c j - 2,uj + [I+ 2 ]uJ +i a h2fj (43)

where Rc j

hbj is the cell Reynolds number,

Equation (43) can be generalized by introducing operator format, i.e.,

rj' u

j-i + rj U j + r^ uj+i M h2(g j f j-1 + qj f + q j f j+1) (44)

where the superscripts (-) minus, (c) center, and (+) plus indicate the difference

weight that multiplies the variable evaluated at the (j-1), (j) and (j+1) grid points

respectively, and where the r j I s and qj 's for grid point j are functions of h,

b . l ._l, bi t bj+l , cj_l, 

C 

and 

cj+l. 

Comparing Eqs. (43) and (44) we can identify

the rj 's and qj 's, viz.,

rj I- Rcj P2 q j 0

rI rt h2 c j

 - 2 q^ (45)

r j - I + Rc j /2 qi 0

We now define the tridiagonal difference operators Q and R

Rjuj^ `

rj 

u j-i + rju + rj uj+i

fjJ ° g j f j-i +g j

fj + gJfj+i

Noting that L(u) = f and substituting Eq. (46) into Eq. (44), we obtain

R[uj] s h2Q[fj ] = h2 Q[ L(u)j]

(46)

(47) 

Alternatively by employing the inverse operator Q -1 and expression, for L(u) j can

obtained

I L( U) j = o" 'RU h2 J

25

.^....^..^.^—.,.,..^^.^.. _.. ^.+ .i. .,*^.' •arrt ."`^, +^ a -°-1 ^^.^+>?n►^;i. F y±'1r.^^,. -... „
For standard central finite differences Q = Q -1 = I, the identity matrix,

(the spatial operator is given explicitly in terms of U j _l , U and Uj*1 ) so that

nothing would appera to have been gained in obtaining Eq. (48). However, in general,

for higher order methods whereas Q is tridiagonal Q -1 is a full matrix. Hence,

Eq. (48) provides us with a convenient expression for the spatial operator for a

wider class of difference approximations. The formalism in Eq. (48) is also applicable

for first and second derivatives appearing alone (cf. Ref. 19). It should be pointed

out, however, that Eq. (48) is not the most general formulation since the compact

implicit formulas cannot be combined to yield a single scalar equation relating the

spatial operator *o the function values (Ref. 19).

In Refs. 3 and 20 a technique due to Berger, et al is described for constructing

fourth order tridiagonal methods which pos;,ess a monotonicity property as the cell

Reynolds number is increased, R c -­. A brief description of the method for

deriving generalized OCI schemes is given in Appendix C and the resulting Q and R

coefficients are presented in Table II.

Another family of schemes that can be expressed in QR operator notation are the

so-called exponential methods. The idea, originally due to Allen (Ref. 21)

(independently derived by Il'in (Ref. 22) and McDonald (Ref. 23)), and employed by

Dennis (Ref. 24) is to set the difference weights so that the numerical solution is

equated to the analytic solution for the locally frozen constant coefficient

equation. The QR coefficien"s of this exponential scheme is given in Table III.

This method is second order accurate for Rc = 0(1) and becomes first order accurate

as Rc -} - where the scheme reverts to first order upwind differencing.

Another exponential scheme which is uniformly second order accurate was 

developed by E1-Mistikawy and Werle (Refs. 25 and 26). The "exponential box scheme"

which is incorporated in their solution of the boundary layer equations with strong

blowing, is based on a spatial operator of the form given in Eq. (39). Berger, et al

(Ref. 27) derived the counterpart for an operator of the form given in Eq. (30), 

but with c = 0. The Q and R coefficients are presented in Table IV. Although this

scheme reverts to second order upwind differences as Rc + -, it does not possess a

discrete maximum principle while the exponential scheme of Allen (Ref. 20) does.

In Table V a centered finite difference scheme is presented which permits the

addition of artificial dissipation to the spatial operator when lRcl>Re so that

max

lRcl never exceeds Rc.max'

26
ORIGINAL PAC'ju 6

OF POOR QUALITY

Application to Coupled Nonlinear Parabolic Equations

Before considering the LBI technique, we discuss some of the limitations placed

on the QR operator scheme in solving a system of nonlinear parabolic equations.

Given a system of m nonlinear parabolic equations in m unknowns,

m

n +1 n

1 (un+I U 1 ) n+

I^ jQI *Q+R

At

Ni Q(u^,U21""Um,X1,X2,X3,t) s 0

j-1,2, ... d+1

where N

n+S is a quasilinear spatial operator, the QR formalism carries directly over

provided that for any equation only one independent variable is operated upon by the

differential operator. For example,

I a(U,W v) u t - uXX + b(u,v,w)u X + c(u,v,w)

is allowed since x derivatives of u only appear, while

I

a(u w v u t ' U XX + b(u ) v,w)u X + c(u,v,w) + d(u,v,w)wx

is not allowed since x derivatives of both u and w appear. The approximate form of

unsteady Navier-Stokes equations used here, when written in quasi-linear form, falls

within the class of allowable differential operators. Thus, for the problem being

addressed in the present study, the OCT schemes are applicable. Note that within the

splitting approach, non allowable terms in the OCI scheme such as dw x above, may be

split off and treated by a special implicit sweep. Provided care is taken and for

instance the Douglas-Gunn formalism is adhered to, no particular problem arises

other than the cost of an additional implicit sweep which is incurred.

Thus, multidimensional problems and/or more general equation forms can usually

be accomodated by a splitting procedure, which reduces the differential operator to a

sequence of one-dimensional problems which have the appropriate allowable form.

27
OF POOR QUALITY

However, as with standard finite differences, to avoid the cost of additional

implicit sweeps, special procedures must be applied to cross derivative terms, e.g.,

extrapolation or explicit treatment.

Linearized Block Implicit Scheme

Consider a system of nonlinear partial differential equations

Aft = + (49)

where ^ is a vector of unknowns and V is a source term vector which is a function of

x1 , x 2 , x3 and t. Extension to source terms which are functions of I are discussed

in Ref. (8). D is a three-dimensional nonlinear differential operator and the

matrix A appearin4 in the momentum equations is equal to pI where p is the density

and I the unity matrix.

Equation (49) may be centered abut the n+s time level, i.e. t o+s = (nW At =

nAt+$At = to+$At, and written

An+^3 L CP - (n

I /At ' 

fin +Q 4) 

+ it (50)

L

where 0 6 4 1 is a parameter allowing one to center the time step, i.e., a = 0

j

corresponds to a forward difference, s = 1/2 to Crank-Nicolson and S = 1 to a backward

difference.

ii After linearizing Equation (50) by Taylor series expansion in time about the nth

j time level by the procedure described in Ref. 8 to give a second-order linearization,

rr we obtain

Anr n+i— 

n

Ot = 1n[ n+9 n1 - yn5n + ^-n+/3

L I (51)

where 11 is the linearized differential operator obtained from D by Taylor series

expansion in time.

The difference between the nonlinear operator V and the linear operator Z is

defined as tin = D n _ Z n. At the intermediate level n + , (D n+S is represented as

n+/3 ' e 5"+ ( i-Q)"^,n (52)

Using these relationships and dropping the vector superbar for convenience a two

level hybrid implicit-explicit scheme is obtained

28
ORIGINAL PAGE 18

OF POOR QUALITY

A n ( ,j) n+I- CDn V A t _ 31n(4)n

+ I-CI)n) +,n 4)n + Mn fin + ,,n 

+/3

(53)

The vector ^n+$ represents all of the terms in the system of equations which are

treated explicitly. More about this will be said later, but for the moment note that

n+s may be approximated to the requisite order of accuracy by some multilevel linear

explicit relationship, or approximated by Vin with a consequent order reduction in

temporal ac,,uracy.

The i-erator ,44 is now expressed as a sum of convenient, easily invertible suboperators „^ < l + .41 2 + m. In the usual ,ADT, framework these suboperators

are associated with a specific coordinate direction. Further, it is supposed that 

these suboperators can be expressed in the QR notation introduced earlier. Writing

^n+S and Mn4) n as a single source term Sn+a , Equation (45) is written as

A [ -' ] /ot = a [ 1, + 1z +13

 ] [' - ^' ] + [1, +1 2 +1 3] d' + s (54) n n+1 n n n n n+i n n n n n n +,8

To solve this system efficiently it is split into a sequence of easily invertible

operations following a generalization of the procedure of Douglas and Gunn (Ref. 19)

in its natural extension to systems of partial differential equations. The DouglasGunn splitting of Eq. (54) is written as the following three-step procedure

An [ (D-- (pn]

/At. pj ( (D* 

-

4

)n ) + [1; + +I 1n ] ^n + Sn+/3

An[Cp*.- 4)n]/At

= ,.2

[Cl) CD n]_+ ^_1 [,

)«*_ (n] + [lin..,I2+1 

a]`Dn + S

n'+R (55)

An [ 4)* ** 4)n] 

/,t 

= p1^ [4). -d)n] + j [

q)* m _ 4) n] + 

Qja [

(DM» * 4)n]

+ p +1? 2 +,03 ]CD n + Sn +^6

29

_. _ _
OF POOR QUALITY

which can be transformed to the alternative foam

[An _A,,,n[^*_,)n] 

_ At[.,, n +1z +,L°3](D n + AtSn+p

[An otQ12] [ex—CDn] = An[ e_ CDn] (56)

[An - At 913 ] [cb*xx-(Dn] = An(4)**_Cl)n]

If the intermediate levels are eliminated, the scheme can be written in the so-called

factored form

( An -QAt-o XAn) 1(An-(3 .At22nXAn)'I(A'LPAt1°3n)(^n +I -(,n)

(57) At( I,n+12n+23n)^n+AtSn+p 

At this point it becomes necessary to consider the structure of the operators

it ,` 2 and Z 3 . It will be recalled from the one-dimensional scalar problem

that use of the QR format greatly facilitated the introduction of a wide variety of

spatial difference formulae. It follows that in the extention to multidimensions

undertaken here, the use of the QR formulation results in the appearance of the

inverse operator Q-1 with the sub-blocks of the Z 1' Z2) Z 3 operators. In

order to implement the scheme the inverse operator Q-1 must be cleared. Accordingly,

the scalar operator Q is generalized to the vector operator Q i with (diagonal) subblocks Qji . In this generalization 3 = 1, 2 apply to the momentum equations and j = 3

applies to the continuity equation. The i subscript is associated with the coordinate

directions of the Z i operators. The discretization results in one diagonal sub-block

for each grid point for each of the three Qi . Each intermediate step of the algorithm

is now premultiplied by the Qi associated with the Zi implicit operator. Writing the

product operator Qi Z i as Li , the inverse operators are thus removed and the scheme

is written, once again dropping the vector superscript for convenience

10,An - AtPL 1 n][ *- Vin] =AtL 1 (pn + AtQl[12n+13](Dn+AtQ1Sn+'B

[Q2An - At QL 2] [ 4)** = QQ2 An [SD * - fin]

[Q3 An - AtPL 3 ] [0***- 4)n]

= Q3An[0* x _ fin] (58)

12^n

= ^2 -I-^2^jn 13n(jn _ Q3 -1 R3 (Dh

4)n +I_ (D***-+0(At3)

30
OF POOR QUALITY

With the removal of the inverse operator Q -1 , the question of the proper intermediate level solution boundary conditions can be addressed. As is pointed out by

Briley and McDonald (Ref. 2), the proper intermediate level boundary conditions may be r^ 

derived by running through the intermediate steps in reverse order. Defining a

Y

boundary condition operator Bn after linearizing the appropriate physical boundary

condition by Taylor series expansion in time as

Bln ( (,pn +I

9 tt, on

and applying this operator to the algorithm defines the boundary conditions as

f,

Ban o3 An L ^^* _ ^

n J w L 83 ^3 qn "

At J3 B Lin [^

_ ^n^

(59)

82 02 qn [ * - fin

]

= ^ 82 ^ 2A n At G B2 L2n 1[ 4) **_ ^n1

Note that unless B3^Li commute (an unlikely event except with Dirchilet boundary

conditions, where B = I) the exact boundary conditions cannot be derived. A number

of possible strategies are possible at this point aimed at various levels of approximation to BiL. For the present, the term A ts BiLi[(p - ,n) is neglected. This

introduces an error of order 0[At(^ -fin)) into the solution but note that this error

disappears at steady state where 0*** = 0** = 0. Neglect of the AtsBiLi [0 - (P n I term

is, of course, equivalent to applying the physical boundary conditions on the intermediate level va i,bAes.

This completes the general derivation of the algorithm and attention is now

given to the specific forms of the Li operators including the rather special form

of the component operator for the continuity equation.

It is worth noting that the operator D or ,. can be split into any number of

components which need not be associated with a particular coordinate direction. As

pointed out by Douglas and Gunn (Ref. 19), the criterion for identifying sub-operators

is that the associated matrices be "easily solved" (i.e., narrow-banded). Thus,

mixed derivatives and the complicating terms which might inhibit the use of OCI can be

treated implicitly within such a framework, although this would increase the number

of intermediate steps and thereby complicate the solution procedure.

n 11 An inspection of Eq. (58) reveals that only the linearized operators L 1 , L2 an

L3 appear. Indeed, the computer code employs this feature by evaluating these thre

31

.. a°"a k .. ,,'^ '3 ..-.;...;r 7: '__ .../se.tr oarl^,,..+r^.s.^ _ +.!:'.K"^c '•?'z- » _ n:^ii
operators before the first sweep, storing them and accessing them as needed in the

subsequent three sweeps. In addition, the terms arising from the nonlinear terms are

immediately absorbed into Sn+a as they appear, allowing for an efficient evaluation of

the terms in the differential equations.

'l The spatial operators appearing in the differential equations are Z i, 2 and

3 must be identified at least formally in order to isolate the coefficients that

are to be used in the construction of the Q and R operators. The operators

1 2'

3 can be represented'in standard form at each grid point, i.e.,

1"^ " = U u ^^,^^ +U2^^,i + a

3 ^^ + f

O 14 R + q 5^3 (60)

a1 In Eq. (60) the first subscript of 0 indicates the velocity component (associated

with the corresponding direction and '" , " indicates a derivative. The subscripts of

the aij refer to the direction (i) and the term in the equation (j) respectively.

Note that the equation is in quasi-linear form, since the coefficients of the

derivative operators need to be identified, for use with the QR operator technique

employed here. Alternate schemes have been proposed by Leventhal (Ref. 28) for

equations in conservation form but are not considered here. In the following section,

a description will be given of how this entire operator is discretized by employing

the QR operator format, and how the discretization is incorporated into the LBI

framework in order to solve the system of equations (58).

The continuity equation is considered first. Since it is a first-order partial

differential equation it does not have the standard form of Eq. (59). Furthermore,

in the linearization process p has been eliminated in favor of the u velocity

components so that the continuity equation has become an equation for the three

velocity components, and not density.

f

1i An inspection of the system of equations under consideration reveals that

?j substantial savings can be realized if the equations are partioned appropriately.

f

{ This is in keeping with the observations of McDonald and Briley (Ref. 8) who noted

that skillful partioning of the resulting matrix can lead to significant decreases in

computation time. Due to the use of a boundary layer coordinate system, the normal

velocity appears only in conjunction with terms associated with the normal "3"

direction in the two momentum equations. Hence, in the first two sweeps one is

required to solve only for the two corresponding velocity components the streamwise

and spanwise momentum equation without the need of considering the continuity equation.

However, on the third sweep where all 3 velocity components appear, one must solve

all 3 equations. This strategy reduces the solution procedure to the inversion of

32
two 2 x 2 block matrices and one 3 x 3 block matrix rather than three 3 x 3 block

matrices which lends to a substantial reduction in computation times. If the full

Navier-Stokes equations were considered (including a normal momentum equation)

the aforementioned partioning could not be applied since the normal velocity would

appear in all three sweeps.

The question that arises is how to appropriately split the continuity equation,

since it is solved only on the third sweep. Here again the Douglas-Gunn formulation

leads to the appropriate choice. The continuity equation written in conservation

form is,

ap + i aX [

^pu i = o

After linearizing and eliminating p, the increment form is obtained

A" Au"*' + BnAWn+I+ Ato 

a s 

[VnAnAun+1 + vnBnAwn+1+ pnAVn+i}

-- Of c^X Jpu i ^n + At)9 a r (pn + n An ? dun+i+ (unBn )1^wn+i^

+ 

AJ^ axe C(pn + wnBn)Awn+1+(wy)Aun+1l

where all, the veloicty components are the contravariant components u = u l , w = u2

and v = u3 . J. is the Jacobian and

n

An

Tn [ 9 1 1 Un + 912wnJ

B = n _ P

-rn 1922wn + 912 

n

J

(61)

(62)

By employing the Douglas-Gunn procedure, Equation (62) is represented as a tl

sweep equation, and a consistent approximation is obtained to the continuity

equation, i.e., the x1 derivative term is evaluated at the * level and the x 2 der!

33
OF POOR QUALITY

term is evaluated at the ** level. The values of the intermediate derivative terms

are obtained after the solution of the first two sweeps of the two momentum equations.

Note that these terms do not contain the normal velocity, The equation can thus be

written in s,pnbalie form

n n 4 n n+i At/3 a

A Au + B Ow + [ J n n n+, n n n+i n n+i 11

{A v Au + v e ^w +p Ov J 6 X3 IJ

(63) 

n S

^Q At ar a" r^^t

aZ [J l J^ [ J { }, J ax 

ax J

Since the only term involving v is in the x 3 derivative term, one can directly

integrate the equation with respect to x 3 , i.e.

Ot f [ An u n+I + Bn L1wn+^ dx3 + At [v n. A n 0u n+i + v n B n Own+1 n

+ p Ov n+I

X 

1

dot [*

Sn d 3

`

X^ J J

X

The next section describes how this is done very easily via the QR operator scheme. 

The concept of integrating directly the continuity equation is not now. Davis (Ref. 29) 

in his coupled procedure for the solution of two-dimensional steady boundary layer 

equation used a trapezoidal rule to integrate the continuity equation. Weinberg

(Refs. 30 and 31) also used a fourth order Simpson integration scheme to solve the 

compressible boundary layer equations. Such procedures are stable and offer a viable 

alternative to approximating the derivatives by finite differences. Note that

conceptually the continuity .equation in integrated form is treated on each sweep

of the Douglas-Gunn splitting, although in actuality this can be viewed as having

the same form as each sweep and the integration operator can be incorporated into

the Z and D difference operators, and as a result the stability and consistency

of the original splitting is retained.

Implementation of the LBI Scheme Employing 

the QR Operator Technique 

Consider the third sweep of Equation (57) in which both momentum equations and 

the continuity equation are solved. The momentum equations are in the form

34
ORIGINAL. PAGE 19

OF POOR QUALITY

j An _ 86t.1

41A CD

l ,J (64)

where p9)*** is the column vector of unknowns, u, v, w. Here it has been implicitly

assumed that the equations have been appropriately normalized and that the Contra

variant velocity components have been suitably transformed into their physical

components. Employing physical. components, (cf. Ref. 32) leads to a better behaved

solution since these components are not unduly influenced by geometrical variations.

For the streamwise momentum equation one obtains

^ A^"(*= A 

u ► 33+ U 23 AU, 3 + 033 A U + Q43 AW + a AV (65a)

while for the spanwise momentum equation one obtains

In Ad)= Aw + b Aw + b Aw + b Av + b Au (65b) 3 >33 23 13 33 43 'S3

where superscript *** has been omitted from pu, Av and 4w. Now in Equation (65a),

the first three terms on the right-hand side are approximated by the operator

equivalent 

AU,33 + ^23 AU,3 + a AU =

01

1 

R2

AU (66)

33 AX3

so that

n l *% QI ^RI^U

13 + 043AW + 0536V LAX (67) 

Similar approximations are made for Equation (65b). After substituting Equation (66)

into Equation (64), and multiplying thru by Q, one obtains for the streamwise

momentum equation 

^QIPn _ jGXR 1 1 AU - J8AtQ1 a 43Aw^/3AtQ1053 6V =. QlpnAu,*' (68)

where a = At/Ax 32

Similarly for the spanwise momentum equation, one obtains 

IQ. ^ BX R 2] AW — QA1C 436V - Rat "7"53 ^` U - Q2 pnAW*'^ (69)

The same type of procedure is employed for the continuity equation. Since the

continuity equation involves only first derivatives, they can be represented as

d QcRc 

dx 3- Ax 3

(70) 

35

# s `t YL' x ^y^CF ^' .^ Qt^^^TwN^'3 rt-- -m R`^" .,.w^ _,... 
ORIGINAL PACE 19

OF POOR QUALITY

The operators Qc and Re are constructed to approximate the weights associated with

either a second order trapezoidal rule or a fourth order Simpson's rule, i.e,

Trapezoidal rule _

c

+

q c - ^ qo- 2 r GGr^ = 0 r C - -I , rG+ = I

Simpson's rule r

c 4 + _ I

qc - 9c 3 q 3 

r- -1 , rC - 0 r+ - , I

The discretized continuity equation thus becomes 

JAnOu+J©new+ OxR Q ^

1 Rc^JAa v nAu+JDo VnAw +JpAV] = RHS (71)

3

where RHS contains all the terms due to the linearization procedure and the terms 

evaluated at the * and ** levels, Multiplying thru by Q and setting w - At/Ax31

the equation reduces to

^Q e d e+ PWRc JAn Vn L1u + ^ Q c J Bn + BwRc JBn V n ] tw + QwFidlp n] dV = QC(RHS) (72) 

The resulting matrix derived from Equations (67), (68) and (71) becomes 

3 tridiagonal matrix (Q and R are tridiagonal operators) with each sub block t

on the form` 

C Q, pn -PXR i -/3,&t Q,a 43 , [- aL1tQ,a,,] Au Q (A Li**)

I

-/3^t Q2b53 Q2Pn -^^ 1 RAtQ2 b431 Qw = Q2(Ow

[Q JAn +awRcJAnv nJ CQc1Bn+f4wRcJ01IpwRcJ^ ] 6V c QC

This matrix is inverted by standard LU decomposition.

Boundary Conditions and Initial Conditions

The type of boundary conditions employed in the solution of the approximate form

of the Navier-Stokes equations are described in this section. On the body surface

36

.. a+ .t^' '^ +^ si*.fit ^Iiv .ISr^/^c !T+ t'ea1F ^. -wX .y
slip is prescribed for all the velocity components. At the outer edge of the viscous

layer the magnitudes of the streamwise and spanwi.se velocity components are also

prescribed, However, the value of the normal velocity component is not set, but

rather computed as part of the numerical solution as is the practice in standard

boundary layer procedures.

At the inflow boundaries, (upstream) velocity prattles are fixed, while extrapolation conditions are employed at the outflow boundaries (downstream). Further

discussion of this matter is given in the section of numerical results.

The intermediate boundary conditions employed on the first two sweeps are the

physical ones. For steady multidimensional problems, the imposition of physical

intermediate boundary conditions did not impair the quality of the solutions

obtained. These results are in keeping with the analysis of McDonald and Briley

(Ref. 2) for second order spatial schemes. For the unsteady cases considered 

physical, intermediate boundary conditions have also been used without any apparent

difficulty.

The question of proper intermediate boundary conditions for fourth order

methods until recently has not been resolved. Fairweather and Mitchell (Ref, 33)

developed nonphysical intermediate boundary conditions for a fourth order solution

of Laplace's equation, and showed that, in general, the use of noncorrected, i.e.,

physical boundary conditions leads to a loss in steady state accuracy £or'their

method. As pointed out by Fairweather and Mitchell (Ref. 33), their scheme is

inconsistent. It is this inconsistency that requires one to use appropriately derived

intermediate boundary condition6 in order to recover a steady state solution independent of time. However, if a consistent scheme were to be used, e.g. Douglas-Gunn,

then physical boundary conditions can be applied without any loss in steady state

accuracy, including fourth order generalized OCI schemes. These conclusions

generalize the results obtained by Briley and McDonald (Ref. 2) for second order

Finite difference methods to higher order schemes and to those schemes that can be

cast into a QR operator framework. However, using physical intermediate boundary

conditions is expected to decrease the overall temporal accuracy to first order.

The Computer. Code

The type of numerical algorithm employed as well as its formulation has a

marked impact on the structure of the computer code. One needs to consider both

the number of CPU operations as well as the memory requirements. Usually, the

number of operations can be reduced at the expense of increasing the amount of

storage. However, for three-dimensional problems the accessible fast (small core

37

o

,a .. ,^"' ?!...:aim.. ^ .«.J.v ..:.:®...T^ F+!nrs^< ++e «.
memory becomes a severe limitation oven in the case where. the coda hate not been

modified to optimize the operation count.

The storage requirements for the solution ov the approximate form of the threedimensional Navier- Stokes equations for even modest size grids (e.g. 30 x 30 x 30)

exceed the available small core memory of a machine like the CDC 7600, One must then

resort to either mass storage devices such as disks or slow access aemory (large core).

In using such devices both accese time and transfer rates must be considered,

When small amounts of data are being transferred frequently then access time becomes a

significant factor. Therefore, a combination of strategies must be employed in raider C

to optimize both access time and transfer rate.

An investigation of the operation count of the LBI scceme in conjunction with

the QR operator technique leads to thr, conclusion that the most significant fraction

of time is spent in computing the matrix coefficients, i.e. the linearization coefficients and difference weights, This amount far exceeds the time required for the matrix

inversion. Hence, it is worthwhile to optimize the calculation of these coefficients,

and if possible store their values. This procedure was accomplished by storing the

operator coefficients Z 2a nd Z3 as they were computed in the first sweep on the

right-band side of the differential equation. On the second and third sweeps, n2

and Z 3 were accessed respectively and were not recomputed. It was for this reason

that the formulation of the LBI scheme referred to the linearized operators Z 2's

instead of the D's on the right.--hand side of the equation.

In order to minimize data transfer, the code was constructed to have in memory

one plane of data at any given time. Recall that 1-2 planes are parallel to the

surface while the 3 direction is normal to the surface. During the solution procedure,

which is described below, the first and second sweeps are conducted on 1-2 planes,

evaluated in sequence in the 3 direction from the surface to the outer edge. The

third sweep is conducted on 1-3 planes evaluated in sequence in the 2 direction.

A problem arises in two dimensions during the first sweep, in that only the '1'

or streamwise direction is required yet 1 1-2' planes must be solved sequ,,ntially. This

means that only '1' lines are solved in the first sweep since the 1 2' direction is

essentially passive, and thus leads to unnecessary data transfer. In order to

alleviate this inefficiency the code was modified so that in two dimensions the

equations would be solved 'in core' in the 'l-3' plane on both sweeps. At the same

time, additional modifications were introduced that resulted in a significant speed

up of the cede. These include specialized hard wired 2 x 2 and 3 x 3 block matrix

multipliers and inverters, and an efficient addressing routine for spatially

deper'Qnt variables. As a result of including these modifications, the CPU time for

38
r

two-dimensional problems was reduced by nearly an order of magnitude, 'file run time

for a two-dimensional turbulent calculation is .00063 sec per grid point per time step

and for a laminar two-dimensional calculation is .00055 sec per grid point per time step.

Although there are still additional modifications remaining to be incorporated into

the code, it is believed that the major gains in efficiency have already been obtained.

The general structure of the computer code will now be described. After the

input section and the initialization of data e.g. geometry, grid transformations,

flowfield,, etc. the actual construction of the difference operators is begun. The

first derivatives of the velocity components and viscosity are obtained for the entire

flow field and stored for ready access when needed for the computation of the appropriate terms in the governing e quations. Thereafter the terms that are to be treated

explicitly are evaluateu and ab,, rbed into the function Sn.

The operators ,441, Z2 

and Z 3 are then computed. These are used to

evaluate the appropriate Q and R coefficients which are then stored for easy retrieval

during each of the ADI sweeps.

In the first sweep the matrix resulting from the application of the Z,

operators for the streamwise and spanwise momentum equation is solved as a 2 x 2

coupled system. The solution of this system, the * level quantities, are then used

to construct the right-hand side of the second sweep equations and to evaluate the

appropriate * level term in the continuity equation. At this point the ^„^ 2 operator

is accessed and again a 2 x 2 system of equations for the streamwise and spanwise

momentum equation is solved. The ** level quantities are then used to construct the

right-hand side of the third sweep equations as well as the appropriate terms in the

continuity equation. For the third sweep equations which consist of the two momentum

equations and the continuity equation, the Z3 

operator is accessed from memory. The

resulting 3 x 3 system of equations is solved for the three velocity components.

After the primary variables are evaluated, the theromodynamic quantities, density,

temperature and viscosity are computed. The procedure is then repeated at the following

time steps.

It is noteworthy that the scheme just described operates on vectors, i.e. lines

of data. Therefore, it could show promise for vectorized machines.

39

14 ' ^ f*4^ ^ { '•^'Y ^ I^iZ::IO^ 1fjR+J^' ^ aV^TR ,! -er.. Y ^ ^ _ i
Numerical Results

In this section we describe resu7,ts of numerical computations that were obtained by exercising the computer co0e discussed previously. The objectives of

these calculations are; < demonstrate the viability of the code in performing two

dimensional steady turbulent flow calculations as well as in computing two dimensional

unsteady flows, to compare second order finite differences with fourth order methods,

to discover any limitations of the solution procedure and to indicate what modifications

would be recommended to improve the code's performance.

In meeting these goals initial calculations were performed with the fourth order

generalized OCI scheme, Model one and two dimensional problems were considered,

progressing to the solution of the laminar boundary layer equations. These calculations which included both the Blasius and Howarth flows were used to validate the

fourth order generalized OCI option as well as the code modifications described in

the previous section.

At this point the choice of the specific calculations to be considered was

addressed.. The first set of calculations consisted of two turbulent cases; zero

pressure gradient (Weighardt Flat Plate (Ref, 34)) and adverse pressure gradient

leading to separation (.Newman Airfoil (Ref. 34)). These cases were chosen due to

the reliable experimental data that are available with which comparisons can be made.

The second set of calculations treat the laminar two dimensional oscillating flow

over a flat plate. Here again there is abundant information (theoretical as well as

computational) with which to make meaningful comparisons.

Two-Dimensional Steady Turbulent Flow

Weighardt Flat Plate

The first case we consider is the turbulent flow over a plat plate with zero

pressure gradient. Flow conditions were specified to match with the experimental data

of Weighardt (Case 1400 in Ref. 34), viz.

ue = 108 ft/sec.

Re = 661,190/ft.

The computational domain was chosen with the inflow boundary located at

x = .6135 ft and the outflow boundary located at x = 13.1232 ft, while the outer

edge was set to a constant value of y = .25 ft for all x. In order to resolve the

sublayer two types of grid transformation were employed, a Roberts type (Ref. 11)

40
•

which is essentially a hyperbolic tangent function and an Oh type (Ref. 35) which

is constructed from a series of complimentary error functions. In the streamwise

direction a uri,.urm grid distribution was employed. The total grid distribution

consisted of 36 points in each direction.

At the inflow boundary the velocity profiles are required to be specified.

These were obtained directly from Weighardt's data by using a Coles' law profile

(Ref. 34 ) i.e.

= K In^ 

YU T ^ + C + 2 K^x) sin2 2 8^ Wall-Wake Law

T (73)

y+ ° u+ laminar sublayer

where

Y + = LUT u+ ° uUT Ur ° Tw /P

and Tr(x)/K is evaluated from the condition that u = u. at y16 = 1. In Eq. (73) the

constants K and C are set at .41 and 5.0 respectively. This leaves two free parameters,

T (or C f ) and 6, which must be chosen to completely specify the profile. In (Ref. 34)

tabulated values of C and 6 are given, which are evaluated from curve fits of the

experimental data as a function of streamwise position. The values corresponding to

x = .6135 ft are

Cf = .004138

d = .015472 ft.

Once the streamwise velocity profiles are obtained the normal velocity can be

determined by approximating 8u/ax and then integrating the continuity equation.

In the calculation procedure the flow was assumed incompressible and in the

streamwise direction the boundary layer option was employed, i.e. streamwise diffusion terms were dropped and a backward difference approximation was used for the

streamwise convective terms. In the direction normal to the wall both second order

finite differences and the fourth order G/OCI schemes were exercised. The full form of

the dissipation function is used in the evaluation of the turbulent viscosity since

in this case it did not have any perverse behavior. In the Newman airfoil calculation

to be discussed shortly, however, the use of the full dissipation function adversely

influenced the behavior of the solution. The solutions were obtained by marching in

time until a steady state was reached.

41
The second order and fourth order calculation required 42 and 43 iterations

respectively to converge, where a convergence criterion of c-10

-5

 was employed and e

corresponds to the maximum change in a primary variable (u,v) over a range of time steps.

The total run time for the fourth order G/OCI calculation was 35.28 sec on a CDC 7600

computer which reduces to .00063sec per grid point per time step.

The skin friction distribution, C as a function of distance along the plate is

presented in Figure 1. A comparison between the experimental data and both computations shows good agreement over most of the plate, with the major discrepancy appearing

near the upstream boundary. There are two reasons for this. First, there is an

r

incompatibility between the inflow profiles (obtained indirectly from experimental

data) and the velocity profiles downstream of the inflow plane, which are obtained

by the numerical solution of the governing equations. There is no direct method of

avoiding this inherent error other than refining the mesh in the streamwise and

normal directions in the vicinity of the inflow plane. This would allow any errors

that are generated at the boundary to damp out several grid points downstream.

Another method that could be used relies on completing two computations; the first

beginning upstream of the inflow plane and extending downstream of it, and the second,

the one of interest, on the actual computational domain. In the second calculation

for the inflow profiles, one uses the profiles computed from the first calculation at

the corresponding x location. Although such a procedure is not always applicable, it

was used successfully for the Newman Airfoil Case (described below).

The second cause of the discrepancy, and probably the more serious, was the lack

of resolution of the velocity profile at the upstream plane. Recall, that the top

boundary was prescribed at some fixed constant height for all streamwise locations.

The extent of the normal distance is based upon the consideration that at the downstream boundary the velocity profle should be totally contained within the computational domain. As a result at the upstream boundary the "boundary layer" occupies a

small portion of the entire normal extent, leading to a lack of resolution there.

This problem can easily be remedied by introducing a y16 transformation which accounts

for boundary layer growth and assures equal resolution in the normal direction for

all x locations. In view of the benefits of this type of transformation it is

recommended that the computer code be modified to include a y/d transformation

on option.

42
43

Despite the shortcomings discussed above, the errors near the inflow boundary

die out rapidly and for the remainder of the flow field computations are in excellent

agreement with the data. The velocity profile at x - 11,44 ft. (of Figure 2), also

shows that good agreement was obtained between the computations and experiment.

In Figure 1, the increased accuracy of the fourth order generalized OCI scheme

is evident. However, there were qualitative differences between the second order and

fourth order solutions. The streamwise velocity profiles obtained by second order

finite differences were smooth throughout the entire boundary layer. In contrast to

this behavior the fourth order generalized OCI solutions tended to have slight wiggles

near the outer edge of the boundary layer even though in theory the generalized OCI

scheme should not if certain inequality constraints are satisfied. Unfortunately,

near the outer edge these conditions are violated allowing for the observed behavior.

We believe that the wiggles are due in part to high order numerical differentiation

of the turbulent viscosity which ;possesses a sharp knee in the wake region of the

layer. It is well known that higher order derivatives of rapidly varying functions

introduce noise, and it is this noise which we believe we are witnessing. Laminar

flow calculations have demonstrated that oscillatory behavior is not exhibited near the

outer edge even for very coarse meshes. It is, therefore, felt that in order to obtain

turbulent solutions commensurate with laminar solutions additional investigations of

'smoothing' procedures for the turbulent viscosity should be undertaken, and appropriate

grid distributions be considered.

As a final note, the computation of C is discussed. Second order finite differences employ a three-point second-order one sided difference. A comparable formula

accurate to fourth-order requires a five-point one-sided difference. As a result of

the grid transformation employed, the formula encompasses points relatively far from

the wall, and hence the computed C will be in error. Instead, for the reported

results in Figure 1, a three-point, one-sided formula was used which gives more accurate

values than the five-point scheme. Such behavior is not uncommon in computing onesided derivatives (cf. Ref. 36).
Newman Airfoil Case

The Newman airfoil flow (case 3500 - series 2, Ref. 34) is considered next. it is

more interesting and more difficult than the flat plate case due to the adverse pressure gradient present which leads to flow separation near the trailing edge. Indeed,

t the only meaningful numerical results obtained were those computed by second-order

finite differences. The fourth-order calculations computed on the same grid were

inaccurate and exhibited excessive oscillations near the outer edge of the computational

i domain. An explanation is suggested for the observed behavior, and a discussion of the

computed results is presented.

{ In the computations a 36 x 36 grid was employed. In the streamwise direction

which extends from x = 2.958 ft. to x = 4.926 ft. the spacing was uniform, while in the

normal direction which extends from the wall to a fixed outer edge at y = .35 ft.,

a Roberts type stretching was used.' The stretching parameter was chosen so that at

the inflow plane the first grid point would be located within the sublayer (y * 1.75).

The Reynolds-number for the flow was 769,231 per foot.

Initial flow computations revealed that the results were sensitive to both the

external velocity distribution and the upstream velocity profile. Differentiating the

velocity data given in Ref. 34 directly for use in the computation of the pressure

gradient was found to work best. Numerically differentiating analytical curve fits

to the external velocity distribution introduced errors in the pressure gradient and

were thus discarded. In order to eliminate incompatabilities caused by a mismatch in

the inflow conditions and the numerical solution, a preliminary calculation was performed on a 21 x 36 grid in the domain 2.759 <- x <_ 3.2, and the computed profiles at

x = 2.958 were used as upstream conditions for the primary calculation. Since the

pressure gradient in the upstream region was mild, there was no difficulty in obtaining

a converged calculation for the 21 x 36 grid (the calculation converged in 21 iterations)

However, the "boundary layer" form of the dissipation function was used since with that

form the solutions behaved better near the outer edge of the computational domain.

In Figure 3 the computed skin friction distribution is presented and compared to

the experiments,a. data. The agreement is good for most of the airfoil with the exception

of the trailing edge region where the second-order calculation predicts separation upstream of the actual separation point. However, for the shape factor, H = S /6 (cf.Fig.4)

the agreement between data and computations is not as good, with the numerical results

x

underpredicting the data. A lower value of H signifies that d is too small, 6 is too

large or a combination of the two. In any event, the lack of agreement in H means that

the velocity profile shape is incorrect. Two effects were considered to account for

44
this discrepancy; lack of spatial resolution and turbulence modeling. As discussed in

the previous section, the inability to account for boundary layer growth leads to a

spatial resolution problem. This effect manifests itself particularly near the upstream

boundary where the boundary layer is thinnest. Since the major discrepeneies in this

case lie near the downstream boundary, where the boundary layer is thickest, poor grid

resolution caused by the lack of a y16 transformation is not a predominant factor for

the observed behavior. Of course, grid resolution could be a real concern, but only in

how it affects the numerical scheme in resolving regions of large gradients, i.e. the

knee in the viscosity profile. In order to investigate the effect of turbulence modeling, a series of calculations were performed using a modified version of Cebeci's

turbulent boundary layer code (Ref. 37) which uses ICeller's Box scheme. This code permitted us not only to investigate turbulence modeling effects, but afforded us with a

means of comparing our computational, procedure with a different method, and running a

series of mesh refinement studies.

The modified Cebeci code was checked out by computing the Weighardt case. After

successful completion of this calculation, the Newman Airfoil case was attempted. In

running the Newman Airfoil anamolies appeared which were not present in the Weighardt

case. The computed skin friction coefficient exhibited oscillations in the streamwise

direction. Since the "mean" C curve compared well with our calculations for x > 3.5 ft.

and since our intent was to investigate the effects of turbulence modeling and grid refinement, we did not pursue the anamolous behavior. Even though oscillation in C were

observed, we believe that they do not invalidate the conclusions to be drawn. The

calculations were performed with 36 points in the streamwise direction (matching our

calculations), and either 36, 60 or 100 points in the normal direction, which was

stretched by a logarithmic transformation. It is important to point out that these

calculations were performed in similarity variables so that the growth of the boundary

layer was taken into account.

The results of the Box Scheme calculations were rather surprising. For the

calculations using either Cebeci's two layer mixing length model (Ref. 14)

or McDonald's model, the C distribution varied insignificantly from our second-order

finite difference solution shown in Figure 3. For the comparison, the oscillating

part was neglected and a mean C curve was considered. Neither was there much

improvement with the use of a finer mesh indicating that 36 point is probably sufficient

to resolve the flow. The major differences between Cebeci's model and Equation 37

are that in Cebeci's model, the outer layer mixing length is proportional to 6 , while

in Equation 37 it is proportional to 6, and in the Van Driest damping Cebeci includes

a term that accounts for the pressure gradient.

45
The effects of these differences were considered next. In order to investigate

the effect of pressure gradient terms in the Van Driest damping, the correction term

was included in McDonald's model (Eq. 37). The results of the calculation with 36

points in the normal direction revealed that there was a significant shift of t

curve upward toward the data, with much better agreement obtained with the data.

' However, the oscillation still persisted. Why that term did not have a similar effect

in Cebeci's mixing length model we cannot answer, and therefore, we believe further

investigation is warranted.

In comparing the shape factor, the H distribution exhibited different behavior.

These results are shown in Figure 4 which presents H as a function of x. It is apparent

that the calculations using Cebeci's mixing length give better agreement with the data.

Since the skin friction distribution did not compare as well, the good agreement of the

shape factor with data must be related to the overall shape of the velocity profile.

Hence, the velocity profiles obtained from Cebeci's model and from Eq. 37 were compared

with data at x = 4.509 ft., (cf Fig. 5) to discover how the turbulence model effects

the shape and thus H. Although the profile using Cebeci's model fits the data more

} closely, the shape of the velocity profile appears to be in error. Therefore, the

good agreement for H shown in Figure 4 may be fortuitous. Further investigations

beyond the scope of the present effort would be necessary to draw additional conclusions.

The following conclusions can be drawn from the calculations employing the Box

Scheme:

(a) Mesh resolution in the normal direction was not a significant factor in

0 the observed results.

(b) Choice of turbulence model can have a. significant effect on the solutions.

In view of the above, the inclusion of a one-equation and/or a two-equation turbulence

model in the computer code is recommended in any future effort.

Until now, we have neglected to say anything about the fourth-order calculations.

The results obtained for a 36 x 36 grid exaggerated the mild oscillations in the

velocity profiles observed in the Weighardt case. The inequality constraints required

for nonoscillatory behavior of the generalized OCI scheme were violated. In order to

satisfy these constraint conditions, more grid points, better transformations or both

are necessary. It is, therefore, recommended that further investigation be performed

in this area.

Two-Dimensional Unsteady Flow

The unsteady flow case considered is that of a nonzero mean flow with a sinusoidal

unsteady component superimposed upon it, i.e.

46
47

_^

u e ( t) = u o( I +A cos wt) (74)

where u is the mean flow velocity, A is the dimensionless amplitude of the periodic

part and w is the circular frequency. For the case of laminar flow over a flat plate,

Lighthill (Ref. 38) obtained expressions for the skin friction coefficient in the limit

of low and high reduced frequency.

0.332 + A(0.498coswt - 0.849wX/uo sln wt ) wX /up <<

Of r- ^ (75)

2 x

0.332 4 A(wX /u O ) I/2 cos(wt + 7r/4) wX /uO » I

It is convenient in the flow analysis to calculate the phase angle between the

external velocity and a boundary layer property (i.e. skin friction coefficient) as a

function of streamwise location, x - x o . The derivation of the phase angle is

presented below.

Consider an unsteady flow with an external velocity field that is a function of

x and t

ue(x,0 =uOW (I+Acoswt) (76)

and which has associated with it a boundary layer property f(x,t). Denoting the

average value of a function by a super bar (-), the average external velocity at

point x is defined as

t +T t +T

`u(X(

T f ^ uO (x O ,t)df = T f Iuo(XO )] 11 +A coswd di (77)

t i ti

or

G ko) = u o (x o ) (78)

where T = 2fr/w is the period of oscillation ant t 1 is a reference time from which

averaging begins. Similarly, the average value of the function f becomes

t ^

+T (79) f ( X O ) = T I f(XO,0dt ti
Here we have tactly assumed that the flow has reached a pseudo-steady state and

initial, transients have died out. For the cases considered here, such a state is

reached after one period so that t l T.

The Function f(x,t) is now expressed as a periodic function with frequency w, Le o

f (x, t) = f (xo) + B coswt + 0(xp)^ (8Q)

where ¢(x o ) is the phase angle between u e and f at point x - x0 . 

In order to determine the phase angle, the following integrals are evaluated

T Jt ^uQ (xo,t) ., uo( xo)) • ^f X 0 , t) - f(xo)) dt (81) { 

cos [c^(x o )^ =

B A uo(xo) 

and

^ +T _

B = T J (fko ,t) - f(xo)) z dt (82)

t^

The character of the flow field is governed in part by the amplitude of oscillation, A. Ackerberg and Phillips (Ref. 39) have shown that for A sufficiently small

A < .3 no backflow will occur. However, for larger values, A > ,3 regions of reverse

flow will appear over the plate, predominantly near the downstream boundary. The

occurrence of reverse flow at the downstream boundary has serious implications on the

numerical solution of the governing equations, and will be discussed in greater

detail below.

For the case considered in this report, A was chosen so as to avoid reverse

flow from occurring, while the other input quantities were specified by numerical

considerations (see below). Hence, the input data employed in the calculation are

u = .10

0

A = .125

w = 5fr/2

XU = , l

XD = 3.6

Re 666,667

48
Before describing the particulars of the calculations, a general discussion

of the solution procedure, its features, and its associated constr4int6 will be given.

The solution procedure for the unsteady case follows along the same lines as for the

steady case, with the essential difference being that in the unsteady case, tile

solution is advanced in real time with the time step chosen from considerations of

temporal accuracy rather than rate of convergence. Whereas in the steady case, all

boundary conditions are time invariant, for the oscillating flow the velocity distributions at the upstream and outer edge boundaries change with time. Tli,rc is no dif -

ficulty in applying the outer edge velocity boundary condition since the velocity

distribution is simply that which is given in Eq. (74). The upstream boundary is

more troublesome, If the upstream boundary is located at x = 0, then at that point 

the boundary layer is of zero thickness and hence in transformed boundary layer 

coordinates remains fixed. In that case, the upstream inflow boundary condition remaind 

Blasius independent of time. However, if the upstream boundary is Located at some

small., but finite, distance downstream of the leading edge ) the velocity profile will

change with time. Since the calculations were conducted in the physical plane rather 

than in the transformed plane, the upstream boundary was required to be placed at a

small finite value of x. Hence, in order to account for the varying upstream boundary

layer profile a method used by Singleton and Nash (Ref. 40) was employed, viz. the

upstream boundary was scaled by the new edge velocity at each new time level. This

procedure fixes the values of the flow variables at the upstream boundary, and permits

the use of function conditions (necessary for well posedness) but as a result also

introduces errors there.

By solving in the physical plane, boundary layer growth could not be very satisfactorily accounted for. This resulted in having the outer edge fixed at y = 0.037 ft.;

the choice of this distance being predicated upon the necessity to accommodate the

boundary layer at the downstream boundary. In Figure 6 the computed phase angle ^ is

presented and compared with Cebeci's calculations (Ref. 41), and to the low and high

frequency predictions of Lighthill (Ref. 38). Both second-order and fourth-order

calculations are shown, for S = 1. The agreement of our computations with the other

predictions is very good for w > .8. For small x there is a discrepancy due to the

implementation of the upstream boundary condition described above, and the lack of a 

y/d transformation. This effect was studied in greater detail at a lower frequency,

w = 7/2 and were compared to the calculations of Murphy (Ref. 42). 

49

.. y.. s ^,It''* :,...u.rrr' ^=,^Y^'±L=rs^.r e..w:,uq%'+^ 'j- "tf^'-» ^•„- `* r 'a'.:

--
i

In Figure 7 second-order solutions with b w 1 and b - 112 (Crank. Nicolson) arc

compared to Cebeci's results. The Crank Nicolson calculation givers larger values

of ^ than both the results for d - 1 and Ref. (41). Since Murphy's and Prenter's

higher order calculations (Ref. 43) show the same trend, this would appear to

indicate the improved accuracy of the Crank Nieolson scheme.

As A is increased, Phillips and Ackerberg (Ref. 39) show that reverse flow

will occur during the cycle over some portion of the flat plate, including the

downstream boundary. Although unsteady boundary layers permit regions of reverse

flow, there is a difficulty in applying appropriate boundary conditions if at the

downstream boundary fluid is entering the computational domain rather than exiting it.

In an attempt to remedy this situation, Ackerberg and Phillips argue that since

disturbances travel at some finite speed then at a distance sufficiently far down-- 

stream the flow will not have felt the disturbances generated at the leading edge

and would, therefore, appear to be that which would exist on an infinite flat plate,

The solution to that problem is well known, i.e. Rayleigh flow. Hence, they Ruggest 

that the Rayleigh solution be set at the downstream boundary (function condition)

consistent with the corresponding x and t. 

We investigated the behavior of our numerical scheme for a problem where reverse

flout occurs, A - .3. The boundary conditions described in Ref. (39) were used for

the so;(,ution of the full equations and central differences with artificial. dissipation was tmployed in the streamr+ise direction. The calculations were run for two

cycles, but the transients had not as yet died out completely. The results which

we are not presenting at this timr: (since they are of a preliminary nature) indicate

that there was no difficulty in the computations even when there were large regions 

of reverse flow. Future work will be aimed at conducting a more comprehensive study

in this very important area. Additional effort will be extended on turbulent 

unsteady flows.

50
CONCLUSIONS

In this report we have presented results obtained by exploiting the computer cod

developed under a previous phase of the research effort to solve nn approximate form of

the time-dependent Navier-Stokes equations. Both two-dimensional turbulent and unstea

laminar cases are considered. The governing equations that are solved are more genera

than the conventional boundary layer equations, notably in the inclusion of streamwise

and spanwise diffusion terms, although the pressure is still imposed by the external

flow, ^s in conventional boundary layer theory. The computer code incorporates the sp

LBI scheme in conjunction with OR operator scheme that permits a variety of spatial di

ference schemes, including standard second-order finite differences, exponential type

methods and fourth-order OCI techniques. In the split LBI scheme, an implicit sweep i

performed in each spatial coordinate direction. A careful ordering of these sweeps

permits an uncoupling of the continuity equation from the system in the first two

implicit sweeps. Thus, on the first two sweeps the (tridiagonal) system block size is

reduced from 3 x 3 to 2 x 2 with a resulting cost savings. On the last sweep of each

time step all the equations in the system are linearly coupled and 3 x 3 blocks must b

eliminated.

Results of computation indicate that the procedure is viable „br more complex

problems of interest. Higher order methods can yield more accurate results although

care rust be taken fo turbulent flows when coarse grids are employed. For unsteady

flows, the method is extremely efficient as a result of the noniterative nature of the

algorithm. Future efforts will be aimed at incorporating a y16 transformation, which

will eliminate some of the difficulties encountered and will make the code more robust.

Further investigations of appropriate turbulence models is also recommended with the

one equation K -k model given priority.

51
ACKNOWLEDGEMENT

The authors would like to acknowledge the valuable assistance given by Dr. John 

Murphy in making available computer printouts of unsteady flow calculations conducted

by him.

52

1 =.S. ♦ r ^R4v. -_ +tii-! ^ /OCR ^Yar/.
APPENDIX A

Linearization Technique

A number of techniques have been used for implicit solution of the

following first-order nonlinear scalar equation in one dependent variable

^(x,t);

dq5 /d1 = F(0) dG(0)/dx (Al)

Special cases of Eq. (Al) include the conservation form if F(¢) - 1, and

quasi-linear flow if G(^) _ ¢. Previous implicit methods for Eq. (Al) 

which employ nonlinear difference equations and also methods based on two- 

step predictor-corrector schemes are discussed by Ames (Ref. 44, p. 82) and

von Rosenburg (Ref. 45), p. 56). Qne such method is to difference nonlinear

terms directly at the implicit time level to obtain nonlinear implicit 

difference equations; these are then solved iteratively by a proced'ire such

as Newton's method. Although otherwise attractive, there may be difficulty

with convergence in the iterative solution of the nonlinear difference 

equations, and some efficiency is sacrificed by the need for iteration. An

implicit predictor-corrector technique has been devised by Douglas and Jones

(Ref. 46) which is applicable to the quasilinear case (G = ¢) of Eq. (A1).

The first step of their procedure is to lineari ,^.e the equation by evaluating

the nonlinear coefficient as F(^ n ) and to predict values of ^n+1/2 using either 

the backward difference or the Crank-Nicolson scheme. Values for n+1 

are

+1/2 )

then computed in a similar manner using F(¢ n and the Crank-Nicolson scheme.

Gourlay and Morris (Ref. 47 ) have also proposed implicit predictor-corrector

techniques which can be applied to Eq. (Al). In the conservative case (F = 1), 

their technique is to define G(¢) by the relation G(0) _ 06(Q) when such a

definition exists, and to evaluate G(,

n+l )

using values for n+l computed by

an explicit predictor scheme. With G thereby known at the implicit time level,

the equation can be treated as linear and corrected values of ^ n+l are computed

by the Crank-Nicolson scheme.

A technique is described here for deriving linear implicit difference

approximations for nonlinear differential equations. The technique is based

53
on an expansion of nonlinear implicit terms about the solution at the known

time level, t n , and leads to a one-step, two-level scheme which, being linear

in unknown (implicit) quantities, can be solved efficiently without iteration.

This idea was applied by Richtmyer and Horton (Ref. 15, p. 203) to a scalar

nonlinear diffusion equation. Here, the technique is developed for problems

governed by R nonlinear equations in k dependent variables which are functions

of time and space coordinates. The technique will be described for the threedimensional, unsteady equations.

The solution domain is discretized by grid points having equal spacings

in the computational coordinates, Ay l , Ay 2 and Ay 3 in the 

y l , y 2 and y3

directions, respectively, and an arbitrary time step, At. The subscripts i, j,

k and superscript n are grid point indices associated with y l , y 2 , y3 and t,

respectively, and thus 

0 i

,j,k denotes ¢(yi, y^, y3, t n ). it is assumed that

the solution is known at the n level, t n , and is desired; at the (n+l) level,

to+l. At the risk of ^n occasional ambiguity, one or more of the subscripts

is frequently omitted, so that On is equivalent to Oni,j, k .

The numerical method employed is quite general and is formally derived for

systems of governing equations which have the following form:

dH(0)/3r=."b(fl+SM (A2)

where 0 is a column vector containing k dependent variables, H and S are

column vector functions of 0, and 2 is a column vector whose elements are

spatial differential operators which may be multidimensional. The generality

of Eq. (A2) allows the method to be developed concisely and permits various

extensions and modifications (e.g., noncartesian coordinate systems, turbulence

models) to be made more or less routinely. It should be emphasized, however,

that the Jacobian aH/aO must usually be nonsingular if the ADI techniques as

applied to Eq. (A2) are to be valid. A necessary condition is that each

dependent variable appear in on= or more of the governing equations as a time

derivative. An exception would occur if for instance, a variable having no

time derivative also appeared in only one equation, so that this equation could

be decoupled from the remaining equations and solved a posteriori by an alternate method.

54
The linearized diffcrvncc- upprovimation is derived from the following

implicit time-difference replacement of Eq. (A2):

(H n+I -H^)/fit =Q[2 ( fin+I)+ Sn +II+(^-Q)[_b `')+Sn) (A3)

J

a

where, for example, H

n+1 

= 

N(n

+1). The form of .2 and the spatial differencing are as yet unspecified. A parameter 0(0 s < 1) has been introduced

so as to permit a variable centering of the scheme in time. Equation (A3)

produces a backward difference formulation for S = 1 and a Crank-Nicolson

formulation for a = 1/2. 

The linearization is performed by a two-step process of expansion about 

the known time level t o and subsequent approximation of the quantity

(a, /at) n At, which arises from chain rule differentiation, by (0n

+l

_ fi n)• The

result is 

Hn +I = H n+(aH /dc,)n ( o

n+I - o n

) + 0(At)2

(A4a)

S n+1= S n

+(aS1,3^) n 

((p n+I_O n

) +0 / jnt) 2

(A4b)

.2(O n+I ) = .2)(,pn)+(d.^)/d¢,)(q5n+l_on)+0(At)2 (A40

The matrices aH/Dq and DS/ 10 4, are standard Jacobians whose elements are defined,

for example, by (aH /a¢,) qr = aH q /air . The operator elements of the matrix

a2 /a4: are similarly ordered, i.e., (a2 /aQ,) gr = a.2 q /ai r ; however, the

intended meaning of the operator elements requires some clarification. For

the q th row, the operation (a2 /a ^ ) n (4 n+1 _ fi

n ) is understood to mean that

g

{a/at q [¢(x,y,z,t)))nAt is computed and that all occurrences of (a¢ r/at)n

arising from chain rule differentiation are replaced by (^r +1

 - ^ )At.

After linearization as in Eqs. (A4), Eq. (A3) becomes the following linear

implicit time-differenced scheme:

55
(aNn/a¢)(0 n

+t _ ^,n ) /At = 2> (on) +sn + 13 (a..b 160 + a S n /d o )(O n+i_0n ) (A5)

Although Hn+l is linearized to second order in Eq. (A4), the division by At

in Eq. (A3) introduces an error term of order At. A technique for maintaining

formal second-order accuracy in the presence of nonlinear time derivatives is

discussed by McDonald and Briley (Ref. 8), however, a three-level scheme

results. Second-order temporal accuracy can also be obtained (for a = 1/2) by

a change in dependent variable to = H(4), provided this is convenient, since

the nonlinear time derivative is then eliminated. The temporal accuracy

is independent of the spatial accuracy.

On examination, it can be seen that Eq. (A5) is linear in the quantity

(¢n+1 _ 0n ) and that all other quantities are either known or evaluated at

the n level. Computationally, it is convenient to solve Eq. (A5) for

(0n+l _ mn ) rather than ¢n+l. This both simplifies Eq. (A5) and reduces 

roundof£ errors, since it is presumably better to compute a small O(At) change

in an 0(1) quantity than the quantity itself. To simplify the notation, a

new dependent variable $ defined by

9n (A6)

is introduced, and thus q)n+1 = 

¢ n+1 _ ^n, and ^ = 0. It is also convenient

to rewrite Eq. (A5) in the following simplified form:

(A+ Al Z) q n+t = D1 [.'2) (fin) +Sn) (A7a)

where the following symbols have been introduced to simplify the notation:

A= 6H n /a cP--/3Z1(aS n/ap) (A7b)

It is noted that .fl(^) is a linear transformation and thus .L(0) = 0. Furthermore if 2(t) is linear, then 2(^ ) = -62)(q)),

56
Spatial differencing of Eq. (A7a) is accomplished simply by replacing

derivative operators such as 8/8y i , o 2 /ay i ay i by corresponding finite

difference operators, D i , D1. Henceforth, it is assumed that .2^ and 2 have

been discretized in this manner, unless otherwise noted.

Before proceeding, some general observations seem appropriate. The

foregoing linearization technique assumes only Taylor expandability, an assumption already implicit in the use of a finite difference method. The governing

equations and boundary conditions are addressed directly as a system of coupled

nonlinear equations which collectively determine the solution. The approach

thus seems more natural than that of making ad hoc linearization and decoupling

approximations, as is often done in applying implicit schemes to coupled

and/or nonlinear partial differential equations. With the present approach,

it is not necessary to associate each governing equation and boundary condition

with a particular dependent variable and then to identify various "nonlinear

coefficients" and "coupling terms" which must then be treated by lagging,

predictor-corrector techniques, or iteration. The Taylor expansion procedure

is analogous to that used in the generalized Newton-Raphson or quasilinearization methods for iterative solution of nonlinear systems by expansion

about a known current guess at the solution (e.g., Bellman & Kalaba, Ref. 48).

However, the concept of expanding about the previous time level apparently

had not been employed to produce a noniterative implicit time-dependent scheme

for coupled equations, wherein nonlinear terms are approximated to a level of

accuracy commensurate with that of the time differencing. The linearization

technique also permits the implicit treatment of cou p led nonlinear boundary

conditions, such as stagnation pressure and enthalpy at subsonic inlet

boundaries, and in practice, this latter feature was found to be crucial to

the stability of the overall method (Ref. 49).

57
A

n N

ORIGINAL PAGG F,n

APPENDIX B OF POOR QUALITY

GEOMETRIC PROPERTIES OF AIRFOILS

In this section we consider the vector properties of curves lying on a surface

and relate them to their tansor equivalents. In particular a NACA four-digit airfoil

section is considered, and its geometric properties are computed in vector and tensor

form.

BASIC CONCEPTS

Consider a radius vector r drawn from an inertial reference frame 0 to a point

P lying on a surface ^ (cf. Figure B-1). The tangent vector to the curve on the surface

passing through P in the x i direction is given by

e = 67 (B-1)

i ax,
ORIGINAL PAGE V^

OF POOR QUALITY

The unit tangent vector is

(B-

where s i is the arc length in the x i direction. But from the chain rule, we obtain

ar ar' asL as, _

_

ax i

dsia x i ti dXi ei

By definition 8s i/8xi

hi , the metric, so that

(no sum on i)

t h (B-3) ax I 

The curve passing through point P in the direction e i has a curvature Ki

associated with it which is given as

atI

at ax, I ati

KI as i

axi as i

hi anti

I a I ar

h i art

hi axi

(B-G)

Ki

ahiar I a2 r

h i ax i

axi + hi 2

axi2 

The unit principal normal to the curve point at P is in the direction of and

denoted by nit while the unit normal to the surface N, at point P, is given by the

cross product of the two tangent vectors at point P

At x t2 = N sin

(B-5)

where ^ is the angle between t and t 2 , i.e.

cosh = ti •

I a I ar

f2 =

r' 

h l ax, h2 ax,

so that

sin= I - cos? = L 9 11 922 - 9122}

(B-6)

h h

t 2

59
where gll , 922 and gig are the metric coefficients.

and

h hi z 9 9i^ zz

The vector normal to both e and N is called the binormal, B i . Hence the curvature

K can be represented by two components, one in the direction of N , and one in the

direction of Bi (cf. Figure B-2). These components are called the normal curvature

Kn , and the geodesic curvature B g respectively

Kn = K - N s K cos 9

(7)

Kg = K- 8 =Ksine

where 6 is the angle between n and N.

n

n

n N

Kn

e

B=e,XN

6

K9 e out of page

Figure B-2

DESCRIPTION OF NACA FOUR-DIGIT AIRFOIL

Consider an inertial reference frame in Cartesian coordinates

xl = (xl , x 2 , x3 ) or equivalently (x, y, z)

Attached to the airfoil is another Cartesian coordinate system given by

xi

(xl , x2 , x3 ) or ( x , y , z)

(cf. Figure B-3).

60
YT

Y

yL

XLx 

Figure B-4

61

^, .>ti" '.b.^.r ^ ..rORtiG:` s^v..::..:,1% — ^e^-uc... -•,; —

^s ^

ORIGINAL PAC". 

C

lL

,)F POORQUALIV

Figure B-3

The coordinate x is inclined at an angle n (the angle of attack) with respect

to x. In Figure 4 a two-dimensional airfoil is shown where (D and Q denote the leading

and trailing edges, respectively.
ORIGINAL PAGE CS

I OF POOR QUALITY

The chord ITT has length

f

C (XT 

-XL)2 + (YT _yL)2)1/2

and in general can be a function of z the spanwise coordinate.

The thickness distribution for a NACA four-digit airfoil is given by

i

+ It = f 2CC Iap X +QI X + a 2+03X3 + a 4 41

where x=x/c (cf Ref. 50). For a cambered airfoil the thickness distribution is

added onto the mean camber line. In the following, we will assume a symmetric

uncambered airfoil with the mean line lying on the x axis.

Hence the designation of the airfoil which is considered is NACA OOXX,

where XX refers to the thickness ratio t/c. In order to obtain a single valued

y = (yt /e) function, the polar angle 0 will be introduced. (cf. Figure B-5).

(B-8)

(B-9)

Figure B-5

62
OF POOR QUALITY

The following relationships hold

x = C/2 COS +c/2 = 2.(I - Cos (P)

(B-10)

X= — = 2(i — COS c)

0 leading edge 0 ¢ w upper surface

R w trailing edge w ¢ < 2v lower surface

Hence, Equation (9) becomes

- y

=

1.20 g (

C

) (B-1) C

where g(y) is Lhe polar representation of the term in the brackets in Equation (B-9)

and t and c can be functions of z. Defining g(^) g(^)/,2 and t(z) t(i)/c(i)

yt becomes

"y = t(i) 9(o) (B-12)

The coordinate x can also be represented as a function of z and ¢, i.e.

x = c(2) (I - Cos h) = CUMO)

and conversely $ can be represented as a function of x and z

Cos- ^ C I - ^xz

The relationship between the (x,y) and (x,y) coordinate systems is

x

{ (-x -' L ) cos a - (y - yt ) si n a (8-13)

C c

Y

_ I

- ^ {(-x -XL

)sina + (y-yt)cosa C

alternatively

x ={x L +xcosa + ysina)

(B-14)

-c = c {yL+ycos a -x sins }

Given a radius vector from the origin of the (x,y,z) frame,

r= xiy+ jk+z

63
ORIGINAL. PACE 19

OF POOR QUALITY

for the particular system under consideration

"Z XL+ c(7) f (0) cosa + 1(z)g(¢.)sina

(B-l5)

y R yt + t(i)g(o)cosa - c(2) f (0) sina

Hence 79 = K ( z ^)

z z (B-16)

T =7(y' 0)

and the surface coordinates will be lines of constant ¢ and constant z.

In general on the surface we can write in symbolic Form

r -x(0,z)i +y(^,z`)1+zk (B-17)

The tangent vectors along the lines ¢ - constant and z - constant can be computed

by differentiating the radius vector with respect to the coordinate line, i.e.

ar ax A + ay h

t 1 = a^ - a4) ► ao j

t_ ar = a7 A + d—y A

2 a'z a zz az`

The unit tangent vectors are

A t1 X^i +yoj

( 8-18) t

i (t 1 I

^xy

A = t2 - Xi i + -i

j +k

Tt2 1

I

+X^ +Y^" (8-19)

The metric coefficients become

2 2 2

h 1 91 1

x0 + Y (B-20)

h 2 = 922= 1 +7Z 2 +7. 2

64
ORIGINAL '4 t"

OF PGOR QUALITY

and

912 -7 4, 7

2 +

if ti and i2 are orthogonal,

gig - 0

The normal to the surface is obtained from the cross product of the surface unit

tangent vectors

A n n

t i x t2 A N sin q

where ^ is the angle between ti and t2 . 

In order to obtain sink, the dot product of the unit tangent vectors is employed 

tl I "t2 -- cos 41

or

h2 ^x i + Yz

 A

. { k ] hi [ 

X^ I + 70 , - J f, ' ^`x + y`y 701

I h2 =

yielding the following relationship

I r

cos tP

h l h2 

I XI, X^+ z -y l

Employing the expression

sin q1 cost

and with home algebraic manipulation the desired result is obtained.

I

f 1 112

sin = h l h z l g li g ^2 g1 2 I

(8-21)

In the case where rp and z are orthogonal to each other, tp "/2 and g12 C 0, or

X 2 x0 + 'y y' p

or 

X i - yz = 0

The normal to the surface can now be computed

A n 1 N = —y0i + X0 _ (X^yz - -y^ Xz ^ A ] z (B-22) 9119zz-gig

65
ORIGINAL PAGE 19

OF POOR QUALITY

From Equations (B-4) and (B-20), the curvature can be obtained in the ¢ direction

(I directit.,i)

I kj=

ah l 17' 62F- a 2 ] 1/2 a 4 A A

—a + 73

+yo [7 4.3

-^ I j +7k]

+

h,2

12 [RA + Y A + f A

2

k

k

h 

1 4

(B-23)

7o 700 + 707001110, - 

1 +

T

+

h-I

Similarly K2 is given as

i ah 2 6F + aer

k

2 h 23 H az h2 aZ2

or.

k 2

4 + x., z

1+ (B-24)

h 2 h 2

The geodesic curvatures can now be computed

A

K K 8 -K (t X N)

9 1 C

where tl , K, and N are given by Equations (B-18), (B-22), and (B-23).

After some algebra, we obtain

K - 'x

9 h 1 4 h2S14

Similarly K becomes

92

K

02 4

z 7,17 ) + (To

+-go 7,1 , (B-26)

hi 

4 h 2 sin qi

66
67

We will now consider a NACA 4 digit airfoil and further assure X L is a constant

(non swept trailing edge). From Equation (15) we can obtain the appropriate

derivatives

R n XL + c('z) f(0) cos a + t(z) g(0) sina

x^ = c(-x)focosa + tg0 sin a

X OO = c cos y f00 + tsin agog

(n-z7)

Xz = ci f cos a + tz gsina

R

if ciz fcosa + tzz gsina

y = -y

L

+ twg(0)cosa - ca)f(0)sina

yo = gotcosa - focsina

-yoo = gootcosa - foocsina (B-28)

yZ = g Cosa tz - f sinaci

yf = gcosatzz - fsina cH

For a rectangular plan form c and t are constant so that

Z = xx if 0

yz = yiz = 0

Substituting Equations (B-27) and (B-28) into Equations (B-18) to (B-24) the appropriate

geometric quantities can be obtained specialized for a NACA OOXX airfoil.
ORIGINAL PAGE 13

OF POOR QUALITY

For the present code, the curvature terms are not used directly. Instead, the

metric coefficients, g ij , the Jacobian, ;7, and the Christoffel syr'-ls, 

rij 

(cf Ref . 4)

are used. These geometric functions are presented below.

Metric Coefficients:

9„ • hl 2 , 1 

2 +Y0

9220 h2 2 - Xi 2 + YZ 2 + I

91 2 X^X q + Y^Yz

9 33 6

I

J2 = 702 + y02 + W2

W = X0YZ YOXZ

Basis vectors

A n

ei x`^a +

e2 Xi s + yij +^C

e3 = -y

01,

+ W k)

and

(B-29)

(B-30)

el = el /h l , e2 = e2/h2 , e3 e3

Christoffell Symbols (27 components)

1

,3 c 4, 3 = r 3 = r 3 = o

13 '23 31 32

I 2

r33 = r33 = O

r, 3 = O, i = 1, 2, 3

68

(B-31)

(,B-32a)
ORIGINAL PAta -- G'_j

OF POW QUALI TY

r,i = ,I2 

+ WOO ) + (V1100

r1 2 = r2l _ J ^^ ^y1 ^^i - yz yZ) + tx^kOz + 70 yo

i )t

1,22 = ^2 It ony l + YO Yqq) + ( RO VE X Z yo) (j xi2 - yf y Z i )J

r1 3 - rai = 3 {(^^Xo XoY ) + [(0¢ (xO xZ + 4,Y, - ^(XZ xOO + y = 40)1}

r2a rat = - llwf (z O x f + yy^) ^(yf y + z= X i )^ +

r „ 

2 = WJ2

r2 = r22 JZ ^^Yz yo Xoz) (B-32b)

r22 : J2 `X^yzZ - y^XZZ)

r2 = r32 = J3 jW q0 yoo+ XoX oo ) - ^^(Xo2 ,} y0 2 )j

r23 = r32 = J3 {41(xX^Z - y^y01 ) - wZ (x + yo 2)1

r^ - 5 ^X ^ y^^ - yo x00 )

r1 2 = r2; _ \X Vol

1,22 = J, t X^ zZZ - y0 y Z Z )

These terms can be specialized for a NACA four digit airfoil by substituting

Equations ( B -27) and (B-28) into Equations (B-29) - (B-32).

69

T ._ .....___ ., ,_,-
ORIGINAL IMAGE is

APPENDIX C

OF POOR QUALITY

Generalized 02erator Compact Implicit Schemes

In this section a procedure for g-±nerating generalized OCI schemes is reviewed.

Given

L (U) =U XX +G(x)U X+CWU

an expression relating L(u) and u is sought in the form

Pa RUj = Q(LU)j + T j

(C-1)

where T. is the truncation error and Q and R are tridiagonal displacement

operators. The maximum accuracy attainable is fourth order, i.e., T^ 'L 0(h4).

Expanding Eq. (20) in terms of q^' c ' + and ri I

we obtain,

y I

T j = h RUj - Q(LU)j - he rj U j-, + r^ U j + rj i)j+i^

[q^ (Lu)j _ 1 + q^ (Lu) j + gj(Lu)

c2

j+I^

A Taylor series expansion *Melds for T.

Tj = T u(xj) +TJ U ( ' ) (xj ) + T2 u(e)(Xjj + Tjau(a)(xj)

4 (4) 5 (5) s (6) 5

(C-3)

+Tj u (x j ) +Tj u ( xj ) + Tj u (xj ) + O(h )

70
URIOINAL PAGE IZI

E OF POOR QUALITY

where superscripts im parentheses denote derivatives with respect to x, and

where

E

T^ = hz ^(rl + r^ + r^) + h 2 ( q ^ c j _ 1 + gjcj + gjcj+, )^

t

TJ = h ^(r^ - rj)-h(g^ b J-i + gJ bJ+gjbJ+i)-h2(gjcJ+i-gJcj-1)]

i TJz

2 (r1 + r q- + qs + q+) (c-4)

- h(q^ bj+l -q ) bj-1)' 

h2 

(gjcJ+,+gjcj-t)

h

TJ hY -2

y! (r J+ + (-I)vr J ) (v-I)!(ql bj+1 + (-I)Y-Lq, bj-1)

I

I - (v

-

2)i (q^ + ( - I)Y q^) + h2 (q^c j+^ + (-I) vg jc J- 1 ) "° 3,4,5,6

f1For second order central finite differences we set TO = T1 = T2 = 0. This

yields, when qc 1 and q = q^ = 0, the following relations

i 

r c - (r" + r +) + h2 ( g j c + qj c j + q+c,+, )

r^ - r J- = hbJ = RcJ

rJ+ + rj 2

which recovers Eq. (45), i.e.,

r` _ - 2 + hzc

r^ = I - RC /2

r t = I + RcJ/2

1

71

1n 
OF POOR QPAGEIs UALITY

To obtain the fourth order operator compact implicit scheme we again set

TO = T1 . T2 , 0 to obtain three expressions for r^' c ' + in terms of the

qi' c ' + , (note that q - , q+ V 0 and q c is not necessarily unity), i.e.,

r^ _ -(r^ +r-) +h2 (gj c 1- +gjcj +gjcj+i ) (C-5a)

r^ - r^ = Rc j-i qj + Rc j q^ + Rcj+iqj+ + h2(g j+cj+i + q. cj-i) (C-5b)

ri +rl = 2(qj +qj +qj) +2I Rc j +i q j - Rcj-i q i I + h2 (g j+ cj+ , + q j- c j _ i ) (c-5c)

Now T3 and T4 must be set to evaluate q-' c ' +. The standard Swartz OCI

method requires T 3 = T4 = 0,

I _ 1

6 ( r + - r - ) - —21 Rc j+i q j + Rcj-iqj ] - [ q j -qj J =0

- 

h2(g1 

cj+, - q, cj-I) (C-5d)

2I4 ( r+ + r-) - —6' Rcj +I q j - Rc j-f q j ] _ 2 q j + q j = 0

- 

h2(g j+cj+i + q j cj -I) (C-5e)

and results in a leading truncation error of the form (Au (S) + Bu ( ' ) )h4 . Substituting (C-5b) and (C-5c) into (C-5d) and (C-53) r+ and r^ can be eliminated and a

system of two equations in q^ and q^ with q^ as a parameter is obtained. The

parameter q

c

is proportional to the determinant of the system. The values of

- c

' +

^- q^' and r c +

j

'' are presented in Table I.

As shown in Ref. 20 a cell Reynolds number stability condition exists

for the Swartz OCI scheme, i.e., for Rc 12 nonrealistic or oscillatory

solutions will be obtained. In order to eliminate this restriction one can

relax the conditions T3 = T4 '= 0, and allow the coefficients of u (3) and u(4)

to be of 0(h4).

72
ORIGINAL PAGE i3

OF POOR QUALITY

By expanding q 3 ' c '

+

 in a series in Rc

3

F`

M=O

Ry 12 parameters, am' c ' + m - 0, 1, 2, 3, are introduced. The equations for

T3 = 0(h 4 ) and T4 = 0(h4 ) yield 5 linear relations, leaving at one's disposal

1 '6 free" parameters plus a factor.

These parameters can be set according to some criteria that would yield

certain desirable properties for the difference equations. The following

constraints are prescribed

qj+ >_O, q J >O , q c >0

C

q I bl 

>_ q l b )-I + q j b)+t

(C-7)

r,+ > r,- >_ 0

- r. r++rand h is sufficiently small so that

10b j - b I _ i - b j+1 >0 and 2+ hcj+,/bj41 >0

for J-2, R R -,J and c j <0

These conditions assure that R is diagonally dominant and Q is invertible for

j all Rc . Further details are given in Ref. 27. The significance of this

approach is that cne can construct a scheme possessing certain desired properties by employing a set of preassigned rules. This is contrary to usual practice,

>I

in which a scheme is chosen, and then its properties are determined. Although

the computational effort in computing the q and r coefficients is not cheap,

the actual percentage of the total cost is minimal. This has been borne out

f

in actual computations.

jThe q and r coefficients for the generalized OCI scheme described in

i Ref. 27 are given in Table H.

73

,. ^ ate• - ...
REFERENCES

1. Kendall, R.M., Bonnett, W.S., Nardo, C.T. and Abbett• , M.J.: Computational

Procedure for Three-Dimensional Boundary Layers on Aircraft and Aerospace Vehicles.

Proceedings AIAA 2nd Fluid Dynamics Conference, Hartford, 1975, pp. 113-121, 

2. Briley, W.R. and McDonald, H.: On the Structure and Use of Linearized Block ADI

and Related Schemes. J. Comp. Physics, Vol. 34, No. 1, 1980.

3. Weinberg, A.C. and McDonald, H.: Solution of Three-Dimensional Time-Dependent

Viscous Flows, Part 1: Investigation of Candidate Algorithms. SRA Report 

R79-90004, Final Contractor's Report,(NAS2-10016), 1979.

4. Weinberg, B.C. and McDonald, H.: Solution of Three-Dimensional Time-Dependent

Viscous Flows, Part 2: Development of the Computer Code, Final Contractors Report, 

(NAS2-10016), 1980. 

5. Nash, J.F. and Patel, V.C.: Three-Dimensional Turbulent Boundary Layers.

Scientific and Business Consultants Inc., 1972.

6. Lin, T.C. and Rubin, S.G.: Viscous Flow over a Cone at Moderate Incidence,

Part 2. Supersonic Boundary Layer, J. Fluid Mech., Vol. 59, Part 3, 1973, pp. 593-620.

7. Briley, W.R.: Numerical Method for Predicting Three-Dimensional Steady Viscous

Flow in Ducts. J. Comp. Physics, Vol. 14, No. 1, 1974.

8. McDonald, H. and Briley, W.R.: Three-Dimensional Supersonic Flow of a Viscous

or Inviscid Gas. J. of Comp. Physics, Vol. 19, No. 2, 1975. 

9. Blottner, F.G.: Investigation of Some Finite Difference Techniques for Solving the

Boundary Layer Equations. Computer Methods in Applied Mechanics and Engineering.

Vol. 6, 1975, pp. 1-30.

10. Howarth, L.: The Boundary Layer in Three-Dimensions-Part 1. Derivation of the

Equations for Flow Along a General Curved Surface. Philosophical Magazine, Vol. 7,

1951, pp. 239.

11. Roberts, G.O.: Computational Meshes for Boundary Layer Problems. Proceedings

of the Second International Conference on Numerical Methods in Fluid Dynamics,

Springer-Verlag, New York, 1971, p. 171.

12. Flugge, W.: Tensor Analysis and Continuium Mechanics, Springer-Verlag, 1972.

13. Eiseman, P.R.: The Numerical Solution of the Fluid Dynamical Equations in Curvilinear Coordinates. Air Force Weapons Laboratory Report AFWL-TR-73-172.

14. Cebeci, T., Kaups, K. and Ramsay, J.A.: A General Method for Calculating

Three-Dimensional Compressible Laminar and Turbulent Boundary Layers on Arbitrary

Wings, NASA CR-2777, January 1977.

15. Richtmyer, R.D. and Morton, K.W.: Difference Methods for Initial Value Problems. 

Second Edition. Interscience Publishers, New York, New York, 1967.

74

+..'.,.-......,...,.n
16. Favre, A.: Equations des Gaz Turbulents Compressibles. J. de Mecanique, Vol. 4,

pp. 361-392, 1965.

17, Kreskovsky, J.P., Shamroth, S.J. and McDonald, H.: Application of a General

Boundary Layer Analysis to Turbulent Boundary Layers Subjected to Strong Favorable

Pressure Gradients, J. of Fluids Engineering, pp. 217-224, June 1975.

18. McDonald, H. and Camarata, F.J.; An Extended Mixing Length Approach for Computing

the Turbulent Boundary Layer Development, Proceedings Computation of Turbulent

Boundary Layers - 1968 AFOSR-IFP Stanford Conference, Vol. 1, pp. 83-98, 1969,

19. Douglas, J. and Gunn, J.E.: A General Formulation of Alternating Direction Methods.

Numerische Math., Vol. 6, 1964, p. 428.

20. Berger, A.E., Solomon, J.M., Ciment, M., Leventhal, S.H. and Weinberg, B.C.:

Generalized OCI Scbemes for Boundary Layer Problems, Math Comp., Vol. 35,

No. 151, 1980.

21, Allen, D.N. and Southwell, R.V.: Relaxation Methods Applied to Determine the

Motion, in Two Dimensions, of a Viscous Fluid Past a Fixeu Cylinder, Quart.

J. Mech. Appl. Math., 8, 1955, pp. 129-145.

22. Win, A.M.: Differencing Scheme for a Differential Equation with a Small

Parameter Affecting the Highest Derivative. Mat. Zametki, 6, 1969, pp. 237-248 -

Math. Notes, 6, 1969, pp. 596-602.

23. McDonald, H.: A Novel Finite Difference Scheme Particularly Suited for Fluid

Flow Problems Containing Discontinuities, Scientific Research Associates Report

P'76-1, August 1976.

24. Dennis, S.C.R.: Finite Differences Associated with Second-Order Differential

Equations, Quart. J. Mech. App1. Math., Vol. 23, 1960, pp. 487-507.

25. El-Mistikawy, T.M. and Werle, M.J.: Numerical Methods for Inviscid, Viscous

Fluid Flows, Report No. AFL 77-9-34, University of Cincinnati, March 1977.

26. El-Mistikawy, 'T.M., and Werle, M.J.: Numerical Method for Boundary Layers with

Blowing - the Exponential Box Scheme. AIAA J., 16, 1978, pp. 749-751.

27. Berger, A.E., Solomon, J.M. and Ciment, M.: An Analysis of a Uniformly Accurate

Finite Difference Method for a Singular Perturbation Problem. Paper presented

at Second International Conference on Computational Methods in Nonlinear Mechanics.

March 26-30, Austin, Texas.

28. Leventhal, S.J.: The Operator Compact Implicit Method for Reservoir Simulation,

Proceedings of Fifth SPE Symposium on Numerical Resevoir Simulation, 1979.

29. Blottner, F.G.: Computational Techniques for Boundary Layers in "Computational

Methods for Inviscid and Viscous Two and Three-Dimensional Flow Fields". Fluid

Dynamics Institute, Hanover, New Hmapshire, 1975.

30. Weinberg, B.C., Leventhal, S.H. and Ciment, M.: The Operator Compact Implicit Scheme

for Viscous Flow Problems. AIAA Paper No. 77-638. Presented at the 3rd AIAA Comp.

Fluid Dynamics Conf., Albuquerque, N.M., 1977.

75

, A , Fig.. ,<-:...v 

•,.•'^^^^^i^^^.;

^

_ _
31. Weinberg, B.C.: Viscous Flow Calculations Employing a Fourth Order Generalized

Operator Compact Implicit Scheme, Paper 79-1468, Presented at 4th AIAA Comp.

Fluid Dynamics Conf., Williamsburg, Va., 1979. 

32. Gibeling, H.J., Shamroth, S.J. and Eiseman, P.R.: Analysis of Stang-Interaction

Dynamic Stall for Laminar flow on Airfoils. NASA CR-2969, April 1978.

33. Fairweather, G. and Mitchell, A.K.: A New Computational Procedure for ADZ

Methods, SIAM J. Numer. Anal., Vol. 4, No. 2, 1967.

34. Coles, D.E. and Hirst, E.A. (Editors): Proceedings Computation of Turbulent

Boundary Layers - 1968. A70SR-IFP Stanford Conference, Vol. I, Compiled Data, 1969.

35. Oh, Y.H.: An Analytical Transformation Technique for Generating Uniformly Spaced

Computational Mesh. Old Dominion University Research Foundation, October 1978.

36, Carnahan, B., Luther, H,A, and Wilkes, J.O.: Applied Numerical Methods,"

John Wiley & Sons, pp. 130, 1969.

37. Cebeci, T. & Bradshaw, P,, "Momentum Transfer in Boundary Layers," Hemisphere

Publishing Co., 1977, Chapter 8,

38. Lighthill, M.J., "The Response of Laminar Skin Friction and Heat Transfer to

Fluctuations in the Stream Velocity," Proc. Roy. Soc. 224A p. 1-23, 1954.

39. Phillips, J.H., and Ackerberg, R.C., "A Numerical Method for Integrating the

Unsteady Boundary-Layer Equations When There are Regions of Back Flow,"

JFM Vol. 58, Part 3, pp. 561-579, 1973.

40^ Singelton, R.E. and Nash, J.F.: A Method for Calculating Unsteady Turbulent

Boundary Layers in Two and Three Dimensions, AIAA J., Vol, 12, No. 5, pp. 590-595,

1974.

41. Cebeci, T. and Carr, L.W.: A Computer Program for Calculating Laminar and

Turbulent Boundary Layers for Two-Dimensional Time-Dependent Flows, NASA TM-78470,

March 1978.

42. Murphy, J., private communication.

43. Murphy, J.D. and Prenter, P.M., "A Hybrid Computing Scheme for Unsteady Turbulent

Boundary Layers," 3rd Symposium Shear Flows, 1981, pp. 8.26, 8.34,

44. Ames, W.F.: Numerical Methods for Partial Differential Equations. Barnes &

Noble, Inc., New York, New York, 1969.

45. von Rosenberg, D.A.: Methods for the Numerical Solution of Partial Differential

Equations. American Elsevier Publishing Co., Inc., New York, New York, 1969.

46. Douglas, J., and Jones, B.F.: On Predictor-Corrector Methods for Nonlinear

Parabolic Differential Equations. Soc. for Indust. Appl. Math., Vol. 11, 1963,

pp. 195-204.

76

.. ,..i ^`.^v '+ -L.ay.Y mod= - .,..00r:.^ e.^...r...wlh
47. Gourley, A.R., and Morris, J.L.: Finite-Difference Methods for Nonlinear Hyperbolic Systems, Math. Comp,, Vol. 22, 1968 ) pp, 28-39.

48. Bellman, R.E., and Kalaba, R.E.t Quasilinearization and Nonlinear BoundaryValue Problems. American Elsevier Publ. Co., Inc., New York, 1965.

49. Briley, W.R.,, and McDonald ) H.: Solution of the Multidimensional Compressible

Navier-Stokes Equations by a Generalized Implicit Method. J. Camp. Physics,

Vol. 24 0 No. k, August, 1977, p. 372.

50. Abbot, 1.11., and VonDoenhoff ) A.E.: Theory of Wing Sections, Dover Publishing

Co., 1959.

77

+ ».; .+k,. +` .Ry

^k .i.! ,,.00M°fase^
78

6 ,. s+F ^*^ ti y

TABLE I, OPERATOR COEFFICIENTS FOR STANDARD

OPERATOR COMPACT IMPLICIT SCHEME

q J ` 6 - 5P1 + 2Pj+i - pJ P J+i

q j 64 + I6p)+l - 16pj- 1 - 4P1- 1 p j+i

q j6 + 5Pj "

 2p)-i - PjPj-1

I ^

rj- II q-( _ 3

2 Pj-i) + q c 

pj) + q +) 

(I 4 2 pj+,) + h 2

 qj Cj-i

rJ -

(ry+ + rj-) + h2 (g j`cj-t + q ;cj + gj+Cj+!)

rj+ qj" - z P1 -I) + q )l 1 + 2 Pj) + q j+ (I + Pj+i) +h2g j Cj+l

where

PJ - hb^
ORIGINAL. PAGE 19

OF DOOR QUALITY

TABLE I.1. OPERATOR COEFFICIENTS FOR CE14ERALIZED

OPERATOR COKPACT lMf'L7CIT SCHEME 

q^ - 6+ L p, -3] Rcl +Lpz]Rci2

q j • 60+1 lop j Rc) + l p 3 JRc1 + [T1+Ip,jRcJ

q • 6 f 1 p, +3) Rc1 + [P,+P,.,)Rc?-1 + [P,)Rc 1^ 

where 

pl V 3 1 P2 - 0 V P3 - max [ 7r, 70

• (T +'r „) r W _ - 16 - 2 + (a -0p 3(T +a- 1 + 7r 

J+1 J I p2 + 1 +) pl + 1 2 p2 2 I r 1+1 2 z 

_ 0 a1?0 0

2p -o-Z 2:0 

rr1 8 oI1(10-r1+1- Tj-1) 0-1 <0 Tr2 (2p,-o•Z)2/13 2p1 - crz <0

c P / 3 + (TI -1 c _Tj,l^ 2Q'2 3Tj _1 — Tj+1 + 10 + 2hT1-1

1 

_1

I I f0 

—T1 +1 T1

-I -i

I c) -1

pq [I+Tj+11- 1 7r 3 W3 ` P3 - 7r 1 + W1 + 2Tj

-I C2 + h b -1 ) pZ

with h sufficiently small so that

10bj - bj-1 -- bj+1 > 0 and 2 + he j+1 /b J41 >0 for J- 2, • • • ,J and cj 50

where

bj-I/bj Tj+1 - bJ+1 /bj and Rc j n hbj

r

i . r j . r+ given in TABLE I

79
TABLE III. OPERATOR COEFFICIENTS FOR

ALLEN SOUTHWELL EXPONENTIAL SCHEME

rj` = Rcj a 

-Rcj/ (I - e-Rcj)

rj + = Rcj /(1-e-Rcy)

rjc = — Rc j + c j

qj-=0

qj° = I

qj+= 0

where Rc j

h b j

80
ORIGINAL. PAGE 19

OF POOR QUALITY

TABLE IV. - OPERATOR COEFFICIENTS FOR EL-MISTTYJM WERLE

EXPONENTIAL BOX SCHEME

r)- p-exp(-p-

) / D exp(-P-),

rj+ • p+ / 11 - exp('-r,

rj _ ( r j+ + r,- )

q j

- (I -.rj) /(2p—)

q

J+ • ( rj+– {) /(2p+)

q j q^ + qj

where

P = 2 (Pj-1 +Pj)

P+ (Pi +pj+1)

and

P) ` hb j

81
ORIGINAL PAGE 1.9

OF POOR QUALITY

TABLE V. - OPERATOR COEFFICIENTS FOR SECOND ORDER

FINITE DIFFERENCES WITH ARTIFICIAL VISCOSITY

r S

rj = — 2+ ct

r^ = I + S

q^ - 0

qj

=t

q^ = 0

where S = Rcmax /2 for Rcj > Rcmax

S = — Rc max /2 for Rcj < — Rcmax

t = Rcmax

IRcjl

and Rc j = h bj

f IRc j I < Rcmax; r ^, r j^ r ^ ►

q ^, ^ q^^q

reduce to standard finite differences.

where Rcmax is the maximum allowable cell Reynolds number

82
ua

^M

xM

to

CA

wd

M

^`--,

q

0

t?

3:

a

•

30

0

q

0

oaM

q

^

a (1)

O•

•

q

0

00•

q

•D

00•

CIDO

DD

m •

E • m •

D OO

D

•

co•

4

N co qt

O (D

t!'

M

M

M Cv

£01 X ;^

F

It

x

O_zO

LL

N

a

LLzU

O

Oza

"'

N_m

OD

LL

crauHOza

Lo

ILNYZ

v

It

-+ 

O

CVh

83

`s1
0.16

0.12

y - ft

0.08

0.04

84

1.0 002

0 0.2 0.4 u/um 0.6 0.8
ORIGINAL PAC42 G

OF POOR QUALITY

1,6

1,4

1,2

1.0

0

" 0.8 wv

0,6

0.4

0.2

0

0

DATA — NEWMAN (Ref. 34)

SECOND ORDER FINITE DIFFERENCES

0

\ 01

0 3 4 5

x — ft

^IG.3. NEWMAN AIRFOIL SERIES 2. SKIN FRICTION DISTRIBUTION

AS A FUNCTION OF x

85
M

w CT

CL

W

-OO

SH

U

M

wmw

O

U

UwmWV

M

d'

O

I zw3aza

d) w

O

OQa

I

i^

N

OOam

F--

O_z . w

U<U

wa

I

N

x

V) Oz0

UM

LLO

,a. UOZa

(r

t

fi.

Lq

M

— a) ti to

N

NN———

86

4
ORIGINAL PAGE gg

OF POOR QUALITY

1.0

O

O

O

r

r

/ O

O O DATA — NEWMAN (Ref. 34)

x = 4.493

CE8EC1

O x = 4.51391

/ ---— Eq. 37

O

^p

0.8

0.6

U

ue

0.4

0.2

0 L

0 .2 .3

y — ft

FIG. 5. COMPARISON OF VELOCITY PROFILES

87
J

qI Q

^ ^1

O0 I vX `

(fj

LL

^8aO

U

cr

cr

0

cr

Wx

13

00

O°

I

ol

wU

Q: HOZ

LCU

O

I —1iy

C3,

UWLLO

(nW

O

q

4

J

q

4O

j I - 0a J i

v ' V O II

\ ^Q

OCI

q 41-0

iiIi

^WQNHU

N

z

J .w'7. 3 w

cr

°w

Cz°

O

°

~3

U

a

n

O

13

>

0w

ww0

3w^-

WmH

Q

w

N

J

Q

U.

cfl a2QNWQz

C9

0

s

R

T

0

0

0

to

aM

O

o

O

N

?

31

p

88

9
Nxwaaaw

zb

J

w

J

°w^

3

a

a

o LL Qz wuo 0

0

of

a

>w

N.

ow

3

°w0w

u

o

13

w

zu

3w

umH a w W J Q

u- 6

tia2QNwQz

0

0

0

0

0

0

0

U

-) It

M cv

—

89

^•

cti

s,

fi