A self-adaptive computational method for transonic turbulent projectile aerodynamics


Material Information

A self-adaptive computational method for transonic turbulent projectile aerodynamics
Physical Description:
vii, 108 leaves : ill. ; 28 cm.
Reed, Christopher William, 1960-
Publication Date:


Subjects / Keywords:
Fluid dynamics   ( lcsh )
Numerical grid generation (Numerical analysis)   ( lcsh )
Curvilinear coordinates   ( lcsh )
Projectiles -- Aerodynamics   ( lcsh )
Aerodynamics, Transonic   ( lcsh )
Engineering Sciences thesis Ph. D
Dissertations, Academic -- Engineering Sciences -- UF
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )


Thesis (Ph. D.)--University of Florida, 1987.
Includes bibliographical references (leaves 105-107).
Additional Physical Form:
Also available online.
Statement of Responsibility:
by Christopher William Reed.
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 - 021625499
oclc - 18145334
System ID:

Full Text







To Shelley and Drew


The support of Dr. Chen-Chi Hsu has been most encouraging during the

course of this work. His guidance and willingness to provide it are

appreciated. In giving me the freedom to experiment, discover, fail and

succeed I have learned as much about research as I have about the topic.

For this I am grateful.

My appreciation is extended as well to the committee members, Dr.

Leadon, Dr. Kurzweg, Dr. Mikolaitis and Dr. Wilson.

This work was partially supported by the AFOSR. Computer time was

provided by NASA Ames, NSF and the Pittsburgh Supercomputer Center.


ACKNOWLEDGEMENTS . . . . . . . . . . .

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


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


11.1 The Curvilinear Coordinate System . . . .
11.2 A Variational Approach . . . . . .
11.3 Inversion of the Grid Generation Equations .


111.1 The Newton-Raphson Iterative Method . . .
111.2 Modified Solution Procedure ............
III.3 The Adaptive Grid Generation Code .........

IV CONTROL FUNCTIONS . . . . . . . . .

IV. 1 The Basic Form . . . . . . . .
IV.2 Sxmoothing . . . . . . . . . .
IV. 3 Other Control Functions . . . . . .

V FLOW SOLVERS . . . . . . . . . .

V.1 The Governing Equations .................
V.2 Solution Algorithms for the Governing Equations
V.3 Examples on Fixed Grid Networks . . . .


VI.1 The General Procedure . . . . . . .
VI. 2 Orthogonality . . . . . . . . .
VI. 3 Curvature . . . . . . . . .
VI.4 Internal Grid Boundaries . . . . . .








* 37






* 62



VIII RESULTS AND DISCUSSION . . . . . . . . . 81

VIII. 1 Choices for the Control Functions . . . . 81
VIII.2 Results Using the Self-Adaptive Computational
Technique: The Beam and Warming Scheme . . .. .84
VIII.3 Results Using the Self-Adaptive Ccaputational
Technique: The TVD Scheme . . . . . .. .90
VIII.4 Applications to the Projectile Base Flow Problem 97

IX CONCIJDING REMARKS . . . . . . . . . .. 102

APPENDIX . . . . . . . . . . . . .. 104

REFERENCES . . . . . . . . . . . . .. 105

BIOGRAHIICAL SKETCH . . . . . . . . . . .. 108

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



Christopher William Reed

August, 1987

Chairman: Chen-Chi Hsu
Major Department: Engineering Sciences

An adaptive grid generation procedure is developed for viscous flow

problems. The equations governing the adaptation are based on a

variational statement resulting in a set of elliptic equations in which

adaptation can occur independently in each coordinate direction. The

method allows for explicit control of adaptation and orthogonality while

grid smoothness is implicit in the elliptic equations. The adaptive grid

generation equations retain a simple relationship between the control

functions and the grid point spacing, are capable of efficiently

providing the extremely refined mesh in the boundary layer regions and

can maintain a smooth grid network in complex geometries.

A self-adaptive computational technique has been developed by

coupling the adaptive grid generation procedure with a thin layer

Navier-Stokes code and a TVD code. In solving an axisymmetric turbulent

transonic projectile flow problem, the grid network was adapted to the

pressure distribution in the streamwise direction and to the velocity

distribution in the direction normal to the projectile surface. Results

obtained for Mach 0.91, 0.96 and 1.2 indicate that the procedure can

provide good adapted grid networks provided good choices are made for the

control functions.




The study of boundary fitted grid systems has become an important

part of ccmpitational fluid dynamics and, in fact, their development has

been cited as a major pacing item in the progress of computational fluid

dynamics [1]. The use of such systems makes the application of finite

difference techniques versatile and effective in solving complex fluid

dynamic problems in arbitrary domains. Boundary fitted grid generation

can be viewed as the development of a general curvilinear coordinate

system in which the coordinates coincide with the boundaries on the edges

of the physical domain. This feature of the boundary fitted grid systems

facilitates the application of finite difference approximations at the

boundaries and eliminates the need for interpolation between grid points

near the boundary. This is particularly advantageous when the boundaries

are highly curved or have slope discontinuities since the use of

interpolation near these boundaries may have a significant effect on the

solution [2].

In the domain interior, the intersection of each family of coordinate

lines denotes each grid point and therefore defines the computational

mesh. The distribution of the coordinate lines in the physical domain may


be concentrated to provide high resolution in certain areas but will

always form an equally spaced rectilinear coordinate system in the

transformed plane. Thus, once the governing equations are transformed

onto the curvilinear coordinates the finite difference algorithms used to

solve these equations may be developed on the equally spaced rectilinear

coordinate system which is inherent in the transformed plane. Certain

characteristics of the general curvilinear coordinate system have been

shown to affect the accuracy of the finite difference approximations in

the transformed plane [3]. These include orthogonality of intersecting

coordinate lines and the smoothness of the grid spacing. It is also known

that a refined mesh is needed in regions of large physical gradients to

obtain accurate approximations. In fact, poor resolution in these high

gradient regions can be detrimental to solution accuracy and to the

convergence of finite difference algorithms for fluid flow problems [4].

Thus, the development of a good adapted grid network is important for the

successful application of finite difference methods to solve complex

fluid dynamic problems.

The self-adaptive grid generation technique has become an important

area in computational fluid dynamics since it has been shown to provide

good grid networks for the complex flow fields occurring in transonic and

supersonic flows [5,6,7]. The use of boundary fitted curvilinear

coordinate systems with transformed governing equations leads naturally

to the concept of solution adaptive grids. As constraints of computer

storage and CPU time prohibit uniform mesh refinement throughout the

entire domain, the coordinate spacing in the physical domain is varied to

increase resolution in only those regions in which the solution is

changing rapidly. For simple flow problems, when the position of

important solution gradients is known, good adapted grids can be obtained

with conventional techniques. However, in more complex problems the

position of these important regions is not known a priori and then the

development of a proper adaptive grid network is difficult. Solution

adaptive grid generation addresses this problem by continuously updating

the grid network as the solution evolves such that the important physical

gradients are sufficiently resolved as they develop. It is thus an

attempt to efficiently reduce the truncation error by only refining the

mesh in those regions where the error may be large. Also, since the grid

characteristics of smoothness and orthogonality also affect the

truncation error, a good adaptive grid generation technique should

include the enhancement of these characteristics as well.

The general approach of most adaptive grid generation schemes is

based on minimization techniques. A measure of each desired grid

characteristic is defined, and the governing equations for the grid are

obtained by minimizing the integral of these measures over the domain. In

many instances adaptation is required in only one coordinate direction.

Dwyer et al. [8], for instance, have developed an equidistribution law to

control the grid spacing along one coordinate direction. The spacing was

adapted to a flame front in a two-dimensional combustion problem and the

improved resolution eliminated spacial oscillations in the flame front.

Gnoffo [9] has developed another one dimensional adaptive scheme in which

tension springs are assumed to connect the grid points along a

coordinate. The spring constants became large when increased resolution

was required, and the points clustered in those regions in an attempt to

minimize the potential energy of the spring system. The application of

this technique to supersonic flow over a Viking Aeroshell successfully


resolved a separated shear layer by adapting to the velocity gradient.

One dimensional adaptation can be extended to higher dimensions by

successive adaptation in each coordinate direction. Nakahashi and Deiwert

[7] have used the spring analogy in such an approach and added torsional

springs that connect intersecting coordinate lines to control

orthogonality. They solved a transonic viscous airfoil problem in which

the grid was adapted to the density gradient in the streamwise direction

to resolve shocks and adapted to the velocity gradient normal to the

airfoil surface to resolve the boundary layer. The one-dimensional

approach to multi-dimension adaptation has the advantage of efficiency

since the grid can usually be obtained using a marching solution

algorithm. Another advantage is the independence of adaptation in each

coordinate direction. However, it is difficult to maintain smoothness in

such approaches, and highly skewed grids may result [10].

Rai and Anderson [11] have used a different approach for two

dimensional adaptation. In their scheme each grid point induces a

velocity, either repulsive or attractive, on every other point. By

choosing some flow gradient to indicate the need for point clustering,

each point with a gradient larger than the average gradient attracts

points, and those with values below the average repel points. By summing

the induced velocities from each point, the velocity of each point is

determined and can be integrated over a time step to determine the new

point location. In applying this scheme to the two dimensional linearized

viscous Burger's equation, they adapted the grid to the derivative of the

dependent variable in each direction and showed that such an approach

will increase solution accuracy.


Saltzman and Brackbill [5] have developed an adaptive grid generation

scheme based on a variational approach. A functional is defined to

measure each grid characteristic and the minimization of these

functionals results in a set of partial differential equations that

govern the grid point spacing. With the proper choice of parameters the

equations are elliptic, which guarantees a one-to-one mapping and results

in a smooth grid. They solved an inviscid supersonic flow problem by

adapting the grid cell size to a function of the pressure gradient to

capture shocks. However, it is not possible in this approach to adapt

the grid spacing independently in each coordinate direction, a feature

that is prohibitive if the spacing in different coordinate directions

vary by a large amount.

The adaptive grid generation technique presented here is similar in

approach to that of Brackbill and Saltzman but is developed for

applications to viscous transonic flow problems. The solution of these

problems usually contains shock structures aligned with one coordinate

direction and boundary shear layers parallel to a streamwise coordinate.

Also, the grid point spacing can vary by orders of magnitude along

different coordinates and it is thus necessary to adapt the grid

independently in each coordinate direction. Therefore, a functional is

defined for independent adaptation in each coordinate direction such that

the resulting equations are elliptic. The use of elliptic equations in

grid generation offer some important advantages that motivate their use.

They constitute a boundary value problem and can thus be used in any

domain. They guarantee a one to one mapping and also, grid smKothness is

inherent in the equations. In addition to the equations governing the

grid, there are two other topics that require attention when developing


an adaptive grid generation technique. First is the definition of the

control functions, which are the part of the grid generation equations

that dictate the need for adaptation. They are dependent upon the

solution and thus provide the link between the grid generation and the

governing equations. Also, a method is needed to account for the effect

of the moving grid on the solution. Since the solution is defined at

each point, any movement will alter the solution. The proposed adaptive

grid generation technique is used to adapt the grid for a transonic

axisynmmetric projectile flow problem which is solved using both an

implicit factorized algorithm and a TVD scheme for the thin layer Navier

Stokes equations. The constraint of axisynmmetric flow reduces the domain

to two independent spatial variables, thus the grid generation equations

are developed in two dimensions. The technique, however, can be readily

extended to three dimensions.


11.1 The Curvilinear Coordinate System

The use of boundary fitted curvilinear coordinate systems in solving

equations on arbitrarily shaped domains was originally motivated due to

problems in applying finite difference approximations at curved

boundaries. Because curved boundaries do not in general coincide with the

rectilinear Cartesian coordinates, the grid points defined along these

coordinates do not usually lie on the bounding curve. It therefore

becomes difficult to accurately represent the boundaries [2] and special

procedures must be used in those regions [12]. A solution to this

problem, which has become an effective and versatile approach, is to map

the arbitrarily shaped physical domain (x,y) into a rectangular

computational domain ( Y7). Such a mapping for a typical projectile flow

domain is shown in Figure 2.1. As indicated, the curved C and 17

coordinate lines coincident with the boundaries in the physical domain

map onto the edges of the rectangular computational domain. The boundary

grid points, being defined by the coordinate intersections, are

coincident with the boundary in the computational plane and thus no

special treatment for the boundary conditions is required. Another

advantage is that virtually any set of nonuniformally spaced curvilinear

coordinates in the physical domain map into an equally spaced regular


(1.JM) j (,.,1M;


________ _________________________ _________i

___ ___ ____ ________ __ __ j__ [___'


Figure 2.1 Curvilinear coordinate system

'1M.! )



mesh in the computational domain. Furthermore, once a successful finite

difference algorithm is developed on the computational domain, any

physical domain may be mapped into this region and solved with the same

numerical algorithm.

The task of numerical grid generation is to define the transformation

between the Cartesian and computational coordinate systems which can be

viewed as the development of a curvilinear coordinate system in the

physical domain. The general mapping between the two coordinate systems


C = (x,y) (2.1.1)

1= 7 (x,y)

Owing to the discrete nature of finite difference approximations, it is

sufficient to define a finite set of coordinate lines, the intersections

of which define the grid points in the physical domain. Thus for each

grid point (x,y) in the physical domain, there are two curvilinear

coordinates ( t) that intersect at that point and designate its

location. It is convenient to view this mapping in the computational

domain, where for each pair (C,t7) there corresponds a point (x,y) in the

physical domain. This correspondence can be conveniently expressed as the

inverse mapping

x = x(C,v) (2.1.2)

y = y(,?I)


It is advantageous to designate the curvilinear coordinates with integer

values. Finite difference expressions in the computational plane are
simplified since A and At become unity and the mapping of eqs. (2.1.2)

will appear as a two dimensional array with integer subscripts, creating
a naturally ordered way to store the grid points. Using the inverse
mapping, the location of any coordinate line in the physical domain can
be determined simply as the set of points obtained by holding one of the
arguments in the mapping constant. The integer indicies i and j are used
here to designated each grid point. Thus the ith point along a
coordinate is i and the jth point along an n coordinate is 1j. For
simplicity the notation (xi,j,yij) is used to represent the point
(x(i,qj),Y(i,?j)) and the total number of points in each direction
and r are IM and JM, respectively. For the projectile problem
considered, the and i coordinate lines will always run in the
streamwise and normal directions as indicated in Figure 2.1.
Once a mapping is defined, in a discrete sense, the governing
equations can be obtained in the computational plane by transforming them
onto the curvilinear coordinate system. Derivatives with respect to the
Cartesian coordinates can be expressed using eqs. (2.1.1) and the chain
rule as

[ax a
{ } I { } (2.1.3)
| a |y a |7
Iay Jy [Qr

where x, y, nx, and ?7y are the metric coefficients of the
transformation. They appear as constants in the transformed equations and
depend solely on the mapping of eqs. (2.1.1). Substitution of eqs.


(2.1.3) into the governing equations makes and n the independent

spatial variables and now a numerical algorithm can be developed using

finite difference approximations on the equally spaced computational

domain. Since the metric coefficients depend only on the mapping,

virtually any domain may be mapped into the computational domain and

solved using the numerical algorithm.

The metric coefficients can be evaluated for a numerically generated

transformation by considering the metrics of the inverse transformation.

The relationship between the metric and inverse metrics is, for two


{x y i = 1 Y -Y J (2.1.4)
y 17y -x X

where J is the Jacobian of the transformation,

J = x y7 xy Y (2.1.5)

The derivation of these equations can be found in reference [13]. Since

and r are the independent variables for the inverse metrics, the

inverse metrics can be evaluated numerically in the computational domain

and subsequently used to determine the metric coefficients.

The primary disadvantage of using boundary fitted curvilinear

coordinate systems is the appearance of additional truncation error terms

in the finite difference approximations on the computational plane.

Mastin [3] has shown, for instance, that the truncation error T for a

second order finite difference approximation of the first derivative is,

in the ccxrputational domain

T = -x2 f xxx fxx (2.1.6)
T=6x fxx 2

The first term is the usual truncation error term containing the grid

point spacing x and the third derivative of the function. The second

term is due to the varied spacing along the $ coordinate lines in the

physical domain. It is thus possible to increase the truncation error if

the coordinates spacing is not smooth. They have also shown that

truncation error for a first derivative varies inversely as the sine of

the angle between the curvilinear coordinates and therefore a lack of

orthogonality can also increase the truncation error. Orthogonality is

also important near solid boundaries for turbulent flow problems because

turbulence models often are based on information in the normal direction.

In the context of a coordinate transformation, the adaptive gridding

can be viewed as an attempt to reduce the truncation error by adjusting

the transformation metrics. This is accomplished by judiciously defining

the curvilinear coordinates such that they vary smoothly, are nearly

orthogonal, and concentrate in regions where the solution is changing


11.2 A Variational Approach

The choice of a method to generate the curvilinear coordinates is

arbitrary in that it has no dependence on the partial differential

equations that are to be solved. Of the many methods available, however,

the variational approach is attractive for purposes of adaptive grid

generation because it can provide a simple, explicit means of

incorporating desired properties into a grid generation scheme and can

result in a set of elliptic equations which govern the grid. To develop
an adaptive grid generation scheme using variational techniques, a
measure of each grid characteristic is defined and its integral is taken
over the physical domain. If the measures are defined such that they
measure a deviation from the desired property then the coordinates and
,7 that minimize the integral will yield a grid with the desired
characteristics. As a minimization problem, the Euler-Lagrange equations
can be applied to obtain a set of partial differential equations, the
solution of which defines the curvilinear coordinates.
The first of three functionals used here is

Ip = J v8dxdy (2.2.1)

which measures the adaptation of the grid spacing in the coordinate
direction to the control function P. When P is small, the quantity V
must also be small in order to minimize the integral; a small value of v7
corresponds to large grid spacing. Consequently, when P is large the grid
spacing along the coordinate will be small. Similarly, the integral

Iq= J Q Vndxdy (2.2.2)

is a measure of the adaptation of the grid spacing in the P direction to
the control function Q. The third functional

Io = (v.vi1)2 dxdy (2.2.3)


measures the orthogonality of the and q coordinates. As will be shown

later, grid smoothness can be controlled through the definitions of the

control functions. These three integrals can be combined into one total

functional IT

IT { V'V + Vn.Vn + (V7.V?7)2 } dxdy (2.2.4)

where x is a parameter that weighs the relative importance of adaptation

and orthogonality. Before applying the Euler-Lagrange equations for and

q, it is beneficial to scale each term in the equations so that they

remain of the same order of magnitude. An ordering analysis shows that

the first term in eq. (2.2.4) is of order (C2/P), the second is of order

(C2/Q) and the orthogonality term is of order (C4/L2) where C is a

computational length scale and L is a physical length scale. It is ccmnnon

practice to use highly stretched grids in transonic flow problems to

provide extremely small grid spacing in the viscous sublayer region. The

spacing along the normal coordinate can increase by five orders of

magnitude, consequently the terms in eq. (2.2.4) may vary in a similar

way. It is therefore beneficial to scale the terms locally rather than

globally since a global scale cannot reflect the size of the term

throughout the domain. The local scales Ap, Aq, and Ao used here are

Ap=PL Aq=QL Ao=JL (2.2.5)

where PL and QL are the local values of the control functions and JL is

the local Jacobian and is of order (L2/C2). Although they vary throughout

the domain they are considered constants during the application of the


Euler-Lagrange equations. This violates the variational principal, but it

has been shown [13] that the use of local scaling will produce better

grids when there are large changes in the grid spacing. After

substituting in the scales, the total functional IT becomes

IT= J ( PL v + QLVn .V + AJL (V'.Vq)2 )dxdy

After applying the Euler-Lagrange equations for and P to the total

functional IT while holding the scales constant, the following set of

partial differential equations is obtained:

xx + yy + Vp.V + A(n2 xx+2qxqy~xy+?yy +

(2Cxx+Cyqy) Y7xx + (t x y+ y x) xy+(x x+2yry) yy)=O


'lxx + 'lyy + V.Vn + ((2xtx+yty)xx+( x~y+qy~x)xy +
(Qx'x+2Cy'y) xy + x'2xx+2x~ytjxy + C217yy)=0

The subscript L has been dropped from the local scales in the

differential equations since they are no longer needed to indicate their

assumed invariance. The and ri coordinates that satisfy these equations

will yield a grid with the desired properties defined in the functionals

of eqs. (2.2.1), (2.2.2) and (2.2.3). It should be noted that there was

no explicit functional to maintain grid smoothness. There are two types

of smoothness that should be differentiated. One type, which can be

controlled by the control functions, is smoothness of the grid point


spacing along a coordinate line. If the control functions change rapidly

along a coordinate, the grid spacing will consequently change rapidly.

This smoothing thus can be governed through the definition of the control

functions. Another type of smoothness, which is inherent in the elliptic

equations, is that between the grid point spacing of adjacent

coordinates. If the grid point distribution along two adjacent

coordinates varies, the elliptic equations will tend to dampen the


Equations (2.2.7) govern the grid coordinates and n and will be

elliptic as long as the parameter A is not too large. They constitute a

boundary value problem and it is therefore necessary to define the

coordinates along the boundaries. It is also necessary to adapt the

boundary points in a consistent manner with the interior points so that

the boundary points remain aligned with the interior points. For this

purpose, a one dimensional functional analogous to equation (2.2.1) or

(2.2.2) is defined to measure adaptation along a boundary coordinate,

IS= f s ds (2.2.8)

where s is the arclength along the boundary coordinate. There is no

consideration for orthogonality and scaling is unnecessary. Applying the

one dimensional Euler-Lagrange equation yields a second order

differential equation,

Sss-Cs Ps/P=O (2.2.9)

Note that a similar equation can be obtained for boundary adaptation

along a r coordinate by replacing in eq. (2.2.9) with ? and P with Q.

This one dimensional equation, together with eqs. (2.2.7), are the


equations that govern the grid coordinates in the physical domain.

11.3 Inversion of the Grid Generation Equations

In order to solve the grid generation equations, they are transformed

onto the curvilinear coordinate system. This is done by inverting eqs.

(2.2.7) and (2.2.9) to make x, y, and s the dependent variables. The

curvilinear coordinates then become the independent variables and the

grid generation equations are solved on the computational domain which is

equivalent to defining the inverse mapping of eqs. (2.1.2). All the

advantages of using the computational domain to solve the governing

equations discussed previously will be realized here as well.

There are two approaches to inverting the equations. One is to

transform the functional IT onto the computational domain and then apply

the Euler-Lagrange equations for x and y. Another approach, used here is

to derive the corresponding expressions on the computational domain for

each term in eqs. (2.2.7) and (2.2.9) and then to substitute them into

the equations. The first term in eq. (2.2.7) can be written using eq.

(2.1.3) as

xx = ( x s-- + "X -- ) X (2.3.1)

Using eq. (2.1.4) to transform the metrics this expression


Y,18 yC a Y'i
Cxx =( - ) (2.3.2)
J a9 J -9i J

and after differentiating, the relationship is found to be

-Y__ 2 x 2 2
1XX- J 2(y2x-2ygyx -y+2x)+-(y~y -2ygyEy7+yy2


Similar expressions can be derived for each of the other second

derivative terms in eq. (2.2.7), they are listed in Appendix A. All the

first derivative terms, of course, can be transformed using eqs. (2.1.4).

By substituting these expressions into eqs. (2.2.7) and collecting like

terms the equations for adaptive grid generation become, in the

cciputational domain,

alxE + a2x n + a3x, + blyE + b2Yn + b3Yq = S1

clxC + C2X, + c3Y, + dly + d2Yn + d3Ynq = S2

in which

a, = -y7(c( + A'#2) 2A'yeaf

a2 = 2y (# + A'f-y) + A'y (4,2+J2)

a3 = -yi (I + A'72) 2A'y af3

bI = x,7 (a + A',2) 21yct(

b2 = -2x,, ( + A'7-y) A'x (4#2+J2)

b3 = x,(7-y + A'-y2) + 2A'y73fi

CI = yE(a + A'a2) + 2A'ynap

C2 = -2y (6 + A'), ) A'y, (4,32+J2)

C3 = y$(-y + A',32) + 2A'yn7p

di = -xa (Q + A'a2) + 2A'yaf3

d2 = 2xE(62 + A'ap) + A'y(462+J2)

d3 = -Y (7 + A'\2) + 21'yr7


S= x2 + Yb2

Y = x2 + y2

= x Y + xy

A' = A/J


S1= +--
P p

Q Jf QJ7
S2 -= +

It is important to note that in inverting the equations the derivative of

P in the Y direction appears even though P was intended to control the

spacing in the direction. This occurs because the contravariant base

vector v actually indicates the spacing in the direction normal to lines

of constant rather than in the direction. If the grid is orthogonal,

the two directions, v1 and V, are equivalent and the ) derivative will

not have any influence since 6 will vanish. It has been found in

numerical experiments that eliminating this term from the equations has a

negligible effect and consequently it is dropped from the equations. The

source terms S1 and S2 in eqs. (2.3.5) then become

S1 =-J

S2 -=Q


The one dimensional equation for the boundary adaptation, eq.

(2.2.9), is also transformed in a similar manner. The second derivative

term can be shown, using the chain rule, to transform as

Sss - 3--


Substituting eq. (2.3.7) into eq. (2.2.9) yields the transformed one

dimensional equation for adaptation along C boundaries


s6 + s PC/P =0

and for adaptation along boundaries coincident with 7 coordinates

snn + SnQn/Q = 0


These equations can be solved directly for the grid point spacing along a

coordinate. For instance, integrating eq (2.3.10) once yields the grid

point spacing

S = C/Q


where c is a constant of integration and can be shown to be

C ---
f da?


where ST is the total arclength along the coordinate in the physical

domain. However, since the solution in the interior will require an

iterative algorithm, the differential form, eqs. (2.3.9) and (2.3.10) are


retained and solved iteratively so that the boundary points remain

consistent with the interior points as the grid evolves.


III.l The Newton-Raphson Iterative Method

Equations (2.3.3), (2.3.9) and (2.3.10) constitute a set of nonlinear

elliptic partial differential equations. In order to solve these

equations numerically they are approximated with second order central

difference expressions which results in a set of coupled algebraic

equations for each grid point. The finite difference approximation for

the equations at grid point, (xij, Yi,j) become


b, (yi+l, j

cl (xi+l,j

dl (yi+l, j

-2xi,j + xi-l,j) + a3(xi,j+1 -2xi,j + xi,j-I) +

-2yij +Yi-l,j) + b3(Yi,j+l -2Yi,j + 2Yi,j-l)-FI=O0


-2xi,j + Xi-i,j) + c3(xij+l -2xij +xi,j-l) +

-2Yi,j +Yi-l,j) + d3(Yi,j+l -2Yi,j +Yi,j-l)-F2=0


F1 = 2 aT1 b2T2

Q J7
Q?2 1-- a i d2T2
F2 a2T1 d2T2


T1 = (Xi+l,j+l + Xi-l,j-l1 Xi+l,j-l xi-l,j+l)/4


T2 = (Yi+l,j+l + Yi-l,j-l Yi+l,j-l Yi-lj+l)/4

where the source terms and the coefficients all contain first derivatives

which when approximated with central differences do not contain the point

xij, Yi,j explicitly. This set of equations is solved using a Newton-

Raphson iterative scheme. The initial guess for the grid point locations

will not satisfy the finite difference representation of the equations

and therefore a non-zero residue will exist

(Rl)i,j = eq (3.1.la)

(R2)i,j = eq (3.1.1b)

where 2 indicates the iteration number. An estimate of the residue at the

next iteration can be obtained by expanding the residue in a Taylor

series about the current residue and retaining only the linear terms,

2+1 2 a(Rl) i,j a(RI)i,j
(Rl)i,j = (Rl)i,j + Axi,j -- + Ayi,j --


S2 2
2+1 1 a(R2)i,j 8(R2)i,j
(R2)i,j = (R2)i,j + AXi,j3 + Ayij,

Where Ax and y are the changes in the grid point locations
where AX and Ay are the changes in the grid point location


AYi,j =Yi,j-Yi,j


The derivatives of the residues with respect to the variables Xji,j and
Yi,j can be obtained by differentiating eqs. (3.1.1) and are

(R )i ,j

a (R2)i ,j

An approximate
residue at the

= -2(al+a3),

= -2(c1+c3),

value of Axij and Ayi, j
next iteration to vanish

I(Rz) ]
.(R2) l,j

a (Rl)-l,j

a (R2)1,j

a (Ri,j


= -2(bl+b3)

= -2(dl+d3)

can be obtained if we require the


a (R2)i ,j



Substituting eqs. (3.1.5) into eqs. (3.1.6)
Axij and Ayij yields

f Axi,j 1 2 r (dl+d3)
I AYi,j J D -(c+c3)

and solving for

-(bl+b3) 1*(Rl)i,,j
(al+a3) (R2) il,j


D = 4((al+a3)(dl+d3)-(cl+C3)(dl+d3))



Knowing Axij and Ayij, the new point location can be calculated by
2+1 2+i
solving eqs. (3.1.4) for xij and Yij.

This same procedure can be applied to the one dimensional equations

for adaptation, eqs. (2.3.9) and (2.3.10). The residue for eq. (2.3.9) is

(Rs)i = si+1 2si + Si-1 + sP/P

and the equation for the new point location becomes

1 1i/2i




The numerical algorithm consists of updating the value of each point

on the boundary with eq. (3.1.11) and of sweeping through the domain and

updating each interior point using eq. (3.1.7). The criterion used here

to indicate convergence is

Rmax < 6


where Rmax is the maximum residue in the interior weighted by the


AX2 j + AYiJ )1/2
Rmax= max( '--' )1/2


Currently, the criterion 6=0.005 is used.


In the original iterative procedure, the most recently updated values

for each grid point location were used in updating the remaining points

and thus the procedure resembles a Gauss-Seidel iteration. This approach,

however, is not suitable for vector processing and thus cannot benefit

from the increased processing speed. Therefore, another approach is used

in which all the grid points along a coordinate line are updated in a

Jacobi iterative fashion. The new grid point values along that coordinate

are then used in updating the points along the next coordinate line. Thus

the technique can be viewed as a half Jacobi, half Gauss-Seidel

iteration. This technique, which is suitable for vector processing,

reduced the computer processing time for each iteration by 40 percent.

Furthermore, only a few more iterations were required for the same

convergence criterion in comparison to the Gauss-Seidel iteration, thus

the effective speed up in processing time is approximately 40 percent.

III.2 Modified Solution Procedure

As pointed out earlier, the grids used for viscous transonic flow

problems contain highly refined grid spacing in the direction normal to

the surface (P? direction) in order to resolve the viscous sublayer. This

spacing results in grid cells with aspect ratios, up to 105 in the

present applications, near solid boundaries. Iterative solution

procedures for elliptic equations become inefficient for such large

aspect ratios since the motion in one direction will be limited by the

small spacing in the other direction. The effect of the aspect ratio on

an iterative process can be ascertained by considering the simple grid

cell shown in Figure 3.1 which for convenience is rectangular and aligned

with the Cartesian


LX wu)__ )


Figure 3.1 High aspect ratio cell

f *1~1

________ 1 _____ -- -. ______

- __ __ ii-
irni - - ___ -~
I ~ i iii--,

Figure 3.2 Initial grid network

y,1 4,

O, b )





coordinates. For this example the aspect ratio AR is equal

to (1/b). Assume for purposes of demonstration that the two

terms (Q /Q) and (P/P) are non-zero at the center point and

will therefore cause the point to move. We can obtain the

changes in the point location by evaluating eqs. (3.1.7) for

this simple case. Many terms for this example are zero


yC =0
y =0

y 17 =0
x =0


yr =0

The coefficients become














and upon substitution into eqs. (3.1.7) yield

P b
Ax'j 2P 1 l+b

Q1 1
Ay, ( )
-2Q l+b


Writing these in terms of the aspect ratio AR, its effect on the grid

point motion becomes clear. If we compare the solution for a square cell

, i.e. AR=1,


AX (.)
)xi 2P 2

Aij=- ( )
Yi,' 2Q 2

with that for a cell of large aspect ratio AR >> 1,

Pe( 1
Axjij 2P l+AR

oi7 AR
Ay, =1 ( AR )
Yij 2Q 1+AR

it becomes evident that the reduction in the grid point motion in the x

direction for large AR is proportional to (1/AR). For aspect ratios of

the order (105) this result becomes prohibitive, rendering the iterative

procedure extremely inefficient.

An example has been formulated to demonstrate this deficiency. In an

initial grid network, as shown in Figure 3.2, the points along the

coordinates are clustered in one region and the points along each v

coordinate are clustered near the bottom surface to form the highly

stretched spacing typical of grid networks used in transonic flow

calculations. The spacing at the surface is approximately 10-5 times

that in the upper region of the grid. In the adaptive grid generation

equations, the control function P is set equal to a constant so that the

points will redistribute themselves with equal spacing along the

coordinates. Using a technique described in Chapter IV, the control

function Q is defined so that the highly stretched spacing along each


coordinate will be retained in the solution to the adaptive grid

generation equations.

Using the Newton-Raphson iterative solution procedure, the solution

is advanced 20 iterations which results in the grid shown in Figure 3.3.

As can be seen, the point spacing along the coordinates is

approximately equal for the coordinates sufficiently above the surface.

However, in the region of the high aspect ratio cells, the movement of

the points in the direction is retarded by the extremely small spacing

in the i direction.

A modification has been made to eq. (3.1.11) which updates the points

along the bounding coordinate in formulating this example. The

movement of points using the one dimensional equation is independent of

the interior points and thus is not restricted by the high aspect ratio

cells. This creates an inconsistency at the boundary since the points

adjacent to the boundary are constrained by the high aspect ratio cells.

Equation (3.1.11) therefore has been modified as

A SA 1
ASi = (Rsi)i (-- ) (3.2.6)

where AR* is the aspect ratio of the cell adjacent to the ith point along

the boundary coordinate. This change only effects the rate of

convergence of the iterative solution procedure for the boundary points

so that their movement will remain consistent with that of the interior

points adjacent to the boundary during the convergence process.

In order to alleviate the deficiency of the above solution procedure,

a technique has been developed which is based on solving for a reduced

grid in which many of the points in the 7 coordinate direction have been


\ \ \ I I


\L/ \ II


Figure 3.3 Grid network after 20 iterations

-- __ -_
- 7M
- 7=

Figure 3.4 Reduced aspect ratio cell



removed to decrease the maximum cell aspect ratio. Figure 3.4 shows a

diagram of a typical grid for viscous flow problems containing large

aspect ratio cells near a solid boundary. A reduced grid can be formed by

removing all the coordinate lines denoted with the dashed lines,

resulting in a grid with fewer points along the r7 coordinate lines which

are designated by the intersection of the solid lines. Currently, enough

points are removed along the q coordinate lines such that the minimum

spacing between coordinate lines is greater than 0.04. The typical

spacing along the boundary is of order (0.1) which results in cell aspect

ratios of order (1). Upon convergence of the solution for the reduced

grid, all the points removed are reinserted along r coordinate lines

using the same one dimensional equation that controls the spacing along

the boundaries.

It is necessary in this process to adjust the control function Q at

the remaining points to account for the spacing requirements of the

points removed. A simple procedure for this purpose can be derived by

using the equation for one dimensional adaptation, eq. (2.3.10). The

procedure will also be adequate for the two dimensional adaptive grid

generation equations. If we write explicitly the finite difference

approximation for eq (2.3.10),

sj+-j+1j-l Qj+I-Qj-l
Sj+l-2sj + sj-i + (- 2 )( 2 )/Qj=0 (3.2.7)

we can regroup the terms in the form

QJ + (--2- -)
(Sj+l-Sj)=(sj-sj-l) (- Qj+---- ) (3.2.8)
Q (Qj+I-Qj-l


which, after using forward difference expressions to indicate the spacing

between two adjacent points can be written

ASj = ASj-1 Cj

Qj+i-Qj -i
QJ + ( 2 )
Cj= 2)
Qj ( --- -)Q



Now, as shown in Figure 3.5, if the ith point is removed it will be

necessary to adjust the control function such that the new relationship

is satisfied

ASj+i = AS*-1 C3-I


where Cj-1 is the modified control function and AS3-i is the increased

grid point spacing

ASjl = ASj-|- + ASj


The correct form of Cj-l can be obtained in the following derivation

ASj+1 = ASj_- Cj-i

ASj Cj+1 = (ASj-l+ASj) Cj-1

Asj-1 Cj Cj+1 = ASj-1 (1+Cj) Cj-1


cj511 c- +
1+ _1
I+Cj J-


Once C3-_ is calculated using eq. (3.2.13), the modified value of Qj-

can be obtained by solving eq. (3.2.10) for Qj-i

Q*- = (- ) (Qj+l-Qj-2) (3.2.14)

This procedure is applied for each point along an j coordinate that is

removed. It should be noted that this procedure does not affect the final

grid, the points remaining in the reduced grid will obtain the same

location as if all the points were used in the iterative solution

procedure. The removed points are, once the solution for the reduced grid

is obtained, inserted back along the r7 coordinates lines. This is done

by approximating the r? coordinate with straight line segments between the

points of the reduced grid and then reinserting the removed points along

each segment with spacing dictated by eqs. (2.3.10) using the original

control function Q. If necessary, the use of a higher order approximation

to the r coordinate lines, such as a cubic spline, can be used to obtain

a better representation of the n coordinate lines. In using this

procedure with the two dimensional grid generation equations, the

solution may vary somewhat since the equations are coupled. Also, the

effects of orthogonality may alter the results in some areas.

In order to show the advantage of this procedure, the initial grid

network of the previous example is advanced 20 iterations using the

reduced grid technique. Figure 3.6 shows the initial reduced grid in

which the points responsible for the high aspect ratio cells are removed.

After modifying the control function Q, the solution to the adaptive grid

generation equations is advanced 20 iterations resulting in the grid

shown in Figure 3.7. As can be seen, the points along all coordinates

Figure 3.5 Corresponding grid point spacing

___Figure 3.6 Initial reduced grid network_ __
Figure 3.6 Initial reduced grid network


Figure 3.7 Reduced grid network after 20 iterations

Figure 3.7 Reduced grid network after 20 iterations

Figure 3.8 Final grid network after reinserting
grid points


are approximately equally spaced and the spacing along the r coordinate

lines in the initial reduced grid is retained which indicates that the

modified control function is correct. Upon inserting the removed points

along the ri coordinates using eq. (2.3.12), the final grid network,

Figure 3.8, is obtained, which is the same grid that would have been

obtained in Figure 3.3 after a large number of iterations.

III.3 The Adaptive Grid Generation Code

A computer code base on the iterative solution procedure described in

the preceding section has been developed and is referred to here as the

Adaptive Grid code. As the solution process is iterative, an initial

guess for the grid point locations is required and the flow variables for

which the grid network is to be adapted are needed for calculation of the

control functions. These are considered as input to the Adaptive Grid

code. The initial guess for the grid network can be a grid network

obtained by any conventional grid generation technique or a previously

adapted grid network. Values for the flow variables at the grid points

are obtained using the flow solvers which are described in Chapter V.

The sequence of calculations in the Adaptive Grid code proceeds as


1. Input an initial grid network and the flow variables at each
grid point.

2. Calculate the control functions P and Q.

3. Remove points to form the reduced grid network and
calculate the modified control function Q*.

4. Solve the grid generation equations for the reduced grid network
using P and Q* .

5. Reinsert the removed points along the n coordinate lines using
the control function Q which then results in the adapted grid


The methods used to calculate the control functions from the solution

variables is described in Chapter IV.


IV. 1 The Basic Form

It is the role of the control functions in the equations governing

grid adaptation to dictate the need for increased resolution and, as

indicated in section 11.1, they should be judicially chosen so that the

resulting grid resolution reduces what would otherwise be large

truncation error. For instance, Pierson and Kutler [14] chose for the

control function, when solving a one dimensional inviscid Burger's

equation, the third derivative of the dependent variable which is also

the leading truncation error term in their finite difference

approximation to Burger's equation. Their one dimensional adaptive grid

scheme was thus an attempt to minimize a measure of the truncation error

over the domain. In more complex problems, however, with more dependent

variables and in higher dimensions, the truncation error expressions are

complex and cumbersome to calculate. Therefore, in practice, an

understanding of the physical problem and experience guide the choices

for the control functions. The general control function considered here


P = 1 + Jpf
Q = 1 + 7qg


where 7p and -Yq are parameters and f and g are sane derivatives of the

dependent variables scaled to range between zero and one. In the

following developments, the terms W,-y and h will be used to indicate

either P,7p and f or Q,-q and g. By evaluating numerically the solution

to the one dimensional equation for adaptation, the following expression

equating the control function to the grid point spacing between points

(i) and (i+l) along a coordinate line is obtained,

ASi = +-yhi (4.1.2)
j l+-yhj

After evaluating this expression for the maximum and minimum spacing

corresponding to h = 1 and h = 0 respectively, and forming their ratio

( m 1 + -y (4.1.3)

it can be seen that the specification of 7 dictates the ratio of maximum

and minimum spacing along a coordinate line. The parameter 7-y also

influences the smoothness of the resulting grid point distribution. If it

is chosen as zero, equal spacing as indicated by eq. (4.1.2) would

result; as 7 is increased, the spacing must change more rapidly to obtain

the varied spacing. Following the work of Nakahashi and Deiwert [7],

another parameter a has been introduced to obtain more control over the

spacing. The general control function now becomes

W = 1 + 7ha


By prescribing the minimum spacing and 7-y, eq (4.1.2) can be written as

1 ST (415)
(ASi)min (l+jh ) -1+ (4.1.5)

in which a is the unknown quantity and can be determined by a root

finding technique. Currently, the Newton-Raphson finding method is

employed. Thus, both the minimum and maximum grid point spacing along a

coordinate can be specified by specifying the two quantities, (Asi)min

and -Y. It should be noted that the prescribed values may not be realized

in practice due to the coupling effects of the two dimensional adaptive

grid generation equations. Also, the constraint of orthogonality may

alter the results in same areas.

IV. 2 Smoothing
The choices for the functions f and g will of course depend on the

problem being solved. However, the functions will in general, contain

derivatives of the solution which are approximated numerically and it is

therefore useful to smooth the control functions to eliminate any

spurious data. Also, since smoothing is most effective when the function

changes most rapidly, it will help reduce any rapid variations in the

control function and subsequently smooth the grid point spacing. The

method of smoothing used here, which is similar to that used in reference

[5] is

Wi+ij+Wi-ij ^J i=2,IM-I
Wij= WiJ + 2 Wij j=I,JM


f Wi,j+l+Wi-lj I i=lIM
Wij= WilJ + -W j=2,JM-I
'j=W~ 2 1,I

where the value of W used here is 0.4. The control function at a point is

replaced by a weighted average of itself and the two adjacent points

along a coordinate line. Currently, two consecutive sweeps of this

algorithm are made over the domain for each of the control functions, P

and Q. In this algorithm there is no accounting for the varying spacing

between points and thus this smoothing is in the computational plane.

IV. 3 Other Control Functions

In many instances a certain point distribution in one coordinate

direction is known to yield good results and it would be advantageous to

maintain this prescribed distribution in that direction while allowing

adaptation to the solution in the other directions. For instance, in some

of the viscous transonic flow problems solved here, it is known that a

distribution in the direction normal to the surface based on an

exponential function will yield good results provided enough points are

used. A procedure to incorporate specified point distributions into the

adaptive grid scheme is easily developed by using eq. (2.3.10).

Consider some given grid point distribution s in the n coordinate

direction which is to be reproduced by the adaptive grid generation



S = S(7) (4.3.1)

After solving for the control function, eq. (2.3.10) becomes

Qi = (4.3.2)

The control function then can be determined from the prescribed

distribution s by using central difference approximations to evaluate

S,. As noted previously, the resulting distribution may vary somewhat due

to the coupling of the two dimensional adaptive grid equations and

effects of orthogonality. This approach, which is used in some of the

cases solved here to specify the spacing in the normal coordinate

direction, is similar to the technique of Thompson et al. [15] to obtain

boundary point distribution in the domain interior. The equation for the

grid point spacing used to obtain the arclength distribution is

Asj = As*(l+e)n (4.3.3)

where e is determined so that ST matches the total length of the

coordinate line (see reference 16), As* is a specified spacing at the

projectile surface, and n is equal to the index j.

The actual choices for the functions h and g in eqs. (4.1.4) are

discussed in Chapter VIII.


V.1 The Governing Equations

The governing equations for viscous transonic flow are the Navier-

Stokes equations and an equation of state. The Navier Stokes equations,

which represent the conservation of mass, momentum and energy form a set

of parabolic/hyperbolic partial differential equations when the time

dependent terms are retained. To solve these equations numerically on an

arbitrary domain the governing equations are scaled and then transformed

onto the curvilinear coordinates. For problems in three dimensions, three

curvilinear coordinates are required, denoted here as , q and .

However, for bodies of revolution, such as the axisymmetric projectiles

at zero angle of attack considered here, the solution is assumed

invariant in the circumferential direction, designated as the

direction, and the equations can be reduced since only two independent

spatial variables are required. A further reduction in the equations

ccmes from the thin layer approximation. In the thin layer approximation

all the viscous derivatives in the direction along a solid surface (the

direction in the present applications) are dropped from the equations.

Due to restrictions on computer storage and CPU time, it is not usually

feasible to resolve the viscous terms in both the directions normal to

the surface and along the surface. For high Reynolds number flows, such


as the transonic flow considered here, large velocity gradients are

expected to occur normal to the surface and most of the grid points are

used to resolve these gradients. As larger spacing is then used along the
surface and the viscous derivatives in that direction are expected to be

small in comparison, they are dropped from the equations.
The conservative form of the transformed thin layer governing
equations for axisynmetric flow, written in the compact vector form are

atq + a8 + a ? + t =Rela



q p
q=J* ( pv
| PW

A pUU+ xp
, E=J pvU+ yp
P pwU+ zp
S(e+p) U

] A
S, F=J

p uV+n xP
( pvV+yP
pwV+7I zp

6=J* ( -pW2y p/y
[ pW(yU+yV)
I 0


uRu, + (u/3)T'7x
uRv 7 + (u/3)T'y
uRw, + (u/3)TVz
A( +(u2+v +w2) +uPr1(-1)-i(a2)2
+(u/3) ( 7xu+, 2yv+t2xw)T


R = 'x2 + ny2 + nz2

T = ?,u + n yv, + nxW,

S = J


Here, p is the density, p is the pressure and E is the total specific

energy. The velocities u,v and w correspond to the x, y and z directions,

respectively, as indicated in Figure 5.1 and U, V, and W are the

contravariant velocity components corresponding to the , t7 and

directions respectively. The coefficient of viscosity is u, Pr is the

Prandtl number, Re is the Reynolds number, -y is the ratio of specific

heats and a is the local speed of sound. The details of the

transformation and application of the thin layer approximation to obtain

eqs. (5.1.1) are available in reference [13]. The assumption of a

calorically and thermally perfect gas was made in closing the system of

equations. The coefficient of viscosity pa is considered variable and it

is also necessary to account for the effects of turbulence for the

present applications. Under the assumption of invariance, only a two

dimensional grid network containing and 7 coordinates in a plane

through the axis of symmetry will be required. It should be noted,

however, that the metric coefficients represent the three dimensional

coordinate transformation and are therefore different from the

expressions of eqs. (2.1.4) and the Jacobian J* now represents the grid

cell volume.

V.2 Solution AlQorithms for the Governing Equations

The implicit spatially factored numerical scheme of Beam and Warming

[17] is used to solve the transformed governing equations. The scheme is

implemented for three dimensional problems in a code due to Steger and

Pulliam [18]; the code used here is a modified version which was

developed by Nietubicz et al. [19] for the efficient calculation of

axisymmetric projectile problems and is referred to here as the Thin

Figure 5.1 Coordinate systems for projectile configuration


Figure 5.2 Boundary configuration for projectile
with sting


Layer code. The code employs Sutherland's law [20] to relate the

viscosity and heat conductivity to the temperature and the effects of

turbulence are modeled using the eddy viscosity model of Baldwin and

Lomax [21].

The numerical algorithm is a tine-marching scheme in which steady

state solutions are obtained as the time asymptotic solution of the

unsteady equations. Spatial derivatives are approximated with central

differences and a forward difference approximation is used for the

unsteady terms. The time marching scheme of the finite difference

approximation to the governing equations are written in delta form as

(I + h 6. eID1) (I + h6n hRelJkJ1 ejD2A

= h(dC + 8, Re1 6nA 6) ED3

where A and B are the Jacobian matrices

a A A a
B (5.2.2)
aq a

and M is obtained as the time linearization of S. Here, 6 is the central

difference operator. Details of this scheme can be found in Warming and

Beam's [17] work, as well as in that of others. The terms D1, and D2 are

implicit dissipation operators and D3 is an explicit dissipation term

D1 = J-1 A v7 J

D2 = J-1 An VI J (5.2.3)

D3 = J-1 ((As V7)2 + A, V7)2) j


where D and v are forward and backward difference operators,

respectively. Since the scheme is not naturally dissipative these terms

are added to the finite difference scheme to prevent nonlinear

instabilities and high frequency spatial oscillations. The constant

coefficients fI and EE used here are 4At and 2At respectively which are

consistent with values used in other applications of this code [19]. In

general, these coefficients should be chosen large enough to dampen the

instabilities but kept small enough to prevent any degradation of

solution accuracy.

In Chapter VII, same examples using the Thin Layer code with the

Adaptive Grid code show that there are problems in obtaining a converged

solution for the thin layer equations when the grid is adapted. Therefore

another numerical scheme for solving the thin layer equations is also

used. This scheme is the TVD scheme developed by Yee [22] for the

multidimensional Navier-Stokes equations. The theory and development of

TVD schemes is beyond the scope of current discussions. However, for

purposes here it is sufficient to view the TVD scheme as an improvement

in the artificial damping terms in the Beam and Warming scheme. In fact,

the TVD scheme developed by Yee for multidimensions can be implemented

into the existing Thin Layer code based on the Beam and Warming scheme by

modifying only the added dissipation terms. The original Thin-Layer code

has been modified to provide an option for the TVD scheme [23] and is

referred to here as the TVD code. These 'sophisticated' dissipation terms

switch the scheme from second order to first order accuracy at points of

extrema and result in a more dissipative scheme in those regions. It

should be pointed out that the primary purpose of the current research is

to develop an adaptive grid generation technique. The TVD scheme is used


here because of problems occurring when using the Beam and Warming scheme

in conjunction with the adaptive grid generation technique. These

problems are discussed more fully in Chapter VII.

The first flow problem considered is a the axisynmetric projectile

flow problem with an attached sting to eliminate the base flow region. In

both the Thin Layer code and the TVD code the boundary conditions are

implemented explicitly by updating the flow variables at the boundary

after each time step. Referring to Figure 5.2 which shows the basic C-

type grid configuration for the axisymmetric projectile flow problem,

there are four different sets of conditions applied at the four bounding

coordinates. Along the outer boundary, line AB, free stream conditions

are assumed. Thus the outer boundary is required to be sufficiently far

from the projectile so that this assumption is valid. Along the upstream

line of symmetry, line AD, the variables are obtained by requiring

symmetry about the x axis. Numerically this is done by setting a second

order forward difference approximation to the first derivative along the

coordinate for each variable equal to zero and solving for the values

at the points on the axis in terms of known values in the domain. Along

the dowmstream boundary, line BC, the variables are obtained by

extrapolation along the coordinate lines. However, if the free stream

velocity is subsonic the pressure is fixed. This is done to prevent

oscillations in the pressure distribution that would otherwise occur

[24]. Along the projectile surface, line CD, the no slip condition holds

for viscous flow, thus the velocity components are zero. The density is

obtained by zero order extrapolation along the 7 coordinate lines and the

pressure is obtained by satisfying the momentum equation along the



Another problem considered is the projectile base flow problem. In

Figure 5.3, which shows the basic 0-type grid configuration, the

coordinate lines bend around the sharp corner at the base and end on a

downstream axis of symmetry which is now an q coordinate. The existing

codes are used for this problem by simply modifying some of the boundary

conditions. The boundary conditions along the upstream axis of symmetry,

line AD, and the projectile surface, line CD, are the same as those used

in the projectile problem with sting. For the downstream axis of

synnetry, line BC, the same conditions along the upstream axis are used;

however, a backward finite difference formula is required. The freestream

conditions are applied along the outer boundary along the line AA' and

the downstream boundary conditions for the projectile with sting problem

are applied along the portion of the outer boundary between A'B, except

that the extrapolation is done along the r coordinate lines.

V.3 Examples on Fixed Grid Networks

In order to show the general characteristics of the Beam and Warming

and TVD schemes, and for purposes of comparison with the solutions

obtained with the adaptive grid generation technique, the two schemes

were used to solve the projectile flow problem with sting on a fixed grid

network. The fixed grid network was obtained using a method developed by

Steger et al. [16] which is an effective technique for external flow

problems. Two equations, one enforcing orthogonality and one specifying

the grid cell volume, form a set of hyperbolic partial differential

equations which are solved numerically to generate the grid network. In

the marching solution algorithm, the initial grid points are specified

along the projectile surface. The solution is marched outward in the 7



Figure 5.3 Boundary configuration for projectile base
flow problem


3.017 -- ----- 2.0


Figure 5.4 SOCBT projectile configuration


coordinate direction, defining each successive coordinate line such

that the spacing between each coordinate line is consistent with the

prescribed grid cell volume. At the bounding r coordinate lines, the

upstream axis of symmetry and the downstream boundary of Figure 5.2, the

coordinate lines are required to intersect the t coordinates at right


The projectile used in all the calculations is a 6-caliber, secant-

ogive cylinder boattail configuration which is shown in Figure 5.4.

There are experimentally determined pressure coefficients available at a

series of points along the projectile surface for a range of Mach numbers

in the transonic range [25]. The data will serve as a means of comparing

the computed solutions. A grid network for this projectile problem with

sting was generated using this technique with 70 points in the

coordinate direction (IM=70) and 60 points in the 7 coordinate direction

(JM=60). As shown in Figure 5.5, most of the points in the t? coordinate

direction were clustered near the projectile surface to resolve the

boundary layer. The spacing at the projectile surface in the i coordinate

direction is 0.00002. The outer boundary is extended to four times the

projectile length. In the streamwise direction, the points were clustered

near the secant-ogive/cylinder and cylinder/boattail junctures in

anticipation of the expansion waves that will occur in transonic flow.

The first solution on this grid network is obtained with the Thin

Layer code for Mach 0.96. The Reynolds number for this and all examples

is 76,000. For all the results obtained, the time step used for a

converged solution is 0.1. Figures 5.6 and 5.7 show the convergence

process of the solution algorithm for the indicated time spans. It is

quite obvious that the solution oscillates in the region containing the

Figure 5.5 Initial fixed grid network for projectile
with sting (70 x 60)

MACH 0.96



-0.2 -


-0.6 -

Figure 5.6



0.0 -

o experiment
- computed

Converging solution for the Beam and Warming

MACH 0.96

;,, *..,,4

-0.2 /

-0.4 /

g i o experiment
-0.6 1 7 computed

Figure 5.7 Converged solution for the Beam and Warming

MACH 0.96

o experiment
- computed

Figure 5.8

R 1


Converged solution for the TVD scheme

Beam and Warming

1 2 3 4 5 N (1000)

Least squared residue




Figure 5.9


shocks. However, after approximately 5000 time steps the solution does

converge and agrees well with the experimental data as indicated by the

comparison of the experimental and computed pressure coefficient along

the projectile surface. Figure 5.8 shows the results obtained using the

TVD code. The solution also agrees well with the experimental data;

moreover it converged within 1600 time steps. This difference in the

convergence ofthe two schemes is shown in Figure 5.9 which plots the

least square residue R, versus the number of time steps N. Here the

slower, oscillatory convergence of the Beam and Warming scheme is evident

while the TVD scheme appears to converge monotonically.

The projectile base flow problem is also solved on a fixed grid using

the TVD code. The grid network for this calculation was obtained using a

modified version of the technique of Steger et al. [16] in which the

boundary condition on the downstream bounding j coordinate line has been

changed so that the coordinate lines intersect with the x axis at right

angles. The grid network obtained with IM=70 and JM=60 is shown in Figure

5.10. As can be seen, the coordinate lines now bend around the sharp

corner, and end at the downstream axis of symmetry forming an 0-type grid

network. The computed pressure coefficient obtained for this problem

using the TVD code is shown in Figure 5.11. Again, the solution

converged in approximately 1600 time steps. It compares well, as before,

over the secant-ogive and cylinder portions but then differs from the

experimental results at the end of the boattail. No experimental data is

available for the base region. There are a number of reasons that may

cause the difference between the computed and experimental results near

the corner. One reason may be the presence of a narrow sting used to

support the experimental model that is not present in the computational


configuration. The difference could also be due to poor resolution in the

projectile corner region or an inadequate turbulence model near the flow

separation region.

Figure 5.10 Initial fixed grid network for projectile
base flow problem (70 x 60)

MACH 0.96


S Or0) U i ed

MACH 0.96

- b




Figure 5.11 Converged solution for projectile base flow
problem using the TVD scheme











VI. 1 The General Procedure

Prior to coupling the Adaptive Grid code with the flow solvers to

solve the projectile flow problem, it was found beneficial to modify the

adaptive grid generation method. The following sections include an

example showing the control over orthogonality and then describes two

modifications made to the adaptive grid scheme to improve the general

technique. It should be noted that each of the modifications is made to

improve the method for the present flow problem and geometry and may not

be necessary or useful in other problems.

For each example the initial grid network used was that of Figure

5.5, or that of 5.10 if the base flow configuration is used in the

example. In the adaptive grid code, the control function P was set equal

to a constant so that the clustered points along each C coordinate line

in the initial grid would tend to spread and equal spacing would result.

The control function Q was defined using eqs. (4.3.2) with the

distribution of eq. (4.3.3) so that the same point distribution along the

r7 coordinates in the initial grid would be retained in the grid obtained

using the Adaptive Grid code. Each of the figures used in the examples

show only the points in the reduced grid network.


VI.2 Orthocronality

In order to demonstrate the effect of enforcing orthogonality, two

cases were considered, one with the parameter A = 0.0, the second with

= 0.5. For the present projectile, the most difficult area to obtain an

orthogonal grid is over the projectile nose. This is due to the sharp

nose configuration which results in the two boundaries, the projectile

surface and the upstream axis of symmetry, meeting at an obtuse angle.

The two sets of coordinates, which run parallel to the boundaries, then

tend to intersect at similar angles. This result is shown in Figure 6. la,

which shows an enlarged view of the grid in the nose section that

resulted when A = 0.0. To improve the orthogonality in this region,

another case for A = 0.5 was run, and as the results of Figure 6. 1b

indicate, the orthogonality of the grid coordinates near the nose was in

general increased. For all the solutions to the transonic flow problem

using the adaptive grid technique reported in Chapter VIII, A was set

equal to 0.5.

VI.3 Curvature

One inherent result of grids obtained as solutions to equations based

on the Laplacian operators is the effect of curved boundaries on the

spacing of interior coordinate lines. In general, boundaries that are

convex will attract coordinate lines and those that are concave will

repel coordinate lines. In the context of the adaptive grid generation

scheme presented here, this effect will alter the grid spacing dictated

by the control functions, either increasing or decreasing the spacing

that would result in the absence of curvature effects. The severity of

the boundary curvature's influence on the grid is best seen in solving

Figure 6.la Grid network using X=0.0

Figure 6.1b Grid network using X=0.5


the adaptive grid generation equations for the base flow grid network.

Recall that thel control function Q is defined so that the original grid

point spacing in the normal direction would be retained. However, the

point spacing resulting as a solution to the adaptive grid equations

differs along each 7 coordinate line as shown in Figure 6.2a. The

attraction of points is most prominent near the sharp corner at the

projectile base where the effect is large due to the high curvature at

the corner.

A simple and effective technique to eliminate the effect of curvature

can be developed by considering the solution of Laplace's equations

between two concentric circles C1 and C2 with radii R1 and R2


V2 = 0
V2 = 0

Solving Laplace's equation is equivalent to solving the adaptive grid

equations with constant control functions and neglecting orthogonality

(i.e. A = 0.0). In using constant control functions, we expect the

coordinates to be equally spaced in each direction. It is advantageous to

write the equations in polar coordinates (r, ), and it is sufficient here

to seek a solution in which q is a function of 0 only and 7 is a function

of r only. The coordinate lines then align themselves with the 0

coordinate direction and the ri coordinate lines become aligned with the

radial coordinate r. The problem can then be formulated as

Figure 6.2a Grid network obtained without corrections
for curvature

Figure 6.2b Grid network obtained with corrections
for curvature


d2 i 1 do
+ - =0
dr2 r dr


= 0 on C1 and C2

S= R1 on C1

= R2 on C2

The boundary conditions corresponds to uniform spacing along the

coordinate lines coincident with the bounding concentric circles. The

spacing along each coordinate direction is the inverse of the first

derivatives O and Tr which by integrating eqs. (6.2.2) are found to be

8 = a
b (6.2.3)
Tr r

where a and b are constants of integration. These results indicate that

the spacing in the direction is uniform while the spacing in the ?I

direction is smaller near the inner circle C1 and larger near the outer

circle C2 The term in the governing equation responsible for this

result, designated here as K, is

1 dT7
K = - (6.2.4)
r dr

which can be viewed as the curvature of the intersecting coordinate

line, (1/r), multiplied by the inverse of the spacing along the tI


coordinate line. The effects of the curved coordinates on the spacing

along the r coordinate lines can be eliminated for this problem simply by

subtracting K frame the governing equation for Y7,

V2 = 0


S1 d?
V2 =0
r dr

This technique can be extended to the general adaptive grid equations

by developing the curvature term K for each curvilinear coordinate line

in the computational domain and then subtracting them fram eqs. (2.2.7).

The term needed to cancel the effects of a curved rP coordinate line on

the spacing along a coordinate line is comprised of the curvature of

the i) coordinate line and the spacing along the C coordinate line. The

curvature p of a n coordinate line is defined as

p = I (6.2.6)

where r is a unit tangent to the rI coordinate and is in the computational


-_ ) (x,,y)
--Z -(6.2.7)

where a is defined in eq. (2.3.5). By using the two relationships

S = / d_ i d (6.2.8)
ds s7 dt7

the curvature p can be shown to be


=( yx -- xLI y-17t (6.2.9)

By multiplying this result with the inverse of the spacing along a

coordinate (1/7-y), the following expression is obtained

Sy~x x~y,,
K1 =f? a? 17 (6.2.10)
j7 ct3/ 2

which is the term that is subtracted from the adaptive grid equation for

in order to cancel the effects of curved t7 coordinate lines. An

expression for the influence of curved coordinates on the spacing along

P7 coordinate lines can be obtained in a similar manner and is

y~x^ x~y^
K2 = (6.2.11)
j/a 'Y- 3/2

Thompson et al. [15] have derived similar expression, but do not

implement them in the same way. It should be noted that in the

computational plane, the terms K1 and K2 contain second derivatives of

the variables x and y, and when eqs. (2.2.7) are transformed onto the

computational plane, the second order central difference approximations

to these terms depend explicitly on the point (xi j, Yi,j). However, in

implementing the Newton-Raphson iterative solution procedure, the terms

K1 and K2 are considered part of the source terms. Thus, the right hand

side of eqs. (2.3.7) becomes

S1 K1J3
-'- P Klj

aJ Q r,
S2 K2J3

This is equivalent to neglecting the dependence of the finite difference

approximation on the point (xi,j, Yi,j) in obtaining eqs. (3.1.5) for the

Newton-Raphson iterative procedure. The J3 term in eqs. (6.2.12) appears

in the transformation of eqs. (2.2.7) into the ccaputational plane.

The adaptive grid generation equations, modified in this manner, will

produce the grid shown in Figure 6.2b. In contrast to Figure 6.2a, it is

clear that the effects of the curved boundaries have been successfully


VI. 4 Internal Grid Boundaries

One other problem that arose in the preliminary testing of the

adaptive grid generation technique is the rapid upstream bending of the ri

coordinate lines emanating from the secant portion of the projectile

surface. Figure 6.3a shows a grid network obtained using the adaptive

grid generation equations with the curvature terms included and with x =

0.5. As can be seen, the Pi coordinate lines emanating frcm the secant

portion of the projectile quickly bend forward to fill the upstream

region of the domain which results in a very skew grid. This

characteristic of the grid does not appear as a problem until a solution

using the Beam and Warming scheme is sought on such a grid. Figure 6.3b

shows a time history of the solution obtained on such a grid for Mach

0.96 and as can be seen the pressure coefficient oscillates in time over

Figure 6.3a Grid network upstream of the projectile

MACH 0.91






o experiment
C- computed

Figure 6.3b

Lack of convergence over secant ogive with
Beam and Warming scheme


the secant portion of the projectile. The solution did not converge and

instead the oscillations

continued and eventually effected the solution downstream.

In order to prevent the upstream bending of the ri coordinates, one

coordinate line in the initial grid network is designated as an internal

boundary. This rt coordinate line remains fixed in space during the

iterativesolution procedure, thus no other rl coordinate lines will pass

through it. However, the points defining the internal boundary coordinate

line are allowed to move along the coordinate so that they remain

consistent with the spacing of the points along the adjacent r, coordinate

lines. This is accomplished by updating the arclength distribution along

the internal boundary after each iterative sweep in which the new

distribution is defined as the average of the arclength distributions

along the two adjacent rY coordinate lines. Figure 6.3c shows the grid

obtained when the 12th q coordinate line from the projectile nose is

designated as an internal boundary and as can be seen it prevents the

coordinate lines from bending upstream. It has been chosen to lie in a

region in which no major grid clustering is expected so that it will not

effect the adaptation of the grid network.

The length and number of points used along the C coordinate lines in

each section may result in an abrupt change in the point spacing along

the coordinate lines passing through the boundary as can be seen in

Figure 6.3c. Therefore the control function P, controlling the spacing

along the coordinate lines is locally modified near the internal

boundary so that the grid point spacing varied smoothly through the

internal boundary. This can be accomplished by requiring the normalized

rate of change of grid point spacing along a coordinate line (sd/sC)

Figure 6.3c Grid network with internal boundaries

Figure 6.4 Grid points near internal boundary

be identical at the two points adjacent to the internal boundary.
Referring to Figure 6.4, if the h coordinate line designated as an
internal boundary contains the ith point along a coordinate line, we

then require

( --] =s-- ] (6.3.1)
s I i-i= s" I i+1

To implement this requirement into the solution procedure, consider the
one dimensional equation for adaptation along a C coordinate, eq. (2.3.8)
which is written here as

P s
PC __
P s

The left hand side of eq. (6.3.2) appears in the source terms of the two
dimensional equations for adaptation, eqs. (2.3.7). Substitution of eq.
(6.3.2) into the source terms for a grid point just ahead of the internal
boundary results in the modified source terms at the point i-1

f s 1 633
(Si)i-l.'j = ai-ljJi-lj I (6.3.3)

The modified source terms for the point just behind the internal boundary
are obtained by replacing i-1 with i+1 in eq. (6.3.3). In order to apply
the iterative solution technique a value for the terms (sc/sc) in eqs.
(6.3.3) for each point adjacent to the internal boundary will be needed.
Their final values are not known a priori since they will depend on the
final position of the points in each section, therefore they are set to

the current value

[s 1_ (si-2 -2si + si+2) (6.3.4)
s Ji-l,j (si+2-si-l)/2

The grid that is obtained when using this technique is shown in Figure

6.5 in which it can be seen that the abrupt change in spacing across the

internal boundary of the grid in Figure 5.3c has been reduced.

The use of internal boundaries provides a simple and effective way to

control the grid points in the domain. In fact, it was found helpful to

also designate the Y coordinate line which separates the region over the

sting from that over the projectile as an internal boundary to keep the 7

coordinate lines that originally began on the projectile surface from

moving over the sting.

Figure 6.5 Smoothed grid spacing near the internal


The adaptive grid generation technique can be incorporated naturally

into the solution process of the flow solvers since the flow solver

algorithms are time-marching. After a specified number of time steps, the

grid network can be adapted to the control functions based on the current

values of the flow variables. Thus the time-marching algorithm of the two

flow solvers used here, the Thin Layer and TVD codes, proceeds as usual

with intermittent calls to the Adaptive Grid code to update the grid

network. In this approach, however, the relocation of the grid points in

the physical domain due to the grid adaptation must be considered. The

values of the flow variables at any time step are defined at the current

grid point locations and any movement of the points will therefore

distort the solution. This effect is critical for unsteady flow problems

since the grid point movement will effect the solution accuracy. For

steady flow problems, the grid point movement may effect the convergence

rate and if it is large enough, it may cause the solution to diverge.

There are two basic approaches to account for the grid point movement

resulting when adapting the grid network. The first is to view the grid

network as a moving grid and include the time metrics in the coordinate

transformation of the Navier-Stokes equations. These metrics can be

written in the computational plane as


Ct = CxXr + y Yr
I?t = 7xXr + 'yYr

The quantities x, and y, can be approximated by using a backward

difference expression

x- = (Xn+l-xn)/At
Yr. (yn+l-yn)/At (7.2)

two in which xn and yn are the grid point locations of the previous grid

network and xn+1 are yn+l are the values of the adapted grid network.

This approach requires the grid to be adapted every time step which is

inefficient for steady flow problems. Furthermore, it has been shown that

special procedures must be used in evaluating the time derivative of the

Jacobian [20] in order to maintain the conservative property of the

finite difference schemes.

Another approach, used here, is to interpolate the flow variables

from the previous grid network onto the adapted grid network. This

approach can be applied to both steady and unsteady flow problems and

does not require the grid to be adapted every time step. Thus this

approach can be used for both steady and unsteady flow problems. For the

two- dimensional grid networks considered here, a linear interpolation is

currently used. Each grid cell in the previous grid network is divided

into two triangles. Once a grid point in the adapted grid network is

located within a triangle of the previous grid network, the value of the

flow variables at that point can be determined by interpolation.

Referring to Figure 7.1, the value of a flow variable q at the point in

the adapted grid network can be determined from the values of the flow


. A,


0 Previous grid network

* current grid network

Figure 7.1 Flow variable interpolation


variables at the vertices of the triangle as

(q1A1 + q2A2 + q3A3)
q = (7.3)
(A1 + A2 + A3)

where A,, A2 and A3 are the areas of the triangles formed with the grid

point in the adapted grid network.

The basic algorithm for incorporating the adaptive grid generation

technique into both the Thin Layer code and the TVD code consists of the

following steps:

1. Input an initial grid network into the flow solver code and
begin the time-marching algorithm.

2. Advance the solution a specified number of time steps, and
submit the current grid network to the Adaptive Grid code.

3. Obtain an adapted grid network using the Adaptive Grid code.

4. Interpolate the flow variables onto the adapted grid network
which now becomes the current grid network.

5. Repeat steps 2 through 5 until the solution converges for steady
flow problems or repeat the steps until a specified time is
reached for unsteady flow problems.

This process which couples the flow solver and the adaptive grid

generation technique is referred to here as the self-adaptive

computational technique. For the two problems considered here, the

axisymmnetric transonic flow over a projectile with sting and the

transonic projectile base flow problem, only steady flow solutions are

sought. Since the adapted grid network depends on the flow variables via

the control functions, the convergence of the solution implies also that

the grid network ceases to change with time. In all the examples of the

following chapters the grid networks of Figures 5.4 and 5.10 are used as


the initial grid.

For the purpose of adapting the grid network, the source terms

corrected to eliminate curvature effects, eqs. (6.2.12) were used in the

adaptive grid generation equations. In solving the projectile problem

with sting two internal t7 coordinates were designated as internal

boundaries. The first is over the secant ogive as shown in Figure 6.5,

the second separates the grid network over the projectile from that over

the sting and is located at approximately x=7.0. In solving the

projectile base flow problem, no internal boundaries were used.


VIII. 1 Choices for the Control Functions

In applying the self-adaptive computational technique it is necessary

to define the control functions for adaptation. For the transonic

projectile problem considered here both shocks and expansion waves will

occur approximately normal to the projectile surface and require refined

grid point spacing in the coordinate direction for adequate resolution.

The original choice for the function f in the control function P to

obtain this clustering was the pressure gradient. However, it was found

in the numerical experiments that the expansion waves occurring at the

secant ogive/cylirnder and cylinder/boattail junctures required smaller

spacing than the shocks for good resolution. The pressure gradient was

largest in the shocks and thus smaller spacing occurred in the shock

regions rather at the expansion waves. A better choice for the control

function was found to be the curvature of the pressure distribution along

the streanwise coordinate direction,

f = ----- I 8S I(8.1.1)
( 1 + (an )2)3/2


where the derivatives are with respect to the arclength s along each

coordinate line. This function is largest at the expansion waves where

the pressure gradient changed sign. At the base and peak of a shock eq.

(8.1.1) also increased but it is not as large as the value near the

expansion waves. Thus the smallest spacing will occur near the expansions

and yet the grid points will also cluster, to a lesser extent, in the

vicinity of shocks.

In the normal direction, the important area requiring increased

resolution is the boundary layer and the logical choice for the function

g in the control function Q is the velocity gradient. However, it was

found that the control function based on the velocity gradient resulted

in too many points being placed in the boundary layer and led thus to a

rapid stretching of the grid point spacing along the Y7 coordinate lines.

Another choice for the function g, the exponential of the velocity


g = exp as (8.1.2)

where u is the velocity magnitude and s measures along the r, coordinate

line, was found to yield better results. Figure 8.1 shows a plot of the

control function Q along a typical S,, coordinate line. The control

function resulting from the exponential of the velocity gradient is

similar to the control function corresponding to the grid point

distribution in the fixed grid networks which is known to give good

results provided enough points are used. The control function obtained

using the velocity gradient however is large for many grid points and,

o velocity gradient

* prescribed distribution

o exponential of the
velocity gradient

20 40

77 60

Figure 8.1

Comparison of control functions

MACH 0.96

o experiment
- computed

Figure 8.2 Lack of convergence using a self-adaptive
computational technique with Beam and Warming

Q (1o)









thus, results in too many points being clustered into the boundary layer

region at the expense of adequate resolution in other areas. The minimum

and maximum grid point spacing along the t coordinate lines is 0.00002

and 6.0 respectively. In this direction it is important to have the small

spacing at each point along the projectile surface in order to adequately

resolve the viscous sublayer region. The parameter a of eq. (4.1.4) was

determined for each j coordinate line so that the small spacing would

occur for each point on the projectile surface. However, in the

coordinate direction, the shocks and expansion waves are not expected to

reach the outer boundary and the small spacing needed along the

projectile surface is not required near that boundary. With the minimum

and maximum grid point spacing prescribed as 0.02 and 0.22, respectively,

a is determined only for the coordinate line coincident with the

projectile surface and that value is used for all the remaining

coordinate lines. This had the effect of increasing the minimum spacing

as the magnitude of the shocks and expansion waves diminished near the

outer boundary.

VIII. 2 Results Using the Self-Adaptive Computational
Technique: The Beam and Warming Scheme

The first application of the self-adaptive coputational technique is

with the Beam and Warming scheme for the calculation of the projectile

flow with sting. The Mach number is 0.96 and the Reynolds number is

76,000, a case for which experimental data is available. The time step

used in the calculation is 0.1 and the explicit and implicit dissipation

terms are 2 and 4 times the time step respectively. The initial grid

network used is that of Figure 5.5 with IM=70 and JM=60. In the first


attempts to use the adaptive grid technique, the grid network was adapted

every 200 time steps. A series of calculated pressure coefficients

obtained with this approach are shown in Figure 8.2. As can be seen, the

pressure coefficient oscillates widely over the projectile surface and

shows no indication of converging. In a series of attempts to obtain

better results, the time step was decreased to 0.05, the dissipation

terms were doubled and the number of time steps between adaptation were

both increased to 500 and reduced to 50. However, similar results were

obtained in each case.

In the next case tried, the number of points in the streamwise

direction was increased to 90 points, thus IM=90, and the grid network

was adapted every time step. The results for this case are nuch better,

however a converged solution was not obtained. Figure 8.3a shows the

computed pressure coefficient at the times indicated, and as can be seen,

the solution appears to have converged everywhere except in the shock

boundary layer interaction regions. The shocks propagate upstream and

downstream in a cyclic fashion. Figure 8.3b shows the same phenomenon

occurring at a later time. In fact, the solution was continued for 8000

time steps and the propagation of the two shocks showed no indication of


Figure 8.4a shows the reduced grid network which resulted at 5000

time steps. As is evident, the grid point spacing in the q coordinate

direction is clustered near the projectile surface with a distribution

similar to that of the hyperbolic grid in Figure 5.4. The spacing varies

somewhat in the shock and expansion wave regions, most likely as a result

of enforcing orthogonality. The clustering in the streamwise direction

MACH 0.96

5a a o
Us wes-

.^ ite,

o experiment

Figure 8.3a Computed solution between 4000 and 5000 time
steps with Beam and Warming scheme


MACH 0.96

0.2 h



o experiment
-0.6 computed

Figure 8.3b Computed solution between 5000 and 6000 time
steps with Beam and Warming scheme







Figure 8.4a Adapted grid network at 5000 time steps
(90 x 60)

Figure 8.4b Adapted grid network at 6000 time steps
(90 x 60)


clearly indicates a response to the control function P, which is based on

the curvature of the pressure function. Clustering occurs at both

junctures on the projectile surface and the adaptation to the shock

structures on the cylinder and boattail are also quite evident. Note that

the smallest spacing appears to be located in the shocks well above the

projectile surface. This occurs because there is more competition for the

small grid spacing on the projectile surface since both the shocks and

expansion waves are strong in that region. Figure 8.4b shows the grid

network corresponding to a later time. In this Figure and those

following, the coordinate lines are not shown in order to gain a better

view of the grid spacing near the projectile surface. The grid

configuration is basically the same as that of Figure 8.4a, however, the

position of the shocks as indicated by the adaptation has changed. The

shock on the cylinder has propagated downstream and the boattail shock

has moved upstream. Referring to the shock positions evident in Figures

8.3a and 8.3b it can be seen the clustering in the grid network is

successfully following the shock patterns.

In an attempt to improve these results the time step was reduced to

0.05 and the dissipation terms were again increased, however, similar

results were obtained. Also the convergence criterion 6 of eq. (3.1.12)

used in the iterative solution procedure for the grid generation equation

was varied. At the value of 6=0.005, only one iteration was required for

the grid network to converge. The rapid convergence occurred because the

grid was updated every time step and the solution, and consequently the

control functions, did not change much between grid adaptations. When the

criterion was reduced to 6=0.0025, 1 to 3 iterations were required for


the grid to converge at each time step, however it had no effect on the

cyclic propagation of the shocks. Although a converged solution was not

obtained using the adaptive grid generation technique with the Beam and

Warming scheme, there are still some encouraging results with respect to

the adaptive grid generation procedure. The maximum and minimum grid

point spacing occurring along the projectile surface in the adapted grid

network were within 5 percent of their prescribed values. The average

minimm spacing along all the rj coordinate lines was within 3 percent of

the prescribed value of 0.00002 with a maximum deviation of 10 percent.

The adapted grid networks shown in Figures 8.4a and 8.4b show that the

adapted grid generation equations are capable of responding reliably and

predictably to the chosen control functions and yield the prescribed

maximum and minimum grid point spacing.

The adaptive grid generation equations, thus, have been shown to

provide a well adapted grid network for the transonic flow problems

provided that a good choice is made for the control functions. However,

the convergence of the solution was effected by the use of the adaptive

grid generation procedure with the Beam and Warming scheme. One reason

for the failure of the solution to converge may be due to interpolation.

Error introduced due to interpolating the flow variables onto each

successive adapted grid network may continuously perturb the solution,

causing it to oscillate. Another cause of this phenomenon may be the form

of the dissipation terms in the Beam and Warming scheme. The added

dissipation terms contain derivatives with respect to the ccuputational

coordinates and were added to the scheme after the governing equations

were transformed. Since these terms are a function of the transformed


coordinates any change in the grid network used will change the

dissipation terms and consequently alter slightly the transformed

governing equations. Thus the continuous adaptation of the grid network

continuously alters the dissipation terms and may be the reason the

solution does not converge.

VIII. 3 Results Using the Self-Adaptive Cocmputational
Technique: The TVD Scheme

Since the solutions obtained using the Beam and Wanrming scheme did

not converge, another scheme, the TVD scheme was also used. Recall that

the primary difference between the two schemes is the dissipation terms

added to the governing equations. The initial grid network used is that

of Figure 5.4 with IM=70 and JM=60. The time step used is 0.1 and the

grid was adapted every 200 time steps. For Mach 0.96, the solution

converged easily within 2000 time steps. Only one iteration was needed

for the grid network to converge at 1800 and 2000 time steps which also

indicates a converged solution. Figure 8.5a shows the computed pressure

coefficient using the adapted grid generation procedure compared to the

solution obtained on the fixed grid network of Figure 5.4. The solutions

agree well with experimental data over most of the projectile surface,

except for the predicted position of the shock over the cylinder. The

position predicted using the adaptive grid generation technique is

centered on the cylinder whereas the position predicted using the fixed

grid network is slightly upstream of the center. Figure 8.5b shows a

shadowgraph for corresponding flow conditions in which it can be seen

that the experimentally predicted shock location is centered over the

cylinder portion of the projectile. Thus, the use of the adaptive grid

generation technique improved the solution accuracy. Figure 8.5c shows


the adapted grid network corresponding to the converged solution. Again,

clustering near the secant ogive/cylinder and cylinder/boattail junctures

is evident as well as the adaptation to the shocks. In comparison to the

fixed grid network of Figure 5.4 it is evident that the grid spacing

along the projectile surface is refined at the center of the cylinder

portion due to the adaptation and is responsible for the increased

solution accuracy.

Two other cases were run using the identical procedure used for the

case of Mach 0.96. Figure 8.6a shows a comparison of the solution

obtained using both the fixed grid network and the adaptive grid

generation technique for Mach 0.91, another case for which experimental

data is available. The only improvement obtained using the adaptive grid

generation procedure is an increase in the peak of the shock occurring

over the cylinder portion of the projectile. However, the peak still does

not reach the experimentally predicted value. In an attempt to increase

the accuracy in that region, 10 more points were added in the

coordinate direction, thus IM=80, however no further improvement in the

solution obtained with the adaptive grid generation technique was

evident. The results thus may be the best that the TVD scheme can yield

given the current grid configuration and the shortened pressure peak is

probably due to the larger dissipative error associated with the TVD

scheme. Figure 8.6b shows the grid network corresponding to the converged

solution. Again grid point clustering in response to the curvature of the

pressure is evident, and as can be seen by the adaptation of the grid

network, both shocks occur directly downstream of the expansion waves.

MACH 0.96



I K7


- adopted grid
. fixed grid

Figure 8.5a Comparison of computed solutions on fixed
and adaptive grid networks with TVD scheme



-B k.. ..

Figure 8.5b SOCBT shadowgraph for Mach 0.96
(reprinted from reference 24)




Figure 8.5c Adapted grid network for converged solution
for Mach 0.96

MACH 0.91
0.2[ o


adapted grid
.-- fixed grid

Figure 8.6a Comparison of computed solutions on fixed and
adapted grid networks with TVD scheme