HTL-51, CFD-21

Project 1928

ISU-ERI-Ames-90030

S. G. Gegg R.H. Pletcher J.L. Steger

A Dual Potential Formulation

of the Navier-Stokes Equations

Final Report

Funds for the support of this study have been allocated

by the NASA-Ames Research Center, Moffett Field, California

under Interchange No. NCA2-174

Heat Transfer Laboratory

Department of Mechanical Engineering

Computational Fluid Dynamics Center

enginee, r.in.g

researcn lnslltule

iowa state university
,4,

In

TABLE OF CONTENTS

LIST OF TABLES ................................ vii

LIST OF FIGURES ............................... ix

ABSTRACT ................................... xvii

ACKNOWLEDGMENTS ........................... xix

NOMENCLATURE ............................... xxi

lo INTRODUCTION ............................. 1

1.1 Background ................................ 1

1.2 Primitive Variable Methods ....................... 3

1.3 Non-primitive Variable Methods ..................... 4

1.3.1 Vorticity/velocity approach ................... 6

1.3.2 Dual potential approach ..................... 7

1.4 Scope of the Present Study ........................ 13

1.5 Organization ............................... 15

2. DUAL POTENTIAL FORMULATION ............... 17

2.1 Introduction ................................ 17

PRECEDING PAGE BLANK NOT FILMED

F_NTENTIONALL¥ 8LAN5
iv

2.2 Governing Fluid Dynamics Equations .................. 17

2.3 Derivations ................................ 21

2.3.1 Velocity decomposition ...................... 21

2.3.2 New dependent variables ..................... 24

2.4 Derivation of Non-Primitive Variable Equations ............ 26

2.4.1 Two dimensions .......................... 26

2.4.2 Three dimensions ......................... 32

o PRELIMINARY ANALYSIS ...................... 35

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

3.2 Boundary Conditions ........................... 35

3.2.1 Scalar potential boundary conditions .............. 36

3.2.2 Vector potential boundary conditions .............. 38

3.2.3 Vorticity boundary conditions .................. 38

3.2.4 Dilatation boundary conditions ................. 46

3.2.5 Density boundary conditions ................... 48

3.2.6 Temperature boundary conditions ................ 49

3.2.7 Velocity boundary conditions .................. 49

3.3 Cartesian Grid Clustering ........................ 50

3.3.1 Clustering near boundaries .................... 53

3.3.2 Clustering near an interior point ................. 55

3.3.3 Arbitrary user defined clustering ................ 57

3.4 Poisson Equation Solvers ......................... 57

3.4.1 Vectorized point Gauss-Seidel with SOR ............ 58
V

,

3.5

3.5.1

3.5.2

3.5.3

3.4.2 AF scheme .............................

3.4.3 ADI scheme ............................

Poisson Solver Comparisons .......................

Scalar potential test problem ..................

Vector potential test problem ..................

Summary of Poisson solver experience .............

3.6 Time Marching ADI Solver for the Transport Equations ........

NUMERICAL RESULTS ........................

4.1 Introduction ................................

4.2 Solution Strategy .............................

4.3 Two-dimensional Cases ..........................

60

62

62

64

70

71

75

4.4

79

79

79

82

Incompressible flow ........................ 83

Compressible flow ......................... 121

Three-dimensional Cases ......................... 172

4.4.1 Incompressible flow ........................ 172

4.4.2 Compressible flow ......................... 180

5. DISCUSSION AND CONCLUSIONS ................ 181

6. RECOMMENDATIONS FOR FUTURE WORK ......... 185

7. REFERENCES ............................... 189

8. APPENDIX A: THREE-DIMENSIONAL VORTICITY AND

DILATATION TRANSPORT EQUATIONS ............ 197
vi

9. APPENDIX B" ALTERNATE NON-PRIMITIVE VARIABLE

EQUATIONS ................................ 203
vii

LIST OF TABLES

Table 3.1:

Table 3.2:

Table 3.3:

Table 3.4:

Table 4.1:

Table 4.2:

Poisson solver comparison on a 21 x 21 uniform grid ..... 66

Poisson solver comparison on a 41 × 41 uniform grid ..... 69

Poisson solver comparison on the scalar potential test problem 69

Poisson solver comparison on the vector potential test problem 70

Speed-up affect of limiting the vorticity solution domain . . . 115

Fully-developed centerline velocity and skin friction ..... 174
ix

LIST OF FIGURES

Figure 3.1:

Figure 3.2:

Figure 3.3:

Figure 3.4:

Figure 3.5:

Figure 3.6:

Figure 3.7:

Figure 3.8:

Figure 3.9:

Figure 3.10:

3-D Cartesian coordinate system with example boundary conditions for the potentials on a solid impermeable surface in the

x-z plane ............................. 39

An initial uniform velocity profile along a viscous boundary

for an impulsive start ...................... 40

2-D channel boundary conditions ................ 51

Typical 2-D channel grid. (a) Physical plane using _ = 1.05, a =

0.5,/_ = 1.2; (b) Computational plane ............. 55

Typical 2-D grid for an isolated airfoil at y = h. (a) Physical

plane using r = 5.0, c_ = 0, t¢ = 0,/3 = 1.5; (b) Computational

plane ............................... 56

Scalar potential test problem .................. 64

Convergence behavior of the Poisson solvers on the scalar potential test problem ....................... 67

Convergence history for the scalar potential test problem . . 68

Vector potential test problem .................. 71

Convergence behavior of the Poisson solvers on the vector potential test problem ....................... 72

PRECEDING PAGE BLANK NOT FILMED
X

Figure 3.11:

Figure 4.1:

Figure 4.2:

Figure 4.3:

Figure 4.4:

Figure 4.5:

Figure 4.6:

Figure 4.7:

Figure 4.8:

Figure 4.9:

Figure 4.10:

Figure 4.11:

Figure 4.12:

Convergence history for the vector potential test problem . . 73

Boundary conditions for laminar incompressible irrotational

flow over a bump ........................ 83

Parabolic arc and sine wave bumps .............. 87

Pressure coefficient for Moo = 0.6 flow over a parabolic arc

airfoil ............................... 88

Pressure coefficient for Moo = 0.6 flow over a sine wave arc

airfoil ............................... 89

Centerline velocity development for a 2-D channel inlet . . . 95

Skin-friction development along a 2-D channel ........ 96

Basic wall heat transfer boundary conditions ......... 97

Local Nusselt number for flow between parallel plates at constant temperature. Pe = 105. Asymptotic Nu = 7.541 .... 100

Local Nusselt number for flow between parallel plates at constant temperature. Pe = 210. Asymptotic Nu = 7.541 .... 101

Local Nusselt number for flow between parallel plates at constant temperature. Dual potential code with boundary-layer

assumptions ........................... 104

Profiles of the temperature parameter at various z* for the

constant wall temperature case ................. 105

Local Nusselt number for flow between parallel plates with

constant wall heat flux. Pe = 210. Asymptotic Nu = 8.235 106
xi

Figure 4.13:

Figure 4.14:

Figure 4.15:

Figure 4.16:

Figure 4.17:

Figure 4.18:

Figure 4.19:

Figure 4.20:

Figure 4.21:

Figure 4.22:

Profiles of the temperature parameter at various x* for the

constant wall heat flux case ................... 107

Local Nusselt number for flow between parallel plates with

mixed wall conditions. Pe = 210. Asymptotic Nu = 5.385

for one wall insulated and the other at constant heat flux.

Asymptotic Nu = 4.861 for one wall insulated and the other

at constant temperature ..................... 108

Profiles of the temperature parameter at various x* for the

case of one wall insulated and the other at constant heat flux 109

Profiles of the temperature parameter at various x* for the

case of one wall insulated and the other at constant temperaturell0

Boundary conditions for laminar incompressible flow over a

flat plate ............................. 112

Laminar incompressible boundary-layer velocity profile at Rex =

312,500 .............................. 117

Skin friction along a flat plate ................. 118

Velocity vectors for laminar flow over a flat plate with an irrotational inlet condition .................... 119

Velocity vectors for laminar flow over a flat plate with a rotational inlet condition ...................... 120

Boundary conditions for lanfinar compressible irrotational flow

over a bump ........................... 122
xii

Figure 4.23:

Figure 4.24:

Figure 4.25:

Figure 4.26:

Figure 4.27:

Figure 4.28:

Figure 4.29:

Figure 4.30:

Figure 4.31:

Figure 4.32:

Figure 4.33:

Figure 4.34:

Figure 4.35:

Pressure coefficient for Moo = 0.6 flow over a sine wave arc

airfoil with nlaximum thickness = 4% of chord ........ 129

Dilatation contours for Moo = 0.6 flow over a sine wave arc

airfoil ............................... 130

Density profiles from inlet to exit for the sine wave arc airfoil 131

Dilatation profiles from inlet to exit for the sine wave arc airfoil132

V velocity profiles from inlet to exit for the sine wave arc airfoil133

Comparison of the dilatation profiles for the stream function

and Euler mode of computation by the DP code. Sine wave

arc airfoil ............................. 134

Pressure coefficient for Moo = 0.6 flow over a parabolic arc

airfoil with maximum thickness = 6% of chord ........ 135

Dilatation contours for Moo = 0.6 flow over a parabolic arc

airfoil ............................... 136

Density profiles from inlet to exit for the parabolic arc airfoil 137

Dilatation profiles from inlet to exit for the parabolic arc airfoil138

V velocity profiles from inlet to exit for the parabolic arc airfoil139

Comparison of the dilatation profiles for the stream function

and Euler mode of computation by the DP code. Parabolic

arc airfoil ............................. 140

Thickening airfoil ........................ 141
Xlll

Figure 4.36:

Figure 4.37:

Figure 4.38:

Figure 4.39:

Figure 4.40:

Figure 4.41:

Figure 4.42:

Figure 4.43:

Figure 4.44:

Figure 4.45:

Figure 4.46:

Variation of the pressure coefficient with time at x/c = 0.525

for a parabolic arc airfoil. Moo = 0.785, maximum thickness

= 0.1 in t = 20 chord lengths .................. 143

Pressure coefficient for a thickening parabolic arc airfoil. Moo =

0.785, maximum thickness = 0.1 in t = 15 chord lengths . . . 144

2-D solution domain for divergence theorem application . . 149

Temperature contours for a variable property flow with constant wall temperature. Twall = 1.1, Re = 40 and M = 0.1 152

U velocity contours for a variable property flow with constant

wall temperature. Twall = 1.1, Re = 40 and M = 0.1 ..... 153

Dilatation contours for a variable property flow with constant

wall temperature. Twall = 1.1, Re = 40 and M = 0.1 ..... 154

Temperature contours for a variable property flow with constant wall temperature. Twall = 1.1, Re = 150 and M = 0.1 155

U velocity contours for a variable property flow with constant

wall temperature. Twall = 1.1, Re = 150 and M = 0.1 .... 156

Dilatation contours for a variable property flow with constant

wall temperature. Twall = 1.1, Re = 150 and M = 0.1 .... 157

Nusselt number comparison for constant wall temperature and

variable properties ........................ 158

U velocity contours near the inlet for a variable property flow

with constant wall temperature. Twall = 2.0, Re = 150 and

M -- 0.1 ............................. 159
x.iv

Figure 4.47:

Figure 4.48:

Figure 4.49:

Figure 4.50:

Figure 4.51:

Figure 4.52:

Figure 4.53:

Figure 4.54:

Figure 4.55:

Figure 4.56:

Figure 4.57:

Figure 4.58:

U velocity contours near the inlet for a variable property flow

with constant wall heat flux. Q = 25.0, Re = 150 and M = 0.1 161

Temperature contours for a variable property flow with constant wall heat flux. Q = 0.5, Re = 150 and M = 0.1 ..... 162

U velocity contours for a variable property flow with constant

wall heat flux. Q = 0.5, Re = 150 and M = 0.1 ........ 163

Dilatation contours for a variable property flow with constant

wall heat flux. Q = 0.5, Re = 150 and M = 0.1 ........ 164

Boundary conditions for laminar compressible flow over a flat

plate ............................... 165

Boundary layer profile at Rex = 100,500 and Moo = 0.5... 166

Temperature distribution for adiabatic flow over a flat plate

at Moo = 0.5 ........................... 167

Velocity profile in an adiabatic laminar boundary layer in supersonic flow ........................... 171

3-D duct geometry ........................ 172

Centerline velocity development for a 3-D square duct .... 175

Fully-developed velocity profile comparison for a square duct

at Re = 50. The Dual Potential (DP) code used a 29 × 29 × 60

uniform grid ........................... 176

Transverse velocity vectors at a plane near the exit of a square

duct with Re = 50 ........................ 177
XV

Figure 4.59: Transverse velocity vectors at a plane near the exit of a 0.50

aspect ratio duct with Re = 50 ................. 178
xvii

ABSTRACT

A dual potential formulation for numerically solving the Navier-Stokes equations

is developed and presented. The velocity field is decomposed using a scalar and

vector potential. Vorticity and dilatation are used as the dependent variables in the

momentum equations. Test cases in two dimensions verify the capability to solve

flows using approximations from potential flow to full Navier-Stokes simulations. A

three-dimensional incompressible flow formulation is also described.

An interesting feature of this approach to solving the Navier-Stokes equations is

the decomposition of the velocity field into a rotational part (vector potential) and

an irrotational part (scalar potential). The Helmholtz decomposition theorem allows

this splitting of the velocity field. This approach has had only limited use since it

increases the number of dependent variables in the solution. However, it has often

been used for incompressible flows where the solution scheme is known to be fast

and accurate. This research extends the usage of this method to fully compressible

Navier-Stokes simulations by using the dilatation variable along with vorticity.

A time-accurate, iterative algorithm is used for the uncoupled solution of the

governing equations. Several levels of flow approximation are available within the

framework of this method. Potential flow, Euler and full Navier-Stokes solutions are

possible using the dual potential formulation. Solution emciency can be enhanced

PRECEDkI_G PAGE BJ,2_/_ l',_OT FILMED

P_INTENTIONALL¥ BLAI_
XVlll

in a straightforward way. For some flows, the vorticity and/or dilatation may be

negligible in certain regions (e.g., far from a viscous boundary in an external flow).

It is possible to drop the calculation of these variables then and optimize the solution

speed. Also, efficient Poisson solvers are available for the potentials.

The relative merits of non-primitive variables versus primitive variables for solution of the Navier-Stokes equations are also discussed.
xix

ACKNOWLEDGMENTS

The authors gratefully acknowledge the support provided for this study by the

NASA/Ames Research Center through Interchange No. NCA2-174. We express special thanks to the Allison Gas Turbine Division of General Motors Corporation for

supporting the first author on a GM fellowship.
xxi

NOMENCLATURE

Roman Symbols

A

A

B

C

c/

cp

Cy

e

Et

F

h

h

I

i,j,k

i, j, k

_,j,k

k

L

M

n

Nu

P

P

Pr

vector potential

cross-sectional area

dilatation

airfoil chord length

skin-friction coefficient (1⁄2 poourW2)

pressure coefficient _,1⁄2Poou_ ]

specificheat at constant pressure

specificheat at constant volume

internalenergy

totalenergy per unit volume

body force

enthalpy

non-dimensional height of a 2-D grid

identity matrix

index in the x, y, z direction respectively

indices in tensor notation

unit vector in the x, y, z direction respectively

thermal conductivity

non-dimensional length of flowfield domain (e.g., channel length)

Mach number

outward unit normal

Nusselt number

exponent in the viscosity power law variation

pressure

Prandtl number (-_)

p_.CE_i_ FAGE ELANK NOT FtL,_D
xxii

Q

Q

R

Re

T

U

U

W

X

Y

Y

g

X

Y

Z

volume flow rate

finite-difference approximation to the exact solution

gas constant

Reynolds number

temperature

x component of velocity

average u velocity at an z cross-section

y component of velocity

z component of velocity

streamwise coordinate direction

spanwise coordinate direction

streamwise direction in the computational plane

spanwise direction in the computationM plane

normal coordinate direction

conservative body force in the z direction

conservative body force in the y direction

conservative body force in the z direction

Greek Symbols

o/

F

5ij

-g

62

52

P

Or

O

¢

IM

weighting parameter in hybrid differencing scheme

parameter in a 2-D stretching transformation

parameter in a 2-D stretching transformation

ratio of specific heats

rate of shear deformation in two dimensions, Uy + vz

Kronecker delta function

hybrid difference operator for convective terms

central difference operator, subscript gives direction to difference

e.g., 5xu = ui+ 1 - ui_ 1

central difference operator, subscript gives direction to difference

e.g., 52u = ui+ 1 - 2u i + ui_ 1

momentum thickness

parameter in a 2-D stretching transformation

density

parameter in a 2-D stretching transformation

temperature parameter

scalar potential

vector 3-D vorticity
**0

XXlll

60 2-D vorticity

Subscripts

hyd

m

n

r

t

t

W

X

y

Z

OO

based on the hydraulic diameter

mean value

normal direction

reference

time

tangential direction

wall

derivative in x direction

derivative in y direction

derivative in z direction

freestream value

Superscripts

k

n

non-dimensional quantity

provisional value

unit vector

iteration level

time level

Other Symbols

V

V 2

V

A

A

[]

gradient operator

Laplacian operator

backward difference operator, subscript gives direction to difference

e.g., Vmu = u i - ui_ 1

forward difference operator, subscript gives direction to difference

e.g.,/Xmu = ui+ 1 - u i

indicates an increment

used to enclose the units of measurement
1. INTRODUCTION

1.1 Background

The topic of this final report is a particular computational approach for solving

the Navier-Stokes equations. The Navier-Stokes equations are usually associated

with the field of fluid mechanics. Solutions to these equations with the appropriate

boundary conditions model fluid motion.

An analysis of fluid motion requires the solution for the physical laws of nature:

1. Conservation of mass

2. Newton's second law of motion

3. Conservation of energy

These laws can be formulated mathematically, with the help of some assumptions,

to become the Navier-Stokes equations. Formally, the Navier-Stokes equations refer

to the mathematical representation of Newton's second law. It will be more convenient for the purposes here to let the term "Navier-Stokes equations" include the

representation of all three physical laws above. Assumptions in the development are

that the coefficients of viscosity are related by a factor of -(2/3) according to Stokes'

hypothesis and that the fluid is Newtonian (Schlichting 1979). A Newtonian fluid
is one in which the fluid shear stress is linearly proportional to the rate of strain.

Additional relationships are included as necessary to describe certain processes or

fluids. Examples are the equation of state for a perfect gas, Sutherland's formula for

viscosity and Fourier's law of heat conduction.

As one can imagine, the successful solution to the Navier-Stokes equations can

help immensely in engineering design and optimization. The numerical solution of

the Navier-Stokes equations can be a complement to experimental and theoretical

fluid mechanics. Unfortunately, the Navier-Stokes equations are coupled and highly

non-linear. Only a few exact analytical solutions are available for simple conditions.

In most configurations of practical interest, numerical techniques must be used to

obtain a solution.

Much progress has been made in obtaining numerical solutions to the NavierStokes equations. Several mathematical formulations for the Navier-Stokes equations

have been developed. They can be divided into two classifications:

1. Primitive variable methods

2. Non-primitive variable methods

As the name suggests, primitive variable methods solve the Navier-Stokes equations using the primary variables as the unknowns. The primary variables are velocity,

total energy (or a variable related to the energy) and pressure or density. One way

to think of the primitive variables is that they are physical quantities which one can

measure in the laboratory. Non-primitive variables, on the other hand, are mathematically derived variables. They are derived from the primitive variables. The nonprimitive variables used in this report will replace the primary variable of velocity.
The replacements used will be vorticity and dilatation. One can devise techniques

to measure vorticity and dilatation experimentally, but the direct measurement of

such quantities is uncommon. Both primitive and non-primitive variable methods

have been used to obtain solutions to the Navier-Stokes equations by numerical techniques. Primitive variable methods are the most widely used for three-dimensional

simulations. Either primitive or non-primitive variables are used for two-dimensional

flow solutions although most applications of non-primitive variables have been for

incompressible flows. The following sections in this chapter will discuss some primitive and non-primitive variable solution methods. The focus of this study wiU be on

a particular non-primitive variable method that is extended to compressible viscous

flOW.

1.2 Primitive Variable Methods

Numerical methods of solving the Navier-Stokes equations using the physical

variables have attracted much attention. Several popular techniques will be mentioned here. The solution method depends on whether the flow is incompressible

or compressible, because the Navier-Stokes equations have a different mathematical

classification depending on the compressibility. For an unsteady incompressible flow,

the governing equations are elliptic/parabolic in time. For an unsteady compressible

flow, the equations are hyperbolic/parabolic in time.

The most common primitive variable solution method for incompressible flow

problems involves the use of a Poisson equation for pressure in place of the continuity

equation. An algorithm which employs this solution method is the SIMPLE (Semi-
4

Implicit Method for Pressure-Linked Equations) procedure (Patankar 1975, 1981).

Another primitive variable solution method is the artificial compressibility approach which modifies the continuity equation to include an unsteady term related to

pressure (Chorin 1967). The resulting equations are a mixed set of hyperbolic/parabolic

equations which can be solved using a time-dependent approach. This approach applies in two and three dimensions and can be modified to compute unsteady flows. An

available computer code that uses this method is INS3D (Kwak et al. 1986; Rogers

et al. 1987; Rogers and Kwak 1988).

A compressible flow solution is often obtained using a time-dependent or timedependent-like approach. Most schemes utilize implicit methods, for example, the

Beam and Warming (1978), Briley and McDonald (1977), or MacCormack (1981)

methods. An available code for these applications is F3D (Steger et al. 1986). Additional discussion and references on primitive variable solution methods can be found

in Hoist (1987).

1.3 Non-primitive Variable Methods

Methods which in some way replace the velocity with derived variables will be

discussed here. At the highest level of approximation, potential flows are typically

solved using either the velocity potential or the stream function. Examples of their

use are found in most fluid mechanics textbooks (Currie 1974). By definition, a

potential flow is irrotational so that the velocity field can be defined by the gradient of a scalar function. This scalar function is called the velocity potential. It is

analogous to the electric field potential. For an incompressible potential flow, the
only additional constraint is mass conservation. The equation of mass conservation

is a Laplace equation for the velocity potential which is easily solved. The velocity

potential is applicable in two and three dimensions. If the stream function is used

it is defined to satisfy mass conservation and the Laplacian of the stream function

satisfies the irrotationality condition. The stream function as defined here only exists

in two dimensions. In both of the above schemes, the momentum equation (vorticity

transport equation) is satisfied automatically since the vorticity is zero everywhere

for the assumption of irrotational flow.

The vorticity/stream function approach is widely used for solution of the twodimensional incompressible Navier-Stokes equations. This method is also discussed

in most fluid mechanics textbooks and is treated extensively in the book by Roache

(1972).

Many would consider this the limit of practicality for non-primitive variable

methods. However, there are at least two other noteworthy approaches to solving

the Navier-Stokes equations in non-primitive variables. Both are valid for two- and

three-dimensional unsteady flows. These methods are known as the

1. Vorticity/velocity approach

2. Vorticity/vector potential approach

These two schemes will be briefly described and then the focus will be placed on

the vorticity/vector potential method. The topic of this thesis will cover the vorticity/vector potential method. This method is also referred to by the the aliases

scalar/vector potential, vorticity/potential, or dual potential method. Since both
vorticity and dilatation are used in this work to replace the primitive variable momentum equations, it does not seem appropriate to identify the method as the vorticity/vector potential approach. Instead, the term dual potential will be used here

following Chaderjian and Steger (1985). This terminology identifies the method as

one which uses two potential functions in a velocity decomposition.

1.3.1 Vorticity/velocity approach

In this method for incompressible flow, the momentum equation is replaced by the

vorticity transport equation. Derivatives of the vorticity definition then yield Poisson

equations for the velocity when the continuity equation is used to make appropriate

substitutions. A more general derivation of the Poisson equations is to take the curl

of the vorticity and substitute in the vorticity definition from velocity. The identity

for the vector triple product then yields Poisson equations for the velocity. The

earliest use of this method was by Fasel (1976). He studied the stability of twodimensionM boundary layers using a coupled and iterative algorithm. Dennis et al.

(1979) used the vorticity/velocity method in the calculation of the cubical driven box

problem. Orlandi (1987) solved high Reynolds number flows over a backward facing

step. Other works using the vorticity/velocity formulation are Osswald et al. (1987),

Guj and Stella (1988), Gatski et al. (1982), Fasel and Booz (1984) and Farouk and

Fusegi (1985). There have been no reported compressible flow applications of this

method. However, the dilatation could be used as a dependent variable, as was done

in this research, to extend the vorticity/velocity method to compressible flow.
1.3.2 Dual potential approach

Derived variables which can be used to represent the three-dimensional continuity

and momentum equations for an incompressible flow are vorticity, a vector potential,

and a scalar potential. This is one possible three-dimensional extension of the more

familiar two-dimensional vorticity/stream function approach. This approach and others that use the vorticity as a dependent variable are appealing because vorticity is

generally located near boundaries in high Reynolds number flows and subsequently

diffused and convected away. For a three-dimensional incompressible flow the usual

procedure in the dual potential method is to solve the vorticity transport equation, a

vector Poisson equation and a scalar Poisson equation. These equations are derived

from the continuity and momentum equations where the velocity is defined as the

curl of a vector potential plus the gradient of a scalar potential. The existence of

these potentials is easily shown for an incompressible flow since the velocity field is

divergence free (Aziz and Hellums 1967).

There has been only one reported formulation of the dual potential method for

three-dimensional compressible, viscous, unsteady flows (Morino 1985). He derived

a set of equations for density, vorticity, entropy and the potentials. There have been

no reported calculations using Morino's formulation.

1.3.2.1 Applications of the dual potential method The dual potential

method has been applied to inviscid and viscous flow problems. Inviscid flow applications include the work of Rao et al. (1989) and Giannakoglou et al. (1988). Rao et al.

(1987) developed a three-dimensional inviscid rotational flow solver based on the dual
8

potential method. They incorporated a boundary layer interaction scheme for viscous

flow problems. Giannakoglou et al. (1988) compute two-dimensional steady rotational

transonic flows in arbitrarily shaped ducts and plane cascades. They decomposed the

mass flux vector into two potentials.

In the viscous regime, the dual potential method has been applied to problems

of three-dimensional natural convection in enclosures (Mallinson and De Vahl Davis

1973) and three-dimensional incompressible flows in ducts (Wong and Reizes 1984).

External viscous flows have been computed by Davis et al. (1989). No attempts

have been reported on the use of this method to solve three-dimensional unsteady

compressible viscous flows.

The dual potential method was first applied to natural convection problems by

Aziz and Hellums (1967). They used the dual potential method to transform the

Navier-Stokes equations. The transformed equations were solved using an alternating

direction implicit (ADI) scheme for the parabolic part of the problem (temperature

and vorticity transport equations) and a successive over-relaxation (SOR) method

for the elliptic portion (vector potential equations). They tested their technique by

applying it to the classical problem of convection in fluid layers bounded by solid

walls in both two and three dimensions.

Aziz and Hellums showed the dual potential method to be faster and more accurate than solutions obtained using the primitive variable approach. In fact, though

the equations are fewer in number for the primitive variable approach, Aziz and

Hellums report that they are much harder to solve than the equations in the dual

potential method. The difficulty arises from the highly non-linear nature of the pres-
9

sure equation and the coupling due to pressure in the momentum equations (as in

incompressible flow problems).

The technique developed above was used by Ozoe and co-workers (Ozoe et al.

1976, 1977, 1979, 1985) in solving a variety of natural convection problems. In the

1985 paper, the problem of three-dimensional turbulent natural convection in a cubical enclosure was solved using a two-equation model for turbulence.

Applications of the dual potential method to incompressible duct flow (throughflow) have not been wholly successful due to confusion over the appropriate vector

potential boundary conditions. The earhest work in this area was by Aregbesola and

Burley (1977). They presented a numerical finite-difference solution for the equations

of motion of a steady laminar incompressible flow in two and three dimensions using

the dual potential method. Wong and Reizes (1984) presented a dual potential formulation for unsteady incompressible flows in ducts of constant but arbitrary cross

section. They showed that the method is capable of handling flows over a wide range

of Reynolds numbers and imply that it can deal with flow situations in which other

models become inadequate. The dual potential method guarantees a zero divergence

of velocity while the usual primitive variable method can at best approximate global

continuity. That formulation was limited to simply connected domains. In a later

paper Wont and Reizes (1986) showed how to use the dual potential method to solve

for the three-dimensional flow in multiply connected regions such as annular geometries. Yang and Camarero (1986) used body fitted coordinates with the dual potential

method to simulate incompressible laminar flows in a square elbow and in a twisted

square elbow. The dual potential method in this paper is shown to be applicable to
10

general duct flow situations in simply connected regions. Hafez et al. (1987) used a

finite element method to solve the steady two-dimensional Navier-Stokes equations

in a dual potential formulation for subsonic and transonic flows. They computed

laminar and turbulent flow cases.

Some viscous external flow solutions were obtained by Rao (1987) and Davis et al.

(1986) for flow over two- and three-dimensional troughs. Rao (1987) used interacting

boundary layer theory to supply the vorticity to a dual potential code for the inviscid

rotational part of the flow. Davis et al. (1986, 1989) use a viscous dual potential

method for the entire flow field. Matching between the outer inviscid flow and the

inner viscous region is automatic in their case.

One possible extension of the dual potential method to three-dimensional compressible viscous flow has been formulated by Morino (1986). There are no reported

results in the compressible viscous regime.

Compressible viscous unsteady flows have been solved by a closely related method,

however. E1-Refaee et al. (1981) used a non-primitive variable method that replaced

the momentum equation with vorticity and dilatation transport equations. The velocity field was obtained from the vorticity and dilatation field by an integral representation. They solved compressible unsteady flows and demonstrated that the solution

field for vorticity and dilatation can easily be limited in their integral representation.

In this report a similar equation set is used, but the velocity is decomposed into two

potentials and solved completely by finite differences. The proposed extension of the

dual potential method to compressible flow would be directly applicable to the vorticity/velocity method. That is, the dilatation variable would be included to account
11

for compressibility.

In view of the short list of references on the dual potential method it is evident

that this method has not been widely implemented in computations. The main

reasons for this have been the need for the solution of several additional variables as

compared to the primitive variable approach and the inability to analyze the numerical

solution process for convergence. There are no good model problems to guide the way.

Three-dimensional flow solvers are computationaily demanding and the introduction of additional variables inevitably increases the computer memory requirement.

However, with the increasing memory of today's computers, and since Aziz and Hellums (1967) have shown that the dual potential approach can lead to faster and more

stable convergence than for primitive variable formulations (for certain problems),

the vector potential will perhaps play an increasing role in the solution of complex

three-dimensionai fluid dynamics problems (Wong and Reizes 1988). Certainly this

kind of formulation deserves continued investigation.

1.3.2.2 Advantages and disadvantages Relative advantages of the primitive variable method and dual potential method are cited in Morino (1985) and

Richardson and Cornish (1977). The major advantage of working in primitive variables is the relative simplicity of the equations and the fact that the primitive variables

have direct physical meaning.

The advantages of the dual potential method are:

1. The vorticity (and dilatation for compressible flow) need only be resolved in

distinct regions.
12

2. Continuity is automatically satisfied for incompressible flow.

3. The equations are weakly coupled (at least for inviscid flow).

4. Good numerical solution routines exist for Poisson equations.

5. Matching between an inviscid region and viscous region occurs automatically

because of the velocity decomposition into rotational and irrotational parts.

Disadvantages of the dual potential method applied to a three-dimensional compressible unsteady flow are:

1. The dual potential method involves ten dependent variables whereas the primitive variable method involves only five to represent conservation of mass, momentum and energy. (In two dimensions the number of dependent variables are

six for the dual potential method and four for the primitive variable method.)

2. The equations for the dual potential method are more complex than the equations associated with the primitive variables (or, they are simply unfamiliar).

3. The potentials do not have direct physical significance.

In addition to the natural disadvantages of the dual potential approach listed

above, there is a lack of available soKware as compared to primitive variable solution

methods. The extension of this approach to unsteady compressible viscous problems

is uncharted territory.
13

1.4 Scope of the Present Study

It has been the goal of this research to extend the capability of the dual potential method to compute unsteady compressible viscous flows. An algorithm has

been developed to provide two-dimensional full Navier-Stokes simulations. A threedimensional algorithm has been developed for incompressible flow only. Test cases

were computed to verify the ability to compute flow fields ranging from full potential

flow to flow fields requiring the full Navier-Stokes equations. It has been demonstrated

in this work that the calculation region can be limited for vorticity and dilatation,

thus providing a speed advantage for certain flows.

Several two-dimensional test cases will be presented to test various aspects of the

dual potential method. Both incompressible and compressible flows will be computed.

Incompressible flows will be studied for steady, irrotational, inviscid conditions

and for steady, rotational, viscous conditions. The steady, irrotational, inviscid test

case is that of flow over a biconvex airfoil (or a bump on a wall). Steady, rotational,

viscous conditions are simulated for a channel inlet and laminar boundary-layer case.

Heat transfer calculations will be made for the channel cases with constant wall

temperature and constant wall heat flux boundary conditions.

For compressible flow, steady and unsteady, irrotational, inviscid flows will be

computed and also steady, rotational, viscous flows. The irrotational, inviscid flows

are biconvex airfoil cases. The steady, rotational, viscous flows are channel inlet

and boundary-layer cases. The channel inlet flows are computed with constant wall

temperature and constant wall heat flux boundary conditions at a Mach number of

0.1. Calculations of the flow over a flat plate are made for a subsonic and supersonic
14

freestream.

In three dimensions only steady, viscous, incompressible channel inlet solutions

were obtained. A summary of the test cases to be presented is given below.

I. Two-Dimensional Cases

A. Incompressible flow

I. Steady irrotational inviscid flow

a) bump cases

2. Steady viscous flow

a) channel inlet with and without heat transfer

b) boundary layer

B. Compressible flow

1. Steady irrotational inviscid flow

a) bump cases

2. Unsteady irrotational inviscid flow

a) bump cases

3. Steady viscous flow

a) variable property channel flows

b) boundary layer

II. Three-Dimensional Cases

A. Incompressible flow

1. Steady viscous flow

a) channel inlet

B. Compressible flow

Progress in this research area has not been easy. There is very little guidance

in the literature on how to proceed with a full Navier-Stokes implementation of a

non-primitive variable method. The governing equations in non-primitive variable

form are unfamiliar. Non-linear terms were simply lumped into the source term and
15

the systemwassolveduncoupledin an iterative manner. As a result, it was necessary

to employ rather simple test cases to check various aspects of the formulation.

The work completed here is primarily in the development and evaluation of the

dual potential method as a flexible approach to solving the Navier-Stokes equations.

Several features of the method have been highlighted. For example, the solution domain for vorticity and/or dilatation may be limited to certain regions. Also demonstrated is the flexibility of the method to accommodate several approximations of

the full Navier-Stokes equations. This effort has advanced the understanding of the

dual potential method in the viscous compressible regime. It represents the first

application of this method to compute throughflow problems with heat transfer. Several basic problems are solved to check out aspects of the algorithm and computer

code. Only Cartesian grids are used for the test problems. Further evaluation and

optimization of the method reported herein are left for future work.

1.5 Organization

The main body of this report consists of six chapters and two appendices. The

presentation follows the logical development of the method from equation derivation

to boundary condition determination, grid generation, numerical algorithm selection

and, finally, flow simulations. In Chapter 2, the dual potential equations are derived

from the velocity decomposition and non-primitive variable dependent variables are

selected to represent the usual primitive variable (or pressure-velocity) form of the

Navier-Stokes equations. In Chapter 3, the numerical representations of the boundary

conditions are derived and the numerical algorithms are presented. The Cartesian
lfi

grid stretching is presented also. A comparison is made of the Poisson equation solvers

since the Poisson equation solution for the potentials can dominate the computation

time. In Chapter 4, the solution strategy is outlined and then numerical results are

reported for two- and three-dimensional test cases. The two-dimensional results represent cases from potential flow to situations requiring the full Navier-Stokes equations.

The three-dimensional results are for incompressible cases only, but are representative

of the speed of this approach for incompressible problems. Chapter 5 includes the

overall assessment of this method. Chapter 6 gives some incentives for future work

on the dual potential method.

The appendices contain equations for the full three-dimensional Navier-Stokes

implementation of this method. Also, alternative non-primitive variables are introduced which could be useful for some problems.
17

2. DUAL POTENTIAL FORMULATION

2.1 Introduction

In this chapter the mathematical equations which model fluid flow are presented.

These equations are the Navier-Stokes equations. They state the conservation of mass,

momentum and energy for a Newtonian, Stokesian fluid. The usual form of these

equations has the primitive variables (p, V, Et) as the primary unknowns. Using the

Hehnholtz decomposition theorem, the velocity field can be split into a rotational

part and an irrotational part. Each part is represented by a potential function.

This decomposition yields a non-primitive variable formulation for the Navier-Stokes

equations.

2.2 Governing Fluid Dynamics Equations

The following equations apply to a continuum fluid.

The conservation of mass is stated

Dp -_ --*

D--'_+ pV. V = 0 (2.1)

The conservation of momentum (Newton's second law), with the assumptions
18

that Stokes' hypothesis holds and that the fluid is Newtonian, is written

(2.2)

Where F = Xi + Yj + Zk is the body force and 6ij is the Kronecker delta function:

6iJ= { 1 ifi=j 0 ifi¢j

In the energy equation only internal and kinetic energy will be considered important. The conservation of energy is then written

De ---+ --_ 0Q _ .-_ + ¢ (2.3) p._+pV. V- Ot

where represents heat energy production by external agencies, q is the heat

conduction and _ is dissipation. Fourier's law of heat conduction will be assumed so

..._4

q =-k VT

The dissipation function for a Newtonian fluid in a Cartesian coordinate system beCOlnes

:

(2.4)

The ideal gas equation of state and a viscosity law are used to close the system for

laminar flow. Constant specific heats are assumed throughout. Reference conditions

are selected to non-dimensionalize the equations. Reference quantities will be denoted

by the subscript r. Fluid properties for air will be used in the calculations.
19

For an ideal (or perfect) gas the following relationshipsexist:

p = pnT (2.5)

Cp

e = cvT h = cpT "_ = _ R = Cp - Cv (2.6)

CV

Sutherland's law of viscosity is used in the form

T3/2

I-t = C1 T + C2 (2.7)

where C 1 and C 2 are constants for a particular gas. For air at moderate temperatures

(approximately 200K-1000K), C 1 = 1.458 × 10 -6 [kg/(m sx/K)] and C 2 = 110.4K.

Power law variations for # were also available in the computer code:

- 0.5 < p < 1.0 (2.8)

A constant Prandtl number is assumed and thermal conductivity is obtained

from the definition, Pr = -_. Typical values of the fixed quantities chosen for air

are:

_, = 1.4 (2.9)

R = 287[_] (2.10)

Pr -- #rCpr - 0.7 (2.11)

kr

A reference length will be designated by Lr [m] and a reference velocity by Ur [n_A].

The reference length is taken to be a characteristic length of the problem such as the

hydraulic diameter for internal flow or chord length for external flow. The reference

velocity is taken to be the magnitude of either the inlet velocity for internal flows or

the freestream velocity for external flows.
2O

The following reference conditions were set for incompressible flow calculations

with the Reynolds number and reference length specified:

Tr = 288.15K (59°F)

#r = C1T :[,,] r + C2

_ W

kr p.rCp Pr [m---_]

1.

Ur - pr Lr

u,.

Mr -- _

Pr = pr RTr

Reference conditions for compressible flow calculations with the Reynolds number, Mach number and reference length specified are as follows:

Tr = 288.15K (59°F)

T3/2 [kg] gr = C1Tr + C2

?r - UrLr [,,] -J
21

The following non-dimensional variables are then obtained:

, x , y , z t*

Lr Lr Lr

t

(Lr/Ur)

, u , v , w , tt

u -Ur v -Ur w -Ur # gr

*_ P * P T* T ,

P Pr p -- -- e

e

R* --

R , Cv , Cp k* k

(U2/Tr) c v _ (U2/Tr) Cp = (Ur2--]-Tr) = k--_

where the variables distinguished by an asterisk are non-dimensional. The nondimensional variables will be used throughout, so the asterisk will be dropped in

the following. The non-dimensional gas constant above is equivalent to

1

R3,M 2

so that

p = pRT

in dimensional or non-dimensionM variables.

2.3 Derivations

2.3.1 Velocity decomposition

The basis of the dual potential method is a splitting of the velocity field into

rotational and irrotational parts. In this section, the impetus for splitting the velocity
22

field in this way is presented. The Helmholtz decomposition theorem formally permits

the splitting.

A useful classification of vector fields is possible using the divergence and curl

operators (Ames 1977). For V • E = 0 at every point of a region R, the vector field

---k

E is said to be solenoidal or divergence free. Physically this means that there are

--+

no sources or sinks in R. If V × E = 0 at every point in R, the field is said to be

irrotational. The following classification of vector fields can then be made:

Class I

Class II

Class III

Class IV

Solenoidal and irrotational:

xE=O .E=O

Irrotational but not solenoidal:

--4

V×E=O V.E_O

Solenoidal but not irrotational:

xE 0 V.E=O

Neither solenoidal nor irrotationah

×E#0 V.E#0

An important theorem in vector field theory is called the Helmholtz decomposition theorem. It states that any vector field can be split into a curl free and divergence

free part. Using the above classifications, it can be observed that the velocity field

of an incompressible fluid is in Class III and the velocity field of a rotational compressible flow is in Class IV. Applying the Helmholtz decomposition theorem to the

velocity vector one obtains:

V=V Cq- V × A (2.12)

----4

It is obvious that the curl free part of the velocity is V ¢. The divergence free

part (recall from the above that another word for divergence free is solenoidal) is the
23

-----+

vector field A. The vector A has been assumed to be solenoidal by design. This is

an arbitrary but appropriate choice to fit the Helmholtz decomposition theorem and

to remove the redundancy of describing a three component vector (V) using another

--->

three component vector (A) plus the gradient of a scalar (¢).

For the classification of the velocity field then, one can compute the divergence

and curl of Equation 2.12 to obtain:

V.V = V2¢=B (2.13)

--+ _ _ _ -+ V2 _ Vx V = w = V (V. A)- (2.14)

.--.+

where ¢ is the scalar potential, B is the dilatation or rate of volumetric strain, w

is the vorticity and A is the vector potential. The vector potential, A, is chosen to

be divergence free. The Laplacian operator in Equation 2.14 is the vector Laplacian.

Throughout this report only rectangular coordinates are used so each component of

the vector Laplacian is similar to a scalar Laplacian.

The vector potential and vorticity will be represented in three dimensions as

follows:

A = AI_+ A2j+A3k (2.15)

w = Wl{ +w2J+w3k (2.16)

The components of the vorticity are obtained from:

w=VxV

w l=wy-vz, w 2=-wz+uz, w 3=vx-uy
24

The velocity components are then:

/

V= v = ey + Alz A3x

w ez + A2 x Aly

(2.17)

----+ _..4

In two dimensions only one component of A and w exist. For the standard

two-dimensional geometry shown in Figure 3.3 the single components are A 3 and

w 3. To simplify things in two dimensions, the subscripts will be dropped on the

vector potential and vorticity so that A and w refer to the two-dimensional case. The

velocity components in two dimensions are then:

u ¢z + Ay I

v= = (2.1s)

v ev - Ax

2.3.2 New dependent variables

In the well known incompressible application of the above decomposition, the

momentum equations become the vorticity transport equation. Continuity is satisfied

--+ .-.+

by the solution of a Laplace equation for ¢, since V. V = 0 in Equation 2.13. Finally,

the potentials are used to compute the velocity field. Any other governing equations

remain unchanged (energy equation, equation of state, etc.).

For a compressible flow, however, V. V # 0. In this case, V. V = V2¢ =

B # 0 in Equation 2.13. An additional equation is required to give the dilatation,

B, for the solution of the scalar potential. By counting the number of equations and

unknowns, one can see that the three-dimensional momentum equations represent

three equations with three unknowns (or two equations and two unknowns in two
25

dimensions). The vorticity transport equations formed by taking the curl of the threedimensional momentum equations yield just two independent vorticity component

equations (see Appendix A). In two dimensions, the vorticity is a single scalar quantity

so again there is one remaining usage of the momentum equations permitted. This

means that an equation (hopefully for the dilatation, B) can be derived by some

operation on the momentum equations in either two or three dimensions. Looking at

possible operations on the momentum equations, one choice is to take the divergence

of the momentum equation.

The divergence of the momentum equation yields for possible dependent variables

either the pressure (as in the pressure Poisson equation) or the divergence of the

---+ --4.

velocity, V • V. The divergence of velocity will be represented by the scalar variable,

B, known as the dilatation or rate of volumetric strain.

Another possible combination of the momentum equations gives a scalar variable

which is the rate of shear deformation (or shear strain rate). Let the shear strain

rate be represented by the symbol, F. In two dimensions, F = uy + vx. This

dependent variable is formed by taking _ of the V momentum equation + _0v of

the z momentum equation. The wall shear stress is simply #F. The variable set

.....4

of F and _v can form the basis for an interesting computational procedure in two

dimensions. Unfortunately, one usually hopes to compute the skin friction, not give

it as a boundary condition. However, this could be a useful inverse type calculation

procedure (see Appendix B).
26

2.4 Derivation of Non-Primitive Variable Equations

The dual potential equations will be derived starting from the governing equations written in primitive variable form. For purposes of comparison, the number of

unknowns required in a flow simulation using either primitive variables or the dual

potentials will be given. The number of unknowns is computed by considering the

continuity, momentum and energy equations only.

2.4.1 Two dimensions

The derivation of the dual potential method in two dimensions will contain the

fewest simplifying assumptions except, of course, that it only considers two space

dimensions. The solution capability will be for flows that require the full NavierStokes equations (i.e., unsteady compressible viscous flow).

In the primitive variables, the continuity, momentum and energy equations represent four equations for four unknowns. These four unknowns may be p, u, v and T.

These variables are solved using:

• continuity:

pt + (p=)=+ (p_)y= o (2.19)

• x momentum:

_v._ )]+_ . _+_

• y momentum:

PD-'t- = Y - _ -_ Re1 Oyo[ # (2Ov _, Oy

2-_. +---- # -4- (2.21)

3 Re 0¢
27

• energy:

--+

DT --+ _ 7 _____P___ _

Dt - (7- 1)(V. V)T + pRePr V. V T + pRecv (2.22)

The pressure may be obtained from the ideal gas law, p = pRT. Any additional

variables (#,k,7, Pr , etc.) must be accompanied by their own equation of state or

constitutive equation.

For the non-primitive variables of the dual potential method, the same twodimensional flow requires the solution for 4, A, w, B, T and p. The potentials give the

velocity field. The potentials themselves are determined by the solution of Poisson

equations derived from the divergence and curl of the velocity (cf. Equations 2.13 and

2.14).

To see what is needed for the dual potential formulation, consider the primitive

variable equations above. They are the governing equations for fluid flow, but now

it is desired to solve not for velocities directly but rather for the potentials. The

velocities are subsequently determined from the potentials by solving Equation 2.18.

The momentum equations must be recast to generate a solution to be used by the

potentials. Equations 2.13 and 2.14 for the velocity field splitting suggest that the

divergence and curl of velocity be sought as dependent variables. These are obtainable

from the x and y momentum equations above by taking the divergence and curl of

the prinfitive variable momentum equation.

The curl of the velocity (vorticity) is obtained as a dependent variable by taking

the curl of the momentum equation. The group, vz - uy = w, is retained as the

dependent variable. The divergence of the velocity is obtained as a dependent variable

by taking _0 of the z momentum equation and _y of the y momentum equation and
28

summing. The grouping, uz + vy = B, can be retained as the dependent variable.

The remaining governing equations can be left unchanged from the primitive variable

equation set.

The continuity equation can be used to compute the density. The energy equation

can be used to solve for T or enthalpy, or any other variable that is related to energy.

The ideal gas equation of state can then be used to compute the pressure if it is

needed.

To summarize: The six variables of the dual potential method corresponding to

the solution of the continuity, momentum and energy equations are determined as

described below:

1. The continuity equation is used to compute the density, p.

2. The curl of the momentum equations gives the vorticity transport equation

for w. Conservative body forces are eliminated by this operation. Also, for

incompressible flow, the pressure is eliminated. For compressible flow, however,

the pressure derivatives remain but can be expressed in terms of other variables

by using the equation of state.

3. The divergence of the momentum equations gives the dilatation transport equation to be solved for B.

4. The energy equation is solved for T, or enthalpy, or a related variable.

5. A Poisson equation is solved for _ with B as the source term.

6. A Poisson equation is solved for A with -w as the source term.
29

The equations for a two-dimensional dual potential formulation will be presented

below and subsequently solved in non-conservative form. Body forces will be neglected.

In two dimensions the dual potential representation of the compressible NavierStokes equations for constant specific heats are:

pt+upz+vpy+pB = 0 (2.23)

wt + uwz + VWy - (p_e) V2w = SI(p,B,w,T ) (2.24)

B t+uBx+vBy- (4.) _-_ V2B

T t + uTx + vTy - pRePr V2T

= S2(p,B,w,T ) (2.25)

= S3(p,B,w,T ) (2.26)

V2¢ = B (2.27)

V2A = -,_ (2.28)

The source terms for the w, B and T transport equations are given below. For

the vorticity transport equation, S 1 contains the compressibility and c 1 contains the

variable viscosity terms.

4 B 4

" [p_(-_y+ _ _)- p_(_ +_By)] S 1 = -Bw+ p'P'_Re

R (2.29)

C l

+ (p_Ty- pyT_)+ Re
3O

c 1 =

(2.30)

For the dilatation transport equation, S 2 contains the "compressibility" and c2

contains the variable viscosity terms.

$2 4 B 4

_nV2T nTv2p-pR (T_p_+Typy)+-_tp_+p_)+ RT. 2 Rec2 (2.31)

c2 2[._y(_y +._)+.=-_+,,yy_y] __7v2B-2.

P

+ #x--[_2Wy + 8 B 3p p - -_ +_)]

+#Y [2wx+_ sBy+ 2-PyB- 2PYvy- P_(u ]

(2.32)

The energy equation source term, $3, contains the compressibility and variable

thermal conductivity.

$3

= -(7- I)BT + pRec----"_

7

pRe Pr [k_T_ + kyTy]

+ w 2 + 4(v_cuy - VyU_c)]

(2.33)

The solution strategy for this system of equations is outlined in Section 4.2.

Briefly, the equations are grouped into an "incompressible" and compressible part.

An incompressible solution is obtained by computing among the equations in the
31

"incompressible" set of equations given below. A compressible solution requires a

pass through both the "incompressible" and compressible sets of equations. In that

case the "incompressible" equations actually contain terms representing compressible

effects. This grouping is used to allow an incompressible solution to be the starting

solution for a compressible problem.

vorticity transport equation

Incompressible V2¢ = 0

V2A = -w

dilatation transport equation

energy equation

continuity equation

Compressible

V2¢= B

ideal gas law,p = pRT

property updates : #, k

The transport equations above are solved using an ADI scheme. Source terms

of the dependent variable (undifferentiated) are treated implicitly. Derivatives of the

dependent variable in the source term are treated explicitly so as not to weight the

off-diagonals. This did provide a slightly faster solution than treating all dependent

variable source terms implicitly.

One will immediately notice the many derivatives introduced in the governing

equations by this method. Even though most of the test cases to be presented are for

subsonic flows, it was necessary to handle some terms conservatively. In particular,

the pressure related terms in the dilatation transport equation and the conduction
32

terms in the energy equation are best handled conservatively for the heat transfer

test cases computed. By treating these terms conservatively it is easier to obtain an

accurate representation than by expanding out in the chain rule form and differencing.

2.4.2 Three dimensions

Only the incompressible equations will be derived here. The three-dimensional

compressible set is given in Appendix A.

2.4.2.1 Incompressible flow The governing equations for viscous incompressible flow in non-dimensional vector form are:

• continuity:

V. V = 0 (2.34)

• momentum:

O--7-+(V. OV _ _V) --*V=- _V P+ ___ V2V-, +F4 (2.35)

There are four unknowns to be determined for the primitive variable solution of a

three-dimensional incompressible flow. Usually these unknowns are u, v, w and p. It

will be seen that the dual potential method requires the solution for seven variables

to satisfy the continuity and momentum equations.

By taking the curl of the momentum equation the pressure is eliminated yielding

the vorticity transport equation:

"-4 _ ---+

0 o_ --, --, (_ _ _ 1 V2 _ +V × F (2.36) o--v+(v, v) • v) v _ Re
33

.--+

If F is a conservative body force (such as gravity) then it too is eliminated by the

curl.

Decomposing the velocity vector according to the Helmholtz decomposition theorem as in Equation 2.12 one obtains:

---+ _ ---+ .-.+

V=V ¢ + V × A (2.37)

-_+ _

Since V • (V × A) = 0, the substitution of Equation 2.37 into Equation 2.34

leads to:

v2¢=0

-.-+

The relation between the vector potential, A, and vorticity is obtained by taking

the curl of Equation 2.37. This becomes:

V2 _ --_ _ -----4' A V (V A)

The vector potential is chosen to be solenoidal so that the above reduces to:

.....+

V 2 A= ----+03

Therefore, for viscous incompressible flow, the momentum and mass conservation

equations and the vorticity definition can be solved using:

-.-+

003ot +(v._ _v) --, 03- • v) v= 1 v2_ ' +v x F (2.3S)

V2¢ = 0 (2.39)

.--+

.--+

V 2 A = -03 (2.40)
34

Then, the velocity is decoded from the potentials according to the Helmholtz

decomposition:

V =V_+ V x A

If the pressure field is needed it can be obtained by solving for pressure from

one of the primitive variable momentum equations or by solving the pressure Poisson

equation. It is also possible to solve for the pressure by computing a force balance on

an appropriate fluid element since the velocity field and hence the shear stress field

is already determined by the solution strategy above.
35

3. PRELIMINARY ANALYSIS

3.1 Introduction

In this chapter the essential parts for assembling a dual potential code are gathered together. First the boundary conditions for the dependent variables are presented. Next, the Cartesian stretched grids are described and, finally, the necessary

solvers are explained.

3.2 Boundary Conditions

A general presentation of the boundary conditions will be given here. The boundary conditions are applicable for two- or three-dimensional problems. Boundaries have

been classified as one of the following:

1. Solid -- impermeable boundary (slip or no slip)

2. Throughflow boundary -- boundary crossed by the streamwise velocity

3. Far-field boundary -- a freestream boundary which may be modeled as an

impermeable boundary, a porous boundary or some other freestream condition

according to the problem.
36

The boundary conditions are, of course, also identified with a partial differential

equation. For ease of presentation, only the boundary types as listed above are

discussed here. To determine what conditions are imposed in a particular problem

it is necessary to know the classification of the governing partial differential equation

(hyperbolic, parabolic or elliptic) and the number and type of derivatives of the

dependent variable. Specific boundary conditions for each model problem will be

given in figures in the discussion of the results. The solid boundary and throughflow

boundary conditions for the potentials are thoroughly derived in Wong and Reizes

(1984) and Hirasaki and Hellums (1970). For the transport variables ( _, B,p,T),

a fully-developed exit condition is specified by dropping second order streamwise

derivatives (except for density) and upwinding other streamwise derivatives at the

exit. Example boundary conditions for a two-dimensional channel throughflow case

are illustrated in Figure 3.3 for reference.

3.2.1 Scalar potential boundary conditions

The scalar potential is obtained from the solution of the Poisson equation, V2qb =

B. This is an elliptic equation so a condition on q_ or its derivative must be given on

all boundaries.

Since the velocity is decomposed into two potentials, it is useful to ascribe certain

of the velocity boundary conditions to each potential. It has already been demonstrated by Hirasaki and Hellums (1970) that if the scalar potential were used to deal

with possible throughflow velocities, then simple boundary conditions on the vector
37

potential are possible. Following Hirasaki and Hellums,

---4

0._____= _ . V

On

on the boundaries. Thus, the scalar potential has Neumann conditions all around.

The normal derivative of the scalar potential is the normal outflow velocity at the

boundary. For a pure Neumann problem such as this, existence and uniqueness of the

solution are concerns. Existence is ensured if the compatibility condition given by

Green's theorem is satisfied. Uniqueness is ensured by prescribing the scalar potential

at some point. For the test cases to be discussed later, it was possible to make one

boundary a Dirichlet boundary and still satisfy the Neumann boundary condition by

virtue of Green's theorem. The positive n direction in the following is the outward

normal,

3.2.1.1 Solid boundaries

on the scalar potential is

For no flow through the boundary, the condition

3.2.1.2 Throughflow boundaries An inlet or exit is the best example of a

throughflow boundary. The condition on ¢ is then

0¢

n

On streamwise velocity, say u i

where,

1//4

u i - Ai V. _ ds

and A i is the cross-sectional area of the throughflow boundary.
38

3.2.1.3 Far-field boundaries A derivative condition on the scalar potential

is also used at a far-field boundary as follows:

0_¢_= normal component of velocity through the far-field boundary, usually O.

On

3.2.2 Vector potential boundary conditions

The vector potential is obtained from the solution of the Poisson equation,

V 2 A= - w. This is an elliptic equation so a condition on the components of

-..4

A or its derivatives must be given on all boundaries.

With the above choice of scalar potential boundary conditions, the boundary

----b

conditions on A for a simply connected region may be shown to be

OAn

At- On -0 (3.1)

where the subscripts t and n denote the tangential and normal components of A

respectively. An example of the potential boundary conditions for a solid surface in

the x-z plane are shown in Figure 3.1.

3.2.3 Vorticity boundary conditions

Vorticity is only needed for rotational flow computations. It is generated, for

example, at no-slip boundaries and diffused and convected away. Vorticity may also

be specified as part of the inlet or initial conditions. Vorticity can also be generated

by shocks, but such flows will not be considered here. The transport equation for

vorticity states the conservation of vorticity. Therefore, the boundary conditions are

extremely important in defining the flow field.
39

/

Z

y

t

A 3 = A 1

=0

I

X

Figure 3.1: 3-D Cartesian coordinate system with example boundary conditions for

the potentials on a solid impermeable surface in the z-z plane
4O

AY 1

y,j index

I

!

|

|

I

2- 7 ---2

\\-,_ ,_ ,%,,\\_ _\\', ,_ ,-x,\\\\\

z, i index

Figure 3.2: An initial uniform velocity profile along a viscous boundary for an impulsive start

3.2.3.1 Solid boundaries The vorticity at a solid boundary is obtained

from the no-slip conditions. Several different formulations are possible. The wall

vorticity can be computed from the velocity derivatives or from a Taylor series expansion of the vector potential. Consider a wall in the (z, z) plane. Using the vorticity

definition from velocities, the wall vorticity (at y = 0) is

w 1 = wy (3.2)

w 2 = 0 (3.3)

w3 = -Uy (3.4)

This method was used by Aziz and Hellums (1967). It was only used in this work

to compute an initial wall vorticity based on an impulsive start. For example, the

initial wall vorticity for a two-dimensional impulsive start is given by w = -uy (see

Figure 3.2. The finite-difference initial wall vorticity is then:
41

u(i,2)

.,(i,1) = y(i,2)- y(i,1)

For subsequent global loops, a better approach (Roache 1972) was used to obtain

the wall vorticity. This better approach uses the vector potential. Consider a vorticity

producing surface at y = 0 (j = 1) in a two-dimensional flow. The grid may be

stretched, so let the transverse spacing to the first node point be Ay I as in Figure 3.2.

For the i index in the z direction and j index in the y direction, write the Taylor

series expansion for the vector potential at (i, 2) about the wall values as follows:

A(i,2) Ayl+ (AYl + (AYl) 3

A( i, 1) + -_y i,1 20y 2 i,1 -6 _ i,1

+ O( Ayl)4 (3.5)

From the velocity decomposition in two dimensions (Equation 2.18), u = Cx + Ay.

At the wall this is

u(i,1) = _i,1 _Y i,1

but, u(/, 1)is zero for a no-slip boundary. Therefore, _yAli,1

can be replaced by the scalar potential derivative along the boundary.

in the Taylor series above

(3.6)

Again using the velocity decomposition in two dimensions, u = Cz + Ay, the y

derivative of this equation gives: uy = Cxy + Ayy. At the solid impermeable wall

this is:

_yyli,102¢i,1 02A+ - OxOy _y2 i,1
42

_A

From the scalar potential boundary condition at a solid boundary, _Y i,1 =

O. Hence,

the second order y derivative of A at the wall is simply:

02A

The wall vorticity in two dimensions is the component w 3 = w = vz - Uy. Along the

wallv=0, sothatw(i,1)= -_gyOui,1. Therefore, the wall vorticity is introduced into

the Taylor series for the Ayy term as:

02A = _yy (gu I = -w(i, 1) (3.7) i,1 i,1

Equation 3.7 is also obtainable from the Poisson equation for the vector potential,

V2A = -w. At the wall, A = 0, so Azz = 0 and Ayy = -w. Substituting for Ay

and Ayy in the Taylor series of Equation 3.5 one obtains:

0¢ I AYl+2(-w(i'l))( Ayl)2+6-_Y 1 03A3 i,1 A(i,2) = A(i, 1)- _-_x i,1

w(i, 1) -

becomes:

2 [A(i,l) -A(i,2)-o_ " l)¢(i+ l'l)-¢(i-I I) ] ( Ayl)2 O-_x('' 7 ' Ay 1 (3.10)

Solving for w(i, 1) yields the following first order approximation for the wall vorticity

¢0(i, 1)- 2[A(i, 1)-A(i,2)-¢zi,IAYl]+O(AYl) (3.9)

(AYl) 2

The vector potentialin two dimensions iszero at a solidimpermeable wall,so A(i, 1) =

0 in the above. Using a second order centraldifferencefor ez with possiblestretching

in the x direction (Section 3.3), the finite-differenceformula for the wall vorticity
43

In a similar fashion the Taylor series and above substitutions could be carried out

in the computational plane ($,_). The unequal grid spacing is accounted for by the

metrics to yield the following formula:

w(i, 1) = 2 _ (i, 1) A(i, 1)- A(i,2)

'

The above first order form was used most frequently in the results to be reported

here. A second order formulation also based on the Taylor series is obtainable from

Equation 3.5 by retaining the (/_yl) 3 term. The second order boundary conditions

will be derived here for a three-dimensional ease. Consider the x-z plane again. The

boundary condition on the w 1 component will be derived first by starting with the

Taylor series expansion of the vector potential component .41 . The Taylor series

expansion for .41 at a mesh point adjacent to the boundary, (y = /_Yl at j = 2), is

given by

oA 1 102Al l (/,yl) 2

Al(i ,2,k) = Al(i ,1,k)+ Oy i,l,k/_Yl+ 2 0y2 i,l,k

+60y103A13 i,l,k ( Ayl)3 +O( Ayl)4 (3.12)

where the indices i, j, k denote the x, y, z positions, respectively.

___+

The first term of the expansion is zero since the tangential components of A

vanish on an impermeable boundary. Hence

Al(i,l,k)=O (3.13)
44

To evaluate the second term, use the velocity component which contains Aly in

its decomposition. From Equation 2.17 note that

w = Cz + A2 z - A1 v

Since w = 0 on the boundary, the derivative in the second term of the expansion can

be written

OA 1 0¢ OA 2 (3.14)

From Equations 2.14 and 3.1, the second order derivative in the third term may easily

be identified as

02A1 l" = -wl(i,l,k ) (3.15)

OY2 i,l,k

Here is where the second order method departs from the first order approach given

earlier. To obtain an expression for the third order derivative of A1, a linear distribution of vorticity and Alyy over the first mesh interval is assumed (Wong and Reizes

1984), so that

Wl

02A1

Oy _

Y

wl(i,l,k ) + --_yl [wl(i,2, k) -wl(i,l,k)] (3.16)

02AI't 0y-_( ", l, k) + _ Y[O20_yl (i,2, k)- 0---_y 02Ali ( ,1, k) ] (3.17)

Combining the above linear distributions in order to write w 1 as a function of Alyy

yields

_1 = ---_y202A l + _ Y [_1(i,2,k)+ 02A1" -_-_-y2(,,2,k)]

Differentiating Equation 3.18 with respect to y yields

Owl 03A1 + l(i,2, k) + 2, k) 0y- 0y - -ff -y2c,

(3.18)

(3.19)
45

Solving Equation 3.19for Alyyy at y = 0 gives

03Al(i,l,k) = wl(i,l,k)-wl(i,2, k)

Oy 3 Ay 1

+-- wl(i,2,k)+ (i,2,k) (3.20)

Ay I (Jy- J

Finally, substituting Equations 3.13, 3.14, 3.15 and 3.20 into Equation 3.12 gives

a second order relationship between the boundary vorticity and the potentials, that

is_

wl(i, 1,k ) = 3Al(i,2, k ) 3 [¢z(i,l,k) + A2x(i,l,k)] (Ayl)2 +

1

+ _Alyy(,,2, k ) + O(Ayl) 2

(3.21)

Similarly, the z component of vorticity at this boundary is found to be:

=

3A3(i,2, k)

( Ayl)2 + _1 A2z(_,l,k" ) + -_A3yy(,,2, 1 . k ) + O( Ay 1 )2 (3.22)

The third and final vorticity component at this boundary is computed from the

no-slip condition and the definition of vorticity to be

w2(i, 1, k)=0 (3.23)

as already stated at the beginning of this discussion. In actual use, the second order

vorticity boundary condition increased the CPU time with no noticeable improvement in accuracy. Roache (1972) also reports that second order vorticity boundary

conditions can be less stable and less capable than the first order boundary conditions.

3.2.3.2 Throu_hflow boundaries The only throughflow boundaries used

here are the inlet plane and exit plane for two- and three-dimensional problems. The

inlet flow field may be specified as either rotational or irrotational. This is controlled
46

by the inlet plane transverse velocity. For two-dimensional flow, the inlet v velocity

is controlled to provide rotational or irrotational flow as follows:

rotational ; set v = 0

inlet condition (3.24)

irrotational ; compute v from the decomposition

The inlet plane vorticity is then computed from the definition, _ = V × V.

For a two-dimensional rotational inlet condition, the inlet transverse velocity is

set to zero and used throughout the code. The irrotational inlet condition is achieved

by computing the inlet plane transverse velocity from the decomposition and using

that non-zero inlet transverse velocity throughout the code. The streamwise velocity

component is unchanged since it fixes the inlet mass flux.

At the exit the second order streamwise derivative is dropped and other streamwise derivatives are upwinded.

3.2.3.3 Far-field boundaries The far-field vorticity must be specified or

determined from the specified upstream flow conditions. It is not safe to extrapolate

the vorticity to a far-field boundary, or any boundary since that does not account for

the physics of vorticity production, convection and diffusion. Far enough away from

a surface generating vorticity, the vorticity should be zero. If the velocity is known

at the far-field, vorticity may be obtained from the definition.

3.2.4 Dilatation boundary conditions

The dilatation, B, is only needed for a compressible flow. Boundary conditions

are developed from the expected velocity boundary conditions and the definition of
47

the dilatation,

B= V. V =uz+vy+wz

or by applying the continuity equation,

Pt + upx + vpy + wpz + pB = 0

The treatment of this variable is the most uncertain. There is little guidance

in the literature. E1-Refaee (1981) used the dilatation in his non-primitive variable

formulation. He extrapolated to obtain the dilatation at the boundary. The boundary

values of dilatation were relaxed and set to zero as the solution approached steady

state. El-Reface solved external flow problems only.

3.2.4.1 Solid boundaries For a viscous flow along an impermeable wall, all

velocities are zero. This reduces the continuity equation to

Pt + pB = 0

For steady boundary conditions on density, or at steady state it is clear that

B=0

Otherwise, the density time derivative or the velocity derivatives must be evaluated

to compute

B - Ptp wall ' from continuity (3.25)

or,

B = uz + Vy + Wz , from the definition of B. (3.26)

In any case, B = 0 at the boundary at steady state.
48

3.2.4.2 Throulghflow boundaries The only throughflow boundaries used

here are the inlet plane and exit plane for two-dimensional compressible flow problems.

The inlet flow may be specified as either rotational or irrotational as stated in the

section on vorticity boundary conditions. From the definition of dilatation, B =

ux + vy, an inlet Dirichlet condition can be computed. For the rotational flow inlet

condition (which corresponds to a uniform parallel inlet flow), v = 0 at the inlet.

This gives B = ux for the inlet boundary condition. For the irrotational flow inlet

condition, v _ 0 and in general v y _ 0. This gives/3 = ux +vy for the inlet boundary

condition. It has been found that an inlet boundary condition of B = 0 may be used

for either inlet condition. By setting B = 0 at the inlet one avoids the ambiguity of

specifying the dilatation inlet condition using velocities which in turn depend on B.

At the exit, the second order streamwise derivative is dropped and other streamwise

derivatives are upwinded.

3.2.4.3 Far-feld boundaries Generally, the velocity field is unchanging at

the far field, so B is given by the unchanging velocity derivatives. A typical far-field

condition is uniform parallel flow so that B = 0.

3.2.5 Density boundary conditions

The density boundary conditions are set using one of the following approaches,

depending on the problem:
49

1. By applying the boundary-layer assumption for viscous compressible flow at a

solid boundary and the ideal gas law,

py = 0

so,

p(i,1) = p(i,2)T(i,2)

T(i, 1)

2. By using the Bernoulli equation for inviscid irrotational compressible flows,

1

p=[1+7_2 1M2(l_u2_v2_w2)}_-':-l--

3. By application of the continuity equation written along the boundary,

Pt + upx + vpy + pB = O

3.2.6 Temperature boundary conditions

The temperature boundary conditions are either set as Dirichlet conditions or

derivative conditions based on a prescribed wall heat flux. First through fourth order

polynomial derivative conditions are included in the computer code as options. The

inlet temperature field is user specified. At the exit plane, second order streamwise

derivatives are dropped and other streamwise derivatives are upwinded.

3.2.7 Velocity boundary conditions

The velocity is not a primary variable in this method. Velocity boundary conditions are dictated by the flow physics and are used to develop boundary conditions on
50

the derived variables, such as the dilatation boundary conditions above. The conditions to use are straightforward--e.g., no-slip at solid boundaries for a viscous fluid.

It is important to point out that the velocity boundary conditions are used, but not

necessarily enforced by the dual potential formulation. For example, a no-slip condition is used to develop boundary conditions for the vorticity at a solid viscous wall.

However, there is a small slip velocity computed by the velocity decomposition. The

slip velocity goes to zero as the grid is refined. To see how the slip velocity arises, consider the velocity decomposition for a two-dimensional flow over a flat plate oriented

as in Figure 3.2. The velocity decomposition for this case is given by Equation 2.18,

repeated here for convenience.

V= -- (3.27)

v Cy - Az

At the flat plate surface, the boundary conditions on the potentials are:

Cy=O

A = 0

It is obvious that the v component of velocity will be zero both analytically and

numerically (provided the same difference formula is used to compute v as was used

to enforce the boundary condition on ¢). However, the u component will only be zero

if ¢;r = -Ay. This is not exactly satisfied numerically.

3.3 Cartesian Grid Clustering

Simple independent variable transformations are used to allow for stretching of

the two-dimensional Cartesian grids (Anderson et al. 1984). The stretching trans-
51

y, j index

inlet:

W---_ O_ or

w = _X -- Uy

B = O, or

B = ux + Vy

T=I

p=l

A=O

Cz = u(1,j)

///////////////////////////////////////

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

w = f(A,¢)

B=O

walls: T = Tw or _ given

p = fixed by py = 0 and T distribution

A=O

Cy = 0

exit:

O_XX _ 0

Bzz = 0

Tzz = 0

pz upwinded

Az = 0

z, i index

Figure 3.3: 2-D channel boundary conditions
52

formations may be readily applied to a third coordinate direction for usein a threedimensionalCartesian grid. The transformations used here are from the family of

generalstretching transformationsproposedby Roberts (1971).

The coordinate transformation mapsthe physicalplane(x, y) into the computational plane (_, _):

where, _ = f(x) only and y = f(y) only. Also, the transformations are scaled so that

/_ = /X_" = 1. The computational grid coordinate values then correspond to the

grid indices like (_, y) = (i - 1,j - 1). This simplifies the coding and avoids divisions

in the numerical algorithm for a slight speedup advantage. In the following, NI is the

largest z index and NJ is the largest y index in the domain. NI and NJ correspond

to the maximum dimensions of the physical grid L, h (see Figure 3.4).

Applying these transformations to the governing fluid dynamic equations requires

the following partial derivatives: (These have been simplified since _ and _ are only

functions of the respective coordinate directions.)

8 OF 8

8z Oz 0_

0 O_ 0

m

Oy Oy O_

82 I 0_ _ 2 02 _ c32 -__ c3
53

0 2 O_Oy 0 2

OzOy Ox Oy O_Og

The above metric terms are computed numerically using second order accurate

differences for the unequal mesh spacing. Analytical calculation of the metrics are

possible for the stretchings to be presented. It was found, however, that the numerical

calculation provided better flow results.

3.3.1 Clusterinl$ near boundaries

Take the z direction to be streamwise. The transverse direction is then y with

"walls" at y = 0 and y = h. The following transformation allows packing near the

inlet and near one or both walls.

Y

= NI(1- In {[¢r + 1 - (z/L)]/[a - 1 + (z/L)]}) 1 < a < o¢ (3.28)

In [(o" + 1) /(tr - 1)]

ct+

c01n ({/_ + [y(2_ + 1)/h]- 2c_} / {/3- [y(2_ + 1)/h] + 2a})(3.29) (1

In [(fl -4- 1) / (/_ - I)1

(a + 1) - (a'- 1) [_1] (1-_-1)

= (3.30)

x L (1___1)+ 1

is readily found as:

This is designed so that 0 < _ < (NI- 1) and 0 < y < (NJ - 1) for 0 <: z < L,

0 < y < h with /_= /k_ = 1. Equation 3.29 forypacks near y = h for a = 0

and near both walls equally for a = 0.5. The inverse of equations 3.28 and 3.29 are

needed to construct the physical grid (x, y). The inverse for the above transformation
54

y = h(- + sgn(t¢)

ioot

-/3 + 2c_)

)

1--0_

(3.31)

where,

1 for _>0

sgn(t¢) =

-1 for t¢ < 0

The inverse for y has been modified using t¢ to direct the y clustering to either

wall as described below. The stretching parameters in the above transformation have

the following affects:

x direction clusterin6

1 < _r < c¢ -- The stretching parameter o" clusters more points near x = 0 as a

approaches 1. The grid becomes more uniform as z _ c¢. An

essentially uniform z grid is generated for _r = 10.

y direction clusterin$

1 < fl < oo -- The stretching parameter/_ controls the y direction clustering

(spacing ratio). Maximum clustering is achieved as/3 approaches 1.

An essentially uniform y grid is generated for _ = 10.

The stretching parameter a is either 0 or 0.5.

_=0

clusters near y = h only.

clusters near y = 0 only.

o_ = 0.5 Cluster points near y = 0 and y = h equally, t¢ makes no difference

when ot = 0.5.

An example of the above clustering is shown in Figure 3.4.
55

]

I

i

I 1

L

= X

(NJ-1) ,, '",'I,

II ' II III

II II III I

II III

, l I Ill II

II III I

II I III

I I I I

,, I ,, I"l i II I III I

0 (NI-1)

(a) (h)

Figure 3.4: Typical 2-D channel grid. (a) Physical plane using _r = 1.05, a =

0.5,fl = 1.2; (b) Computational plane

3.3.2 Clusterin$ near an interior point

This clustering technique was used for the bump problem (isolated airfoil). It

can also be used to cluster points about an obstacle located within the grid. The

equations given will work in either the x or y directions and are designed to cluster

near a single point or to symmetrically cluster about an object by reflecting the

generated grid about the line of symmetry.

The z coordinate scheme for the bump problem will be given here. This requires

an odd number of i points and assumes that the bump or airfoil is always placed in

the center of the x grid. (The bump in the test cases is simulated by the blowing

condition rather than occupying x, y space.)

As stated previously, ¥ is simply given by the grid point index, i - 1, so that
56

Zc L

Ca) (b)

Figure 3.5: Typical 2-D grid for an isolated airfoil at y = h. (a) Physical plane using

r = 5.0, a = 0, t¢ = 0, fl = 1.5; (b) Computational plane

/X_ = 1. The inverse transformation is:

_ -W }

sinh[ )] (332,

x = xc 1+ sinh(rW)

where

1

W - In

2r

0 < r < oo (3.33)

The stretching parameter, r varies from zero (no stretching) to large values which

produce the most refinement near x = Xc. An example of this stretching is shown in

Figure 3.5 where the grid has been refined near the line zc and reflected about the

line of symmetry at x = L/2.
57

3.3.3 Arbitrary user defined clustering

The dual potential computer program will accept any user defined clustering

so long as _ = f(x) only and _ = f(y) only. These may be input as arbitrary

x,g points or as equations relating x,g to x,g as in the above examples. The code

automatically scales the metrics to form a two-dimensionaJ computational space with

A_ = A_ = 1.

3.4 Poisson Equation Solvers

Three different methods have been used to solve the Poisson equations for the

potentials. One method is the vectorized point Gauss-Seidel scheme with successive

over-relaxation (SOB.). The other two are alternating direction implicit (ADI) type

schemes (Mitchell and Grifllths 1980). The two ADI type schemes will be distinguished as follows:

1. A scheme formed by factoring and then splitting into a two step formula similar

to the D'Yakonov (1963) approach. This will be referred to as the approximately

factored (AF) scheme in this discussion.

2. A Peaceman-Rachford type scheme with a half-time step level (Peaceman and

Rachford 1955). This scheme will be referred to hereafter as the ADI scheme.

In summary, the three schemes to be used to solve the Poisson equations are:

1. Vectorized point Gauss-Seidel with SOR

2. AF scheme
58

3. ADI scheme

The model equation V2q _ = S will be used to demonstrate the above algorithms.

For the AF and ADI schemes, the Poisson equation is written with a fictitious time

derivative: q_t + V2q_ = S. The iterations represent "time" levels with time step h.

The time step may also be thought of as a relaxation parameter. At convergence the

time derivative term is negligible.

For a Cartesian grid with possible stretching according to the coordinate transformations in Section 3.3, the equation to discretize is:

The computational plane grid spacing is unity so/k_ and/k_ do not appear in

the algorithms below.

3.4.1 Vectorized point Gauss-Seidel with SOR

The finite-difference form of the equation to solve is:

Since Gauss-Seidel is a point iterative method, the exact application of the above

algorithm will depend on the mesh point ordering. The Gauss-Seidel method is based

on immediate use of the most recent values. Therefore, the _most recent used in

the above equation is either _k4-1 or _k. Solving the above equation iteratively

by points will not vectorize due to data dependency. This can be illustrated by a

simple example. Consider the five-point formula finite-difference scheme for the twodimensional Laplace equation, V2U = 0, on a uniform Cartesian grid with Dirichlet
59

boundary conditions:

Ui+ 1,j - 2Ui,j + Ui- 1,j + Ui,j+1 - 2Ui,j + Ui,j- 1

(/X_)2 ( Ay)2

As a further simplification take Az =

would code the following:

4o

-0

Ay = 1 so that solving this by points one

DO 40 J=2,N3-1

DO 40 I=2,NI-1

U(I,J) = 0.25 * (U(I-I,J) + U(I+I,J) + U(I,J-I)

CONTINUE

+ U(I,J+I))

The dependency of U (I, J) on U (I- 1, J) inhibits vectorization. Notice, however,

that for a fixed J, the even I indexed values of U on the left hand side depend only

on the odd I indexed ones on the right hand side. A vectorization strategy is now

apparent. The data dependency is removed from the computation by "coloring" the

grid as a checkerboard and updating the U in two sweeps. At a starting J, the odd

I index points may be thought of as black squares on a checkerboard and the even

I index points may be thought of as red squares. At the next J, the odd I points

are then red squares and the even I points are black squares. This red-black coloring

continues in J until the grid is patterned like a checkerboard. In one sweep the black

points are updated using only red points and in another sweep the red points are

updated using only black points. This is easily implemented by incrementing the I

loop by 2. Some initial work is required to determine the starting and ending I indices

for each J. Therefore, in two sweeps the solution is iteratively updated and the code

will vectorize. The compiler, however, will not recognize that the data dependency

has been removed. The programmer must direct the compiler to vectorize the loops.

This strategy is coded below. Note that the starting and stopping I indices are a
60

function of J. It is a simple alternating function, such as 2,3,2,3,... for IBSTRT(J) and

then 3,2,3,2,... for IRSTRT(J).

C...BLACK POINTS

DO 40 J=2,NJ-I

C...COMPILER DIRECTIVE: IGNORE VECTOR DEPENDENCIES IN THE I LOOP

DO 40 I=IBSTRT(J),IBEND(J),2

U(I,J) = 0.25 * (U(I-I,J) + U(I+I,J) + U(I,J-1) + U(I,J+I))

40 CONTINUE

C...RED POINTS

DO 41J=2,NJ-1

C...COMPILER DIRECTIVE: IGNORE VECTOR DEPENDENCIES IN THE I LOOP

DO 411=IRSTRT(J),IREND(J),2

U(I,J) = 0.25 * (U(I-i,J) + U(I+l,J) + U(I,J-l) + U(I,J+l))

41 CONTINUE

This may be combined with SOR for a further speed advantage. The extension of this vectorization concept to three dimensions is straightforward. The threedimensional problem may be solved as a stack of two-dimensional checkerboards, as

a four-color point method, or by extending the idea of colored points to colored lines

and solving by lines rather than points. The red-black strategy and other vectorizable

structures are discussed in Gentzsch and Neves (1987).

3.4.2 AF scheme

Using first order temporal differencing and second order spatial differencing on

the model Poisson equation one obtains:

(3.34)
61

It is convenient to use the "delta" form of the dependent variable so define

Ak_ = ffk+l ffk

To form a "delta" in the difference equation above subtract the spatial derivatives of

ffk from each side. Multiplying through by the time step (or relaxation parameter),

h, gives

The "delta" form allows easy implicit handling of the boundary conditions.

Steady Dirichlet conditions are automatically satisfied by the fact that Akq ' = 0

on the boundary. Steady Neumann conditions are easily handled by reflecting the

Akff at the boundary. The derivative function cancels so there is nothing to be added

to the right hand side. The actual Dirichlet values or derivatives are input by the

source term spatial derivatives of q,k.

An approximately factored form of the equation above is:

where the factors can be denoted as L 1 and L 2 so that the AF representation is

LIL 2 A k ,I_ = RHS.

As it is written, the algorithm is implemented by sweeping implicitly in the z

direction then in the y direction. The solution (Akq_) is attained in the two steps:
62

Step 1:

Step 2:

L 1 TEMP = RHS

L 2 A k (I, = TEMP

3.4.3 ADI scheme

The starting form of the difference equation is the same as for the AF scheme

(Equation 3.34). The ADI scheme splits the calculation for ff into two steps. In the

first step, the z derivatives are treated implicitly at the iteration level k + 1⁄2, while

the y derivative terms are lagged. The provisional solution is denoted as ffk+1⁄2.

The second step obtains the solution, _k+l by discretizing the y derivative terms

implicitly at the iteration level k + 1 and using the x derivative terms at the k + 1⁄2

level.

[1+h/0 \22 k+1⁄2hs j [I- h(Oy 262h[O2y - ul k

[

I +

h_(OY_22 -_\-_y,] _5_+ h(02y_ 2i_0y2) _]q_kd-. 2 J t,3 1 = 2Si,J h k +[

I - g_zz)_- h[0-£\22 7_-_z2)-Yj h{O2"Z_'_lk+ i,j 1⁄2

3.5 Poisson Solver Comparisons

An assessment of the two-dimensional Poisson solvers was made on test problems for the scalar and vector potentials. The test problems are from incompressible

channel throughflow cases. The three solvers were tested on stretched and uniform

grids with various aspect ratios. The L 2 norm of the error is used to compare the
63

convergence history for each method:

I E(ok+I _ ok)2 L2 = number of points

However, the requirement used for a converged solution is

Iok+l - Ckl < e (3.35)

qk+l -

max

where t} represents either potential and k is the iteration level. The tolerance used for

the comparisons below is _ = 10 -6. This rather strict tolerance can yield a decrease

in the L 2 error by as much as 6 orders of magnitude. It was found that the scalar

potential had to be obtained to this accuracy for reliable flow solutions.

All of the Poisson solvers were coded with the boundary conditions incorporated

in the solution algorithm. The relaxation parameter used for the vectorized SOR

method will be designated by w. The range for w is 0 < o; < 2. The relaxation

parameter used for the AF and ADI schemes will be designated by h, the fictitious

time step. The time derivative term which was added to the Poisson equation simulates a parabolic problem, Og + V20 = S. The way this is written, marching is

only permitted in the negative "time" direction. Therefore, a negative h was used

to march the AF and ADI solvers. The need for a negative h is also evident in the

numerical representation of these schemes since a negative h will add to the diagonal term of the coe_cient matrix. The AF and ADI schemes were solved using a

vectorized tridiagonal solver.
64

¢=0

Y ey = 0

v2¢=0

ey = 0

Figure 3.6: Scalar potential test problem

3.5.1 Scalar potential test problem

The scalar potential is given by the Poisson equation, V2¢ = B. For an incompressible flow the dilatation, B, is zero. The test problem for incompressible

flow through a channel is then, V2¢ = 0, with boundary conditions as shown in

Figure 3.6. The exact solution for this case is ¢ = x (notice that the problem is

actually one-dimensional). This is a difficult problem to solve numerically due to the

many Neumann boundary conditions. For incompressible problems, whether steady

or unsteady, the equation for the scalar potential can be solved once and for all. An

efficient solver may not seem important for the incompressible case. For a compressible flow, however, the dilatation field will be computed by a time marching method

so the scalar potential will have to be solved as often as every global iteration. The

cases to be studied in this report are mostly subsonic and the dilatation may be expected to be small. Therefore, the incompressible solution may be used as a starting

solution for a compressible problem. Also, in some cases the dilatation may be so
65

small as to be treated as a perturbation of the solution to the incompressible problem. For these reasons and the fact that the exact solution is easily obtained for the

incompressible case, the incompressible problem will be used to assess the solvers for

the scalar potential.

The first comparison on the scalar potential test problem is with uniform grids.

The aspect ratio, defined as a - _, is varied for each case. The results for a 21 × 21

grid are shown in Table 3.1 and for a 41 × 41 grid in Table 3.2. The SOR method

has trouble with this problem and requires an unusually high optimum relaxation

parameter, w. It can be seen from these results that the SOR method cannot compete

with ADI or AF when the grid aspect ratio is far from 1. Stretched grids will be

necessary in the solution of viscous problems so this immediately excludes the SOR

method for use in solving for the scalar potential, especially when B is non-zero.

Notice that the AF and ADI schemes solve the problem in the same number of

iterations regardless of the cell aspect ratio. Also, the optimum relaxation parameter

for AF and ADI can be reasonably predicted from the results shown in Tables 3.1

and 3.2 for 0 < a < 1. This is the most likely range of a for channel type viscous

flow geometries. The optimum h for the AF scheme behaves like

The optimum h for the ADI scheme behaves like

where ha 1 denotes the optimum h for a uniform grid with a = 1. This can be used to

get a reasonable estimate for the optimum relaxation parameter to use for a stretched

grid.
66

Table 3.1: Poisson solver comparison on a 21 x 21 uniform grid

Uniform 21 × 21 grid

method

SOR

AF

ADI

iterations to convergence

a =110 0.1 0.01

109 611 7253

3 3 3

138 138 138

CPU time (s)

1.0 0.1 0.01 MFLOP

0.0634 0.349 4.11 39.7

0.001645 0.001626 0.001667 54.1

0.1234 0.1265 0.1236 64.6

Optimum relaxation parameters (w or h)

method

SOR

AF

ADI

a =1.0

1.898

-160

-0.0089

0.1

1.9854

-10240

-0.89

0.01

1.9991

-1310720

-89

Testing these solvers for the scalar potential on a realistic grid with stretching gives the convergence behavior shown in Figure 3.7. The grid is stretched as

shown in Figure 3.4. The cell aspect ratios for this stretched grid range from cr =

0.4044-0.01265. Notice that the AF scheme converges over a very wide range of the

relaxation parameter, h. In contrast, the SOR method has a very limited range of

relaxation parameter which gives convergence in the thousands of iterations at the

very fastest. The extremely good performance of the AF scheme on this test problem

is incidental. It is explained by the fact that the exact solution for this test case

is ¢ = x, and the factorization error introduced in the AF scheme contains cross

derivative terms. The cross derivatives and hence the factorization error therefore go

to zero quickly. The optimum conditions for the solvers on this problem are shown

in Table 3.3. The convergence history at the optimum conditions for the solvers is

shown in Figure 3.8.
67

COE©>

__90

gOC0

(I)

CO

._)0

[3_

(])

_]2

.p

q_00

._)

(I)

._)

(_

(])

+.) o0

[_

C]

cC__0C

O)

03LC0

C..D

-0

O_

_JE0

t-O

-0

fl)

t-OgmEo

t_

cl

.o

-d

...._

Q. g0

i@

e@ @,EJ @ _, @1

c- w__CI

Q0 _CI O_

U0 L_..,nJ

_.Iu-_

[] o,,_I

OC .......

._D

o6o

o"

d

!I

Jill I II I I Jill I I I I I Jill I I I I

_0[ _0[ _0I

eoue6de^uoo oq SUOR1Oda17

0

CDo

r

mer0

i_o

•_.0

0 (-

- .,..)o

-m• ._ L _ 0 ,._

"oa-_

g-0

-0

0

r

(23
68

Convergence HLstor4

-o

scoLor potentLot test problem

o

C3 ,

v--4-

i

0

CD

1

1

I

I \

\

O

L

[.

• "

10_

I

T

0

C3

J

C3

10 °

I

t

1

' /

I

t

!

' i

I

I

I

I

' /

I

Figure 3.8:

LEGEND

SOR scheme

RF scheme

ADI scheme

I I I I Ilil I I I I I III/ I I ! ! I iill I

101 102 I0_

LterotLon number

Convergence history for the scalar potential test problem
69

Table 3.2: Poisson solver comparison on a 41 x 41 uniform grid

Uniform 41 × 41 grid

method

SOR

AF

ADI

iterations to convergence

=1.0 0.1 0.01

195 1260 8105

3 3 3

258 258 258

CPU time (s)

1.0 0.I 0.01 MFLOP

0.306 2.01 12.7 55.3

0.00468 0.00475 0.00483 72.8

0.719 0.713 0.701 82.1

Optimum relaxation parameters (w or h)

method c_ =1.0 0.1 0.01

SOR 1.9482 1.9925 1.99932

AF -160 -10240 -1310720

ADI -0.0045 -0.45 -45

With stretching, AF and ADI outperform the vectorized Gauss-Seidel easily on

this problem. With non-zero dilatation and more complicated geometry the AF

scheme is not expected to display such a tremendous advantage over ADI as in the

example here. Also note the high relaxation parameter required for the fastest convergence by the SOR method.

Table 3.3: Poisson solver comparison on the scalar potential test problem

41 × 41 stretched channel grid

method:

SOR

AF

ADI

optimum

relaxation param.

1.99644

-80000

-3.0

iterations to

convergence

1954

142

CPU time (s)

3.33

0.00509

0.4059

MFLOP

53.7

70.9

80.8
7O

Table 3.4: Poisson solver comparison on the vector potential test problem

41 × 41 stretched channel grid

optimum iterations to

method: relaxation param, convergence CPU time (s) MFLOP

SOR 1.721 39 0.0629 54.1

AF -0.003 35 0.0531 70.1

ADI -0.00071 105 0.2143 82.2

3.5.2 Vector potential test problem

The vector potential test problem has been taken from an incompressible channel

throughflow case. The equation and boundary conditions are shown in Figure 3.9.

For these conditions, the SOR method did not exhibit extreme sensitivity to the grid

aspect ratio and compared favorably with the AF and ADI methods. The Poisson

equation for the vector potential, with three of the four boundaries having Dirichlet

conditions, is easily handled by SOR using optimum relaxation parameters in the

approximate range 1.2 < ,, < 1.8. This equation must be solved as often as every

global iteration for a rotational flow.

The vector potential equation was solved here on a 41 × 41 grid stretched as

shown in Figure 3.4. A converged incompressible channel flow at Re = 300 provided

the vorticity source term. The behavior of the solvers over a range of relaxation

parameters is given in Figure 3.10. The optimum conditions are given in Table 3.4

and the L 2 error for the optimum convergence of each solver is shown in Figure 3.11.

In this case the AF and SOR schemes compare favorably in CPU time. In an actual

flow solution, however, the ADI scheme has achieved convergence faster than AF or

SOR when the vector potential solver is called at each global iteration.
71

A=0

Y A = 0

A=O

Ax=O

X

Figure 3.9: Vector potential test problem

3.5.3 Summary of Poisson solver experience

This was not meant to be an exhaustive study of methods for solving Poisson's

equation numerically. Two very specific Poisson equation problems were studied

to assess the numerical solvers. The results for the scalar potential test problem

are particularly perplexing. The point SOR method performs poorly while the AF

scheme is unbelievably fast. This is at least partly due to the one-dimensional nature

of the scalar potential test problem. Point SOR cannot take advantage of the onedimensional nature of the solution because it solves the domain pointwise. In the first

iteration, point SOR creates a two-dimensional distribution. This happens because

the solver cannot sense the information from opposite boundaries simultaneously.

This is unlike the AF and ADI schemes which can immediately detect the influence

of opposing boundaries in the implicit direction. The AF and ADI performance on the

scalar potential test problem could be dependent on the sweep order, especially since

the solution is one-dimensional. This was not tested, however. The biggest advantage
?2

L

>

-->

0 "L.

C •

o E

....)

0 -o

O)

O_ r

U

--_ c

,--,

o ,_

0 .,

----) E

..-) 0

L c_

4_) ..,->

u

0

L

L" .-_

0 o..

C

-,->

L o

>

C

0

I I I I I I I I I I I I I

Oi _0[

eoue_jeAuoa oq suoT_oJe_ 7

I,)

0

el

° _,,,I

0

o

o

0

0

o

0

0

0

c_

,,"4
73

Convergence HL stor _1

vector potentLot Lest problem

0

LESEND

SOR scheme

RF scheme

RD] scheme

T

C)

i

o

\

\

\

\

\

\

\

\

\

\

\

Figure 3.11:

I I I 1

25.0 50.0 75.0 lO0.O 125.0

LterotLon number

Convergence history for the vector potential test problem
74

for the AF scheme is that the factorization error is fortuitously zero because of the

one-dimensional solution to the scalar potential test problem.

The vector potential problem doesn't offer the anomalies that the scalar potential

problem does. Point SOR and the AF scheme are very competitive. The poorer

performance of the ADI scheme is perhaps due to the fact that it is not written in

delta form.

The apparent overall best solver for the Poisson equations in the dual potential

method is the AF scheme. At the optimum convergence condition it converges in

the fewest iterations and in the least CPU time. The AF scheme converges rapidly

over a large range of relaxation parameter. This is important in getting started with

a solution for a new problem. Almost any reasonable relaxation parameter, h, will

work whereas some effort is required to get a fast solution by the SOR method. The

ADI scheme is a consistent performer for the very different problems tested, but its

optimum is not as fast as AF at its optimum for the conditions tested above. However,

in actual flow calculations where the solver may be called at each global iteration the

ADI scheme is always competitive. SOR is competitive when the boundary conditions

are Dirichlet. It has a narrow band of relaxation parameter giving fast convergence

as compared to either AF or ADI. This narrow band and steep slope (see Figures 3.7

and 3.10) is undesirable for predicting a good relaxation parameter to use on a new

problem. One rule of thumb, however, is that the optimum relaxation parameter for

SOR increases as the number of grid points increase on the same problem.

The above results and discussion suggest that the scalar potential Poisson equation is best handled by AF or ADI. The vector potential problem is best handled
75

by AF or SOR, though the ADI scheme was found to be competitive on actual flow

problems. In actual use, the AF scheme was the easiest to code, followed by the ADI

scheme and then the vectorized point SOR. All three schemes were fastest when the

boundary condition information was incorporated into the algorithm. Not indicated

in the simple test cases above is the fact that SOR loses any advantage on large grids.

The computational effort for SOR is greater than for the AF or ADI schemes. Beyond

some grid size then, the SOR scheme should be dropped in favor of the more efficient

AF or ADI schemes. Pointwise iterative schemes tend to bias the solution depending

on the sweep direction. The effect of sweep direction was not studied for the SOR

scheme above.

For the three-dimensional (incompressible) problems the vector potential equation was solved using point SOR. The exact solution of the scalar potential was used

so that the Poisson equation was not solved numerically.

3.6 Time Marching ADI Solver for the Transport Equations

The ADI scheme presented here is based on the one proposed by Douglas (1962).

The transport equations for -_, B, p and T in three dimensions may be written as

( o2s vN Oz2 + + +cs-e

.....+

where S represents any of the dependent variables o.,, B, p or T and l{'-eerepresents

the coefficient of the diffusion terms. The term 0 includes all remaining terms. The

source term cS may be lumped into 0, but here it will be treated implicitly.

The above is written in a uniform Cartesian grid. A stretching transformation can

easily be introduced, but would only clutter the development here. The stretching
76

transformations were used in Section 3.4 on the Poissonequations to provide an

example of their use. Here, the complete three-dimensional algorithm on a uniform

grid will be sufficiently complicated without metric terms included.

In the development of the ADI scheme for the transport equations, the following

definitions will be used

x i = i Ax, yj = j Ay, z k = k Az (3.36)

t n = n At (3.37)

_z u + lul V¢ u- lul (3.38)

_ = _uy-_+(1-_) 2 A_ + 2

_y _y +lvlVy v-lvl

= aV2Ay+(l_a) _ V 2 Ay + 2 _y (3.39)

5z = aw2_ z+(1-a) [w+2 IwIVzAz +w-2 Iwl_Az ] (3.40)

where 5x, 5y and 5z are hybrid finite-difference expressions for the convective terms

o, <

c_ < 1. A central difference formula is obtained for a = 1 and an upwind formula

for _ = 0. A different weighting parameter may be used for each direction. Note

that this formula is written for the physical grid with uniform spacing. If there is

stretching, the metrics are required.

The notation in Equations 3.38-3.40 could be confusing. It is important to note

that the V and A in the numerator are operators and the A in the denominator

represents increments.

Let Q denote the finite-difference approximation to S and define:

Qn = Qn(i,j,k ) = Q(zi, Yj,Zk,t n) (3.41)
77

6xQn = (sRe( 62 -_x) Qn (3.42)

]

6YOn = Re (/Xy)2

6zQn = ( l_es ( _-_z)2 62 _z) O n (3.44)

solution of S may now be developed.

estimate of the solution for time level

the differences with respect to x at a

An alternating direction method for the

From a known solution at time level n, a first

n + 1 may be obtained by evaluating half of

provisional time level denoted by n + 1" (Aziz and Hellums, 1967). The source term

cS will also be evaluated at the implicit time level, so that

16 x (Qn+l* +Qn)+ cQn+l* + 6yQ n +6zQ n

= [(on+l*-on)/At]+O (3.45)

can be evaluated at n,n + 1⁄2 or n + 1. In the algorithm described here, 0

where

was evaluated at the n time level only. In the following, one asterisk denotes the first

approximation with more asterisks for successive estimates. Each successive estimate

(3.46)

16 (on+l* +0 + 1 (on+l** 16 (On+ 1 cOn+ 1

= [(O n+l -O n )/At] +0 (3.47)

The equations 3.46 and 3.47 may be simplified by subtracting Equation 3.45

from 3.46 and Equation 3.46 from 3.47, respectively. The resulting algorithm, which

is made by evaluating half of the differences in one direction implicitly.

16 (Qn+l** _ cQn+l** +

y + Qn/+ +

6zQ n

= [(on+i**-On)/At]+O

-_16 x (Qn+l* + Qn)
78

is second order accurate in time and up to second order accurate in space, is:

(_a: + 2c - 21At) Qn+l*

(_y + 2c- 2/At) Qn+l**

(6z + 2c- 2/At)Q n+l

= - (_ +26y+ 26z+ 2/At)Qn+ 2e (3.48)

= _yQn +(2c_2/At) Qn+l* (3.49)

= 8zQn+(2c-2/At)Q n+l** (3.50)

This is the algorithm used for the three-dimensional solution. A tridiagonal

system of linear algebraic equations for Q is solved three times to obtain Qn+l. The

first pass treats the z direction implicitly, the second treats y implicitly and the final

pass treats z implicitly. The alternating direction implicit method used here reduces

to the one proposed by Douglas if 6x, 6y and 6z are replaced by (52x,(Sy2 and _fz2. Also,

the non-linearities in the method of Douglas are included as part of 0 only.

In two dimensions the provisional time level is n + 1⁄2. The time increment is Tzxt

and the entire differences in one direction are taken at the implicit time level. In the

(x, y) plane, the algorithm to solve the transport equations can be written

Qn+1⁄2_ Qn

At2

Qn+l_ Qn+1⁄2

At/2

: _zQn+l + _yQn + cQn+1⁄2 _ 0

= _fzQn+1⁄2 + 6yQ n+l + cQ n+l - 8

In a form ready to code, the two-dimensional algorithm becomes:

(3.51)

(3.52)

: (6y + 2//_t) Qn _ 0 (3.53)

= (_x+2/At) Q n+1⁄2 -0 (3.54)
79

4. NUMERICAL RESULTS

4.1 Introduction

In this chapter the test cases will be discussed. It will be seen that the dual

potential method can be used to solve the full potential equation, the Euler equations

and the Navier-Stokes equations. The results generated from the dual potential code

will be designated as DP in the comparisons.

4.2 Solution Strategy

The solution strategy has been broken down into three parts: an incompressible

segment, a compressible segment and an iterative update. Certain steps are common

to both the incompressible and compressible procedure. The incompressible segment

and iterative update are always needed. The transport equations are solved uncoupled

and iterated only once at each time level. The ADI scheme for the transport variables

(-_, B,T,p) provides up to second order accuracy in time and space. The Poisson

equations may be iterated to convergence at each time level for a time accurate result.

It is possible to limit the iterations on the Poisson equations to speed up the solution

for a steady state case. The iterations may also be limited for a time accurate solution,

but numerical experimentation is required.
8O

The step by step calculation procedure is:

"Incompressible" segment:

___} ---}

1. Using known values (or assumed initial values at startup) for p, a_, B, T, V at all

points at a given time level, n, solve the vorticity transport equation by an ADI

--}n+l

method to give o_

2. Solve the vector potential equation using the n+ 1 level vorticity. If the calculation

is for a constant property, incompressible flow without heat transfer, go to step

7. Otherwise, continue with step number 3.

Compressible segment:

3. Solve the dilatation transport equation for B n+l using n level source terms.

4. Solve the energy equation for T n+l using n level source terms.

5. Solve the continuity equation for pn+l using n level source terms.

6. Update the properties #, k and p using n + 1 level quantities.

Iterative update:

7. Compute the scalar potential, _.

8. Update the velocity to time level n + 1.

9. Solve for the vorticity boundary conditions at time level n + 1.

10. Check for steady state convergence and stop if the solution has converged. Otherwise, transfer the just computed n + 1 level results to the n level and repeat
81

the sequence from step 1 until steady state convergence or to some point of

interest.

Though the segments above may be somewhat misnamed (since the equations

listed under the "incompressible" segment must be solved for a compressible solution),

the idea is to suggest possible solution strategies. One strategy which may be used

to obtain a compressible flow solution is to start with a converged incompressible

solution. The governing equations of the dual potential method are neatly segregated

by incompressible/compressible or constant property/variable property criteria.

The convergence requirement is (usually) applied to all dependent variables. For

example, a compressible case requires that -_, B, T, p, A and ¢ all satisfy the convergence requirement. The requirement is

I+n+l - +hi < e (4.1)

¢n+1

max

applied at each point, where (I>representsany dependent variable and n isthe iterationlevelor time levelaccording to the nature of the variable.The toleranceisusually

e = 10-5. This type of convergence requirement ismore stringentthan an L 2 norm

and more helpfulfor a self-interrogatingscheme to determine where to concentrate

computational eiTortand where to avoid computation. In particular,forthe variables

used in this method, the vorticityand dilatationcan be negligiblysmall. For example,in external boundary-layer flows the vorticityand dilatationcan be expected to

approach zero at some distance away from a disturbance. This point by point error

checking, normalized on a representativemaxhnum fieldvalue, can indicate that the

largesterrors are near the disturbance and initiatea check to see iffunction values,

--->

such as B or w, are small at some distance away. If this is the case, the code can
82

automatically reduce the computational field for that variable. The more common L 2

measure of error that other investigators use has been found to be somewhat deceiving. Convergence in the L 2 sense may not satisfy convergence in the absolute sense

above. Since the L 2 measurement averages out the error contribution of each point,

it doesn't identify particular locations of large error. That is, the L 2 measurement

may satisfy its convergence criterion and yet there may be isolated areas of large and

unacceptable error. Hence L 2 convergence may have little to do with the accuracy

attained by the solution.

Regarding the actual tolerance used in calculations, it has been observed that

the vector potential error tolerance may be loosened to speed the solution with no

negative effects. The scalar potential, on the other hand, requires a tight tolerance

in all cases. The potentials, given by solutions of the Poisson equations, can be very

time consuming to calculate.

4.3 Two-dimensional Cases

The dual potential formulation is rather thoroughly tested in two dimensions.

Example problems will be computed for incompressible and compressible flows. A

potential flow solution will be computed for the flow about a thin biconvex airfoil (or

bump on an inviscid wall). An unsteady calculation for the compressible flow about a

thickening airfoil will also be discussed. Steady viscous flows will be computed for the

channel inlet case and flat plate boundary layer. The results with discussion follow.
83

y, j index

inlet:

¢=0

u=l

freestream CY = _z

and bump: v =

!

bump

freestream: CY = 0

u=l

exit:

Cx = 1

u=l

Figure 4.1: Boundary conditions for laminar incompressible irrotational flow over a

bump

4.3.1 Incompressible flow

4.3.1.1 Steady irrotational inviscid flow

Bump cases This test case and the corresponding compressible test

cases for the flow over thin biconvex airfoils were included to demonstrate the calculation of simple irrotational inviscid flows. By considering only very thin airfoils, the

small disturbance theory boundary conditions may be used. The geometry is therefore

simple so that the dual potential formulation can be tested with fewer complications.

Also, velocity potential and stream function solutions are readily obtainable for com-
84

parison. The biconvex airfoil may also be considered to be a thin bump on an inviscid

wall. Airfoil and bump will be used interchangeably in referring to the bump test

cases.

Problems of this type were solved extensively in the 1970s as models for the

flow over helicopter blades (Ballhaus and Steger 1975; Beam and Ballhaus 1975).

The bump is generated by blowing according to the small disturbance theory (Ashley

and Landahl 1965). The Mrfoil in this case is very thin so that the freestream is only

slightly perturbed by the presence of the body. Tangency is then imposed by requiring

the resultant flow velocity at the body to be tangent to the thin body. Since the body

is assumed to be very thin, the tangency condition can be applied at the airfoil chord

line. This permits solutions to this problem to be made on a simple Cartesian grid.

The assumptions for this case are steady, irrotational, inviscid, isentropic flow.

A sketch of the problem is shown in Figure 4.1. This problem may be solved in

many ways using potential methods. Only one of the potentials in the dual potential

method is required to solve this flow. Either one of the potentials can be used with

little modification to the code. Small perturbation approximations will be used in

this problem. The velocity components will be represented by

u = Uoo+fi

v = _

(4.2)

(4.3)

where fi and _ are perturbation quantities. These are all non-dimensional quantities normalized on the freestream velocity so that Uoo = 1.0 and the perturbation

quantities are assumed to be much less than 1.

The scalar potential, _b, may be used without modification for a potential flow
85

solution to this problem. The definition of the scalar potential is designed to satisfy

the irrotationality condition automatically:

¢ = ux + vy (4.4)

Conservation of mass is then satisfied by the solution of

v2¢=0 (4.5)

The boundary conditions for a solution of this case using the scalar potential are

shown in Figure 4.1.

If the vector potential, A, is used it is actually treated as the stream function for

irrotational incompressible flow. This is a redefinition of the vector potential that is

used elsewhere in this development. The similarity of the stream function with the

vector potential is evident in the Poisson type of equation and Dirichlet boundary

conditions. This permits the single two-dimensional dual potential code to solve

this problem using either potential. Let A be the stream function in the following

development. The stream function is designed to satisfy continuity automatically:

Ax = u (4.6)

Ay = -v (4.7)

The irrotationality condition is then written

V2A = O (4.8)

The boundary conditions on the stream function are all Dirichlet with A = y. For a

compressible flow solution the vector potential requires a boundary condition modification to emulate the stream function.
86

The scalar potential is perhaps the most natural potential function to use for

this case since there is no need to change its definition. It simply assumes the role of

the traditional velocity potential.

Two different bump geometries were computed by the dual potential code. Solutions have been computed for the parabolic arc and the sine wave arc bumps. The

parabolic arc has stagnation points at the leading and trailing edges. The sine wave

bump was included to avoid stagnation points in the early development of the compressible code. The bumps are described using the small disturbance theory boundary

conditions. That is, the bumps are generated by blowing through the boundary. The

= (4.9)

PP

2 [1 - cos(27rz)] (4.10)

parabolic and sine wave arcs are described by;

parabolic arc:

Y

sine wave arc:

y -

for 0 < z < 1. r/2 is the height of the parabolic arc and PP (for peak-to-peak) is the

height of the sine wave arc. The sine wave amplitude is -_, but the sine wave arc

is used so that the bump meets the boundary with zero slope. The bump geometries

are indicated in Figure 4.2. The computational domain was chosen to be three chord

lengths in the streamwise direction and two chord lengths away from the airfoil to the

far field. The airfoil is centered in the streamwise direction. A uniform Cartesian grid

was used for the results presented here with 61 streamwise and 41 transverse points.

The airfoil surface will be defined by the flow tangency condition imposed along

the mean surface of the body. Flow tangency at the airfoil surface can be represented
8T

W

W

0_

J

_J

Z

C

01,0"0

c_

z _I_,

r...Oa- __I i

T,T "_'10,

_I m_

6_II_ i

Zln_,

rJ')IO_ _

I

I

I I I

S£0"0 0£0"0 5_0"0

I I I I

0_0"0 5tO'O 0[0"0 500"0000"0

I

tf_

-_X

c_

c_

i •

c_

0

o

,.o

o_,,_

P_

ImO
88

CO

D

C)

II

>--

-p0

-pC

(D

..D

L)

..)

(+_

(+_@0

C_)

(I) L

CO

O_

(I)

[_

0_

bO

CD6I

_E

"@

.COCOO

.JOL0

cl

£'0

Q_

O_

O_

O _LO

,_O

,_O

O

.qO

':3OOO

o

,_o

cco

0r z_ 'g uD L U O

._I Q)C

.-$ J

C

-J "O O O

Q_O

Om

O ,q

o_

o_

c_

Q_

I

_'0

I

I'0

I

0"0 d3

-

I!

I"0- 8"0-

C3

O0

-6

to

c)660

-6

-6IdI

£'0-

e_

°,.._ 0

c@

e_

e_

;.4 00

II8

_J

QJ O

Q;

d4

°_,,i
89

r,..D

a

C3

I!

0

C

©

._)

0

._)

qqO_

0

C_

L

CO

[I_

0

0

I

_C

T

.C

0

L

0

0

C

.J

_o

S_'0

O

o

o

O

O

I I I I I

0_'0 SI'O 01"0 SO'O 00"0

d 3 -

x

I-4

t_

e_

0

0

II

8

_e

0

4_

°,.m
9O

by

( fi U-oo U-oo (4.11) surface Uoo + fi Uc¢ 1 +

fi is small compared to 1 so it may be dropped leaving the In the above equation,

following approximation for the tangency boundary condition:

= ¢y for Uc_ = 1. (4.12)

surface U-oo surface surface

Finally, the tangency condition is applied at the mean surface of the airfoil so that

surface Uc¢ surface grid boundary

This small perturbation type of boundary condition does not require the use of a

body-oriented grid and is therefore very easy to implement. The above development

has been for the scalar potential. As used in the test problem, the boundary condition

on the scalar potential is

O<_z<_l

otherwise

(4.14)

where the solution domain is -1 < z < 2.

A similar development for the stream function gives the following small perturbation boundary condition

A = _ y + local airfoil height 0 _< z _< 1 (4.15)

[ y otherwise

The airfoil height referred to here is the half thickness of the airfoil at a particular z

position.
91

The pressure coefficient is defined as

(4.16)

where V 2 = u 2 + v 2. The pressure coefficient for the bump cases is computed by

an approximation to the definition

-- po¢ -2fi

-- (4.17)

Uc¢

This approximation is obtained by dropping squares of the perturbation velocity components. In terms of the scalar potential or stream function, the pressure coefficient

is approximated by

Cp _ 2(1 - @x) = 2 (1 - Ay) (4.18)

where U¢¢ = 1.

The pressure coefficient for an inviscid and incompressible calculation of M = 0.6

flow over a 6% thick (v = 0.06) parabolic arc airfoil is shown in Figure 4.3. This

result was computed solely to check the DP code in an incompressible calculation.

The linear theory result plotted in Figure 4.3 is obtained for the parabolic arc airfoil

by a method given in Ashley and Landahl (1965). That method gives the pressure

coefficient for thin airfoils in incompressible flow. It can be shown that the pressure

coefficient for a subsonic compressible flow is related to the pressure coefficient for

incompressible flow by the factor f_ = <1 - M2:

1

Cp = _ (CP)incompressible (4.19)

Hence, the DP incompressible flow result in Figure 4.3 is scaled to a compressible

result using Equation 4.19 just as the linear theory result is scaled from Ashley and
92

Landahl's incompressible calculation. A compressible stream function solution on the

same 61 x 41 grid lies right on the scaled incompressible solution for the parabolic

arc airfoil case. A similar result is shown for a 4% thick (PP = 0.02) sine wave arc

airfoil in Figure 4.4. The incompressible solution using the dual potential code scales

to the compressible result of a stream function solution. The computed incompressible solutions here have been scaled to the compressible result for comparison with

calculations obtained from available compressible stream function codes.

4.3.1.2 Steady viscous flow

Channel inlet The two-dimensional channel inlet flow has been computed by many researchers (Wang and Longwell 1964; McDonald et al. 1972; Tenpas

and Pletcher 1987). This case provides the first test of the vorticity transport equation

and the vector potential solver. Also, stretched grids must be used in this problem

for the first time to adequately resolve the inlet features. This problem will be a

stepping stone to the heat transfer cases presented subsequently. It will be seen that

the dual potential formulation can solve this flow problem for either a rotational or

irrotational inlet condition.

The developing flow in a two-dimensional inlet has been computed for Re = 10

to 7500. The Reynolds number is based on the hydraulic diameter,

4 x cross-sectional area

Dhy d = wetted perimeter

For a two-dimensional passage the hydraulic diameter is twice the wall spacing, h.

The physical distances, z and y, for the two-dimensional channel cases are nondimensionalized by the hydraulic diameter. The non-dimensional wall separation
93

distanceis then ymax = 1⁄2. Note that y for a two-dimensional channel problem will

then be in the range 0 < y < 1⁄2 for the results plotted here from wall to wall. The

geometry and boundary conditions are given in Figure 3.3. A uniform inlet velocity

profile was prescribed. The inlet conditions can be made to correspond to a rotational

or irrotational inlet flow. The irrotational inlet may be more realistic for experimental

models with a rounded entrance (Van Dyke 1970).

The centerline velocity development is shown in Figure 4.5 for Re = 300. Good

agreement is obtained for both the irrotational and rotational inlet conditions. Note

that the z used in the axis label, z/(Re* Dhyd)* 16, in Figure 4.5 is dimensional. This

was done to form the dimensionless group used by other researchers. The multiplier

of 16 is used to give the same z axis range as others who have used the channel half

height instead of the hydraulic diameter in the Reynolds number and dimensionless

axial length.

For the results shown in Figure 4.5, all convective terms have been central differenced. Refining the grid for the irrotational inlet case shows that the numerical

algorithm is clearly better than first order accurate but is not second order accurate

on a stretched grid. The truncation error is estimated to be approximately O (/kz) 1"6

where /_x represents the z and y grid spacing here. (The solution for the 41 × 41

grid is approximately one-third of the way vertically between the 21 × 21 grid solution

and the 81 × 81 grid solution). The need to upwind difference the convective terms

becomes necessary as the Reynolds number increases in order to maintain stability.

The skin-friction development is shown in Figure 4.6 for the rotational and irrotational inlet conditions. Both computational results asymptote to the expected
94

fully-developed GfRe = 24.0 for the two-dimensional channel.

The complete Navier-Stokes equations for incompressible flow were used for these

calculations. It is permissible to drop the streamwise second derivatives for high

Reynolds number flows. Dropping the _-'zz terms causes the flow to develop slightly

faster (i.e., in a shorter z distance).

Constant property heat transfer Heat transfer cases have been

computed for constant properties. The results are in good agreement with published

computational results (Schade and McEligot 1971; Hwang and Fan 1964). Both

references above obtained solutions for this problem utilizing boundary-layer approximations.

The heat transfer cases can be classified as in Figure 4.7. As shown in the figure,

there are three basic types of problems for specified wall temperature and/or wall heat

flux for the parallel plates geometry. Test cases will be computed for wall conditions

that are constant with x. The most detailed presentation will be for the constant

wall temperature case.

Any combination of the wall boundary conditions in Figure 4.7 can be solved by

the dual potential code. The boundary conditions may be functions of x.

The notation used in this heat transfer section is the same as used by Shah

and London (1978) and Kays and Crawford (1980). The x distance in the following

heat transfer results is referred to the Peclet number as in Shah and London (1978).

Therefore,

, x

z - Pe = RePr

Pe'

A quadratic fit to the computed temperature profile was used to compute the local
95

0

0

.._

C

Ill

E

0._

0

>

[b

CD

...)

0

0

Ill

[b

C

._9

J

C

O_

C

C.3

r_

...._

cc

0

..C

0

O

I

_g

._ _._

%.

z

W

U

0

°

I:D 010 O:

n," n_ : ._._

r_ r,,._lr,_q C.9!

f i _ --J --_

'_ _: C

• _ .' :0 0 0

_i_I "< C C

0_01 .: _ 0 0

I_ I---.: C O r'_

I_1_ o:o o o

x,",x,

\ "'b,^

w

X \

\

S'I t,'I £'I _'I I'l O'I

uoqn/n

!

0

CO

"O

.32

C3

)1[

0

n_

X

0

,,?

0

t-I

_o

E_

o

o

I::I

r..)

I-,I

o,-¢
gfi

___) ©CC0

dED

CDI

C_0

O] C0

_b0C0

._)0

.J

L

q_C

.J

U]

@

DO!0

n_

__D

q-

.JL

Q_E0U

.SC

.J

E0

_J

00!

._o. o_

.,o .___4•

o cc_

-J, .J c

n

0•

9_ _i 0 0

_

_C .J (.,_i 0 -o fJ'

_1 "J0

00U

.o t.0

0LX

[]C .J0

•i

IO

[]

[]

[]0

[]

[]O

[]OO00

[]0

[]

[]OOO0OOt

o noo o °°oo I I

ff

O'OO o'og 0"0_ O'OE

e_l_JO

000

-o_0

O0000

,.4000o06

0"0

X

oo

.o

o_,-_

I:::1

°_,_

I...i

0e.._
97

' o_.T_"" , ' qwl ' ' qwl !

I I I ! I I

I I I I I I

'

, ', ', qw2 ', ', ,,

Figure 4.7: Basic wall heat transfer boundary conditions

Nusselt number, where

Nu- hDhyd

k (4.20)

The Nusselt number was computed from the temperature solution. The notation of

Kays and Crawford (1980) will be used in the following derivation. Specifically, an

overdot indicates a time derivative, the subscript o will mean at the wall and a double

prime (tt) will indicate per unit area. The wall heat flux can then be written as

qtor = h(Twall- Tmean) (4.21)

_tot = k OT

- -_Y wall (4.22)

The wall heat flux is given by either the convection heat transfer equation or the wall

conduction equation as written above. These two must be equal, so

h(Twall - Tmean) = -k OT

Oy wall

(4.23)

The definition of Nusselt number from Equation 4.20 can then easily be written in

terms of the computed temperatures by rearranging Equation 4.23. The result is

OT

Nu - hDhyd --_ wallDhyd

k - (Twall- Tmean) (4.24)
98

whereTmean is the local mixed mean fluid temperature. It is also sometimes referred

to as the bulk fluid temperature or the mixing cup temperature. Tmean is computed

as the single (mean) temperature for the axial convection that equals the integrated

axial convected thermal energy rate. That is,

fitcpTmean = (p Ac V ) cpTmean (4.25)

= fAc pucpT dAc (4.26)

where fit is the mass flow rate, Ac is the cross-sectional area at the z position for

which Tmean is being computed and V is the mean velocity. Therefore,

Tmean - f Ac uT dAc

fAc U dA c (4.27)

The constant property cases to be presented have fully-developed temperature

solutions which are self-similar. The self-similar temperature function is

O = Twall- T(i,j) (4.28)

Twall - Tmean

The variable, O, will be referred to as the temperature parameter in the following

results. It is defined and used in the dual potential code to indicate when the temperature solution has become fully developed. This is in addition to an error tolerance

on the temperature at each point.

It is important to pack many points close to the channel inlet, so the grid was

stretched until grid independent solutions were obtained for a 41 x 41 grid. Those

stretching parameters were then used for further grid refinement by adding points.

The length of channel was also chosen so as not to interfere with the natural flow

development. The fully-developed skin friction for these two-dimensional cases is
99

CfRe = 24.0 as computed analytically and by the dual potential code. It is interesting to note that all of the results for the heat transfer cases were computed for

the rotational inlet condition. This is the inlet condition used by others reporting

computational results. Van Dyke (1970) reports that the irrotational inlet condition

may be more realistic for comparison with experimental models having a rounded

entrance. The dual potential code can easily handle either a rotational or irrotational

inlet condition.

Constant wall temperature Some comparisons of constant property heat transfer results for constant and equal wall temperatures are

shown in Figures 4.8 and 4.9. The published results of Schade and McEligot (1971)

are represented in each figure as a solid line. The present results using the dual

potential code are indicated by symbols. The grid size for the dual potential results is

indicated in the figure legend. The Reynolds number based on the hydraulic diameter

is 150 for the present results in Figure 4.8 and 300 for the present results in Figure 4.9.

Schade and McEligot use the boundary-layer assumptions and drop their numerical

results for approximately the first 20 points in z. The dual potential results for the

161 × 81 grid have only the first 4 points dropped. There is nothing wrong with those

4 points, they are just well ahead of the data presented in the literature. The dual

potential incompressible code makes no assumptions other than constant properties.

This accounts for the Nusselt number discrepancy between the two solutions near the

inlet. The dual potential code predicts a higher heat transfer rate near the inlet than

the results of Schade and McEligot. Figures 4.8 and 4.9 show that there is a Peclet

number dependence for the Nusselt number. For the range of results presented, a
100

(])

C

.-4.3

0

C

(])

Q_

E

(])

[.._._,

E

0

O3

E

0

(_3

CO

o)

4-)

0

J

(3_

O)

_3

_._3

0

f._

©

Q__

0

Lf3

m,

C

..3

__2

0

..2

0

C

0"05

U') LO I.f')

000

! ! I

O_ O. O_

O" • • •

°z,.d, , ,

W U

j "13 0 0 0

E

0 _ [-_ [-_

-0 _4 03 O0

0 X X X

0 _l_ O0 OO

CO ,--4

0 'q +

+

O'Q

+

+

Jr

,q

I

0"0_,

I

0"0£

I

O'OE

l'XnN

I

0"01

_0

_0

9

_0

T

i 0

0"0

X

X

¢@

¢@

0

@@

c@

o_:

'= II

_z

II

Q _
101

(1)

(

0

(_

(])

O_

E

(])

C

0

C

0

CD

_J

0

__)

O_

__)

_D

__)

0

C

0

O_

0

0

rY

C

.J

._)

0

c

o

0

f_

O'OS

CD 0

! !

• •

13_ O-

.j _,-i v-4

_,,

!.1 u

_ .J _J

_.i "13 0 0

0_

"O_CO

0 % %

O,_OO

O3

0

0

I I I I

0"01, 0"0£ 0"0_ 0"0!

L'xnN

0"0

"7

_0

,I.--,4

_CD

,,?

_CD

%

X

I-I

I...4

o

I=1

,.o

II

_._

_ .

z_

_ H
102

Peeler number of about 210 is required to give good agreement with the boundarylayer assumptions. Lower Peclet numbers give a higher heat transfer rate near the

inlet. Profiles of the temperature parameter,/9, are plotted in Figure 4.11 for several

streamwise positions. The fully-developed profile is parabolic.

To match the boundary-layer assumptions used by Schade and McEligot, the

streamwise second derivatives must be dropped. The affect is similar to assuming a

high Reynolds number flow. Schade and McEligot report that Pe > 100 is sufficient

for the boundary-layer assumptions to be valid. Also, Schade and McEligot neglect

viscous dissipation. For incompressible flow, the only streamwise second derivatives

in the dual potential method are wmx, A_:m and Txm. When these derivatives and the

viscous dissipation are dropped, the results approach those of Schade and McEligot,

as can be seen in Figure 4.10. It is clear that the boundary-layer assumptions cannot

be used to get accurate results near the inlet at lower Reynolds numbers or Peclet

numbers. The boundary-layer assumptions neglect the important elliptic effects at

the inlet and underpredict the heat transfer.

The computation rate for the dual potential code is 80-100 MFLOPS overall

using a single processor on a Cray X-MP. The computation rate can be stated in a

5.5/us

more useful manner as _,- global iteration × grid point on a single processor of the

Cray X-MP. The number of global iterations is on the order of the number of grid

points. A conservative estimate of Cray X-MP cpu time is:

cpu time (s) = 5.5 × 10 -6 × (number of grid points) 2

The two-dimensional dual potential code is a compressible code. The compressible

terms are computed for this calculation even though they are zero. If the code was
103

designed purely for incompressible solutions the code would be much fasternat least

twice as fast based on a rough operation count. The error tolerance used for these

calculations was ¢ < 10 -6.

The solution to the constant property case is independent of wall temperature.

This was checked using Twall = 2.0, 1.1, 1.01, 1.001, 1.0001 with the same constant

property results obtained for each case. Shah and London (1978) also present data

for the constant wall temperature case. For the present results, the grid was stretched

using cr = 1.05, c_ = 0.5 and fl = 1.2.

Constant heat flux The constant temperature case has

revealed that Pe = 210 gives good agreement with the boundary-layer approximations

used by Schade and McEligot. Just one case will be presented here in Figure 4.12. The

results compare well to solutions obtained using the boundary-layer approximations,

but are more accurate near the inlet. A second order polynomial curve fit to the

temperature distribution was used to obtain the wall temperature. The temperature

parameter profiles for this case are presented in Figure 4.13.

Mixed wall boundary conditions Solutions were computed for two cases with one wall insulated and the other wall at either constant

heat flux or constant temperature. Shah and London (1978) refer to these boundary

conditions as the fundamental boundary conditions of the second kind and third kind

respectively. The asymptotic Nusselt number for the case of one wall insulated and

the other at constant heat flux is computed analytically to be Nu = 5.385. This case

will be designated by H2 in the results. The asymptotic Nusselt number for the case
104

G) LD

.+_) 0

£_

[b

O_E

-pE0

O] E00

CO

(b

-4-> 0

__)

0_.

__)

[D

__.)

_.)0L0

O_

0"05 O'OE

I- m

0"0[

0

0

%

0"0

o

---7,o

×

t-I

2:; _"

.°

._.._
i

!

'/!I

!

!//

105

\

i

£'0

o

_o

_0

°_
106

X

__D

h_0

G)

n-C0

£0C0

C_3

CO

(D0

__)

[3_

_.i)

(])

___)

_i)0C0

EL

o

c37

(D

fY0JC

i_D0

Co

.-)0

O'OS

oi

h.11--*

o I z 70 w c ID

_._o Qbaj ID -_ i0o O0

.EX

0 ,-'

03 CO0 0

0/

IIII

0"0_ 0"0£ 0"0_ O'Ol

Xl7N

,l"I4 ?

_C3

,,2

_0

"o

0"0

XX

4,a, o

,.c:l

_._

°_,,_
107

o

,..)

,,c:l

H

0

_>

I.=I

I=4

I-I

¢)

.=-i

q,) I-4

o,.=_
O_

C

0

..D

.+J

._D

-0

C

0

CD

(._

0

-0

C

0

O3

_J

0

Z

-0

fl)

X

._b

7

O)

fl)

J

0

__)

fl_

__)

fl)

__)

_D

0

L

0

(]_

0

0

V

@

(]d

2

C

.J

_J

0

.J

0

0

L

O'OS

108

I I e_l Drl i

-vii

@ • .

__o_-_l-_i

0 I0:

_ -r 1- .j .j:

,., .o g4 J _J:X X 0 O:

,.it ,.._ r:

[]

O

[]

0

o

0

0

D

O

[]

0

o

0 0

O 0

0 0

O 0

O 0

[] 0

0

0

O

0

I

O'OY

I

0"0£

XflN

I

O'OE

I

0"0!

0

0

_CD

0"0

x

x

s--4

II

P_

0_

_; o O

_-_

0 0

°_,,i °_,,i

_ e

_ II II

.,., Z Z

0 0

4_

°_,,_
109

O)

co g

(]) .j

._.)

_ c

o o

L_

o

(1)

O_ a::

L.

C m

-_ ..J

0 j

L o

g

f") .j

E .o

C

........................... •.

"....,o.

! • :

o o,_'£i_

Z_gw 0,01 £1919i 9

wj 00lOI IIIII o!o ?

X

×

S'O

I I I I

I_'0 £'0 _'0 ["0

0

o

0

0"0

_J

0

0

_o

_8

_O
110

0

0

0

_8

°_,,q
111

of one wall insulated and the other at constant temperature is computed analytically

to be Nu = 4.861. This case will be designated by H3 in the results. The computed

Nusselt number development for these two cases and the analytical asymptotic Nusselt number are shown in Figure 4.14. Temperature parameter profiles for these two

cases are plotted in Figures 4.16 and 4.15.

Boundary layer This case is considered to provide yet another check

on the accuracy of the dual potential scheme since the laminar boundary-layer solution is well known. Despite the simplicity of this flow, it is not always easy to obtain

good comparisons with numerical solutions to the full Navier-Stokes equations.

The boundary conditions for the laminar incompressible flow over a flat plate are

shown in Figure 4.17. Since this is an incompressible constant property flow, B, T

and p need not be computed.

Concerning the vorticity boundary conditions, an irrotational or rotational inlet

flow may be specified. The best comparison with Blasius results was obtained for

a zero inlet plane vorticity (irrotational condition). This gives a non-zero v velocity

at the inlet. Both u and v are fixed to be zero at the stagnation point. The wall

vorticity is computed as in all viscous wall cases. The freestream is assumed to be

irrotational. The exit condition is based on the usual fully-developed assumption that

second order streamwise derivatives are zero.

The vector potential inlet and viscous wall conditions are the usual A = 0. In this

case the vector potential is used to account for the throughflow velocity (i.e., leakage

velocity) at the top boundary. It was stated in the boundary condition section that

only the scalar potential would be used to give throughflow velocities. This is the only
112

y,j index freestream:

w=0

Az = Cy - v

Cy = 0

u:l

vy = 0

p ...........

inlet:

03 : 0 or_

03 _ Vy -- UX

A=0

¢=0

u=l

v#O

_\\\\\\\\\\\\\\\\\\\\

03 = f(A,¢)

A=O

flat plate: Cy = 0

u=0

v-=0

exit:

03XX _- 0

Azx = 0

UXX _ 0

"Vxx _ 0

Figure 4.17: Boundary conditions for laminar incompressible flow over a flat plate
113

exception. The v velocity is given at the top boundary by assigning the same velocity

at the boundary as at the point just inside. This implies that 0v = 0. Hence, Az

along the top boundary is computed from the assumed v leakage velocity using the

velocity decomposition for v, v = _by - Az. More elaborate calculations can be used

to compute the v leakage velocity from the displacement thickness. For example, the

leakage can be accurately written as

dg* v

dz U_

Thus the leakage can be controlled based on a known or computed displacement thickness, 6*. These ideas worked no better than simply assigning v(i, nj)= v(i,nj- 1)

at the freestream boundary. At the exit ua: = 0 was acceptable, but uxz = 0 gave

a slightly better boundary-layer profile as compared to the Blasius (1908) solution.

Note that the exit condition on the vector potential, Azz = 0, is equivalent to an

extrapolation technique. The vector potential at a point just outside the exit plane

boundary is computed by Azx = 0 for use in the Poisson equation for A. The Poisson

equation is solved at the exit plane. This assumes that the flow is fully developed and

removes the need to specify the velocity at the exit. The scalar potential boundary

conditions are the same as for the incompressible channel case. Therefore, the solution of the scalar potential Poisson equation is known a priori and _b = x is simply

assigned and not re-computed. The same boundary conditions as the channel case

were used to avoid a solution for the scalar potential in this case. The only difficulty was a need to change the top boundary condition on the vector potential to

account for a leakage flow. This leakage flow is small and could be neglected with the

boundary far enough away. If it is neglected and the top boundary is too close, the
114

boundary-layerprofile tends to overshoot the freestream value by several percentage

points. Allowing leakage at the top wall permits the solution to be obtained on a

smaller grid.

The laminar boundary-layer computation agrees well with the Blasius solution.

A boundary-layer profile comparison is shown in Figure 4.18. An 81x81 grid was

used with stretching in the y direction. This puts about 22 points in the boundary

layer at Rex = 312,500. The solution is attained in 54.5 seconds at 85 MFLOPS

overall on a single processor of the Cray X-MP. The data marked with triangles

were computed by limiting the solution domain of the vorticity transport equation to

31 transverse points. This example was specifically used to demonstrate the ability

to reduce the computational effort by dropping the calculation of vorticity where

it is negligible and still obtaining a full field solution for the problem. Even faster

solutions are possible using this method. Normally the Poisson equations dominate

the solution time. That does occur in this problem also, but iterations on the vector

potential calculation were limited with no deleterious effect on the final solution. With

the vector potential limited, the vorticity transport equation becomes the dominant

time consuming computation for this problem. It is this condition that was used to

demonstrate the ability to reduce computation time for the vorticity. This speed-up

affect is shown in Table 4.1.

For problems that one knows have vorticity approaching zero at a boundary,

the point-wise convergence check used in the code can be helpful in deciding where

and how much to limit the solution domain for vorticity. For example, if vorticity

goes to zero at the freestream boundary as in this problem, one will note that for the
115

Table 4.1: Speed-up affect of limiting the vorticity solution domain

81× 81 grid. Re L = 312,500

Transverse pts for vorticity solution CPU time (s) MFLOP

81 54.5 85

41 40.2 77

31 36.9 73

convergence checking used here (which is point by point and normalized on the largest

value in a plane), the biggest error will be in the region where vorticity is important.

If the point of maximum error is far away from the known _,, = 0 boundary, then

limiting of the solution domain is indicated.

The MFLOP rate in Table 4.1 drops as the solution domain for vorticity is limited

because the vorticity transport subroutine had a high MFLOP rate. The solution was

speeded up by limiting the vorticity solution domain and thereby reducing the CPU

time spent in computing the vorticity in regions where it is nearly zero anyway. The

vorticity transport subroutine is a 107+ MFLOP code. Reducing the computation

time on vorticity reduced its high MFLOP rate in the weighted average to compute

the overall MFLOP rate for the code.

Solutions were obtained for Rez up to 2,000,000 using both rotational and irrotational inlet conditions. The rotational inlet condition gave a higher skin friction and

a slight velocity overshoot. The computed skin friction along the flat plate for the

two inlet conditions is shown in Figure 4.19. The rotational inlet condition produces

a higher Cf than the irrotational inlet by the conservation of vorticity. Since the

inlet vorticity is non-zero in the rotational case, there is an additional amount of vorticity in the rotational as opposed to the irrotational inlet case. The boundary-layer
116

calculation for the rotational inlet (w # 0) has vorticity in excess of that assumed in

the Blasius solution. The extra vorticity is the inlet amount. The Blasius solution

is obviously for an irrotational inlet. The results presented here were obtained using

first order upwinding of the streamwise convective term. Weighting toward central

differencing works but requires more CPU time to converge. The equations used here

do not use the boundary-layer assumptions and so can accurately model the near

leading edge flow. The skin friction is computed as:

Ou

C f _ = "walll _ wall2 t/pc_ Vc_ x (4.29 )

_pc_V_ V _c_

The velocity derivative is computed using a second order polynomial fit to the u velocity component. The computed asymptotic value for the skin friction is C f Rx/Rx/Rx/Rx/Rx/Rx/R_x =

0.664 for the irrotational inlet case and Cf Rv/-_x = 0.704 for the rotational inlet.

The $ and y distances in Figures 4.19-4.21 are non-dimensionalized by the length

of plate required to give Re L = 312,500. Velocity vectors for this boundary-layer flow

are plotted in Figures 4.20 and 4.21. The most easily observed difference in these two

figures is the flow at the inlet plane. The irrotational inlet flow condition produces a

v velocity component at the inlet plane.
117

0

e:l 0

t,=l

p,, 0

Q;

Q;0

.,_

.,_

e_

66

,@

Q;

_O

o_._
118

00

IF)

(-',3

["0

II

_....1 0

__)

0_00

03C00C0

._)0

._)

C_C

.D

(/3

C

._) _.J c J ._) o

c--_ o o

._C

•OO

NO .--)

•O .O L .O

.-)L -_

0_ 0_ --J

¢'_O 0CIi [] O

OD

Nn

OO°O IOO

°l OOO°O IoO

°; OOO°O I00

O0

OO

OO

OO

O0

OO

OO

OO 0

O

0 []

0 []

I 'IIII

0"[ 0"0 9"0 _'0 [ xe_ ] l_lOSx3O_'0

0

C)

O0 0S

tO

C)

C)

c_

('g0

c_00

0"0

X

cm0

"-a

L_

a_

I-i

°,,._
119

1

I

0

l l lllll 

N.

P

0

0

11

l

11

II

1

llllll 

llll 

11111 

1,1,11111I 

qlD In ,q' _ N

0 0 0 0 0 O.

0 0 (_ 0 0 0

0

g

(;x

0

0

0

0

O

0

0

>..

0

ov,q

0

"4.-)

0w,q

0

,g.q

0

0

0

f_

0

"V

0

oo

¢,q
120

I

q@

0

0

I

!

0

0

I

IIII 

IIIII 

IIIII

IIllII 

II,rIT T,TTTT

0 0 0

0 0 0

w

Clt

0

m

P.

c;

g

0

0

0

N.

0

"Z

0

0

0

O

O

>,

¢:

0

0

0

0

4,)

i:=

f_

O

0

0

¢.)

0

_>
121

4.3.2 Compressible flow

4.3.2.1 Steady irrotational inviscid flow

Bump cases The flow assumptions for this problem are the same as

in the previous bump test case:

1. inviscid

2. irrotational

3. isentropic

4. no body forces

The flow can be incompressible or compressible and steady or unsteady. The boundary conditions for this case are shown in Figure 4.22. The same boundary conditions

apply for the unsteady problem which is discussed in the next section. The assumptions for this test case are the same as in the bump cases previously discussed on

pages 83-92 except that the compressible solution is computed directly rather than

scaled from an incompressible result.

There are many possible well-tested solution methods for this irrotational problem. The dual potential method is easily adapted to solve the flow using any of the

following combinations of governing equations. Either potential may be used. The

scalar potential may be used as the traditional velocity potential for a full potential

solution. The vector potential may assume the role of the stream function for a compressible stream function solution. In addition, an Euler mode of solution is available
122

y,j index

inlet:

B=0

¢=0

p=l

T=I

u=l

v=0

freestream:

B = Vy initially, then by continuity

p from Bernoulli

T from constant total enthalpy

bump

B from continuity

Cy = 0

p from Bernoulli freestream:

exit:

B=O

Cx = 1

p-1

T=I

u=l

v = Cy

D

z, i index

T from constant totalenthalpy

U=¢X

v=0

Figure 4.22: Boundary conditions for laminar compressible irrotational flow over a

bump
123

using the dual potential method with viscous terms and thermal conductivity set to

zero.

The traditional potential methods simply require the velocity potential or stream

function with a version of the Bernoulli equation. These two methods will be presented in a form only applicable to steady flows. The full potential equation written

in conservation law form is

0 0

(pCy) = 0 (4.30)

This may be solved by the scalar potential solver in the dual potential code with a

redefinition of the scalar potential derivatives to be the parenthesized quantities in

the equation above. This requires a modification of the boundary conditions also.

Since the flows to be computed here are subsonic, it is permissible to write the full

potential equation in non-conservative form:

V2 ¢ = PxCx + PyCy (4.31)

P

This is easily solved by the scalar potential solver with no modification to the boundary conditions. The steady compressible Bernoulli equation is then solved for density,

1

Using the stream function, only the velocity decomposition and Bernoulli equation are required. For the stream function approach, solve the following:

1--+ V = V×¢

P

1
124

T=

u 2 + v 2

2Cp

The stream function, ¢, above is defined to satisfy continuity automatically by:

Cy

U --

P

V --

P

The above would require some code modification to define the vector potential,

A, as the stream function, ¢. An alternate approach is to define

CX

Az -

p

¢_EY

Ay =

P

then the vector potential solver remains unchanged for a stream function solution and

only the boundary conditions must be modified by p. This also permits easy use of

the non-linear form of Bernoulli's equation which becomes

1

The energy equation is not specifically needed if the Bernoulli equation is solved since

it is a statement of conservation of momentum and energy.

The dual potential formulation may be easily converted to an Euler solver by

dropping the viscous terms and thermal conductivity. All that is required is to set

# = k = 0. For the bump problems, the flow is irrotational so that the vorticity

transport equation is also dropped. The DP code in the Euler mode solves the

governing equations in the following order:

1. V2¢=B
125

2. v=v¢

3. the energy equation for T

4. an equation for B

5. an equation for p

6. the equation of state, p = pRT

The assumptions for this problem allow the energy equation to simplify to a

statement of constant total enthalpy. There are several equations to choose from to

compute density and dilatation. Using the DP code in Euler mode, three possible

equation sets are given below. These are listed in order of increasing complexity, or,

in other words, from fastest to slowest.

--.+

v=v¢

V2¢= B

u__v 2

T = Tt- ZCp

p by Bernoulli

B from continuity

p by Bernoulli

B by its transport equation

p from continuity

B by its transport equation

All these equation sets give the same results. Obtaining p from the Bernoulli

equation and B from continuity is as fast as the compressible stream function ap-
126

proach. Using the continuity equation to obtain p and the dilatation transport equation for B takes about 10 times more CPU time than the fastest methods. Only the

last equation set will be applicable to an unsteady flow, provided the energy equation

is also time accurate.

For the compressible flow over a bump, the dual potential code uses the following

initial conditions:

1. u = 1, v = 0 except at the bump where v = _ (the slope of the boundary)

2. B = 0 everywhere except at the bump, where B = Vy (Euler version only)

3. p = 1 everywhere

4. ¢ = 0 everywhere (Euler and full potential versions only)

5. A = 0 everywhere (Stream function version only)

6. T = 1 everywhere

7. # = k = 0 everywhere for all time.

There were initial difficulties solving the bump test cases due to abrupt function

changes at the boundaries. It was necessary to smooth the variables (T,p,B) into

the boundaries either by extrapolation or by solving for the variable at the boundary

using available information and a governing equation--usually a simplified equation

such as the Bernoulli equation. The far- and near-field boundaries then have values

that are near the expected steady condition but differ slightly to provide smooth

derivatives into the boundary. For example, the density is not exactly 1.0 at the

far-field boundary when it is computed by the Bernoulli equation at that boundary,
127

but the derivatives of density which occur in the dilatation transport equation and

elsewhere are much smoother using this approach. The y=constant boundaries were

the most critical.

The pressure coefficient can be computed from the definition

cp

or, if the scalar potential is used,

p - poo

- 1 2

2poo V_

% = -2¢ t-2¢x-¢2

whereV 2 =u 2+v 2.

The flow over a 4% thick sine wave arc airfoil and a 6% thick parabolic arc airfoil

was computed for a freestream Mach number of Moo = 0.6. The dual potential code

was used in a compressible stream function mode and in an Euler mode. Results

for the sine wave arc are shown in Figures 4.23-4.28 and results for the parabolic

arc in Figures 4.29-4.34. The DP code is solved on a uniform 61 × 41 grid. The

computational domain is 3 chords in the streamwise direction and extends 2 chords

away from the airfoil. The airfoil is centered in z. The sine wave arc was the first

successful compressible calculation. There are no stagnation points to contend with

on the sine wave arc. Note the fair agreement between the stream function and Euler

solution for this case shown in Figure 4.23. A grid refinement study shows that the

61 x 41 grid gives an Euler solution independent of further refinement. The solutions

have not been obtained on a larger domain to examine whether the solution is grid

independent in that respect. The discrepancy between the stream function result and

the Euler solution is caused in part by the inconsistency of the small perturbation
128

boundary condition with the stream function and Euler solvers. The intent here is to

show the capability of the solver, not ultimate solutions. The Euler solution requires

about 3-10 times more CPU time than the potential solution using the same DP

code--depending on which equations are used to solve for the density and dilatation.

In the Euler code the convective terms are weighted toward central differencing by

0.9. The Euler solution is not exactly symmetric as can be seen in the contour plots

of dilatation (Figures 4.24 and 4.30). This has been observed by others for inviscid

solutions.

4.3.2.2 Unsteady irrotational inviscid flow

Bump cases The flow assumptions for this problem are the same

as in the previous bump test case, except that the flow will be unsteady. The dual

potential code will be used in an Euler mode. Constant total enthalpy has been

assumed for the energy equation. Formally this requires a steady flow, but it will be

applied here at each time step for the unsteady problem. The particular unsteady

problems to be computed here have primarily low frequency disturbances (Ballhaus

and Steger 1975; Beam and Ballhaus 1975). The validity of using the energy equation

in a quasi-steady fashion was evaluated by computing the unsteady flow with the

complete energy equation. The pressure coefficient results were within 4%. This

problem was studied mainly to test the dilatation transport equations and the density

determination from the continuity equation. The constant total enthalpy (or constant

total temperature since cp is constant) is written

u 2 + v 2

T t = T + (4.34)

2cp
129

I

Ol"0

I

SO'O

dOI

00"0

X

° ,,,,,t

,,--i

0

o

0

II v

8 °

d_
130

¢,_

/_

I I I I

k") 0 I/I

r,. Ln N

0 0 0

>,,

0

m

Wl

t.,l

t_

P.

0

0

I0

6 x

0

0

0

!

0

0

I

0

I

7

0

f-i

.,,_

f-t

;>

p.

0

0

II

8

N

o

0

¢)

of,,_
131

cq

cD

• m

__)

._) -_

__ _co3

0 .-)

o

t_ _c

CL

__ o

._)

C o

--D @

c

._)

t_O'I

S .'

t' ' :T[

' I ,I 7["-, (:10,0"1 0 i

t_,t ¢) '_'i

I

I I I

£0"I _0"I I0"I

I

O0 [ 66" 0

g_Tsuep

I

86"0

0

o

N

o

i

0

o

U')

o

0

×

o

m

-c;

0

o

m

0

I

L_

0

!

-6

I

o

(El

, T'

Z6"O 96"0

I-I

¢I

0

0

_tJ

0

t-I

,.;5

1.4
132

J

J

0

:3 rJ

E_

EL o

0 i

f'- .j

0

CO

0 •

-o 5

_...)

.D

r-7

_..-.

I I I :l

I I I I I I I

5['0 OI'O £0"0 00'0 £0"0- Of'O- £1"O-O_'Ou07q.o3oq7p

0

0

Q

o

I,n

0

_. x

o

N

o

Q

o

0

to

I

0

o

I

-c;

I

0

C3

I...I

u

4=1

,.£o

0

'-a

0

0

I,-i

0

° ,.,.._

°,-i

I-i
133

cq

c3

CD

!

___) .-)

(+_

o G

C L

EL o

> >

o

c

...)

5LO'O

I

OSO"0

o

o

I

C)

I

X

t_

e_

I-4

0

,---4

0

Q)

0

0

_0
134

_.b

X

0

o

c_

(1) .o

_.--) ..C

C

_c

E .,

0 o

(_ c_

o

c+._

C >

o

0

_.D

.D

C7

Og'O

I

J

I

:f

,%

C

0

.J

0

tr_

Z _-

bJ

___ E

bJ 0

_J •

t_

_J

I

SI "0

! ! !

"-'_ "7 "7

¢_

@

:d

i£

Q

+ × 0

I

OI "0

I I

SO"0 00"0 SO'0-

UOTqOlo77P

I I

01"0- Sl "0-0_'0-

0

o

o

s_

Q

0

i

0

o

I.N

o

o

o

c;

Lfl

I

0

-c;

I

c_

I

0

o

"7

0

F_

c@

O

O

_._

_ .,-_

t_

o _3

d6

°,i_
135

__..] 0

n_

(Ii Z

TE-,

(Zl

r_

r,100JL

(_

D

0

['000

II

C-

.J

L00

._)O

£'0

O_

O_

O_

O_

0

0

°f

b _- <3 _o---_

6.1L

(.90

U (9

--]C

.J

_.J

(_ ._J

@W

O_ OL

0£3 0 <_

0

0 ,Q

0'_

0"_

O_ I

_'0 I

I'O I0"0 dO- II'0- I_'0-

(N] 0

000

tO6

tN]

-c;o

-c;

CN] !0!

£'0-

X

_P

° _,,,i 000•

0

8°

_D

8_

°,,_
136

_t

0.00

_ o _.

p_

O

O

mr1

O

3_

O

N

in

mr_

N

O

O

Q

N.

O

!

O

O

I

O

I

"7

O

O

>

O

O

II

8

_o

S-4

O

O

o

°°

c_

t_

b0

°,.q
137

0

o

I

Z6"O 96"0

X

,.--M

¢J

f-4

_J

0

0

o

0

G)

_0

°_._
138

_._)

_b

0

c:)

D__ o

0 i

C" .-_

0

C

...) _.9

o

0 o

_.)

.__)

t, 0

---_-_"__ _

• • ..... .

¢¢'° _

.. i I

%- - "..........:.2,'-

 .loiCrJ!I C:_ i ---_

 lT',7'T:T

I !

fi'O Z'O

! I

1"0 0"0 I "Ouo_1olol_p

I

C'0-

!

C'O- _'0-

C_

Q

N

Lt)

m •

o

u_

u_

12)

c_

CD

u_ X

CD

u_

-c_

o

CD

-c_

U3

_XJ

"c_

!

0

L_

I

L_

c_

I

0

0

T

°_

0

..o

0

o

P_

0

.2

°_=_

C_

_JD

°_
139

Dr)

(2)

C]

!

(D o_

___) .__

c+_.

o G

C_ c

O_ o

0

._)

o

_o

o

[_

0

O_

SI"0

I l

01"0 SO'O

I

00"0 SO"OR3?oO]OA A

0

o

c_

u';

o

0

D

r,.

o

0

Lt)

-- • X

C_

0

0

C_

C_

c_

-c;

!

0

-c;

I

-c;

!

0

0

T

SI "0-

0

c@

0

0

0

0

;>
140

.__)

X

(D

o

(D ._

....9 _.c

C ..Dm

---) •

..C

o

L

0

(3

._)

._)

o

(_

E

0

(_

q_

C

0

._b

0

+)

0

__9

__)

C_

9O

0

@

0.' :

,.Y

J

! I

.I' _

"_ '6----" '

C I

0

.J

0

C

r_

Z _-

ILl

C__ E

_,oII)

t.

I I I

I:--

n

I

!

I

I

t,'0

I

_'0

FX :.f

J

I

!./°

)

I

0"0 _'0-

uo7qoqo77P

I

_'0-

0

c)

0

(:3

0

c3

(:3

X

o

0

c3

(:3

£x]

0

I

o

1.13

!

b..

-6

I

0

o

7"

9"0-

0

I,)

0

0

•£ _

_8
141

Voo

Tm a.x [ to

Vmax t > to

Till ax

P

to t

Figure 4.35: Thickening airfoil

The unsteady case computed here is that of an airfoil thickening from zero thickness to 10% of chord. The airfoil grows in 20 chord lengths of time measured with

respect to the freestream velocity. The airfoil growth is given in Figure 4.35. For this

problem rmax = 0.1 and to = 20. This test case may be used as a model for the flow

about the advancing rotor of a helicopter. The dual potential results are compared

to an Euler solution and transonic small disturbance equation solution of Beam and

Ballhaus (1975) in Figure 4.36. The pressure coefficient as a function of time is plotted at the position z/c = 0.525 measured from the leading edge of the airfoil. The

dual potential code in Euler mode computes the unsteady solution for the thickening

parabolic arc airfoil in best agreement with the small disturbance computation. All

results in the figure are for a domain that is three chord lengths long and extending

two chord lengths out from the bump. The grids are nearly alike. The DP code had
142

61 streamwise by 41 transverse points all uniformly spaced. With a 121 x 81 grid on

this domain the DP solution is only slightly lower (no more than 1.3%) with a more

pronounced overshoot. Beam and Ballhaus used 25 points on the airfoil. No claim

is made that this solution is grid independent. In fact, further calculations using the

DP method on refined grids and on a domain that is five chords in the streamwise

and four chords to the far field have given pressure coefficient results that asymptote

about 8.5% below the results in Figure 4.36. Though the DP Euler solution is short

of the overshoot of Beam and Ballhaus, it is close to the transonic small disturbance

(TSD) equation result for this case.

The solution for this unsteady problem converges slowly at first as the scalar potential is computed from an initial field of $ = 0. Later, during periods of slow change,

the solution is obtained rapidly. When the bump grows the fastest the solution speed

is slowest.

Another result is shown in Figure 4.37. For this test the parabolic arc airfoil

grows to 10% thickness in t = 15 chord lengths. The dual potential solution is

compared to a solution of the linearized transonic small disturbance equation. The

pressure coefficient is plotted for z/c = 0.525 as measured from the leading edge

of the airfoil. Early unsteady calculations made it to steady state by reducing the

tolerance on the dilatation, B. The boundary conditions for the unsteady problem

are actually fixed for all time except at the bump. Thus, unsteady terms are not

required in most boundary conditions at the far- and near-field. A tight tolerance on

the scalar potential is important here.
143

LD

_-)o

C:)×

o

O

(DE

...)

__C:

...)

C1_

C.D

¢+_

C) C

C)

..9O

...)

¢_

C_

IEO

.+.>

'N

t_

00

C_I

..DC

..D

<.,__

r-

._>

l"0 'I

I | .'

"OI Q)

,._' ol

f.'CI

,,_41 "_L:

Xl-_l

O '-" ffJI Z _ .-_"

.J O_ I_l

_.- c: oi C II1, .-J.

-_, Cl

\\ o, OI

\\L W, _I

\'.Gi CI

\\,,: a., cl I I I I I I

9"0 5"0 _'0 £'0 _'0 I'0 0"0 d Ooo

-t_

I'No

m. o

•

••>0

m..o

ffl

rm- DO

f-O

m [---, o

u-; oo

•"_ II _ 8

o_,,_

0_,_
144

00

O)

zy)

7D00

O90C0E

O_0

.__) 00

7:3

O_

C3

00000

.°

IL8

°F,,4

°_,,q0

19 ,,_ _ m

t-I0 "_

°_,_
145

4.3.2.3 Steady viscous flow

Variable property channel flows Only a few results will be shown

to demonstrate the capability of the dual potential code. Clearly, an infinite number of different wall boundary conditions exist. The full Navier-Stokes solution for

variable property flow in a two-dimensional channel could be computed for any wall

conditions of interest. The computation of variable property channel flows has not

been wholly successful. High wall temperatures and high wall heat fluxes lead to a

computed flow separation near the inlet and subsequent reattachment and flow development. The addition of energy to the flow gives a lower "fully-developed" Nusselt

number than the constant property cases. Also, the Nusselt number does not asymptote as it does in the constant property case. Rather, the Nusselt number continues

150 will be presented for the constant wall

case with _ = 1.1. The higher Reynolds number is more difficult temperature

to compute. Contours of temperature, u velocity and dilatation are shown in Figures 4.39, 4.40 and 4.41 for the Re = 40 case. Similar plots for Re = 150 are shown

in Figures 4.42, 4.43 and 4.44. The Nusselt number development for these eases is

plotted in Figure 4.45 along with results from Schade and McEligot (1971) and a

solution from the code of Nelson and Pletcher (1974). Both of these references used

the boundary-layer equations.

For a wall temperature of _ = 2.0 the flow separates near the inlet and redevelops. The separation is evident in the contour plot of u velocity in Figure 4.46.

This case does not converge. A high wall heat flux will also not yield a converged

to drop for wall heating.

Results for Re = 40 and Re =
146

solution. A non-dimensional heat flux, Q, is defined as

II

Q_ q°Dhyd (4.35)

koTo

A near inlet separation develops as shown in Figure 4.47 for Q = 25.0. Lower wall

heat fluxes do yield converged solutions. A non-dimensional wall heat ftux of Q = 0.5

was computed for Re = 150 and M = 0.1. Contours of the temperature, u velocity and

dilatation are shown in Figures 4.48, 4.49 and 4.50 respectively. The Nusselt number

was again below constant property predictions and does not asymptote but continues

to fall. This behavior is expected for solutions which include viscous dissipation (Shah

and London 1978).

The difficulties in solving this problem have been isolated to the density and

dilatation determination. Let us consider the density first. Density is computed from

the continuity equation. Since that equation has first order time, a: and y derivatives,

only one condition on time, x and y are needed. In Figure 3.3 for the boundary

conditions of this problem, however, note that density is constrained by an inlet and

two wall conditions. An initial condition is also assumed. This is one condition too

many. The specification of the wall density at both walls seems to be the problem.

Though many methods of determining the wall density are possible, none seemed to

be completely satisfactory. The wall density has been determined alternatively from

1. the x momentum equation solved along the wall,

2. continuity solved at the line of points just above the wall followed by a statement

that py = 0, or

3. density equals a constant, p = _.
147

Of these, number two gives a constant pressure on a transverse plane, but a realistic

looking negative streamwise pressure gradient has not been obtained. Mass is conserved for all of these approaches. All approaches give the same Nusselt number for

the conditions computed. So, at least for these low Mach number, and "low" Reynolds

number flows, the density doesn't seem to interfere much with the temperature field

determination.

Next, consider the dilatation variable, B. The dilatation transport equation

evidenced early on that it would be difficult to solve. Many approaches were tried:

1. Divergence theorem constraint on the B field

2. Under-relaxation

3. Weighting the B field using the transport equation solution and the definition,

B = uz + v v from the velocities.

4. Smaller time steps.

A useful check is to make sure that the solution for the dilatation satisfies the

divergence theorem (also called Green's theorem or Green's identity). Only twodimensionM compressible problems in Cartesian grids were studied so the following

derivation will be useful for all the compressible test cases.

From the Poisson equation for the scalar potential, Equation 2.13:

V2¢ = B (4.36)

Thus;

ffD V2¢dA = ffD B dA (4.37)
148

Applying the divergence theorem yields:

V2¢dA= D V ¢. n ds = D

f fDBdA=ffD dA=/oD (4.39)

where D is the dom_n, c3D is the boundary of the domain and n is the outward

unit normal. The implication of the above is that the scalar potential should be used

to handle the throughflow velocities at domain boundaries. This is exactly how the

scalar potential boundaries are imposed (except for the boundary-layer calculations

where the vector potential is used for the throughflow velocity). A representative

two-dimensional geometry is shown in Figure 4.38. The outward normal derivative is

indicated on the figure. Proceeding counterclockwise around the domain and applying

the divergence theorem, the area integral of B must satisfy:

ffDBdA = - foXlvdx + foYludy

+ -v dx + u dx (4.40)

1 1

Writing the above integrals in the positive coordinate directions gives:

ffDBdA = - lvdx+ ludy

+_lvdx-_ludx (4.41)

The area integral of dilatation and the boundary integrals of the throughflow velocities

are easily computed to provide a check on the dilatation field.

This constraint was observed to be satisfied automatically to acceptable tolerance

for all the compressible problems except the viscous variable property channel flows.
149

Y

Yl

ny=l

ds = -dz

nx = -1

ds = -dy

D

.__... OD

nz=l

ds : dy

D

0 $ z 1 z

ny = --1

ds = dz

Figure 4.38: 2-D solution domain for divergence theorem application
150

In that case,convergedsolutions can best be obtained by enforcingthe divergence

theorem constraint on the computedB field.

In spite of these problems, solutions were obtained without applying any simplifying assumptions (i.e., no terms were dropped in the governing equations). The

results here were obtained by method 2 for the density boundary condition and the

divergence theorem constraint applied to the solution for B.

Some trends that were observed will be reported next. For variable property

flow, the Nusselt number does not asymptote but continues to fall for wall heating.

This is in agreement with the findings of others. Shah and London (1978) indicate

that the limiting Nusselt number is zero for variable property flow with wall heating.

This is reason for confidence in the trends computed by the dual potential method.

Dilatation has no choice but to conform to the numerical satisfaction of Green's

theorem. This constraint must be used throughout the calculation. The correction

becomes smaller and smaller as the steady state is approached, but without the

correction the convergence is slowed and seemingly halted. At this point the only

remaining problem appears to be getting a believable negative streamwise pressure

gradient. The finest grid solution presented here is for 41 × 41 points. For this coarse

grid the skin friction asymptotes to CfRe = 23.6. The fully-developed skin friction

for an incompressible channel case would be CfRe = 24.0.

The code runs at 83 MFLOP for full Navier-Stokes calculations. The low wall

32 #s

heating cases converge at _ global iteration × grid point"
151

Constant wall temperature Converged solutions were obtained for _ = 1.1. Results for Re = 40 and Re = 150 will be presented. The

flow Mach number is M = 0.1. Contours of temperature, u velocity and dilatation are

shown in Figures 4.39, 4.40 and 4.41 for the Re = 40 case. Similar plots for Re = 150

are shown in Figures 4.42, 4.43 and 4.44. Recall that the physical distances, z and

y, for the two-dimensional channel cases are non-dimensionalized by the hydraulic

1

diameter. The non-dimensional wall separation distance is then ymax = _ so that

Y is in the range 0 < y < 1⁄2 for the results plotted here from wall to wall. The it

velocity contours in Figures 4.40 and 4.43 are similar to results for an incompressible channel case. Notice that the core velocity attains u = 1.5 which is the exact

value for the incompressible case. The temperature contours in Figures 4.39 and

4.42 are realistic considering that the wall temperature is a constant _ = 1.1.

The dilatation contours indicate that most of the compressible effects for this flow

are concentrated near the inlet. The Nusselt number development for these cases is

plotted in Figure 4.45 along with results from Schade and McEligot (1971) and a

solution from the code of Nelson and Pletcher (1974). Schade and McEligot used the

boundary-layer assumptions and neglect viscous dissipation in the energy equation.

They compute an increase in Nusselt number for heating when properties are variable.

The computations of Pletcher are also based on the boundary-layer equations but include the dissipation terms in the energy equation. His results confirm the behavior

predicted by the dual potential solution of the Navier-Stokes equations. Shah and

London (1978) also confirm that the Nusselt number does not asymptote for variable

property conditions with heating.
152

(_ o o 6 o

>-

0

X

O

0

0

0

0

N

0

0

0

0

0

o

°,-_

o

;>

o,_
153

0

0

0

u_

o

oq,q

0

t-c "

II

o

o_.._

>.
154

x

_. _ N, -. o

0 0 0 0 0 0

I

(/J

0

°1.-i

0

o

bO

°_
155

0

(#3

v

o_,._
156

,,--4

0

i¢

0

•._ II

H

°_

o

o_

&-i
157

O

u_

,4

I@

X

N

N

O

0

0

0

°_

0 •

_ v

°_
158

(9

.4_) •

o ?

(l) .,

E •

E-_ .J

c- c

.+.) ._

C o

0 (-

E__ .,

r,,o ._

(1) -o

._ t_

o X

_ 0

O_ c. O_

._)

___)

___)

0

(_

0

n

0

0

>

0"05

01 "-I _ I

I BI • •

I1,'

,-,,-, L 8

I ...), .-_ to,o-_T'

o -',

_N_2

Z (_ ._),X X -,Jl

K..,..JQ3 L,-lt..-.",','-' O)

o o,_r _r -_I

r,l c T- I ...I (_ 't .0..I

C "10, ",0 731

C,O Lr) Cl

0,_I _ "_ Ol

"C t I I I

.,-,-0,0_ O,e r,-" @ (_I

OlO- Z

ff')lO 0

I

,_ +

I

I

i d

I

/

i

O'Oi,

i

0"0£

i

O'OE

l'XnN

I

0"0!

0

0

?

_0

T

0

)K

('M

0"0

)K

X

.=

o

o

7:

o

f.4

.o

°_

_0
159

00

N

\\\\

0

In0

ut_

'2,

X

mt_ o0o0

00

0_

0

0,.q
160

T ,,

Figure 4.46 is from the results of the case with _ = 2.0. This case does not

converge. These results are shown mainly to show the separation that is computed

by the DP code. For convenience, the constant wall temperature is given as Twall in

the figures. It is to be understood that it is actually the ratio of _.

-lnle_

Constant wall heat flux Figure 4.47 is for the case with

Q = 25.0. This case also does not converge. These results are shown mainly to show

the separation that is computed by the DP code. The high heat flux and high wall

temperature cases exhibit this same feature.

A converged solution is obtained for lower heat fluxes. Results are plotted in

Figures 4.48-4.50 for Q = 0.5, Re = 150 and M = 0.1. The non-dimensional wall

separation distance is Vmax = 1⁄2 since V is non-dimensionalized on the hydraulic

diameter which is simply twice the wall spacing.

Compressible boundary layer

Subsonic freestream The compressible subsonic flow over

a flat plate was computed to test the code on an external, viscous, compressible flow

case. The conditions at the plate are adiabatic. The normal pressure gradient at the

plate, py, is assumed to be zero. This is consistent with a boundary-layer assumption.

The boundary conditions for this case are shown in Figure 4.51.

The full dual potential method and a variation were used to compute the solution.

The variation on the dual potential method was used as a self check on this problem.

It uses dilatation computed from the continuity equation and density computed from
161

0

I

0

I

e4

0

>..

x

url

o

o

o

0

0

0

c_

_8
162

0

0

U

ol, q

0

Q::

0

o_

°_..q
163

in

o

X

N

N

mtj

o

0

0

0

o

o_

H

o_

Q
164

0

q)

o

°,.'4

Q)

0

Q) LL

_qX

°_

°°

°,.._
165

V, J index

inlet:

w=O

B=O

T=I

p=l

A=O

¢=0

u=l

v#O

freestream:

w=O

B - -[upz+vPV) =_ 0

IX

P

T=I

p=l

Az = Cy - v

Cy = 0

u=l

Vy : 0

exit:

tt;XX = 0

Bxz : 0

Txz = 0

Px upwinded

Azz = 0

_XX = 0

VX x = 0

_\\\\\\\\\\\\\\\\\\\\\

flat plate:

w = f(A, c#) x, i index

B=O

OT = 0

p = fixed by py = 0 and T distribution

A=O

CV _ 0

u=0

V=0

Figure 4.51: Boundary conditions for laminar compressible flow over a flat plate
166

CD

CD

LD

CD

CD

T---4

G)

C_

O Ln

___) T-

._)

(_+_

O

l_

0_

D')

+)

._)

o

O

G)

0Z 0"5 0"I, 0"£ 0"_

( x/nu/ju717 ) "_Jbsx_

0

II

8

II

0

=

0
167

L

0

3q

0

C --J

0 :n

-J b

D E

_ g

..D n

L

03 c

.D _

-13 o

J

(1) •

L t.

D_ g

E o

©

O

I

90" I

I'-"

O,OD

.a.i r,o_

r7 ._i ._Jlr7

Z

C9 01 r-.,r-.

-J .--_OlC3

I III

L. LIL

@_ Q_lO_

I

SO" I

I

I

t

/

I I

7 I

/,

till_

//

,i ,

EO" I

1

I,O"

I

£0" 1

I

I0"I

I

00" I 66 0

o

t,,,.

if')

o

to

(_

u")

If)

o. x

"_ •

n"

U')[--,

(D

U')

o X

• X

If)

o

u')

o

u')

o

o

II

8

t,_

o

o

°,"4

.£

o,,,_

,--j

LO

I:W)
168

the ideal gas law as follows. From continuity, and assuming steady state conditions,

B = -(upx + Vpy)

P

It is possible to assume that the pressure is constant and uniform for this case, so

that density can be computed from the ideal gas law by

p 1 where p

P-RT-T _=1

This variant is possible because of the simplifying assumptions of this problem. This

gives a much faster solution than the transport equations. When the transport equations are solved for dilatation and density, the boundary conditions shown in Figure 4.51 are used. The wall density can actually be obtained by other methods as

mentioned in the variable property channel flow discussion on page 146. A comparison of the two solution methods and results from a boundary-layer finite-difference

scheme (Christoph and Pletcher 1983) are shown in Figure 4.52 for a M = 0.5 flow.

Pletcher's data were computed using the above mentioned finite-difference boundarylayer scheme with approximately 50 points in the boundary layer. For reference, the

Blasius profile and a solution at M = 1.0 from Schlichting (1979) are included in Figure 4.52. It can be seen that the compressible boundary layer thickens as the Mach

number increases.

The temperature profile for this case is shown in Figure 4.53. The results from

the boundary-layer finite-difference scheme are not plotted but coincide with the

Pr = 0.7 theoretical curve. More grid points in the dual potential solution could

be required for better agreement. However, the results are within 1% of theory. In

the present calculations, the wall temperature is computed by using the zero wall
169

derivative condition for temperature as a boundary condition on the energy equation.

Pletcher reported setting the adiabatic wall temperature in his calculations. The

wall temperature is 1.042 by theory and 1.039 as computed using the dual potential

method.

This was the first compressible viscous case to be computed successfully. There

was initially a problem computing this flow using the full dual potential equation

set. The dilatation transport equation could only be converged by dropping the

terms which originated from the pressure terms in the primitive variable momentum

equations. (Dropping those terms is equivalent to assuming that the pressure is

everywhere constant.) This led to the determination that those terms (second order

derivatives in the B equation source term) must be differenced conservatively. This

is necessary because the pressure gradient is zero, yet numerically it is non-zero for

p = pRT substituted into the pressure derivative terms.

These cases were computed by the full dual potential method without using

boundary-layer assumptions other than py = 0 at the plate. The full energy equation

was used. Streamwise convective terms were upwinded. The time step limitation is

controlled by the dilatation variable. The skin-friction development is nearly indistinguishable from the incompressible boundary layer solution for an irrotational inlet.

The skin friction asymptotes to a value slightly lower than for incompressible flow

with the difference in the fourth place behind the decimal point (cf. Van Driest 1952).

The asymptotic value computed here was Cfv/-Re = 0.6646 at Rex = 100,500. At

the same Reynolds number the incompressible solution computed Cf_ = 0.6647.

With 81 transverse points it is observed that the vorticity magnitude is less than
170

10 -5 from j = 35-81 and the dilatation magnitude is less than 10 -5 from j = 54-81.

The solvers for vorticity and dilatation can be shut off beyond these points with no

adverse affect on the solution and a speed up of approximately 14%.

Supersonic freestream A supersonic flow over a flat plate

was also computed. The conditions at the plate are adiabatic. The result in Figure 4.54 was obtained by computing the dilatation from the continuity equation. A

zero streamwise pressure gradient was assumed. Then, by the boundary-layer assumption, a zero normal pressure gradient is also assumed so that pressure is constant for

this problem. The dilatation may then be computed from the continuity equation

B= - + py) (4.42)

P

and density can be computed from the ideal gas law and the constant pressure assumption

p 1 for a unit pressure field.

It was not possible to compute this flow using the transport equation for B.

(4.43)

The abscissa in Figure 4.54 is _2 where _2 is the momentum thickness defined

by

_2 = =0 pcvUoo f7c¢ dy

Experimental results for this case were obtained by O'Donnell (1954). The theoretical

calculation is by Chapman and Rubesin (1949). The achievement of good results here

may be due in part to the use of the fixed pressure assumption.
171

CDcD

CD

CD

c_o

fl -oCo OC

t_

O .._E

__.)

.-_ __,_ _ ..o 0 .$

0_•L

_,_)o

...)_

0 ,--, o

t

C@

..0

-0C0

QC

Z0

WE

00.

6.)0

-J .C

_r

(.

J

CC --_ 0 • Q C • .J C0 CE .0

.-_L 03

0@

@ [I. ", 3S X

_._@ r7

0

I

0"!

fI I I

8"0 9"0 I,'0 _'0 0"0

Ju7rT/n

0

0o

I--4

c_00

c_

0

o -0J

0o•

-0

o0oo

c_

@

_m@@

,.o

,.m

_3

_Q@@

_mO

o_-_
172

a

Figure 4.55: 3-D duct geometry

4.4 Three-dimensional Cases

Only a first step to three-dimensional calculations using the dual potential method

has been taken here. For incompressible flow, solutions were computed for the developing flow in rectangular ducts of constant cross section. This may be considered as

a stepping stone to more complex geometries. Compressible flow solutions in threedimensions have not been obtained. Time did not permit further development of the

dual potential method, but the capability of this method has been demonstrated for

three-dimensional incompressible flow calculations in a simple geometry.

4.4.1 Incompressible flow

4.4.1.1 Incompressible channel inlet A three-dimensional, laminar, incompressible flow code has been programmed using the dual potential method. The
173

formulation (Wong and Reizes 1984) is applicable to ducts of arbitrary but constant

cross section. The code has been used to calculate the developing flow in rectangular

ducts of various aspect ratios for Reynolds numbers of 10 and 50. The duct geometry

is shown in Figure 4.55. Uniform Cartesian grids were used to compute the solutions

to be presented.

For the three-dimensional duct geometry shown in Figure 4.55, the aspect ratio is

b/a. White (1974) and Shah and London (1978) give a formula to compute the fullydeveloped velocity profile for constant property incompressible flow in a rectangular

duct of constant cross-section.

For the streamwise direction z and a rectangular section with -a < y < a and

-b < z < b:

#1r3 -_--_ _ (-1)( n-l)/2 [1 - cosh (nTrz/2a)] (nTry_

n=l,S,5,.., cosh(n_rb/2a)J cos \ 2a }

4ba3 (-_-_) 1 (b) _ 1 5 _ _ _-_ tanh n=1,3,5,... J

Q __

Where Q is the volume flow rate, Q = ureA and A = 2a × 2b.

Note that the viscosity It does not matter for the fully-developed velocity profile

of u

_-mm" The pressure gradient may be eliminated in the equation for u so that u_u___ Um

may be determined readily from the above equations.

The computed flow development and fully-developed profiles agree well with

the known results. Better agreement could be obtained using stretched grids. The

centerline velocity development for a square duct is shown in Figure 4.56. Computed

results using the dual potential method are indicated by symbols. The centerline

velocity development as computed on stretched grids by Wong and Reizes (1984)
174

Table 4.2: Fully-developedcenterlinevelocity and skin friction

Umax/Um C/Re

Aspect ratio Re[dual pot. grid ref.

1.00 10 15 × 15 x 45 2.096

1.00 50 29 × 29 × 60 2.096

0.75 50 21 x 15 × 30 2.077

0.50 50 31 × 15 × 30 1.992

0.25 50 61 x 15 x 30 1.774

I dual pot. ref. dual pot.

2.065 14.22708 14.07

2.084 14.22708 14.23

2.051 14.47570 14.38

1.968 15.54806 15.40

1.752 18.23278 17.90

is shown as the Re = 10, 50 and 200 results. Experimental results of Goldstein

and Kreid (1967) would be very near the Re = 200 results. The fully-developed

velocity profile for Re = 50 is compared to the analytical solution in Figure 4.57. The

transverse velocity vectors near the exit of a square duct are shown in Figure 4.58 and

near the exit of an aspect ratio = 0.50 rectangular duct in Figure 4.59. The velocity

vectors are magnified 100 times to show the persistent vortical flow in the corners

even at the fully-developed condition. A similar symmetric pattern is exhibited in all

the rectangular ducts. The transverse flow is toward the corner along the walls and

away from the corners along the corner bisector. Note that this flow pattern happens

to be the opposite of the Reynolds stress driven secondary flow for a fully-developed

turbulent flow in a rectangular duct (Demuren and Rodi 1987; Speziale 1987a). The

flow patterns here probably result from initial disturbances that have not yet decayed.

Stretched grids are needed for more efficient calculation of this flow.

A summary of the fully-developed centerline velocity and wall skin friction coefficient for various aspect ratio ducts is given in Table 4.2. The tolerance used in

computing the dual potential results was 0.0001. These results were computed on a

uniform Cartesian grid. For the square cross-section duct, the skin friction coefficient
175

.+3

C

(1)

E

OO

._.3

>

(9

r-)

:::o

..,o

...3

O

O

.__)

(1)

(1)

C

_._.)

(__

C

C._

I

0"_

0

I I

I I I I

O'I 9"1 t,'l _'[

uoewn/×o_n

0"!
176

S'0

0

I I I I

I,'O Z'O Z'O ["0

z

II

hO

X

X

o_,,_

u_

oo

°_,,_
177

\\1 i/.......,.,,. \ _.

,,\_

------'_

..: ,,

I

;:¢y.,';,2zI': ',t

",'//"1" el ' . ..../z l

/,.----',

i

i

_ \ I / / t'., t

- /

/

,,_t "/1 p

/ /

_ \ \.,. _II

\

Figure 4.58: Transverse velocity vectors at a plane near the exit of a square duct

with Re = 50
178

Figure 4.59: Transverse velocity vectors at a plane near the exit of a 0.50 aspect

ratio duct with Re = 50

asymptotic value is computed to be CfRe = 14.23 on a 29 x 29 x 60 grid. This is the

exact answer to two places behind the decimal point.

The shear force was computed by integrating the product of the first order velocity derivative and the local wall area over the entire wall surface. A second order

polynomial fit of the u velocity component is used to compute the velocity derivative

in the shear stress equation. The local wall shear stress is computed from:

rxy = g +

rzz = # (0z+ 0z)

ryz = . N+

The shear stresses vxy and rxz are used in this computation. The integrated wall

shear force is then used to compute CfRe.

The three-dimensional code has not been vectorized and uses point Gauss-Seidel

with SOR for the vector potential Poisson equation. That particular scheme becomes
179

more and more inefficient as the number of grid points increases. The computation

time can be reduced somewhat for steady state calculations by tightening the tolerance gradually as iterations increase. For a uniform Cartesian grid, 15 x 15 x 30

points is a good compromise of accuracy and solution speed. It has been observed

that a solution tolerance of 0.001-0.0001 is adequate. For more accurate results on a

uniform Cartesian grid at the Reynolds numbers reported here, the spacing as in the

square duct case with 29 x 29 x 60 points should be used. This gives fully-developed

flow solutions within 1% of the exact. Note that for incompressible flow the scalar

potential can be solved once and for all for a given fixed geometry. This feature is

used in the solution presented here. For the rectangular duct geometry, ¢ = x is the

analytical solution to the Laplace equation V2¢ = 0 with the boundary conditions

discussed in Section 3.2.1.

74Fs

The computation rate for this unvectorized code is _ global iteration x grid point

on a single processor of the Cray X-MP. The complete flow field is solved by the current code. The symmetry of this problem was not used to speed the solution. This

would immediately reduce the computation time by about a factor of four. The

MFLOP rate is 12 MFLOP overall, without any enhancements.

The three-dimensional incompressible code should be extended by adding the

energy equation, generalized coordinates (so a curved duct case can be computed)

and a rotating coordinate capability for centrifugal compressor modeling. A more

efficient Poisson solver is needed and other enhancements to improve the computation

speed.
180

4.4.2 Compressible flow

The transport equations for w and B in the three-dimensional compressible

formulation are given in Appendix A. No solutions were attempted.
181

5. DISCUSSION AND CONCLUSIONS

A dual potential procedure has been developed and evaluated for flow calculations

ranging from potential flow to full Navier-Stokes solutions in two dimensions. A threedimensional incompressible formulation has also been presented.

Compressibility has been handled by the dilatation variable. A dilatation transport equation is obtained from the momentum equations. Workable dilatation boundary conditions have been presented.

From this study the following conclusions can be made:

1. The dual potential formulation of the Navier-Stokes equations is very flexible

in that it is easy to compute subset equations (potential, Euler).

2. The dual potential method can simulate irrotational and rotational inflow easily.

In fact, it is interesting to note that the test cases for which experimental results

are available agreed best with the computed results for an irrotational inflow

condition--boundary layer and three-dimensional duct inlet flow. Yet many

of the computations had to be compared with the rotational inlet condition as

this is what other computational investigators have used, namely in the channel

heat transfer cases. Van Dyke (1970) reports that the irrotational inlet condition

may be more realistic for experimental models with a rounded entrance. The
182

above statement and the small sampling of experimental results suggests that

the irrotational inlet condition is an important feature to be able to simulate.

The dual potential method handles either condition easily.

3. Computational effort can be reduced in irrotational regions (_= 0) and in

incompressible regions (B = 0).

4. The method appears to be very accurate for incompressible calculations.

5. It appears possible to extend the dual potential method to compressible flow,

but the dilatation transport equation may need more work. It seems that some

problem still persists here--perhaps self consistency of the inflow B boundary

condition--that is most noticeable for viscous compressible problems. Also,

the dilatation can undergo much variation in certain regions of the flow field

making it difficult to resolve. For example, in stagnation regions B changes sign

as shown in Figure 4.32 for the flow over a parabolic arc bump. This behavior

of the dilatation variable may be undesirable in practical use of this method.

6. The scalar potential requires a tight convergence tolerance as used in the test

cases here. This may be explained by the fact that the scalar potential is the

largest component in the velocity decomposition for the streamwise velocity.

The convergence tolerance on the vector potential can be 1-2 orders of magnitude less restrictive than for the scalar potential. This observation may be

problem specific since the streamwise velocity component was the dominant

component for the test cases studied in this work.

7. For incompressible flows, the dual potential approach can be very competitive
183

with primitive variable methods. As an example, the original INS3D code

100 tts

(Kwak et al. 1986) performs at _ global iteration x grid point on the Cray

X-MP. The non-optimized version of the three-dimensional incompressible dual

74 #s

potential code runs at _ global iteration × grid point on the X-MP. The only

significant point to make about this comparison is that the dual potential code is

written for uniform Cartesian grids while the INS3D code can handle generalized

three-dimensional coordinates.

8. It is possible to solve the dual potential equation set in an iterative, uncoupled

way.

A number of capabilities of this formulation have been demonstrated, some for

the first time. The success seems to be limited for compressible viscous flows with a

pressure gradient. Determination of the wall density is awkward for these cases and

perhaps causes some problems for the dilatation transport equation. The dilatation

transport equation seems to work fine for a known density and temperature field as

was the case for the boundary-layer solutions. More remains to be done before a

trouble free compressible formulation is available.
185

6. RECOMMENDATIONS FOR FUTURE WORK

The key advantages and disadvantages of the dual potential method have been

stated in Section 1.3.2.2. Speziale (1987b) describes further advantages of nonprimitive variable methods using vorticity. He shows that these methods have simpler

boundary conditions than primitive variable formulations for problems in a rotating

reference frame. A review of the literature (Ozoe et al. 1985) suggests that a dual

potential formulation is the method of choice for confined flows. The solution speed

of the dual potential type methods for incompressible flows is well documented (Aziz

and Hellums 1967). However, the disadvantages have precluded solutions in complex

geometries.

It seems best then to concentrate research efforts on the known advantages of this

method and solve problems of practical interest. A few problems will be suggested

that exploit the advantages of the dual potential method. Global weather forecasting

is one such interesting practical problem. The geometry is simple, the flow is confined

and density variations can be obtained from simple correlations. Incompressible flows

inside curved, twisted ducts have been solved using the dual potential method (Yang

and Camarero 1986). It should be possible to build on these solutions by allowing

the duct to rotate. Many practical fluid flow problems could then be modeled such

as the flow through an automotive water pump, a centrifugal compressor, a flow

i_E.CEDiI_ PAGE BLANK NOT FILMED

I_.__][N]'iO/IJti_[ BLANK
186

meter, etc. Further interesting practical solutions could be obtained for problems that

resemble the bump cases studied here. The "walls" in that two-dimensional problem

can be allowed to move and the bump can have a time-dependent growth. Also,

the flow can be treated as viscous. Since solutions have been obtained for subsonic

compressible flows, the dual potential method may be well suited for computing flows

with embedded compressible regions. Examples of practical flows with these features

are the flow under an automobile and flow through a carburetor. Heat transfer effects

can be included in all of these simulations.

In other research areas the dual potential method may also be of importance.

Direct simulations of channel flow and a flat plate boundary layer have been obtained

by Rai and Moin (private communication, NASA/Ames Research Center). They

used a non-conservative primitive variable method on a staggered grid to simulate

incompressible turbulent flow. The dual potential method could be used to compute

direct simulations of these flows with little modification. If the speed advantage

of the dual potential method for incompressible flow calculations extends to direct

simulations, then interest in the formulation will surely increase.

It is always wise to seek out similarities in the other disciplines. The Helmholtz

decomposition theorem was developed and used extensively by scientists working in

electrodynamics. There may be some analogous physical problem in physics, mathematics or other field that could help in the understanding and use of the dual potential

method. A search of the literature quickly turned up the paper by Miiller (1987) on

the topic of vector splitting applied to a physics problem.

Since the dependent variables in the dual potential formulation are vorticity and
187

dilatation, one should use these variables to his advantage not only in computing a

flow, but in the understanding of the physics. For example, vorticity is generated

by such natural phenomenon as walls and shocks. A code which solves for vorticity

directly, such as the dual potential method, may be useful in problems in which one

wishes to track the vorticity. The dilatation variable can also be useful in understanding some flow situations. The dilatation changes sign in the vicinity of stagnation

points and shocks. This physical behavior can help give insight into the problem

being studied (and cause problems numerically!).

Finally, it is important to have a code that is a good Navier-Stokes solver, but

also is easily and efficiently used in other modes. It has been demonstrated that the

dual potential method can be used in a potential, Euler and Navier-Stokes mode.

Further work can improve the performance even more and address the issues of shock

capturing, turbulence modeling, direct simulations, etc. If this method is ever to

prove useful in applications, it must be able to run in all these modes efficiently and

reliably. In addition it must be faster than existing methods. Unfortunately, the

ideas behind the dual potential method are not as familiar as the primitive variable

approaches. This unfamiliarity is a handicap since the details of application codes

should be easily understood. To increase the familiarity of this method it must be

tested on some problems of practical interest.

The two- and three-dimensional dual potential codes written as a part of this

research are still under development. A copy of a version of the codes may be obtained

from the first author or through the Department of Mechanical Engineering, Iowa

State University, Ames, Iowa 50010 (Attn: Professor R. H. Pletcher).
189

7. REFERENCES

Ames, William F. (1977). Numerical methods for partial differential equations. 2d

ed. Orlando: Academic Press.

Anderson, D. A., Tannehill, J. C. and Pletcher, R. H. (1984). Computational fluid

mechanics and heat transfer. New York: McGraw-Hill.

Aregbesola, Y. A. S. and Burley, D. M. (1977). "The vector and scalar potential for

the numerical solution of two- and three-dimensional Navier-Stokes equations."

J. Computational Physics, 24, 398-415.

Ashley, Holt and Landahl, Marten (1965). Aerodynamics of wings and bodies.

Reading, Massachusetts: Addison-Wesley.

Aziz, K. and Hellums, J. D. (1967). "Numerical solution of the three-dimensional

equations of motion for laminar natural convection." Phys. Fluids, 10, 314-325.

Ballhaus, W. F. and Steger, J. L. (1975). "Implicit approximate-factorization

schemes for the low-frequency transonic equation." NASA TM X-73,082.

Beam, R. M. and Ballhaus, W. F. (1975). "Numerical integration of the

small-disturbance potential and Euler equations for unsteady transonic flow."

NASA SP-347, 2, 789-809.

Beam, R. M. and Warming, R. F. (1978). "An implicit factored scheme for the

compressible Navier-Stokes equations." AIAA J., 16, 393-401.

Blasius, H. (1908). "Grenzschichten in Fliissigkeiten mit kleiner Reibung." Z.

Math. Phys., 56, 1-37. English translation in NACA TM 1256.

Briley, W. R. and McDonald, H. (1977). "Solution of the multidimensional

compressible Navier-Stokes equations by a generalized implicit method." J.

P_J_CEDII_G PAGE BLANK NOT FILMED

I,_HN'[10N/CMI IIL_
190

Computational Physics, 24, 372-397.

Chaderjian, N. M. and Steger, J. L. (1985). "The numerical simulation of steady

transonic rotational flow using a dual potential formulation." AIAA Paper

85-0368, Reno, Nevada.

Chapman, Dean R. and Rubesin, Morris W. (1949). "Temperature and velocity

profiles in the compressible laminar boundary layer with arbitrary distribution

of surface temperature." J. of the Aeronautical Sciences, 16, 547-565.

Chorin, A.J. (1967). "A numerical method for solving incompressible viscous flow

problems." J. Computational Physics, 2, 12-26.

Christoph, G. H. and Pletcher, R. H. (1983). "Prediction of rough-wall skin friction

and heat transfer." AIAA J., 21,509-515.

Currie, I. G. (1974). Fundamental mechanics of fluids. New York: McGraw-Hill.

Davis, R. L., Carter, J. E. and Hafez, M. (1986). "Three-dimensional viscous flow

solutions with a vorticity-stream function formulation." AFOSR-TR-86 (see

also AIAA Paper 87-0601, Reno, Nevada).

Davis, R. L., Carter, J. E. and Hafez, M. (1989). "Three-dimensional viscous flow

solutions with a vorticity-stream function formulation." AIAA J., 27, 892-900.

Dennis, S. C. R., Ingham, D. B. and Cook, R. N. (1979). "Finite-difference

methods for calculating steady incompressible flows in three dimensions." J.

Computational Physics, 33, 325-339.

Douglas, J., Jr. (1962). "Alternating directions method for three space variables."

Numerische Mathematik, 4, 41-63.

D'Yakonov, Ye. G. (1963). "On the application of disintegrating difference

operators." Zh. vflchisl. Mat. Mat. Fiz., 3, 385-388.

E1-Refaee, M. M. (1981). "A numerical study of laminar unsteady compressible

flow over airfoils." Ph.D. dissertation, Georgia Institute of Technology, Atlanta,

Georgia.

E1-Refaee, M. M., Wu, J. C. and Lekoudis, S. G. (1981). "Solutions of the

compressible Navier-Stokes equations using the integral method." AIAA J., 20,
191

356-362.

Farouk, B. and Fusegi, T. (1985). "A coupled solution of the vorticity-velocity

formulation of the incompressible Navier-Stokes equations." Int. J. ]:or Numer.

Meth. in Fluids, 5, 1017-1034.

Fasel, H. (1976). "Investigation of the stability of boundary layers by a

finite-difference model of the Navier-Stokes equations." J. Fluid Mechanics, 78,

355-383.

Fasel, H. and Booz, O. (1984). "Numerical investigation of supercritical

Taylor-vortex flow for a wide gap." J. Fluid Mechanics, 138, 21-52.

Gatski, T. B., Grosch, C. E. and Rose, M. E. (1982). "A numerical study of the

two-dimensional Navier-Stokes equations in vorticity-velocity variables." J.

Computational Physics, 48, 1-22.

Gentzsch, W. and Neves, K. J. (1987). "Computational fluid dynamics: Algorithms

and supercomputers." In AGARDograph No. 311, ed. H. Yoshihara.

Giannakoglou, K., Chaviaropoulos, P. and Papailiou, K. D. (1988). "Computation

of rotational transonic flows using a decomposition method." AIAA J., 26,

1175-1180.

Goldstein, R. J. and Kreid, D. K. (1967). "Measurement of laminar flow

development in a square duct using a laser-doppler flowmeter." Trans. ASME

Ser. E, J. Appl. Mech., 34, 813-818.

Guj, G. and Stella, F. (1988). "Numerical solutions of high-Re recirculating flows

in vorticity-velocity form." Int. J. for Numer. Meth. in Fluids, 8, 405-416.

Hafez, M. M., Habashi, W. G., Przybytkowski, J. L. and Peeters, M. F. (1987).

"Compressible viscous internal flow calculations by a finite element method."

AIAA Paper 87-0644, Reno, Nevada.

Hirasaki, G. J. and Hellums, J. D. (1970). "Boundary conditions on the vector and

scalar potentials in viscous three-dimensional hydrodynamics." Q. Appl. Math.,

28, 293-296.

Hoist, T. L. (1987). "Numerical solution of the Navier-Stokes equations about

three-dimensional configurations--A survey." Supercomputing in Aerospace.
192

NASA CP 2454, 281-298.

Hwang, Ching-Lal and Fan, Liang-Tseng (1964). "Finite difference analysis of

forced-convection heat transfer in entrance region of a flat rectangular duct."

Appl. Sci. Res., Sect. A, 13, 401-422.

Kays, W. M. and Crawford, M. E. (1980). Convective heat and mass transfer. 2d

ed. New York: McGraw-Hill.

Kwak, D., Chang, J. L. C., Shanks, S. P. and Chakravarthy, S. (1986). "A

three-dimensional incompressible Navier-Stokes flow solver using primitive

variables." AIAA J., 24, 390-396.

MacCormack, R. W. (1981). "A numerical method for solving the equations of

compressible viscous flow." AIAA Paper 81-0110, St. Louis, Missouri.

Mallinson, G. D. and De Vahl Davis, G. (1973). "The method of the false transient

for the solution of coupled elliptic equations." J. Computational Physics, 12,

435-461.

McDonald, J. W., Denny, V. E. and Mills, A. F. (1972). "Numerical solutions of

the Navier-Stokes equations in inlet regions." Trans. ASME Set. E, J. Appl.

Mech., 39, 873-878.

Mitchell, A. R. and Griffiths, D. F. (1980). The finite difference method in partial

differential equations. New York: John Wiley & Sons.

Morino, L. (1985). "Scalar/vector potential formulation for compressible viscous

unsteady flows." NASA CR 3921.

Morino, L. (1986). "Helmholtz decomposition revisited: Vorticity generation and

trailing edge condition." Computational Mechanics, 1, 65-90.

Mfiller, E. E. (1987). "Scalar potentials for vector fields in quantum

electrodynamics." J. Math Phys., 28, 2786-2790.

Nelson, R. M. and Pletcher, R. H. (1974). "An explicit scheme for the calculation

of confined turbulent flows with heat transfer." Proceedings of the 1974 Heat

Transfer and Fluid Mechanics Institute. Stanford University Press. 154-170.

O'Donnell, Robert M. (1954). "Experimental investigation at a Mach number of
193

2.41of averageskin-friction coefficientsand velocity profilesfor laminar and

turbulent boundary layersand an assessmentof probeeffects." NACA TN 3122.

Orlandi, P. (1987). "Vorticity-velocity formulation for high Reflows." Computers _J

Fluids, 15, 137-149.

Osswald, G. A., Ghia, K. N. and Ghia, U. (1987). "A direct algorithm for solution

of incompressible three-dimensional unsteady Navier-Stokes equations." AIAA

Paper 87-1139.

Ozoe, H., Yamamoto, K., Churchill, S. W. and Sayama, H. (1976).

"Three-dimensional, numerical analysis of laminar natural convection in a

confined fluid heated from below." J. Heat Transfer, Transactions of the

ASME, 98, 202-207.

Ozoe, H., Yamamoto, K., Sayama, H. and Churchill, S. (1977). "Natural convection

patterns in a long inclined rectangular box heated from below. Part II.

Three-dimensional numerical results." Int. J. Heat Mass Transfer, 20, 131-139.

Ozoe, H., Yamamoto, K. and Churchill, S. W. (1979). "Three-dimensional

numerical analysis of natural convection in an inclined channel with a square

cross section." AIChE J., 25,709-716.

Ozoe, H., Mouri, A., Ohmuro, M., Churchill, S. and Lior, N. (1985). "Numerical

calculations of laminar and turbulent natural convection in water in

rectangular channels heated and cooled isothermally on the opposing vertical

walls." Int. J. Heat Mass Transfer, 28, 125-138.

Patankar, S. V. (1975). "Numerical prediction of three-dimensional flows." In

Studies in convection: theory, measurement, and applications, ed. B. E.

Launder, 1, 1-78. New York: Academic Press.

Patankar, S. V. (1981). "A calculation procedure for two-dimensional elliptic

situations." Numev. Heat Transfer, 4, 409-425.

Peaceman, D. W. and Rachford, H. H. (1955). "The numerical solution of parabolic

and elliptic differential equations." J. Soc. Indust. App. Math., 3, 28-41.

Rao, K. V. (1987). "A three-dimensional dual potential procedure with

applications to wind tunnel inlets and interacting boundary layers." Ph.D.

dissertation, Iowa State University, Ames, Iowa.
194

Rao,K. V., Steger,J. L. and Pletcher, R. H. (1989). "Three-dimensional dual

potential procedure for inlets and indraft wind tunnels." AIAA J., 27, 876-884.

Richardson, S. M. and Cornish, A. R. H. (1977). "Solution of three-dimensional

incompressible flow problems." J. Fluid Mechanics, 82, 309-319.

Roache, P. J. (1972). Computational fluid dynamics. New Mexico: Hermosa

Publishers.

Roberts, G. O. (1971). "Computational meshes for boundary layer problems."

Proc. Second Int. Conf. Num. Methods Fluid Dyn., Lecture Notes in Physics, 8,

171-177, New York: Springer-Verlag.

Rogers, S. E. and Kwak, D. (1988). "An upwind differencing scheme for the

time-accurate incompressible Navier-Stokes equations." AIAA Paper 88-2583,

WiUiamsburg, Virginia.

Rogers, S. E., Kwak, D. and Chang, J. L. C. (1987). "INS3D--An incompressible

Navier-Stokes code in generalized three-dimensional coordinates." NASA

Technical Memorandum 100012.

Schade, K. W. and McEligot, D. M. (1971). "Cartesian Graetz problems with air

property variation." Int. J. Heat Mass Transfer, 14,653-666.

Schlichting, H. (1979). Boundary-layer theory. 7th ed. Translated by J. Kestin.

New York: McGraw-Hill.

Shah, R. K. and London, A.L. (1978). "Laminar flow forced convection in ducts:

A source book for compact heat exchanger analytical data." In Advances in

Heat Transfer, ed. Thomas F. Irvine, Jr. and James P. Hartnett. New York:

Academic Press.

Speziale, Charles G. (1987a). "On nonlinear K-1 and K-e models of turbulence." J.

Fluid Mechanics, 178, 459-475.

Speziale, Charles G. (1987b). "On the advantages of the vorticity-velocity

formulation of the equations of fluid dynamics." J. Computational Physics, 73,

476-480.

Steger, J. L., Ying, S. X. and Schiff, L. B. (1986). "A partially flux-split algorithm

for numerical simulation of compressible inviscid and viscous flow." Proceedings
195

of the Workshopon Computational Fluid Dynamics. Institute of Nonlinear

Sciences.University of California, Davis, CA.

Tenpas,P. W. and Pletcher,R. H. (1987). "Solution of the Navier-Stokes equations

for subsonic flows using a coupled space-marching method." AIAA Paper

87-1173.

Van Driest, E. R. (1952). "Investigation of laminar boundary layer in compressible

fluids using the Crocco method." NACA TN 2597.

Van Dyke, M. (1970). "Entry flow in a channel." J. Fluid Mechanics, 44, 813-823.

Wang, Y. L. and Longwell, P. A. (1964). "Laminar flow in the inlet section of

parallel plates." AIChE J., 10,323-329.

White, F. M. (1974). Viscous fluid flow. New York: McGraw-Hill.

Wong, A. K. and Reizes, J. A. (1984). "An effective vorticity-vector potential

formulation for the numerical solution of three-dimensional duct flow

problems." J. Computational Physics, 55, 98-114.

Wong, A. K. and Reizes, J. A. (1986). "The vector potential in the numerical

solution of three-dimensional fluid dynamics problems in multiply connected

regions." J. Computational Physics, 62, 124-142.

Yang, H. and Camarero, R. (1986). "An improved vorticity-potential method for

three-dimensional duct flow problems." Int. J. for Numer. Meth. in Fluids, 6,

35-45.
197

8. APPENDIX A: THREE-DIMENSIONAL VORTICITY AND

DILATATION TRANSPORT EQUATIONS

The momentum equations for a three-dimensional flow of a Newtonian Stokesian

fluid will be presented here. They will be written in terms of vorticity (actually _) and

dilatation as the dependent variables. First a point will be made about the number

of transport equations which must be solved for the vorticity in three dimensions.

In three dimensions it is only necessary to solve two of the three component

vorticity transport equations because the third vorticity can always be obtained by a

linear combination of derivatives of the other two. To see this, consider the vorticity

definition,

where,

---'> -..+

0,=V× V

0,1 = Wy - Vz

0,2 = UZ -- WX

0,3 = vx - Uy

therefore,

0,3z = -0,1z - 0,2y

_¢___,?_it_ PACE. BL._K NOT FILMED
198

Hence, the third vorticity component can always be obtained from the other two

components. In the above, the numbers 1,2,3 represent directions and the letters

x, y, z are derivatives.

The vorticity transport and the dilatation transport equations follow. They

represent the momentum equations for a three-dimensional flow. To compute a flow

field using the dual potential method, additional equations are required to satisfy

mass and energy conservation. Also, the Poisson equations for the potentials must

be solved (Equations 2.13 and 2.14):

V2¢ = B

V 2--+A -- -- O,}

Equations 8.1-8.3 display the vorticity transport equations in the dependent

variable, (_). The three-dimensional dilatation transport equation is presented in

Equation 8.4. The ideal gas law (p = pRT) has been used to replace the pressure

gradient term in the momentum equations. Body forces are neglected. The following

equations have not been coded. Their numerical behavior is unknown.
199

Vorticity transport: x-component

uz + -_ Ree #zPz + #ypy + #zPz

ll[p Wl Wl ]

+.(?)+uz(?)

R

11

10+

(_)_ + + <7) 0 Pz 20-+ -_

+ g_zv • v- gNv. v

tte +(7)y(_. +_,z)

+(-;)_. #z (2Wz 2-_,_,._)

-9-)z++V=)

-("_>o,z (_+=_)] (8.1)
200

Vorticity transport: y-component

(?)_+_(?)_+v(?)+w(?)z=

[vy + p211Re {#V2P + #xpx + pypy + #zPz

+- #x pRee •

+- #V 2

P

R [Txpz - Tzpz]

+p---g

1 1

(?)

+ _y(7)y o'2 +#z(7)z _v2 ]

10_

+ gv.v-(7-)sgv.v

#z (2uz 2 -_

(8.2)
201

Vorticity transport: z-component

[WZ+p21Rel {pV2P+pxPx+/_YPY-}-PzPz}] (_p)

+p-_

+_ _ v2, + _v.

+ ._v.v- 'Y o _v.v

Pz +(7)_(_z+_)

-sv._)

(s.3)
202

Dilatation transport:

B t+uBz+vBy+wBz- (4.) _-_ V2B=

-(u 2 + v 2 + w 2 + 2uyv_ + 2UzWx + 2VzWy)

I_ 4 B

4B (w3-Wlz+_ Y)

_ 4

_RV2T RTv2 p_p R (T*pz+Typy+TzPz)+--_ RT. "(pz+p2+p2) 2

pRe p

2

+7_ ("_u_ +"yy'' +"=wz) 2 Z V2" (8.4) 3 pRe
203

9. APPENDIX B: ALTERNATE NON-PRIMITIVE VARIABLE

EQUATIONS

An approach different from that presented in the main body of this report is to

form a different combination with the z, y momentum equations.

In two dimensions the governing equations are:

1. continuity

2. x momentum

3. y momentum

4. energy

The cross product applied to the momentum equations yields the vorticity transport

equation. This leaves one with still another usage of the momentum equations. In

the main body of this report the dilatation variable, B = uz + Vy, was selected as the

final usage. Another possible operation on the momentum equations is to compute

0

(x momentum) + _z (Y momentum). Then, choosing the dependent variable to

be (Uy +vx), an equation for the rate of shear deformation is derived. Let I" = Uy +vz.
204

The two-dimensional equation for F is

p [F t + uFx + vFy + FB] ttV2F _---_VF V# R_V2# = -2pxy + _ + • +

1 0 2 Du Dv

+ Re OxOy #B - PY-_ - Px Dt

(9.1)

Then, the momentum equations provide the solution for w and F. Recalling the

definitions,

¢,, = vx - uy (9.2)

F = vx +uy (9.3)

one obtains the following compact formulas for some velocity derivatives:

w+F

vx - (9.4)

2

F-w

- 2 (9.5)

This adaptation to the dual potential method could be used as an inverse solution

procedure since the wall shear stress is simply #F.

NOTICE: At a no-slip impermeable boundary, vx = 0. There is no trouble

getting the vorticity at the wall. Therefore, F should be simple to obtain once the

vorticity is computed! At the wall F = -w. Also, the potentials are not needed for

this use of the momentum equations.