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 picture

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 som

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

submission.

Produced by the NASA Center for Aerospace Information (CASI)
NASA CONTRACTOR REPORT 166565 - PZ^' • ?..

^r ^.1

CC

cr

X1111 ±,^^

WSJ, FAoUTY

T

z, ^` ^ES$ pEP. i

Solution of Three.-Dimensional Time-Dependent Viscous Flows

Part 2; Development of the Computer Code

(NASA -CR- 166565 -, Pt- 2) SOLUTTCN OF N84-27002

3-DIMENSIONAL TIME-DEPENDENT VISCOUS .LOWS.

PART 2; DEVELOPMENT OF THE CGMPUTER .CODE

(Scien.t.Liic Researc.h Associates, inc,) 82 ,p Unclas

HC h05/M,F A01 CSCL 20D G3/34 13669

Bernard C. Weinberg

Henry McDonald

CONTRACT NAS2- 10016

December 1980

NNISA
NASA CONTRACTOR REPORT 166565

Solution of Time-Dependent Viscous Flows

Part 2: Development of the Computer Code

Bernard C. Weinberg

Henry McDonald

Scient:'_fic Research Associates, Inc.

Glastonbury, Connecticut

Prepared for

Ames Research Center

under Contract NAS2-10016

NASA

PJational Aeronautics and

Space Administration

Ames Research Center

Moffett Field, California 94035
0

TABLE OF CONTENTS

PAGE

SUtLMARY .

. . . . . . . . . . . . . . . . . . . . . . . . . . 1

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

LISTOF SYMBOLS . . . . . . . . . . . . . . . . . . . . . . . . . 4

ANALYSIS. . . . . . . . . . . . . . . . . . . . . . 6

Background, . . . . . . . . . . . . . . . . . . . 6

Coordinate System . . . . . . . . , . . . . . . . . . . . . 10

Governing Equations . . . . . . . . . . . . . . . . . . 12

Continuity Equation. . . . . . . . . . . . . 12

Momentum Equations . . . . . 14

Energy Equation. . . 19

Equation of State. . . . . . . . . . . . . . . . . 19

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

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

Spatial Difference Approximations . . . . . . . . . . . . . 21

QROperator . , . . . . . . . . . . . . . . . . . . . 21

Linearized Block Implicit Scheme . . . . . . . . . . . . . 26

Implementation of the LBI Scheme Employing

the QR Operator Technique . . . . . . . . • • . . . . . . 32

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

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

NUMERICALRESULTS . . . . . . . . . . . . . . . . . . . . . . . . . 39

Oscillating Flat Plate. . . . . . . . . . . . . . . . . . . 39

Three-Dimensional Flat Plate Boundary Layer . . . . . . . . 41

Two-Dimensional Howarth Flow . . . . . . . . . . . . . . . . . 42

CONCLUSIONS. . . . . . . . . . . . . . . . . . 44

APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . 45

Appendix A - Linearization Technique . . . . . . . . . . . . . 45

Appendix B - Tensor Notation and General Coordinates. . . . . 50

REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . 72

i
SUMMARY

Y

The flow over a helicopter rotor blade in forward flight is an important example

of three-dimensional time-dependent flow. The b!)undary layers on the rotor blade set

loss levels and control retreating blade stall. Ae a consequence there is considerable

interest in developing a numetl,.^al scheme for solving the time-dependent viscous compressible three-dimensional, flow equations to aid in the design of helicopter rotors.

In the present report, the development of a computer code to solve a threedimensional unsteady approximate form of the Navier-Stokes equations employing a

Linearized block Implicit technique in conjunction with a QR operator scheme is

described. Results of calculations of several Cartesian test cases are presented.

These results indicate that the computer code can be applied to more complex flow

fields such as these encountered on rotating airfoils.

1
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 losses 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, some progress has been made in two space dimensions but little to date

` R 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 problem

for steady three-dimensional be ,ndary layers in a very illuminating fashion. This

particular problem arises when the spanwise component of velo/;ity changes sign and will

be discussed in detail subsequently. Because of the interest by external aerodynamicists

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 coordinates

in generalized tensor form. The rationale for the choice of this approach is discussed

ir, detail in fief. 3.

2
The basic a ,)sumpt;ions 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. A novel method is

i

employed for solving the continuity equation r.n conjunction with the reduced

Navier-Stokes equations. The continuity equation is split by employing the DouglasGunn procedure to obtain a consistent approximation to the full equation which is

then solved as an integral. Results of computations on model problems in Cartesian

coordinates are presented.

3

P,
a, b, c, d

n

aij

ds

ei

fi , f (x)

giJ ! 9 i

h

hl' h2

J

L(u)

21

Li

Din

P

Q

C +

q s q , q

4.

y

R

C +

r , r , r

Re

R

c^

S

t

T

T

0

LIST Or SYMBOLS

coefficients of differential equation

coefficient of differential operator

nonlinear three-dimensional differentia], operator

differential arc lexigth.

covariant basis vector

right-hand side of differential equation

components of metric tensor

step size

scale factors

Jacobian

linear differential operator

linear three-dimensional differential operator

product operator Q i i

difference between 'Z n and `P n

pressure

tridiagonal difference operator

components of Q

velocity vector

tridiagonal difference operator

components of R

Reynolds number

cell Reynolds number

source term

time

static temperature

stagnation temperature

4

A.. t',r Y ^'.o/^9x ^YJ .. .,^iL.^I^^: t^/^ '^ ^ F. 4,IT{C ^1LL. ^Rn .p^ . 
5

V7

V

u

U, V, w

x

Greek Symbols

0

Y

r k

A

At

Ax

C 

jk

P

ik

a

W

W

Subscript

e

p

CO

Superscripts

n

*, *V?, ***

contravariant velocity component

physical velocity component

coordinate direction in i direction

factor for centering time step

ratio of specific heats

Christoffel symbol

Kronecker delta function

velocity divergence

time step

step size in i direction

strain tensor

boundary layer coordinate

viscous stability parameter

molecular viscosity

density

stress tensor

vector of unknowns in LBI scheme

source term in LBI scheme

frequency of oscillation

CFL stability parameter

edge

physical

free stream

time level

intermediate time level in -ADI procedure
5

ANALYSIS

Background

In this section the requirements of a three-dimensional unsteady viscous flow

computer code for flow over airfoils are discussed.

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 time-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 purpose:

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 (Ref. 4),

that very siGnificant transient boundary layer effects can be observed.

In examining the flow problems of practical 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 three-dimensional.

In rotary wir g 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. As a consequence it might be supposed that the

usual three-dimensional thin boundary sheet approximations (Nash and Patel, Ref. 5)

could be used to produce a valid set of governing equations. Fortl^nately some improvements in thin boundary sheet approximations are possible as a result of having to

eliminate the negative cross flow problem mentioned earlier.

The negative cross flow problem is best explained in a somewhat intuitive manner,

and a good physical description of the problem is given by Kendall, et al (Ref. 1).

Looking at the suction surface of a conventional swept back wing the boundary layer

cross flow, w, is usually outward in the z positive direction along the span from

root to tip. Thus 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 view of the

physics of the problem, the spanwise marching scheme does not normally encounter

6
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 twodimensional boundary layer equations into a region of separated flow and its obvious

relationship to the negative cross flow problem, it is not surprising that spanwise

marching into a negative cross flow region is not accomplished without special treatment.

Recently conventional boundary layer developers have been turning to performing an

implicit spanwise construction to remove the restriction of only positive cross flows

(Kendall, et al, Ref. 1). Lin and Rubin (Ref. 6) in their predictor-corrector boundary

region solutions for flow over a yawed cone at moderate incidence also show that allowing diffusion in the spanwise direction not only eliminates the problems associated

with negative cross flow, but improves upon the solutions obtained by three-dimensional

boundary layer techniques.

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 the rotary wing applications and at least desirable for

fined wing applications, especially as it can be had for a very modest increase in

code computational labor. Based on the experience in Refs. 7 and 8, the spanwise

implicit sweep would only result in about a 20% increase relative to the explicit

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

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 timedependent flows a similar structure is to be had so that at each time level one

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

in the spanwine and normal directions to give the desired unconditional lability.

► 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 some cost as the explicit sweep does arise. If

one does perform a streamwise explicit sweep, then the linearization of nnlinear

terms is performed about the known spatial marching level. If an implicit streamwise

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, it is easier to obtain a consistent spatial-temporal order

accurate linearization by marching in time than in space (in time the nonlinear

marching derivatives have the form u 4

 whereas in space marching they have the form

u

i

 U). Further by structuring implicitly in the space marching direction, (small)

regions of axial reverse flow would be permitted. As a result of these combined

benefits of linearization and separation, -, stroarr q ise implicit structure is advocated

and has been implemented in this effort.

Transient calculations mean that, in essence, a full 3-D spatia3 integration is

carried out at each time step. Thus, spatial accuracy is very important to minimize

the spatial grid point density for efficiency since many time steps are contemplated

in a given cycle. In order to get the most out of a given spatial difference formula,

the errors from representing nonlinear terms by linear combinations of terms should be

less than or equal to the spatial discretizaLion 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 riot recommended since this only reduces

the linearization error and computationally costs as much as a complete Zime step.

Cutting back the time step would be preferable 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 . U n+i W n + Un

 W 

no

- U 

n

W n + 0( At?)

8
and clearly one cannot lag wn+l at the old time level n without introducing a first

order time o:?or in order to get an uncoupled system, Lei, wn*l not appearing in the

streamwise momentum equation. Thus formal linearization and consideration of the

resulting errors indicate the coupled system ought to be treated from the accuracy

point of view, ^nis 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 approximatioii arises unconnected with linearization

but arising from basic coupling terms in the original equations and if indeed some

terms 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 (Ref, 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 difference molecule

accuracy. The linearization technique is described in Ref. 8, together with its application to block coupled splitting schemes. Schemes of this ge-neral type are here termed

"split linearized block implicit" or split- LBT 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 CrankNicolson centered time implicit scheme for instance, although second order in time,

shows quite a dispersion 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 dependenc y is continuously

input through initial and boundary conditions and relatively the concern is with "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 of a wave in a shear after a long propagation time. The

interest is in forced 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.

9 
The governing equations that are considered here are the Navier-Stokes equations,

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

form for a bo'ay oriented coordinate system (boundary layer coordinates), In 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 herd only for 

purposes of simplification, In the analysis that follows, the full energy equation

could equally well have been usad. Employing the equation of state which relates the

pressure p to the velocity components .0 and w by an algebraic egwntion, 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 oi: .::Ra.s effort is to solve for the flow over airfoils an understanding of tb, tyx,^ ;f geometries to be considered is essential to guide the choice of

the coordinate s;stwit. and the structure of the computer code.

Consider a typical finite span swept wing airfoil as shown in Fig. 1. 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-Stakes 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 while the other coordinate directions must lie on the body surface. Furthermore,

the coordinate lines normal to the surface are straight. 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.

10

- ^ !,. .' +. ,,'6+ . s .• ?^....^ ...rat> i..rr..iru!^I+---^ .t!^as .^, .... w.. __^^_:
A

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

gained by the Judicious choice of coordinates. The most obvious advantage is that the

physical boundaries of m flow region can uo. represented by coordinate surfaces. This

removes the need for fractional cellar in general; hence, the complications and loss of

accuracy associated with a boundary 'interpolation are removed. Another advantage is

that a uniform numerical method can be used. The solution can teen be performed with

a fixed number of cells in any given direction and with a uniform mesh spa,'ing.

In Fig. 2 we can see the advantages of a nonorthogonal grid which conforms with

the boundaries of the swept wing and covers the entire airfoil over a Cartesian grid

whi , does not. In addition the coordinate transformation can be constructed to contain

distributions for pnysic:al space mesh points. In this context, the uniform mesh of

computational space is simply mapped into a suitably distributed mesh in physical space.

The resolution of large solution gradients is the major objective in the selection of a

coordinate mesh distribution, as in the resolution of an attached boundary layer,

Another ;:,,ore subtle example is the resolution oL large gradients in computational coordinates due to regions of high curvature on the bounding surfaces. When the transfor=tion contains the mesh point distribution there is no need to construct the apparatus for the discrete approximation of derivatives on a nonuniform mesh. This results in

a savings in both computer logic and storage.

Therefore, 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, in order to suitably distribute grid points in regions of

large gradients, provisions are made for analytical grid transformations (Ref. 11) in

each coordinate direction.

In view of the type of geome,ries to be considered and the assumptions made to

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

coordinate system is advocated where the metric tensor which has four independent

components is given by

gU 912 0

9 ij 912 922 0

0 0 933

11
Tile subscripts 1 and 2 refer to the directions oil surface of the body while subscript 3 refers to the direction normal to the body. Furthermore, the metric data in

the coordinate directions on the airfoil surface do not vary with tine normal direction,

i.e, the metric data in a 1 - 2 surface above the body are evaluated on the

the body surface (Ref. 10). Since we will be dealing with nonorthogonal coordinates

it is advantageous to derive the equations in general nonorthogonal coordinates employing generalized tensors. In Appendix B a brief description of tensor notation is given.

Further details can be found in Refs. 12 and 13.

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

—ordinate dependent). The grid transformations are absorbed into the geometrical

coefficients, leaving the equations unaltered in form. This point has a considerable

effect on the development of the computer code. Since the only geometric information

that must be input is the definition of the metric data and their derivatives, it can

be contained in one subroutine without modifying the remainder of the code.

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. We will show that this notation aids in the ordering of the various

terms in the equations and in many respects simplifies the construction of the computer

code.

In the following derivation the governing equations are nondimensionalized as

fos.".ows, xi with respect to the characteristic length L, the velocity with respect to

U., density, pressure and temperature with respect to p . , p.Um2 and U. 2/c

p 

respectively

and time with respect to L/Uw . Viscosity is nondimensionalized with respect to uo.

Continuity Equation

Consider the continuity equation written in vector form so that it is independent

of coordinate system i.e.,

a

ap +o - pq =0

12

(1)

ti
Joe

13

where p is the density and q is the velocity vector. L,:pressing the velocity vector

in a covariant basis

q = u F^ (2)

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

vector in the x i direction. The velocity vector may be expressed in a number of different forms, each with certain attributes. Here for the moment the velocity vector is

expressed in a covariant basis, for simplicity. Later the velocity components will be

transformed into physical components for numerical solution. The contravariant basis

exhibits variation in its components for instance in slug flow, 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 se.

As such, it is felt that the actual computations are better performed on the physical

components.

The divergence of a vector (cf Appendix B) is given by

0 ' P q = pukt^ _ (puk ), 9 + pub rk

= {Jpuk),k (3)

where pukIk is the covariant derivative, pu k , k is the partial derivative in the xl`

direction, J is the Jacobian and P ik is the Christoffel symbol (cf Appendix B).

In Bquatiun 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 jmoothness of the geometrical

quantities. However, herein we use either form solely from the point of convenience.

For the continuity equation we use the form involving the Jacobian while in the momentum

equations the form involving the Christoffel symbols is employed for the evaluation of

the explicit (lagged) diffusion terms. Thus the form of the continuity equation which

is used can be expressed as

at + J (Jnu k ), k = 0 (4)
OF POOR QUALITY IVIk

Momentum Equations

The momentum equations in vector form can be written as

Pq = P[ dQ (V.=

V

' a, (5)

where v is the stress Censor.

In generalized tensor notation Equation 5 becomes

P (

t

uk ul I I e l _

(,

a-iklkei (6) a t

The stress tensor is defined as

Ik = rr 2 μ' 1

* ik li

a- - - [

P+ ^J' +

elk

3 Re Re (7)

where )i is the viscosity, p is the pressure, Q is the velocity divergence attd c ik is

the strain tensor. The Reynolds number, Re, is defined as p.V,,'L/u.. The strain

tensor is defined as

Elk = U lmg mk + u klm Sim

($)

where gmk and gmi are the components of the metric tensor. Employing the fact that 6 k = 9 k and substituting the definition of the strain into the stress tensor

we obtain

Pik = _ 

\P 

+ 3 R , gik+ R ui ^ lm *mk R u k Imgmi

(9)

Substituting Equation 9 into the momentum equation and employing the relationship

9kiI = 0 (10)

k

14

a T-_•`ms s\^^ ^ -a^!FP ^rrp

' _..:aro.*r o^.^. .^^ _ ,.
we obtain for the i-th momentum equation in the ei direction

P[ dti + ukUflk J = _g1k(p + 3 Re &))k

(11)

+ 

g' L R U' LI I k + 9m'1 R uk lmI I 

In Ref. 3 it was pointi,d out that the QR Operator scheme requires that derivatives

in any direction operate on only one variable. In the momentum equation this requirewent 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. However, the use of

the QR Operator scheme requires these terms be treated explicitly. This together with

the quasi-linear form of the governing equations are the major limitations that arises

in the treatment of the approximate form of the Navier-Stokes equations considered from

the use of the QR Operator scheme. 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 quasi-linear and, for instance; the full conservative form of the differential equations, are equivalent. In discrete form, various formulations of the 

governing equations exhibit different properties (Ref.29). 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 auxiliary equation relating the pressure gradient term to the derivatives

of the velocity components and would increase the block size of the system. An 

assessment to 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, we will first split the governing equations into

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

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

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

15
(12)

Diffusion Terms

Consider the term

(μUml,)I,

ORIGINAL PAGE 19

OF POOR QUALITY

Since mi:ad partial derivatives are commonly treated explicitly in orthogonal

coordinate systems, we will do likewise in generalized nonorthogonal coordinates and

extend this concept 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 construction of thee computer code,

and is thus adopted here. Furthermore, by not splitting up the covariant derivative

for the purpose of making ih implicit, but rather retaining it as a unit could prevent

instabilities that may arise due to time splitting. This occurs when two portions cf

one term should cancel identically but cannot due to their being split between two sweeps.

If j = k the term is retained as implicit, and if j ^ k then it is lagged. We will

consider the case j ^ k first. Upon expanding the explicit part of the diffusion

term it becomes

( /J. U1

- ^ μ(U! j+ 

Un rijn)1,k + μL U m tj 

+Unr Jln]rim.,

(13)

-μ[Ul , m + UnrImn] r;

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)

i i )

L

n ( + U (μu ri. UI T^ ^ + )U.Si I (U + (μu ), j =μ

rjn μ n ^j in -N ^n' jj N 1nl

n+ T i (14)

Jj

where

Sijni = 

rmjn rl mj - rImn rm jj +rl jn,j (15)

16
t

ORIGINAL. PAGE 18

and OF POOR QUALITY

Tljj ° 2f` l

U ^'j rMi un', j

r mj^ + f ► j^ um ^mj +unrjln^ (lb)

m--i+l, n=5-2i and no sum on m and n

Since T velocity components and derivatives in directions other 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

quasi-linear form as

E l CC 1

1 9 jj /, lu l , jj + IC1 9 (^u,j+21ir lji) — Z Clkgkkμr kk]u^,J

j =1 k=1

+ E Clk9k ilS1kj^ u j + ^^i

 [ Z ^I gkk^^k rk, 1 u(1 (17)

k = 1 k=1 l J

+ Z Tj^ Cljgjj + z f z 9nk[ μ ul l

n 

]I k + Z 

gn1 

LN' uk I n ] IkJ ):I krl n=1 n=I

On n^i

0 0

where s1j

I 0 I 0

10 0 1 (18)

and 2 1 1

CI j = I 2 I

(19)

1 1 2

and repeated indicies do not indicate summation. The last two summations can be

combined into one explicit term, so that the diffusion term for the i-th momentum

equation becomes

17

,p
ORIGINAL. PAGE 13

0P POOR QUALITY

jZ l ^ ^I g)) μ u'il7+ [

C ) 9 )) ( μ0 + 2,.F' ) 

- k ' C^ g kk^r Tk k ] u 1 ^)

+ a1 ^ikgkkμ^k Tkii ] ui +

ck gkk^ Sik)k] 

uJ

k-1 k`•I (20)

+ l e i gkkThk + Z 9nk[ N ui ^n^ Ik+

gni 

C ,' uk In ^ I k k=1 n=1 n=1

kin 0i

and there is no summation for repeated indices. Note that the diffusion operator

terms for a given direction have been cast into standard quasi-linear form i.e.,

au xx + bu x + cu + d

Convective Terms

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

PUj ui = Pu j Iu i , ) + uminIji (21)

which equals when expanded

i 3 J i

2 1 m i 1 3 n J i

^" u,u I1 ^^ { p u u+m^ pu u Pmj + b3 Z1u

^^ } pu (22)

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

by substituting Equations (20) and (22) into Equation (11) and treating the pressure

gradient and velocity devergence 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

".., a. 9 ys,'. w.^ +r + -.r.^^ r".. _ANR7 ^sAsor.
ORIGINAL PAGE E9

OF POOR QUALITY

Energy Equation

The energy equation employed here states that the stagnation temperature is 

constant throughout

To-- T+ 2 q2 (23)

The generalized tensor notation q2 is given by

q2 = U l u j g .

I^

where u i ane u are the contravariant velocity components. Incorporating the assumptions made concerning the coordinate system we employ, i.e.

=p

g13 g23 

we obtain

q2 = ^U1^2 911 + 2u U2912 (u2) 2922+(U3)2 933

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

physical velocity components, i.e.

Up =U 1 h 1 WP = U2h2

we obtain

To = T + 2 ( up2 + W 2)) + h

12 

u pW p (24)

1 2

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-1

Y PT

(25)

19
ORIG WAL PAGE 19

OF POOR QUALITY

Linearizations

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

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

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

In order to overcome this difficulty we employ the linearization procedure (described

in Ref, 8 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

(p.e)n+9 = (pngq)en+A+ (pnen)4n+A

+

 p T

n nen

n

( (u' )n + 9 1 2 (u 2 ) n I (ui)n +a

L

h 

1 2

+ pn Ten n r(u2)n+ 

hl n

(UI)n1(u2)n+p (26)

l iz 11

pn^nen (l n+/3

+ n L 2 ( T n - •T o ) — 7 2L / f ^ 2pn nen

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

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

a velocity component. In the case of a term containing pu i we set ^n_,n+S = 1.

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

velocity components are used. However, as noted in Ref. 14 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.

20
The Turbulence Model

We treat the set of three-dimensional ensemble averaged turbulent reduced NavierStokes equations. Ensemble averaging permits the appearance of low frequency (relative to the turbulence) time dependent "mean" flow. It is, therefore, necessary to

specify a turbulence model suitable for this problem.

The approach taken in the present effort assumes an isotropic Lurbulent viscosity,

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

Reynolds stress = 

cr Ray Re L EIk 3 s^^ ^J (27)

Using Favre averaging (Ref. 15$ 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, )1, and the turbulent

viscosity, uT.

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 will be considered. The very versatile QR Operator notation will be introduced, whifrh allows

as special cases a variety of schemes such as standard second order finite differences, first order upwind differences, fourth order operator compact implicit

(OCI), fourth order generalized OCI and exponential type methods. Since all these

schemes are of the same form, a single subroutine which defines the difference

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

structure of the program unaltered. Subsequently, 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 detail in Ref. 3.

The QR formulation alloo•,, 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 mesh formulations.

21

-_s-

*
Q

The general concepts and notation for two poilit boundary value problems will be

introduced and then the methodology entended to more general linear and nonlinear

parabolic partial differential equations in one dimension. The extension to multidimensional problems will also be indicated.

Consider the two point boundary value problem

MU) - a(x)u XX + b(x)u X + C,Wu r i(x) (28)

with u(0) and u(l) prescribed. Derivative boundary conditions, although not described

here, can easily be incorporated into the framework of the Q-R operation notation.

Let the domain be discretized so that xi . (J-1)h, j - 1, 2,..., J + 1, and U 

u(x1 ), F ax (x1 ) S^ "j uxx (xi ) and h n l/J is the mesh width, The numbering convention was chosen here to be compatible with FORTRAN coding.

Without loss in generality for a(x) y0 , Eq. (29) can be divided by a(x) so that

we may treat instead the following equation

E(u) - u xx + b(_x)ux + c(x)u f(x) (29)

where

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

The spatial differential operator is identified as

L(u) - uXX + b( x)u X + c(x)u

(30)

Substituting the finite difference approximations to the first and second

derivatives

Do

2h U - U + 2hj FU' ) - u X ( x ) ) + 0( h2 ) (31)

D+D_ U 

= 

u j _ 1 -2 2j +Uj+I e S

uXX(xj) + 0( h2) (32)

h 1 h j

t

22
ORIGINAL FACE 19

OF POOR QUALITY

into Eq. (29) And rearranging, we obtain

J 1J Ch 2h ^ J [ J h^ IT 2h ^ J )

or

C

R c (' Rc

I - Uj.^ + Lhz c j • 21UJ + I + 2 ]UJ+^ h2 j (33)

where Rc bb^ is the cell Reynolds number.

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

rj U. I + rc U + r* -up - hz(q" f - + `f + ' f )

J .^ 1 J J j 1 J q1 1 q J J+1 (34>

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

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

respectively, and where the r i 's and qj 's for grid point 3 are functions of h, bJ_1,

bi, bj+1, cj _1 , cj and cj+1• Comparing Eqs. (33) and (34) we can identify the r3 's and

qj 's, viz.,

rj (- Rcj /2 qj = 0

rj . hzcJ _ 2 qj = I

(35) 

rj I + Rc J /2 q^ n 0

We now define the tridiagonal difference operators Q and R

` RLUjl rJ Uj-1 + r U j + ry Uj+1

- Q[ fj] ` q J fj - 1 + qj fj + q) fj+l (36

23 
^i

C

Noting that L(u) w f and substituting Eq, (36) into Eq. (13) we obtain

R L Uj j - h?Q [ L(u) j

` h2Q

[

t 1] (37)

Alternatively by employing the inverse operator Q- an expression for L(u),

J 

call be

obtained

L(U)j R V Q" RUj

For standard central finite d3.fferences Q - Q71 - I, the identity matrix,

(the spatial operator is given explicitly in terms of UJ„10 U

i 

and U

J1

) so that

nothing was gained in obtaining Eq. (38). However, in general, for highter order

methods Q is tridiagonal and Q

- 1 

is a full matrix, Hence Eq. (38) gives us a means

of expressing the spatial operator for a wider class of difference approximations.

The formalism in Eq. (38) is also applicable for first and second derivatives appearing alone (ci. Ref. 28). It should be pointed out, however, that Eq. (38) is not tdie

most general formulation since the compact- implicit formulas cannot be combined to

yield a single scalar equation rebating the spatial operator to the function values

(Ref. 28).

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

fourth order tridiagonal methods which possess a monotonicity property as the cell

Reynolds number is increased, R c-* O', We will not repeat it here, However, the resulting Q and R coeffio,tents are given in Table 11,

Another family of schemes that can be expressed in Q-R operator notation are

the so-called exponential methods. The idea, original..Ly due to . Allen (Ref. 17)

(independently derived by I1'in (Ref. 18) and McDonald (Ref. 19)) and employed by

Dennis (Ref. 20), 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 Q and R coefficients 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 i - where

the scheme reverts to first order upwind differenci.ng.

Another exponential scheme which is ujj;iformly second order accurate was developed

by El-Mistikawy and Werle (Refs. 24 and 25). 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. (29). Berger, et al (Ref. 23)

derived the counterpart for an operator of the form given in Eq. (30), but with c = 4),

The Q and R coefficients are presented in Table IV. Although this scheme reverts to

24
Y

OF POOR QUALITY

second orfer upwind differences as Rc - r­, it does not possess a maximum principle

analogous to the ordinary differential equation it is approximating as does the

exponential scheme of Allen (Ref. 17). In Table V a second order c n tral difference

scheme is presented which adds artificial viscosity tothe spatial operator when

jRcj>2 so that jRcj never exceeds 2. This scheme is employed in the solution of

several model problems and will be d;tscussed in greater detail in the section on

numet.;ical results.

Application to Gougled Nonlinear Parabolic Equations

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

on the Q-R operator scheme in solving a system of nonlinear parabolic equations.

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

T

(UI `] -U ) n+Q

^T o n+p ^ ^t

NI (U^,U2, .... U m , X 1 x 2 ,X 3 ,t) 0

Ij

j•112, ... 0J+I

where Nn

+$

 is a quasilinear spatial operator, the Q-R formalism carries

directly over provided that for any equation only one independent variable is operated

upon by the differential operator. For example,

' a(u,w,v) U t * UXX + b(U,V,W)UX + G(U,V,W)

is allowed since x derivatives of u only appear, while

a(u,w,V) U t ` U9X + 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 OCI schemes are applicable. Note that within the

25

f

'..,•_Oil
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 tRken and for 

instance the Douglas-Gunn formalism is achieved too, no particular problem arises

other than the cost of an additional implicit sweep is incurred,

aL 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.

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

A ^ f = . + (39)

where ' is a vector of unknowns and I is a source term vector which is a function

of x l , x2 , x3 and t. Extension to source terms which are functions of 4) are discussed in Ref. (8). „b is a three-dimensional nonlinear differential operator and `

of the matrix A acting on the momentum equations is equal to pI where p is the

density and I the unity matrix.

Equation (39) may be entered about n+p time level, i.e. t

o+$ (n+S)At =

nAt+OAt = t o+CAt, and written

(40) An+Q C cl) _ ^

n ^ ^^t n+^^n+^+ Tn+^

0 < 3 1 is a parameter allowing one to center the time step, i.e., S 0 corresponds

to a forward difference, B = 1/2 to Crank-Nicolson and B = 1 to a backward difference.

Equation (40) can be linearized by Taylor series expansion in time about the n th 

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

An rn+^-

n

At

 = 

1 n [;5 n +Q ^n ] _ ^nc^ n +

n +/3 (41) 

where , is the linearized differential operator obtained from _2) by Taylor Series 

expansion in time.

26
ORIGINAL PAGE 19

OF POOR QUALITY

The difference between the nonlinear operator .!b and the linear operator is

definedas M

n `^ n _ n. The intermediate level n + S is defined as

n +/3 = Q ;51 +1 + — Q^ ;$n

(42)

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

hybrid implicit explicit scheme is obtained

n n +l n

n n +l n In n n n n +Q

(43)

The vector n*O 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 

^n with a consequent order reduction in

temporal accuracy.

The operator .1 is now expressed as sum of convenient, easily invertible sub- 

operators.1 + .1+ ,^ In the usual ADI framework these suboperators l 2

M.

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

suboperators are expressed in the QR notation introduced earlier. Writing ,fn+s and

Mn(Pn as a single source term S n+S the algorithm is written

An[ (Dn+l _4)

n ]

/At = ^3[1 1n f .1 +.I n ,[^n+1 ^^n] +[1; +12+13] n + Sn (44) 

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. 28)

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

An [ 4)*-(Dn ] /At = /31i ((D"'-4)n ) + [1"+1z+I n cDn + Sn +Q

(45)

An [(,*._ (,nl ^,t 

= Q1 [cl'* _ CD n] + Q1 z [ ^*

»_ n + [1 1n +1 n +1 3] CD +S n +JB

An [,,**r 4)nl /fit = ^1;[^* _^nl +Q12 [^*w - n1 + R13 [^M** ^n1

+ [1n +.1 +,03 ]^n + Sn+/3

27

114:
OF POOR QUALITY

which can be written in the alternative form

An _ At /31' n] [ 4) 4)n = At C1°# n +,p2 +1 n n + A t S n+p

[An 

^t p12 ] [e*_ (,n = An 

C 

e_ (,)n] (46)

CAn- At p13 ] C^**x_n^ = An C^ — (pn^J 

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

factored form

(An - RAt.Zen)An-(An-/3At^2n)An-'(A^/3At'P3n)( n+1 -fin ) -

(47) At (

n + 1 ^" +1 3n) 

(,n+ 

A t S n+a

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

2 and 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 it follows that use of the QR formultion results in the appearance of

the inverse operator Q-1 with the sub-blocks of the -.^ 1, -P , . 0 3 operators. In

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

the scalar operator Q is generalized to the vector operator Qi with (diagonal) subr

blocks Qji . In this generalization j = 1,2 apply on the momentum equations and j = 3

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

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

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

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

product operator Q .P i as Li, the inverse operators are thus removed and the scheme

written, once again dropping the vector superscript for convenience

[Q i An Ot(3Lin]C4)'i- 4)n AtL i (pn + AtQ, [12+13]4)n+OtQISn+Q.

l

(48)

C Q,An — At (3 L n][4)** J C ^iFlk_ ^j n^ _ Q? An [^ * _ ^nl

CQ3 An - At (3L a 

CCD Q3 A

n r - C15n^

12 ,) n 

= Q2 -1 R2 ) n .13n4)n = Q3 -1 R3 cDh

cl) n+i = (D* * *-+ of Ot 3 )

28
With the removal of the inverse operator Q -1 

, the question of the proper intermediate level solution boundary conditions can be addressed. As is rointed out by

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

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

condition operator Bi after linearizing the appropriate physical boundary condition by

Taylor series expansion in time as

BIn «n+1 ) = 9 ( t , 4'n)

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

Ban Q3 An[ (p** -

4

)n ] = [ B3 03 An _ Lit p,3 L3n ] [ 4)* **_ cDn1

BZ 02 

An[ V - CV]_ [ BQ p2An — Lit P13 L2 J L (D

** -

(Pn] (49)

and note that unless BiLi com;nute (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 B nLn . For the present the term At$BiLn ,, - Vi n] is neglected. This introduces

an error of order 0 [At (4, -4+ n )] into the solution but note that this error disappears

at steady state where 0*** _ 0** _ 0*. Neglect of the atsBiLi ,n] term is of course

equivalent to applying the physical boundary conditions on the intermediate level

variables.

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 ` ,1 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. 28), 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.

An inspection of Eq t (48) reveals that only the linearized operators Li, L2 and

L appear. Indeed, the computer code employs this feature by evaluating these three

operators before the first sweep, storin g them and accessing them as needed in the

29

____

t317.""_"^
subsequent three sweeps. In addition, the terms arising from the nonlinear terms are

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

the terms in the differential equation

The operators n ll n ,^ n

p 1^1, .z, 3 can be represented in standard form at each grid point,

i.e.

ln^ 4) n 

° an `^'i,ii + a^2`^'i,i + an (^ + 

a n

4`p + a n^a (50)

In Eq. (50) the subscripts of V indicate the velocity component (associated with the

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

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

the equation is in quasi-linear form, a necessity of the QR operator technique employed

here since the coefficients of the derivative operators need to be identified. Alternate schemes have been proposed by Leventhal (Ref. 33) for equations in conservation

form but are not considered here. In the following section we will describe how this

entire operator is discretized by employing the QR operator format, and show how

the discretization is incorporated into the LBI framework in order to solve the system of

equations (48).

We consider first, however, the continuity equation. Since it is a first order

partial differential equation it does not have the standard form of Eq. (49). Furthermore, p has been linearized and 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.

As pointed out by McDonald and Briley (8) skillful partioning of the resulting

matrix can lead to significant decreases in computation time. An inspection of the

system of equations under consideration reveals that substantial savings can be

realized if the equations are partioned appropriately. Due to the use of a boundary

.layer coordinate system, in the two momentum equations the normal velocity appears

only in conjunction with terms associated with the normal "3" direction. Hence, in

the first two sweeps, for the streamwise and spanwise momentum equation one is required

to solve only for the two corresponding velocity components 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 two 2 x 2 block matrices and one 3 x 3 block matrix

rather than three 3 x 3 block matrices which accounts for the reduction in computation

time.

If we were to consider the full Navier-Stokes equations which include a normal

momentum equation the aforementioned partioning could not be applied since the

normal velocity would appear in all three sweeps.
ORIGINAL PAGE 13

OF POOR QUALITY

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

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

gives us the answer. The continuity equation written in conservation form becomes,

at + J a 4 d P ul , ° xT (51)

linearizing p we obtain in increment form

A" AU"+t + Bn A Wn+t + 'AJQ aX 3 w"A" Aun+t + 

v n B Awn+t + P" Av"+il

At a [ dPU I ]n + Q ^ ( (Pn + un e) Aun+t + (U"en)^Wn+t] (S2)

+ At,8 a2 [ (Pn + wy^Awn+t+(Wne)AUn+tl

i ax J

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

and v _ u3 . J is the Jacobian and

n Pn ( n n]

A - Tn 9 11 + 912 W J

B - n Pn Tn I 9 22Wn + Un 9t2

Approximating equation (52) by employing the Douglas-Gunn procedure as a third

sweep equation, we obtain a consistent approximation to the continuity equation, i.e.,

the x derivative term is evaluated at the * level and the x2 derivative 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. We can thus write the

equation in symbolic form

31
OF POOR QUALITY

n Oto d n +i n n+111 1

A Dun+ ^+ B

n ^wn+ ^+ J 

dx3 

C J {A v n Qun+^ + v

n Bn ^w n + p Qv I J

Sn _ Al 

._a i_ [ i J 1 f * -011, 

a 2 

J dx x L J l (53)

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

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

J3 

L An A u n+l + Bn AWn+

^

dx3 + pt At (vnAn Dun+i + vn BnAWn+, + pnAv +j J

r^ I Sn - /3^t ^3^t r * }

dx 3

(54)

x J L 

In the nc:t section we describe how this is done very easily via the Q-R operator

scheme. The concept of integrating directly the continuity equation is not new.

Davis (Ref. 34) in his coupled procedure for the solution of two-dimensional steady

boundary layer equation used a trapezoidal rule to integrate the continuity equation.

Weinberg (Ref. 35) (Ref. 36) 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 the2 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 Eq. (47) in which we solve both momentum equations

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

[ 

n - 

,,,n3 ] ^ *** ' 

AnOCD (55)

32
OF POOR QUALITY

where A(P*** is the column vector of unknowns, u ) v, w. Here we have implicitly

assumed that the equtions have been appropriately normalized and that the contravariant velocity components have been suitably transformed into their physical components. Employing physical components, (cf. Ref. 14) leads to a better behaved solution

since these components are not unduly influenced by geometrical variations.

For the streamwise momentum equation we obtain

.13n A U133+ U 23 AU 13 + 0

33 A u + a43 LAW + a 53AV (55a)

while for the spanwise momentum equation we obtain

I

3 .^ik*Ik= AW, 33+ b23 '6^ 3 + b330w + b430V + b53Au (55b)

where we have omitted superscript *** from Au, Av and Aw. Now $.n equation 55a, we can

approximate

Q I^ RI 1 38 + a23AU13 + a 33AU by

AX 

2 dU

3 (56)

the operator equivalent, so that

13n A V**= AX 2 + a43AW + a 53 AV

Q - ^R Liu

3

(57)

Similar approximations are made for Equation (55b). After substituting Equation (56)

into Equation (54) and multiplying thru by Q we obtain for the streamwise momentum

equation

I Q i Pn -

Q X R 1 1 Du - j8At Qi a43Aw - j3At0ja53 AV = Q 1 PnAu** (58)

where A = At/Ax 32

33
34

Similarly for the spanwise momentum equation we obtain

CQ.pn -- PXR2^ Ow , PAtQP43AV — AAtOP53Qu = Q2 p^^w * (59)

We employ the same type of procedure for the continuity equation. Since the continuity

equation involves only first derivatives, so that they can be represented as

aQcRc

d 3 = QX3 (60)

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

qr =0 , 4^= 2 , q^ _

2

r

c 

= 0 , r^ = -1 , rc+= I

Simpson's rule

qc

I c

3 q c r

4 +_ I

'

q c - 3

r^ _ - I , r^ = 0 r^

The continuity equation thus becomes

JAn ^u +JBn ^w+ QXa Q c ^ R jjeV"Au+ JBn vnAW + J^nAV, = RNS

3

(61)
CRICINAL PAGE 19

OF POOR QUALITY

where R11S 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/Ax3

we obtain

IQcJAn •+ PWRcJAn vn lAU + LQcJBn + PwRcJBn v n IAW + [ PwRcdpn]AV = Q jRHS) 

(62)

The resulting matrix derived from Equations 57, 58 and 61 becomes a block 3 tridiagonal

matrix (Q and R are tridiagonal operators) with each sub block taking ors the form

r

ii

Qi Pn —/3XRj J '^iAtQ,a,, [-,BAt0,a5, Au

it

Y

er

.4

ai

., l r `)9At Q 2 b 53 ] [ Q2Pn -(3XR2 

J 

C j3AtQ2b43 , 6 

ti

^j

k' [QcJe +,3w Rc JAnV n] CQc J6n +PwRcJ01 C13wRc Jp J AV

x

Qi(ALi**)

= Q2(AW**)

Qc(RHS)

This matrix is inverted by standard LU decomposition.

35
ti

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 no

slip is prescribed for all the velocity components. At the outer edge of the viscous

layer the values of streamwise and spanwise 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 profiles 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. Since all the multidimensional problems considered here to da; =_ have

been steady, their has not yet been any need for a high temporal accuracy solution,

and 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. (The type of schemes

used to date in the present study). Fourth order methods have not yet been applied

to any of the test problems in this report.

The question of proper intermediate boundary conditions for fourth order methods

until recently has not been resolved. Fairweather and Mitchell (Ref. 37) 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 sccuracy for their method. As pointed out

by Fairweather and Mitchell (Ref. 37) their scheme is inconsistent. It is this inconsistency that requires one to use appropriately derived intermediate boundary conditions in order 1.o recover a steady state solution independent of time. However, if a consi

tent scheme were to be used, e.g. Douglas-Guni, then physical boundary conditions can

be applied without any loss in steady state accuracy. A more complete discussion of the

implementation of ADl schemes with the appropriate intermediate boundary conditions

for QR operator schemes (including the fourth-order generalized OCl.scheme) will appear

in a forthcoming report. 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.

36
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. We need to consider boat 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) memory becomes a

severe limitation even without attempting to optimize the operation count.

The storage requirements for the solution of 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 CAC 7600. One must

then resort to either mass storage devices such as disks or slow access memory

large core).

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

When small amounts of data are being transferred frequently (what we are considering)

then access time becomes a significant factor. Therefore, a combination of strategies

must be employed in order to optimize both access time and transfer rate. As a result

of these necessary manipulations, the resulting code was far more complex than the

one that would have been developed specifically for a machine with "unlimited" storage

resources.

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

the QR Operator technique leads to the conclusion that the most significant i fraction a ction 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 on 2 and 3 as they were computed in the first sweep on the

right-hand side of the differential equation. On the second and third sweeps I

and 3 

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

n, 

that the formulation of the LBI scheme referred to linearized operators ,^ i s

instead of I n , s on the right-hand side of the equation. 

The general structure of the computer code will not 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

37

Z. 4 ,. ,a^^; ;^ _^,.....^!`.^ ...• ^ - - -.mss
flow field and stored for ready access when needed for the comput;)Vion of the appropriate terms in the governizig equations. Thereafter the terms that are to be lagged

(treated explicitly) are evaluated and absorbed into the function Sn

The operators I1,J02 and .13 are then computed. Those 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.11 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 1 1

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. The

third sweep equations consist of the two momentum equations and the continuity equation, with the -03 operator being accessed from memory. The resulting 3 x 3 system of

equations are solved for the three velocity components.

After tho ADI procedure is completed, the thermodynamic quantities, density,

temperature and viscosity are computed. The procedure is then repeated at the following

time step.

It is noteworthy that the scheme just described operates on vectors, i.e. lines

of data. Therefore, it could show promise for vectorized machines.

38

*4.
NUMERICAL RESULTS

In this section we describe the numerical results obtained by using the computer

code described in the previous sections on several Cartesian test problems. The goal

of these calculations is to validate the computer code. No specific timing experiments have yet been conducted to evaluate the code's efficiency. As mentioned in the

previous section, the code was structured for Oe calculation of three-dimensional

problems so that the disk writing and large core data allocation are necessary to run

the code efficiently in that mode. For two-dimensional problems these artifices,

although not necessary except for very large mesh problems, are still employed. Hence

as an observation, it .3 noted that the code ie not as efficient in the two-dimensional

mode as it could be.

Since the primary goal in this portion of the effort was to obtain a working

computer code, the following test cases are considered:

1. One-dimensional unsteady oscillating flat plate;

2. Three-dimensional boundary layer without pressure gradient; and

3. Two-dimensional steady boundary layer with and without pressure gradient.

These cases check out three basic features of the code, viz., time dependent behavior,

pressure gradient effects and three-dimensional effects.

All the calculations employ second order finite difference techniques. Although

a fourth order generalized OCI option is incorporated into the code (and the code is

structured particularly to encompass that scheme), to date none of the cases considered

were run employing this scheme. In subsequent work fourth order calculations will be

compared to second order finite difference as they relate to efficiency gains and

convergence rates.

Oscillating Flat Plate

The first case considered is that of an oscillating flat plate in a unifo,14 stream.

It is a one-dimensional unsteady problem and is thus only a function of the normal

distance to the wall and time.

Schlich r '::g (Ref. 38) gives the exact solution to this problem, which he terms

'Stokes second problem', as

u(77,t) U CO e-'1coS(wt -r7)

39
where n . y 3 w and w . the frequency of oscillation,

2v

The boundary conditions are

at the wall

oq = 0 U(0 1 t) = Uapcos wt

v(O,t) = 0

at the outer edge

7-- a0 U(O,t) = 0

Since we must specify the velocity r; a finite value of n we set u at n = ne to the

•' exact value of u(ne ,t). As initial Londitions, the velocity field is set to the exact

value at t = 0, Since the solution does not take into account any initial conditions

a

of the fluid, the relationship for u(n,t) represents a "steady state solution".

The computer is run in the two-dimensional mode, with 5 grid points in the

streamwise direction, so that additional boundary conditions are required at the inflow

and exit planes. This is accomplished by setting the first derivative of the velocity

components to zero there, which eliminates any variation in the streamwise direction.

For the problem considered, the following parameters are employed.

w = 100 rad/sec

Re = 100

ne = 5.0

A uniformly distributed grid of 11 points in the normal direction is used as well

as a fixed time step of wet = 5°. Second order central differences are employed for

the calculation which includes all the viscous terms appearing in the governing equation,

Since there is no variation in the streamwise direction and the normal velocity is zero

throughout, the governing equations reduce computationally to their one-dimensional

counterpart. In Fig. 3 the computed skin friction coefficient is compared to the

analytical result over one and one-half time cycles. The analytical value is

40
Cf - Re -. V a 1 (.14- - wfl

Even with a relatively coarse mesh there is good agreement between theory and the

numerical results. Although the Crank-Nicolson results (s - 1/2) are more accurate

than the fully implicit calculation (S = 1) as can be seen in Fig. 4 where the L2

M 

error of each calculation is compared, the skin friction results appear much closer.

We can account for the discrepancy by looking at the evaluation of the skin friction

coefficient.

The error in the skin friction coefficient is made up in part from the error ;gin

the solution of the differential equation and in part from the approximation of the

first derivative in the computation of the shear. The error in the evaluation of the

first derivative is given as

3

error (, 7,1,, = 

0 

) - at l u 

'1 ` a2 7 2v )

where al and a2 are bounded constants

For our case v2 = 50 so that the total truncation error is dominated by this term.

The results of this calculation indicate that accurate time dependent calculations

can be obtained.

Three-Dimensional Flat Plate Boundary Layer

The second problem considered is the three-dimensional viscous flow over a flat

plate skewed at an angle of 45° to the flow direction. In Fig. 5 the computational

domain is shown in which the solution is obtained. An equally spaced mesh of 11

points in each of the three directions is employed.

Two cases are considered. In the first case boundary layer assumptions are

made, and the only viscous terms retained in the momentum equation are those that

appear in the standard three-dimensional boundary layer equations. Diffusion terms in

the streamwise and spanwise directions are ommitted. In addition, first order approximations to the streamwise and spanwise first derivatives (in the marchingdirections)

41
are made. The second case considered retains all the viscous terms. The derivative

approximation to the spatial operators employed entered differences with artificial

viscos i ty (cf Table V).

The boundary conditions applied are no slip at the surface and velocity specified

at the outer edge of the viscous layer. There the velocity components are set to lhr-2—.

Note that the vector sum of the two components is 1. At the inflow boundaries 1 and 2

(cf Fig. 5) the velocity profiles are set while at the outflow boundaries 3 and 4

second derivative extrapolation conditions are employed.

Both cases compared well with one another and to the Blasius solution. In Fig. 6

the displacement thickness b* and the momentum thickness t, along the diagonal of the

computational domain are compared to the theoretical result. Although symmetry was

retained across the diagonal, at a constant ^ value, there was some variation of the

integrated properties as a function of g from their values on the diagonal.

Two-Dimensional Howarth Flow

The third case considered is a two-dimensional steady flow problem; the laminar

boundary layer on flat plate in the presence of an adverse pressure gradient, the

Howarth flow. Here again we consider two cases. The first employs the boundary layer

approximations while the second retains all the terms in the streamwise momentum

equations.

A 21 x 21 mesh is used in both cases. The boundary conditions employed are

no slip at the wall, streamwise velocity specified at the outer edge, velocity profile

specified at the inflow plane and second derivative extrapolation at the downstream

plane.

The external velocity field is linearly retarded, i.e. u = 1 - Kx, so that

there is an adverse pressure gradient. Due to the adverse pressure gradient separation will occur at some point downstream. Howarth (Ref. 39) computed the separation

point at a value of x* = Kx = .1200. In the calculations considered here we choose

our domain to span .05 x* .12 which terminates at a lorar;nn --cnn M no f-n

Howarth's predicted separation point.

42

.—:..i ^_
The initial conditions for the problem were constructed from the local Falkner

Skan similarity solution with the appropriate pressure gradient parameter ,imposed

at the corresponding streamwise location.

In Fig. 7 the computed skin friction coefficients for both cases are compared to

the results of Howarth. The case considered corresponds to Briley's cases (Ref. 40),

i.e. Re - 62,500 and k = 3. There is good agreement over most of the domain, with an

overprediction of the shear as the separation point is approached. This diserepency

is not surprising due to the coarse mesh employed. It is interesting that the

"boundary layer" case gives somewhat better predictions near the separation point

than the "full equation" case. This is probably due to the lack of upstream influence

in the boundary layer solution which did not allow the flow to adjust upstream of

separation. Both calculations converged within 30 - 40 iterations, with a maximum

change between two (,Large) time steps being less than 5*10-4.

43
44

.. ^ .rw ,+^'ry.-rya _^.^_ "' ••.r' ice'' •^'+^''r"' s, . . ^• ^^

CONCLUSIONS

In this report a computer code is described that can be applied to the solution

of an approximate form of the time-dependent Navier-Stokes equations over airfoil

sections. The governing equations are more general 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, as in conventional

boundary layer theory. The computer code which treats the governing equations incorporates the split LBI scheme in conjunction with QR operator scheme that permits a

variety of spatial difference schemes, including standard second order finite differences,

exponential type methods and fourth order OCI techniques. In the split LBI scheme, an

implicit sweep is 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 be eliminated. Results of several Cartesian problems employing second order

finite differences indicated that the proposed method is viable and has application to

more complex flows. Future efforts will aim at exercising the fourth order OCI schemes

option, treating actual airfoil sections and considering turbulent flow. Comparisons

will also be made among the various types of spatial schemes to assess the overall

efficiency gains made by employing higher order methods. The treatment of inflow and

outflow boundary conditions will also be of primary interest.
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

0

^(x,t):

d¢/dt =F(0) dG(OWX (Al)

Special cases of Eq. (Al) include the conservation form if F(^) - 1, and

quasi-linear flow if G(0) = ^. Previous implicit methods for Eq. (Al)

which employ nonlinear difference equations and also methods based on twostep predictor-corrector schemes are discussed by Ames (Ref. 41, p. 82) and

von Rosenburg (Ref. 42), p. 56). O.ne 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 procedure 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. 43) which is applicable to the quasilinear case (G = 0) of Eq. (Al).

The first step of their procedure is to linearize the equation by evaluating

the nonlinear coefficient as F(On) and to predict values of ¢n+112 using either

the backward difference or the Crank-Nicolson scheme. Values for ^n+1 

are

then computed in a similar manner using F(¢ n+1/2) 

and the Crank-Nicolson scheme.

Gourlay and Morris (Ref. 44) 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 6(0) by the relation G(¢) = rG^ (¢) when such a

definition exists, and to evaluate G(¢n+l) using values for 0n+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 0n+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

45

,_ w. ^ 21^.e wat 

.h* _ I :—++^t _^VK, _ id. ^^..^^^y►

f' F' .t!„t4s.., ,!^. „,,
on an expansion of nonlinear implicit terms about the solution at the known

time level, tn , 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 Morton (Ref. 29, p. 203) to a scalar

nonlinear diffusion equation. Here, the technique is developed for problems

governed by R nonlinear equations in t 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, Ayl, Ay2 and Ay in the y 1 , y` and y 

3

directions, respectively, and an arbitrary time step, At. The subscripts i, J,

k and superscript n are grid point indices associated with y1 

y2, 

y3 and t,

respectively, and thus¢i ,J,k denotes O(yi, yj, y 3 , tn). It is assumed that

the solution is known at the n level, t n , and is desired at the (n+l) level,

t o+l . At the risk of an occasional ambiguity, one or more of the subscripts

is frequently omitted, so that ¢n is equivalent to 0ni,j,k.

The numerical method employed is quite general and is formally derived for

systems of governing equations which have the following form:

JH (0)/dt =-2)(0) +S($) (A2)

where ¢ 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/a^ 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 one 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.

46
The linearized difference approximation is derived from the following

implicit time-difference replacement of Eq. (A2):

(Hn+l-Hn)/.&f=p'i-2)( ,n+I)+Sn+Ij+(I(¢,n)+Sn] (A3)

where, for example, H

n+1 

r 

H(On+l). The form of 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 (0)

produces a backward difference formulation for S - 1 and a Crank-Ni.colson

formulation for 3 = 1/2.

The linearization is performed by a two-step process of expansion about

the known time level to and subsequent approximation of the quantity

(a¢/at) net, which arises from chain rule differentiation, by (0n+1 - 0n). The

result is

Hn+I = H +(aH/d,,)n ( ,bn +I _,bn) + 0 (Lt)2 (A4a)

S n+l = S n +(dS/d¢ )n (on+I _.p n) +0 ( 00 2

(A4b)

b(O n +I) _

n) +(a.^/a^)(^n+I_4,n)+0(Qt)2 (A4c)

The matrices aH/a^ and aS 180 are standard Jacobians whose elements are defined,

for example, by (aH/3¢) qr = aHg /a¢r . The operator elements of the matrix

D2 /ao are similarly ordered, i.e., (a.2^ /a¢) qr 5 3.2 q /a¢ r ; however, the

intended meaning of the operator elements requires some clarification. For

the q th row, the operation (a.`^ q/9W(On+1 _ 

¢n) is understood to mean that

fa/at2q (O(x,y,,Z,t))) net is computed and that all occurrences of (a¢r/at)n

w arising from chain rule differentiation are replaced by (0n+1 - 0r)/At.

After linearization as in Eqs. (A4), Eq. (A3) becomes the following linear

implicit time-differenced scheme:

47
(dHn/a¢)(On +t_on)/At --.2)(On) + Sn + /3 ( a `h/ d¢+aSn /ao)(On+I—on ) (A5)

Although H

n+l is linearized to second order in Eq. (M), 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 8 1/2) by

a change in dependent variable to ¢ _- H(¢), provided this is convenient, since

the nodlinear 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

(

0n+l _ 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 - Vin) rather than ¢n+l. This both simplifies Eq. (A5) and reduces

roundoff errors, since it is presumably better to compute a small 0(At) change

in an 0(l) quantity than the quantity itself. To simplify the notation, a

new dependent variable ^ defined by

, = ¢_0n (A6)

is introduced, and thus ^n+l = ^n+l _ ^n, and ^n = 0. It is also convenient

to rewrite Eq. (A5) in the following simplified form:

(A+ At .1') 4 n+I = At [-',6 (00 )+S"] (A7a)

where the following symbols have been introduced to simplify the notation:

A—= aH n /a(k —Rot(as"iao) (A7b)

.1 = - 13 (a .rb /do) (A7c)

It is noted that .1(*) is a linear transformation and thus _P(0) = 0. Furthermore if 2(^) is linear, then 2(^ ) =

48

s

.e t wt ♦ +b:ik " it
Spatial. differencing of £q. (A7a) is accomplished simply by replacing

derivative operators such as D /Dy i, 82 /8yiay i by corresponding finite

difference operators, Di , Di. Henceforth, it is assumed that 2Z and I 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. 45).

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 coupled 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. 30).

49
APPENDIX h

Tensor Notation and General Coordinates

In tho space R with coordinates (x l, x2 , x3 , ,xn) let the position

vector be denoted by r. A basis can be formed by constructing the covariant

vectors e 

ar

e, _ ax

i

The inner (dot) product of these vectors forms the function g ij , which we will

refer to as the metric tensor

91j ° e, ' ej

The reason for this terminology will be described subsequently. In threedimension g ij has the form

911 9 12 913

g l j g21 9 22 923

931 9 32 933

The metric tensor g ij is symmetric in the indices i and ,j so that there are 6

independent components. For an orthogonal coordinate system the off diagonal terms

g ip i#j are zero.

The differential arc length ds is defined as

(ds) 2 = dr • dr

a i a d,,_axer dxj

a X7

= dx I dxj e e j

ds 2 = 9 ij dxi dxj

(al)

(B2)

(B3)

50

,.R,.._, ^^.^.^---•-_ . _ ^. ^ . _ mow „

;^ -^'ri,,,^ ^, t ;,.^:;,r.-_ -^-^_^-----_____ 

... a .t y^it row .^k" . Ys ..All ./^.
OF POOR QUALITY

Expanding we obtain

ds 2 = q l l dx l dx l + g22dx 2 dx 2 + 933 dx 3 dx 3 04)

+ 2 gl2dxl dx ?- + 2g13 dx i dx 3 + 2g23dx2dx3

For orthogonal curvilinear coordinates the arc length reduces to

A ds2 gll dx l dx l + g22 dx 2 dx2 +g33dx3 dx 3 ($5)

so that gll , 9221 933 can be identified with the scale factors as hl, h2 and h3

respectively.

The Kronecker delta function is defined as

if1 I = j

e f ej ^ sj -- to 10 j (a6)

Notes on raring and lowering indices

(a) multiplying by 6 1J changee the index but keeps it in same position

j

u j SI

= ui

U I s1 = uj

(b) multiplying by g i3 and gij

ui g lj = U lowers index

U i glj ° Uj raises index

(r_) multiplying 6 1J by the metric tensor gjk

SI gjk s ki raises index

1

s1 

9 j = 

g1k substitutes index

51
OF POOR QUALITY

Hence we obtain the following relationship

slk g1k

Note that 6 i is not the Kronecker delta function in general. Only in the case of

Cartesian coordinates where gik R 1 and i - k is 6 i the Kronecker delta function.

In order to obtain the contravariant basis ei' we write

-- ^- k

e l ' e j g lj = glkSI

Employing equation (B-6)

r^

r{

el, ej = 9Ik Ce k'ej J

we obtain

e)91kek

f

(B7)

r

Hence the metric tensor alk relates the covariant to the contravariant

base vectors.

We may now define g

jk 

as

9jk = e j e k (B8)

We wish to relate g jk to gjk . Dotting equation (B-7) by ej we obtain

j k j

e = 91k -e

gi = glkg

j kj (B9)

Now from the definition of g ij (Equation B-8) we obtain the desired result

sj el =gij= 91k8 j=91keje k k

52
(Blo)

ORIGINAL PAGE IS

OF POOR QUALITY

Therefore, g also relates the covariant and contravariant basis vectors, i.e.

e I = 91k' k

The metric tensors gij and gi3 can also be W ised to raise and lower indices.

For example, consider the product of the tensors Cjk, C3k C3 k , and Ci k with gij

or gij

Cjk 

gij

CI

Cjk gi j = Clk

C.j k g lj = Cki

C kj gij = elk

Derivatives

The derivative of a basis vector ei is a vector in Rn so that it can be

represented ar a linear combination of other basis vectors, e s. Let e i e b a

vector, then its derivative with res p ect to xi is

aX, = D, e i = 

rjk ek

Dotting this expression with e^ we obtain

Dj e^ • et = r^ ek .eA = r ^k = rj'

R,

where rji is the Christoffel symbol.

r^= 

Djei 

.e

53

MW
ORIGINAL PAGE 19

OF POOR QUALITY

R

Note that the Christoffel symbol Pii is not a third order tensor.

Representing rii in terms of the g i3 's we obtain

a

9r„j = emIej

i

DI 9mj= Dl em• ej = D i em -e j + em , D I ej

r 

rlmek 

ej 

+em

1 j ek

i 

i

k k

t p 9mj

rm9kj 

+ 

9mkrj

Similarly we can show that

k k

Dj 9mi

= 9

k rjm + 9mkrjl

and

kk

DM91J = - 9kjrmi - 9I^., ! mj

Employing the symmetry property of the P's in their lower index, i.e.

k k

r

j 

= rj I

and adding the three derivatives we obtain

k 9 k

Iij = 2 DI 9 mj+ Dj 9ml - Dm9lj

which relates the Christoffel symbol to the derivatives of the components of the

metric tensor.

54
(b)

ORIGINAL, PAGE 04

The Kronecker delta function is defined as

OF POOR QUALITY

s k - e ...k

j `ej

Taking its derivative, we obtain

DI 8jk=

0= Dlek-ej

. 0 = D, ek ej + ek , DI e j

Substituting in the expression for the derivative of a basks vector

M

D i ej = rij em

and rearranging,

D^ek . e j = - ek - r. gm = — rj 8 = —rij

we obtain

rjk - —Diek•ej

which relates the Christoffel symbol to the derivative of a contravariant basis

vector.

Other properties of Christoffel symbols are

(a) r gk.C= r 

z = 

I, k 

se

t}k i} Ij k

k

rij 9kt - r je

There is also symmetry of lower indices, i.e.

rljk = rjik

r

f r "
ORIGINAL PACE IS

OF POOR QUALITY

Covariant Derivative

Consider the derivative of a vector u, expressed in a covariant basis

r

r, U =Uiei

_ Dj 

U i ei 

T Dj 

Ui ei + Ui 

Djei

Now the derivative of a basis vector is

Dj ei _ rj ek

so that the last term in the expression becomes

r

Ui Dj ei U i I j k ek

w

Since k and i are dummy indices, we may interchange indices, i.e.

U i Dj ei U i rj -6k = U k r' j ei

Combining, we obtain the desired expression

d

( D jU i + Ukrkj)ei

The term in the parentheses is called the covariant derivative and is denoted

as follows

du = Ui l ei

dx1 j

F

In a Cartesian coordinate system the covariant derivative reduces to the partial

i

derivative, i.e.

Ui 1 j = Ui,j

i

' Consider the derivative of a vector in a contravariant basis which is given by

D i Uj e = Di Uj -d

i

+ UjDie

i

56

i
OF POOR QUALITY

The last term can be expressed in terms of Christoffel symbols, i.e.

uj Di6 u jriikek=ukrjel

Thus we obtain the desired result

Di u jej = (D i uj - uk ij e j

= (ujri - u k ri j) ej

Di u jej = ujliej

Now in Cartesian coordinates the covariant derivative of the metric tensor is zero,

9ij,k = Sij'k 0

Since this is a tensor equation it is valid not only in a Cartesian system but in

every coordinate system, i.e.

9ijl k = 0

Consider now the vector uR.

_

U _ u k 91k

Taking the covariant derivative of u R and using the previous relationship, we obtain

the desired result

us11 = (uk90) ( i = u k l i 9Ak+ uk9`^i

uL i = u

k

Ii 9Qk

57
An alternate expression for the derivative of a covariant vector is given by

D i u U111 e i = uj 11 gzje4

U 11 ej

Note that the dummy indices are interchanged.

Vector Operators

In this section we consider the vector operators, gradient and divergence

operating on scalars, vectors and tensors.

Gradient

Scalar (^)

The gradient of a scalar is given as

Vu = ek®Dk^ = gkmem@Dk^

where @ denotes a tensor product.

Vector (-T = u -4,i )

The gradient of a vector becomes

0u = e k ©DkuI e^ eke l Dk ufef + ui Dkei^

• uIIkekC^ei

• gmk u i Ik em Or,i

58
ORIGINAL. PACE 19

OF POOR QUALITY

which is a second order tensor.

Tensor

Consider a tensor of the form v = aj el ©ej . The gradient of this tensor

can be obtained as follows

M

0u = ek ®Dk olj el ®ej

=ek 8{Dkaj (TI ©ej ) +alj ^(Dk el ) ®ej +oljel 0Dkej^

=ek © {o alj(e 8e )+aij rrye 8e + alje ®I''^e} 

k I j

'kl L j i kj

which yields a third order tensor

_ 

0u = { D

k° Ij + ° jrkc + o'` rk,j e k ©ei

 ®ej

ll

f Dka ij + a j i , + 

aitrtl 

11 

gmem O ei De

^ 

= (aij Ik) 9'n M ©el Oej

59
ORIGINAL PAGE 19

OF POOR QUALITY

Divergence

The divergence of a vector of the form u^e^) can be obtained as follc •s

0.7 e k • Dku^ei

=e k •uilk e^

 - ul Ike' •el

ui lk slk

D • u =uklk

From the definition of the covariant derivative we obtain

uk l k Dkuk + u1 J k

It can be shown (Ref. 13) that

k 

* 1k J Di d

Substituting this expression into the definition of the covariant derivative we

obtain an alternate form of the divergence.involving the Jacobian

k) k 1 i U k = Dk u + u D i d = D k uk + uk i Dkd

U 

Ik - Dk ( d u k )

Note that this relationship is in conservation form.

60

,...

,r

k ..coi - ' i

,I^'..+^'n. F *'^!^°'. -"•s' '* __;L4!^M-
0 • u = [ Di a il + ai1 ^L + aitr^j ej

which reduces to

ORIGINAL, PAGE 19

OF POOR QUALITY

Consider now the tensor of the Form (u = a tj et © ej

Again we obtain its divergence as follows

V ' U = V'a ij e ®e t j

e k• Dk ( a ijei ©ej )

=:e k ID k a^ ij e Oej + a ij De ki ©e + a @D el tje kejl

= DO e k • et© ej + aiJ 

r. 

ek ' e^ ©ej + a'^ek • ei © rkj ec

= D k

al j Sk i}e, + a ij rt sk e + ai1 rt sk ki Q j kj i e

by i

nterchanging indices j and k we obtain

V- =u = aijli ej

61
OF POOR QUALITY

or

j 0 • u = l J Di(ai J )+allrii^aj

We will consider nowthe tensor of the following form (u = aj ei ©ej) ,^ j

The divergence of u becomes

D . u = 6 k

1 DO ei ®ej + aj Dkei 0-0 + aj e i ®Dk ej^

= ek-[D a e Oe j + a i rkN Oe - aj a' F, @I kj e^] k j i j i G j i kb

Rearranging and interchanging dummy indices we obtain

V • u= IDk aj +ajrkl — at

rkj

^ rk ^ei)©ej

= [ D I aj + al rE - °yL rj ^ ej

Now the contravariant basis vector is defined as

e j = ek gjk

Substituting into the expression for V-u and with some manipulation we obtain

the divergence of the second order tensor as

0 • u Di aj + ° rid aE rj J gjk ek

0 • u - 

L J Di (a 1 J) - l 

Ij J g,kek

62
63

._

In order to obtain D •u in the e coordinate system we employ the property of the

metric tensor to raise and lower indices

0'u ^Dk aj + a^ rk1^_ aI Ff' ek e} ©ef

ek. e = gk^

and obtain the desired result,

0 o u l L Dk a^ + a^ 114 - a^ 

rk, gkj} el
REFERENCES

1. Kendall, R.M., Bonnett, W.S., Nardo, CT,, 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. Braley, W.R., and McDonald, Ii.; On the Structure and Use of Linearized Block

ADI and Related Schemes. J. Comp. Physics, Vol. 34, No. 1, 3.980.

3. Weinberg, B.C., and McDonald H.: Solution of Three-Dimensional Time Dependent

Viscous Flows, Part 1: Investigation of Candidate A13orithms. SRA Report

R79-90004, Final Contrnctor's Report (NAS2-10016), 1979.

4. McCroskey, W.J.: Some Current Research in Unsteady Fluid Dynamics. J. Fluids 

Engineering, Vol. 99, March 1977, pp. 8-38.

5. Nash, J.F., and Patel, V.C.: Three-Dimensional T=;rbulent Boundary Layers. 

r^ 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,

P p • 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 I. 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,

1973. 

14. Gibeling, H.J., Shamroth, S.J., and Eiseman, P.R.: Analysis of StrongInteraction Dynamic Stall for Laminar Flow on Airfoils. NASA CR-2969,

April 1978.

64

4 ^, c

_. ^.A/.t^k'r •^ ,w,::^ it ,!aR ^,^,.: ^ _ _ ..^ .r...r=

'^r
15. Favre, A.: Equations des Gaz Turbulents Compressibles. J. de Mecanique,

Vol. 4, pp. 361-392, 1965,

16. Berger, A.E., Solomon, J.M., Ciment, M., Leventhal, 8.H., and Weinberg, B.C.:

Generalized OCI Schemes for Boundary Layer Problems, Math Camp., Vol. 35,

No. 151, 1980.

17. Allen, D.N., and Soutiwell, R.V.: Relaxation Methods Applied to Determine the

Motion, in Two Dimensions, of a Viscous Fluid Past a Fixed Cylinder, Quart.

J. Mech. Appl. Math., A, 1955 0 pp. 129-145.

18. I1"in, 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 0 1969, pp. 596-602.

19. McDonald, H.: A Novel Finite Difference Scheme Particularly Suited for Fluid

Flow Problems Containing Discontinuities. Scientific Research Associates

Report P76-1, August 1976.

20. Dennis, S.C.R.: Finite Differences Associated with Second-Order Differential

Equations, Quart. J. Mech. Appl. Math., Vol. 23, 1960, pp. 457-507.

21. E1-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.

22. El-Mistikawy, T,:.. and Werle, M.J.: Numerical Method for Boundary Layers with

Blowing - the Exponential Box Scheme. AIAA J., 16, 1978, pp. 749-751.

23. 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.

24. Peaceman, D.W., and Rachford, H.H.: The Numerical Solution of Parabolic and

Elliptic Differential Equations. Soc. for Indust. Appl. Math., Vol. 3, 1955,

pp. 28-41.

25. Douglas, J.: On the Numerical Integration ofu = u by Implicit Methods.

Soc. for Indust. Appl. Math., Vol. 3, 1955, pp.^ i2-65. t

26. Mitchell, A.R.: Computational Methods in Partial Differential Equations.

Wiley, New York, 1969.

27. Yanenko, N.N.: The Method of Fractional Steps. Springer-Verlag, 1971.

28. Douglas, J., and Gunn, J.E.: A General Formulation of Alternating Direction

Methods. Numerische Math., Vol. 6, 1964, p. 428.

29. Richtmyer, R.D., and Morton, K.W.: Difference Methods for Initial Value

Problems. Second Edition. In'tarscience Publishers, New York, New York, 1967.

30. Briley, W.?..,, and McDonald, H.: Solution of the Multidimensional Compressible

Navier-Stokes Equations by a Generalized Implicit Method. J. Comp. Physics,

Vol. 24, No. 4, August, 1977, p. 372.

65
31. Gibeling, H.J., McDonald, it., and Briley, W.R.: Development of a ThreeDimensional Combustor Flow Analysis, Volume I: Theoretical Studies. Technical

Report AFAPL-TR-75-59, Vol. 1, July 1975.

32. Gibeling, H.J., McDonald,^ i,to, and Briley, W.R.: Development of a ThreeDimensional Combustor Flo, Analysis, Voluma 11: Theoretical Studies. Technical

Report AFAPL-TR-75- 59, Vol. 1I, Octob er, 1976.

33, Leventhal, S.H.: The Operator Compact Implicit Method for Reservoir Simulation,

Proceedings of Fif:', SPE Symposium on Numerical. Resevoir Simulation, 1979.

34. 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 Hampshire, 1975.

35. 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

3,:d AIAA Comp. Fluid Dynamics Conf., Albuquerque, N.M., 1977.

36. 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.

37, Fairweather, G., and Mitchell, A.K.: A New Computational Procedure for ADI

Methods, SIAM J, Numer. Anal. Vol. 4, No. 2, 1967.

38. Schlichting, H,: Boundary Layer Theory, Fourth Edition, McGraw-Hill Book Co., Inc.

New York, 1960.

39. Howarth, L.: On the Solution of the Laminar Boundary Layer Equations, Proc. Rey.

Soc. A. 164, 1938.

40. Briley, W.R.: A Numerical Study of Laminar Separation Bubbles Using the NavierStokes Equations, J. Fluid Mech., Vol. 47, Part 4, pp. 713-736, 1971•

41. Ames, W.F.: Numerical Methods for Partial Differential Equations. Barnes &

Noble, Inc., New York, New York, 1969.

42. von Rosenberg, D.A.: Methods for the Numerical Solution of Partial Differential

Equations. American Elsevier Publishing Co., Inc., New York, New York, 1969.

43. Douglas, J., and Jones, B.F.: On Predictor-Corrector Methods for Nonlinear

Parabolic Differential Equtions. Soc. for Indust. Appl. Math., Vol. 11, 1963,

pp. 195-204.

44. Gourlay, A.R., and Morris, J.L.: Finite-Difference Methods for Nonlinear Hy'eerbolic Systems. Math. Comp., Vol. 22, 1968, pp. 25-39.

45, Bellman, R.E., and Kalaba, R.E.: Quasilinearization and Nonlinear BoundaryValue Problems. American Elsevier Publ. Co., Inc., New 'York, 1965.

66
ORIGINAL PAGE 13

OF POOR QUALITY

F

f

^ F

A

,w

t

f ^^

t

$`t

TABLE I. - OPERATOR COEFFICIENTS FOR STANDARD

OPERATOR COMPACT IMPLICIT SCHEME

q j ` 6 - 5'pj + 2pj+i pjpj+i

q jc 60 + 16pj +i - 16pj- i - 4pjl pj +i

q j 6 + 5p - 2pj-i - pjpj-i

rj- q ^ \ I - 2 `pj-1) + q l ( - I I pi) + qJ ^I + 2 pj+i) + h2q 1 cj-i

rj _ ( rj+ + rj ) + h2(gj cj_i + qcj + qj+cj+')

rj qj pj-s +

^I — z ) + q j ( I + 2 p' ) + q j

+(I + 2 pj+i) + h2g j c j +I

where

pj•hbj

67

n

AMP
ORIGINAL PAGE 19

OF POOR QUALITY

TABLE II. - OPERATOR COEFFICIENTS FOR GENERALIZED

OPERATOR COMPACT IMPLICIT SCHEME

q 6+1pt-3^Rcj +rpZ]Rc1t l

q^ n 60 +^IOpI^Rc1 

+ L P 3) Rcj + LT1*1P4)Rc1

q1+ . 6 + I pt + 3] Rcj + [p. +p.]Rc

12 + [p4]Rc 1 3

where

Pt ' 3 , PZ ' O 9 P3 ' max I VI 

, 7r2 1

7rt - (T) +I+TJ-t)p2

 + Tj+t pi + W, irZ . 15 2pZ + ( o-z - I)pt - 3 ( Tj+t+a-Z) + Wr

0 0-1 2:0_ J 2pt- v'2 >O

rt a a'(10 - Tj+t - Tj-t) 0-1:50

^Z (2pt- o.2)2 /8 2p,-c,<O

vt ` Pt 3 + ^G Tj T+t) 2a-? - 3Tj+t - Tj-t + 10 + 2hTj-t \ b _t

1 +t jt 1

t cJ-t )

P4

2 (+ Tj+,1 ^3 n3 ' P3 7ri + t + 2^rj-t 2 + h bJ

 -1 p 

with h sufficiently small so that

10bj- bj- 1

- bj+1 > 0 and 2 + hc j+,,IbJ41 > 0 for j-2, • ••,J and cjs0

where

F ; Tj_t " bl_ t /bJ s T)+t " bj+t /bj and Rcj - hbj

ri , r j, r j given in TABLE I

68
TABLE III.- OPERATOR COEFFICIENTS FOR

ALLEN SOUTHWELL EXPONENTIAL SCHEME

r — = Rcj a 

-Re 

J / 0 - e-Rcj)

rj+

Rej /( I-a -Rol)

r jc = — Rcj + cj

qj-= 0

qj°=i

qj+= 0

where Rc j = h b j

1

69

I
TABLE IV. - OPERATOR COEFFICIENTS FOR EL-MISTIKAWY WERLE

EXPONENTIAL BOX SCHEME

where

r J- - p-exp (-P-) / 11 - exp(-p-)]

r J+ - p+/ ^I -exp(-P+)] 

r^ _ -(r i + + rj-)

q j- _ (1-rJ

-

) /(2p-)

qj - (rJ+ I) / ( 2p+)

q j ' Qj + qj

p- _ ( PJ-i + PJ ) P+ ` 2 ( PJ + PJ+i)

P) = hbJ

and

70
TABLE V,- OPERATOR COEFFICIENTS FOR SECOND ORDER

FINITE DIFFERENCES WITH ARTIFICIAL VISCOSITY

N

^a

r

r ) - = t S

Jc = -2 + ct

rj+ = I + S

Cl - 0

c

=t

qj

q^ = 0

where S = I for Rc j > 2

S = - I for Rcj < -2

t = -21Rcj^

and Rci = hbi

for ,Rcj ) <2, r3 , r^, r^, q^, q^, qj

reduce to standard finite differences.

71
G

O

.rj

u

a

N

O

U

n

p

O

.r-j

41

U

u

cn

rl

.H

O

W

r-I

cC

u

.H

a

H^

H

GO

W

W

G

cz

r-4

PW

r1

\J

N

M

1

ORIGINAL PAGE 19

OF POOR QUALITY

72
73

ORIGINAL PAGE 19

OF poop QUALITY T

CN

PIN

2

4J

ro

0

8

0

0 'rj

b6D

0

H

O 1J

0

z

-H

0

9

U)

'^

al

td

0

4-4

0

U

a)

4.

4-4

C14

F34

.H

M

P

0

0

O

Cl)

0)

4j

cU

C14 

"r,

-1
ORIGINAL PAGE 18

OF POOR QUALITY

M

N

w1

G1 u

a

N k'

bq

t0

rJ

r-1

.H

u

N

O

N

N w

u

.H

w

u.,

v

OU

O.rj

LJ

u

.H

$4

W

G

x

I

ri

t0 -H

W

J

NO O

r—I r I O O O ri r{

Cl!

W

u

74
CJ

N tdH

N

-W

4

0

H

0

44

FX4

H0

CY)

co

O

Is

0

C-1

49

H

0

N

-H

P4

f

ORIGINAL PAGE 19

OF POOR QUALITY

75

0
ORIGINAL PAGE It

OF POOR QUALITY

W

r•7

M

76

6oNW3OWNN

tO

r-1

Iq

r-q

cdCO

.HN

C3 Ad

a^H4HNOW

.HoOA

r-1

cd

C+OL

fd

MW

rzOUI

Ln

00

.r{ W
j0

as

ro A

OqRO

.0N

1JNANu

cs

rz 0zbNyq

Uq

r^U

.r{

.0H

tJG

U1EaU

r-ia

r

ORIGINAL PAGE 19

OF POOR QUALITY

pp

+D

0o

RO

uj

00

p

00

p

0 ^^ 00 N

bOJ

'^ L

N ^

o a w

x O

H u

0

0

110 W
ORIGINAL PAGE 19

OF POOR QUALITY

0 Full Equation

1.2

<> Boundary Layer

1.1

Howarth Solution

1.0

.9

c .8

.7

.6

.5

.4

.3

.2 0 0 0 oIX

0

0 .01 02 .03 .04 .05 .06 .07 .08 1 09 .10 11 .12

X*

Fig. 7 - Skin Friction Disti-11,uHon for 'r-,,o-D! ronsional Hovarth Flow.

78

Ik