Numerical simulation of the aerodynamics of turbulent transonic flow past a projectile


Material Information

Numerical simulation of the aerodynamics of turbulent transonic flow past a projectile
Physical Description:
vi, 134 leaves : ill. ; 28 cm.
Shiau, Nae-Haur Eric, 1958-
Publication Date:


Subjects / Keywords:
Aerodynamics, Transonic   ( lcsh )
Projectiles -- Aerodynamics   ( lcsh )
Numerical calculations   ( lcsh )
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )


Thesis (Ph. D.)--University of Florida, 1987.
Includes bibliographical references (leaves 132-133).
Statement of Responsibility:
by Nae-Haur Eric Shiau.
General Note:
General Note:

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 001061709
notis - AFE5622
oclc - 18763017
System ID:

This item is only available as the following downloads:

Full Text








The author wishes to express his sincere gratitude for

Dr. Chen-Chi Hsu, chairman of the committee, for his

guidance, support and encouragement during the time of


Also, his many thanks are extended to the committee

members, Dr. Malvern, Dr. Leadon, Dr. Kurzweg and Dr. Shih

for their interest and time.

The calculation in this study was performed by a Cray

system at Bell Laboratories and by a Cray X-MP/48 at

Pittsburgh Supercomputing Center with all CPU time provided

by National Science Foundation through Grant ECS-8515761.


ACKNOWLEDGMENTS ........................................

ABSTRACT ...............................................





I INTRODUCTION ...................................

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

Navier-Stokes Equations ........................
Transformed Navier-Stokes Equations ............
Thin-Layer Navier-Stokes Equations .............
Axisymmetric Thin-Layer Navier-Stokes
Equations ....................................

III NUMERICAL ALGORITHM ............................

Time-Differencing ...................... .......
Approximate Factorization ......................
Numerical Dissipation ..........................
Finite Difference Equations ....................
Finite Difference Scheme for Axisymmetric
Thin-Layer Navier-Stokes Equations ...........


Flow Domain ....................................
Boundary Conditions .............................
Initial Condition ..............................
Viscosity ...... ...............................
Forces and Moments .............................
Rate of Convergence ............................

V HYPERBOLIC GRID GENERATION .....................

Governing Equations ............................
Finite Difference Form .........................

VI FLOW PROBLEM ...................................



Inviscid Flow .................................. 59
Flow Domain Determination ..................... 59
Effect of Artificial Dissipations ........... 60
Viscous Flow .................................. 62
Effect of Artificial Dissipations .......... 65
Effect of Eddy Viscosity Distribution ....... 67
Effect of Grid Resolution in Streamwise
Direction ......... ........... ............. 67
Effect of Grid Resolution in Normal Direction 71
Importance of Viscous Drag to Pressure Drag 76


Flow Cases at 100 Angle of Attack .............. 90
Flow case at M =0.96 ....... ................. 90
Flow case at M =0.91 and 1.20 ...............103
Flow Case at 450 Angle of Attack ...............110
Comparison of Flow Cases at Different Angles of
Attack ........................................126

IX CONCLUDING REMARKS ..............................128

APPENDIX ................................................ 130

REFERENCES ................................................ 132

BIOGRAPHICAL SKETCH ..................................... 135

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



Nae-Haur Eric Shiau

December 1987

Chairman: Chen-Chi Hsu
Major Department: Engineering Sciences

Steady, transonic flow over a projectile was numerically

investigated to gain insight and understanding of the

thin-layer Navier-Stokes approximation and to assess the

relative importance of viscous effects on the aerodynamic

forces acting on a projectile at transonic speeds. The

turbulence model used in the computation is an algebraic

two-layer eddy viscosity turbulence model proposed by

Baldwin and Lomax. The projectile model considered was a

secant-ogive-cylinder-boattail configuration for which

accurate surface pressure measurements are available to

assess the numerical results. The sensitivity of the

numerical scheme to the grid resolution and added numerical

dissipation terms was investigated. The surface shear stress

distribution has also been computed to aid in identifying

regions of the reversed flow.

The computed results, based on a grid system generated by

a set of hyperbolic partial differential equations, were

verified by comparing with experimental measurements. Good

agreement between the numerical results and measured data

indicates that the numerical scheme can predict the

essential features of the flow field. The calculated

results obtained at zero angle of attack and Mach number

0.91 show that the viscous drag is as important as the

pressure drag (no base drag). However, the relative

importance of viscous drag decreases with increasing Mach

number and angle of attack. At 10 and 45 degrees angle of

attack, the viscous forces had a negligible effect on both

the lift force and the pitching moment.


An accurate prediction of the aerodynamic forces is

important in the design of flight vehicles. For complex flow

fields, an accurate simulation of the flow field is crucial

for the successful prediction of such forces. Considerable

experimental efforts have been undertaken to enhance our

understanding of these complex flows and, recently, with the

increase in computational capabilities, there has been a

major effort to numerically simulate such flows. The most

challenging features of transonic flow to be numerically

simulated are the viscous/inviscid and shock/boundary layer

interactions. There are two basic approaches used to

numerically simulate these flows. One method is a patching

technique in which the solution is obtained from inviscid

equations and boundary layer equations separately and then

patched together in some iterative process. Another

approach is to solve a suitable subset of the Navier-Stokes

equations for the entire flow domain. The primary advantage

of the patching method is their computational efficiency.

However, it cannot solve the separated flow. The latter

approach has the advantage that it provides a continuous

treatment of the solution through the interaction regions.

- 1 -

2 -

In the transonic flow regime, critical aerodynamic

behavior has been observed which can be attributed to the

complex shock pattern and local flow expansions that exist

over the projectile at these speeds. Shadowgraphs which can

be found in reference 7, show the shock pattern existing on

a typical projectile configuration for transonic flow. The

formation of shock waves imbedded in the flow field produces

severe changes in the aerodynamic coefficients such as the

drag force and the pitching moment. For example, the

nondimensional drag force for a typical projectile shape has

been found to change by as much as 100% through a Mach

number range of 0.95 to 0.97 [9]. Similar results have also

been evident in this study. For such large changes in the

magnitude of the aerodynamic forces, it is essential to

understand and be capable of computing the features of the

flow field that contribute to this phenomenon. In the

transonic flow regime, the flow field is characterized by

strong viscous/inviscid, shock/boundary layer interactions.

It is, therefore, advantageous to use the approach which

solves the suitable subset of Navier-Stokes equations since

this method considers these interactions in a fully coupled


A three dimensional code was developed in 1978 by Pulliam

and Steger at NASA Ames to solve the thin-layer

Navier-Stokes equations for an ideal gas in a transformed

boundary-fitted coordinate system [14]. In this code, which

is applicable to high speed compressible flow problems, the

3 -

governing equations are approximated by the spatially

factored implicit scheme of Beam and Warming [3,4]. To

control the numerical stability of the solution algorithm,

second order implicit and fourth order explicit dissipation

terms are added. Consequently it results in a block

tridiagonal system of equations to be solved at each time

level. The turbulence closure model implemented in the code

is an algebraic two-layer eddy viscosity model of Baldwin

and Lomax [2]. The unsteady thin-layer Navier-Stokes code

has an option for calculating inviscid flow problems and a

steady flow solution can be obtained as the converged

solution of a time asymptotic flow problem. This code has

been widely applied to various flow problems.

The three-dimensional thin-layer Navier-Stokes code has

also been simplified for axisymmetric flow cases. Both codes

have been investigated to some extent by the U.S. Army

Ballistic Research Laboratory on transonic projectile

aerodynamic problems [9,11]. The planar grid network

provided to a Navier-Stokes code is obtained from a grid

generation code GRIDGEN which can give either an elliptic

grid or a hyperbolic grid [10]. Moreover, a three

dimensional grid system is formed by a sequence of planar

grids around the axis of a projectile. For a

secant-ogive-cylinder-boattail projectile with sting at zero

angle of attack, the published result showed that the

computed surface pressure coefficient over the secant-ogive

portion and the boattail portion of the projectile agrees

4 -

rather well with the measured data but the agreement on the

cylinder portion of shock/boundary layer interaction region

is not very satisfactory for some cases considered. For the

projectile model at 2 degrees angle of attack the

Cp-distribution computed agrees qualitatively with the

measured data but quantitatively the agreement over the

cylinder portion and boattail portion is not satisfactory

[18]. The published results seem to indicate that the

thin-layer Navier-Stokes codes can give acceptably accurate

surface pressure for the complex transonic projectile

aerodynamic problem if a good adaptive grid system is

provided. However, more precise causes for the

unsatisfactory results reported are yet to be investigated;

moreover, no result on the skin friction coefficient Cf has

been reported and discussed.

The GRIDGEN code obtained from the U.S. Army has been

investigated and modified [6]. The modified hyperbolic grid

is generated using predetermined adaptive boundary points.

In order to gain the smoothness property, the hyperbolic

grid obtained from GRIDGEN also has been modified to bend

the grid lines in the nose region of the projectile, which

has been found to be important in the convergence process

and for the solution accuracy. The modified hyperbolic grid

generation code is efficient, requiring less computer time

than elliptic grid generation techniques; moreover, it can

provide a good grid network and is therefore used in this


5 -

An accurate prediction of a complex viscous flow field via

thin-layer Navier-Stokes codes depends upon many factors.

These include grid resolution, turbulence models, artificial

dissipation terms, flow characteristics, etc. For a

successful application and/or modification of a

Navier-Stokes code, one must be familiar with certain

relevant features of the code. This motivates the present

investigation of the computer codes. The main objectives of

this study are to further advance our understanding and gain

in-depth insight on the application of thin-layer

Navier-Stokes codes implemented with the Baldwin-Lomax

turbulence model for accurate computation of transonic

projectile aerodynamics.


Navier-Stokes Equations

In three dimensions, the unsteady compressible

Navier-Stokes equations referenced to Cartesian coordinates

without body forces and external heat addition can be

written in the following vector form:


q 3E 8F aG aEV aFV aG
- + + + -- + +
St ax Sy az ax Sy 8z

in which q is an unknown vector, E, F, G are the inviscid

flux vectors and EV F GV are the viscous flux vectors.

The explicit expressions for these vectors can be written as

p pu pv pw E ]T
pu pu +p puv puw

pv puv pv2+p pvw

pw puw pvw pw2+p

0 rxx 'xy 'xz x

0 yx Tyy Tyz y

0 tzx Tzy 'zz z

(E+p)u IT

(E+p)v ]

(E+p)w ]T



- 6 -

q = [

E = [


EG = [

FG= [
G =


7 -

Here, p is the density; u, v and w are the velocity

components in x, y and z-directions, respectively; p is the

pressure; and E is the total energy per unit volume which is

related to the internal energy e and the kinetic energy by

the formula

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

For a Newtonian fluid with zero bulk viscosity, the

relationships between stresses and velocity gradients are



TZ -- _

xy =

Tyz = 1

Sx =

2 ( 2--

li (2-

V (2--

au av

ay ax

av aw

az ay

aw au

ax az









) = zy

) XZ

and the B's in equations (2) are related to viscous

dissipation and heat conduction by

- 8 -

0x =Uxx + Vxy + xz + k -

y = Uyx + VTyy + Wyz + k (5)

z = UTz + VT + WT + k -

where p and k are the coefficients of dynamic viscosity and

thermal conductivity, respectively, and T is the


The complete set of time-dependent Navier-Stokes equations

(1) includes one continuity equation, three momentum

equations and one energy equation, but there are seven

unknowns p, u, v, w, E, p and T. In order to close the

system of equations, the fluid is assumed to be a perfect

gas. The equation of state and other relations for a

perfect gas are

p = pRT

e = cT h = T p(6)

where R is the gas constant, h is enthalpy, and cp and c

are the specific heat at constant pressure and volume,

respectively. I is the ratio of specific heats, which for

air at standard conditions is 1.4. Now, by using equations

9 -

(3) and (6), one can obtain two additional equations

relating p and T to the dependent variables q

p = (T 1)[E p (u2 + v + w2)]
(7-1) E 1
T [-- (u2 + v2 + w2)]
R p 2

Hence, equations (1) and (7) are the governing equations for

a class of unsteady compressible flow problems.

The number of parameters for flow problems of the same

class can be reduced by introducing dimensionless variables.

Once the equations are put into nondimensional form, the

solution can easily be used to correlate data for succinct

presentation using the minimum possible amount of

experimental testing. The dimensionless quantities used here

are obtained by choosing characteristic scales in the

following manner:

S x y z t
x y z t

S u v w T
u v w T (8)

E *
E p p =
V. P.. PV 2

10 -

where the nondimensional variables are denoted by an

asterisk and the freestream conditions by o. For instance,

V is a characteristic velocity of the fluid at freestream

and L is the characteristic length scale.

The substitution of the dimensionless quantities, equation

(8) into the governing equations (1) and (7), yields

+ + + -Re (- + + ) (9)
at ax 8y 8z ax ay az


*2 *2 *2
p =( 1)[E p (u + v + w)]
E 1
2 *2 *2 *2
T = (T-1) M [ (u2 + v + w )]
p 2

where the Reynolds number Re and the Mach number at

freestream are defined as

p V L V
Re M (11)

The unknown vector q and flux vectors E F G E ,
F G* have the same form of equation (2) as well as the

stresses in equations (4) except that all the quantities are

nondimensional quantities. The components of P become

- 11 -

= *
x Uxx

* ** *
xy "xz

B = U T + V T + W yz
zy yx yy yz

S =U t + v T + W
z Zx zy zz

+ -
(T-l)Pr Me 8x
(T-1)Pr M2 y*
+ ----M2 -z*
2 *
(r-1)Pr M az

Here, the Prandtl number Pr is

cP v
Pr p

The typical values of the Prandtl number for air at standard

conditions are 0.9 and 0.72 for turbulent and laminar flow,

respectively. In order to simplify the notation, the

asterisks are removed from the nondimensional equations

hereafter except where noted.

Transformed Navier-Stokes Equations

To enhance numerical simplicity, accuracy and efficiency,

a generalized boundary-fitted curvilinear coordinate system

is used. The boundary-fitted coordinate transformation maps

all boundary surfaces in the physical domain onto coordinate

surfaces of the computational domain. The boundary-fitted

curvilinear coordinates not only facilitate implementation

of boundary conditions accurately, but also allow all



12 -

computations to be performed on a fixed square grid in the

computational domain. The numerical solution technique,

thus, becomes more versatile.

Consider a generalized transformation of the form

S= ( ( x,y,z,t )

n = n ( x,y,z,t )

5 = 5 ( x,y,z,t )

= t

where (t,n,T,r) are the coordinates of the computational

domain and (x,y,z,t) are those of physical domain. It is

known, e.g. reference 1, that the Jacobian of the

transformation J-1, which represent the volume of the grid

cells in physical space, is given by

J- xty z + xqyCZt + x;yzn xtyZz xCyRzt (14)

x ytzc

and the metrics in the computational domain can be related

to those in the physical domain by using a chain rule

expansion of x,, yt, etc., and then solving for tx, ty'

etc., to give

- 13 -

tx = J ( y zC yzlI )
ty = J ( zqxC z x1 )

tz = J ( x Iy x Cy )
nx = J ( y cz Ytz )

y = J ( z x zx )

"z = J ( x yt xty )
Cx = J ( ySz ynz )

y = J ( zgx zxt )

Cz = J ( xy, xy )

tt = -xSx YjSy z Jz
qt = -Xqx YXIy zZ z

Ct = -X Ix YXCy z TC

Note that the subscript in the metrics denotes

differentiation. The metric identities of the

transformation, e.g. reference 1 or 20, are

a 1 a St a t a t
(-)+ (- ) +- (- ) +-) =0
aa J at J al J at j

a t a IX a C
x x x
( ) +- (- ) +-( ) =o
at J an J at J
a y. a y a y
( ) +- ( ) +- ( )=0
at J an J at J
a 2z a lz a
zz z
-( )+- ( )+- )=
at J an J ac J

14 -

By applying the generalized transformation to the

compressible Navier-Stokes equations (9) in vector form,

dividing by the Jacobian and then rearranging into

conservation law form with the use of metric identities,

equation (15), the following form of the transformed

Navier-Stokes equations is obtained:

aq 8E aF aG aEv F aG
+ + + = Re ( + -- + -- ) (16)
at as an as as an as

in which

q = J-q = [ p pu pv pw E ]

a8 a8 a8 as
E = J- ( -q + -E + -F + -G )
at- ax ay az
SJ-1 [ pU pUu+Sxp pUv+Syp pUw+Szp (E+p)U-Stp ]T

an an a8 n aa n
F = J- ( --q + -E + -F + --G )
at ax ay az
SJ-1 [ pV pVu+nxp pVv+n p pVw+nzp (E+p)V-ntp T

as as ac as
G = J- ( -q + -E + -F + --G )
at ax ay az
= j-1 [ pW pWu+xp pWv+yp pWW+Cp (E+p)W-etp ]T


1 x xx y yx + z zx
EV j tx xy + tyyy + tz zy

SXX + SyTyz + SzT
Zx xz y yz + z zz
xOx +tyy z z

- 15 -


S x xx + yTyx + zTzx
Fv Ax xy Ay yyy Az zy

x xz y yz A+ zzz

x 0 + y y + zBz


1 x xx +y yx +z Tzx
v xTxy + Cyy z Tzy

x Txz + yIyz + z zz

Cxx + yy y + -zoz

The stresses r and temperature gradients in the B terms

contain the metrics of the transformation; e.g.,

au av au 3u au av av av
xy=1r(-----)=V(-y +y y y-- x-- x ---cx )= xy
ay ax at aA ac a& an aS

etc. U, V and W are the contravariant velocity components

along the three curvilinear coordinates A, 1 and i,

respectively. They are defined in terms of the Cartesian

velocity components u, v and w as

U = tt + U +x + Vy + wz

V = At + unx + Vy + wnz (17)

W = 4t + ucx + Vy + wz

16 -

Thin Layer Navier-Stokes Equations

In general, a very large amount of computer time and

extensive storage are needed for solving a complex

,aerodynamic problem governed by the complete Navier-Stokes

equations (16). At high Reynolds numbers, however, the flow

around the body can be considered to be mostly inviscid,

except in a thin-layer region close to the body surface

where the viscous effects become important. With this in

mind, the thin-layer approximation is used to simplify the

Navier-Stokes equations in generalized coordinates.

The thin-layer approximation to the Navier-Stokes

Equations is based on an assumption that the diffusion

processes along the body surface are negligibly small in

comparison with the diffusion process in the direction

normal to the body surface. This is generally true for high

Reynolds number flow. 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. Assume that a body surface is mapped into a

4I-plane in the transformed space. Then and T derivatives

of the viscous terms can be neglected. Hence, the thin-layer

Navier-Stokes equations in curvilinear coordinates are

17 -

3q 3E aF 8G 3S
+ + + = Re (18)
at ac aa c



mlu + m2Cx

S = -1 mlv + m2 y

ml1w + m2Cz

mlm3 + m2( xu + yV + z w)

S= (x2 + y2 + z2)

m2 = P(xuC + yV + 4zw )/3

m3 = (u2 + v2 + w2) /2 + Pr-1(-l)-l(c2 )

and c= J/RT is the local speed of sound. Since the

thin-layer Navier-Stokes equations (18) retain the normal

momentum equation, the pressure can vary through the viscous

layer. Hence, the thin-layer model is more general than the

boundary layer equations and is capable of predicting flow


Axisymmetric Thin-Layer Navier-Stokes Equations

For an axisymmetric flow in which the flow variables are

invariant in the n-direction, the three-dimensional

thin-layer Navier-Stokes equations can be further simplified

by imposing the following restrictions: (i) the body

18 -

geometry is axisymmetric; (ii) the flow variables and the

contravariant velocities do not vary in the n-direction.

Hence, the 8F/ai term of equation (18) can be reduced to a

source term H. Based on the relationships between the

coordinate systems, the Cartesian coordinates (x,y,z), the

cylindrical coordinates (x,O,R) and the curvilinear

coordinates (t,n,,) shown in figure 1, one has

x = x(.S,C,' )

y = R( t,C,T ) sing

z = R( t,(,T ) cos#

By evaluating the metric terms with the use of the above

relations and substituting them into transformed thin-layer

Navier-Stokes equations (18), the resulting unsteady

axisymmetric thin-layer Navier-Stokes equations become

8q 3E 8G as
+ + + H = Re (19)
at at at 8

With the metrics and flux vectors evaluated at plane 0 = 00,

the source vector H [9,11] can then be expressed as

- 19 -

Projectile-type body


x = constant plane


# = constant plane

Figure 1. Axisymmetric body and coordinate system

- 20 -

H = J-1


pV[R (U-Vt)

-pVR (V-nt)



+ R (W-Ct)]

- p/R


which has replaced aF/8n of equation (18).

Equation (19) contains only two spatial derivatives but

does retain all three momentum equations by allowing

non-zero velocity components in the invariant n-direction.

Thus, flow over axisymmetric spinning projectiles can be

solved with equation (19).


There have been various numerical algorithms developed

for approximating the Navier-Stokes equations. For the

thin-layer Navier-Stokes codes obtained from the U.S. Army

Ballistic Research Laboratory and used in this study,

equations (18) and (19) are approximated by the Beam and

Warming scheme [3,4,21], which is an implicit approximate

factorization finite difference scheme. This effective

finite-difference scheme is described and discussed in the

following. Note that for a fixed grid system, the T equals

to t.


Consider three different temporal differencing

approximations given by

(1) Backward Euler implicit formula

n+l n n+1 2
q = At (-)l + O[(At)]

- 21 -

- 22 -

(2) Trapezoidal formula

At aq +q
q+l qn [(-) + ()n+l] + O[(At)3]
2 at at

(3) Three-point backward formula

3qn+l 4qn + qn- = 2At (--)n+ + O[(At)3]

By defining Aqn = qn+ qn, the-above formulas can be

combined into a single generalized time-differencing formula

in delta form

qn+l 8qn
rAt = (1 + s)Aqn sAqn- At (1 r)
at at

+ 0[(r-s-1/2)(At)2 + (At)3] (20)

With the appropriate choice of the parameters r and s,

equation (20) reproduces many familiar two- and three-level

explicit and implicit schemes which are listed in the

following table.

23 -

Table 1. Partial list of schemes contained in equation (20).

r s Scheme Truncation error

0 0 Euler, explicit O[At]

0 -1/2 Leapfrog, explicit O[(At)2]

1 0 Euler, implicit O[At]

1/2 0 Trapezoidal, implicit O[(At)2]

1 1/2 3-point-backward, implicit 0[(At)2]

After substituting equation (18) into equation (20), we have

rAt 3E aF
Aqn + --(- + --
l+s at a3i

l-r aE
-At --(-
1+s 8a


+ +

as s
-1 n+1 n-1
Re --)n Aq
;c 1+s

aG as
- Re-1 n
a4 a

which is either first order or second order in time,

depending on the choice of r and s as indicated in the above


The flux vectors E, F and G at the advanced n+1 time step

are nonlinear functions of the unknown vector q. It is

possible to remove the nonlinear nature and still maintain

the order of accuracy of the finite difference approximation

by using Taylor series expansion about qn


- 24 -

En+1 = En + AAqn + O[(At)2] where A = -

F+1 = Fn + BnAqn + O[(At)2] where B -

Gn+1 = n + CnAq + O[(At)2] where C -

Matrices A, B and C are so called Jacobian matrices. The

viscous flux term can also be linearized by using a method

suggested by Steger [16]. In order to apply this

linearization method, the coefficients of viscosity p and

thermal conductivity k are assumed to be locally independent

of q. As a result of this assumption, elements of S have the

general form

s = J-1

where ai is independent of q = J- q given in equation (16),

and Bi is a function of q. These elements are linearized in

the following manner by Taylor series expansion

sin+1 = sin + -1 an- [ ( n Aqmn] + O[(At)2
a; m aq

25 -

These relations can be written in a vector form as

Sn+1 = Sn j-lnJAqn

which has a truncation error of second order in time. The

explicit expression of Jacobian matrices A, B, C and M is

given in the Appendix. Equation (21) now becomes

rAt aA aB aC a
[I + --(- + + Re--(J-IMJ)]nAqn
l+s at an ac ac

s At 3E aF aG as
Aq1 --(- + -- + Re-) (22)
l+s l+s at an ac a

where the matrix I is the identity matrix.

Approximate Factorization

To increase computational efficiency, approximate

factorization is employed for the solution algorithm of

equation (22), which reduces the solution process to three

one-dimensional problems at a given time level. The solution

process for Aqn then involves the inversion of three

matrices with smaller bandwidth, for which there exist

efficient solution algorithms. Consider the following

factorized equation approximating equation (22).

26 -

rAt 8A rAt 3B
[ I + [ n + -- -
l+s a l+s an

rAt aC a
[ I + Re- -(J_-MJ) ] Aqn = RHS(22) (23)
l+s ay at

Here RHS(22) is used to indicate the right-hand side of

equation (22). The error introduced by the factorization

procedure, equation (23), is the leading third order term

At3 aA aB Ba aC aC aA At4 3A aB sC 8qn
[--(- -+- -)+-- -]n -+O[(At)5]
4 a~ an an ac ac a 8 a3 a3 a8 at

which does not affect the second order accuracy of the

difference approximation.

Numerical Dissipation

For a complex aerodynamic problem, a direct application

of equation (23) often experiences convergence problems.

Hence, dissipation terms are added to the equations to damp

out high frequency oscillations and to control numerical

instabilities. Several different numerical dissipation terms

have been investigated, e.g. reference 13. Currently, the

numerical dissipations proposed by Beam and Warming [3,4]

and Steger [16] are employed in the thin-layer Navier-Stokes

codes. The second order implicit dissipation operators

- 27 -

D = -eI j-1V A J

Dn J-1V A (24)

D = -E J-1V A

are, respectively, added to the left hand side implicit

one-dimensional 5, n and C operators of equation (23) while

an explicit fourth order dissipation operator

DE = -E J-1[(V A)2 + (VA )2 + (V A)2]J (25)

is added to right hand side of equation (23). The symbols A

and V are the standard forward and backward difference

operators, that is,

Afi = (fi+ fi) + O(h)

Vfi = (fi fi-) + O(h)

Finite Difference Equations

The spatial derivatives appearing in equation (23) have

to be represented by finite difference quotients. In the

computer codes, central difference approximations are used.

Thus, the equation (23) is rewritten as

- 28 -

LE LL L Aqn = Rn (26)


L + 6 An + D

L I + 6 Bn + DI

rAt rAt
L I + -- 6 Cn Re-16 (J-1nJ) + D
1+s l+s

s At
Rn Aq- (6E + 6 F + 6 G Re-1S) + DE
l+s l+s

At this point, equation (26) can be first or second order

accurate in time depending on the choice of parameters r and

s. For the finite difference approximation of spatial

derivative terms, a second order central differencing

6f= (fi fi-) + O(h2)

is employed on the implicit part of 6 A, 6 B and 6 C for

maintaining the block tridiagonal matrix structure, while a

fourth order central difference approximation

6f (-f.+2 + 8fi+1 8fi- + f ) + O(h4

29 -

is used on the convection terms 6 E, 6 F and 6 G of the

RHS(26) to keep the convection truncation error from

exceeding the magnitude of the truncation error of the

viscous terms 6 S which is approximated by second order

central difference. In this study, the Euler implicit

two-level scheme is used; i.e., r=l and s=O. Therefore, the

solution algorithm is first order accurate in time and

second order accurate in space.

The unknown vector q at n+l time level governed by

equation (26) can now be computed effectively as follows:

** Rn
Lq Aq = R
LE Aq = Aq

L Aqn = Aq*

qn+l = q + Aqn

where Aq and Aq are intermediate solutions. The

numerical process thus results in a block tridiagonal matrix

inversion for each coordinate. The size of the blocks

depends on the number of elements in the unknown vector q.

Finite Difference Scheme for Axisymmetric
Thin-Layer Navier-Stokes Equations

In a manner similar to the three-dimensional case, the

application of Beam and Warming noniterative implicit scheme

30 -

to the axisymmetric thin-layer Navier-Stokes equations (19)

results in the following finite difference equations:

L Lt Aqn = Qn (27)


s At
Qn Aqn-1 ---(6 E + 6 F+ 6 G Re-16 S)n
1+s 1+s

Hn E J-[(V A)2 + (VA ) n

and Lt and L, are given in equation (26).

Again, the Euler implicit scheme, i.e., r=l and s=O, is

used in the axisymmetric thin-layer Navier-Stokes code which

results in first order accuracy in time and second order

accuracy in space.


An axisymmetric thin-layer Navier-Stokes code and a

three-dimensional thin-layer Navier-Stokes code obtained

from the U.S. Army Ballistic Research Laboratory (BRL) are

employed in this study for the numerical simulation of

transonic turbulent flows past a projectile model with a

sting. The axisymmetric code is a product of BRL and NASA

Ames joint effort [11]; it is a simplified version of the

three-dimensional code developed at NASA Ames [14]. The

application of these codes to transonic projectile

aerodynamic problems has been investigated to some extent by

the U.S. Army Ballistic Research Laboratory.

The three-dimensional code and the axisymmetric code are

based on the transformed thin-layer Navier-Stokes equations

(18) and (19), respectively, which are valid only for high

Reynolds number ideal gas flows. In these codes, the

characteristic velocity used in the transformation, equation

(8), is the speed of sound at freestream condition c

instead of the freestream velocity V ; consequently, the

freestream Mach number M in equations (10) and (12) becomes

unity. The solution algorithm adopted for the transformed

governing,equations (18) and (19) is a noniterative,

implicit, approximate factorization, finite difference

31 -

32 -

scheme of Beam and Warming which has been discussed in the

previous chapter. It is noted that both codes have an option

for solving inviscid flow problems while a steady state

solution is obtained as the converged solution of the

unsteady impulsive flow problem.

Both the axisymmetric code and the three-dimensional code

have been studied in detail. Necessary modifications to the

computation of turbulent eddy viscosity have been made to

the axisymmetric code and are discussed in the following;

furthermore, the code has been vectorized for Cray X-MP

supercomputers to improve the efficiency of computation.

Moreover, essential subroutines for the computation of

aerodynamic forces have been written and implemented into

both axisymmetric and three-dimensional codes.

For a successful application and/or modification of a

Navier-Stokes code, one must be familiar with certain

relevant features of the code. The essential features of the

Navier-Stokes codes considered and the subroutines

implemented are described and discussed in the following

subsections for the transonic projectile aerodynamics

problems investigated.

Flow Domain

In this study, flows past an axisymmetric projectile

model with a sting are considered for the investigation. For

axisymmetric flow problems, a schematic planar flow domain

33 -

and the corresponding computational domain, ABCD, are shown

in figure 2. Here, the line AB is the solid boundary of the

projectile model, the line BC is the exit or downstream

boundary and the line CD is the outer boundary while the

line AD is called the upstream center line. For a projectile

model at an angle of attack, that is, a three-dimensional

flow problem, the flow domain in the physical space is

generated by the revolution of the planar domain about the

axis of the projectile model. Hence, a grid network for a

three-dimensional problem can be formed by a sequence of

positions of the planar grid ABCD around the axis of the

projectile model. The corresponding grid system in the

computational domain is then formed by stacking the

rectangular grid domain ABCD. Accordingly, as shown in

figure 2, the plane ABB'A' is the solid boundary of the

projectile model, the plane BCC'B' is the downstream

boundary, the plane CDD'C' is the outer boundary, while the

plane DAA'D' is a singular plane representing the upstream

center line AD.

Boundary Conditions

For a selected boundary-fitted transformation, such as

the one discussed above, suitable boundary conditions must

be specified. The boundary conditions in the thin-layer

Navier-Stokes codes are applied explicitly and require the

- 34 -



*l 0







35 -

specification of all variables p, u, v, w and E in unknown

vector q on the line ABCD.

The application of solid boundary conditions, line AB, is

simplified through the boundary fitted curvilinear

coordinate transformation since the body surface has been

mapped onto the plane of 4=0. The no-slip boundary condition

for viscous flow is enforced by setting the contravariant

velocities to zero, i.e., U=V=W=0 on the projectile surface.

For inviscid flow problems, a linear extrapolation of

tangency contravariant velocities U and V is used while W=0

for impermeable wall. The velocity components u, v and w

are then obtained by inverting equations (17). Densities at

the surface are assumed to be equal to the values calculated

on the first grid line adjacent to the boundary, which was

always well within the viscous sublayer. The pressure along

the body surface is obtained from the three transformed

inviscid momentum equations. By combining [ xx(&-momentum) +

Syx(n-momentum) + z x(C-momentum)] with the use of the

continuity equation, metric identities (15) and neglecting

viscous terms, one finds

(txCx + y y + z z)P + (1x4x + ,yy + +zz)p,
+ (2 + y2 + z) (28)

= P(at + uax + va y + wa )

pU( xU + cyvt + zw) pV(Cxu + yv + CZW)

36 -

Since viscous terms have negligible effect on surface

pressure, the same relation (28) has been used in viscous

flow computation [14,16]. Equation (28) is solved at the

surface by using second order central differences in the (

and T directions and second order forward differences in the

C direction. The total energy E is then obtained from the

known pressure, velocities and density at the surface by

using equation (10).

Along the upstream line of symmetry, line AD, which

emanates from the projectile nose, symmetry conditions are

imposed by setting the second order finite difference

expression for the first derivative equal to zero. In

three-dimensional problems, the flow variables are

determined along the line AD independently from each plane

containing line AD. The averaged values obtained from all

planes are then used on the singular plane DAA'D'. Along

the far field outer boundary, line DC, freestream conditions

are imposed. On the downstream boundary, line CB, linear

extrapolation is used for all flow variables as the exit

plane flow conditions for M 21. For M <1, all flow

variables are obtained from linear extrapolation except the

total energy E which is obtained from equation (10) by

assuming the freestream pressure p. fixed at the exit plane.

37 -

Initial Condition

A steady flow problem solved by the Navier-Stokes code is

treated as an impulsive flow problem. The computational

procedure is then started marching in time until an

asymptotic steady state is reached. Hence, for the

projectile flow problem, the freestream conditions are taken

for the initial conditions. To ease the sudden jump in

solid boundary conditions when flow impacts the projectile

impulsively, the boundary conditions at the surface are

turned on slowly over the first 30 time steps. Thus, the

boundary conditions of the unknown vector at the body

surface are scaled by

q = qbc x SL + (1 SL) q4


SL = (10 15 x r + 6 x r2) r3

r = (NC 1)/30

where NC is the number of iterations, qbc is the correct

boundary conditions of the unknown vector which has been

discussed in the previous section.

38 -


The effective coefficient of viscosity i in the

Navier-Stokes equations is comprised of two parts, the

laminar viscosity py and the turbulent eddy viscosity pt.

In the thin-layer Navier-Stokes codes, the evaluation of

laminar viscosity pt as a function of temperature is

obtained from Sutherland's theory of viscosity and given by

V T T + S 1 +S
S = ( )15 = T*5 ( )
O T T + S T +S

where the constant S is 198.60R [15] and S =S/T .

The turbulent eddy viscosity model employed in the codes

is the one proposed by Baldwin and Lomax [2]. It is a

two-layer algebraic model in which an eddy viscosity is

calculated for an inner and an outer region. The inner

region follows the Prandtl-Van Driest formulation

(Pt)inner = pt21wl (29)

with e = ky[l exp(-y+/A+)]

where 7 is the normal distance to the surface and Iwl is the

magnitude of the vorticity

|w = [(a y-axv)2 + (azv-ayw)2 + ( zu)21/2

39 -

and y is the law of the wall coordinate

S(Pw w) 1/2
y =

where w, P and v are the local shear stress, density and

laminar viscosity evaluated at the wall, respectively.

The eddy viscosity for the outer region is given by the

Clauser formulation

(ut)outer = pKCcpFwakeFkleb() (30)

where K is the Clauser constant, Ccp is an another constant

and Fwake is the smaller value of the following:

Fwake = YmaxFmax
Fwake = CwkYmaxdif /Fmax

The quantities ymax and Fmax are determined from the


F(Y) = ylwl[l exp(-y+/A+)]

where Fma, is the maximum value of F(y) and yax is the

value of y at which it occurs. The Klebanoff intermittency

factor Fkleb(Y) is given by

40 -

Fkleb() = 11 + 5.5(CklebJ7/yax)6 -1

and the quantity Udif is the difference between maximum and

minimum total speed

dif = (u2 + V2 + W2)max1/2 2 + + 2)min/2

where the minimum total speed is zero, unless the spinning

projectile is considered.

Rather than attempting to match (Vt)inner and (Vt)outer

values in the overlap region, both (lt)inner and (Vt)outer

are calculated at every grid point in the flow field and the

smaller value is used as turbulent eddy viscosity. The

constants used for this model are A+ = 26, C = 1.6, Cwk
cp Wj
0.25 and Ckleb = 0.3 Note that in this section, all

quantities are dimensional except those with a superscript

In the original axisymmetric thin-layer Navier-Stokes

code, the turbulent eddy viscosity was calculated only for

the first 25 grid points off the projectile surface in order

to increase computational efficiency. However, in the

current grid networks used, the 25th point is well within

the boundary layer and subsequent zero values for the

turbulent eddy viscosity in the remaining region of the

boundary layer led to instabilities in the numerical

algorithm. The code has been modified to include the

41 -

'turbulent eddy viscosity at every grid point in the entire


Forces and Moments

The motion of the fluid relative to a solid boundary

exerts a dynamic force on the surface which consists of both

frictional and pressure forces. These surface forces can

create the pitching, yawing and rolling moments relative to

the center of gravity. In the cases we considered, the

axisymmetric projectile without spinning, there is no

contribution to the side (lateral) force, nor the yawing and

rolling moments. Note that both pressure and shear stresses

can cause the pitching moment for the projectile at an angle

of attack.

The forces and the pitching moment acting on an arbitrary

shape can be obtained by integration of the surface pressure

and shear stress distributions and moments about the center

of gravity, respectively, as follows:

(1) Force components in axial direction

FpA = I p sin( dA
A (31)
FfA = I$ T cos9 dA

42 -

(2) Force components in normal direction

FpN = If p cos8 coso dA
FfN = /f T sin8 cos0 dA

(3) Pitching moment about the center of gravity

MC.G. = II [FN (XC.G. X) + FA Z cos#] dA (33)

Here, p and T are the pressure and shear stress

distributions on the body. The subscripts p and f refer to

the forces due to pressure and shear stresses, respectively,

and subscripts A and N indicate the force components in the

axial and normal direction, respectively. The angle *

starts on the leeward side and measures in the

circumferential direction and 0 is the inclination of local

slope of the surface, as shown in figure 3. XC.G. is the

axial position of the center of gravity measured from the

nose of projectile which is 3.6 for this particular

projectile. X is the local axial position and Z is local

position in the radial direction, measures from the axis of

the projectile. The pitching moment is considered positive

when it is in the nose-up sense.

With the same dimensionless quantities in equations (8)

and characteristic velocity c forces in the equations (31)

and (32) can be put into

- 43 -




0) r

C- W

-U U
^! c

,4 )

a. ,




44 -

Fp = L2 p co F* = F *

1 *
Ff = L yi c Ff= Ff

The ratio of the forces due to pressure and shear stress is
therefore F to Ff /Re.

The pressure and shear stress coefficients are expressed


p p. P p.
p ( l / 2 ) p mV 2 ( 1 / 2 ) M .(

T /Re
f (1/2)p V (1/2)M2

and the pitching moment coefficient as

Cm (36)
(1/2)p V .L (1/2)M. L

Because the moment has dimensions of a force-length product,

the additional length scale L in the denominator of moment

coefficient is introduced.

The component of the resultant force in the direction of

the freestream velocity V that passes the body is the drag,

which is one of the most important parameters in aerodynamic

designs. It can be expressed as

- 45 -

FD = FN sina + FA cosa (37)

where a is the angle of attack. The corresponding lift force


FL = FN cosa FA sina (38)

For axisymmetric flow cases, i.e. a=00, the drag force will

be the axial force and lift force will be zero.

The computation of forces, moment and shear stress

distribution has been added to both thin-layer Navier-Stokes

codes. Thus, the computer codes have become more complete in

evaluating the flow features.

Rate of Convergence

The numerical algorithm is a time-marching scheme in

which a steady state solution is obtained as a time

asymptotic solution of the unsteady equations. The

convergence of the process can be observed from the average


(R 2 1/2
Average residual = fm (39)
At /(JMAX-2)(KMAX-2)(LMAX-2)

46 -

where R is either the vector of RHS(26) or RHS(27). The

total numbers of grid points in the n and C direction are

JMAX, KMAX and LMAX, respectively. And summations over j, k,

i range from 2 to JMAX-1, KMAX-1, LMAX-1, respectively. The

summation of m is from 1 to the total number elements of

unknown vector q, which for this case is 5.

For steady state solutions, from equation (26) or (27), we

know that the vector R is zero since the solution will no

longer vary. This implies then that the average residual

becomes zero when the steady state solution is reached. In

the numerical process of this study, the convergence

criterion for the steady state solutions being reached has

been set to the average residual less than 10-4.


In the applications of thin-layer Navier-Stokes codes, a

grid network must be provided. The grid network used in this

study is generated from a grid generation code, named

GRIDGEN, obtained from the U.S. Army Ballistic Research

Laboratory. This code can provide a grid network generated

either by elliptic or by hyperbolic equations. For external

flow problems, hyperbolic grid generation is feasible.

Since specifying the precise positions and shape of the

outer boundary is not important, the grid around the

geometry of interest can be obtained by the specification of

grid positions at one boundary, and then marched outward

from this boundary until the far field is reached. The

hyperbolic grid generation is a noniterative scheme, so the

required computer time is almost equivalent to that of one

iterative sweep in solving elliptic grid generation

equations. Furthermore, this hyperbolic grid has the

orthogonal property, and grid smoothness can be controlled

by a selected grid distribution function.

- 47 -

48 -

Governing Equations

Construction of hyperbolic planar grid (x,z) proceeds by

selecting two differential equations which serve as

constraint equations that are solved to obtain the

coordinates of grid points [17]. It is known [8] that an

orthogonal grid will reduce the truncation error terms

associated with finite difference approximation formulated

in generalized coordinates; thus one equation selected is

that which enforces orthogonality

V *0 V4 = 0 or &x

For a two-dimensional general coordinate transformation,

the metric relations and Jacobian J of the transformation


z x z& xt
S x z (41)

J = xz4 xcz (42)

The equation (40) can then be written in the computational

space (&,4) as

xx + zz4 = 0


49 -

The second equation is chosen to specify a finite grid

cell area to insure a nonsingular mapping through the

transformation. In numerical implementations, AE=AC=l, the

Jacobian in equation (42) is approximately equal to the

physical cell area. A set of hyperbolic grid generation

equations is then formulated by equations (43) and (42).

Equations (43) and (42) are nonlinear equations which can

be locally linearized by expanding about a known state

(xo,z0). Neglecting the high order term, e.g.,

(x-xo")(x-x*),, we have

x0x = [x + (x x [ + ( )]

= x 4x + xox x x

etc. where the superscript 0 denote the known state. Upon

substituting these relations into equations (43) and (42),

one finds

x ox + z oz + x ox + z z = 0
z 0x x oz z0x + x 0z = J + J0

where the grid areas J are to be specified and JO is

determined at the known state. Therefore, all quantities of

RHS(44) are known. Equations (44) can-be put into the

compact form

- 50 -

AW + BW = f (45)


x 0 z C B x z f = 0 x X
A = B = f= W =
z o -x 0z o J + 0 z

Assume that the C coordinate runs in the direction away

from a specified boundary. Then the governing equations for

the propagation problem can be written as

CW + W = g

where C = B-1A g = B-1f

The eigenvalues of C are two distinct real roots which

implies that the governing equations (45) form a hyperbolic


Finite Difference Form

In the GRIDGEN code, the finite difference approximations

used for the hyperbolic grid equations (45) are second order

central difference for t derivatives and first order

backward difference approximation for derivatives in the

marching direction C; i.e.,

- 51 -

( ),t ( W jlt+ Wj-l,t+ ) + o( A2 )

( WW )1 = ) + O( A )

Then, an implicit finite difference scheme can be formed

1 1
----A e W'+1 + + Bj wj,' --A2 'W]-1l,+1

= ft+ + B, eWi,

which forms a 2x2 block-tridiagonal set of equations to be

solved on each successive constant C-line. Note that the

column matrix f ,e, on RHS(46) is assumed known, since it

only involves the quantities J+JO. The grid areas J on the

next e+ constant C-line can be specified in any way. In

the GRIDGEN code, it has been specified as

Jj,+= AL x ASt

where AL is the total length of two adjacent cells, as shown

in figure 4, and ASe is the specified height of two

successive constant 4-lines. The distribution AS used in the

code is an exponential distribution, and is determined by

the relation

- 52 -

I -.f

Figure 4. Specification of grid area in GRIDGEN code.

53 -

ASZ = ASO (1 + )i i = 1,IMAX-1 (47)

where IMAX is the total number of grid points along the

given arc length, and ASO is the specified grid spacing at

one end of arc. The parameter e is determined by a

Newton-Raphson iteration process so that the sum of the

above increments matches the known arc length. The grid

spacing ASO next to the boundary considered in this study is

0.01 for inviscid flow problems; however, for viscous

turbulence cases, first grid spacing ASO of 2x10-5 is


Several different distribution functions have been

discussed and investigated in detail [19,20,22], and, in

comparison, a hyperbolic tangent distribution has the better

overall distribution [19,20]. For extensive study of the

effectiveness of the grid network on the flow solution, a

hyperbolic tangent distribution is implemented in the

computer code. A hyperbolic tangent distribution is

constructed by

6 1
tanh[-( 1)
2 IMAX 1
Si = Stot {1 + tanh(

ASi = Si. Si = 1, IMAX-1

54 -

where Stot is the total length of the arc and the parameter

6 is determined by

6 AS0
S = (IMAX 1)
sinh(6) Stot

Two computational planar grids for the projectile

configuration with the same boundary distribution and

ASO=2x10-5 grid but different normal distributions,

equations (47) and (48), for the projectile configuration

are shown in figure 5. We observed that the grid with

hyperbolic tangent clustering provides a smoother grid away

from the projectile than the one with exponential clustering

in the radial C-direction. However, near the body surface,

the grid with exponential clustering is smoother than

hyperbolic clustering.

A hyperbolic grid generation is a propagation problem; it

requires the initial grid points. The grid point

distribution on the body surface, in this study, is adapted

to the flow solution in the t-direction and then propagated

outward in (-direction. Furthermore, a three dimensional

grid for the axisymmetric projectile geometry considered is

generated by rotating a two dimensional planar grid about

the axis of symmetry and maintaining a uniform angular space

between the rotated planes.

55 -




-16e e 18 /

(a) Based on exponential clustering.




-18 18i X-'

(b) Based on hyperbolic tangent clustering.

Figure 5. 90x40 hyperbolic grids with different clusterings
in the normal direction.


The objectives of this study are to gain further

understanding of complex flow simulations using the

thin-layer Navier-Stokes approximation and to investigate

the relative importance of the viscous force to the pressure

force on a projectile at transonic speed.

The geometric configuration considered is an axisymmetric

secant-ogive-cylinder-boattail (SOCBT) projectile. The

projectile model has a 3-caliber secant-ogive part followed

by a 2-caliber cylinder and 1-caliber 7-degree boattail as

shown in figure 6. The model was supported by a base-mounted

sting in the experiments. Experimental data are available

for the pressure distribution along the projectile surface

for a range of Mach numbers in the transonic regime, 0.91,

0.94, 0.96, 0.98, 1.10 and 1.20; at angles of attack of 00,

20, 40, 60 and 100; and at circumferential positions around

the model in 22.50 increments [7]. The computational model

of the projectile, however, is formed by extending the

boattail part of the projectile for another 1.77 calibers to

meet the sting. This modification is made to avoid the

difficulty of simulating the base flow region.

For the external flow problems, which are considered here,

the outer boundary needs to be specified. Theoretically, the

56 -

- 57 -




58 -

outer boundary should be as far as possible from the

projectile, so that freestream conditions can be imposed.

The wind tunnel experiments supply the freestream stagnation

temperature TO=580R; moreover, the projectile body is

assumed adiabatic and impermeable.


For the purposes of evaluating the accuracy and

applicability of the thin-layer Navier-Stokes codes

described previously for complex flow problems, comparisons

of measured and predicted aerodynamic characteristics and

associated flow fields for various flow conditions have been

made. The finite difference scheme for the Navier-Stokes

equations requires the additional artificial second order

implicit and fourth order explicit dissipation terms,

equations (24) and (25) respectively, to control numerical

instabilities. According to a linear stability analysis, the

scheme is unconditionally stable for linear equations when

EI=2eE [3,4,21]. This weight is used in all subsequent

investigation of axisymmetric flow problems.

Inviscid Flow

Flow Domain Determination

For SOCBT projectile model, the experimental data are

available in a range of M =0.91 to 1.20. A flow case of

M =0.96 is selected for detailed investigation. In order to

59 -

- 60 -

determine adequate flow domains, the axisymmetric

Navier-Stokes equations (19) were solved on a 90x40

hyperbolic grid with exponential clustering, equation (47),

in the normal direction. The two domains used, Stot=18 and

24 are about 3 and 4 times the projectile-length,

respectively. The results presented in figure 7 show that

the computed surface pressures are very similar and give an

accurate prediction of the experimental data except on the

boattail region. The results indicate that both domains 18

and 24 seem to be sufficiently large for imposing freestream

conditions on the outer boundaries. Nevertheless, in theory,

the better solution is obtained with a larger domain when

enough grid resolution is provided.

Effect of Artificial Dissipations

Since the flow domain of 18 can provide accurate results,

the 78x28 hyperbolic grid with domain 18 was used with two

different sets of dissipation terms: one is EI=2E=4At;

another is EI=2EE=8At for the inviscid flow case. The

influence of artificial dissipation terms on the flow

solutions is presented in figure 8. When more dissipation is

added, smearing of the numerical solution occurs near the

shock/boundary layer interaction' region, about at the middle

of cylinder part, since the dissipative nature tends to

smear the predicted shock over several mesh intervals. In

order to prevent the deterioration of the solution, the

- 61

M =.96

2. 4. 6. XD
Figure 7. Computed surface pressure of inviscid solution
with a 90x40 grid but different domain.
: Sot=18 ........ ot24
tot'l tot2

M = .96

2. 4. 6. WD
Figure 8. Computed surface pressure of inviscid solution
with a 78x28 grid but different dissipations.
-- : E-I=2EE=4At ------: EI=2EE=8At

- 62 -

added dissipation terms should be as small as possible;

however, the solution algorithm became unstable when

CI=2EE=2At is employed.

Viscous Flow

In turbulent flow applications, the closure turbulence

model used in the codes is a two-layer algebraic eddy

viscosity turbulence model of Baldwin and Lomax. The

computed results obtained show that a 90x60 hyperbolic grid

with exponential clustering in C-direction and flow domain

Stot=24, which is about 4 times the projectile length, can

give very accurate solutions. Figure 9a shows that the

computed surface pressure coefficient Cp is in excellent

agreement with measured data while the computed shear stress

T, presented in figure 9b indicates that there is no flow

separation on the projectile surface. The corresponding

pressure and Mach contours are shown in figure 10. The dense

regions in the contour plots indicate that large gradients

occur; i.e., it implies the corresponding flow property

changes rapidly. The coalescence of the Mach lines to the

right of the supersonic region represents the position of

the shock. Two shock locations, one in the middle of

cylinder part and another around the middle of boattail

portion, agree quite well with a shadowgraph in reference 7.

The expansion waves also are shown at the ogive-cylinder and

- 63 -

M = .96

2 4 6 X/D

(a). Pressure coefficient;


0: measured data.


2 4

6 X/D

(b). Surface shear stress TC

Figure 9. Computed results with a 90x60 grid.






- 64 -

0 2 4 6
(a) Pressure contours
p: 0.48+0.86, Ap=0.02

u 2 4 6
(b) Mach contours
M: 0.88-1.0, AM=O.01; M: 1.0-1.26, AM=0.02

Figure 10. Contour plots of M =0.96.

65 -

cylinder-boattail junctures. The computed nondimensional

drag force, D=FDx10- due to surface pressure is Dp=40.28

over the 6-caliber projectile, excluding the base drag,

while the drag force resulting from skin-friction is

Df=20.53. This shows the importance of viscous flow

computation for accurate aerodynamic force prediction.

Effect of Artificial Dissipations

The effect of artificial dissipation terms upon the

accuracy of numerical solutions of turbulent flow cases was

next considered. The same grid network 90x60 of domain

Stot=24 was used but with EI=2EE=8At. The computed C as

shown in figure lla seems to agree rather well with the

accurate results, while the quality of the solution is

degraded near the shock/boundary layer interaction region

owing to the nature of dissipation. The resulting pressure

drag is Dp=32.91 which is 18% less than that of the accurate

solution. Figure llb shows the difference between the

calculated shear stress distributions; however, the

resulting shear drag is about the same, Df=20.12. The

difference in the drag forces due to pressure is mainly due

to the difference in the predicted pressure on the boattail

part, since the pressure intensity on the cylinder portion

does not contribute to the drag force in axisymmetric flow


p- 66 -

M = .96

Station 23 39 46

2. 4. 6. oD
(a). Pressure coefficient; 0: measured data.

M = .96


iL 2

S2. 4. 6. XoD
(b). Surface shear stress C
Figure 11. Computed results with a 90x60 grid but
different dissipations.
---- : I=2E=4At ------: e=2E =8At

- 67 -

Effect of Eddy Viscosity Distribution

The effect of an averaged-technique originally

implemented in the axisymmetric thin-layer Navier-Stokes

code for computing eddy viscosity also has been investigated

on 90x60 grid at M,=0.96. Figure 12 shows the distribution

of eddy viscosity with and without the averaging technique

at three different boundary point stations identified in

figure lla. It is clear that the averaging scheme yields

improper zig-zag eddy viscosity distribution above the

viscous sublayer; nevertheless, it has negligible effect on

surface pressure. In fact the computed C -distribution is

almost exactly the same as that of the accurate solution

shown in figure 13; however, the shear stress computed is

consistently less than that of the accurate solution over

the entire projectile surface. The corresponding pressure

drag force is Dp=41.41 and skin-friction drag is Df=18.99,

which are different by 3% and 8% respectively.

Effect of Grid Resolution in Streamwise Direction

The grid networks 78x40 and 90x40 with domain Stot=18 are

used for comparison. Of the 12 additional points in the

t-direction of the 90x40 grid network, 2 points were added

in the ogive part and 10 points were added in cylinder part.

From the results presented in figure 14, we observed that

the solutions are generally quite the same except in the

- 68 -

.. Stat ion
M) ... ------ 23

o .. /

c ^ "* ~ '- ^
L ..5 -..-. .....

......... ...... -.--- ---

.. ... ... ....
S*** .' ...

-^-'-'-'--^^ *^'^-'----^*

5 8 1 0 1 5 0
Eddy Viscosity
(a) Sample eddy viscosity distributions based on
an averaged-scheme implemented in the original
axisymmetric thin-layer Navier-Stokes code.


S- 23

S\ ............ 39
S\ ---- 46

. -
5o 100 158

) ,**. /

50 1800 158
Eddy Uiscosity
Figure 12. (b) Sample eddy viscosity distributions computed
without averaging.

- 69 -

M = .96

(a). Pressure coefficient; 0: measured data.

I 2. 4. 6. WD
(b). Surface shear stress T4

Figure 13. Computed results with and without an
--: without averaged-scheme ----: with averaged-scheme




70 -

M = .96


2. 4. 6. X/D
Figure 14. Computed surface pressure with domain S t=18
but different grid resolution in t-direction.
......: 78x40 grid -- 90x40 grid.

71 -

cylinder portion which is a shock/boundary layer interaction

region. In shock/boundary layer interaction regions, the

solutions are more sensitive to grid resolution. The results

show that the 90x40 grid, with more points in the cylinder

portion, has indeed improved the solution characteristics,

especially right after the ogive-cylinder juncture. However,

the computed Cp-distribution of both cases do not reach the

peak value of experimental data on the cylinder part. This

could be attributed to insufficient domain or grid

resolution on the shock/boundary layer interaction region,


Effect of Grid Resolution in Normal Direction

To further test for the required grid resolution in the

normal direction, the 90x40 grid with domain Stot=24 was

solved. The results obtained and shown in figure 15 are

compared with the most accurate solution from the 90x60

grid. Both grid networks have the same grid distribution on

the body surface and the same exponential clustering

function but different numbers of grid points, 40 and 60, in

the -direction. As is evident in figure 15 by the poor

accuracy on the cylinder and boattail portions, the 90x40

grid did not provide enough resolution. The computed drag

force due to pressure is 23.95 which is about 41% less than

that of the 90x60 case. However, the drag force due to skin

friction is 20.28 which has no significant difference from

72 -

M = .96

.4 -

2. 4. 6. X/D
(a). Pressure coefficient; 0: measured data.

M = .96


i 2. 4. 6. /D
(b). Surface shear stress rT
Figure 15. Computed results with Stot=24 but different
grid resolution in C-direction.
----: 90x40 grid --- : 90x60.grid.

73 -

20.53 of the 90x60 case, even though the shear stress

distributions are quite different, as shown in figure 15b.

The flow problem with domain Stot=24 has been solved again

with a 90x40 hyperbolic grid based on hyperbolic tangent

clustering, equations (48), in the normal direction. The

characteristics of exponential and hyperbolic tangent

clustering functions have been discussed in chapter V; the

corresponding grid networks are shown in figure 5. The

computed Cp presented in figure 16a shows that the grid

resolution using 40 points with exponential clustering in

the normal direction is not sufficient; moreover, a grid

which resulted from the use of hyperbolic tangent clustering

yields a better result; it yielded Dp=27.38 while the other

grid resulted in Dp=23.95. These pressure drag forces are,

respectively, 32% and 41% less than the accurate value of

Dp=40.28. The computed shear stress given in figure 16b

indicates that the shear stress distribution is very

sensitive to the pressure field, particularly in the

shock/boundary layer interaction region; however, the

resulting shear drag forces are about the same, 20.05 and

20.28, respectively. It should be pointed out that a grid

with hyperbolic tangent clustering does not always give more

accurate results. In fact the flow problem with a domain

Stot=18 also has been solved on the 78x28 grids. The results

obtained and presented in figure 17 show that a grid with

hyperbolic tangent clustering gives better result on the

boattail portion but the other grid with exponential

- 74 -

M = .96

2. 4. 6. XID
(a). Pressure coefficient; 0: measured data.



2. 4. 6. )WD
(b). Surface shear stress T

Figure 16. Computed results with a 90x40 grid but different
clustering in C-direction.
S : with exponential clustering
-: with hyperbolic tangent clustering.

- 75 -

M = .96

S4. 6. IWD
Figure 17. Computed results with a 78x28 grid but different
clustering in C-direction.
---- : with hyperbolic tangent clustering.
---- : with exponential clustering

- 76 -

clustering yields a much more accurate solution on the

cylinder portion of the projectile.

Importance of Viscous Drag to Pressure Drag

In order to gain further insight as to the relative

importance of the viscous drag to the pressure drag, the

sequence of transonic flows M =0.91, 0.94, 0.98, 1.10 and

1.20 past the SOCBT projectile model at zero angle of attack

have also been solved with the same 90x60 hyperbolic grid

using exponential clustering in the normal direction. We

observed by comparing computed C -distribution with

experimental data, shown in figures 18, 19, 20, 21 and 22,

that the grid which is a very good grid for M =0.96

apparently is not the one best suited for the other flow

cases; however, the results computed are acceptably

accurate. The shear stress distributions computed have shown

that there is a small region of the reversed flow at the

boattail portion for M.=0.91 and 0.94; surprisingly the

computed C on the boattail portion is in excellent

agreement with measured data. It seems to imply that the

Baldwin-Lomax turbulence model can give an accurate solution

for small separated flow regions. The drag forces due to

pressure and shear stress are listed in the following table.

- 77


2 4 6 X/D
(a). Pressure coefficient; 0 : measured data.


a 4



(b). Surface shear

stress rTC

Figure 18. Computed results with a 90x60 grid.

- 78


(a). Pressure coefficient;


C 1

0 2

0: measured data.

2 4 6 X/D

(b). Surface shear stress 14

Figure 19. Computed results with a 90x60 grid.

- 79

(a). Pressure coefficient;


* : measured data.

6 X/D

(b). Surface shear stress rT

Figure 20. Computed results with a 90x60 grid.

, 4


- 80

Cp M = 1.1


2 4
(a). Pressure coefficient; 0 : measly

I M = 1.1






2 4

(b). Surface shear stress Ir

6 X/D

ured data.

Figure 21. Computed results with a 90x60 grid.

6 X/D

- 81

(a). Pressure coefficient;



0 : measured data.

M = 1.2

6 X/D

(b). Surface shear stress rT

Figure 22. Computed results with a 90x60 grid.

- 82 -

Table 2. Drag forces due to pressure and shear stress
at different Mach numbers.

M 0.91 0.94 0.96 0.98 1.10 1.20

Dp 20.81 32.22 40.28 63.78 78.95 98.24

Df 18.35 19.33 20.53 21.88 27.67 32.56

The data in the table indicate that the forces due to

pressure as well as due to skin friction are increased with

increasing Mach number. The ratio of the computed pressure

drag Dp to skin friction drag Df is, respectively, 1.13,

1.67, 1.96, 2.92, 2.85 and 3.02 for flow cases M =0.91,

0.94, 0.96, 0.98, 1.10 and 1.20. These computed data

indicate that the skin friction drag is rather important in

the low transonic regime, and is not as significant in the

higher transonic regime.

The corresponding pressure and Mach contour plots of

M =0.91, 0.94, 0.98, 1.10 and 1.20 are shown in figure 23.

From a series of comparisons, from figures 10 and 23, we

found that initially at M =0.91 only a small disturbance is

felt and expansion waves are observed at ogive-cylinder and

cylinder-boattail junctures. Right after the expansion

waves, the boundary layer thickness decreases due to the

increasing speed and decreasing pressure. As the Mach number

increased to 0.94, shocks are seen to occur on the middle of

- 83 -

M =0.91

0 2 4 6
(a) Pressure contours
p: 0.46-0.84, Ap=0.02

2 4 6

M: 0.85-1.0,

(b) Mach contours
AM=O.O0; M: 1.0-1.25,

Figure 23. Contour plots.


- 84 -

M =0.94

2 4
(c) Pressure contours
p: 0.46-0.86, Ap=0.02

Figure 23. (continued)

(d) Mach contours
M: 0.88-1.0, AM=O.01; M: 1.0-1.2, AM=0.05

- 85 -

M =0.98

0 2 4
(e) Pressure contours
p: 0.46-0.86, Ap=0.02

8 2 4 6
(f) Mach contours
M: 0.88-1.0, AM=0.01; M: 1.0-1.26, AM=0.02

Figure 23. (continued)

- 86 -

M =1.10

2 4
(g) Pressure contours
p: 0.52-0.92, Ap=0.02

2 4
(h) Mach contours
M: 1.0-1.32, AM=0.02

Figure 23. (continued)

87 -

M =1.20

2 4
(i) Pressure contours
p: 0.5-+0.9, Ap=0.02

2 4
(j) Mach contours
M: 1.08-+1.5, AM=0.02

Figure 23. (continued)

88 -

the cylinder and boattail parts. Right after the shock

waves, the boundary layer thickness increases, which is

clearly shown in the Mach contour plot, due to the

increasing back pressure and decreasing velocity. With a

further increase in Mach number, the regions of supersonic

flow have grown and the shocks have moved downstream. At

M =0.98 speed, the computed results show that the shock in

the middle of the cylinder part has disappeared, and the

only shock occurring is at the end of the boattail portion.

This is also evident via the shear stress E, plot in figure

20. A further increase in the freestream condition, to

M =1.10, shows the beginning of a bow shock at the nose and

expansion fans are evident at the surface junctures. At

higher speed M =1.20, the contour plots show the shock as

well as expansion waves are more inclined toward the flow

direction than the case of M =1.10.

From the validity of the Baldwin and Lomax turbulence

model reported in reference 2, we believe that the above

computed results for viscous cases are accurate and good.


To obtain further understanding on the application of the

thin-layer Navier-Stokes approximation, equation (10), to

three dimensional flow problems, the projectile at angles of

attack was investigated.

The three dimensional grid in this study is obtained by

rotating a hyperbolic.planar grid around the axis of the

projectile with a uniform angle increment, A0, in the

circumferential direction. The angle 0 is the angle measured

from the leeward side. It is assumed that the flow problem

considered is symmetric with respect to xz-plane and thus

the grid system generated is only on one half of the

projectile. Two extra planes are included to facilitate the

calculation of the central difference approximation of the

spatial derivatives on the windward and leeward planes.

Based on experience in running the axisymmetric code, the

90x60 hyperbolic planar grid from GRIDGEN with domain 24 is

used to provide a good grid and ensures obtaining accurate

results. This planar grid was rotated to form the three

dimensional grid. The total grid points used becomes 90, 19

and 60 in the streamwise circumferential n and radial (

direction, respectively. The transonic flows of M,=0.91,

- 89 -

90 -

0.96, 1.20 past the projectile model at a=100 were chosen

for numerical simulation.

Flow Cases at 100 Angle of Attack

Flow case at M =0.96

For the flow case M,=0.96, the computed results are in

excellent agreement with experimental data for the surface

pressure distribution. Figures 24 shows the pressure

distribution in the longitudinal direction at different

circumferential planes. The pressure on the ogive portion is

seen to be higher than on the windward side, thus generating

a upward force on the projectile. The boattail region,

however, shows a reversal with the higher pressure now

occurring on the leeward side, thus causing a downward

force. The resulting moment about center of gravity causes a

nose-up behavior of the projectile. This effect can also be

seen in the contour plots of figures 25. The computed

contour plots also show that the shocks on the leeward side

are stronger than those on the windward side. In the Mach

contour plot, it is seen that the boundary layer thickness

on the leeward side is much thicker than on the windward

side, especially after shock which is imbedded in the middle

of boattail part of leeward side. The shock positions on the

windward side appear to be slightly downstream of those on

- 91 -

M= .96

1 -0
4 -------- X
7 -----. -

x/ 2-0
S.. x.*' -------- X

/* /
?4 v

2. 6. WO

Figure 24. Computed surface pressure at various #-stations
of SOCBT projectile at a=10.

- 92 -


3 0
"4 ;''\ M = .96 / 3

6 6

x .4. 3 -
-- x /6 -------- X
d .. 9l ------ -
...-. --.. ."_ .1 A / _0

o 4f .: ^ -

Figure 24. (continued)

- 93 -


* II
0 04

&t *


- 94 -


O +-


-96 .C

'._ _.--
t t
06' 0


1 O/1 /


Full Text
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 EHUCE6PG2_TG0438 INGEST_TIME 2012-02-17T15:56:18Z PACKAGE AA00003797_00001