Citation
Numerical simulation of transonic turbulent flows over a complex configuration

Material Information

Title:
Numerical simulation of transonic turbulent flows over a complex configuration
Creator:
Yang, Shu-Cheng, 1954-
Publication Date:
Language:
English
Physical Description:
vi, 144 leaves : ill. ; 28 cm.

Subjects

Subjects / Keywords:
Aircraft wings ( jstor )
Boundary conditions ( jstor )
Grid generation ( jstor )
Jacobians ( jstor )
Navier Stokes equation ( jstor )
Simulations ( jstor )
Solid surfaces ( jstor )
Transonic flow ( jstor )
Turbulent flow ( jstor )
TVD schemes ( jstor )
Aerodynamics, Transonic ( lcsh )
Aerospace Engineering, Mechanics and Engineering Science thesis Ph. D
Dissertations, Academic -- Aerospace Engineering, Mechanics and Engineering Science -- UF
Turbulence -- Mathematical models ( lcsh )
Genre:
bibliography ( marcgt )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1990.
Bibliography:
Includes bibliographical references (leaves 139-143).
General Note:
Typescript.
General Note:
Vita.
Statement of Responsibility:
by Shu-Cheng Yang.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. §107) for non-profit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact the RDS coordinator (ufdissertations@uflib.ufl.edu) with any additional information they can provide.
Resource Identifier:
026233036 ( ALEPH )
25034739 ( OCLC )

Downloads

This item has the following downloads:


Full Text









NUMERICAL SIMULATION
OF TRANSONIC TURBULENT FLOWS OVER
A COMPLEX CONFIGURATION














By

SHU-CHENG YANG


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


1990













ACKNOWLEDGEMENTS


The author wishes to express his sincere appreciation to his dissertation advisor, Dr. Chen-Chi Hsu, for his guidance and encouragement throughout this endeavor. Special thanks are expressed to Dr. W. Shyy for his invaluable criticism of the draft of this dissertation and his stimulating discussions. Many thanks are also extended to the other supervisory committee members, Dr. B.M. Leadon, Dr. U.H. Kurzweg and Dr. A.K. Varma for their interest and comments on the research project.

Thanks are also due to Dr. M.H. Chen and Dr. N.H. Shiau for their assistance on their previous work of the TVD scheme and computer codes. Special thanks must go to Mr. J. Ordonez for his tireless proofreading of the manuscript. Also, the author has a debt of gratitude to Dr. N.C. Yao and Dr. C.G. Tu of the Chinese Navy, without whose encouragement and help he would not have had a chance of pursuing the advanced study here at the University of Florida. Finally, the author is indebted to his family for their continuous support.

This work was partially supported by a NASA-Ames research grant, NAG2473. Most of the calculations were performed with a Cray-2 of the NAS System at NASA Ames Research Center.


ii














TABLE OF CONTENTS



Page

ACKNOWLEDGEMENTS,-.................-.--........-- .-.--....... ii

ABSTRACT ... - ---................................................. v

CHAPTERS

I INTRODUCTION ----......................................... 1

II GOVERNING EQUATIONS ................................... 7

2.1 Navier-Stokes Equations and Equation of State ................ 7
2.2 Non-Dimensionalization -..-................................ 12
2.3 Transformed Navier-Stokes Equations ........................ 14
2.4 Thin-Layer Approximation.--.............................. 17
2.5 Turbulence Model ........................................ 19

III NUMERICAL ALGORITHM - -............................... 25

3.1 Beam and Warming Scheme .--............................. 26
3.2 TVD Scheme .---........................................ 32
3.3 Flow Solver .-- - - --........................................ 43

IV THE FLOW PROBLEM AND SOLUTION STRATEGY .......... 45 4.1 The Flow Problem and Preliminary Work ................... 45
4.2 Solution Strategy .---..................................... 51
4.3 Six-Block Grid Topology - --............................... 54

V GRID GENERATION ----..................................... 62

5.1 Surface Grid Generation ..--............................... 63
5.2 Volume Grid Generation .................................. 66
5.3 The Six-Block Grid ....................................... 74




iii










VI SOLUTION METHOD .-.-..-........................-........ 88

6.1 Jacobian and Metric Coefficient ............................. 89
6.2 Boundary Condition .---................................... 95
6.3 Interface Boundary Condition -- --........................... 105
6.4 Convergence Criterion - --.-................................ 108
6.5 Pressure and Shear Stress Coefficient . -..................... 109

VII RESULTS AND DISCUSSION -.-............................. 111

7.1 The Result from Base Grid -............................ 112
7.2 The Result from Refined Grid.-.......................... 115

VIII CONCLUDING REMARKS -- --............................... 130

APPENDICES

A JACOBIAN MATRICES .................................... 132
B EIGEN-FUNCTIONS ..--..................................... 134
C CLUSTERING FUNCTIONS ---............................... 136

REFERENCES ................................................ 139

BIOGRAPHICAL SKETCH ...................................... 144


iv













Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy NUMERICAL SIMULATION OF TRANSONIC TURBULENT FLOWS OVER A COMPLEX CONFIGURATION By

Shu-Cheng Yang

December 1990


Chairperson: Chen-Chi Hsu
Major Department: Aerospace Engineering, Mechanics and Engineering Science

A multi-block computational method is developed for three-dimensional high speed turbulent flows over complex configurations. In this method, the flow domain is divided into several contiguous blocks such that each block is partially bounded by a solid surface. The chosen grid topology has advantages in computing eddy viscosity as well as in applying a thin-layer Navier-Stokes code. In the solution process, each block is then treated as an independent flow problem governed by the same set of time dependent thin-layer Navier-Stokes equations. The interface boundary conditions between blocks are updated at every time step of the solution algorithm. The application of the method with distance weighted interpolation for updating the interface boundary condition is rather simple and straightforward; however, special measures in grid topology and


v








generation are required. The flow solver employed is a total variation diminishing thin-layer Navier-Stokes code implemented with an algebraic eddy viscosity model of Baldwin-Lomax. In this study, the proposed method is investigated for the simulation of a transonic turbulent flow past a wing-fuselage configuration in which the flow field is properly divided into six blocks. A coarser grid was first used for flow simulation, followed by grid refinement studies. With the use of the refined grid, the calculated wing surface pressure distribution agrees well with measured data. Numerical results obtained also show that the computed shock location and pressure distribution on the wing can be in good agreement with the experiment data if the block grids used are adaptive to important solution characteristics. It is concluded that the proposed method is indeed a very promising approach to be developed further for complex configuration flow simulation.


vi













CHAPTER I
INTRODUCTION


For a high speed flight vehicle, numerous complex three-dimensional fluid flow phenomena appear, including the formation and shedding of vortices, the shock-wave boundary-layer interaction and induced flow separation, and the shock on shock interaction. Since the performance of a flight vehicle is strongly affected by these complex flow phenomena, accurate prediction of them and associated aerodynamic forces and moments are critical for advanced design of a modem flight vehicle. At present an accurate wind-tunnel simulation of the complex flow problem in flight conditions can be very difficult, if not impossible. With the advent of supercomputers and the advancement of solution techniques for nonlinear problems, computational fluid dynamics (CFD) is under extensive development and has been applied to complement the wind-tunnel experiment and field test in the design process of a flight vehicle.

The use of CFD in the aerodynamic design of flight vehicles has

progressed rapidly over the last few years. In the early days, CFD was used to support the validity of a design that developed in the wind-tunnel by trial and error. The state-of-the-art of the CFD method has progressed to the point of being regarded as an important design tool for many flow problems [1]. The reduction of wind-tunnel occupancy is only a direct benefit of CFD. Detailed


1








2

flow information that is unattainable from the wind-tunnel test often can be obtained by numerical simulations. This information about the underlying flows can lead to a better understanding and ultimately to a better design. In fact, the implementation of CFD has fostered a revolution in the design process of flight vehicles.

Current CFD methods have only demonstrated an ability to simulate flows about complex geometries with simple physics or about simple geometries with more complex physics. The panel methods are unique in that they provide a capability for solving the flow about completely general configuration. Their principal limitation is that they are restricted to simple physics as modeled by the linear equation. The Euler code provides more flow information than the panel method; however, the important viscous effects are excluded in the approach. For more realistic flow simulations, the application of Navier-Stokes equations is required. However, the use of full Navier-Stokes equations is simply too expensive to be practical for today's computers. An alternative is to use the thinlayer Navier-Stokes equations. In general, the governing equations are approximated by a numerical scheme and then solved in either a structured or an unstructured grid network. The latter is quite flexible in gridding but it requires more work in bookkeeping and computer coding. In the structured grid approach, the governing equations are first approximated by either a finitedifference or a finite-volume scheme, and then solved in a boundary-fitted curvilinear coordinate system. This approach has advantages in boundary condition treatment and computer coding; however, it may have difficulty in








3

generating good grid system for complex configurations. As evidenced by papers published in journals and presented at conferences, e.g. [2-4], the flow simulations involving simple geometries have progressed considerably; however, the capability for turbulent flow problems involving complex 3-D configurations is still in a development state. Hence, further research and development effort is needed to advance the CFD methods for more complex geometry flow simulations.

For flow past a complex configuration, such as a flight vehicle, the

generation of an acceptable single structured grid for the entire flow field is rather difficult. It is much easier to divide the flow field into several subregions, termed blocks, and then generate a good grid for each subregion. Hence, it is natural to develop a solution method based on the multi-block grid network. There have been a number of multi-block solution techniques proposed for the Euler simulation of aircraft models and other complex configurations [5-7]. However, the topology of block gridding for viscous turbulent flow simulation is more involved than that for inviscid flows and requires a lot more grid points to resolve the thin viscous layer. In general, a good block grid topology for Euler simulation may not be a proper one for Navier-Stokes simulation.

There are also zonal methods developed to use different equation sets in different regions of the flow field. In general, the flow field is divided into inviscid and viscous zones, and then the Euler equations are applied in the inviscid zone, while the boundary-layer or Navier-Stokes equations are used in the viscous zone. A zonal approach that was based on the Euler/thin-layer NavierStokes equations was proposed and tested successfully for transonic wing flow [8].








4

Another zonal algorithm that applied the Euler and thin-layer Navier-Stokes equations for simulating the flow field of isolated wings was reported in [9]. The flow field was divided into four zones; two inviscid zones with coarse grids and two viscous zones with clustered grids. The computed results were in good agreement with experimental data for the surface pressure, except in the immediate vicinity of the tip and in the shock-induced separated region. The zonal approach also has been used to simulate flow problems for more complex geometry. A transonic viscous flow past an F-16A wing-fuselage configuration has been simulated by a zonal approach [10]. The flow field was divided into as many as 16 zones; in the inner zones adjacent to no-slip surfaces the thin-layer Navier-Stokes equations are solved, while in the outer zones the Euler equations are used. The prediction of wing surface pressures was quite good, except at the leading edge and the aft-shock position. The zonal approach has been shown to be rather efficient. However, the difficulties are to locate properly the interfaces for the inviscid and viscous zones and to patch the solutions that are obtained from different equation sets at these interfaces. Recently, a zonal method with Chimera overset grid scheme was investigated for subsonic turbulent flow about the F-18 fuselage forebody and the combined wing-fuselage [11]. Special coding was required for the overlapped grid region, but the computed results have been shown to be in good agreement with flight-test flow visualization and surface pressure measurements.

In this study, a novel multi-block computational method is proposed and explored for the simulation of transonic turbulent flows over complex








5

configurations. In this method, the flow domain is divided into several contiguous blocks such that each block is partially bounded by a solid surface. For ease in applying thin-layer approximation and the Baldwin-Lomax turbulence model, the solid surface is mapped onto an entire boundary plane in computational space. In the solution process, each block is then treated as an independent flow problem governed by the same set of time dependent thin-layer Navier-Stokes equations. The interface boundary conditions between blocks are updated at every time step of the solution algorithm. The solution algorithm with distance weighted interpolation for updating the interface boundary condition and the computer coding are rather simple and straightforward. The multi-block method offers several distinct advantages over the single block computation. In fact, the use of block grids makes the problem of simulating flow fields about complex geometry more tenable. If a geometric feature is changed, only the related block requires to be modified without completely redoing the basic grid topology. Also, the solution procedures are adapted to the modern computer architecture; in parallel processing, each block can be solved at different processors; or in case short of main memory, each block can be solved at a time while the other blocks remain on extended storage.

A transonic turbulent flow over a wing-fuselage configuration is

investigated for the development of the proposed method. The flow field is properly divided into six contiguous blocks. A coarse base grid was first generated and used for the flow simulation, and then a refined grid was generated from the base grid by halving the grid spacing along the wing in the








6

chordwise direction. A flow simulation package that includes grid generation, flow solver, and plotting program, was developed and applied to the six-block flow problem. The method was first investigated for a transonic turbulent flow of Mach 0.8 past the wing-fuselage configuration at zero angle of attack, then simulations for different angles of attack, a=4* and a=-3', were conducted to gain further insight and understanding on the solution algorithm. With the use of the refined grid, the calculated wing surface pressure distribution agrees well with measured data. Numerical experiments conducted also show that special measures and care are required in generating good grid systems involving excessive distortion between the physical and computational domains; however, the results obtained indicated that the proposed method is a very promising approach for the simulation of complex configuration aerodynamics.













CHAPTER II
GOVERNING EQUATIONS


Some theoretical background that required for transonic turbulent flow

simulations are described in this chapter. The Navier-Stokes equations and ideal gas equation of state are first discussed. Then, they are non-dimensionlized and transformed to a curvilinear coordinate system for ease in numerical applications. The thin-layer Navier-Stokes equations valid for high Reynolds number flows are obtained by neglecting the diffusion terms along the direction of the solid surface. The time-averaged thin-layer Navier-stokes equations are used for turbulent flow simulations and the Baldwin-Lomax model [12] is used as the turbulence closure.


2.1 Navier-Stokes Equations and Equation of State


2.1.1 Navier-Stokes Equations


The fundamental equations of fluid dynamics are based on the following universal laws of conservation: (a) conservation of mass; (b) conservation of momentum; (c) conservation of energy. The Navier-Stokes equations that govern most of the compressible flow problems can be obtained by applying those universal laws to a fluid flow [13-14].


7








8

Continuity equation. The conservation of mass law applied to a fluid

passing through an infinitesimal, fixed control volume without sources or sinks of mass yields

8p + apu + 8pV + =pW 0 (2.1)
at ax ay az

where p is the fluid density and u, v, w are the components of fluid velocity in x, y, z direction, respectively.

Momentum equations. The conservation of momentum law applied to the control volume while neglecting body forces reads

apu+ (PU2+p)+ apuv apuw _ a. ar a r
at ax ay az ax ay az

apv + apvu + a( apvw ar ar ar
- y - + -" y+ - I-z (2.2)
at ax ay az ax ay az

apw apwu + apwv + a _ ar, +r ar.
at ax ay az ax ay az

where the components of the viscous stress tensor r, for a Newtonian fluid are given by the constitutive equation

au a w= 2 au av aw
Tx1C = ax ay a ax 3
A(-+ -- + - )+2 - (2 av a)
ax az ax 3 ax ay az

av aw au av 2 a w a
r=A(--+--+--)+2,u - =-(- ---ay az ax ay 3 ay az ax

aw au av aw = 2 aw au av
r= X( -+ -+--)+2 - -p( - --) (2.2a)
az ax ay az 3 az ax ay

au av
r= 1( -+ - ) = r
S ay ax








9


av aw 8z ay aw au
rz = p( + - ) = *r,
ax az

Here, p is the thermodynamic pressure and y is the dynamic viscosity. The fluid is assumed to be isotropic and satisfying the Stokes condition, A = -2/3 p, which is a good approximation for flow problems without relaxation process.

For air, the viscosity p depends mainly upon the temperature and can be estimated by using an interpolation formula based on Sutherland's theory of viscosity, i.e.,

T m T. + S
= T + S (2.2b)

where p. denotes the viscosity at the reference temperature T., and S is a constant which assumes the value 198.6'R [15].

Energy equation. The conservation of energy law applied to the control volume and assuming no external heat addition results aE a a a = #3 ~ 23
+ [(E+p)u]+ [(E +p)v] + [(E+p)w] - Ip" +g -L- (2.3)
at ax ay az ax ay az

where

aT
, = Urx +vr +wr +k -ax

aT
+ = uryx+vryy+wry +k ay (2.3a)
ay

aT
16Z = urz +vrT,+wr. +k -ay








10

Here, E is the total energy per unit volume, T is the temperature and k is the thermal conductivity. The heat loss by conduction qj is related to the temperature gradient by Fourier's law of conduction

aT
-k -T (2.4)
axj

The equations described above, Eqs.(2.1-2.3), constitute a complete set of the time dependent Navier-Stokes equations that is valid for a class of compressible flow problems without heat addition and chemical reaction (or relaxation). The Navier-Stokes equations includes one continuity equation, three momentum equations and one energy equation, but there are seven unknowns, namely, p, u, v, w, p, E, and T. The relationships between the thermodynamic variables (p, p, T, E) are established by the equation of state of the working fluid.


2.1.2 Equation of State


In this study, air is the assumed working fluid and the perfect gas equation of state is employed. For a perfect gas, the thermal and caloric equations of state are

p = pRT, e = CT ( C, = constant) (2.5)

respectively. Here, R is the universal gas constant, e is the internal energy. Some other useful relations are

h = CT, -y = C,/C,, a2 = yRT (2.6)

where h is the enthalpy, C, and C, are the specific heat at constant pressure and volume, respectively, y is the ratio of specific heats, and a is the local speed of sound.








11

If only internal energy e and kinetic energy are considered significant, the total energy E can be written as

E = p[e+(u2+v2+w2)/2]

Consequently, one obtains two specific equations relating p, T to the dependent variables p, u, v, w, and E, i.e.,

p = (-y-1)[E-g 1P(U2+V2+w2)]
p 2 [ (2.7)
T = _[ Ef..12+V2+wlv1
R P 2

If expresses heat conduction q in terms of internal energy through the use of caloric equation of state, one can rewrite Eq.(2.4) as follows.

aT k ae yp ae
q=-k -- - (2.8)
ax C, axj Pr axj

Here, the Prandtl number Pr is defined as Pr = Cp/k, for air Pr = 0.72, approximately, and the internal energy is given by

E 12
e = P- y(u +v2+w2) (2.8a)

and the p's in Eq.(2.3a) can also be rewritten as


fx = urx +vr, +wr, + y ae Pr ax

fi = ur +vr+wr, + / (2.9)
'~Pr ay

fi = ur +vr3,+wr, + -Y ae Pr az








12


2.1.3 Vector Form


The Navier-Stokes equations, Eqs.(2.1-2.3), form a system of five equations and can be cast in vector form as

aq aE aF 8G aE, aF, aG,
-_+ -+ -+ -- - + --+ - (2.10)
at ax ay az ax ay az

where q is an unknown vector, E, F, G are the flux vectors and E,, F,, G, are the viscous flux vectors. The explicit expressions for these vectors are


pu P p2 +p pvU
q pv, E= puv F pv2+p

,E (E+p)u J (E+p)v
(2. 10a)
pw ' O ~~ ''
pwu r
G = pwv - E,= rX, F, = rYY- G, = y
G1pwv+ rx Ez z
(E + p~w. lax 1,. A

in which ru are given by Eq.(2.2a), p's are defined in Eq.(2.3a), and p, T are related to q by Eq.(2.7).


2.2 Non-dimensionalization


To reduce the number of parameters that characterize the same class of problems, a dimensional analysis is employed. The non-dimensional quantities are obtained by choosing characteristic scales in the following manner:

. x . y * = * t
L ' ' L ' L/a,

. U . v . w T
U= -, V - , w T = T - (2.11)
a. a. a.T.








13

E . p . p
E= P= -, p = 2

where the non-dimensional variables are denoted by an asterisk *, the freestream conditions by -, a, is the freestream speed of sound, and L is a characteristic length scale.

By substituting the above relations into the governing equations, Eq.(2.10), the non-dimensional equations in vector form are

aq aE aF* aG* aE aF'* aG
-- + --. + - + -r = Re-'( -+ -+ -i) (2.12)
at ax ay az ax ay az

where the Reynolds number, Re, is defined as

Re =

and the relationships from the perfect gas equation of state, Eq.(2.7) become

1 *
p (--1)[E*-,P(u +v"+w)= (_y-1)p e'


E* 2
T = 7(7-1)[.f* 2.(u2+v2+w")] (2.13)

The unknown vector q* and flux vectors E*, F*, G*, E,*, F,' and G, have the same form as in Eq.(2.10a), except that all the quantities are now non-dimensional.

In summary, the dimensionless equations govern a class of flow problem. Only two parameters, namely, the Reynolds number, Re, and the Prandtl number, Pr, appear in the equations, the Mach number is dropped out due to the use of a, as characteristic velocity scale. In the following, the asterisks * are removed from the non-dimensional equations except where noted.








14

2.3 Transformed Navier-Stokes Equations In the finite-difference approach, to enhance numerical accuracy and efficiency in applying boundary conditions, the governing equations are transformed from the Cartesian coordinates to a boundary-fitted curvilinear coordinate system. Consequently, a general computational code can be developed such that the computation is performed in an uniformly-discretized computational domain [16].

If one introduces a generalized coordinate system

= ( xy,z,t )

17 = ( x,y,z,t ) (2.14)

= ( x,y,z,t )

r =t

and applies the chain rule of partial differentiation, the non-dimensional governing equation, Eq.(2.12), becomes q,+ q,+ q,,, + q,+ (,Ec + ,E,7+ .Ec) + ( F+ +yF,, ++F) + (Gc+,G,,+ G ) = Re"[( .E,+,E,,+ E,,)+( yF,,+ +yF,,+ F,,)+( G, ++,G,,+,G )] (2.15) The metric coefficients appearing in these equations can be determined by the matrix relation as follow


[C '7Y ?z7t yf y,, yC y' (2.16)
x~~ z z z,
0' y7z7 = 1 dy y 1y1y'


or








15


[( 1 x x,, x C
17, 1, = J, y y
. f y . I Iz f z,, z C
[y,,ze-yrz,, -(x,,ze-x.z,,) x,,y,-x~y, -( z-yze) xzz, -x-z -(x"yC-x-yX) yez,,yze -(x z,,-x,,7z ) x~y,,-x,,y
ctl (xt (
17 ~ f 1 y7 , 17 7

where J is Jacobian of the transformation and J1, representing the volume of the grid cell in physical space, is given by

1= a(x,y,z) x x, , x
a(Q,1S) I zc z, z'
= x +yz-y z,)-x,(yz -yrzd)+ x(yez,-y,,ze) (2.17)

In order to put the transformed equations, Eq.(2.15), into conservation-law form, the equation is first divided by the Jacobian and then rearranged to give

(q,)/J + (q + Ee + Fc+ zGe)/J + ( +,+q,, ++,,+ 7F,+ zG,)/J

+ (q&+ .Ec+ YFc+ zG )/J = Re-[((E,,+ X,+ezG,,)

+ ( + +E,+tF,+%G,) + ( .E,+ YF,+ zG,)]/J

and then rewritten into the following conservation law form

(q/J), + ((q+ .E+ YF+ zG)/J)c + ((,7q+t E+yF+nzG)/J),

+ (( q+ .E+ YF+ zG)/J)c = Re-'[((QE,+ YF,+ zG,)/J)c +

+ ((7E,+ ,7F,+ 9,G,)/J),, + ((,E,+ F,+ zG,)/J)j] (2.18)

in which the following metric identities [10] have been applied.

()+ ( )+ ( ) = 0



( )+ ( )+0
a ~~ ~ 0 , J a








16


a


J


(-1 )+ an J


(2.19)


a ( ) = 0 a J


(L)+ a(-)+ a ( I.) = 0
a8 J an J a J

The transformed equation, Eq.(2.18), has the same form as the original equation, Eq.(2.12), and can be written as


AA A A A
aq +E aF aG . E,
-+-+ - + - = Re(--+ 8r a8 an a a


A
aF, +
- + an


aG,
--)


(2.20)


where


A
q = q/J,


F = (ji~q+q.E+ryF+rG)/J,
A
E, (c.E,+ CYF,+ .G,)/J,


A
E = (,q+ ,E+ YF+ ,G)/J,
A
G = ( tq+ ,E+ ,F+ ,G)/J,
A
F, = (q.E,+P,~.F,+n.iG,)/J,


G, = ( .E,+ gF,+ G,)/J

The explicit expressions for the transformed unknown vector, flux vectors and viscous flux vectors are


A pu
q =J-' pv ,
pw I
E

PV )
PuV + 17,P
F = J-P VV+17p -,

(E+p)V-ri7p .


0
X X+ y X+ rzx C.Pr + IYY z)%


pU

E = J pvuU+C p
pwU+ p
(E+p)U-QpJ

AP pUW +
= *1puW+gp
G pvW+ryp -y,
pwW+ p
(E+p)W- tp


I ,


A
F,


[00
=J-1 IlxTXX + 7yryx + 7zrzx
1 7x7rxy +17yryy + nz1zy ??xr z+ 7 yT + nzrT
,7'x +VV +770


(2.20a)


(2.20b)


I ,


A V
F.,








17


0
AX i %D.Y Yryy+ *zzi

=Xp J1+ UY ,,+ ZA

where U, V and W are contravariant velocity components given by

U = 1+c'u+q+ 'w

V = ,+17,u+iyv+,7w (2.20c)

W = r+ Xu+ Yv+ w

It is noted that the velocity gradients in r,, Eq.(2.2a), and the

temperature gradient in p's, Eq.(2.3a), have to be transformed into computational space also, e.g.,

'rx au +av +8w )+A u
r= -+ + -)+2pax ay az ax

au u au +av _av _v
a8 a a a Y an a

+ ( aw-,L + aw nz+ aw .)]+2u[( a L + a ,+ -au
at an a a an as


2.4 Thin-Layer Approximation


The numerical simulation of complete Navier-Stokes equations, Eq.(2.20), is in general rather time consuming and needs a huge amount of computer storage. For high Reynolds number flows, however, the viscous effects are often confined to a thin layer near no-slip boundaries, called the boundary-layer, outside of which the flow can be considered to be mostly inviscid. With this in mind, it is reasonable to assume that the diffusion processes along the body








18

surface are negligibly small in comparison with those in the direction normal to the body surface. As a consequence of this assumption, all viscous terms containing derivatives parallel to the body surface are dropped since they are substantially smaller than viscous terms containing derivatives normal to the surface, while all other terms in the momentum equations are retained. One of the advantages of retaining the terms which are normally neglected in boundarylayer theory is that separated and reverse flow regions can also be computed. Also, flows which contain a large normal pressure gradient can be readily computed.

Another justification of the use of a thin-layer approximation can be given from a computational viewpoint. In practice, a substantial fraction of the available computer storage and time is expended in resolving the cross-stream gradients in the boundary layer since a highly stretched grid is required, and the resolution along the body is treated as what is needed in inviscid flow. As a result, even though the full derivatives are retained in the equations, the gradients along the body are not resolved in an adequate manner. Hence, one can drop those terms that are not being adequately resolved, provided that they are reasonably small.

Since the thin-layer model requires a boundary-layer type of coordinate system, it is assumed that the solid body surface is mapped onto an entire g1plane in the transformed space, then the thin-layer approximation requires that the and t derivatives in the viscous terms be dropped. Upon simplifying the complete Navier-Stokes equations, Eq.(2.20), the system of governing equations becomes








19



aq a E aF aG .1 S
+ - + -+- = Re( --) (2.21)
ar a a7 a a

where

0
J ~ m~u,+ m27.
S =J miv + m2y- (2.21a)
m1w. +W.2~

I = A( .2+ ,2+ .2)

m2 = u(,u, + Yv,+ zw )/3 (2.21b)

m3 (u2+v2+w2),/2 + Pr-1(-11)-(T),

Note that the notation ^ has been removed from the above equations for simplicity and the heat conduction term are replaced by

q _ y ae _ aT (2.22)
Pr a (-y-1)Pr a8

The thin-layer Navier-Stokes equations, Eq.(2.21), in conjunction with the equation of state, Eq.(2.13), constitute the required governing equations for this study. They are valid for both laminar and turbulent flow problems. However, due to the widely disparate scales of turbulence, it is impractical to apply the equations directly for turbulent flow simulation. An alternative is to use the time-averaged equations coupled with a turbulence model.


2.5 Turbulence Model


The time-averaged equation is obtained by decomposing the dependent

variables into time mean and fluctuating components and then time-averaging the








20

entire equation. In the averaging process, new unknowns will be introduced and must be treated by some approximating techniques so as to close the problem. In this study, the mass-weighted time-averaged equations are used for turbulence flow simulation with the eddy viscosity that appears in the averaged equations provided by the Baldwin-Lomax two-layer algebraic model. In this approch, the mass-weighted time-averaged thin-layer Navier-Stokes equations will have the same form as the original equations, Eq.(2.21), while all the dependent variables are time-averaged. The effective coefficient of viscosity I is comprised of two parts, the laminar viscosity p, and the turbulent eddy viscosity p,. The laminar viscosity p, is obtained from Sutherland's theory of viscosity, Eq.(2.2b), and the eddy viscosity p, is provided by a turbulence model [16].

Similarly, the effective thermal conductivity k is also comprised of two parts, k, and k,. The heat conduction is expressed in dimensional quantities as

aT C, & a
q = -(k,+k,) - _( C,+ ) (2.23)
az Pr, Pr, az

The heat conduction term in Eq.(2.22) is then replaced by its dimensionless counterpart

q - ( A' + At ) aT (2.23a)
q - 1 Pr, Pr, ag

Here, the turbulent Prandtl number Pr, is defined as Pr, = pcp/k, and is about

0.9 for air [10].

In this approach, the eddy viscosity pt is the only quantity that has to be modeled. In what follows, the theoretical background about the algebraic model will be discussed first, and then the Baldwin Lomax model is described.








21


2.5.1 Two layer algebraic model


In analogy with the coefficient of viscosity in Stoke's law for laminar flow, Boussinesq [17] introduced a mixing coefficient, pt, or referred to as eddyviscosity, for the Reynolds-stress that appears in the time-averaged equation.

au
t= a (2.24)
ay

where u is the streamwise velocity, y is the spatial coordinate normal to the solid surface, and au/ay is assumed to be the only significant component of strain rate.

With Prandtl's mixing length hypothesis, the coefficient pi, is then related to mean velocity field by

p = P12 -u (2.25)
ay
where 1 is the mixing length.

The algebraic model is then to establish a relation between the mixing length and the characteristic length of the respective flow. For simple wall shear flows, the mixing length in the law of the wall region is proportional to the normal distance from the wall, i.e., 1 = ky, where k = 0.41 is a universal constant determined experimentally. In the viscous sublayer, however, the mixing length theory is not valid. A corrected form of the mixing length valid down to the wall was deduced by van Driest. He introduced a damping factor and wrote

1 = ky[1-exp(-y+/A+)] (2.26)








22

where y' = ypu,/ is the law of the wall coordinate, A' denotes an empirical constant which equals 26 for flat-plate flow, and u, the friction velocity, is defined as u, = (r./p)', where r, is the shear stress on the wail.

In the law of the wake region, the Clauser model is commonly used in conjunction with a factor to account for the intermittent effect. The eddy viscosity is then written as

y, = kCUs*Y (2.27)

where the value k, is obtained experimentally and assumed to be 0.0168, U is the streamwise velocity on the edge of the boundary layer, 6* is the boundary layer thickness, -y is the intermittency factor deduced from a curve fit of measured data and given by Reynolds and Cebeci [18] as

y = [1+5.5(y/6)6]-l (2.27a)

where 6 is the boundary layer thickness.

Based on the above arguments, Cebeci and Smith [19] proposed a well known two-layer algebraic model as follows. In the inner-layer:

t= p12 au (2.28)
ay

1 = ky[1-exp(-y+/A*)] (2.28a)

In the outer-layer:

p= kCU6*-y (2.28b)

-y = [1+5.5(y/6)6]1 (2.28c)








23

This model has been shown to be in excellent agreement with experimental observations for simple shear flows. However, the model requires the boundary layer thickness and the edge velocity which constitute a practical disadvantage in the implementation of the model. In an attempt to overcome this difficulty, Baldwin and Lomax [12] modified the outer layer formulation to eliminate the necessity of finding the boundary layer thickness. In addition, they replaced the velocity gradient by the vorticity, which is invariant under coordinate transformation, so that the model can be easily applied to 3-D problems. The Baldwin-Lomax model is described next. 2.5.2 Baldwin-Lomax Turbulence model


The Baldwin-Lomax model is similar to the Cebeci-Smith two-layer

algebraic turbulence model. The eddy viscosity in a profile is calculated for each inner and outer layer, and then matching both layer as follows [12]:

t I(Ji..,, , y5 , (2.29)
("I)outer , y z y,
where y. is the smallest value of y at which values the inner and outer formulas are equal.

The inner layer follows the Prandtl-Van Driest formulation

('t)ir.', = pl21wI with 1 = ky[1-exp(-y*/A')] (2.30)

where y is the normal distance to the surface and Iw is the magnitude of the vorticity

Jw| = [(u,-v.)2+(v.-w,)2+(w.-u)2 ]' (2.31)








24

y+ = y(pr.)'/p, is the law of the wall coordinate, A+ = 26, and p., r., A, are the local density, shear stress, and laminar viscosity evaluated at wall, respectively.

The eddy viscosity for the outer layer is given by

(pt),.t,, = pk CF.IkFkib(Y) (2.32)

where k. = 0.0168 is the Clauser constant and C, = 1.6. The parameter F., is given by

F~k = min { F"m / where Ck = 0.25 (2.32a)
CAgy.U'djf/F.
and

Udi = (u2+v2+W2)1,. - (U2+V2+W2)",,i. (2.32b)
The quantities y. and F, are determined from the function

F(y) = ylwl[1-exp(-y+/A')] (2.32c)

where F. is the maximum value of F(y) and y. is the value of y at which F. occurs. The Klebanoff intermittency factor Felb(y) is given by

Fkleb(y) = [1+5.5(Cjy/yM)']4 and Ckeb = 0.3 (2.32d)

The Baldwin-Lomax model has been investigated and applied to a variety of 2-D and 3-D flow calculations with satisfactory results [10]. The model is capable of handling the compressible flow problems with modest flow separation [20]. The implementation of algebraic models requires only a small amount of computer time and storage; this is particularly important in 3-D computations.













CHAPTER III
NUMERICAL ALGORITHM


In numerical simulations, the transformed governing equations described in the previous chapter are approximated by finite-difference equations and then solved in a equally spaced computational space. The basic computer code used in this study is a TVD thin-layer Navier-Stokes code which is based on an original version of thin-layer Navier-Stokes code, ARC3D of NASA Ames Research Center. The original ARC3D was based on the Beam and Warming scheme [21], and has been improved for its efficiency, accuracy and convergence [22]; for instance, there are options for using different numerical dissipations for better stability criteria. Alternatively, to improve the robustness of the code, a symmetric Total Variation Diminishing (TVD) scheme was implemented into the original ARC3D. This modified ARC3D has shown to be more efficient and robust than the original one [23-26], and hence is adopted in this study.

The Beam-Warming scheme and the TVD scheme were constructed by totally different philosophies; however, in applications they can be viewed as the same scheme with different numerical dissipation terms [27]. Both schemes will be described and discussed in the following.


25








26

3.1 Beam and Warming Scheme


The Beam-Warming scheme is an non-iterative, implicit approximate factorization finite-difference scheme. The implicit operator is factorized (spatially split) to obtain an ADI-like ( Alternating Direction Implicit ) algorithm such that the three-dimensional inversion process is reduced to three onedimensional processes. The basic algorithm is second-order accurate in both time and space [28]. The implementation of this scheme to the system of governing equations, Eq.(2.21), is described as follows.


3.1.1 Temporal-Differencing


To approximate the time derivative of Eq.(2.21), the well-known secondorder Trapezoidal Rule or Crank-Nicolson Formula is employed to give

Aq" 1 aq"*l 8q"
- = ( + q) + O[(Ar)2] where Aq'=qn+l-qn (3.1)
Ar 2 ar a r

If one rewrites Eq.(2.21) as

aq aE aF aG - aS
ar a ai a as

and then substitutes it into Eq.(3.1), the following equation results
Aq r aE aF .kIGRe(aS
,q" - { --- + - + -Re-'( --)]"+1
2 a an a a

aE aF aG 1aS
+ [ + + --Re-'( )]"} + O[(Ar)3] (3.2)
a a?? a a

In the above equation, the flux vectors E, F, G, and viscous flux vector S at the advanced n+1 time step are in general nonlinear functions of the unknown








27


vector q. They are linearlized while maintaining second-order time accuracy as follows.

The flux vectors are linearlized by using a Taylor series expansion about q" as follows

a E
E"n+1= E"+AnAq" + O[(Aq)2] where A =
aq

F"*1 = F"+BAq" + O[(Aq)2] where B (3.3)
aq

G"*I = G"+CnAq" + 0[(Aq)2] where C = a aq

The above approximation is second-order in time since

Aq" = q"*1-q" = O(Ar) and O[(Aq)2] = 2[(Ar)2]

The viscous flux vector S is linearized by using the procedure suggested by Steger [291. The coefficient of viscosity y and thermal conductivity k are assumed to be locally independent of Jq ( recall that Jq is the unknown vector in physical domain ), so that elements of S can be expressed in the general form

si = J-1, ai where ct is independent of Jq.
a

Each element is then linearized in the following manner by Taylor series expansion

sin+I = s "+J-1{C : )"]}JAg * + 0[(Ar)2
a aJqj

Written in vector form, the viscous flux vector becomes

S"*I = S"+J-'M"JAq" + O[(Ar)2] (3.4)








28

The matrices A, B, C and M are so called Jacobian matrices and the elements of these matrices are given in Appendix A.

After substituting all the linearized terms, Eqs.(3.3- 3.4) into Eq.(3.2), one has

[I+ AT( + + -Re-, a (J-IMJ)"Aqn
2 ac an a a

= -Ar( aE + aF + --Re-, )" + O[(Ar)2] (3.5)
a8 an a a

where I is the identity matrix, and Aq" is the delta form unknown to be solved for. Note also that a first-order time-accurate Euler implicit scheme can be obtained by simply replacing Ar/2 on the left-hand side of Eq.(3.5) by Ar [281.


3.1.2 Approximate Factorization


Solving Eq.(3.5) directly is very inefficient since it involves the inversion of a very large system resulted from three dimensions. Beam and Warming applied an approximate factorization procedure to reduce the inversion to three onedimensional inversions. They rewrote Eq.(3.5) in the factored form

Ar 8A Ar B Ar a
+ ]"[+ ]"[+ A -Re-' --(J'1MJ)]"Aq"
2 a 2 an 2 a a

= -Ar( aE + aF + -G-Re-, as )" + O[(Ar)2] (3.6)
a an a a

Since the error introduced by this factorization procedure has the leading thirdorder term

Ar3 aA aB aB aC +C 8A ) aq"
4 a8 an an a a a ar

the second-order accuracy of the difference approximation is not affected.








29


3.1.3 Spatial differencing


The spatial derivatives appearing on the left-hand side of Eq.(3.6) are approximated by second-order finite-difference formula. e.g.


- 6f + O[(A )2I _ 1-(f-fi.1) + O[(A )2

hence the resulting scheme possesses a block-tridiagonal structure that can be solved efficiently. On the right-hand side, the same second-order differencing is employed for the viscous term, while a fourth-order finite differencing approximation, e.g.

af 6 - + O[(AI)4] (fi+2+8fjar-8fj-1+fi-2) + O[AO)41
at12Aq~

is used for the convection terms to keep the convection truncation error from exceeding the magnitude of that of the viscous term.

After applying spatial differencing into Eq.(3.6), the spatial derivatives are replaced by difference operators, and the finite-difference form of Eq.(2.21) is then

[I+ "[I+B "+ (6,C-Re-'s,(J-'MJ))]"Aq"

= -Ar(6,E+6 1F+6cG - Re-16,S)" (3.7)

This is the basic Beam-Warming algorithm. The implicit operator has been decoupled to produce an ADI-like scheme, and due to the second-order central-difference operators employed, the algorithm produces block tridiagonal systems for each spatial direction.








30


3.1.4 Numerical Dissipation


The Beam-Warming scheme requires the use of numerical dissipations to damp out the spurious oscillations that occurs when dealing with discontinuities. Several different numerical dissipations have been investigated and discussed in [30]. Those employed in the original ARC3D code were proposed by Beam and Warming [28], and Steger [29] and are described as follows:

The second-order implicit dissipation operators

De = -e1J1VIVeJ

D,, = -eJ-Yv A,,,J (3.8)

D. = -e1J-VCASJ

are added to the implicit operators on the left-hand side of Eq.(3.7) in the , and directions, respectively, and on the right-hand side the following explicit fourth-order dissipation term is added:

DE = -EEJ_1V [AvCA2)+ (VA7)2 + (V Ar)2 ]JAq' (3.9)

in which symbols A and v are the forward and backward difference operators, e.g.

af - --Af, + O[Af] = (fi-1f.) + O[Af]


f - VA + O[Af] =-(f1-f1.1) + O[Af]

and eI, EE are constant coefficients used to control the amount of numerical dissipations. It was suggested to set EE=At and eI=2EE so that as the time step increases, the numerical dissipations added relative to the spatial derivatives of convection and diffusion remains constant.








31

After inserting numerical dissipations into Eq.(3.7), the finite-difference representation of Eq.(2.21) becomes


[I+ A6A+D]"[I2+ 6,B+D ,]"[I+ Ar (6 C-Re-'6,(J-'MJ))+D ]"Aq"

= -Ar(6,E+6,,F+6,G - Re-'sS)"+DE (3.10)


3.1.5 Solution procedure


For simplicity, the following notation is introduced

Le = [I+ 6,A+D]


L,7 = [+ .6,,B+D,]"l


Lr = [I+ A(6C-Re-'6 (J-'MJ))+Dc]"Aq"

R" = -Ar(6 E+6 ,,F+ 6 G-Re-16,S)"+DE

where the L's are operators, then the Beam-Warming scheme, Eq.(3.10), can be written as

L L,,LAq" = R' (3.11)

If one further defines intermediate solutions as

Aq' = L rAq"

Aq* = L,,Aq*

the unknown vector q at the n +1 time level can be computed as follows.

LcAq** = R"

L,,Aq* = Aq*

LcAq" = Aq* (3.12)








32


qn+ = q"+Aq"

The Beam-Warming scheme has been shown to be accurate and efficient for simple flow problems (10], however, for more complex flow cases it can be very sensitive to grid resolutions [25].


3.2 TVD Scheme


It is known that finite-difference schemes may produce oscillations when applied across a discontinuity. For instance, the second-order Lax-Wendroff scheme

un+ I= iuj -v( u",.-u,' )-Hv(1-v)( u"+-2u"+uj' ), v = aAt/Ax (3.13) for a linear hyperbolic equation

au au
-- + a-- = 0, a > 0.
at ax

can produce spurious oscillations near the discontinuity. The Total Variation Diminishing (TVD) scheme is one of the methods designed to eliminate those undesirable oscillations. The notion of TVD was first introduced by Harten [31]. He derived a set of sufficient conditions from which the constructed schemes will not generate oscillations across discontinuities. Since the theoretical background of TVD schemes is rather complex and involved, only relevant basic definitions are given here. The construction of a symmetric TVD scheme and its application to thin-layer Navier-Stokes equations is discussed in the subsequent sections.








33

A numerical scheme, either explicit or implicit, can be written in the following operator form [32]

L U"*1 = R u" (3.14)

where u"+1 is the mesh function at the advanced time step to be solved, u" is the known mesh function at current time step, and L and R are some finitedifference operators designed for the particular scheme, e.g., for the second-order Lax-Wendroff scheme given in Eq.(3.13)

L 1 R = 1
R =1 -vAj+4-Mu2(1-v,)Aj.,A Aj

Here and in the following, ey stands for the central difference operator, i.e., Aiu = uj+i - uj. The total variation of a mesh function u' is defined to be

TV(u") = .z I u'y1-u"3I = . I A+ u" (3.15)
J=- _Wi =Here, the notation - is used to emphasize that the summation is taken over the entire domain. If a numerical scheme applies to a mesh function u, the total variation of the mesh function is non-increasing as a function of time, i.e.

TV ( u"+1 ) s TV ( u" ) for all n

the scheme is said to be TVD. In practice, the following sufficient conditions derived by Harten [31] are served for constructing TVD schemes:

TV ( R u" ) TV(u )

TV ( L u"*' ) TV ( u"') (3.16)

That is, if the numerical scheme shown in Eq.(3.13) satisfies the sufficient conditions, Eq.(3.16), the scheme is TVD.








34


It has been proved that a TVD scheme is monotonicity preserving [32], that is, for any monotone mesh function u that increasing (or decreasing) monotonously with time, w = Qu is also monotone, where Q is a finite-difference operator. Since a monotonicity preserving scheme will not generate spurious oscillations, and neither will the TVD scheme.

In what follows, the construction of the symmetric TVD scheme for a onedimension hyperbolic equation will be described first, then the scheme is applied to the Euler equations, and finally it is extended to the thin-layer Navier-Stokes equations.


3.2.1 Symmetric TVD Scheme


The symmetric TVD scheme was developed by Yee [27] and it is a generalization of Roe [33], Davis [34] and Sweby's [35] TVD Iax-Wendroff scheme. To illustrate the construction of this scheme, a one-dimension linear hyperbolic equation of the form

8u au
-- + a - 0 (3.17)
at ax

is considered. Sweby has shown that the second-order Lax-Wendroff scheme for this equation can be constructed by adding an anti-diffusive flux to the first-order upwind scheme [35], i.e.,


uji= u1 - vIgAu" - M.(l-v)A.)Aj+%u", if a > 0.


i4 + = i - Lj+%u" - v(1+v)j+%A. u" , if a < 0.


(3.18)








35

where v = aAt/Ax is the Courant number. Then he introduced a limiting function p and rewrote the above scheme as


11+1 = uj - V 1,+(1-v)[p(rj')/r,*-p(r'j.1)]}ag.gu" if a > 0 (3.19)
114+ = u1 - v{1-%(1+v)[p(rg-,)-p(r-j)/r-j]}Aj,4u" if a < 0
where p is chosen to be a function of the ratio of consecutive cell gradients rj; p, = p(rj), r,+ = Aju/A,.+,u, r; = A+,u/A-u, with the restrictions (sp/rs2, 0sp:2 to satisfy the TVD constraints, Eq.(3.16).

The Sweby's scheme, Eq.(3.19), was reformulated by Davis [34] as a LaxWendroff scheme plus an upwind weighted numerical dissipation term, i.e.


14+1 =j - 1/(1+v)A u"- /hv(l-v)Aj u" +

[K (rj+)+K-(r- 1)]Aj.,u" - [K (r*j-,)+K-(r-)]Aj.u" (3.20)

where

K v(1-v)[1-p(rj+)] if a > 0
1.0 ifas0

K-(r) 0 if a <: 0 (3.20a)
v(l+v)[p(r- )-1] if a < 0
To avoid determining the upstream direction, Davis then modified the dissipation term as follows.

K+(rj+) = jv|(1-|vj)[1-p(r')]/2

K-(rj ) = Jv|(1-|vI)[1-p(r;)]/2 (3.20b)

Shortly after that, Roe [331 further generalized Davis' scheme, Eq.(3.20), and suggested a new limiting function Q. Roe's scheme takes the form








36


uj+7= uI - 2v(l+L)A.4u" - %v(1-V)AjU"

- Vvj(1-|Vj)(1-Qj.4)Lj_%U" + V j (1-jvj)(1-Qj %)Aj 4U" (3.21) where the limiter Q depends on three consecutive gradients given by

QJ+% = Q(Aj.-u/Aj+%u , Ay1,u/Aj+,u) = Q(r*j , r- 1) (3.21a)
and

0 < Q < 2/(1- v)

0 < Qyg/r-,i < 2/Ivl (3.21b)

0 < Q3+,/r*j < 2/1&,l

to satisfy TVD sufficient conditions, Eq.(3.16). There are a variety of limiter Q satisfing the above restrictions, such as

Q(r',r+) = max[O, min(2r-,1), min(r-,2)]

+ max[O, min(2r*,1), min(r*,2)] + 1 (3.21c)

Roe's scheme was later modified and extended to a one-parameter family TVD scheme by Yee [27]. She first rewrote Roe'scheme, Eq.(3.21), into the form

uj*1+ Ae(h,'%-hj'g1) = uj" + A(1-e)(hj,-hj",) (3.22)

hj+ = [a(up+1+uj)-{Aa2Qj+&+IaI(1-Qj+4)}Aj+4u]/2, A = At/Ax (3.23) and then simplified the numerical flux to

hj+4 = [a(uji+uj)-ala(l-Qj, .+ u]/2 (3.24)

Note that the term Aa2Qj+%Aj+%u/2 was taken out from Eq.(3.23), but it makes no difference as long as the limiter function Q is chosen such that the scheme satisfies the TVD sufficient conditions.








37

For nonlinear hyperbolic conservation law of the form

+ af(u) = 0 (3.25)
a t a x

the numerical flux, Eq.(3.24), was modified to

h = [(f ,i+fj)-w(a , )(1-Qjg)Aj 4u]/2 (3.26)

where f=f(uj), A + u=uj,-u+ , a(u)=df(u)/du is characteristic speed and

aj+ { (fj+rffj)/j+%u , - 0 (3.26a)
a(uj) , Ay+u 0

The function F(a+,) that is sometimes referred to as the coefficient of numerical viscosity, is defined as
{ |zI , IzI 3 2 b
'i2(z) = {(2+z2)/2c ,z< (3.26b)

to avoid the entropy violating problem [31], where E is a parameter that is described in Ref.[36]. In order for the above scheme, Eqs.(3.22, 3.26), to be TVD, that is, satisfying the TVD constraints, Eq.(3.16), the following conditions obtained by Yee [37], must be satisfied

0 < Qi+ < 2

0 < Qj+ /r*j < 2/[v(1-9)a,.g,] - 2

0 < Qi+4/r-;+1 < 2/[L.(l-e)Ia,+1I1 - 2 (3.26c)

v Iaj+, I< 2/[3(1-e)]

If one selects the following limiter that satisfies the above constraints

Qi% = minmod(Aj.u , Aj ,u) + minmod(Aj ,u , Aj 1,u) - Au+u (3.27)
then from Eqs.(3.26-3.27) Yee's numerical flux can be expressed as

hj+%= [fj + fj + j+gJ/2 (3.28)

O+ = *(a,+)(gj+gj+1) - 2,&4.,u (3.28a)








38

where the minmod function minmod(xy) is defined as

0 , if x y have opposite sign
minmod(x,y)= 9

and

gj = minmod(t-.,.u , j.,u) (3.28b)

The above scheme, Eqs.(3.22, 3.28), form a class of TVD schemes with the numerical dissipation terms centered and hence are often referred to as symmetric TVD schemes. It is noted here that the scheme is TVD for onedimensional nonlinear scaler hyperbolic conservation laws, and is also TVD for constant coefficient hyperbolic systems. However, in applications, it is formally extended to more complex equations, such as the Euler equations, the thin-layer Navier-Stokes equations, and evaluated by numerical experiments.


3.2.3 Extension to Euler equations


Extension of the scaler TVD scheme to systems of conservation laws is not unique. Here, one follows the method proposed by Yee [38,39]. The method is accomplished by defining at each point a "local" system of characteristic fields and then applying the scheme to each of the scalar characteristic equations and hence is referred to as the "local characteristic approach".

Consider the inviscid part of the thin-layer Navier-Stokes equations, Eq.(2.21),

aq aE aF aG
-- + - + -- +--- = 0 (3.29)
at a an; a








39

Recall that the Jacobian matrices of E, F, G are A, B, and C, respectively. Denote the eigenvalues of A, B, and C as


(al a ai '), a~ia~ia) (4i,agaaa")
and Rt, R, R, as the eigen-matrices whose columns are eigenvectors of A, B, C, respectively, and R-'E, R',, R-, as the inverses of the corresponding eigenmatrices. The explicit form of these eigenvalues and matrices are given in Appendix B.

In the -direction, let q,,, be some symmetric average of %, and q,, e.g.

qj+% -= 0.5(q 4+q,,W)
or Roe's average

% + = (P gjq*+ p 4j1,jqy 1W)/(PWh + Phy ) (3.30)
and denote a',+, R,, R1j+% as the quantities a',, R,, R71, evaluated at qAk. The difference of the characteristic variables in the local -direction are defined as

aj+4 = R+'j g(J ,,Wqj +,-Jgq+)/[(Jj, +J J)/2] (3.31)

One can also define the corresponding parameters for the r7 and directions in a similar manner.

Using these notation and applying Eq.(3.22) to each equations, a oneparameter family of TVD scheme for the system of equations can be written as q + Ae[(Ej";* j-E _* F + " F") +(Gj" -Gj')] = qJk - \(1-9)[(Ej",sW-E _4,)+(F "k -Fjl_%, +(GjL+ 4-Gj"0.%)] (3.32)








40


where A=Ar since A =Ar =A =1 and the numerical flux functions obtained from Eq.(3.28) are given by

Ej+* = [Ej j+Ej+Ryg4Ds]/2

Fj,k+, = [Fgky+Fj,k+lI+Rk+ 4Dk+ ,]/2 (3.33)

Gj4.,.% = [Gjjo+Gj,+RaseggD,]/2
the expression of the elements of Dj.,. are obtained from Eqs.(3.28a-3.28b) and are given as

0 += t(aj+,)(egj, -2aj+,) (3.33a)

g'. = minmod(4x. ,a.+%) (3.33b)

Likewise, one can define Dk+4, II+4*


3.3.6 Linearization and ADI decomposition


In order to solve Eq.(3.32) efficiently, the left-hand side of Eq.(3.32) is first linearized by the procedure described in Section 3.1.1. The resulting scheme written in delta form is

[I+ Ae[(Hj+ %*+-Hjg +(Hj+j-Hjk-.,)+(Hjg g-H p.g)]Aq" = -A[(Ej +,4-Ej.-,1) + (Fj+,.j-Fk,1) + (Gjg,+4-Gj%)]" (3.34)

where I is a 5 X 5 identity matrix and the H operators are given by:

Hj+%,, = [A +n s,]/2

H.,k,, = [Bjk+1 j,k+4,I]/2 (3.34a)

Hjgj+ = [C +~ogj,I+]/2

in which

nj+.k = {R diag[p'-2,y(a')] R-1}+qA+%








41


Qk+m, = {R diag[p'-2k(a')] R-1}k+4Ak+4 (3.34b)

aj4, = {R diag[p'-2,Y(a')] R-1,,.A and

flj+ = *(aj 4)(gi +gjj )/ct'j ,

+= k k+-Jdk+ek+0/aik+A (3.34c)
Si.. =+ *(a',.j,g(g e,ig1a'j ,

Here, the expression diag(z) denotes a diagonal matrix with diagonal elements z'.

To calculate the o's at every time step is quite costly. For steady-state

applications, the temporal accuracy is not important. One can use a spatial firstorder implicit operator (that is, by setting the limiters equal zero); i.e., by redefining the n's as

nj+%,v = {R diag[-*(a)] R-1})+,+

ak+%4 = {R diag[-*(a')] R-1}k+%Ak+% (3.35)

njW+4 = {R diag[-i(a')] R-1}1,A1+

Since the temporal accuracy is not important, the computation can be further reduced by simplifying the n's diagonal matrices, i.e.

n+v.*ji = {diag[-maxt(a')]},A+,

ok+v4 = {diag[-max(a)]}k+%Ak+% (3.36)

nj., = {diag[-maxY(a')]},Asa

If one applies the approximate factorization procedure discussed in Section

3.1.2 to Eq.(3.34), the ADI form is obtained, i.e. [I + Ae(Hj+ %,-H.%,k)][I + ) (Hjk.,.-Hjk-&,)][I+ A(H ,.%-Hj.%)]Aqn = -A[(Ej+%*k-Ej_. ) + (Fjk-+%I-Fj,k-,) + (Gj,+ -G .%)]" (3.37)








42

With the use of Eq.(3.33) and Eq.(3.34a), the above equation can be expressed in delta form as follows

[I+ Ae6,A+D,]"[I+ Ae6,,B+D,,]"[I+ As6,C+D,]"Aq"

= -A(6fE+6,,F+6G)"+DE (3.38)

where

D=

D,= (3.38a)

D=

DE = -,+R-R.i+RR + s-14%4-R1_%oj.%) (3.38c)


3.3.7 Application to thin-layer Navier-Stokes equations


To apply the TVD scheme to the thin-layer Navier-Stokes equations,

Eq.(2.21), is similar to that to the Euler equations [38-39]. The viscous term that not appears in the Euler equations is first central differenced and linearlized as discussed in Section 3.1.1 and then added to Eq.(3.38). The resulting finitedifference equation is

[I+ Ae6,A+D,]n[I+ +e6,,B+D,,]" xe(6,C-Re-'6,(J-'MJ))+D,]"aq"

= -A(6fE+6,,F+s G - Re-s 6S)"+DE (3.39)

It is noted here that A=At since A =Ar =A =1, and e is a parameter in

the range of zero to one; Euler explicit or implicit scheme is obtained by setting 9 = 0 or 1; if e =1/2, it is a second-order implicit scheme. However, the order of accuracy is degraded due to the simplifications that have been made in the derivation of the scheme. It is also noted that if 9 =1/2, Eq.(3.39) is similar to








43

the Beam-Warming scheme, Eq.(3.10), except that the numerical dissipation terms are replaced by the ones derived from TVD constraints. In fact, the TVD scheme can be viewed as a central difference scheme with a "smart" numerical dissipation which is an automatic feed back mechanism used to control the amount of numerical dissipation [38].


3.3 Flow Solver


Applications of the original version of ARC3D and its axisymmetric

version [40] for the simulation of transonic turbulent flow past a projectile model showed that the codes implemented with the Baldwin-Lomax turbulence model can give surface pressure distribution in excellent agreement with measured data if good grid networks are provided to the Navier-Stokes codes for the simulation [25]. Numerical experiments also showed that the Beam-Warming solution algorithm is very sensitive to grid characteristics in complex flow regions. In order to improve the robustness of the codes, a symmetric TVD scheme was first extended for modifying the axisymmetric code and investigated for the projectile aerodynamics simulation; as reported in [25,41], the modified code is indeed less sensitive to the grid characteristics in complex flow regions. Later on, the original ARC3D was modified to implement the symmetric TVD scheme. As have been mentioned before, from the viewpoint of numerical implementation, the only difference between the Beam-Warming scheme and the TVD scheme are the numerical dissipations. Hence, to implement the TVD scheme into the original ARC3D, one simply replaces the dissipation operator, Eq.(3.8), and the








44

dissipation terms, Eq.(3.9), with those designed for the TVD scheme, i.e., Eqs.(3.38a, 3.38b). However, the latter is more complicated than the former. New subroutines for computing the TVD dissipations had been written and the original ARC3D was then modified into a TVD thin-layer Navier-Stokes code [26]. This TVD code has been applied to the projectile flow problem, and the results showed that the code was indeed accurate and robust. Hence, it is adopted in this study for developing the multi-block computational method.













CHAPTER IV
THE FLOW PROBLEM AND SOLUTION STRATEGY


A transonic turbulent flow of Mach 0.8 past a wing-fuselage configuration shown in Figure 4.1, is adopted as the model problem to develop the multi-block computational method. The flow problem considered is indeed a complex one, involving many viscous phenomena, such as viscous sublayer, shock boundary layer interaction, flow separation, as well as complicated boundary geometry. There are wing surface pressure measurements available for assessing the accuracy of the numerical solutions. The flow field is assumed to be symmetric with respect to the symmetry plane of the wing-fuselage configuration and hence only half of the flow domain is required for the simulation.


4.1 The Flow Problem and Preliminary Work


The wing-fuselage configuration shown in Figure 4.1 is a model of a nextgeneration 150-passenger transport. The wing, optimized for Mach 0.77, has an aspect ratio of 5.17 based on the averaged chord length, wing leading edge sweep of 30*, and an advanced supercritical airfoil. This wing-fuselage configuration is used as a model to develop the computational methods for propfan powered aircraft [42]. A rather coarse one block H-grid for the configuration was provided by the NASA Ames Research Center. The flow field is assumed symmetric with respect to the z = 0 plane and hence only the z > 0 half-space is


45








46


Ilk ~ i':


Figure 4.1 Wing-fuselage configuration


.,wi%








47

considered for the simulation. Figure 4.2 shows an overview of the physical domain and its computational domain.

In the following, the integers j, k and 1 are used to denote the grid points in the , 7, and directions, respectively. Here, is along the streamwise direction, t is along the wing spanwise direction, and is about normal to the solid surface. The dimensions of this original grid network are jmax = 101, kmax = 45, and lmax = 29. The solid surfaces are transformed into part of the 1=1 plane with the shape H, as shown in Figure 4.3; k =1 and k = kmax are symmetric planes; while j =1, j =jmax, and I= lmax are artificial outer boundaries.

To further illustrate the configuration quantitatively, one takes the

maximum diameter of the fuselage, D, as a reference length. The fuselage is about 7D long; the wing span is about 3.3D; the chord length at the wing root is about 1.4D; the diameter of the sting is about 0.075D; the distance of the outer boundaries far upstream and downstream are about 15D, while that of the lateral side is about 10D. These relative lengths and the j, k, I triplet in physical space are shown in Figure 4.4. Note also that the first spacing next to the solid surface is in the range of 0.008D to 0.026D.

The 101x45x29 H-grid one-block system was first used for a preliminary study. Because of the specific grid topology, the thin-layer Navier-Stokes code simply can not be applied directly without modification for turbulent flow simulation. As can be seen from Figure 4.3, only part of the 1=1 surface is the solid surface. Extensive modifications were made for proper Jacobians and metrics calculations, boundary conditions, and eddy viscosity computations. Also, 30 grid points that extremely close to the solid surface were added to the original










48












9-2C





cF
E L JK45


z- J-1 D






A






00
Qp
Y GH



. C


Figure 4.2 An overview of the original grid










49


(a) Physical space





0 U v PH


NOSE FUSELAGE TAIL STING
0 R F

WING


K L


WING


E 0 R F

NOSE FUSELAGE TAIL STING


0


S


(b) Computational space


Figure 4.3 The 1=1 surface of the original grid


K L















F


E 0
z S,


T


w












8D

I 1,23.,1)
3.GD ..C(1.,45.,29)
37 ,23,1)
E-15E K L 76 23.,
,29)
0 p J -t,2131, 1
(13,l.,l) 7D
Q
CA(90 1 )10D
Z37D' 15D F D
GH (101, 45,29) y B 1OD 0.075D
x (101,1l,29)








51

grid to resolve the thin viscous-layer; and two mirror image surfaces were also added for ease in symmetric boundary condition treatment. The 101x47x59 Hgrid network was then provided to the modified thin-layer Navier-Stokes code for the flow simulation. As shown in Figure 4.5, the computed results are in poor agreement with experimental data. This seems to indicate that the number of grid points on the wing surface, 18 in the spanwise direction and 40 in the chordwise direction, is not fine enough to resolve the complex flow phenomnia. Also, modifications on the thin-layer Navier-Stokes code may be needed before the one-block H-grid topology can be used for this complex turbulent flow simulation. A flexible and efficient multi-block technique is more desirable to numerically simulate a complex geometry turbulent flow problem, such as the wing-fuselage aerodynamics. It is also noted here that the calculations of this one-block approach were performed with the Cray X-MP/48 of the NASA Ames Research Center.


4.2 Solution Strategy


For flow past a complex configuration, such as the wing-fuselage

combination, to generate an acceptable single grid for the entire flow field is rather difficult. It is much easier to divide the flow field into several subregions and then generate a good grid for each subregion. Hence, it is natural to develop a solution method based on the multi-block grid network by treating each block as an independent flow problem. Each block is then solved by the same thinlayer Navier-Stokes code with the information between blocks transferred via the interface boundary conditions.









52


/P 0
M -


M - 0.8
-000 00 o0


0


o -


0



a.-.4
uym 6-


I


(
/ /
/ / I
/ /
/
/
/
/


/
-I








/1


Figure 4.5 Computed results of the one-block approach


CP.


I








53

For the wing-fuselage flow problem, an effective and yet simple block-byblock computational method is proposed. In this method, the entire flow field is divided into six contiguous blocks and each block is partially bounded by a solid surface. Each block is then treated as an independent flow problem, which is solved by the same thin-layer Navier-Stokes code, with the interface boundary conditions updated at every time step. This method does not need to worry about patching solutions between different equation zone as in the zonal equation approach, and it is easier to deal with the interface boundary conditions than that of the Chimera overset grid approach. However, in the same manner as the zonal methods, poor interfacing can introduce errors at the interface boundaries and consequently propagated into the solution field. The error introduced can depend on the solution behavior and the grid characteristics around the interface as well as on the interpolation technique employed for updating the interface boundary conditions. For a successful multi-block computational method, it is then important and critical to minimize the associated interface error and its effect on the solution accuracy.

To attack the wing-fuselage flow problem by using the proposed method, a flow simulation package has been developed which includes four independent computer codes; a surface grid generation (SG3D); a volume grid generation (GG3D); a multi-block flow solver (WF6B); and a plotting program (PP3D). They will be briefly discussed in the subsequent chapter; however, a detailed description can be found in a reference manual [43]. Most of the figures presented in this dissertation were produced by the plotting program PP3D.








54

4.3 Six-Block Grid Topology


For a multi-block computational method, the definition of a good

composite grid topology for a complex configuration turbulent flow simulation is more involved and difficult than that for an inviscid flow or a laminar viscous flow simulation. For instance, the solid surface has to be transformed onto an entire plane in computational space so that the thin-layer approximation can be applied easily; the Baldwin-Lomax turbulence model requires the normal distance of each grid point to the solid surface; and the grids should be able to cluster points close to the solid surface to resolve the thin viscous-layer. Hence, for ease in eddy viscosity calculation and in applying the thin-layer approximation, a special 6-block grid topology was defined and investigated for the wing-fuselage configuration. The main idea of dividing the flow field is to isolate the wing by a cone-like boundary surface that sits between the wing and the fuselage, as shown in Figure 4.6. The isolated wing was divided into 2 blocks and the rest of the flow field was divided into 4 other blocks. By doing so, the flow field is divided into six contiguous blocks such that each block is partially bounded by a solid surface, and the solid surface of each block is transformed onto = 1 plane in the computational space. An overview of the topology is shown in Figure 4.7 while Figure 4.8 gives side and cross-sectional views. The physical and computational space of each block is shown in Figure 4.9.

The selected grid topology has the clear advantage that each block can be solved by the same thin-layer Navier-Stokes code; however, the choice has its deficiences in the block grid generations. For instance, for the wing blocks, Block








55


zY
dividing surface



















z Y
ividing surface


Figure 4.6 The dividing surface isolating the wing








56


-o


I


NAV


CC /
c0


C-C UtJ


Figure 4.7 Six-block grid topology (overview)









57


Figure 4.8 Side and cross-sectional views


H S
T N





F BOC 1\C /BLOCK 5 x
A B K Z G M
L x U


z

H,S,T,N



b' BLOCK 3 BLOCK 4,


R,Q

BLOCK 2 BLOCK 5
b 0

G, m B,K D,P 0,1










58


H















E A
B GB

BLOCK I


Figure 4.9 The transformation of the six block grids


z


y











F
G Ht

C

AE
DD

B C D


(a) Block 1

H
z

y












N G N



CI

B b
CL
2 -L K


(b) Block 2


G - - - - - - - -B-
M BLOCK









59


S

N H~f
)I/
fitI
\_e
QR
C
C L
bL


(c) Block 3



Y












H

N
T


R C
L
BLOCK 4

(d) Block 4


Figure 4.9 ( continued )


B







BLOCK 3


s






T N






C
L


AH


s







T





(e) Block 5


z

y -- x













x z U y K

L K


(f) Block 6


Figure 4.9 ( continued )


60


0


z


N










U



P K



BLOCK 6


C
L p-- - -- -8125


H
00
N



C

L








61

3 and Block 4, shown in Figure 4.9(c,d), four of the six bounding surfaces in the computational domain, CLRQ, LNTR, RTSQ, and QSHC form a nearly flat surface in the physical space. Therefore, special measure and care in grid generation are required to overcome the difficulty resulting from the drastic change in the domain transformation. Also, for the nose block, Block 1, shown in Figure 4.9(a), the boundary plane AEJF in the computational domain degenerated into a single line in physical space and care must be taken when dealing with the boundary conditions.













CHAPTER V
GRID GENERATION


For the wing-fuselage flow problem, the flow domain has been divided into six blocks, however, it is still a great challenge to generate acceptable grids for each block. To generate grid networks for the six-block wing-fuselage flow problem, surface grids must be obtained first. A special surface grid generation code (SG3D) was developed to utilize the original H-grid network for constructing the required surface grids. The resulting surface grids will be presented in section 5.1. Generating volume grids for the six-block flow domain, it is required to develop an adequate grid generation code. Various grid generation methods [44-51] have been studied in that regards. The main concern is that the code can produce grid networks that satisfy the required grid characteristics, such as: (a) grid lines in different directions must be nearly orthogonal to each other; (b) grid lines in c-direction must be about normal to the solid surface; (c) = constant grid surfaces must be clustered near the solid surface; (d) grid lines must be continuous across the interfaces; and (e) the blocks must be contiguous (no holes or overlapping). A combination of two-boundary algebraic and elliptic grid generation methods will satisfy these requirements best. Accordingly, a general purpose 3-D Grid Generation code (GG3D) that based on these two methods was developed and applied to generate the 3-D grid networks.


62








63

A detailed description of GG3D can be found in [43]. The techniques employed to generate the six block grids are described in Section 5.2 and the resulting grid networks are presented in Section 5.3.


5.1 Surface Grid Generation


To generate grid networks for the six block grids, surface grids must be obtained first. A two-boundary algebraic grid generation method is used to generate volume grids; consequently, only the solid surface grid (1=1) and the cooresponding outer surface grid (1=imax) are required. The other boundary surfaces are then defined by the straight lines that connect the boundary curves of the solid and outer surface grids.

A special surface grid generation code (SG3D) has been developed that uses the original H-grid to construct each surface grids by cubic splining, transfinite interpolation, strentching, shifting and dividing. The resulting surface grids obtained are shown in Figure 5.1; the dimension of these surface grids are 28x23, 40x12, 40x35, 40x35, 40x12 and 26x23 for Block 1-6, respectively.

It is noted here that the outer surfaces are defined according to the grid

topology and will have strong effect on the grid characteristics of the block grid to be generated. Hence, the boundary curves of each outer surface should be defined carefully. For block 1, a cone-like shape that is resemble to the solid surface is desired for the outer surface. As shown in the grid topology, Figure

4.7, the boundary curve GHI is defined such that the bounding surface lies in the middle of the wing leading edge and the fore-body of the fuselage. Similarly, the outer surface of block 6 is defined such that bounding surface of the boundary








64


NA


(a) Solid surfaces


Figure 5.1 The 6-block surface grids


lapr


6








65


fir


(b) Outer surfaces

Figure 5.1 (continued)


N








66

curve MNO lies in the middle of the wing trailing edge and the aft-body. The boundary curves Hb'N of block 2 and Ha'N of block 5 are also defined such that the boundary surfaces Hb'NLC and Ha'NLC lie in the middle of the wing and the fuselage. For the outer surfaces of block 3 and block 4, the boundary curve HN, that is defined as a straight line, is located such that the interface surface of block

3 and block 4, i.e., surface HRCQLN, is not incline to either upper side or lower side of the wing. In other words, the plane spanned by HCLN should pass throught the wing tip, and hence divides the wing into upper wing and lower wing. If this were not the case, one would have difficulty to generate an usable volume grid for either block 3 or block 4.


5.2 Volume Grid Generation


5.2.1 Two-Boundary Grid Generation


The two-boundary grid generation technique used is based on Hermite

transfinite interpolation [46-48]. Two surface grids are specified and the interior grid lines are connected by curves such that boundary orthogonality is enforced. The other four boundary surfaces are then defined using straight lines. The distribution of interior surfaces is defined by a suitable clustering function such that more grid points are clustered near the solid surface as it is usually required. A description of clustering functions is given in Appendix C.

Denoting the position vectors of the two specified boundary surfaces as

rQ,,0) = r1(C,q) and r(Q,r,1) = r2(Q,q) (5.1)

and the derivatives with respective to as








67

DR, = ar(Q, i,0) and DR2 = ar(Q,,1) (5.2)
a a

the following functional form, based on Hermite transfinite interpolation, may be written for the connecting curves:

r(,v, ) = r(4,t)hj( ) + r2(,i)h2() + DRih3(g) +DR2h4( ) (5.3)

To satisfy Eq.(5.1) and Eq.(5.2), the following must be true:

hl(0) = 1, hI(1) = 0, ah1(0) = 0, ahb(1) = 0 a a

h2(0) = 0, h2(1) = 1, ah2(0) = 0, ah2(1) = 0


h3(0) = 0, h3(1) = 0, ah3(0) = 1 ah3(l) = 0 (5.4)


h4(0) = 0, h4(1) = 0, ah4(0) = 0, ah4(1) a . a .

Hence, a cubic polynomial in for hj( ) satisfying conditions given in Eq.(5.4) can be determined. They are

h= 2 3-3 2+1 h2 =-2 3+3 2

h3 = c3-2 2+ (5.5)

h4 = _-2

Here, is a normalized point distribution, and is usually defined by a clustering function as shown in Appendix C.

The expressions for the derivatives, DR, and DR2, appearing in Eq.(5.3) are the key elements used to enforce the boundary orthogonality and are derived as follows.








68

A vector normal to the surface r = 0 that points into the spatial domain is given by

n - ar1(,y) ar1(Q,7)
a a 1

where x denotes the vector cross product. A vector that tangents to the connecting curve at r =0 can be defined as

e = ar(Q,2,O)
a.

Hence, in order for the connecting curves to be normal to the r = 0 surface, the cross product of the vector tangent to the connecting curves e and the surface normal n must be zero, i.e.,

e x n = 0 at = 0 (5.5)

or written in component form

ay ax ay ax ay az az ax az ax a a 1 a a a 1 a a n a a a n

+ ax ay ax ay ax az ay az ay az
a a? ae a a1 a an ae ae a n

+ ax ax az ax az ay az 8) 8z ay )]k
a a7 ae ae a17 a( a7 ae ae an
=0 (5.5a)

From the above equation, it can be seen that the connecting curves will be perpendicular to the surface r = 0 if the following conditions are satisfied at r = 0:

DX= - ( ay az ay az
a at ac ae an








69

DY, - ay = CK,(,o)( ax az ax az (5.7)
a at a a an

az ax ay ax ay
DZ1 = = CK1(,n)(
an a a an

Here, CK, are arbitrary functions introduced to control the shape of the connecting curves, and if set to zero the curves become straight lines.

In a similar manner, one can show that the connecting curves will intersect the - = 1 surface perpendicularly if

DX2 = x -CK2.Q,7)( ay a ya
an a a an

DY = = CK2(ci) --ax az -z) (5.8)
aa ay ax ay

DZ2 _ az = -CK2(, x)(
an a a an

are satisfied at g = 1.

By substituting Eqs.(5.7-5.8) into Eq.(5.3), one obtains the following

expressions for the resulting volume grid with boundary orthogonality enforced:

x(Q,n,,) = x( ,ii)hl( )+x2(,7)h2( )+DXlh3( )+DX2h4( )

y(,n, ) = yl( ,7)hl( )+y2( ,7)h2( )+DYlh3( )+DY2h4( ) (5.9)

z(,q, ) = z,( ,t7)h,( )+ z2(,tI)h2( )+DZh3( )+DZ2h4( )

It is noted here that the derivatives on the right hand side of Eqs.(5.7-5.8) are discretized by proper finite- difference fomula, e.g., central difference

ay1 _ yl(k+ 1)-y(k-1) + O(An2)
an n(k+ 1)-a(k-1)

and e and q are normalized as follows:








70

0(j) = (j-1)/(jmax-1), j=1,2,...,jmax

n(k)= (k-1)/(kmax-1), k=1,2,- - -,kmax (5.10)

The CK's appearing in Eqs.(5.7-5.8) have a strong effect on the shape of

the connecting curves. In fact, they control the depth of the perpendicular part of the curves. Values for the CK's are usually specified as constants by trial and error. A shape function was devised to modify the CK's so that different values can be specified at different j, k locations. The shape function is of the form

CK = CKI*cosw{CPJ[(j-1)/(jmax-1)]+CSJ}

*cosr{CPK[(k-)/(kmax-1)] + CSK} (5.11)

where CKI is a constant to be specified, CPJ and CPK are used to control the frequency while CSJ and CSK are used to control the phase of the relative cosine function. In doing so, values for the CK's become functions of j and k, and can be smoothly distributed on the surface. For example, if one specifies CPJ=0.5, CSJ = 0, then CK = CKI when j =1 and CK decreases as j increases; or if one specifies CPJ =1, CSJ = -0.5 then CK increases from zero to CKI and then decreases to zero as j varies from 1 to jmax. Figure 5.2 shows some examples to illustrate this shape function.


5.2.2 Elliptic Grid Generation


The technique and subroutines employed are based on the elliptic grid generation system of Thompson [49-51]. The volume grid is constructed by forcing the Cartesian coordinates of interior points to satisfy specially devised Poisson equations. In the Poisson equations, three control functions for each curvilinear direction are set up to control the grid point distributions. The system











71


CPJ-0. 5, CSJ-0 ,CPK-0,CSK-8


CPJ-0,CSJ-0,CPK-.6,CSK--.4


CPIJ-.8, CSJ--O. 4, CPK-.5,CSK--.5 CPJ-0, CSJ-e, CPK-8,CSK-0


CPJ-.5,CSJ--.5,CPK-.5,CSK--.5


CPIJ-1,CSJ--.5,CPK-.5,CSK--.5


CPJ-.6,CSJ--.3,CPK-.6,CSK--.3


CPJ-1,CSJ--.5,CPK-1,CSK--.5


Figure 5.2 Shape function for CK's








72

of Poisson equations are solved by point SOR (success overrelaxation) iterations. The user-provided initial volume grid is used not only to start the SOR iterations but also to set up the control functions.

The elliptic grid generation system is based on the Poisson equations of the form

V'= (vs.v )Pi

V = (vn -vn)P2 (5.13)

V2= (v .v)P3

where the P's are so-called control functions which can be specified to control the grid smoothness and/or orthogonality. Since the actual computations are performed in the computational domain where , Y7 and are the independent variables, with r = (x,y,z)T as the dependent variables. The poisson equations, Eq.(5.13), are then transformed to computational domain. The transformed equations are

(V. 2 r+(I.V a2 r V. 2r
a2r8r a~r
(v 7vg) a + (v7.v) + (v vI)- +
agoaano agag

a2r a2r a2r
(v .v-,) + (v7-v-)- + (v .v7) - +

aga ana aga
(v-v)- + (vn-.vg) + (v-v)- +
ar ar ara
(v-v )P Ia +(Vn-Vn)P2 a+(VV,)P3 - 0 (5.14)
a a?7 a

The finite-difference equations of Eq.(5.14) are obtained by replacing the derivatives with central differences and these equations are then solved by point SOR iteration procedure.








73

If the control function P's are evaluated from the initial grid, the three components of the elliptic grid generation system, Eq.(5.14), provide a set of three equations in which the control functions P's are treated as unknowns and the values for the position vectors are obtained from the initial grid. The transformed poisson equations, Eq.(5.14), can be rewritten as

(V -V )Pj a + (V17- V1)P2 ar +(v,-V,)P3 arar a1 ar

-[( V ) a2r + ( V O a') + (V2.V,) + aga alna asa
8~r 8~r a~r

( -V 7) ar + (vi.vE) + (v .v7) - +
ar ar ar

(v .v) + (v97.v7) + (v -v,) + (5.15)
aa a1W aar

which can be solved simultaneously at each point for the three control functions, Pi, i= 1,3. If the control functions are computed solely based on the above equations, and the elliptic system, Eq.(5.14), will reproduce the same grid as the initial grid which is of trivial interest. However, if the computed control functions are smoothed before they are applied into the elliptic system, a smoother grid network can then be produced.

A variation to the evaluation of the control functions is obtained by assuming that the grid lines are orthogonal to each other, and Eq.(5.15) is reduced to

(v -v )P - ar +(v"-v)P2 ar +(vr.vO)P3 ra a17 as

ar ar a~r
-[(v -v) + (v.v) + (v .v ) -] (5.16)
aa8 a017 asa








74

The control functions evaluated from these equations will force the elliptic system to produce a more orthogonal grid.


5.3 The Six-Block Grid


A set of volume grids with 60 grid points in the 1-direction was first

generated by two-boundary grid generation with grid lines set normal to the solid surface whenever possible. The first spacing off of the solid surface was set to about 2x10-5D, which is considered fine enough for turbulence simulations. Figure 5.3 shows the 3-D close-up plots of the six block grids with cut-offs to show the details of the interior surfaces. These grids were generated from the solid and outer surfaces described in Section 5.1. The other four boundary surfaces were then defined by straight lines that connect the boundary curves of the solid and outer surfaces. Sixty grid points were distributed on each of the connecting straight lines by using the hyperbolic tangent clustering function with the first spacing of 2x10-5D and the last spacing of six times of the average spacing, and hence, more grid points were clustered near the solid surfaces. To generate grid network using two-boundary grid generation, one has to face a difficulty in specifying a proper value for the parameter CK, Eqs.(5.7-5.8). In general, the value for CK is specified by trial and error, however, its variation on a surface can be justified by the use of the shape function, Eq.(5.11), according to the geometry of the grid to be generated. For the six block grids, the grid orthogonality on the outer surface was not enforced, that is CK2=0, to avoid the possibility of grid overlapping. On the other hand, the grid orthogonality on the








75








1z


(a) Block 1


- Figure shows the surfaces 1=1, 1=50, k=5, k=15, j=21


z
Y x
















(b) Block 2 -- Figure shows the surfaces
1=1, 1=50, k=3, k=13, j=17

Figure 5.3 Six block grids -- 3D close-up with cut-off








76


z

Lx








1I






(c) Block 3 -- Figure shows the surfaces
1=1, 1=50, k=1, k=15, k=30, j=10, j=18, j=32


z







x 1 M


(d) Block 4 -- Figure shows the surfaces 1=1, 1=50, k=1, k=15, k=30, j=6, j=14, j=32


Figure 5.3 (continued)








77


z

Yj















(e) Block 5 -- Figure shows the surfaces
1=1, 1=50, k=3, k=13, j=17


















Y



(f) Block 6 -- Figure shows the surfaces
1=1,1=50, k=5, k=15,j=1,j=10


Figure 5.3 (continued)








78

solid surface is desirable. From numerical experiments conducted the following was one of the best to set up the CK, distribution by the shape function for each block.

Block 1 - CPJ=0.5 CSJ=0 CPK=O CSK=O

Block 2 - CPJ=1 CSJ=-0.5 CPK=0.5 CSK=O

Block 3 - CPJ=0.8 CSJ=-0.4 CPK=1 CSK=-0.5

Block 4 - CPJ=0.8 CSJ=-0.4 CPK=1 CSK=-0.5

Block 5 - CPJ=1 CSJ=-0.5 CPK=0.5 CSK=O Block 6 - CPJ=0.5 CSJ=-0.5 CPK=O CSK=O

These distributions are shown in Figure 5.4. For block 1, CK, decreases as j increases and becomes zero when j =jmax. This is reasonable since the j =jmax surface incline about 45 degrees to the solid surface. Similar argument can be addressed for the other blocks. To further improve the smoothness of the grid network, these algebraic grids were then used as initial grids for elliptic grid generation. Smoother grids were obtained for block 1, 2, 5, and 6. However, it failed to generate elliptic grid for either block 3 or 4. This may due to their special geometries such as three boundary surfaces form nearly a flat plane. Figure 5.5 illustrates the difference between algebraic and elliptic grids. Notice that the elliptic grids are smoother and more orthogonal than the algebraic grids.

These six block grids were further modified to include mirror image planes for imposing symmetric boundary conditions in the Navier-Stokes simulations. The resulting grid network was then used as a basic grid in this study and is referred to as the "base grid" hereafter. The dimensions of this base grid are 28x25x60, 40x13x60, 40x35x60, 40x35x60, 40x13x60, and 26x25x60 for blocks 1-6,











79


BLOCK I


RLIC f3 ~


RLC rWV


BLOCK 2











j-jmax kka


BLOCK 4


BLOCK 6


Figure 5.4 Shape function for the six-block grid


k-kmax J-JBx


m nrv 3


I


J- jmax k-mkmax


k-kmax j- jmOX x k - k max
Bi nCK 5
BLOCK 6










6 j-jmax


k-kmaxj- jmx








80


(a) Block 1 - k=13


Figure 5.5 Elliptic and algebraic grids


z
















Block 1, k-13, Algebraic Grid


z

















Block 1, k-13, Elliptic Grid








81


z


I


I-


A F/ 7///// / /


Block 6, k-13, Algebraic Grid





z















Block 6, k-13, Elliptic Grid


(b) Block 6 -- k=13 Figure 5.5 (continued)


Z///, F111 / / /








82

respectively. Figure 5.6 shows some detailed plots to illustrate the grid characteristics. Figures 5.6a-d show that the grid lines are about normal to solid surface whenever possible. Figures 5.6e-f show that some k=constant surfaces are quite skew for block 2, 3, 4, and 5. Figure 5.6g shows the grid distribution at k =1, 12, and 18 wing cross-sections.

In order to investigate the effect of grid resolution on the flow solution, a series of grid refinement studies have been performed. Several locally refined grids (near the wing tip or near the solid surface, etc.) were generated and used to study the possible effect due to the grid characteristics. A set of "refined grids" was also generated by doubling the grid points in the wing core-wise direction. The dimensions of these refined grids are 28x25x60, 79x13x60, 79x35x60, 79x35x60, 79x13x60 and 26x25x60. Figure 5.7 shows the solid surfaces of block 2 and 3 of this refined grid. Note that the grid lines were generated by using a cubic spline so that the solid surfaces would not be distorted.











83


z


j-1 or 40 Block 1 k-13


Block 3 or 4
k-35


'h ~//A'$"' K (I// <. /
/ // /






Block S k13


(a) Block 1, 3, 6


Figure 5.6 Some detail plots for six-block grids


J-40 or 1


(b) Close-up to show the boundary orthogonality
















Block 3


84


z Black 4



















Block 5
J-20


(c) Block 2, 3, 4, 5


Figure 5.6 (continued)


B lock 2
i-20


(d) Close-up to show the boundary orthogonality


Z










85


Y


x















k-6 1-1-35


(e) Block 3, 4 -- k surfaces are quite skew


k-18,13 1 3
8 1ock 2 I

z







(f) Block 2, 5 -- k surfaces are quite skew


Figure 5.6 (continued)











86


Block 4 k-I











Y



Block 3 k-I 1-1-35

(g) Block 3, 4 - k= 1, 12, 18 to show the grid distribution
on the leading and trailing edge


Figure 5.6 (continued)


Block 4 k-18











c . kY Block 3 k 18


Block 4 k-12











Y



Block 3 k-12









87


(a) Base grid -- Block 3, 4 -- 40x35
Block 2, 5 - 40x13


x
Y
z


Refined Grid - 79x35 on wing surface


(b) Refined grid


- Block 3, 4 -- 79x13
Block 2, 5 -- 79x13


Figure 5.7 Base grid and refined grid


x


Base Grid - 48x35 on wing surface


I













CHAPTER VI
SOLUTION METHOD


The multi-block computational method considered is rather simple and

straightforward in application. For the wing-fuselage flow problem, the flow field has been properly divided into six contiguous blocks, and the grid network has also been generated. In the solution process, each block is then treated as an independent flow problem that is solved by the same thin-layer Navier-Stokes code. The interface boundary conditions between blocks are updated with a distance weighted interpolation before the computations march for another time step. A computer code has been developed and investigated for the six-block wing-fuselage aerodynamics in which the unsteady 3-D TVD thin-layer NavierStokes code is made to a subroutine subprogram. Since the multi-block computer code has to deal with several blocks of different configurations and their interfaces, the code has to be modular and flexible. Also, the code has to adapt to the architecture of the host computer system. A six-block wing-fuselage Navier-Stokes code (WF6B) was then written for the available computer system, Cray-2 of NAS system at NASA Ames Research Center. The special feature of the supercomputer, Cray-2, is its huge main memory, however, the computing speed is a little slower than Cray Y-MP. In order to take advantage of the large main memory while not sacrifice the computing speed, WF6B uses part of the


88








89

main memory for storing the intermediate data. The other main memory is then used as active area for executing computations. A set of common arrays is then allocated at the active area and used for the computations of each blocks. A block diagram to show this data arrangement is given in Figure 6.1. A detailed description of WF6B can be found in [43].

In order to illustrate the multi-block solution method a flowchart of WF6B is given in Figure 6.2. As shown in the figure, two do-loops run over each blocks and the number of time steps requested. In the inner loop, for each block, BC is used to set up the boundary conditions and then TLNS is to formulate the governing equations and perform the three inversions. A value of L2 residual for each block is computed by RESID every 10 time steps. After the computations run over each block, INTERBC updates the interface boundary conditions, and then the computation continues to the next time step. Finally, the pressure and the shear stress coefficients are computed by CPP and CSP, respectively. In what follows, the techniques used to calculate the metric coefficients and Jacobians of transformation are described first. Then the boundary conditions imposed on the boundary surfaces of the six-block flow problem are discussed while the interface boundary conditions are given in a separate section. Finally, the equations used to compute the L2 residual and pressure and shear stress coefficients are given.


6.1 Jacobian and Metric Coefficient


The metric coefficients and Jacobians of transformation can be calculated from Eqs.(2.16-2.17), if the covariant metrics, x,, y,, z,, etc, are known. For a second-order approximation, the covariant metrics are calculated at interior grid









90


Block 1

28x25x60 rl(3,j,k,l) ql(j,k,1,6)


Block 2

40x13x60 r2(3,j,k,l) q2(j,k,1,6)


Block 3

40x35x60 r3(3,j,k,l) q3(j,k,l,6)


Block 4

40x35x60 r4(3,j ,k,l) q4(j,k,1,6)


Block 5

40x13x60 r5(3,j,k,l) q5(j,k,1,6)


Block 6

26x25x60 r6(3,j,k,l) q6(j ,k,l,6)


Note: r's are the Cartesian coordinates
q's are the unknown vectors Block data sweep in


Active Area

40x35x60 Common Blocks


/varl/ /vars/


/btrid/ a(3,j,1,5,5) /vis/ turmu(j,k,l)
/var3/ xx(9,1)
etc.
r-Cartesian coordinates q-unknown vector s-RHS vector a-5x5 matrices turmu-eddy viscosity xx-metric coefficients


results sweep out


read in grid data and/or restart file


write out results


Figure 6.1 Data arrangement in WF6B


r(3,j,k,1) s(j,k,l,5)
q(j,k,1,6)


Permanent Disk stores r & q for
each block









91


do N = 1,NMAX


grid data


restart data


do BLOCK = 1,6


BC


TLNS


L E. U


INTERBC


residual


write out


r 9


restart data


pressure coefficients


shear stress coefficients


Figure 6.2 The flowchart of WF6B


-1


1


H-


+ -


if CPP = yes
CPP


if CSP = yes
CSP


read in grid data


compute Jacobians &
generate initial condition or read in restart data --


I I


---








92

points by a central-difference formula and second-order one-sided differencing at the boundaries, e.g.,

X = (x,+1,,-xj.1,)/2 for 1< j
xe = (-3x,k+4xj+j-xj+2,)/2 for j = 1 (6.1)

xf = ( 3xj4-4x._W+xj_.2)/2 for j = jmax

The Jacobian, J, of each grid cell is then computed from Eq.(2.17). The contravariant metric coefficients are defined in Eq.(2.16), and can be calculated accordingly.

The Jacobians and the metric coefficients evaluated by the above finitedifference approach are rather sensitive to grid characteristics, and can produce a large truncation error if the grid network does not have good properties of smoothness and orthogonality. It is known that they can also be approximated from geometric relations such as cell volume and cell-face area. In the discretized sense, the Jacobian defined in Eq.(2.17) is equivalent to one over the cell volume, 1/Vol, and the contravariant metric coefficients defined in Eq.(2.16) are equivalent to one over the cell-face normal. If one defines a computational cell that is centered at the corresponding grid point, with its cell-faces located half-way to the neighboring points, it is then appropriate to replace the Jacobian, J, by the corresponding cell density, 1/Vol. Furthermore, since ,/J, 4,/J, and c./J defined the x, y, z components of the normals (not unit normals) to the local tangent plane of the constant ( surface, one can replace , , by n/Vol, ny/Vol, and n,,/Vol, where ny, nr and nj are the components of the normals to the local constant j planes. Likewise, one can replace t x, '7, 7. and ., Y, . by








93

nh/Vol, nh/Vol, nh/Vol and nb,/Vol, n1y/Vol, n./Vol, respectively.

A general hexahedron, defined by eight arbitrary corner points, can be partitioned into six tetrahedrons, Figure 6.3a. The volume of a tetrahedron denoted by its vertices (ab,c,d) equals one-sixth of the triple product of the three vectors emanating from one of the vertices and ordered according to the righthand directions, Figure 6.3b. For example,

V'e(a,b,c,d) = [rd -(rbd x rd)]/6 (6.2)

= [(xe-x+(b-Yxa-z -xd)(y-yd)(ze-zd)

+ (xa-xd)(ye-yd)(zb-zd)-(Xa-Xd)(yb-yd)(Zc-Zd)

- (x,-xd(y-Yd)s-zd)- -xXy-yd)(z.zd)]/6
where r, = r - ri, and ri and rj are the position vectors. If the eight vertices of the hexahedron are denoted by numerals 1-8 and are associated with the j,k,l triplets as follows:

1 - j,kl 2 - j+1,k,l

3 - j+1,k+1,1 4 - j,k+1,1

5 - j,kl+1 6 - j+1,kl+1

7 - j+1,k+1,1+1 8 - j,k+1,1+1

the volume of the hexahedron is

Vol(j,k,l) = Vet(1,2,3,7) + Vtel(1,2,7,6) + Ve(1,5,6,7)

+ V'l(1,3,4,7) + Vt(1,4,8,7) + Ve(1,5,7,8) (6.3)

The Jacobian of a grid-centered computational cell, which is bounded by the surfaces at the middle of two adjacent grid points, is computed as one over the average volume of all the associated hexahedrons.








94


Vol(j,k,1) = V"(1,2,3,7) + V'(1,2,7,6) + V'(1,5,6,7) +
V*0(1,3,4,7) + Vm(1,4,8,7) + Ve(1,5,7,8)

(a) Partition of a hexahedron




a ~


a--------------b


C


V"(a,b,c,d) = [r -(rd x rd)]/6

(b) Volume computation of a tetrahedrons

Figure 6.3 Volume computation of a hexahedron


5 6


2 3 /


3




Full Text

PAGE 1

NUMERICAL SIMULATION OF TRANSONIC TURBULENT FLOWS OVER A COMPLEX CONFIGURATION SHU-CHENG YANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTL^L FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1990

PAGE 2

ACKNOWLEDGEMENTS The author wishes to express his sincere appreciation to his dissertation advisor. Dr. Chen-Chi Hsu, for his guidance and encouragement throughout this endeavor. Special thanks are expressed to Dr. W. Shyy for his invaluable criticism of the draft of this dissertation and his stimulating discussions. Many thanks are also extended to the other supervisory committee members. Dr. B.M. Leadon, Dr. U.H. Kurzweg and Dr. A.K. Varma for their interest and comments on the research project. . . » ^ i . : • Thanks are also due to Dr. M.H. Chen and Dr. N.H. Shiau for their assistance on their previous work of the TVD scheme and computer codes. Special thanks must go to Mr. J. Ordonez for his tireless proofreading of the manuscript. Also, the author has a debt of gratitude to Dr. N.C. Yao and Dr. C.G. Tu of the Chinese Navy, without whose encouragement and help he would not have had a chance of pursuing the advanced study here at the University of Florida. Finally, the author is indebted to his family for their continuous support. This work was partially supported by a NASA-Ames research grant, NAG2473. Most of the calculations were performed with a Cray-2 of the NAS System at NASA Ames Research Center. ii

PAGE 3

TABLE OF CONTENTS Pa ge ACKNOWLEDGEMENTS ii ABSTRACT v CHAPTERS I INTRODUCTION 1 n GOVERNING EQUATIONS 7 2.1 Navier-Stokes Equations and Equation of State 7 2.2 Non-Dimensionalization 12 2.3 Transformed Navier-Stokes Equations 14 2.4 Thin-Layer Approximation 17 2.5 Turbulence Model I9 m NUMERICAL ALGORITHM ?. 25 3.1 Beam and Warming Scheme 26 3.2 TVD Scheme 32 3.3 Flow Solver 43 IV THE FLOW PROBLEM AND SOLUTION STRATEGY 45 4.1 The Flow Problem and Preliminary Work 45 4.2 Solution Strategy 5I 4.3 Six-Block Grid Topology 54 V GRID GENERATION 62 5.1 Surface Grid Generation 63 5.2 Volume Grid Generation 66 5.3 The Six-Block Grid 74 11

PAGE 4

VI SOLUTION METHOD 88 6.1 Jacobian and Metric Coefficient 89 6.2 Boundary Condition 95 6.3 Interface Boundary Condition 105 6.4 Convergence Criterion 108 6.5 Pressure and Shear Stress Coefficient 109 Vn RESULTS AND DISCUSSION Ill 7.1 The Result from Base Grid 112 7.2 The Result from Refined Grid 115 Vm CONCLUDING REMARKS 130 APPENDICES A JACOBIAN MATRICES 132 B EIGEN-FUNCnONS 134 C CLUSTERING FUNCTIONS 136 REFERENCES 139 BIOGRAPfflCAL SKETCH 144 iv

PAGE 5

•'f Abstract of Dissertation Presented to the Graduate School • of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy NUMERICAL SIMULATION OF TRANSONIC TURBULENT FLOWS OVER A COMPLEX CONFIGURATION By Shu-Cheng Yang December 1990 Chairperson: Chen-Chi Hsu Major Department: Aerospace Engineering, Mechanics and Engineering Science A multi-block computational method is developed for three-dimensional high speed turbulent flows over complex configurations. In this method, the flow domain is divided into several contiguous blocks such that each block is partially bounded by a soUd surface. The chosen grid topology has advantages in computing eddy viscosity as well as in applying a thin-layer Navier-Stokes code. In the solution process, each block is then treated as an independent flow problem governed by the same set of time dependent thin-layer Navier-Stokes equations. The interface boundary conditions between blocks are updated at every time step of the solution algorithm. The application of the method with distance weighted interpolation for updating the interface boundary condition is rather simple and straightforward; however, special measures in grid topology and V

PAGE 6

generation are required. The flow solver employed is a total variation diminishing thin-layer Navier-Stokes code implemented with an algebraic eddy viscosity model of Baldwin-Lomax. In this study, the proposed method is investigated for the simulation of a transonic turbulent flow past a wing-fuselage configuration in which the flow field is properly divided into six blocks. A coarser grid was first used for flow simulation, followed by grid refinement smdies. With the use of the refined grid, the calculated wing surface pressure distribution agrees well with measured data. Numerical results obtained also show that the computed shock location and pressure distribution on the wing can be in good agreement with the experiment data if the block grids used are adaptive to important solution characteristics. It is concluded that the proposed method is indeed a very promising approach to be developed further for complex configuration flow simulation. vi

PAGE 7

CHAPTER I INTRODUCTION For a high speed flight vehicle, numerous complex three-dimensional fluid flow phenomena appear, including the formation and shedding of vortices, the shock-wave boundary-layer interaction and induced flow separation, and the shock on shock interaction. Since the performance of a flight vehicle is strongly affected by these complex flow phenomena, accurate prediction of them and associated aerodynamic forces and moments are critical for advanced design of a modem flight vehicle. At present an accurate wind-tunnel simulation of the complex flow problem in flight conditions can be very difficult, if not impossible. With the advent of supercomputers and the advancement of solution techniques for nonlinear problems, computational fluid dynamics (CFD) is under extensive development and has been apphed to complement the wind-tunnel experiment and field test in the design process of a flight vehicle. The use of CFD in the aerodynamic design of flight vehicles has / progressed rapidly over the last few years. In the early days, CFD was used to support the vaUdity of a design that developed in the wind-tunnel by trial and error. The state-of-the-art of the CFD method has progressed to the point of being regarded as an important design tool for many flow problems [1]. The reduction of wind-timnel occupancy is only a direct benefit of CFD. Detailed 1

PAGE 8

flow information that is unattainable from the wind-tuimel test often can be obtained by numerical simulations. This information about the underlying flows can lead to a better understanding and ultimately to a better design. In fact, the implementation of CFD has fostered a revolution in the design process of flight vehicles. Current CFD methods have only demonstrated an ability to simulate flows about complex geometries with simple physics or about simple geometries with more complex physics. The panel methods are unique in that they provide a capabihty for solving the flow about completely general configuration. Their principal limitation is that they are restricted to simple physics as modeled by the linear equation. The Euler code provides more flow information than the panel method; however, the important viscous effects are excluded in the approach. For more realistic flow simulations, the application of Navier-Stokes equations is required. However, the use of full Navier-Stokes equations is simply too expensive to be practical for today's computers. An alternative is to use the thinlayer Navier-Stokes equations. In general, the governing equations are approximated by a numerical scheme and then solved in either a structured or an unstructured grid network. The latter is quite flexible in gridding but it requires more work in bookkeeping and computer coding. In the structured grid approach, the governing equations are first approximated by either a finitedifference or a finite-volume scheme, and then solved in a boimdary-fitted curvilinear coordinate system. This approach has advantages in boundary condition treatment and computer coding; however, it may have difficulty in

PAGE 9

3 generating good grid system for complex configurations. As evidenced by papers published in journals and presented at conferences, e.g. [2-4], the flow simulations involving simple geometries have progressed considerably; however, the capability for turbulent flow problems involving complex 3-D configurations is still in a development state. Hence, further research and development effort is needed to advance the CFD methods for more complex geometry flow simulations. For flow past a complex configuration, such as a flight vehicle, the generation of an acceptable single structured grid for the entire flow field is rather difficult. It is much easier to divide the flow field into several subregions, termed blocks, and then generate a good grid for each subregion. Hence, it is natural to develop a solution method based on the multi-block grid network. There have been a number of multi-block solution techniques proposed for the Euler simulation of aircraft models and other complex configurations [5-7]. However, the topology of block gridding for viscous turbulent flow simulation is more involved than that for inviscid flows and requires a lot more grid points to resolve the thin viscous layer. In general, a good block grid topology for Euler simulation may not be a proper one for Navier-Stokes sunulation. There are also zonal methods developed to use different equation sets in different regions of the flow field. In general, the flow field is divided into inviscid and viscous zones, and then the Euler equations are appUed in the inviscid zone, while the boundary-layer or Navier-Stokes equations are used in the viscous zone. A zonal approach that was based on the Euler/thin-layer NavierStokes equations was proposed and tested successfully for transonic wing flow [8].

PAGE 10

Another zonal algorithm that applied the Euler and thin-layer Navier-Stokes equations for simulating the flow field of isolated wings was reported in [9]. The flow field was divided into four zones; two inviscid zones with coarse grids and two viscous zones with clustered grids. The computed results were in good agreement with experimental data for the surface pressure, except in the immediate vicinity of the tip and in the shock-induced separated region. The zonal approach also has been used to simulate flow problems for more complex geometry. A transonic viscous flow past an F-16A wing-fuselage configuration has been simulated by a zonal approach [10]. The flow field was divided into as many as 16 zones; in the inner zones adjacent to no-slip surfaces the thin-layer Navier-Stokes equations are solved, while in the outer zones the Euler equations are used. The prediction of wing surface pressures was quite good, except at the leading edge and the aft-shock position. The zonal approach has been shown to be rather efficient. However, the difficulties are to locate properly the interfaces for the inviscid and viscous zones and to patch the solutions that are obtained from different equation sets at these interfaces. Recently, a zonal method with Chimera overset grid scheme was investigated for subsonic turbulent flow about the F-18 fuselage forebody and the combined wing-fuselage [11]. Special coding was required for the overlapped grid region, but the computed results have been shown to be in good agreement with flight-test flow visualization and surface pressure measurements. In this study, a novel multi-block computational method is proposed and explored for the simulation of transonic turbulent flows over complex

PAGE 11

configurations. In this method, the flow domain is divided into several contiguous blocks such that each block is partially boimded by a solid surface. For ease in applying thin-layer approximation and the Baldwin-Lomax turbulence model, the sohd smface is mapped onto an entire boundary plane in computational space. In the solution process, each block is then treated as an independent flow problem governed by the same set of time dependent thin-layer Navier-Stokes equations. The interface boimdary conditions between blocks are updated at every time step of the solution algorithm. The solution algorithm with distance weighted interpolation for updating the interface boundary condition and the computer coding are rather simple and straightforward. The multi-block method offers several distinct advantages over the single block computation. In fact, the use of block grids makes the problem of simulating flow fields about complex geometry more tenable. If a geometric feature is changed, only the related block requires to be modified without completely redoing the basic grid topology. Also, the solution procedures are adapted to the modern computer architecmre; in parallel processing, each block can be solved at different processors; or in case * short of main memory, each block can be solved at a time while the other blocks remain on extended storage. A transonic turbulent flow over a wing-fuselage configuration is investigated for the development of the proposed method. The flow field is properly divided into six contiguous blocks. A coarse base grid was first generated and used for the flow simulation, and then a refined grid was generated from the base grid by halving the grid spacing along the wing in the

PAGE 12

6 chordwise direction. A flow simulation package that includes grid generation, flow solver, and plotting program, was developed and applied to the six-block flow problem. The method was first investigated for a transonic tiu-bulent flow of Mach 0.8 past the wing-fuselage configuration at zero angle of attack, then simulations for different angles of attack, a =4" and a =-3°, were conducted to gain further insight and understanding on the solution algorithm. With the use of the refined grid, the calculated wing surface pressure distribution agrees well with measured data. Nimierical experiments conducted also show that special measures and care are required in generating good grid systems involving excessive distortion between the physical and computational domains; however, the results obtained indicated that the proposed method is a very promising approach for the simulation of complex configuration aerodynamics.

PAGE 13

CHAPTER n GOVERNING EQUATIONS Some theoretical background that required for transonic turbulent flow simulations are described in this chapter. The Navier-Stokes equations and ideal gas equation of state are first discussed. Then, they are non-dimensionlized and transformed to a curvilinear coordinate system for ease in numerical appUcations. The thin-layer Navier-Stokes equations valid for high Reynolds number flows are obtained by neglecting the diffusion terms along the du-ection of the soUd surface. The time-averaged thin-layer Navier-stokes equations are used for turbulent flow simulations and the Baldwin-Lomax model [12] is used as the turbulence closure. 2.1 Navier-Stokes Equations and Equation of State 2.1.1 Navier-Stokes Equations The fundamental equations of fluid dynamics are based on the following universal laws of conservation: (a) conservation of mass; (b) conservation of momentum; (c) conservation of energy. The Navier-Stokes equations that govern most of the compressible flow problems can be obtained by applying those universal laws to a fluid flow [13-14]. 7

PAGE 14

8 Continuity equation . The conservation of mass law applied to a fluid passing through an infinitesimal, fixed control volume without sources or sinks of mass yields dp dpyx dps apw ^ , , — + + + = 0 (2.1) at ax ay az ^ ' where p is the fluid density and u, v, w are the components of fluid velocity in x, y, z direction, respectively. Momentum equations . The conservation of momentum law appUed to the control volume while neglecting body forces reads i£^4--^(pu^+p)+-^+i^ = zx at ax " dy az ax ay az at ax ay az ax ay az a^w apwu apwv a , , ^ 8t„ dr ar +— +— + (pw^ + p) = S.+ .L_S.+ "'^a (2.2) at ax ay az ax ay az where the components of the viscous stress tensor ry for a Newtonian fluid are given by the constitutive equation au 3v aw. ^ au 2 au av aw '"xx = A( — +— +— ) + 2m— = ±^(2^-^.^) ax ay az ax 3 ^ ax ay az ^ " 'ay sz 5x ' ' ay 3 '^^ ay az ax ' _ , 3w 3u av aw 2 aw au av , au av r

PAGE 15

av aw , aw au , IS is a Here, p is the thermodynamic pressure and fi is the dynamic viscosity. The fluid assumed to be isotropic and satisfying the Stokes condition, A = -2/3 n, which good approximation for flow problems without relaxation process. For air, the viscosity ^ depends mainly upon the temperature and can be estimated by using an interpolation formula based on Sutherland's theory of viscosity, i.e., , T,L5 T„ + S M = f^.(^r g (2.2b) where denotes the viscosity at the reference temperature T„, and S is a constant which assimies the value 198.6"R [15]. ' ' Energy equation. The conservation of energy law applied to the control volume and assuming no external heat addition results -^[(E+p)u]+ i-[(E+p)v]+ 4-[(E+p)w] = ^ (2.3) 01 ax ay az ax ay az where ! . aT ^, = Ur„ + Vr^ + Wr„+k— ^ . aT '' . ^ = Ury^ + VTyy + Wry, +k— " , (2.3a) aT = Ur„+Vr^ + Wr„ +k— ,

PAGE 16

10 Here, E is the total energy per unit volume, T is the temperature and k is the thermal conductivity. The heat loss by conduction qj is related to the temperature gradient by Fourier's law of conduction aT = -k— (2.4) The equations described above, Eqs.(2. 1-2.3), constitute a complete set of the time dependent Navier-Stokes equations that is valid for a class of compressible flow problems without heat addition and chemical reaction (or relaxation). The Navier-Stokes equations includes one continuity equation, three momentum equations and one energy equation, but there are seven unknowns, namely, p, u, v, w, p, E, and T. The relationships between the thermodynamic variables (p, p, T, E) are established by the equation of state of the working fluid. 2.1.2 Equation of State In this study, air is the assumed working fluid and the perfect gas equation of state is employed. For a perfect gas, the thermal and caloric equations of state are p = pRT, e = QT ( Q = constant ) (2.5) respectively. Here, R is the universal gas constant, e is the internal energy. Some other useful relations are h = CpT, 7 = Cp/Q a^ = -yRT (2.6) where h is the enthalpy, Cp and Q are the specific heat at constant pressure and volume, respectively, 7 is the ratio of specific heats, and a is the local speed of sound.

PAGE 17

If only internal energy e and kinetic energy are considered significant, the total energy E can be written as E = p[e + (u^+v^+w^)/2] Consequently, one obtains two specific equations relating p, T to the dependent variables p, u, v, w, and E, i.e., 1 c 1 (2.7) If expresses heat conduction qj in terms of internal energy through the use of caloric equation of state, one can rewrite Eq.(2.4) as follows. , aT k ae yu ae Here, the Prandtl number Pr is defined as Pr = C^ti/k, for air Pr = 0.72, approximately, and the internal energy is given by e = ^-yCu'+v^+w^) (2.8a) and the /9's in Eq.(2.3a) can also be rewritten as iii de = Ur„ + Vr„ + Wr„ + " Pr ax ae = Ur^ + Vr^+Wr^ + — — (2.9) = Ur„ + Vr„ + Wr + Pr ay 7/i ae " ° Pr az

PAGE 18

2.1,3 Vector Form 12 The Navier-Stokes equations, Eqs.(2. 1-2.3), form a system of five equations and can be cast in vector form as aq aE aF aG aEL aF, aG, at ax ay az ax ay az (2.10) where q is an unknown vector, E, F, G are the flux vectors and E^ F^ G^ are the viscous flux vectors. The explicit expressions for these vectors are q = p \ 1 pU pu +p ' pV E = < pUV , F = pW puw [e J L(E+p)uJ v pV pVU pv^+p pVW G = pW pWU pWV 3W^ + pW^ + P Ue+p)wJ 0 xz I9x J F = '0 " ro (2.10a) ' T zy [fir) in which ry are given by Eq.(2.2a), ^'s are defined in Eq.(2.3a), and p, T are related to q by Eq.(2.7). 2.2 Non-dimensionalization To reduce the number of parameters that characterize the same class of problems, a dimensional analysis is employed. The non-dimensional quantities are obtained by choosing characteristic scales in the following manner:

PAGE 19

where the non-dimensional variables are denoted by an asterisk *, the freestream conditions by «, a„ is the freestream speed of sound, and L is a characteristic length scale. " ' By substituting the above relations into the governing equations, Eq.(2.10), the non-dimensional equations in vector form are aq* , aE* aF* aG* „ dE* dF* aG ' -T+ — s+ — s+ — r = Re X — ^ + — ^ + — /) (2.12) at ax ay az ^ ax ay az ^ . ; ^ ^ where the Reynolds number, Re, is defined as v. Re = and the relationships from the perfect gas equation of state, Eq,(2.7) become p* = (^-l)[E*Ap\u*^+v-2+w-2)]= (-y-l)pV C V • E* 1 ' " 1^ = 7(7-l)[— -^(u'^+v'^+w'^)] ''^^ ^ • (2.13) The unknown vector q* and flux vectors E*, F*, G*, E/, F/ and G* have the same form as in Eq.(2.10a), except that all the quantities are now non-dimensional. In summary, the dimensionless equations govern a class of flow problem. Only two parameters, namely, the Reynolds number, Re, and the Prandtl number, Pr, appear in the equations, the Mach number is dropped out due to the use of a, as characteristic velocity scale. In the following, the asterisks * are removed from the non-dimensional equations except where noted.

PAGE 20

14 2.3 Transformed Navier-Stokes Equations In the finite-difference approach, to enhance numerical accvu-acy and efficiency in applying boundary conditions, the governing equations are transformed from the Cartesian coordinates to a boundary-fitted curvilinear coordinate system. Consequently, a general computational code can be developed such that the computation is performed in an uniformly-discretized computational domain [16]. If one introduces a generalized coordinate system e = ?( x,y,z,t ) n = r,( x,y,z,t ) (2.14) r = r( x,y,z,t ) r = t and applies the chain rule of partial differentiation, the non-dimensional governing equation, Eq.(2.12), becomes q. + QeC. + q.-^. + q^ r . + (C^^ + v^fi, + ^^,) + (CyF^ + VyF, + TyF,) + (e,G^ + r,fi^ + r,G,) = Re-^[(^ A, + r,^^ + rxE^) + + VyF^ + CyF^) + (C,G^ + vfi^ + r,G^)] (2. 15) The metric coefficients appearing in these equations can be determined by the matrix relation as follow i^ ' ^? ^f) Vy »?, y? y. Vc yr k '<} . T.f Z„ Zf. . o' 0" \ . or

PAGE 21

15 f/x 'y X ''y ''z fx f y f z . Xj. y? y. yf Zf -1 = J (yfZf-yfZe) /fZ^-x^z -(x^y^-x,y^) (XfZ^-x,Zj) x^y^-x^y^ J yez.-y.Zf '^x Cy ^0 ' »»» '?x '?y '?Z "yr ' I Tx Ty fz J Xj X^ Xj. y? y. yf z. where J is Jacobian of the transformation and J'\ representing the volume of the grid cell in physical space, is given by _ a(x,y,z) a(e,i7.f) = Xe(y.ZryfZ,)-x,(yeZf-yfZ^)+x^(y^z^-y^Zj) . . (2.17) In order to put the transformed equations, Eq.(2,15), into conservation-law form, the equation is first divided by the Jacobian and then rearranged to give (qJ/J + (l.qf + ?xEe + eyF, + e.G,)/J + („,q^ + ^ JE, + + ^^g,)/j + (r.qf + rxE, + fyF, + f ,G,)/J = Re-^[(e + e/,, + ?,G,,) and then rewritten into the following conservation law form (q/J)r + (a.q+e,E+eyF+C,G)/J)^ + (M+»7^+r,yF+r,,G)/J)„ + ((r.q+fxE+fyF+r,G)/J), = Re-^[((eA+^yF,+ ?,GJ/J), + . + ((.7A+'7yF,+ r,,G,)/J), + ((rJE,+ ryF,+ f,G,)/J),] , (2.18) in which the follov^ing metric identities [10] have been apphed.

PAGE 22

16 (2.19) — (7^)+— (f) = 0 a? ^ J ^ at; M ar M ^ The transformed equation, Eq,(2.18), has the same form as the original equation, Eq.(2. 12), and can be written as v V ' .v ar ae a»7 ar a$ a?; ac (2.20) where q = q/J, E = (^,q+|^+C/+|,G)/J, F = (,7.q+»7^ + r7yF+r,,G)/J, A G = (r,q+fJE+fyF+f,G)/J, (2.20a) A = (^A+?,Ft+?zG.)/J, A = („A+'/yFT+'7zG,)/J, = (rA+ryF,+f,G,)/j The explicit expressions for the transformed unknown vector, flux vectors and viscous flux vectors are T-i q = J p pV pVi E F = J-^ E = J -1 /juU+^j) pvU+^p pwU+^j) L(E+p)U-^j5j puV+r/J) pvV+»7yP pwV+r/J) (E+p)V-„j)J 0 Cx'-xy+Vyy + ^z''xy ^x'-xt+Vp'^^^^^ ^x^x+?A+^A G = J-^ puW+fj) pvW+Typ pwW+fj) [(E+p)w-rj)J F, = J-^ 0 '?x'"xx+Vyx+''z'"zx '7x^xy+Vyy + ''z'-zy (2.20b)

PAGE 23

0 where U, V and W are contravariant velocity components given by V = f7,+T7,u+r;yV+»7,w (2.20c) w = r,+f,u+fyV+r^w ' It is noted that the velocity gradients in Ty, Eq. (2.2a), and the temperature gradient in /9's, Eq.(2.3a), have to be transformed into computational space also, e.g., • -k ax ay az ' ax ^r^ 3u au au av^ av av ^ a? a»; ar . aw ^ _ aw aw ^ au au au + (— — »7z+ — rz)]+24(— -Ffx)] ae a»7 ar a^ at? ar 2.4 Thin-Layer Approximation The numerical simulation of complete Navier-Stokes equations, Eq.(2.20), is in general rather time consuming and needs a huge amount of computer storage. For high Reynolds number flows, however, the viscous effects are often confined to a thin layer near no-slip boundaries, called the boundary-layer, outside of which the flow can be considered to be mostly inviscid. With this in mind, it is reasonable to assume that the diffusion processes along the body

PAGE 24

18 surface are negligibly small in comparison with those in the direction normal to the body surface. As a consequence of this assumption, all viscous terms containing derivatives parallel to the body surface are dropped since they are substantially smaller than viscous terms containing derivatives normal to the surface, while all other terms in the momentum equations are retained. One of the advantages of retaining the terms which are normally neglected in boundarylayer theory is that separated and reverse flow regions can also be computed. Also, flows which contain a large normal pressure gradient can be readily computed. Another justification of the use of a thin-layer approximation can be given from a computational viewpoint. In practice, a substantial fraction of the available computer storage and time is expended in resolving the cross-stream gradients in the boundary layer since a highly stretched grid is required, and the resolution along the body is treated as what is needed in inviscid flow. As a result, even though the full derivatives are retained in the equations, the gradients along the body are not resolved in an adequate manner. Hence, one can drop those terms that are not being adequately resolved, provided that they are reasonably small. Since the thin-layer model requires a boundary-layer type of coordinate system, it is assumed that the sohd body surface is mapped onto an entire C»7plane in the transformed space, then the thin-layer approximation requires that the e and f; derivatives in the viscous terms be dropped. Upon simplifying the complete Navier-Stokes equations, Eq.(2.20), the system of governing equations becomes i

PAGE 25

aq , aE , aF , aG _ aS dr where ae dri = Re-\— ) ac (2.21) S = J-^ 0 * miUj. + m25-, miVf + mzfy I mim3 + m2(r^u + ryV+r,w) J (2.21a) ^2 = M(r,Uf + fyVf + r,w^)/3 (2.21b) mj = (u2+v2+w2)^/2 + PT\y-ir\r), Note that the notation has been removed from the above equations for simplicity and the heat conduction term are replaced by jfi ae _ /i aT '^""Pr ar ~ "(7-l)Pr "aT ^^'^^^ The thin-layer Navier-Stokes equations, Eq.(2.21), in conjunction with the equation of state, Eq.(2.13), constitute the required governing equations for this study. They are valid for both laminar and turbulent flow problems. However, due to the widely disparate scales of turbulence, it is impractical to apply the equations directly for turbulent flow simulation. An alternative is to use the time-averaged equations coupled with a turbulence model. 2.5 Turbulence Model The time-averaged equation is obtained by decomposing the dependent variables into time mean and fluctuating components and then time-averaging the

PAGE 26

20 entire equation. In the averaging process, new unknowns will be introduced and must be treated by some approximating techniques so as to close the problem. In this study, the mass-weighted time-averaged equations are used for turbulence flow simulation with the eddy viscosity that appears in the averaged equations provided by the Baldwin-Lomax two-layer algebraic model. In this approch, the mass-weighted time-averaged thin-layer Navier-Stokes equations will have the same form as the original equations, Eq.(2.21), while all the dependent variables are time-averaged. The effective coefficient of viscosity n is comprised of two parts, the laminar viscosity /i, and the tiu-bulent eddy viscosity /i,. The laminar viscosity A*, is obtained from Sutherland's theory of viscosity, Eq.(2.2b), and the eddy viscosity ^l^ is provided by a turbulence model [16]. Similarly, the effective thermal conductivity k is also comprised of two parts, k, and k,. The heat conduction is expressed in dimensional quantities as The heat conduction term in Eq.(2.22) is then replaced by its dimensionless counterpart Here, the turbulent Prandtl number Pr, is defined as Pr, = ^x^c^/k^, and is about 0.9 for air [10]. In this approach, the eddy viscosity /i, is the only quantity that has to be modeled. In what follows, the theoretical background about the algebraic model will be discussed first, and then the Baldwin Lomax model is described.

PAGE 27

21 2.5.1 Two layer algebraic model In analogy with the coefficient of viscosity in Stoke's law for laminar flow, Boussinesq [17] introduced a mixing coefficient, /i„ or referred to as eddyviscosity, for the Reynolds-stress that appears in the time-averaged equation, au ^ = /^.-T— (2.24) ay where u is the streamwise velocity, y is the spatial coordinate normal to the solid siuface, and 3u/ay is assumed to be the only significant component of strain rate. With Prandtl's mixing length hypothesis, the coefficient n^, is then related to mean velocity field by 'vdu dy . . . :r (2.25) where 1 is the mixing length. ' The algebraic model is then to establish a relation between the mixing length and the characteristic length of the respective flow. For simple wall shear flows, the mixing length in the law of the wall region is proportional to the normal distance from the wall, i.e., 1 = ky, where k = 0.41 is a universal constant determined experimentally. In the viscous sublayer, however, the mixing length theory is not valid. A corrected form of the mixing length valid down to the wall was deduced by van Driest. He introduced a damping factor and wrote 1 = ky[l-exp(-yVA^)] . , (2.26)

PAGE 28

22 where y"*" = ypujn is the law of the wall coordinate, A"^ denotes an empirical constant which equals 26 for flat-plate flow, and u^, the friction velocity, is defined as = {rjpf', where r„ is the shear stress on the wall. In the law of the wake region, the Clauser model is commonly used in conjunction with a factor to account for the intermittent effect. The eddy viscosity is then written as where the value is obtamed experimentally and assumed to be 0.0168, U is the streamvdse velocity on the edge of the boundary layer, 5* is the boundary layer thickness, 7 is the intermittency factor deduced from a curve fit of measured data and given by Reynolds and Cebeci [18] as where s is the boimdary layer thickness. Based on the above arguments, Cebeci and Smith [19] proposed a well known two-layer algebraic model as follows. In the iimer-layer: (2.27) 7 = [l + 5.5(y/5)'] (2.27a) (2.28) 1 = ky[l-exp(-yVA*)] (228a) In the outer-layer: A*, = k,U5*7 (2.28b) 7 = [l+5.5(yA)«] (2.28c)

PAGE 29

This model has been shown to be in excellent agreement with experimental ' ; , observations for simple shear flows. However, the model requires the boundary layer thickness and the edge velocity which constitute a practical disadvantage in the implementation of the model. In an attempt to overcome this difficulty, Baldwin and Lomax [12] modified the outer layer formulation to eliminate the necessity of finding the boundary layer thickness. In addition, they replaced the velocity gradient by the vorticity, which is invariant under coordinate transformation, so that the model can be easily applied to 3-D problems. The Baldwin-Lomax model is described next, 2.5.2 Baldwin-Lomax Turbulence model The Baldwin-Lomax model is similar to the Cebeci-Smith two-layer algebraic turbulence model. The eddy viscosity in a profile is calculated for each inner and outer layer, and then matching both layer as follows [12]: where y, is the smallest value of y at which values the inner and outer formulas are equal. The inner layer follows the Prandtl-Van Driest formulation where y is the normal distance to the surface and |w| is the magnimde of the vorticity (2.29) (/^.)in»er = pl'kl with 1 = ky[l-exp(-y VA*)] (2.30) (2.31)

PAGE 30

24 = y(py,^J^/f^ is the law of the wall coordinate, A"^ = 26, and /x^ are the local density, shear stress, and laminar viscosity evaluated at wall, respectively. The eddy viscosity for the outer layer is given by (M,)outer = pk,QpIVakeFkleb(y) ' ' ' • ' " (2.32) where = 0.0168 is the Clauser constant and = 1.6. * . H ' , The parameter F^^^^ is given by Fwak. = mini J^rf /p where = 0.25 (232a) and Udif = (u'+v^+w^)^^ (u^+v^+w^)^^ . (2.32b) The quantities y^ and are determined from the function F(y) = yk|[l-exp(-yVA*)] (2.32c) where is the maximum value of F(y) and y^ is the value of y at which F^ occurs. The Klebanoff intermittency factor Fu,,,(y) is given by Fkid.(y) = [l + 5.5(Cu^/y^']-» and C^.^ = 0.3 (2.32d) The Baldwin-Lomax model has been investigated and apphed to a variety of 2-D and 3-D flow calculations with satisfactory results [10]. The model is capable of handling the compressible flow problems with modest flow separation [20]. The implementation of algebraic models requires only a small amount of computer time and storage; this is particularly important in 3-D computations.

PAGE 31

CHAPTER in NUMERICAL ALGORITHM In numerical simulations, the transformed governing equations described in the previous chapter are approximated by finite-difference equations and then solved in a equally spaced computational space. The basic computer code used in this study is a TVD thin-layer Navier-Stokes code which is based on an original version of thin-layer Navier-Stokes code, ARC3D of NASA Ames Research Center. The original ARC3D was based on the Beam and Warming scheme [21], and has been improved for its efficiency, accuracy and convergence [22]; for instance, there are options for using different numerical dissipations for better stability criteria. Alternatively, to improve the robustness of the code, a symmetric Total Variation Diminishmg (TVD) scheme was implemented into the original ARC3D. This modified ARC3D has shown to be more efficient and robust than the original one [23-26], and hence is adopted in this study. _ The BeamWarming scheme and the TVD scheme were constructed by totally different philosophies; however, in applications they can be viewed as the same scheme with different numerical dissipation terms [27]. Both schemes will be described and discussed in the following. ' 25

PAGE 32

26 3.1 Beam and Warming Scheme The BeamWanning scheme is an non-iterative, implicit approximate factorization finite-difference scheme. The implicit operator is factorized (spatially split) to obtain an ADI-like ( Alternating Direction ImpUcit ) algorithm such that the three-dimensional inversion process is reduced to three onedimensional processes. The basic algorithm is second-order accurate in both time and space [28]. The implementation of this scheme to the system of governing equations, Eq.(2.21), is described as follows. 3.1.1 Temporal-Differencing To approximate the time derivative of Eq.(2.21), the well-known secondorder Trapezoidal Rule or Crank-Nicolson Formula is employed to give ^ = + 3-7^) + 0[(Ar)2] where Aq-=q-*V (3.1) If one rewrites Eq.(2.21) as and then substitutes it into Eq.(3.1), the following equation results 2 " 3{ 3, ac ^ 3f In the above equation, the flux vectors E, F, G, and viscous flux vector S at the advanced n+ 1 time step are in general nonlinear functions of the unknown

PAGE 33

; if 27 4 vector q. They are linearlized while maintaining second-order time accuracy as follows. The flux vectors are linearlized by using a Taylor series expansion about q" as follows E"+^ = E"+A"Aq" + 0[(Aq)^] where pn+i ^ pn^gn^qn ^ o[(Aq)^] wherc G°** = G"+C"Aq" + 0[(Aq)^] where The above approximation is second-order in time since Aq" = q"-'»-q" = 0(Ar) and 0[(Aqf] = 0[(Arf] The viscous flux vector S is linearized by using the procedure suggested by Steger [29]. The coefficient of viscosity ^ and thermal conductivity k are assumed to be locally independent of Jq ( recall that Jq is the unknown vector in physical domain ), so that elements of S can be expressed in the general form Sj = J'^aj -— i where a, is independent of Jq. Each element is then linearized in the following manner by Taylor series expansion S,-^ = S,Vj-l{a,"-^[E(-^r]}jAq/ + 0[(Ar)^] Written in vector form, the viscous flux vector becomes gn+i ^ s"+J-^M"JAq" + 0[(Ar)2] (3,4) A = B = C = 3E aq dF aq aG aq (3.3)

PAGE 34

28 The matrices A, B, C and M are so called Jacobian matrices and the elements of these matrices are given in Appendix A. After substituting all the linearized terms, Eqs.(3.33.4) into Eq.(3.2), one has 2 ^ 3^ dr, af ar ^ ^ , ^ . r ... ^ = .^,(iE^_aF^_aG_j^^.,^ , ^ v a^ ar; ar af ^ ^ ^ ^ where I is the identity matrix, and Aq" is the delta form unknown to be solved for. Note also that a first-order time-accurate Euler implicit scheme can be obtamed by simply replacing Ar/2 on the left-hand side of Eq.(3.5) by Ar [28]. 3.1.2 Approximate Factorization Solving Eq.(3.5) directly is very inefficient since it involves the inversion of a very large system resulted from three dimensions. Beam and Warming appUed an approximate factorization procedure to reduce the inversion to three onedimensional inversions. They rewrote Eq.(3.5) in the factored form Since the error introduced by this factorization procedure has the leading thirdorder term Ar^._aA_ _aB_^ aB_ _aC aC aA ag" 4 ac dr, dr, ar ar ac the second-order accuracy of the difference approximation is not affected.

PAGE 35

3.1.3 Spatial differencing The spatial derivatives appearing on the left-hand side of Eq.(3.6) are approximated by second-order finite-difference formula, e.g. hence the resulting scheme possesses a block-tridiagonal structure that can be solved efficiently. On the right-hand side, the same second-order differencing is employed for the viscous term, while a fourth-order finite differencing approximation, e.g. is used for the convection terms to keep the convection truncation error from exceeding the magnitude of that of the viscous term, ^ ; After applying spatial differencing into Eq.(3.6), the spatial derivatives are replaced by difference operators, and the finite-difference form of Eq.(2.21) is then . .. " , [1+ y-6,A]"[I+ y^5,B]"[I+ ^(5^C-Re-^S,(J-^MJ))]"Aq= -ArC^^E+sZ+SfG-Re-^S)" (3.7) This is the basic Beam-Warming algorithm. The implicit operator has been decoupled to produce an ADI-like scheme, and due to the second-order central-difference operators employed, the algorithm produces block tridiagonal systems for each spatial direction.

PAGE 36

30 3.1.4 Numerical Dissipation The BeamWanning scheme requires the use of numerical dissipations to damp out the spurious oscillations that occurs when dealing with discontinuities. Several different numerical dissipations have been investigated and discussed in [30]. Those employed in the original ARC3D code were proposed by Beam and Warming [28], and Steger [29] and are described as follows: The second-order implicit dissipation operators = -«iJ"\a^J (3.8) are added to the imphcit operators on the left-hand side of Eq.(3.7) in the |, rj and f directions, respectively, and on the right-hand side the following explicit fourth-order dissipation term is added: De = -^EJ-'[(VA)'+(V/+(VfA,)']JAq(3.9) in which symbols a and v are the forward and backward difference operators, e.g. ^ = + = -J^i^rki) + 0[A|] and €j, eg are constant coefficients used to control the amount of numerical dissipations. It was suggested to set €E=At and ei=2€E so that as the time step increases, the numerical dissipations added relative to the spatial derivatives of convection and diffusion remains constant.

PAGE 37

After inserting numerical dissipations into Eq.(3.7), the finite-difference representation of Eq.(2.21) becomes [1+ |l-5,A+D,r[I+ ^6^B + DJ"[I+ y^(5,C-Re-^5,(J-^MJ)) + D,]"Aq" = -Ar(SjE+6^F+5j.G Re"^5j.S)"+DE (3.10) 3.1.5 Solution procedure For simplicity, the following notation is introduced h = [I+y-5,A+D,r L„ = [1+ ^«„B + D 1" Lf = [1+ y^(5fC-Re-^5f(J-^MJ)) + Df]"Aq" R" = -Ar(5^E+5,,F+5^G-Re-^5^S)"+DE where the L's are operators, then the Beam-Warming scheme, Eq.(3.10), can be written as L,L,L,Aq= R(3,11) If one further defines intermediate solutions as Aq* = Lj-Aq" Aq = L,Aq the unknown vector q at the n+1 tune level can be computed as follows. LfAq= R» L. • /q = Aq Vq" = Aq* (3 12)

PAGE 38

q"+i = q-+Aq" The Beam-Warming scheme has been shown to be accurate and efficient for simple flow problems [10], however, for more complex flow cases it can be very sensitive to grid resolutions [25]. . : 3.2 TVD Scheme It is known that finite-difference schemes may produce oscillations when appUed across a discontinuity. For instance, the second-order Lax-Wendroff scheme u;*i = -K uj^i-u? )-V2uil-u){ uj"^i-2uj"+u/!i ), u = aAt/Ax (3.13) for a linear hyperbolic equation au + a =0, a > 0. at ax can produce spurious oscillations near the discontinuity. The Total Variation Diminishing (TVD) scheme is one of the methods designed to eliminate those undesirable oscillations. The notion of TVD was first introduced by Harten [31]. He derived a set of sufficient conditions from which the constructed schemes will not generate oscillations across discontinuities. Since the theoretical background of TVD schemes is rather complex and involved, only relevant basic definitions are given here. The construction of a symmetric TVD scheme and its application to thin-layer Navier-Stokes equations is discussed in the subsequent sections.

PAGE 39

A numerical scheme, either explicit or implicit, can be written in the following operator form [32] L u°** = R u" (3.14) where u""^* is the mesh function at the advanced time step to be solved, u" is the known mesh function at current time step, and L and R are some finitedifference operators designed for the particular scheme, e.g., for the second-order Lax-Wendroff scheme given in Eq.(3.13) L = 1 R = 1 -.^Vm-MMWh ' V:l Here and in the following, Aj+j^ stands for the central difference operator, i.e., Aj^.j4U = Uj+i Uj. The total variation of a mesh function u° is defined to be =iJL I^Vi-u"! =jL ''^ > ' ' (3.15) Here, the notation « is used to emphasize that the summation is taken over the entire domain. If a numerical scheme applies to a mesh function u, the total variation of the mesh function is non-increasing as a function of time, i.e. TV(u"*MTV(u»-^M (3.16) That is, if the numerical scheme shown in Eq.(3.13) satisfies the sufficient conditions, Eq.(3.16), the scheme is TVD.

PAGE 40

34 It has been proved that a TVD scheme is monotonicity preserving [32], that is, for any monotone mesh function u that increasing (or decreasing) monotonously with time, w = Qu is also monotone, where Q is a finite-difference operator. Since a monotonicity preserving scheme will not generate spurious oscillations, and neither will the TVD scheme. In what follows, the construction of the symmetric TVD scheme for a onedimension hyperbolic equation will be described first, then the scheme is applied to the Euler equations, and finally it is extended to the thin-layer Navier-Stokes equations. 3.2.1 Symmetric TVD Scheme The symmetric TVD scheme was developed by Yee [27] and it is a generalization of Roe [33], Davis [34] and Sweb/s [35] TVD Lax-Wendroff scheme. To illustrate the construction of this scheme, a one-dhnension linear hyperbolic equation of the form 8u du + a — at ax + a— = 0 (3.17) is considered. Sweby has shown that the second-order Lax-Wendroff scheme for this equation can be constructed by adding an anti-diffusive flux to the first-order upwind scheme [35], i.e., Uj-** = i.A..j,u" y2i.(l-..)Aj.j4Aj^^,U» , if a > 0. uj*» = ^Aj^^u" Vi«.(l + ..)A.^^Aj.j,u, if a < 0. (3.18)

PAGE 41

where v = aAt/Ax is the Courant number. Then he introduced a limiting function p and rewrote the above scheme as i^*i=uJ-.{l + y2(l-.)[p(r/)/r/-p(rV,)]}Aj.,^^ , if a > 0 ^^^^^ 1^*1= uJ-u{l-y2(l + ..)[p(rj;i)-p(r-p/rj]}Aj,^u, if a < 0 where p is chosen to be a function of the ratio of consecutive cell gradients rj; pj = P(rj). = ^-mu/Aj+mU, rj' = Aj+^^u/Aj.^^u, with the restrictions 0

I 0 if a < ( KYr-) ^^-^ (3.20a) \'^(1 + u)\p(t: )-1] if a < 0 To avoid determining the upstream direction, Davis then modified the dissipation term as follows, = M(i-IH)[i-P(V)]/2 K-(ij-) = IH(l-|H)[l-p(ri-)]/2 . (3.20b) Shortly after that. Roe [33] further generalized Davis' scheme, Eq.(3.20), and suggested a new limiting function Q. Roe's scheme takes the form

PAGE 42

36 ** ^ V V2|H(i-kl)(i-Qj.H)w" + ^2|H(i-M)(i-Qj.h)Vv4u" (3-21) where the limiter Q depends on three consecutive gradients given by Qj+M = Q(Aj.WAj+MU , Aj+ii^u/Aj+i^u) = Q(r*j , r-_j^i) (3.21a) and 0< Qj,H <2/(l-|H) 0 < Qj+^Aj+i < 2/ 1 1/1 (3.21b) 0 < Q^.Jr^j. < 2/IH to satisfy TVD sufficient conditions, Eq.(3.16). There are a variety of limiter Q satisfing the above restrictions, such as Q(r',r'^) = max[0, min(2r",l), min(r",2)] + max[0, min(2r*,l), min(r+,2)] + 1 (3.21c) Roe's scheme was later modified and extended to a one-parameter family TVD scheme by Yee [27]. She first rewrote Roe'scheme, Eq.(3.21), into the form uy*V Ae(hP;^^-h-^) = uf + A(l-e)(hj",^-h;^) (3.22) hj+H = [a(Uj+i + up-{Aa'Qj^,4+ |a|(l-Qj^H)}Aj+Mu]/2, a = At/Ax (3.23) and then simplified the numerical flux to hj+M = [a(Uj+i + Uj)-|a|(l-Qj^^)Aj^^,u]/2 (3.24) Note that the term Aa^Qj^.^^Aj^.^u/2 was taken out from Eq.(3.23), but it makes no difference as long as the limiter function Q is chosen such that the scheme satisfies the TVD sufficient conditions.

PAGE 43

37 For nonlinear hyperbolic conservation law of the form au 3f(u) the numerical flux, Eq.(3.24), was modified to hj+vi = [(fj+i+fjH(Vv4)(l-Qj+w)Aj+^u]/2 (3.26) where ^=f(Uj), Aj^.^4U=Uj+i-Uj, a(u) = df(u)/du is characteristic speed and la(u,) .VhU = 0 (^-^ea) The function *(aj+j4) that is sometimes referred to as the coefficient of numerical viscosity, is defined as :;: *<^) = {!^+z^)/2< ;|z| J (3.26b) to avoid the entropy violating problem [31], where e is a parameter that is described in Ref.[36]. In order for the above scheme, Eqs.(3.22, 3.26), to be TVD, that is, satisfying the TVD constraints, Eq.(3.16), the following conditions obtained by Yee [37], must be satisfied 0 < Qj+H < 2 :^ 0 < Qi^Jx\ < 2/[.(l-e)|aj.J] 2 0 < Qj^^/r^^i < 2/[Kl-e)|aj^m|]-2 ^. (3.26c) Haj^J < 2/[3(l-e)] If one selects the following Umiter that satisfies the above constraints Qj+^ = minmod(Aj.j4U , a^+^u) + minmod(Aj+^u , Aj+j^u) Aj^^^u (3.27) then fi-om Eqs.(3.26-3.27) Yee's numerical flux can be expressed as , ^i^^ = [fj + fj+i + (3.28) fl^j+ii = *(VH)(gj+gj+i) 2A,^i,u ^ (3.28a)

PAGE 44

38 where the minmod function niinmod(x,y) is defined as minmod(x,y) = x , x < y ' 0 , if X y have opposite sign ly, X > y and g. = minmod(Aj+v4U , a,.^u) (3.28b) The above scheme, Eqs.(3.22, 3.28), form a class of TVD schemes with the numerical dissipation terms centered and hence are often referred to as symmetric TVD schemes. It is noted here that the scheme is TVD for onedimensional nonlinear sccder hyperbolic conservation laws, and is also TVD for constant coefficient hyperbolic systems. However, in apphcations, it is formally extended to more complex equations, such as the Euler equations, the thin-layer Navier-Stokes equations, and evaluated by numerical experiments. 3.2.3 Extension to Euler equations Extension of the scaler TVD scheme to systems of conservation laws is not unique. Here, one follows the method proposed by Yee [38,39]. The method is accomplished by defining at each point a "local" system of characteristic fields and then applying the scheme to each of the scalar characteristic equations and hence is referred to as the "local characteristic approach". Consider the inviscid part of the thin-layer Navier-Stokes equations, Eq.(2.21), dq aE aF 3G = 0 at a^ art a^ (3.29)

PAGE 45

Recall that the Jacobian matrices of E, F, G are A, B, and C, respectively. Denote the eigenvalues of A, B, and C as (af.a-fja'j.aj.a^), (a^,a^,a?^,a'^,a^), (apaji^apa^a^), and Rp R^, Rj. as the eigen-matrices whose columns are eigenvectors of A, B, C, respectively, and R"*^ R"*,, R**j. as the inverses of the corresponding eigenmatrices. The explicit form of these eigenvalues and matrices are given in Appendix B. In the ^-direction, let Qj+h^ be some symmetric average of and Qj+i^, e.g. , or Roe's average qj+v4w = (^^jW*Jj,k4+^Viw*Jj+iw)/(^\ia+^Viw) (3-30) and denote R^^.^^, R'\+^ as the quantities a'^, R^, R"*^ evaluated at qj+M,k4The difference of the characteristic variables in the local f -direction are defined as = RVv4(JjMW
PAGE 46

40 where A = At since A^=Ar7=Ar = l and the numerical flux functions obtained from Eq.(3.28) are given by ^ [Pj,k4''"^j,li+14''"^k+V4*k+M]/2 (3.33) the expression of the elements of ^j+j^ are obtained from Eqs.(3.28a-3.28b) and are given as ^j+H = *(a'j,J(g'j+g'jM-2a'j^^) (3.33a) gfj = minmod(aV.v4,a'j+^) (3.33b) Likewise, one can define $^+^4, 3.3.6 Linearization and API decomposition In order to solve Eq.(3.32) efficiently, the left-hand side of Eq.(3.32) is first linearized by the procedure described in Section 3.1.1. The resulting scheme written in delta form is [I + Ae[(Hj,^^-Hj.^^ + (Hj^,H4-Hj^.^^) + (H^w^H-Hjw.v,)]Aq» = -M(Ej.H,k4-%Mw) + (Fj^+M4-FjW + (Gjw+H-G^^^^ . (3.34) where I is a 5 X 5 identity matrix and the H operators are given by: ^+MW ~ tAi+W'""j+K,k4]/2 Hj.k+M4 = [Bjj.+M+"j4c+v4.i]/2 (334a) in which "j+M,M = {R diag[^'-2*(a')] R-%^Ai^^

PAGE 47

... . 41 i nj^+H4 = diag[;9'-2*(a')] R-\^^a^,^ (3.34b) "jW+H = {R diag[y3'-2*(a')] R-'},+mAi+m and = *(a^+H)(g^+g^+i)/«j+M ^ ". . J^;^ ^ " )8k+v4 = *(ak+v4)(g'k+g'k+i)/«k+v4 " ' (3-34C) Here, the expression diag(z') denotes a diagonal matrix with diagonal elements z*. To calculate the n's at every time step is quite costly. For steady-state applications, the temporal accuracy is not important. One can use a spatial first, order implicit operator (that is, by setting the limiters equal zero); i.e., by redefining the n's as ! ; , «l+vum = {R diag[-*(a')] R-*}j+HAj+M ' ' ^ V nj*+M4 = {R diag[-*(a*)] R-*}k^i4Ak^j4 : v "^ ': " ^. (3.35) "iw+M = {R cliag[-*(a')] R-*},+hA,^„ : " Since the temporal accuracy is not important, the computation can be further ; reduced by simplifying the n's diagonal matrices, i.e. "j+vw = {diag[-max*(a')]}j+iiAj+^4 "j,k+vM = {diag[-max*(a')]}k^^4A,,+^4 ;^ (3.36) "ik4+M = {diag[-max*(a')]},+i4A,^^ . If one applies the approximate factorization procedure discussed in Section 3.1.2 to Eq.(3.34), the ADI form is obtained, i.e. * 1 = -M(Ej,Hw-EjW + (%*H,i-%-v44) + (Gjw+v4-Gjw.H)]" (3.37)

PAGE 48

With the use of Eq.(3.33) and Eq.(3.34a), the above equation can be expressed in deha form as follows [I + Ae5 jA+ Dj]"[I + AGS + D +xes^C+ ]"Aq" = -a(5jE + 6„F+6^G)" + De (3.38) where = -^e(nik+wa-nj,k.M4) (3.38a) = -Ae(nj^+^-nj.,kj.v4) 3.3.7 A pplication to thin-layer Navier-Stokes equations To apply the TVD scheme to the thin-layer Navier-Stokes equations, Eq.(2.21), is similar to that to the Euler equations [38-39]. The viscous term that not appears in the Euler equations is first central differenced and Unearlized as discussed in Section 3.1.1 and then added to Eq.(3.38). The resulting finitedifference equation is [I + Ae« D^]"[I + Aes^B + D + Ae(« ^ C-Re'^* ^.(J-^MJ)) + D^]"Aq" = -x(6^E+s^F+s^G Re-^5j.S)"+DE (3.39) It is noted here that A=At since A^=Aj7=Af = 1, and e is a parameter in the range of zero to one; Euler explicit or implicit scheme is obtained by setting 9=0 or 1; if 0 = 1/2, it is a second-order impHcit scheme. However, the order of accuracy is degraded due to the simpUfications that have been made in the derivation of the scheme. It is also noted that if e = l/2, Eq.(3.39) is similar to

PAGE 49

the BeamWarming scheme, Eq.(3.10), except that the nmnerical dissipation tenns are replaced by the ones derived from TVD constraints. In fact, the TVD scheme can be viewed as a central difference scheme with a "smart" numerical dissipation which is an automatic feed back mechanism used to control the amount of numerical dissipation [38]. 3.3 Flow Solver AppUcations of the original version of ARC3D and its axisymmetric version [40] for the simulation of transonic mrbulent flow past a projectile model showed that the codes implemented with the Baldwin-Lomax turbulence model can give surface pressure distribution in excellent agreement with measured data if good grid networks are provided to the Navier-Stokes codes for the simulation [25], Nmnerical experiments also showed that the Beam-Warming solution algorithm is very sensitive to grid characteristics in complex flow regions. In order to improve the robustness of the codes, a symmetric TVD scheme was first extended for modifying the axisymmetric code and investigated for the projectile aerodynamics simulation; as reported in [25,41], the modified code is indeed less sensitive to the grid characteristics in complex flow regions. Later on, the original ARC3D was modified to implement the symmetric TVD scheme. As have been mentioned before, from the viewpoint of numerical implementation, the only difference between the BeamWarming scheme and the TVD scheme are the numerical dissipations. Hence, to implement the TVD scheme into the original ARC3D, one simply replaces the dissipation operator, Eq.(3.8), and the

PAGE 50

44 " '' • dissipation terms, Eq,(3.9), with those designed for the TVD scheme, i.e., Eqs.(3.38a, 3.38b). However, the latter is more complicated than the former. New subroutines for computing the TVD dissipations had been written and the original ARC3D was then modified into a TVD thin-layer Navier-Stokes code [26]. This TVD code has been applied to the projectile flow problem, and the results showed that the code was indeed accurate and robust. Hence, it is adopted in this study for developing the multi-block computational method.

PAGE 51

CHAPTER IV THE FLOW PROBLEM AND SOLUTION STRATEGY A transonic turbulent flow of Mach 0.8 past a wing-fuselage configuration shown in Figure 4.1, is adopted as the model problem to develop the multi-block computational method. The flow problem considered is indeed a complex one, involving many viscous phenomena, such as viscous sublayer, shock boundary layer interaction, flow separation, as well as complicated boundary geometry. There are wing surface pressure measurements available for assessing the accuracy of the numerical solutions. The flow field is assumed to be symmetric with respect to the symmetry plane of the wing-fuselage configuration and hence only half of the flow domain is required for the simulation, v , 4.1 The Flow Problem and Preliminary Work > ... . • •• -j: ' ' The wing-fuselage configuration shown in Figure 4.1 is a model of a nextgeneration 150-passenger transport. The wing, optimized for Mach 0.77, has an aspect ratio of 5.17 based on the averaged chord length, wing leading edge sweep of 30", and an advanced supercritical airfoil. This wing-fuselage configuration is used as a model to develop the computational methods for propfan powered aircraft [42]. A rather coarse one block H-grid for the configuration was provided by the NASA Ames Research Center. The flow field is assumed symmetric with respect to the z = 0 plane and hence only the z > 0 half-space is 45

PAGE 52

Figure 4.1 Wing-fuselage configuration

PAGE 53

47 considered for the simulation. Figure 4.2 shows an overview of the physical domain and its computational domain. In the following, the integers j, k and 1 are used to denote the grid points in the ^, and r directions, respectively. Here, ^ is along the streamwise direction, is along the wing spanwise direction, and f is about normal to the solid surface. The dimensions of this original grid network are jmax = 101, kmax = 45, and Imax = 29. The solid smfaces are transformed into part of the 1 = 1 plane with the shape H, as shown in Figure 4.3; k=l and k=kmax are symmetric planes; while j = l, j=jmax, and l=lmax are artificial outer bovmdaries. To further illustrate the configuration quantitatively, one takes the maximum diameter of the fuselage, D, as a reference length. The fuselage is about 7D long; the wing span is about 3.3D; the chord length at the wing root is about 1.4D; the diameter of the sting is about 0.075D; the distance of the outer boimdaries far upstream and downstream are about 15D, while that of the lateral side is about lOD. These relative lengths and the j, k, 1 triplet in physical space are shown in Figure 4.4. Note also that the first spacing next to the soUd surface is in the range of 0.008D to 0.026D. The 101x45x29 H-grid one-block system was first used for a preUminary study. Because of the specific grid topology, the thin-layer Navier-Stokes code simply can not be apphed directly without modification for turbulent flow simulation. As can be seen from Figure 4.3, only part of the 1 = 1 surface is the sohd smface. Extensive modifications were made for proper Jacobians and metrics calculations, boundary conditions, and eddy viscosity computations. Also, 30 grid points that extremely close to the sohd surface were added to the original

PAGE 54

48 Figure 42 An overview of the original grid

PAGE 55

49 (a) Physical space 0 u u p NOSE FUSELAGE TOIL STING 0 R K UINQ U L 0 R UIMG U NOSE FUSELAGE TOIL STING J 0 S T Q (b) Computational space Figure 4.3 The 1=1 surface of the original grid

PAGE 56

Figure 4.4 Relative lengths in physical space

PAGE 57

, . v; , I ' / grid to resolve the thin viscous-layer; and two mirror image surfaces were also added for ease in symmetric boimdary condition treatment. The 101x47x59 Hgrid network was then provided to the modified thin-layer Navier-Stokes code for the flow simulation. As shown in Figure 4.5, the computed results are in poor agreement with experimental data. This seems to indicate that the number of grid points on the wing surface, 18 in the spanwise direction and 40 in the chordwise direction, is not fine enough to resolve the complex flow phenomna. Also, modifications on the thin-layer Navier-Stokes code may be needed before the one-block H-grid topology can be used for this complex turbulent flow simulation. A flexible and efficient multi-block technique is more desirable to numerically simulate a complex geometry turbulent flow problem, such as the wing-fuselage aerodynamics. It is also noted here that the calculations of this one-block approach were performed with the Cray X-MP/48 of the NASA Ames Research Center. 4.2 Solution Strategy For flow past a complex configuration, such as the wing-fuselage combination, to generate an acceptable single grid for the entire flow field is rather difficult. It is much easier to divide the flow field into several subregions and then generate a good grid for each subregion. Hence, it is natural to develop a solution method based on the multi-block grid network by treating each block as an independent flow problem. Each block is then solved by the same thinlayer Navier-Stokes code with the information between blocks transferred via the interface boundary conditions.

PAGE 58

Figure 4.5 Computed results of the one-block approach

PAGE 59

53 For the wing-fuselage flow problem, an effective and yet simple block-byblock computational method is proposed. In this method, the entire flow field is divided into six contiguous blocks and each block is partially bounded by a solid smface. Each block is then treated as an independent flow problem, which is solved by the same thin-layer Navier-Stokes code, with the interface boundary conditions updated at every time step. This method does not need to worry about patching solutions between different equation zone as in the zonal equation approach, and it is easier to deal with the interface bovmdary conditions than that of the Chimera overset grid approach. However, in the same manner as the zonal methods, poor interfacing can introduce errors at the interface boundaries and consequently propagated into the solution field. The error introduced can depend on the solution behavior and the grid characteristics around the interface as well as on the interpolation technique employed for updating the interface boundary conditions. For a successful multi-block computational method, it is then important and critical to minimize the associated interface error and its effect on the solution accuracy. To attack the wing-fuselage flow problem by using the proposed method, a flow simulation package has been developed which includes foiuindependent computer codes; a surface grid generation (SG3D); a volume grid generation (GG3D); a multi-block flow solver (WF6B); and a plotting program (PP3D). They will be briefly discussed in the subsequent chapter; however, a detailed description can be found in a reference manual [43]. Most of the figures presented in this dissertation were produced by the plotting program PP3D.

PAGE 60

54 4.3 Six-Block Grid Topolog y For a multi-block computational method, the definition of a good composite grid topology for a complex configuration turbulent flow simulation is more involved and difficult than that for an inviscid flow or a laminar viscous flow simulation. For instance, the solid surface has to be transformed onto an entire plane in computational space so that the thin-layer approximation can be applied easily; the Baldwin-Lomax turbulence model requires the normal distance of each grid point to the soUd surface; and the grids should be able to cluster points close to the solid surface to resolve the thin viscous-layer. Hence, for ease in eddy viscosity calculation and in applying the thin-layer approximation, a special 6-block grid topology was defined and investigated for the wing-fuselage configuration. The main idea of dividing the flow field is to isolate the wing by a cone-like boundary surface that sits between the wing and the fuselage, as shown in Figure 4.6. The isolated wing was divided into 2 blocks and the rest of the flow field was divided into 4 other blocks. By doing so, the flow field is divided into six contiguous blocks such that each block is partially bounded by a solid surface, and the solid surface of each block is transformed onto f = 1 plane in the computational space. An overview of the topology is shown in Figure 4.7 while Figure 4.8 gives side and cross-sectional views. The physical and computational space of each block is shown in Figure 4.9. The selected grid topology has the clear advantage that each block can be solved by the same thin-layer Navier-Stokes code; however, the choice has its deficiences in the block grid generations. For instance, for the wing blocks. Block

PAGE 61

55 Figure 4.6 The dividing surface isolating the wing

PAGE 62

Figure 4.7 Six-block grid topology (overview)

PAGE 63

57 Figure 4.8 Side and aoss-sectional views

PAGE 64

58 z (a) Block 1 (b) Block 2 Figure 4.9 The transformation of the six block grids

PAGE 65

(d) Block 4 Figure 4.9 ( continued )

PAGE 66

(f) Block 6 Figure 4.9 ( continued )

PAGE 67

61 3 and Block 4, shown in Figure 4.9(c,d), four of the six bounding surfaces in the computational domain, CLRQ, LNTR, RTSQ, and QSHC form a nearly flat surface in the physical space. Therefore, special measure and care in grid generation are required to overcome the difficulty resulting from the drastic change in the domain transformation. Also, for the nose block. Block 1, shown in Figure 4.9(a), the boundary plane AEJF in the computational domain degenerated into a single line in physical space and care must be taken when dealing with the boundary conditions.

PAGE 68

CHAPTER V GRID GENERATION For the wing-fuselage flow problem, the flow domain has been divided into six blocks, however, it is still a great challenge to generate acceptable grids for each block. To generate grid networks for the six-block wing-fuselage flow problem, surface grids must be obtained first. A special surface grid generation code (SG3D) was developed to utilize the original H-grid network for constructing the required surface grids. The resulting surface grids will be presented in section 5.1. Generating volume grids for the six-block flow domain, it is required to develop an adequate grid generation code. Various grid generation methods [44-51] have been studied in that regards. The main concern is that the code can produce grid networks that satisfy the required grid characteristics, such as: (a) grid lines in different directions must be nearly orthogonal to each other; (b) grid lines in {--direction must be about normal to the solid surface; (c) r = constant grid surfaces must be clustered near the solid surface; (d) grid lines must be continuous across the interfaces; and (e) the blocks must be contiguous (no holes or overlapping). A combination of two-boundary algebraic and elliptic grid generation methods will satisfy these requirements best. Accordingly, a general purpose 3-D Grid Generation code (GG3D) that based on these two methods was developed and applied to generate the 3-D grid networks. 62

PAGE 69

A detailed description of GG3D can be found in [43]. The techniques employed to generate the six block grids are described in Section 5.2 and the resulting grid networks are presented in Section 5.3. 5.1 Surface Grid Generation To generate grid networks for the six block grids, surface grids must be obtained first. A two-boundary algebraic grid generation method is used to generate volume grids; consequently, only the solid surface grid (1=1) and the cooresponding outer surface grid (l=lmax) are required. The other boundary surfaces are then defined by the straight Unes that connect the boundary curves of the solid and outer surface grids. A special surface grid generation code (SG3D) has been developed that uses the original H-grid to construct each surface grids by cubic splining, transfinite interpolation, strentching, shifting and dividing. The resulting surface grids obtained are shown in Figure 5.1; the dimension of these surface grids are 28x23, 40x12, 40x35, 40x35, 40x12 and 26x23 for Block 1-6, respectively. It is noted here that the outer surfaces are defined according to the grid topology and will have strong effect on the grid characteristics of the block grid to be generated. Hence, the boundary curves of each outer surface should be defined carefully. For block 1, a cone-like shape that is resemble to the solid surface is desired for the outer surface. As shown in the grid topology. Figure 4.7, the boundary curve GHI is defined such that the bounding surface lies in the middle of the wing leading edge and the fore-body of the fuselage. Similarly, the outer surface of block 6 is defined such that bounding surface of the boundary

PAGE 70

(a) Solid surfaces Figure 5.1 The 6-bIock surface grids

PAGE 72

66 curve MNO lies in the middle of the wing trailing edge and the aft-body. The boundary curves Hb'N of block 2 and Ha'N of block 5 are also defined such that the boundary surfaces Hb'NLC and Ha'NLC lie in the middle of the wing and the fuselage. For the outer surfaces of block 3 and block 4, the boundary curve HN, that is defined as a straight line, is located such that the interface surface of block 3 and block 4, i.e., surface HRCQLN, is not incline to either upper side or lower side of the wing. In other words, the plane spanned by HCLN should pass throught the wing tip, and hence divides the wing into upper wing and lower ; wing. If this were not the case, one would have difficulty to generate an usable volume grid for either block 3 or block 4. » ' ' . 5.2 Volume Grid Generation 'C * '. -. -" ^ 5.2.1 Two-Boundary Grid Generation The two-boundary grid generation technique used is based on Hermite transfinite interpolation [46-48]. Two surface grids are specified and the interior grid lines are connected by curves such that boundary orthogonality is enforced. The other four boundary surfaces are then defined using straight lines. The distribution of interior stufaces is defined by a suitable clustering function such that more grid points are clustered near the solid surface as it is usually required. A description of clustering functions is given in Appendix C. Denoting the position vectors of the two specified boundary siufaces as r(?,r;,0) = r,(^,,,) and r(e,»?,l) = r2(^,,,) (5.1) and the derivatives with respective to f as

PAGE 73

67 ar the following functional form, based on Hermite transfinite interpolation, may be written for the connecting curves: ra.'/.r) = ri(|,„)hi(f) + r2(^,„)h2(f) + DRih3(r) +DR2h4(f) (5.3) To satisfy Eq.(5.1) and Eq.(5.2), the following must be true: h.(0) = 1, h.(i) = 0, iMOL = 0, iMlI = 0 h,(0) = 0. h,(l) .1.^=0, ^ = 0 h,(o) = 0, h3(i) = 0, ^M2L = 1, iMlL = 0 ^ ar ar h,(0) = 0, h,(i) ,0,^=0, iMl). = 1 af ar Hence, a cubic polynomial in {• for hj(f ) satisfying conditions given in Eq,(5.4) can be determined. They are hi = 2r^-3r^+l h2 = -2r'+3f2 ha = r'-2r'+r (5.5) Here, f is a normalized point distribution, and is usually defined by a clustering function as shown in Appendix C. The expressions for the derivatives, DR, and DRj, appearing in Eq.(5.3) are the key elements used to enforce the boundary orthogonahty and are derived as follows.

PAGE 74

68 A vector normal to the surface r = 0 that points into the spatial domain is given by where x denotes the vector cross product. A vector that tangents to the connecting curve at r=0 can be defined as ar(e,>?.0) e = ac Hence, in order for the connecting curves to be normal to the r = 0 surface, the cross product of the vector tangent to the connecting curves e and the surface normal n must be zero, i.e., exn = Oatf=0 (5.5) or written in component form r .^ci^_-^ \y ar a»7 ai ac ari ar at; a^ "aT ^ af ^ ar; ac ' af at; ^" ar a»; a^ "ai~ "a^^ + f iL. ( IL. .££_ ^ :?y_ / iy_ ar ar? a^ ac ar? ar ar? ae . = 0 (5.5a) From the above equation, it can be seen that the connecting curves will be perpendicular to the surface r = 0 if the following conditions are satisfied at f = 0: DX, = ^ = -CK,(e,.)( JLJL.^JL^ ar ^^^"^^ dr, ar, '

PAGE 75

69 DYi — CKi(e,»7)(-— — --— — ) ' (5.7) ari di arf rT and C and ri are normalized as follows:

PAGE 76

70 CG) = (j-l)/(jmax-l), j = l,2,...,jmax r,(k) = (k-l)/(kmax-l), k = l,2,...,kmax (5.10) The CK's appearing in Eqs.(5.7-5.8) have a strong effect on the shape of the connecting curves. In fact, they control the depth of the perpendicular part of the curves. Values for the CK's are usually specified as constants by trial and error. A shape function was devised to modify the CK's so that different values can be specified at different j, k locations. The shape function is of the form CK = CKI*cosx{CPJ[0-l)/Omax-l)] + CSJ} > «cos^{CPK[(k-l)/(kmax-l)] + CSK} (5.11) where CKI is a constant to be specified, CPJ and CPK are used to control the frequency while CSJ and CSK are used to control the phase of the relative cosine function. In doing so, values for the CK's become functions of j and k, and can be smoothly distributed on the surface. For example, if one specifies CPJ = 0.5, CSJ = 0, then CK=CKI when j = l and CK decreases as j increases; or if one specifies CPJ = 1, CSJ =-0.5 then CK increases from zero to CKI and then decreases to zero as j varies fi-om 1 to jmax. Figure 5.2 shows some examples to illustrate this shape function. 5.2.2 Elliptic Grid Generation The technique and subroutines employed are based on the elliptic grid generation system of Thompson [49-51]. The volume grid is constructed by forcing the Cartesian coordinates of interior points to satisfy specially devised Poisson equations. In the Poisson equations, three control functions for each curvilinear direction are set up to control the grid point distributions. The system

PAGE 77

71 CP j-e. 5 , cs j-e , cPK-e , csK-e cp >e , cs j-0 . cpk. s , csk— . 4 CPU-. 8 , cs J— a . 4, CPK. 5, csK— . 5 CP j-e, cs j-e, cPK-a, csK-9 Figure 5.2 Shape function for CK's

PAGE 78

72 of Poisson equations are solved by point SOR (success overrelaxation) iterations. The user-provided initial volume grid is used not only to start the SOR iterations but also to set up the control functions. The elliptic grid generation system is based on the Poisson equations of the form = (v^.v^Pi v^r, = (v„.v„)P2 (5.13) A = (vr-vr)P3 where the P's are so-called control functions which can be specified to control the grid smoothness and/or orthogonahty. Since the actual computations are performed in the computational domain where ^, rt and $• are the independent variables, with r = (x,y,z)^ as the dependent variables. The poisson equations, Eq.(5.13), are then transformed to computational domain. The transformed equations £u-e a^r a^r d^v (V^.V^) __ + (V„.V„) — + (Vr.7„) -— -H a^OT) dtjdr} d^dr/ a^r a^r a^r acsr a»7a$' a^ar (7e.vOPi^+(Vr,.V„)P2^+(Vf.Vf)P3^ = 0 (5.14) OS or/ af The finite-difference equafions of Eq.(5.14) are obtained by replacing the derivatives with central differences and these equations are then solved by point SOR iteration procedure.

PAGE 79

73 If the control function P's are evaluated from the initial grid, the three components of the elliptic grid generation system, Eq.(5.14), provide a set of three equations in which the control functions P's are treated as unknowns and the values for the position vectors are obtained from the initial grid. The transformed poisson equations, Eq.(5.14), can be rewritten as (vevOPi + (v»7 • v,,)P2 ^ + (vr . vr)P3 d^d^ dfidi arac / X / X , , aV (VC-Vr,)— +(V„.V„)— — +(vr.V„)— + a^drj drjdri d^drj which can be solved simultaneously at each point for the three control functions, Pj, i=l,3. If the control functions are computed solely based on the above equations, and the elliptic system, Eq.(5.14), will reproduce the same grid as the initial grid which is of trivial interest. However, if the computed control functions are smoothed before they are applied into the elliptic system, a smoother grid network can then be produced. A variation to the evaluation of the control functions is obtained by assuming that the grid lines are orthogonal to each other, and Eq.(5.15) is reduced to (ve.vOPi-^+(V,.V,)P2i^^+(Vr.Vf)P3^ = oT) a$" -[("^^'^O 4hc + (^''•^'') -fr + (^^-"^f ) -TT] (5.16) aca? dr}dri acaj^

PAGE 80

The control functions evaluated from these equations will force the elliptic system to produce a more orthogonal grid. 5.3 The Six-Block Grid A set of volume grids with 60 grid points in the I-direction was first generated by two-boundary grid generation with grid hues set normal to the solid surface whenever possible. The first spacing off of the solid surface was set to about 2xlO'^D, which is considered fine enough for turbulence simulations. Figure 5.3 shows the 3-D close-up plots of the six block grids with cut-offs to show the details of the interior surfaces. These grids were generated from the soUd and outer surfaces described in Section 5.1. The other four boundary surfaces were then defined by straight lines that connect the boundary curves of the solid and outer surfaces. Sixty grid points were distributed on each of the connecting straight lines by using the hyperboHc tangent clustering function with the first spacing of 2xlO'*D and the last spacing of six times of the average spacing, and hence, more grid points were clustered near the solid surfaces. To generate grid network using two-boundary grid generation, one has to face a difficulty in specifying a proper value for the parameter CK, Eqs.(5.7-5.8). In general, the value for CK is specified by trial and error, however, its variation on a surface can be justified by the use of the shape function, Eq.(5.11), according to the geometry of the grid to be generated. For the six block grids, the grid orthogonahty on the outer surface was not enforced, that is CK2 = 0, to avoid the possibiUty of grid overiapping. On the other hand, the grid orthogonality on the

PAGE 81

(b) Block 2 Figure shows the surfaces 1=1, 1=50, k=3, k=13, j = 17 Figure 53 Six block grids 3D close-up with cut-off

PAGE 82

76 (c) Block 3 ~ Figure shows the surfaces 1=1, 1=50, k=l, k=15, k=30, j = 10, j = 18, j=32 Figure 5.3 (continued)

PAGE 83

77 (e) Block 5 -Figure shows the surfaces 1=1, 1=50, k=3, k=13, j = 17 (f) Block 6 ~ Figure shows the surfaces 1=1, 1=50, k=5,k=15,j=l,j=10 Figure 53 (continued)

PAGE 84

78 solid surface is desirable. From numerical experiments conducted the following was one of the best to set up the CKj distribution by the shape function for each block. Block 1 CPJ =0.5 CSJ= 0 CPK= 0 CSK= 0 Block 2 CPJ = 1 CSJ = -0.5 CPK= 0.5 CSK= 0 Block 3 CPJ = 0.8 CSJ = -0.4 CPK= 1 CSK= -0.5 Block 4 CPJ =0.8 CSJ = -0.4 CPK= 1 CSK= -0.5 Block 5 CPJ = 1 CSJ = -0.5 CPK= 0.5 CSK= 0 Block 6 CPJ = 0.5 CSJ = -0.5 CPK= 0 CSK= 0 These distributions are shown in Figure 5.4. For block 1, CKj decreases as j increases and becomes zero when j = jmax. This is reasonable since the j = jmax surface incline about 45 degrees to the solid surface. Similar argument can be addressed for the other blocks. To further improve the smoothness of the grid network, these algebraic grids were then used as initial grids for elliptic grid generation. Smoother grids were obtained for block 1, 2, 5, and 6. However, it failed to generate elliptic grid for either block 3 or 4. This may due to their special geometries such as three boimdary surfaces form nearly a flat plane. Figure 5.5 illustrates the difference between algebraic and elliptic grids. Notice that the elliptic grids are smoother and more orthogonal than the algebraic grids. These six block grids were further modified to include mirror image planes for imposing symmetric boundary conditions in the Navier-Stokes simulations. The resulting grid network was then used as a basic grid in this stody and is referred to as the "base grid" hereafter. The dimensions of this base grid are 28x25x60, 40x13x60, 40x35x60, 40x35x60, 40x13x60, and 26x25x60 for blocks 1-6,

PAGE 85

Figure 5,4 Shape function for the six-block grid

PAGE 86

1 1 I 80 2 Block 1, k-13, fllgsbralc Grid Z Block 1, k-13. Elliptic Grid (a) Block 1 k=13 Figure 5.5 Elliptic and algebraic grids

PAGE 87

81 (b)Block6-k=13 Figure 5.5 (continued)

PAGE 88

^ 1 82 respectively. Figure 5.6 shows some detailed plots to illustrate the grid characteristics. Figures 5.6a-d show that the grid lines are about normal to solid surface whenever possible. Figures 5.6e-f show that some k= constant surfaces are quite skew for block 2, 3, 4, and 5. Figure 5.6g shows the grid distribution at k = l, 12, and 18 wing cross-sections. In order to investigate the effect of grid resolution on the flow solution, a series of grid refinement studies have been performed. Several locally refined grids (near the wing tip or near the soUd surface, etc.) were generated and used to smdy the possible effect due to the grid characteristics. A set of "refined grids" was also generated by doubling the grid points in the wing core-wise direction. The dimensions of these refined grids are 28x25x60, 79x13x60, 79x35x60, 79x35x60, 79x13x60 and 26x25x60. Figure 5.7 shows the solid surfaces of block 2 and 3 of this refined grid. Note that the grid lines were generated by using a cubic sphne so that the soUd surfaces would not be distorted.

PAGE 89

83 L j-1 or 48 Block 3 or 4 k-35 J-40 or 1 Slock 1 k-13 Block 6 k-13 (a) Block 1, 3, 6 (b) Qose-up to show the boundary orthogonality Figure 5.6 Some detail plots for six-block grids

PAGE 90

84 Block 2 Block 5 J-20 j-2e (c) Block 2, 3, 4, 5 (d) Close-up to show the boundary orthogonality Figure 5.6 (continued)

PAGE 91

85 (e) Block 3, 4 k surfaces are quite skew (f) Block 2, 5 k surfaces are quite skew Figure 5.6 (continued)

PAGE 92

86 Slock 4 k-18 Block 3 k-18 Block 4 k-12 Block 3 k-1 1-1-35 (g) Block 3, 4 k=l, 12, 18 to show the grid distribution on the leading and trailing edge Figure 5.6 (continued)

PAGE 93

87 Base Grid 48x35 on uing surfaca (a) Base grid Block 3, 4 40x35 Block 2, 5 -40x13 Refined Grid 79x35 on wing surface (b) Refined grid Block 3, 4 79x13 Block 2, 5 79x13 Figure 5.7 Base grid and refined grid

PAGE 94

CHAPTER VI SOLUTION METHOD The multi-block computational method considered is rather simple and straightforward in application. For the wing-fuselage flow problem, the flow field has been properly divided into six contiguous blocks, and the grid network has also been generated. In the solution process, each block is then treated as an independent flow problem that is solved by the same thin-layer Navier-Stokes code. The interface boimdary conditions between blocks are updated with a distance weighted interpolation before the computations march for another time step. A computer code has been developed and investigated for the six-block wing-fuselage aerodynamics in which the unsteady 3-D TVD thin-layer NavierStokes code is made to a subroutine subprogram. Since the multi-block computer code has to deal with several blocks of different configurations and their interfaces, the code has to be modular and flexible. Also, the code has to adapt to the architecture of the host computer system. A six-block wing-fuselage Navier-Stokes code (WF6B) was then written for the available computer system, Cray-2 of NAS system at NASA Ames Research Center. The special feature of the supercomputer, Cray-2, is its huge main memory, however, the computing speed is a Uttle slower than Cray Y-MP. In order to take advantage of the large main memory while not sacrijSce the computing speed, WF6B uses part of the 88

PAGE 95

89 main memory for storing the intermediate data. The other main memory is then used as active area for executing computations. A set of common arrays is then allocated at the active area and used for the computations of each blocks. A block diagram to show this data arrangement is given in Figure 6.1. A detailed description of WF6B can be found in [43]. In order to illustrate the multi-block solution method a flowchart of WF6B is given in Figure 6,2. As shown in the figure, two do-loops run over each blocks and the nmnber of time steps requested. In the inner loop, for each block, BC is used to set up the boundary conditions and then TLNS is to formulate the governing equations and perform the three inversions. A value of L2 residual for each block is computed by RESID every 10 time steps. After the computations run over each block, INTERBC updates the interface boundary conditions, and then the computation continues to the next time step. Finally, the pressure and the shear stress coefficients are computed by CPP and CSP, respectively. In what follows, the techniques used to calculate the metric coefficients and Jacobians of transformation are described first. Then the boundary conditions imposed on the boundary surfaces of the six-block flow problem are discussed while the interface boundary conditions are given in a separate section. Fmally, the equations used to compute the 12 residual and pressure and shear stress coefficients are given. 6.1 Jacobian and Metric Coefficient The metric coefficients and Jacobians of transformation can be calculated from Eqs.(2.16-2.17), if the covariant metrics, x^, y^, z^, etc, are known. For a second-order approximation, the covariant metrics are calculated at interior grid

PAGE 96

90 Block 1 28x25x60 rl(3, j,k,l) ql(j,k,l,6) Block 2 40x13x60 r2(3, j,k,l) q2(j ,k,l,6) Block 3 40x35x60 r3(3,j,k,l) q3(j,k,l,6) Block 4 40x35x60 r4(3, j,k,l) q4(j,k,l,6) Block 5 40x13x60 r5(3, j,k,l) q5(j ,k,l,6) Block 6 26x25x60 r6(3,j,k,l) q6(j ,k,l,6) Note: r's are the Cartesian coordinates q's are the unknown vectors Block data sweep in Active Area = 40x35x60 Common Blocks /van/ r(3,j,k,l) /vars/ s(j,k,l,5) q(j,k,l,6) /btrid/ a (3, j, 1,5, 5) /vis/ tuirmu( j ,k, 1) /var3/ XX (9,1) etc. r-Cartesian coordinates q-unknown vector s-RHS vector a-5x5 matrices turmu-eddy viscosity xx-metric coefficients results sweep out read in grid data and/or restart file write out results Permanent Disk stores r & q for each block Figure 6.1 Data arrangement in WF6B

PAGE 97

91 read in grid data compute Jacob ians & generate initial condition or read in restart data do N = 1,NMAX do BLOCK =1,6 BC TLNS RESID INTERBC write out if CPP = yes CPP if CSP = yes CSP grid data restart data residual restart data pressure coefficients shear stress coefficients Figure 6.2 The flowchart of WF6B

PAGE 98

92 points by a central-difference formula and second-order one-sided differencing at the boundaries, e.g., = (Xj+ijc4-XiW/2 for 1 < j < jmax Xf = (-3Xj,m+4Xj+i^-Xj+2j^/2 for j = 1 (6.1) = ( 3Xj^-4Xj.ij^+Xj.2j^/2 for j = jmax The Jacobian, J, of each grid cell is then computed from Eq.(2.17). The contravariant metric coefficients are defined in Eq.(2.16), and can be calculated accordingly. The Jacobians and the metric coefficients evaluated by the above finitedifference approach are rather sensitive to grid characteristics, and can produce a large truncation error if the grid network does not have good properties of smoothness and orthogonality. It is known that they can also be approximated from geometric relations such as cell volume and cell-face area. In the , ^ discretized sense, the Jacobian defined in Eq.(2.17) is equivalent to one over the cell volume, 1/Vol, and the contravariant metric coefficients defined in Eq.(2.16) are equivalent to one over the cell-face normal. If one defines a computational cell that is centered at the corresponding grid point, with its cell-faces located half-way to the neighboring points, it is then appropriate to replace the Jacobian, J, by the corresponding cell density, 1/Vol. Furthermore, since C^/J, ^y/J, and iji defined the x, y, z components of the normals (not unit normals) to the local tangent plane of the constant | surface, one can replace by n^/Wo\, Uj-y/Vol, and Uj^/Vol, where n^^ Ujy, and are the components of the normals to the local constant j planes. Likewise, one can replace rj^ t)^ and c,, fy, by

PAGE 99

93 n^/Vol, n,jy/Vol, n,„/Vol and n^j/Vol, n,y/Vol, n,^/Vol, respectively. A general hexahedron, defined by eight arbitrary comer points, can be partitioned into six tetrahedrons, Figure 6.3a. The volume of a tetrahedron denoted by its vertices (a,b,c,d) equals one-sixth of the triple product of the three vectors emanating from one of the vertices and ordered according to the righthand directions. Figure 6.3b. For example, V'*(a,b,c,d) = [r.,.(r^ x r,,)]/6 (6.2) = [(Xc-Xd)(yb-yd)(z,-Zd) + (Xb-Xd)(y.-yd)(z,-Zd) + (Xa-Xd)(yc-yd)(Zb-Zd)-(vxd)(yb-yd)(Zc-Zd) (vxd)(ya-yd)(Zb-Zd)-(vxd)(yc-yd)(ZaZd)]/6 where = Tj Fj, and Tj and Tj are the position vectors. If the eight vertices of the hexahedron are denoted by numerals 1-8 and are associated with the j,k,l triplets as follows: l-j,k,l 2-j + l,k,l 3-j+l,k+l,l 4-j,k+l,l 5-j,k,l+l 6-j + l,k,l+l 7 j+l,k+l,l+l 8-j,k+l,l+l the volume of the hexahedron is Voiak,l) = V*'*(l,2,3,7) + V***(l,2,7,6) + V***(l,5,6,7) + V*'*(l,3,4,7) + V*(l,4,8,7) + V*''(l,5,7,8) (6.3) The Jacobian of a grid-centered computational cell, which is bounded by the surfaces at the middle of two adjacent grid points, is computed as one over the average volmne of all the associated hexahedrons.

PAGE 100

94 1 . -» . s • 7 / / / / 1 / 4 s s a VolOMc,!) = V«(U3,7) + V-(1A7,6)+V^(1,5,6,7) + \^(1,3,4,7) + V*«*(1,4,8,7)+V*-(1,5,7,8) (a) Partition of a hexahedron d a ^ c V-(a,b,cd) = [r,,.(rMxrJ]/6 (b) Volume computation of a tetrahedrons Figure 6.3 Volume computation of a hexahedron

PAGE 101

95 The metric coefficients, based on the cell-face area, are calculated in a similar maimer. The area vector of a triangle denoted by its vertices (a,b,c) is A^(a,b,c) = (r^b X r^)/2 = {[(yb-yJ(vzb)-(yb-ya)(Zc-Zb)]» + [(Zb-Za)(VXb)-(Zb-Za)(VXb)]j + [(vxa)(yc-yb)-(vxa)(yc-yb)]k}/2 (6.4) The area vector of a general surface, defined by four comer points (a,b,c,d), may be defined as AREA(a,b,c,d) = A^(a,b,c) + A''(c,d,a) (6.5) Hence, nj^ are the x,y,z components of AREA(a,b,c,d). For a constant j plane, one has %x = {[(yb-ya)(Zc-Zb)-(yb-y.)(Zc-Zb)] + [(yd-yc)(vzd)-(yd-ye)(vZd)]}/2 njy = {[(Zb-Za)(Xc-Xb)-(Zb-Za)(x,-Xb)] + [(Zd-Zc)(vXd)-(Zd-Zc)(vXd)]}/2 /• (6.6) Djz = {[(vxa)(yc-ybHvxJ(yc-yb)] + [(Xd-Xc)(ya-yd)-(Xd-Xc)(ya-yd)]}/2 and similarly, one can compute n^^ n^y, n^, and n^^ n,y, n,^. 6.2 Boundary Condition There are different types of boundaries encountered in nmnerical simulations, such as solid surface, far-field, inflow/outflow, and symmetric boundaries. The type of boundary condition specified at a boundary surface

PAGE 102

96 depends on the physical behavior and numerical requirement at the surface. According to the grid topology, the proper boundary condition has been set up for each block in WF6B. Figure 6.4 shows the boundary conditions applied at the boundary surfaces of each block and the relative interfaces between them. In the figure, SOLIDBC, FSBC, SYMBC, SGLRBC, OUTBC and INTERBC stand for soUd, freestream, symmetric, singular, outflow and interface boundary condition, respectively. The j = 1 surface of block 1 is a special surface that degenerates into a single line in physical space, special treatment is required. Since on the surface, each ri line maps to a point in physical space, the same boundary condition is set at each point on this line, and the value is obtained by distance weighted extrapolation from the interior points. The other boundary conditions are discussed in the following and the interface boundary condition is given in the next section. 6.2.1 Solid Boundary Condition In the computer code, the r = 0 plane is assumed the solid surface. The techniques employed to obtain the velocities, pressure, density, and total energy will be described for both inviscid and viscous flows, and then a slow start mechanism for viscous flows is given. Inviscid flow. At a soUd surface, the flow tangency condition is satisfied. So the contravariant velocity component in f -direction W is set to be zero for non-injecting wall and the other contravariant velocities U and V are foimd by extrapolation, e.g..

PAGE 103

97 -1=1 iS. — X 1 =1 X X X — X ilia A. X
PAGE 104

98 Vjw = 2V^^ V.,^ (6.7) The velocity components u, v and w are then obtained from Eq.(2.20c). Note that ^, , r;, and are zero for a fixed grid system employed here. The total energy E at the solid surface is obtained from Eq.(2.7) in which the pressure is computed from the normal momentum equation. The normal momentum equation is obtained by combining the three transformed inviscid momentimi equations [10], i.e.. a(pu/J) 1 a[(puU+^j))/J] 1 a[(puV+r,J))/J] _ a[(puW+rj5)/J] = 0 at a(pv/J) ^ 3[(pvU + e,p)/J] ^ a[(/.vV+„,p)/J] ^ a[(pvW+r,p)/J] = 0 at ae ar; ar a(pw/J) 1 a[(pwU+ej))/J] 1 a[(pwV+r,j))/J] 1 a[(pwW+rj))/J] = 0 at a€ a»? ar If combining r^'C-momentum + ry •'/-momentum + r^'r-momentum with the use of the continuity equation and metric identities, Eq.(2.19), one has the normal momentum equation Uxrx+ Cyry + i.r,)Pf + (»?,rx+ Vy + "zf z)p.+ (rxrx+ ryr, + rzrz)Pf = -pUCr^u^+ryVf+r.Wf^pVCr^u.+ryV^+r.w,) (6.8) For simplicity, the equation is written as aJ>j + a2P^ + a3P^ = R ^ /^^^ " ^ ' (6.9) with "1 = (^xrx+Cyry+^zrJ " \ . "2 = fex+'J/y+'JirJ •

PAGE 105

99 °3 = (rxfx+ryr,+fzrz) R = -pU(5:^u5+ryVj+rzWe)-pV(f,u^+ryV^+r,w^) (6.10) Note that the left-hand side of Eq.(6.8) equals vp«Vf which is the normal pressure gradient without normalization. The normal momentum equation, Eq.(6.9), is solved at the solid surface by using ADI decomposition in the C and ij directions and second-order forward differences in the f direction. A second-order forward difference is first applied to Pj., i.e., Pf = (-3Pj,M + '*Pjw-Pj,M)/2 + O(Ar') (6.11) The normal momentum equation, Eq.(6.9), can then be written in semi-discritized form as 2ti 2c( 2R Pik,r 3;^ Pr 17P. = ('^Pjw-Pjw— )/3 (6.12) If the ^and r/derivatives are central differenced, then the finite-difference equation of Eq.(6.12) expressed in operator form is 2a, 2a, 2R . Pj.k.r 3—«fPj,k,ri— ^.Pjjci = (4Pj,M-Pjjc,3 )/3 or ^^'^'^'^ '"^Pj'^i = (4p^,K,2-Pj,M^ )/3 (6.13) By applying ADI-de composition to the left-hand side of Eq.(6.13), one has ^ 1?'"^PJ'M = (4Pjw-Pjw^ )/3 (6.14) which is second-order accurate in each spatial coordinate. The pressure at the solid surface can then be solved by the following two one-dimension inversions

PAGE 106

100 (|^Vl)P* = (4pj^-Pi«-|^)/3 (^Vl)Pj.M = P* (6.15) The left-hand side of Eq.(6.15) forms systems of tridiagonal matrices which are solved by the Thomas algorithm. The density at the solid surface is obtained via the ideal gas equation of state when the pressure and temperature are known. For an isothermal wall at constant temperature T^ the density is computed from the dimensionless equation of state P = ypTJT^ (6.16) For an adiabatic wall, the stagnation enthalpy is held constant. By using the equation for enthalpy H„ = E+p/(p,p) = E„+pJp„ (6.17) and the computed velocities and pressure, a value of density is obtained at the solid surface, i.e., _ yppa ' ~ (7-l)[(E.+p,)-,.(u^+v^+w^)/2] (6-18) Viscous flow. For viscous flows, the no-slip condition is enforced by setting all the velocity components equal to zero on the impermeable wall, i.e., UjW = Yiw = = 0 (6.19) The total energy is again computed from Eq.(2.7) in which the pressure and density are obtained as follows. Since the viscous terms have a negligible effect on the surface pressure, the pressure distribution for viscous flows is obtained by the same procedure as

PAGE 107

101 that for inviscid flow. It is worth noting that for viscous flows, U and V are both zero on the solid wall; and hence, the normal momentum equation being solved becomes Uxfx+^yfy+e.rJPe+('7xrx+'7yry+'?,fJp,+(r,f«+ryry+fxrz)Pf = o (6.20) In other words, based on the inviscid normal momentum equation, Eq.(6.9), the solid boundary condition for the viscous flow is an |vr| The density for an isothermal wall at constant temperature is computed from the same dimensionless equation of state, Eq.(6.17). However, for an adiabatic wall, extrapolation is used for viscous flow problems, e.g., PiM = f>w (6.21) since both pressure and temperature gradients are zero at the wall, and dp =0 at a wall an Slow start. In viscous flow problems, the flow field is usually started from an initial condition of uniform free-stream flow. To avoid the sudden jump on the solid surface, the no-slip condition is turned on slowly over the first several time steps. In the code, the boundary conditions are smoothed over the first 30 time steps by setting the dependent variables q at the sohd wall as q = <^ + (l-«)q. (6.22) with u> = (10-155 + 6^')^^ 9 = (NC-l)/30 (6.22a) where NC is the number of time steps and q^^ is the correct boundary conditions

PAGE 108

102 discuss above [28]. Equation (6.22a) is a special function designed for w such that w increases very slowly over the first few time steps. 6.2.2 Far-Field Condition , In order to numerically solve a problem that is posed in an exterior region, the infinite domain is commonly replaced by a finite domain with artificial bound£uies. Stretched grids are usually used to place those boundaries far away from the object such that no spurious reflected waves can be generated. Hence, the free-stream conditions can be specified as the far-field boundary conditions [30]. • 3 . ' . The subroutine FSBC is used to impose the far-field boundary conditions at 1= Imax by setting all the dependent variables q at those boundaries equal to free-stream condition, i.e., . .. ' '-'l i' q = q« ' ' ' " : ' (6.23) 6.2.3 Inflow/Outflow Condition ; \ The inflow and outflow boundaries may be viewed as special cases of farfield boundaries, and hence, the non-reflecting property is also required. Since the incoming and outgoing directions are known at those boundaries, it is easier to specify the boundary conditions based on the wave theory [30]. When the flows at the inflow boundaries are supersonic, the flow conditions are not affected by the downstream flow characteristics and hence can be set equal to the freestream conditions, i.e., q = q.

PAGE 109

103 For supersonic flow at outflow boundaries, extrapolation can be employed since no influence should be felt upstream in the region of interest, e.g., at j = jmax 4jniax,k44 ~ ^jinax-l,k44 (6.24) For subsonic flows the situation is more complicated. More sophicicated boundary conditions may be needed if strong reflection waves are expected at the artificial boundaries. Here, only the techniques used in the original ARC3D are described. In ARC3D, the free-stream condition, Eq.(6.23), is employed at inflow boundaries while at the outflow boundaries, the velocity components u, v, w and the density p are determined by extrapolation and the total energy E is obtained from Eq.(2.7) by assuming the pressure is equal to the free-stream pressure, e.g., at j = jmax ^jinax,k4 ~ Ujmax-lJtJ X)max,k,l ~ yjmax-l,k4 ^jmax,M ~ ^jmax-l,k4 (6.25) E = pV(7-1) + ^p(u2+v2+w2) 6.2.4 Symmetric Boundary Condition The symmetry boundary is commonly adopted to simplify flow problems if symmetric conditions can be assumed. In the computer code, the »? = 0 plane is assumed to be a symmetry plane which is also a constant z plane in physical space as shown in Figure 6.5a. On this plane, the V-velocity and the first

PAGE 110

104 (a) Symmetric boundary Synmetry Plane Y z-0 Mirror image plane k-32 2-0 A. 1 Mirror imoge plane (b) Symmetric boundary with mirror image plane Figure 6.5 Symmetric boundary condition

PAGE 111

105 derivatives of the other flow properties with respect to n are set to zero to enforce symmetry, i.e., qj = 0 %W = (Hmw qj,k+ui)/3 for i ^ 4 (6.26) Recall that qj^y^ is pw at grid point and V=w on this symmetry plane. Another way of dealing with symmetry boundaries is to add one more plane that is the mirror image of the plane next to the symmetry plane. For instance, if k = 2 is the symmetry plane, then k = 1 is the added boundary plane which is the mirror image of k = 3 plane. As shown in Figure 6.5b, if the k =2 plane is also a constant z plane then the flow conditions at k = 1 plane can then be specified as %W = q^Ai for i ^ 4 (6.27) to enforce symmetry. Note that q^y^ is (pw)j^y and equals -(pw)^^ since k = 1 plane is a constant z plane. 6.3 Interface Boundary Condition The success of the multi-block grid method relies on correct communication among blocks. Inappropriate implementation of interface boundary conditions may result in numerical instability or reduced accuracy. There have been various techniques developed for updating the interface boundary conditions with different degree of sophicication [52-53]. A simple and efficient interpolation method has been successfully tested on a projectile flow problem. In this approach, the information at the interface boundaries is obtained by distance weighted interpolation from the neighboring blocks. For

PAGE 112

106 instance, the unknowns q at an interface boundary point can be estimated from the known quantities at an interior point at block 1 and q2 at block 2 as follows: q = (qiSi + q2Si)/(si + &z) (6.28) where s^ and Sj are the distances of the interior points to the boundary point. In other words, the gradient of q in the local one-dimension direction s is assmned constant, i.e.. It is noted that the conservative properties are not conserved in this method and errors may be introduced if the solution variations are large near the interfaces. However, the method is quite efficient and easy to implement and has been tested on a projectile flow problem with satisfactory results. Numerical experiments showed that only negligible errors were introduced if the grid was smooth and orthogonal near the interfaces. The interfaces of the six-block flow domain are identified by the block diagrams shovm in Figure 6.6. Here, some special lines that require special treatment are discussed. As shown in Figure 6.6a, the line j=28, k=13 of Block 1 is a common line among Block 1, 2, 3, 4 and 5, the boundary conditions at this line are obtained by interpolating over all these blocks, i.e., for 1 = 1, 60 qi28434 = q2i,i24 = q3i,ij = q44o,i4 = ^^iw = (Slj27q32,u + S3j2ql27,i34 + ^^i^
PAGE 113

107 Block 1 j=jinax k=l, 13 k=13 k=13,25 I k=l, 13 k=13 k=l k=l k=13 k=i: M j=l j=l j=l j=jmax j=jmax j=jmax Block 2 Block 3 Block 4 Block 5 j=jmax j=jmax j=jmax j=l j=l j=l k=l, 13 k=13 k=l k=l k=13 k=13 , 1 k=l,l3 k=13 k=13,25 Block 6 j=l (a) Block 1 & 2,5 and Block 6 & 2,5 Block 2 Block 3 Block 4 Block 5 k=kmax j=l, jmax k=l j=l, jmax j=l k=l,kmax j=jiiiax k=l , kmax k=l j=l, jmax k=kmax j=l, jmax k=kmax j=l,kmax k=kmax j=jmax, 1 j=jmax k=l , kmax j=l k=l,kmax (b) Block 2&3, 3&4and4&5 Figure 6.6 Interface block diagram

PAGE 114

108 Here, ql28,u4 stands for the unknown q at grid point (28,13,1) of block 1 and slj27 stands for the distance between grid point (27,13,1) and (28,13,1) in block 1, and likewise the other notation can be defined. A similar formula is applied to the line j = 1, k= 13 of Block 6 which is a conmion line of Block 2, 3, 4, 5 and 6. As can be seen from the grid network, the three surfaces j = l, k=35 and j=40 of Block 3, which are also j= 40, k=35, j = l surfaces of Block 4, form almost a flat plane. The boundary conditions at the lines between the surfaces, i.e., j = l, k=35 and j=40, k=35 of Block 3, are interpolated not only from block 3 and 4 but also from the neighboring surfaces, i.e., for 1 = 1, 60
PAGE 115

1 109 ^ "'^'^ ' A,[amJlKS')(lmax-l)l" <"2' Here, jmax, kmax, and Imax are the number of grid points in each curvilinear directions. 6.5 Pressure and shear stress coefficient The pressure and shear stress coefficients which are defined as _ _ r _ r'/Re are computed as follows: (a) pressure coefficient p* is first calculated from Eq.(2.9) and p*, = l/y, Cp is then readily computed from Eq.(6.33). (b) shear stress coefficient The shear stress tensor in Cartesian coordinate system can be obtained firom Eq.(2.2b). Since only the shear stresses on the solid surface are of interest, the shear stress tensor is then expressed in the curvilinear coordinate system, i.e., [r(^,n,^)] = [R][r][R]-» (634) where [R] is the transformation matrix that can be expressed in terms of the normaUzed metrics coefficients. The shear stress components t^^ and t^^ are r,, = (Ae,+Bey+CO/(^,'+c/+e/)^ (6.35) r,, = (Ar,, + Bv,+ CrjJ/(r,,' + v/+V.Y (6,35a)

PAGE 116

110 where B = (r^^ + r^y^ + r^,)/(x,2+y^'+z^')^ C = (',^f + ryj, + r^,)/(x/+y,^+z,^)^ (6.35c) The shear stress coefficients are then computed from Eq.(6.35a), i.e., Cffj = 2r5.j/ReM„ Cf,^ = Ir'yReM^ • (6.36) The drag and Uft forces can be obtained by integration of the surface pressure and shear stress distributions over the body surface.

PAGE 117

CHAPTER VII RESULTS AND DISCUSSION The proposed multi-block computational method is investigated for the simulation of transonic turbulent flows of Mach 0.8 past a wing-fuselage configuration by using the six-block grid network. Since the only available experimental data are the surface pressure coefficient at eight wing cross-sections, the discussion will concentrate on the results computed at these locations. These eight wing cross-sections, expressed in normalized distances from the symmetry plane, are, respectively, 0.15, 0.225, 0.305, 0.405, 0.5, 0.7, 0.85, and 0.95, in which the distance from the symmetry plane to wing tip is one. The multi-block thin-layer Navier-Stokes code (WF6B) provided the options to compute the Jacobian and metric coefficient by either finite-difference approach or volume/cell-face area approach. The use of finite-difference approach for computing Jacobians of transformation results in negative values at the certain regions near wing tip, the leading edge, and trailing edge. This seems to indicate that the grid network is quite skew near these regions. Hence, in order to accurately compute the geometric quantities, it is necessary to use the volume approach for Jacobian calculation. The metric coefficients can be computed by either approach. It is known that the use of the cell-face area approach for computing metric coefficient is more consistent with the volume 111

PAGE 118

112 approach for Jacobian calculation; however, the computed pressure coefficients are found to be insensitive to either method. 7.1 The Result from Base Grid The computational method was first investigated for a transonic mrbulent flow of Mach 0.8 past the wing-fuselage configuration at zero angle of attack with the use of the coarse base grid. The grid sizes of block 1 to block 6 are 28x25x60, 40x13x60, 40x35x60, 40x35x60, 40x13x60, and 26x25x60, respectively. A converged solution was obtained in 2500 time steps. Figure 7.1 shows the comparison of the computed Cp-distribution and the experimental data at four wing cross-sections. The result shows that the calculated wing surface pressure coefficient agrees well qualitatively with the experimental data but the quantitative agreement is not satisfactory at all at certain locations. One can " clearly observe the large discrepancy in the leading edge area and around the shock wave and high pressure gradient regions. In order to investigate the effect of the leading edge grid resolution upon solution accuracy, the base grid was modified by adding 10 more grid points in the regions close to the leading edge of the wing. This revised grid system was then used for a computation that restarted from the afready obtained solution. The computed Cp-distribution obtained in 500 time steps is shown in Figure 7.2, and does not show any improvement. This seems to indicate that around the leading edge, the grid spacing is not the major cause of the discrepancy. It is suspected that the highly skewed grid in that region as well as the turbulence model may be responsible for the disagreement between the prediction and the measurement.

PAGE 119

Figure 7.1 Cp-plot for a=0° based on the base grid

PAGE 120

114 Figure 7.2 Cp-plot based on leading edge refined grid

PAGE 121

' 1 115 In order to gain further insight and understanding on the solution algorithm, computation was made for flow at angle of attack a =-3° by using the base grid. The solution obtained from the case of zero angle of attack was used as an initial guess for restarting the computation. The Cp-distribution obtained after 400 time steps is shown in Figure 7.3. As shown in the figure, the Cpdistribution on the upper wing surface agrees well with the measured data; however, the agreement is not acceptable on the leeward side of the wing. This seems to suggest that the grid resolution is not fine enough to resolve the flow phenomena involving shock waves and/or large solution variations. Hence, it is demanded for a grid refinement study. 7.2 The Result from Refined Grid Limited by computer resources and accessibility, only the grid resolution in the wing chordwise direction was improved. A refined grid was generated by using natural cubic spline to double the number of grid points in the wing chordwise direction. With the use of the refined grid, flow at zero angle of attack was again computed. A converged solution was obtained in 500 time steps after restarting from the result of the base grid. Figure 7.4 presents the calculated Cpdistribution and the experimental data at eight different cross-section along the wing span. As can be seen, the computed surface pressure coefficient indeed agrees more closely to the measured data now, particularly in the shock and high gradient regions. These results clearly indicate that the accuracy of numerical sunulation can be further improved if a further improved grid resolution or a solution-adaptive grid is used in the simulation. Figure 7.5 shows the pressure

PAGE 123

Figure 7.4 Cp-plot for 0=0° based on refined grid

PAGE 124

Figure 7.4 (continued)

PAGE 125

119 (b) Top and side views; Back fuselage Figure 7.5 Cp-contour for a=0"

PAGE 126

120 (c) Top and side views; Block 5 fuselage (d) Top and side views; Block 2 fuselage Figure 7.5 (continued)

PAGE 127

121 0 (e) Upper wing surface (f) Lower wing surface Figure 7.5 (continued)

PAGE 128

122 coefficient contour on the surface of the wing-fuselage configuration with an interval of ACp=0.05. A Mach contour plot on a surface 1=20, a short distance above the wing surface, is given in Figure 7.6, which provides the information of the location, extent and strength of shock wave. The distributions of calculated shear stress components, t^^ and t^^ on the wing surface, as given in Figures 7.77.9, show that there is no reversed flow in the streamwise direction on the wing and that the fluid particles sUghtly above the wing surface are moving away from the wing root toward the wing tip except that some of those near the under surface (block 3) trailing edge move toward the wing root as they are swept downstream. It was observed that the flow characteristics represented by the shear stress distribution are in agreement with those exhibited by the surface pressure distribution. For instance, the variation of r^^ shown in Figure 7.7 clearly indicates a sudden decrease in fluid velocity across the shock. To further investigate the method, the refined grid was then used for the simulation of flows at 4" angle of attack. The result obtained from the base grid was again used as the initial guess for the computation. Figure 7.10 presents the Cp-distribution obtained in 500 time steps. The computed Cp on the windward side is in good agreement with measure data; however, the agreement on the leeward is not favorable. The predicted shock location on the upper surface is downstream of the measured location. It is also noticed that the computed Cp has shown oscillations near the leading and trailing edges of the wing. It is argued that the oscillations could be caused by the improper use of the cubic spline technique in refining the grids, since they appear only in the results

PAGE 129

123 ; * V (a) Upper wing, block 4 (b) Lower wing, block 3 Figure 7.6 Mach-contour at 1=20

PAGE 130

Figure 7.7 r^^-plot for 0=0" based on refined grid

PAGE 131

125 (b) Lower wing surface Figure 7.8 Surface shear stress component r

PAGE 132

I 126 (a) Upper wing surface (b) Lower wing surface Figure 7.9 Surface shear stress component r

PAGE 133

1 I t 127 Figure 7.10 Cp-plot for 0=4"

PAGE 134

128 obtained from the refined grid. The use of cubic spline may result in more skewed grid from the original skew base grid. The calculated shear stress distribution is shown in Figure 7.11, indicating that there is shock induced separation. Finally, some information regarding the CPU time requirement for the problems studied is given to help illustrate the scope as well as the complexity of the present methodology. All the computations were performed at the supercomputer, Cray-2 of NAS system at NASA Ames Research Center. The computer code WF6B has the efficiency of about 7.7x10'^ second per grid point per time step. The base grid has totally 311,400 grid points, and requires about 24 second for one time step. The result for zero angle of attack were computed for 2,500 time steps, that is, about 17 hours of CPU time. The refined grid has 536,040 grid points and requires 41 second for one time step. The result obtained was computed for 500 time steps after restarting from the result of the base grid. The CPU time is about 5.5 hours for 500 time steps.

PAGE 135

Figure 7.11 Shear stress coefficient plot for 0=4"

PAGE 136

CHAPTER Vm CONCLUDING REMARKS A multi-block computational method is developed for three dimensional high speed turbulent flows over complex configurations. In this method, the flow field is divided into several contiguous blocks such that each block is governed by the same set of thin-layer Navier-Stokes equations. The solution algorithm with distance weighted interpolation for updating the interface boundary condition is rather simple and straightforward; however, special measures in the grid topology and grid generation are required before a satisfactory prediction can be made. A flow simulation package which includes grid generation, flow solver, and plotting program, is developed and apphed to a transonic turbulent flow of Mach 0.8 past a wing-fuselage configuration. The flow solver is good for high Reynolds number transonic turbulent flow problems with modest flow separations. A special designed six-block grid topology is used for the development of the method. The flow field is properly divided into six blocks that are solved by the same unsteady TVD thin-layer Navier-Stokes code. A coarse base grid and a refined grid are generated and used for flow simulations. Numerical experiments performed in this study have shown that special measures and care are required in generating good grid system involving excessive distortion between the physical and computational domain; however, the ' .A " V 130

PAGE 137

131 computed shock location and pressure distribution can be in good agreement with the experiment data if the block grids used are adaptive to important solution characteristics. The proposed method has demonstrated its flexibility, efficiency, and capability in the simtilation of complex geometry flow problems. In the demand of the CFD methods for simulating turbulent flows over complex geometry, the proposed method is indeed a very promising approach to be further developed. As experienced in this study, a good grid network is crucial to achieve accurate results. To develop a cost-effective computational method for 3-D complex configiiration aerodynamics, a technique to generate solution-adaptive grid network is desired. However, a good grid that adapts to the important flow characteristics of the present complex flows are not known a priori; hence, it is very desirable to incorporate a self-adaptive grid generation technique that can give a proper grid distribution for yielding accurate solutions with minimum number of grid points to minimize the computing cost.

PAGE 138

APPENDIX A JACOBIAN MATRICES The Jacobian matrices A, B or C equals k. kx k, 0 ' k,^-u-(cf-\)u6 ky^-(y-l)ve kji,-(-r-l)we k^+js . where i> = 7(e/p)-<^ = (7-l)(u^+v^+w^)/2 e = + ky + with k = ^, T7 or r for A, B or C, respectively. The viscous flux Jacobian M is 0 0 0 0 0 «l«f(Vp) a2«f(l/p) «3«f(Vp) 0 "2«f(l/p) «4«f(Vp) 0 «35f(l/p) «5«f(l/p) c^eScii/p) 0 . ^51 mjj a„5f(l/p) , where ni2i = ^(-u/p) + a^S j.(-v/p) + ajSf (-w/p) 132

PAGE 139

133 = a2«f(-u//') + a45f(-v/p) + a5«f(-w/p) ni4i = a3«f(-u/p) + Q5Sf(-v//') + a6«f(-w/p) nisi = ai6f(-uVp)+a25f(-2uv/p)+«35f(-2uw/p) + a^S f (V/p ) + aj6 ^( V/p) + aj* ^ (-2vw/p) °l52 = -ni2i-ao*f(u/p) = -m4i-ao5f(w/p) ao = 7MPr-'(r/+ry'+rz') «i = 4(4/3)r/+fy'+r/] "2 = W3)r.ry «4 = M[f,'+(4/3)r/+f,'] "5 = (M/3)ryrz «6 = /^[r,'+f/+(4/3)r,']

PAGE 140

APPENDIX B : EIGEN-FUNCnONS ; The eigenvalues of A, B and C are {a/,a/,a/,a,^a/} = {U,U,\J,lJ+ci^,'+i/+^,y,V-c(i,'+i/+Cf} {a,^a,^a/,a,^a/} = {V,V,V,V^-c(,,^+,/+,,^)^V-c(,,2+,/+,,^)^} {a/,a,^a,^a,^a,^} = {W,W,W,W+c(f,2+fy2+r/)^W-c(r,'+r/+rz')^} respectively, where c = (yp/p)^ is the local speed of sound and the contravariant velocity components U, V, W are U = C,+l,u + CyV+^^ * ' ^ ^ 'V . V = i7t + >7xU + »7yV+i7,W "t'? w = c,+rxU+r,v+f^ • ^rv. > The similarity transformation for the Jacobian matrices A, B, C are A = R^,R, i , B = R^A^R^ i , C = R,A,R,-^ where the a's are the diagonal matrices A, = D[u,u,u,u+c(^/+e/+e,2)^u-cU/+?/+?,2)«] :., y v A, = D[V,V,V,V+c(„^+„/+„/)^V-c(„/+„/+„,2)V4j A, = D[w,w,w,w+c(^/+f/+^,')^w-c(f/+r/+f/)^] If one defines the notation as follows and k = ^, »; or f * 134 '

PAGE 141

135 e = ICjU + ICyV + lCjW and a = J2p/2c, ^ = /Z/2pC, 7 = 7-1 0 = (7-l)(u^+v2+w2)/2 then = a a a(u + lc,c) q(u-ICjC) lc,v-Jc,p a(v + lCyC) a(v-ICyC) a(w + I^c) a(w-I^c) p(k,w-I^u) p(RyU-Ic^v) c[( + <^)/j + ce] a[(^ + c')/7 -ce] and R^-* = ^(i-^/<): ^S'^u/c^-ic./p ^7v/c^ i^7w/c^+E,/p (lc,w-R,u)/p R,(lWc^)^yu/<^+yp 1c,7v/c^-Vp 1^7w/c^ -1^7/c' ^(^-Ce) ^[Ic,C-7u] ^[lSC-7V] ^[1^C-7W] /37 ^P( + Ce) -p\^C + yu] -^[EyC + 7V] -/3[I^c + 7w] ^7

PAGE 142

APPENDIX C CLUSTERING FUNCTIONS The algebraic grid generation requires some sort of clustering functions or stretching functions to distribute the grid points so that more grid points are clustered where the solution varies rapidly, e.g. near the soUd surface. Four clustering functions have been implemented in GG3D and are discussed as follows: 1. CUBICLU Cubic Function Distribution For a space curve in the 1-direction, for instance, the normalized arc length s is defined, in terms of arc length S(l), as S(l)-S(l) ^® = S(lma.)-S(l) (^1' Let the specified first spacing and the last spacing be denoted by asj and AS2, respectively. The cubic function distribution is obtained as follows. First, define h = L/(hnax-l) Then a, b, c are computed as a = [asi + ASj -2h] b = [ASi h c(h^-h)]/(h^-h) (C2) c = 1 b c The normalized cubic function distribution, v(l), is then given by v(l) = a(l-l) +b(l-l)^ + c(l-l)' (C3)

PAGE 143

137 2. TANHCLU Hyperbolic Tangent Distribution Using the same notation as above, the hyperboUc tzingent distribution is constructed as follows. First, compute A and B A = Vasj/Vasi B = (lmax-l)yAS2ASi (C.4) Then the following nonlinear equation is solved for 6 by the Newton-Raphson method. sinh(5)/5 = 1/B (C.5) The normalized hyperbolic tangent distribution, v(l), is then given by = Arfi-A)u(i) where urn = -a+ tanh[5{(l-l)/(hnax-l)-0.5}] 2 ^ tanh(fi/2) ^ ^^''^ 3. SINHCLU Hvperbolic Sine Distrhutinn If only the first spacing asi is specified, the hyperbolic sine distribution is calculated as follows. First, B is computed B = (hnax-l)ASi (C.8) and 6 is determined by the Newton-Raphson method as the solution of sinh(5)/5 = 1/B , The normalized hyperbolic sine distribution, v(L), is given by ?

PAGE 144

138 smh[s{(l-l)/(lmax-l)}] ^ sinh(6) 4. EXPOCLU Exponential Distribution ' ^ If only the first spacing, asj, is specified, the exponential distribution is calculated as follows. First, 5 is obtained by the Newton-Raphson method fi"om the following nonlinear equation: 1 = ASi[(l + 5)'™"-*-l]/« ' '': (CIO) Then the normalized exponential distribution is given by . ^ . i * _ , i> V(l) = V(l-l) + ASi(l + 5)" V . (C.11) . , ^ — : " ^ . 5. Summary The truncation error is strongly affected by the point distribution, and studies of the effect of different distribution functions on the truncation error have been made in the past. The exponential distribution can cluster more grid points near the 1=1 surface. However, the variation of the spacing is large. The cubic function is good for clustering points near both ends. The hyperbolic sine function gives a smoother distribution in the immediate vicinity of the 1 = 1 surface, while the hyperboHc tangent function gives a smoother overall distribution.

PAGE 145

REFERENCES 1. Goldhammer, M. I. and Rubbert, P. E., "CFD in Design: An Airframe Perspective," AIAA 89-0092, AIAA 27th Aerospace Science Meeting, Reno, Nevada, Jan. 1989. 2. Schiff, L. B., Cummings, R. M., Sorenson, R. L. and Rizk, Y. M., "Numerical Simulation of high-Incidence Flow over the F-18 Fuselage forebody," AIAA 89-0339, AIAA 27th Aerospace Sciences Meeting, Reno, Nevada, Jan. 1989. 3. Abid, R., Vatsa, V. N., Johnson, D. A. and Wedan, B. W., "Prediction of Separated Transonic Wing Flow with Nonequilibriimi Algebraic Turbulence Model," AIAA Journal, Vol. 28, No. 8, Aug. 1990, pp. 14261431. 4. Jameson, A. and Yoon, S., "Lower-Upper Imphcit Schemes with Multiple Grids for the Euler Equations," AIAA Journal, Vol. 25, July 1987, pp. 929935. 5. Rai, M. M., "A Relaxation Approach to Patched-Grid Calculations with the Euler Equations," Journal of Computational Physics, Vol. 66, No. 1, 1986, pp. 99-113. 6. Gatlin, B. and Whitfield, D. L., "An Implicit, Upwind, FiniteVolume Method for Solving the 3-D Thin-Layer Navier-Stokes Equations," AIAA 87-1149, AIAA 8th Computational Fluid Dynamics Conference, Honolulu, Hawaii, June 1987, pp. 516-526. 7. lijewski, L. E., 'Transonic Euler Solutions on Mutually interfering Finned Bodies," AIAA Journal, Vol. 28, No. 6, June 1990, pp. 982-988. 8. Hoist, T. L., Kaynak, U., Gundy, K. L., Thomas, S. D., Flores, J. and Chaderjian, N., "Transonic Wing Flows Using an Euler/Navier-Stokes Zonal Approach," Journal of Aircraft, Vol. 24, No. 1, Jan. 1987, pp. 17-24. 9. Srinivasan, G. R., McCroskey, W. J., Baeder, J. D. and Edwards, T. A., "Numerical Simulation of Tip Vortices of Wings in Subsonic and Transonic Hows," AL\A Journal, Vol. 26, No. 10, Oct. 1988, pp. 1153-1162. 139

PAGE 146

140 10. Flores, J. and Hoist, T. L., Numerical Solution of the Navier-Stokes Equations for Complex Configurations. Lecture Note, Univ. of Tenn. Space Institute, Mar. 1988. 11. Cummings, R. M., Rizk, Y. M., Schiff, L. B. and Chaderjian, N. M., "Navier-Stokes Predictions of the Bowfield Around F-18 (HARV) Wing and Fuselage at Large Incidence," AIAA 90-0099, 28th Aerospace Sciences Meeting, Jan. 1990. 12. Baldwin, B. S. and Lomax, H. "Thin Layer Approximation and Algebraic Model for Separated Turbulent Flows," AIAA 78-257, AIAA 16th Aerospace Sciences Meeting, Huntsville, Ala., Jan. 1978. 13. Malvern, L. E., Introduction to the Mechanics of a Continuous Medium. Prentice-Hall Inc., Englewood Cliffs, New Jersey, 1969, pp. 423-433. 14. Yih, C. S., Fluid Mechanics, West River Press, Ann Arbor, Michigan, 1988, pp. 40-41. 15. Schlichting, H., Boundary Layler Theory, 7th edition, McGraw-Hill Book Co., New York, 1978. 16. Anderson, D. A., Tannehill, J. C. and Fletcher, R. H., Computational Fluid Mechanics and Heat Transfer, McGraw-Hill Book Co., New York, 1984, pp. 479-502. 17. Tennekes, H. and Lumley, J. L., A First Course in Turbulence, 11th printing. The MIT Press, Cambridge, Massachusetts, 1987. 18. Reynolds, W. C. and Cebeci, T., "Calculation of Turbulent Flows," Topics in Applied Physics, Turbulence, Vol. 12, edited by P. Bradshaw, SpringerVerlag, New York, 1978, pp. 193-229. 19. Cebcei, T. and Smith, A. M. O., Analysis of Turbulent Boundary Layer, Academic Press, New York, 1974. 20. Lakshminarayana, B., 'Turbulence Modelling for Complex Shear Flows," AIAA Journal, Vol. 24, No. 12, Dec. 1986, pp. 1900-1917. 21. Beam, R. M. and Warming, R. F. "An Imphcit Finite-Difference Algorithm for Hyperbohc Systems in Conservation-Law Form," Journal of Computational Physics, Vol. 22, Sep. 1976, pp. 87-110. 22. PuUiam, T. H., "Artificial Dissipation Models for the Euler Equations " AIAA Journal, Vol. 24, Dec. 1986, pp. 1931-1940.

PAGE 147

141 23. Shiau, N. H. and Hsu, C. C, "A Diagonalized TVD Scheme for Turbulent Transonic Projectile Aerodynamics Computation," AIAA 88-0217, AIAA 26th Aerospace Science Meeting, Reno, Nevada, Jan. 1988. 24. Chen, M. H., "On Total Variation Diminishing Schemes for Transonic Turbulent Flow Computation," Ph.D. Dissertation, University of Florida, 1988. 25. Hsu, C. C. and Shiau, N. H., "Numerical Simulation of Turbulent Transonic Projectile Aerodynamics," AIAA 87-1237, AIAA 19th Fluid Dynamics, Plasma Dynamics and Lasers Conference, June 1987. 26. Shiau, N. H., Hsu, C. C. and Chyu, W. J., "Numerical Simulation of 3-D Transonic Turbulent Projectile Aerodynamics by TVD schemes," AIAA 890335, AIAA 27th Aerospace Sciences Meeting, Jan. 1989. 27. Yee, H. C, "Construction of Exphcit and ImpUcit Symmetric TVD Schemes and Their AppUcations," Journal of Computational Physics, Vol. 68, 1987, pp. 151-179. 28. Beam, R. M. and Warming, R. F., "An ImpUcit Factored Scheme for Compressible Navier-Stokes Equations." AIAA Journal, Vol. 16, 1978, pp. 393-401. 29. Steger, J. L., "ImpUcit Finite-Difference Simulation of Flow about Arbitrary Two-dimensional Geometries," AIAA Journal, Vol. 16, July 1978, pp. 679-686. 30. PuUiam, T. H. and Steger, J. L., "Recent Improvements in Efficiency, Accuracy, and Convergence for ImpUcit Approximate Factorization Algorithm," AIAA 85-0369, AIAA 23th Aerospace Sciences Meeting, Jan. 1985. 31. Harten, A., "A High Resolution Scheme for the Computation of Weak Solutions of HyperboUc Conservation Laws." Journal of Computational Physics, Vol. 49, 1983, pp. 357-393. 32. Sod, C. A., Numerical Methods in Fluid Dynamics, Cambridge University Press, New York, 1985, pp. 300-355. 33. Roe, P. L., "GeneraUzed Formulation of TVD Lax-Wendroff Scheme," ICASE Report No. 84-53, Oct. 1984. 34. Davis, S. F., "TVD Finite Difference Schemes and Artificial Viscosity " ICASE Report No. 84-20, 1984.

PAGE 148

142 35. Sweby, P. K,, "High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws," SIAM J. Numer. Anal., Vol. 21, 1984, pp. 995-1011. V. 36. Harten, A. and Hyman, J. M., "A Self-Adjusting Grid for the Computation of Weak Solutions of Hyperbolic Conservation Laws," Journal of Computational Physics, Vol. 50, 1983, pp. 235-269. 37. Yee, H. C, "Linearized Form of Implicit TVD Schemes for the Multidimensional Euler and Navier-Stokes Equations," Computation and Mathematics with AppUcations. Vol. 12A. Nov. 1986, pp. 413-432. 38. Yee, H. C, "On Symmetric and Upwind TVD Schemes," Proceedings of the 6th GAMM Conference on Numerical Methods in Fluid Mechanics, Sep. 1985, pp. 399-407. 39. Yee, H. C. and Harten, A., "ImpUcit TVD Schemes for Hyperbolic Conservation Laws in Curvilinear Coordinates," ALAA 85-15 13-CP, AIAA 7th Computational Fluid Dynamics Conference, Cinn,, Ohio, July 1985. 40. Nietubicz, C. J., Pulliam, T. H. and Steger, J. L., "Numerical Solution of the Azknuthal-Invariant Thin-Layer Navier-Stokes Equations," ALAA 790010, AL\A 17th Aerospace Sciences Meeting, Jan. 1979. 41. Hsu, C. C, Shiau, N. H. and Reed, C. W., "Numerical Simulation of Transonic Turbulent Flow Past a Real Projectile," AJAA 88-0218, AL\A 26th Aerospace Sciences Meeting, Jan. 1988. 42. Facey, J. R., "Return of the Turboprops," Aerospace America, Vol. 26, Oct. 1988, pp. 14-23. 43. Yang, S. C. and Hsu, C. C, A Multi-Block Flow Simulation Package Reference Manual, Department of Aerospace Engineering, Mechanics, and Engineering Sciences, University of Rorida, 1990. 44. Thompson, J. F. (Ed.), Numerical Grid Generation, North-Holland, New York, 1982. 45. Eiseman, P. R., "Grid Generation for Fluid Mechanics Computations," Annual Review of Fluid Mechanics, Vol. 17, 1985, pp. 487-522. 46. Yang, S-L. and Shih, T. I-P., "An Algebraic Grid Generation Technique for TimeVarying Two-Dimensional Spatial Domains," International Journal for Numerical Methods in Fluid, Vol 6, 1986, pp. 291-304. 47. Shih, T. I-P., Finite-Difference Methods in Computational Fluid Dynamics to be published, pp. 313-407.

PAGE 149

143 48. Smith, R. E. and Eriksson, L-E., "Algebraic Grid Generation," Computer Methods in AppUed Mechanics & Engineering, Vol. 64, 1987, pp. 285-300. 49. Thompson, J. F., "A General Three-Dimensional Elliptic Grid Generation System on A Composite Block Structure," Computer Methods in Applied Mechanics & Engineering, Vol. 64, 1987, pp. 377-411. 50. Thompson, J. F., Warsi, Z.U.A, and Mastin, C. W., Numerical Grid Generation Foundations and Applications, North-Holland, New York, 1985, pp. 188-330. 51. Thompson, J. F. and Gatlin, B., "Program EAGLE User's Manual, Vol. II and m. Surface and Grid Generation Code," AFATL-TR-88-117, 1988. 52. Berger, M., "On Conservation at Grid Interfaces," SIAM Journal of Numerical Analysis, Vol. 24, No. 5, Oct. 1987, pp. 967-983. 53. Moon, Y. J. and Liou, M. S., "Conservative Treatment of Boundary Interfaces for overlaid Grids and Multi-lever Grid Adaptation," AlAA 891980, AIAA 9th Computational Fluid Dynamics Conference, Buffalo, New York, June 1989, pp. 480-487.

PAGE 150

BIOGRAPHICAL SKETCH Shu-Cheng Yang was born in 1954, at Kimmen, the Republic of China. In 1972, he went to the Chung-Cheng Institute of Technology, Taiwan, where he received his bachelor's in mechanical engineering in 1976. In 1981, he earned his M.S. in mechanical engineering from the University of Cincinnati. In 1986, he enrolled as a graduate student in the Department of Aerospace Engineering, Mechanics and Engineering Sciences at the University of Florida.

PAGE 151

qu^ty, a. a dissenation for the'd^p\ToTDoa^'of%MLo^^^^^^^ " '"^ Chen-Chi Hsu, Chairperson Professor of Aerospace Engiiieering, Mechanics, and Engineering Sciences quahty, as a dissertation for the degree of Doctor of Philosophy ^ Bernard M. Leadon Professor of Aerospace Engineering, Mechanics, and Engineering Sciences accemable^ndti H '"^^^ ^ °Pi^o° " confonns to acceptable standards of scholarly presentation and is fully adequate, in scope and quahty, as a dissertation for the degree of Doctor of PhHosophy ^ Uhich H. Kurzweg Professor of Aerospace Engineering, Mechanics, and Engineering Sciences

PAGE 152

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Wei Shyy Associate Professor of Aerospace Engineering, Mechanics, and Engineering Sciences I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Arun KVarma Professor of Mathematics This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. December 1990 , Winfred M. PhiUips 7^ Dean, College of Engineering Madelyn M. Lockhart Dean, Graduate School


xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EFYZGXKFA_WLDDON INGEST_TIME 2015-04-14T19:00:34Z PACKAGE AA00029914_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES