,q 

0 

t 

e 

- d \ ;; [>/.’JDEPARTMENT OF MECHANICAL ENGINEERING AND MECHANICS 

COLLEGE OF ENGINEERING AND TECHNOLOGY 

OLD - 6 L/--;;A /L, DOMINION UNIVERSITY /P 

NORFOLK, VIRGINIA 23529 

136057 

J CV? 

CONSERVATIVE FINITE VOLUME SOLUTIONS OF A LINEAR 

HYPERBOLIC TRANSPORT EQUATION IN TWO AND THREE 

OIMENSIORS USIffi MULTIPLE GRIDS 

Monchai Kathong 

and 

Surendra N. Tiwari, Principal Investigator 

Progress Report 

For the period ended June 30, 1987 

Prepared for the 

National Aeronautics and Space Administration 

Langley Research Center 

Hampton, Vi rgini a 23665 

’ 

Under 

Cooperative Agreement NCC1-68 

Or. Robert E. Smith, Technical Monitor 

ACD-Computer Applications Branch 

(PAS1-CP-1827 15) CCISEtiVA31VE EIIX‘IB VOLUUE N88-2C934 

SCLUIl.CYS OF A IIbI1B EYPESlECLlC IBAISPOBT 

ECUAllCIR 1Y TbO ELC IB€iEE CIEEESJCUS USIHG 

(Cld Ccriaica Uoiv.) 53 p CSCL 12A 63/64 0136057 

BELTIELE GRID: frcgress Befcrt, 3C Jnn, 1987 Unclas 

I 
DEPARTMENT OF MECHANICAL ENGINEERING AND MECHANICS 

COLLEGE OF ENGINEERING AND TECHNOLOGY 

OLD DOMINION UNIVERSITY 

NORFOLK, VIRGINIA 23529 

CONSERVATIVE FINITE VOLUME SOLUTIOWS OF A LINEAR 

HYPERBOLIC TRANSPORT EQUATION IN TU0 AND THREE 

DIMENSIONS USING MULTIPLE GRIDS 

BY 

Monchai Kathong 

and 

Surendra N. Tiwari, Principal Investigator 

Progress Report 

For the period ended June 30, 1987 

Prepared for the 

National Aeronautics and Space Administration 

Langley Research Center 

Hampton, Vi rgini a 23665 

Under 

Cooperative Agreement NCC1-68 

Dr. Robert E. Smith, Technical Monitor 

ACD-Computer Applications Branch 

Submitted by the 

Old Dominion University Research Foundation 

P.O. Box 6369 

Norfolk, Virginia 23508 

July 1987 
FOREWORD 

This is a progress report on the research project "Numerical Solutions 

of Three-Dimensional Navier-Stokes Equations for C1 osed-91 uff Bodies," for 

the period ending June 30, 1987. 

investigate the 'Tonservative Finite Volume Solutions of a Linear Hyperbolic 

Transport Equation in Two and Three Dimensions Using Multiple Grids." This 

work was supported by the NASWLangley Research Center (Computer Appl ications Branch, Analysis and Computation Division) through the cooperative 

agreement NCC1-68. 

Smith of the Computer Applications Branch, ACD. 

Specific attention has been directed to 

The cooperative agreement was monitored by Dr. Robert E. 

DRECEDWG PAGE BLANK NOT FILMED 

iii 
CONSERVATIVE FINITE VOLUME SOLUTIONS OF A LINEAR HYPERBOLIC TRANSPORT 

EQUATION IN TWO AND THREE DIMENSIONS USING MULTIPLE GRIDS 

BY 

M. Kathongl and S.N. Tiwari2 

The feasibility of 

solving linear hyperbol 

cases. The results are 

the multiple gr 

c equations for 

compared with e 

ABSTRACT 

d technique is investigated by 

simple two- and three-dimensional 

act solutions and solutions obtained 

from the single grid calculations. It is demonstrated that the technique 

works reasonably well when two grid systems contain grid cells of comparative sizes. 

introduce any significant error and that it can be used to attack more cow 

plex problems. 

The study indicates that use of the multiple grid does not 

Graduate Research Assistant, Department of Mechanical Engineering and 

Eminent Professor, Department of Mechanical Engineering and Mechanics, 01 d 

Mechanics, 01 d Dominion University, Norfolk, Virginia 23529. 

Dominion University, Norfolk, Virginia 23529. 

iv 
TABLE OF CONTENTS 

. 

Page 

FOREWORD ............................................................... iv 

ABSTRACT ................................................................ v 

1 . INTRODUCTION ........................................................ 1 

2 . DISCUSSION ON CONSERVATIVE REZONING ALGORITHM ....................... 5 

3 . TWO DIMENSIONAL HYPRBOLIC EQUATION .................................. 9 

3.1 Formulation of 2D Hyperbolic Equation on Curvilinear Grid ..... 10 

Three-Stage Runge-Kutta Integration .................... 18 

Implementation on Multiple Grid System ........................ 19 

3.1.1 Difference Equation .................................... 14 

3.1.2 

3.2 

3.3 Results and Discussion ........................................ 19 

4 . THREE-DIMENSIONAL HYPERBOLIC EQUATION .............................. 24 

4.1 Formulation of 3D Hyperbolic Equation on Curvilinear Grid ..... 25 

4.1.1 Difference Equation .................................... 28 

4.2 

4.3 Results and Discussion ........................................ 39 

Implementation on Multiple Grid System ........................ 36 

5 . CONCLUSIONS ........................................................ 39 

REFERENCES ............................................................. 44 

LIST OF TABLES 

Table Page 

3.1 

3.2 

4.1 

Root Mean Square Errors. 20 Case 1 ................................ 22 

Root Mean Square Errors, 30 Case 1 ................................ 40 

Root Mean Square Errors. 2D Case 2 ................................ 23 

4.2 Root Mean Square Errors. 3d Case 2 ................................ 41 
TABLE OF CONTENTS (Concluded) 

LIST OF FIGURES 

Figure Page 

1.1 Zoning of multiply connected region .............................. 2 

2.1( a) Over1 apped Grids .............................................. 6 

2.l(b) Old Grid with known Qokl ...................................... 6 

2.l(c) New Grid with unknown Qn .................................... 6 

2.2 Overlap areas arising from the superposition of the old mesh 

(solid lines) and the new mesh (dash lines) . elementary overlap area is shaded ................................ 8 

3.1 Physical vs . Computational Domain ............................... 12 

3.2 A typical two-dimensional element ............................... 15 

3.3 Grid systems for a circular geometry using 11 x 11 and 

21 x 21 grids ................................................... 20 

Grid systems for a circular geometry using 21 x 21 and 

21 x 21 grids ................................................... 20 

Grid systems for an elliptic geometry using 11 x 11 and 

21 x 21 grids ................................................... 20 

ij 

A typical 

3.4 

3.5 

4.1 

4.2 

4.3 

4.4 An interface where ho is obtained by interpolating h from 

Physical vs . Computational Domain ............................... 26 

A typical volume element ........................................ 29 

A computational domain arising from two patched grids ........... 37 

Ramshaw's technique ............................................. 38 

A Butler-Wing configuration ..................................... 43 

the points indicated by dotted triangle and hn is obtained by 

5.1 

vi 
1. INTRODUCTION 

In recent years much progress has been made in the solutions of steady 

state equations of motion in both two and three dimensions. For complex 

shapes these calculations are usually based on a body-fitted curvilinear 

grid. For general three-dimensional bodies (for example, an aircraft 

configuration), it is very difficult to construct a body-fitted coordinate 

system [1]*. To simplify this problem, it is becoming more common to use 

several grids (multiple grids) at once, each in a different coordinate 

system [2-41. 

region at the interfaces of the various grids (Fig. 1.1). 

information be transferred from one grid to another accurately, it is important to treat grid points on the interfaces with care. 

nature of the equations of motion permits solutions with discontinuities 

such as shocks and slip surfaces. 

the right strength and physical location, it is imperative that the scheme 

used for the calculations be conservative [5]. In a multiple grid calculation, it is important that the interfaces are also treated in a conservative manner so that the discontinuities can move freely across these 

i nterf aces. 

This approach results in new boundaries within the given 

In order that 

The non-linear 

In order that such discontinuities assume 

The question of conservation when switching between two different grids 

or numerical schemes has been considered by several authors. Warming and 

Beam [6] derived transition operators for switching conservatively between 

MacCormack's method and a second order upwind scheme. 

[7] applied this transition operator approach to derive so-called zonal 

Hessenius and Pulliam 

*The nunbers in brackets indicate references. 
\ ZONE 9 

/ 

/ 

\ *ONES I ZoNE2 

Fig. 1.1. .Zoning of multiply connected region. 

2 
interface conditions, but with a significant loss of accuracy at the zonal 

interfaces. Rai [8] has developed conservative zonal interface conditions 

for zonal grids which share a common grid line, and has nice calculations 

demonstrating the shock capturing ability of zonal grids with a discontinuity crossing zones. Earlier work in the area of multiple grids includes 

that of Cambier -- et al. [9] who analyzed the zonal-boundary problem for a 

system of hyperbolic equations and used the compatibility equations to develop a zonal-boundary scheme. 

channel flow. 

zonal-boundary scheme that is not conservative and, hence, unsuitable for 

problem in which flow discontinuities move from one grid to another. Rai 

-- et al.[lo] present results obtained on metric discontinuous grids; 

integration scheme used is the Osher upwind scheme. 

results obtained on overlaid grids in conjunction with the stream function 

approach. 

201. 

Steger, and Doughtery [3]. 

Good results are presented for transonic 

However, the use of the compatibility equations results in a 

the 

Reference 11 gives 

For more information on multiple grid one should refer to [12- 

the need for conservative grid interfaces is illustrated by Benek, 

Dukowicz [21] described a rezoning method for arbitrary quadrilateral 

meshes in two dimensions. 

ilar to the method of Dukowicz, but is simpler and more direct. 

they defined rezoning as a method for transferring a conserved quantity, Q, 

from one generalized mesh to another when the volunetric density of Q is 

uniform within each cell of the original mesh. A computer program which 

follows Ramshaw's procedure has been written and tested with various types 

of grids and variables. The object i v e o f this study is to establish whether or not this technique is feasible 

for various grid systems while keeping the governing equation simple. 

Ramshaw [22] suggested a procedure which is simNote that 

The program seems to be working well. 

3 
For simplicity, the scalar 20 and 30 hyperbolic equations are chosen 

for this study. Hyperbolic equations have had a history of being used as 

model equations for testing newly developed schemes, for example, see [23- 

301. 

time integration. 

accuracy. 

Kutta to their nmerical calculations [38-421. 

Finite volume approach is chosen along with three-stage Runge-Kutta 

The three-stage Runge-Kutta integration is of 2nd order 

Many authors have applied the finite volume approach and RungeThis study can be divided into two portions. 

hyperbolic equations, qt + aqx + bqy = 0 is considered. 

is solved on a two dimensional grid system which is changed into another 

grid system at some time, say t = tl. The information obtained from the 

first grid calculation is transferred to the second grid by Ramshaw's technique. 

retain the conservative property of the nunerical scheme, the interpolation 

of flux across the interface is needed. q itself is the 

flux. Ramshaw's technique is used to transfer q across the interface. 

The second portion of this study is to examine the technique with the scalar 

3D hyperbolic equation, qt + aq, + bq,, + cq, = 0. 

suitable as the model equation of the equations of motion. Here, the grid 

system is changed from one to another at some x = y plane, say plane z = z 

Plane z1 is the interface mentioned previously. 

is going through this plane, i.e., flux in z direction = h = cq, needs to be 

transferred across the plane. 

doing so. 

equation is given in later sections. The results are compared with exact 

solutions and solutions obtained from the single grid calculations, i.e., 

without applying the Ramshaw's technique. 

First, the scalar 20 

Here, the equation 

Many authors, for example Berger [43], suggested that in order to 

For this problem, 

This equation is more 

1' 

For this problem, the flux 

Again, Ramshaw's technique is implemented in 

Details and procedures of solving both 2D and 30 hyperbolic 

It is important to note that 

4 
the objective of this study is to determine whether the technique is applicable for these simple governing equations. If it is applicable, one may 

have some confidence that it will be applicable for more complicated equations. It is expected that the results from this study will 

ficant contribution in the area of Computational Fluid Dynam 

In Sec. 2, a brief discussion on Ramshaw's technique is 

tions 3 and 4 are concerned about the two and three dimensio 

equation, respectively. Finally, the conclusion is given in 

yeild a signics. 

given. Seca1 hyperbol i c 

Sec. 5. 

2. DISCUSSION OW CONSERVATIVE REZOWIWG ALGORITHU 

Rezoning is a method for transferring a conserved quantity Q from one 

generalized mesh to another when the volumetric density of Q is uniform 

within each cell of the orginial mesh. 

in Sec. 1. 

dimensions has recently been described by Dukowicz [Zl]. 

suggested the procedure in doing so which is similar in spirit to the method 

of Dukowicz, but is simpler and more direct. 

this procedure has been written and is working well for example grids and a 

wide variety of choice of variables. 

The need for rezoning was discussed 

A rezoning method for arbitrary quadrilateral meshes in two 

Ramshaw [22] 

A computer program following 

By far the most common type of generalized mesh is the arbitrary quadrilateral mesh, which is convenient to work with because it has the same 

simple topological and logical structure as a square or rectangular mesh. 

the basic idea behind the rezoning is simple. Here, two grid systems are 

overlapped each other in some fashion (Fig. 2.1). 

denoted by Qokl, is to be transferred from the old grid system (Aokl is the 

area of each mesh in the old grid system) to the new grid system in which 

The conserved quantity, 

5 
1 I 

I 

I 1 

I I 

I 

* 

t---w---+---- - - -- I 

I 1 

I I 

I 1% 

I I 

I I 

- -- - 1 - -- - - - -1 - --- 

I I 

I I 

I I 

1 I 

1 I I 1 

Fig. 2.l(a). Overlapped Grids 

i 

Fig. Z.l(b) Old Grid with known Qokl 

I I 

I I 

I I 

I I 

1 I 

I I 

I I I ---- * ------- +---- I 

I I I 

I I I 

I I I 

I I 

I II I 

I I 

I I 

I I 

I 

I 

I 

I 

1 

F g. Z.l(c). New Grid with unknown - 

Qni 

6 
-- 

Anij is the area of each. 

quantity in each mesh of the new grid system. 

by : 

The quantity Qnij is denoted as the transferred 

Thus, Qnij can be computed 

where An is the portion of the area Anij which is contained in the area 

Aokl, and the sunmation is up to the nunber of the old meshes contained in 

Anij. The quantity qokl = Qokl represents the volumetric density of Qokl 

and is assuned to be uniform within each cell. The task, now, is to find 

Anokl and the number of the old meshes contained in each Anij. 

the polygon P is given by [44] 

Ok 1 

The area of 

where the sunmation is over all the sides of P, and E! is either +1 or -1 

depending on whether P lies to the left or right of side s. Note that the 

endpoint coordinates (xl, yf) and (xs, y 2) are considered to be associated 

with the side s and not with the particular polygon P. 

S S 

Figure 2.2 shows arbitrary overlap grids. 

gons whose sides are segments of the old-mesh lines. 

The overlap areas are polyThe number of sides of 

each type, and total nunber of sides, will be different for different overlap areas. Each side is common to two overlap areas, the one on the left 

(L) and the one on the right (R), and these overlap areas may be considered 

7 
Fig. 2.2. Overlap areas arising from the superposition 

of the old mesh (solid lines) and the new 

mesh (dash lines). 

overlap area is shaded. 

A typical elementary 

8 
to be associated with the side. 

The objective is to apportion a conserved quantity Q, whose volumetric 

density q is considered uniform within each cell of the old mesh, into the 

cells of the new mesh. It would be inefficient and difficult to automate in 

a computer, to naively sweep over the overlap areas directly. Instead, 

Ramshaw suggests to evaluate the same contributions by sweeping over the 

sides or segments s. The side or segment s is any side or segment of the 

polygon (overlap area). 

will be denoted by (xl, yl) and (xz, yz). 

The coordinate of the two endpoints of side s 

ss ss 

If the side s is a segment of the old mesh, then the quantity Q in 

the new-mesh cell containing side s is to be incremented by an amount: 

If the side s is a segment of the new mesh, then the contribution to cell 

N1 

2 

on its left is As = - qo (4, - qR) (xl, S yq) and (x;, y:) while the 

N contribution to the cell on its right is just -As, where q, is the 

volumetric density of the quantity Q of the old mesh cell in which side s 

lies. Adding A: and for each of the new mesh cell yields the quantity Q 

contained in each of the new mesh cells. 

N 

3. TWO-DIMEWSIOWAL HYPERBOLIC EQUATIOn 

Consider the two-dimensional hyperbolic equation, 

(3. la) 

9 
with the initial condition, 

(3.lb) 

The exact solution of Eq. 3.1 can be found as 

q(x,y,t) = f(x - at, y - bt) (3.2) 

This problem is well posed [45] if boundary conditions are specified at 

A the boundaries where the dot product, (a,b)* $ > 0; here (a,b) = ae, + 

b$ and $is the unit normal at the boundaries. 

multiple grid system using the finite-volune approach along with the threestage Runge-Kutta integration in time. Section 3.1 is concerned with the 

formulation of the 2D hyperbolic equation on the curvilinear grid system. 

The implementation on multiple grid system is briefly discussed in Sec. 3.2. 

Finally, the results and discussion on various grid configuration and initial condition are given in Sec. 3.3. 

This problem is solved on a 

Y 

3.1 Formulation of 2D Hyperbolic Equation on Curvilinear Grid 

Recall Eq. 3.1, 

qt + aqx + bqy = 0 

Since a and b are constants, one can write Eq. 3.1 in the conservation 

for as [46]. 

qt + fx + gy = 0 (3.2) 

where 

f = aq, g = bq 

10 
Equation 3.2 could be solved on any (x,y) coord nates, called physica 

coordinates. Instead, one normally solve Eq. 3.2 on the orthogonal 

coordinates system, called computational coordinates [47]. 

3.2 needs to be transformed to the computational coordinates. 

the transformation into the computational coordinates produces some 

additional terms in the transformed partial differential equation (thus 

increases the amount of computation required), it has some advantages. 

Boundary conditions can be treated in the straight forward way. Also, 

Thus, equation 

Even though 

transformed equations are solved in the simple orthogonal and possibly equal 

spaced region. Severe departure from orthogonality will introduce 

truncation error in difference expressions [47]. Figure 3.1 illustrates the 

physical versus computational domain in a two-dimensional space. 

the chain rule, Eq. (3.2) is written as: 

By using 

Also using the chain rule, one can write: 

dx =x + x dn 5 n 

dy = yE& + yndn 

which can be written in the matrix form as: 

11 
Physical Domain 

Fig. 3.1. Physical vs. Computational Domain 

Computational Domain 

12 
The inversion of the matrix yields: 

Simi 1 ar 1 y, 

From Eqs. 3.4 and 3.5, it can be seen that: 

where J= 1 

"Yn - y€xn 

(3.4) 

(3.6a) 

is the Jacobian of transformation. obviously, Eq. 3.6 exists when J f - or 

1-x -[E x .] f 0. 

YE yn 

E = -& Y n From Eq. 3.6 it is found that 5, = Jyn 

(3.6b) 

13 
Equation 3.3 can be written as: 

- 1 - 1 

(3.7) 

q = q/J 

Equation 3.7 is the transformed equation which is solved on the 

computational coordinates. 

3.1.1 Difference Equation 

In this section, the finite volme approach is applied to Eq. 3.7 along 

with three-stage Runge-Kutta integration in time. 

difference equations on the equally spaced in each direction (5 - TI), 

computational domain is given below. 

The formulation of the 

Consider a computational cell illustrated in Fig. 3.2. Recall Eq. 3.7 

qt +- fE + g, = O 

14 
I I 

a 

Fig. 3.2. A typical two-dimensional element. 

15 
Integrating over the cell yields 

qt 

Using the divergence theorem, one obtains 

i 

Applying this to Fig. 3.2, one can write 

2 2 

An 2 2 

(3.8) 

16 
Each term on the right hand side of Eq. 3.8 can be written as 

2 2 

An 2 

2 2 

L 2 2 2 2 2 2 

2 L 2 2 2 

2 2 2 2 31 

17 
t) = f(xl - at, y1 - bt), are 1’ yl’ The exact boundary conditions , q( x 

given where the product (a, b) n > 0, and the extrapolated values of q 

are given at the boundaries where (a, b) n < 0. This is due to the nature 

of the hyperbolic equations. 

3.1.2 Three-S t age Runge- Kut t a I n tegr at ion 

Recall that Eq. 3.8 is written as 

2 2 

where 

2 2 

By applying three-stage Runge-Kutta, one may express: 

qn + 1* = qn + Atf(qn) 

1 n + 1** n n + 1* 

q = q + - 1 Atf(qn) + - 1 Atf(q 

2 2 
3.2 Implementation on Multiple Grid System 

As mentioned previously, for multiple grid system, the interpolation of 

flux will yield a conservative scheme. For the scalar 20 hyperbolic 

equation, Eq. 3.1, q itself is the flux transporting in time. For this 

problem, the difference equations as described in Sec. 3.1.1 are solved on 

one grid system from initial time t = 0 up to some time, say t = Then, 

at t = tl the grid system is changed and the difference equations are solved 

on this new grid system. The initial conditions for the new grid system are 

obtained by the interpolation of the final values in the old grid system, 

i .e., q(x,y,tl). 

implemented in doing so. 

exact solutions and with the solutions obtained by using the new grid system 

for all time, i.e., from t = 0 to t = final without changing the grid 

system. 

5' 

Ramshaw's technique, as described in Sec. 2, is 

The results at t = tfinal are compared with the 

3.3 Results and Discussion 

Various grid systems, as shown in Figs. 3.3 - 3.5, and initial 

conditions are implemented as described in Sec. 3.2. 

compared with the exact solutions and the solutions obtained by using the 

second grid system for all time. 

The results are 

For the case of exponentially stretching, 

k2 -1 

e 

kl -1 e 

the derivative contained in the Jacobian of transformation, J, xs, xn, etc. 

can be computed analytically. Here, results using both analytically 

19 
Fig. 3.3. Grid systems for a circular geometry using 11 x 11 and 

21 x 21 grids. 

Fig. 3.4. Grid systems for a circular geometry using 21 x 21 and 

21 x 21 grids. 

Fig. 3.5. Grid systems for an elliptic geometry using 11 x 11 and 

21 x 21 grids. 20 
computed derivatives and nunerically computed derivatives are shown. 

Before discussing any results, it should be noted that some error 

should be expected since Ramshaw's technique assumes constant volumetric 

density of the conserved quantity in each grid cell. But in this study, the 

analytic initial conditions which vary within grid cells have been selected. 

However, the objective of this study is to determine whether the errors grow 

significantly as Ramshaw's idea is implemented. 

Tables 3.1 and 3.2 show the root mean square errors at final time, 

t = t If errors are defined as the 

difference between the calculated solutions and the exact solutions. the 

root mean square errors can be obtained as 

(both tables are at 200 time steps). 

where rn x n are the number of cells contained in the grid system where the 

errors are computed. 

Tables 3.1 and 3.2 indicate that the technique does not introduce 

significant errors to the solutions. 

can be seen clearly by looking at the top row of Table 3.2. In the case of 

constant mesh size with linear initial conditions, the finite volume calculation should yield the exact solutions. Note that multiple grid calculations even yield smaller errors than the single grid calculation in some 

cases. 

significant errors to the solutions, especially when it is implemented on 

two grid systems which contain grid cells of comparative sizes. 

The errors introduced by the technique 

Thus, it can be concluded that the technique does not introduce 

21 
Table 3.1 Root Mean Square Errors, 2D Case 1 

initial conditions, q = ~cos(lx)cin(By) 

exact 601ut ions, 

2 

Q = ~cos(l(x-at))rin(l(~-bt)) 2 

1 Orid vs. 

e grid 

mxn 

1 Nrid 

const ant 

mesh size 

expo anal. 

kll.5 num. 

rig 3.3a vs. 3.3b 

'ig 3.4a vs. 3.4b 

:ig 3.5a vs. 3.5b 

11x11 

llxll 

11x11 

llxll 

21x21 

11x11 

mxn 

2 grid 

21x21 

21x21 

21x21 

21x21 

21x21 

21x21 

rms. error 

1 grid 

. 0110307 

. oeaess 

,02015 

02 16039 

.002740 

.0241553 

rms. error 

2 grid 

. 0027907 

.OW452 

.0074 14 

. 009214 

.OO92 14 

. 0099474 

rmr. error 

multiple grid 

. 0029005 

.007656 

. 00762 

. oio2a46 

.009668 

. 0 108937 

22 
Table 3.2 Root Mean Square Errors, 2D Case 2 

1 grid vs. 

2 grid 

const ant 

mesh size 

expo anal. 

k-1.S num. 

3g 3.3b vs. 

Fig 3.3a 3.4a 

Initial conditions, q = Wx+y 

exact rolut ions, q 

mxn 

1 gria 

11x11 

11x11 

11x11 

11x11 

21x21 

11x11 

- 

mxn 

2 grid 

21x21 

21 x21 

21x21 

21x21 

21x21 

21x21 

rms. error 

1 grid 

.0000001 

.0067252 

.0064324 

. 0025948 

,00065 18 

0052474 

rmr. mrror 

2 grid 

0000000 

0017874 

. 00 I6984 

.0052553 

. 0052553 

,0042585 

rmr. error 

multiple grid 

~~ 

.0009590 

.0020J91 

.0019984 

.0075165 

.00532 13 

.0043775 

23 
4. THREE-DIMENSIONAL HYPERBOLIC EQUATION 

Consider the three-dimensional hyperbolic equation, 

+ aqx + bqy + cq, = 0 0 <x <1, 0 <y<l, 0 <z<l qt 

with the initial condition, 

The exact solution of this equation can be found as 

q(x, y, z, t) = f(x-at, y-bt, z-ct) 

As mentioned in Sec. 3, this problem is well posed if one specifies the 

-a boundary conditions at the boundaries where the dot product ( a,b,c) n > 

0, where (a,b,c) = a? +b$ +cT and is the unit normal at the boundaries. 

XYZ 

This problem is solved on the multiple grid system using the finitevolume approach along with the three-stage Runge-Kutta integration in time. 

In Sec. 4.1, the formulation of the 3D hyperbolic equation on the curvilinear grid system is described. 

the multiple grid is given in Sec. 4.2. 

from various grid configurations and initial conditions. 

hyperbolic equation is more suitable to test Ramshaw's technique. 

because most real physical computations are performed on a three dimensional 

space and the governing equations are similar to the 3D hyperbolic 

equation. 

The brief discussion'on the implementation on 

Section 4.3 discusses the results 

Note that the 3D 

This is 

24 
4.1 Formulation of 3D Hyperbolic Equation on Curvilinear Grid 

Cons 

which can 

der the three dimension1 hyperbolic equation, 

be written in the conservative form as [46] 

qt + fx + gy +hz=O 

(4.1 1 

(4.2 1 

where 

f = aq, g = bq, h = cq 

and a, b, c are constants. 

This section shows the formulation of the transformed equation of Eq. 

4.2 on the computational coordinates (E-TI-<). 

forming Eq. 4.2 into the computational coordinates were mentioned in Sec. 

3.1. 

The advantages of transFigure 4.1 illustrates the physical versus computational domain. 

By using the chain rule, Eq. 4.2 is written as 

+ +hr)=O. (f<<x + gr<y < z 

After some mathematics manipulations, one can write 

(1 q) + (15 f + 15 g + 15 h) + (hxf + 1n g + 1n h) -x --Y " -2 E - -.Y - -z n 

J J J J J J J 

25 
Physical Domain Computational Domain 

Fig. 4.1. Physical vs. Computational Domain. 

26 
+ (1s f + 1s g + IS h) -x -y -z 5 

J J J 

= 0 

I where 

(4.3) 

5, sy sy 

"x "y "z 

5, cy <z 

Similar to Eq. 3.6, one can see that 

1 

J 

- 

1 

J 

- 

Equation 4.3 can be written as 

27 
where 

A 

q=q 

J 

Equation 4.4 is the transformed equation which is solved on the computat ional coordinate (6 ,n ,c). 

4.1.1 Difference Equation 

Figure 4.2 illustrates a volume element on the computational coordinates. 

approach to Eq. 4.4 with the volume element as in Fig. 4.2. 

the formulation are given below. 

The difference equations are derived by applying the finite volume 

The details of 

Recall Eq. 4.4, 

$t + + Cn +% = 0 

< 

or 

28 
i + 1, 

I 

8 

j + 1, k+ 1 

Fig. 4.2. A typical volume element. 

29 
Integrating over the volune element, it is expressed as 

f ‘W t 5 

dV = -JW (?E +;n +$’ )dV 

Applying the divergence theorem to the right hand side, the above equation 

can be written as 

Utilizing Fig. 4.2, the difference equation for the above equation can be 

written as 

A 2 

-h i+1, - j+l, - k 

2 2 

(4.5) 

30 
X 

rl i, j +f 

Yrl + i, j +L 

2 

- j + 1, k - 'i, j, k 

All 2 

+x -x 1 i, j + 1, k + 1 i, j, k +1 

An 

- '('i, j + 1, k - 'i, j, k 9 L 

Arl 

1 + 'i, j + 1, k + 1 

All 

All L 2 2 I 

-z 1 , j + 1, k+ 1 i, j, k+l + zi 

All 

X 5 i, + 1, 

2 

- 9 

2 

+ 

-x k+l = - j, k + 1 i, j, k 

All 

2 2 

+x -x 1 i, j+1, k+l i, J+1, k 

All 

k+ 1 - = - '('i, j, k + 1 'i, j + 1, k 

Ac 2 2 

1 - + 'i, j + 1, k + 1 

AT; 

'i, j, k 

= l(Zi -z zl si,j+l,k+l - - - n , j, k + 1 i, j, k 

As L 2 2 I 

+z -2 1 i,j+l,k+l i,j+l,k 

31 
A 

qi = [ - (yg5 - z y )f + (x z + - 1,j, k + - 1 55 5 5 - ZSXr)g 

5 X 

2 2 

= l(Xi -x 

i+L, j,k+l - j, k + 1 i, j, k 

2 

2 2 

'C 

-x 1 i + 1, j, k + 1 i, j, k +1 +x 

- 

i +L, j, k+l = - '('i, j, k + 1 'i, j, k 

2 

'2 2 A€ 

) - + 'i + 1, j, k + 1 'i, j, k +1 

AE 

-z ) i + 1, j, k + 1 i, j, k +1 +z 

AE 

+X -x ) i+l, j, k+l i+l, j, k 

AC 

1 - 

+ 'i + 1, j, k + 1 'i + 1, j, k 

32 
5 i + 1, j, k + 1 =L('i, j, k + 1 - 'i, j, k -Q 

YE 

1 A5 

L '2 2 

-z i + 1, j, k+ 1 i + 1, j, k +z 

A5 

- - 

i +1, j +L, k - - '('i + 1, j, k 'i, j, k 3 

2 2 

= l(Xi -x 

%i+l,j+l k - , j + 1, k i, j, k 9 -9 

2 2 

-x - 

j + 1 k +1, j, k i, j& 3 -9 

I AE L 2 .2 

-x 1 i, j + 1, k +x i+l, j+l, k 

) ''i + 1, j + 1, k - 'i, j + 1, k 

A€ 

+z -z 1 i+l, j+l, k i, j+l, k 

A6 

All L I2 2 

+x -x 1 i+l, j+1, k i+l, j, k 

An 

33 
- 

’ni+l,j+l,k l- - - - ‘(Yi, j + 1, k - yi, j, k CI 

- 1

J 

An L ‘2 2 

=[x(yz -zy)-x(yz -2Y1 i+l, - - j+l,k+l -E 05 nt n 55 55 

2 2 2 

+’i 1 + 1, j + 1, k - ’i + 1, j, k 

= l(Zi -2 ‘ni+l,j+l l- -9k - c1 ,j+l,k i,j,k 

L ‘2 2 

+x -x ) i+l, j+l, k i+l, j, k 

f = aq, g = bq, h = cq 

2 2 2 2 2 2 2 2 

2 2 2 2 2 2 2 2 L 

) - 

‘i+l, -j+l, - k i(qi+l, - j+l, - k+l+qi+l, - - j+l, - k-1 - 

2 2 2 2 2 2 2 2 

2 2 2 

34 
X 

5 

5-1 

i 

i+l, 

+ 1, j 

2 

+ 1, k 

2 

- 

+ 1 -1. + 1, j, k - 'i, j, k 

A5 

- 22 2 

+x -x 1 i+l, j+l, k i, j+l, k 

A6 

-x 1 

-x 11 

+L('i + 1, j, k + 1 i, j, k +1 

AE 

2 

+x i+l, j+l, k+l i,j+l, k+l 

AE 

j + 1, 

2 

j + 1, 

2 

= 1 [l(x -x 

k+l -- i, j + 1, k i, j, k 

All 

- 22 2 

+x -x 1 i+l, j+l, k i+l, j, k 

All 

-x 1 + - i, j, k +1 ('i, j + 1, k + 1 

All 

2 

+x -x 11 i+l,j+l,k+l i+l,j,k+l 

All 

= 1 [ l(Xi -x 

k+L - - Ak+l 

A5 22 2 

+'i+l, j, k+l -x i + 1, j, k 1 

1 i, j + 1, k 

-x + 1 ('i, j + 1, k + 1 

2 

+x -x 

i+l,j+l,k+l i+l,j+l,k 

and z . Finally, Yn9 Y5' z 5' zl19 5 

Similar expressions can be written for y 

Eq. 4.5 can be written in the form: 

35 
The three-stage Runge-Kutta integration as described in Sec. 3.1.2 can be 

appl i ed here. 

4.2 IMPLEMENTATION ON MULTIPLE GRID SYSTEM 

Figure 4.3 illustrates an example of multiple grid used in this study. 

Here, the grid system changes from one to another at the plane z = zl. In 

order to solve the difference equations, as in section 4.1.1, on this grid 

system, the extra boundary conditions are needed at the interface, i.e., 

plane z = z1 where the two single grid system meet. 

tions are not the physical boundary conditions. Many authors refer to them 

as interface conditions, for example see Berger [43]. 

mush be given for these boundary conditions. 

interface in this study. 

These boundary condiSpecial treatments 

Figure 4.4 shows the typical 

In Fig. 4.4, ho denotes the flux h coming out from the first grid 

system and hn denotes the flux h going into the second grid system. Here, 

ho is needed for the first grid calculations and hn is needed for the calculations on the second grid. 

mentioned previously. 

between the first and the second grid. 

ferring ho across the interface using the Ramshaw's technique described in 

Sec. 2. The calculations begin from initial time, t = 0, to some time, say 

t = tfinal. 

with the exact solutions and the solutions obtained from the single grid 

calculations for various grid system and initial conditions. 

They are the extra boundary conditions 

In this study, ho is obtained by the interpolations 

The flux h is obtained by transn 

The results from the multiple grid calculations are compared 

36 
Interface 

----------- I 

I I 

z = z1 

I 

I 

I 

I 

I 

I 

-4 

1 

1 

Fig. 4.3. A computational domain arising from two 

patched grids. 

37 
I/ 

I I I 

I 

I 

Fig. 4.4. An interface where ho is obtained by 

interpolating h from the points indicated by dotted triangle and hn is obtained by Ramshaw's technique. 

38 
4.3 RESULTS AND DISCUSSIONS 

Tables 4.1 and 4.2 show the results obtained from 

various grid systems with two different initial condit 

the calculations on 

ons. If the error is 

defined as the difference between the numerical solutions and the exact 

solutions, the root mean square error can be found as 

rms. error ;jaer;orz 

where N = total number of grid cells. 

Both tables illustrate the rms. errors from the multi le rid c 1 Ulations as well as the calculations from each of single grid which compose 

the multiple grid. Table 4.1 gives the results with linear initial conditions, while the results with trigonomitrical initial conditions are given 

in Table 4.2. 

lation does not introduce significant errors to the over all solutions. 

some cases, the multiple grid calculations even give better solutions than 

the calculations from the single grid. 

It can be seen from both tables that the multiple grid calcuIn 

5. CONCLUSIONS 

The results from this study have shown that the technique works reasonably well when the two grid systems contain grid cells of comparative sizes. 

This can be clearly seen when one recalls the assumptions made for the technique, i.e., each grid cell contains uniform conserved quantity. Since, in 

general, this criterion cannot be met in the real physical calculations, 

some effort may need to be made to seek for the best possible types of grids 

in order to apply the technique. It should be pointed out that only simple 

39 
Table 4.1 Root Mean Square Errors, 3D Case 1. 

initial conditions, q = R+x+y+z 

exact solutions, q = R+(x-at)+(y-bt)+(z-ct) 

conf i gur-at ion 

cunst. mesh 

?xp. streching 

with k=3.0 

:*=const. mesh) 

cylinder 

E u t 1 er -w i n g 

-~ 

f, No. of pts. 1 vs. 2’ 

llxllxll vs. 

llx21xll vs. 

17x17~11 vs. 

2lx21x11 vs. 

L9xl9xll VS. 

>lX2lXll VI. 

~1x21xl1 vs. 

1lxi21xll vs. 

2lx21xll 

11x1 1x1 1 

2lx2lxll 

17x17~11 

2lx21xll 

l9x19xll 

2lx2lxll 

21X~lXll 

O-type vs. H-type 

l1x2lx1l vs. 21X2lX11 

llxllxll Vf. 21X21Xll 

17x17~11 vs. 21x21~11 

llx21xll vs. Ellx2lxll 

rms. error 

single grid 

first 

IO00000 

.000000 

t 000000 

,000000 

,000000 

,000000 

, 000000 

,002 106 

000608 

002842 

001495 

00 1275 

second 

. 000000 

. 000000 

0 000000 

,000000 

.000000 

. 000000 

.002 106 

. 000000 

. 0 15802 

,012321 

,012321 

,012321 

mul. grid 

.013660 

.000004 

.005720 

. 003334 

s 002 163 

.001703 

,002389 

,007 139 

,009620 

,012083 

,008672 

,000379 

40 
Table 4.2 Root Mean Square Errors, 3D Case 2 

llxllxll vu. 2ix21x1: 

2lxESlxll VS. llxllxll 

17x17~11 vs. 2lx21xll 

21x21~11 VS. 17xl7xll 

19xl9x11 VI. 21X21Xll 

21X21X11 V6. 19X19Xll 

2 1 x2 1 x 1 P vs. 

22n21x11 vs. 21X21Xli 

21 x21 x 1 1 

O-type VS. H-type 

21n21xll vs. 21x21~11 

llxllxll vs. 21X21Xll 

17x17~11 V8. 21X21X11 

21X21X11 V6- 21X21X11 

rms. error 

first 

.00556C 

.001901 

.002494 

.001901 

,002136 

.001901 

IO01901 

,004424 

,001651 

007848 

003297 

002180 

6econd 

001901 

.005560 

.001901 

002494 

001901 

,002 136 

-004424 

.001901 

,021766 

,021 175: 

,021 1752 

I 021 1752 

mul. grid 

. 0 19420 

.003786 

.006363 

,004219 

.003414 

003070 

.00557B 

008227 

,011515 

.017058 

. 0 12327 

,011819 

41 
configurations are chosen f o r t h i s study, i.e., a single grid system can be 

generated around these configurations. 

multiple grid approach may be the only way to attack the problem if the 

complexity of the geometries is to be maintained. 

For more complex bodies, the 

This study demonstrates that use of the multiple grid does not introduce any significant error to the problems. 

technique to solve the equations of motion over an ideal aircraft configuration such as a Butler-Wing configuration (Fig. 5.1). If the results from 

the Butler-Wing calculations are satisfied, the next step is to 

consider this technique along with the equations of motion over the real 

aircraft configuration. It is expected that results from this study will 

give some contributions to the field of Computational Fluid Dynamics. 

The next step is to apply this 

42 
Fig. 5.1. A Butler-Wing configuration. 

43 
REFERENCES 

I 

I Paper 83-1929, July 1983, pp. 293-302. 

1. Jameson, A. and Baker, T.J., "Solution of the Euler Equations for 

Complex Configurations," AIAA 6th Comput. F1 uids Dynamics Conference, 

2. Atta, E., "Component-Adaptive Grid Interfacing," AIAA 19th Aerospace 

Sciences Meeting, Paper 81-0382, January 1981. 

3. Benek, J.A., Steger, J.L., and Dougherty, F.C., "A Flexible Grid 

Embedding Technique with Application to the Euler Equations," AIAA bth 

Comput. Fluids Dynamics Conference, Paper 83-1944, July 1983, pp. 373- 

382. 

4. Boppe, C., "Calculation of Transonic Wing Flows by Grid Embedding," 

AIAA 15th Aerospace Sciences Meeting, Paper 77-207, January 1977. 

5. Lax, P., "Hyperbolic Systems of Conservation Laws and the Mathematical 

Theory of Shock Waves," SIAM, 1972. 

6. Warming, R. and Beam, R., "Upwind Second-Order Difference Schemes and 

Applications in Aerodynamic F1 ows," AIAA Journal Volume 14, September 

1976, pp. 1241-1249. 

7. Hessenius, K. and Pulliam, T., "A Zonal Approach to Solution of the 

Euler Equations," AIAMASME 3rd Joint Thermophysics, F1 uids, P1 asma and 

Heat Transfer Conference, Paper 82-0369, June 1982. 

8. Rai, M.M., "A Conservative Treatment of Zonal Boundaries for Euler 

Equation Calculations," AIAA 22nd Aerospace Sciences Meeting, paper 84- 

0164, January 1984. 

9. Cambier, L., Ghazzi, W., Veuillot, J.P., and Viviand, H., "Une 

Approache par Domaines pour le Calcu d'Ecoulements Compressibles," Cinquieme Col loque International sur les Methodes de calcu, France, 

December 1981. 

10. Rai, M.M., Chakravarthy, S.R., and Hessenius, K.A., "Zonal Grid 

Calculations Using the Osher Scheme," International Journal of 

Computers - and Fluids, Volume 12, No. 3, 1984 , pp. 16r-TT5-r11. Steger, J.L., Dougherty, F.C., and Benek, J.A., "A Chimera Grid 

Scheme," Mini-Symposium on Advances in Grid Generations, ASME Applied, 

Bioengineering and Fluids Engineering Conference, Houston, Texas, June 

1983. 

12. Beam, R. and Warming, R., "An implicit Factored Scheme for the 

Compressible Navier-Stokes Equations," AIAA Journal, Volune 16, April 1978, pp. 393-402. 

13. Atta, E. and Vadyak, J., "A Grid Interfacing Algorithm for ThreeDimensional Transonic F1 ows and Aircraft Configurations," AIAMASME 3rd 

44 
14. 

15. 

16. 

17. 

18. 

19. 

20. 

21. 

Joint Thermophysics, F1 uids, P1 asma and Heat Transfer Conference, Paper 

82-1017, June 1982. 

Osher, S. and Sanders, R., "Numerical Approximations to Nonlinear 

Conservation Laws with Locally Varying Time and Space Grids," Math. 

Comp., Volume 41, September 1983, pp. 321-336. 

Berger, M.J., "Stability of Interfaces with Mesh Refinement," 

ICASE/NASA/Langley Research Center, Hampton, Vi rgini a, Report No. 83- 

42, August 1983. 

Trefetnen, L.N., "Stability of Finite Difference Models Containing Two 

Boundaries or Interfaces," ICASE NASMLangley Research Center, Hampton, 

Virginia, Report No. 84-11, March 1984. 

Berger, M.J. and Oliger, J., "Adaptive Mesh Refinement for Hyperbolic 

Partial Differential Equations," Journal of Computational Physics, 

Volume 53, March 1984, pp. 484-512. - 

Browning, G., Kreiss, H.-O., and Oliger, J., "Mesh Refinement ,I' Math. 

Comp., Volume 27, January 1973, pp. 29-39. 

Cement, M., "Stable Difference Schemes with Uneven Mesh Spacings," 

Math. Comp., Volume 25, April 1971, pp. 219-226. 

MacCamy, R.C. and Martin, S.P., "A Finite Element Method for Exterior 

Interface Problems," Int. J. Math. Math. Sci., Volume 3, 1980, pp.311- 

350. 

----- 

Dukowicz, J.K., "Conservative Rezoning (Remapping) for General 

Quadrilateral Meshes," Journal - of Computational Physics, Volume 54, 

June 1984, pp. 411-424. 

22. Ramshaw, J.D., "Conservative Rezoning Algorithm for Generalized TwoDimensional Meshes," Journal - of Computational Physics, Volume 59, June 

1985, pp. 193-199. 

23. Harten, A., Lax, P., and Leer, B., "On Upstream Differencing and 

Gudunov-Type Schemes for Hyperbolic Conservation Laws," ICASE 

NASA/Langley Research Center, Hampton, Virginia, Report No. 82-5, March 

1982. 

24. Leer, B., "Multidimensional Explicit Difference Schemes for Hyperbolic 

Conservation Laws," ICASE NASWLangley Research Center, Hampton, 

Virginia, Report No. 172254, November 1983. 

25. Beam, R. and Warming, R., "An Implicit Finite Difference Algorithm for 

Hyperbolic Systems in the Conservation Law Form," Journal of 

Computational Physics, Vol me 22, September 1976, pp. 87-110. 

26. Turkel, E., "Composite Methods for Hyperbolic Equations," SIAM Journal 

of Numerical Analysis, Volume 14, September 1977, pp. 744-759. 

45 
27. Starius, G., "On Composite Mesh Difference Methods for Hyperbolic 

Differential Equations," Numerical Math., 1980, pp. 241-255. 

28. Oliger, J., "Fourth-Order Difference Methods for the Initi al-Boundary 

Value Problem for Hyprbolic Equations," Math. Comp. Volume 28, January 

1974, pp. 15-25. 

29. Oliger, J., "Hybrid Difference Methods for the Initial-Boundary Value 

Problem for Hyperbolic Equations," Math. Comp. Volume 30, October 1976, 

pp. 724-738. 

30. Harten, A., "High Resolution Schemes for Hyperbolic Conservation Laws," 

Journal - of Computational Physics, Volume 49, March 1983, pp. 357-393. 

31. Gottlieb, D., Gunzburger, M., and Turkel, E., "On Numerical Boundary 

Treatment of Hyperbolic Systems for Fi ni te Di fference and Fi ni te 

Element Methods ,I1 ICASE NASA/Langley Research Center, Hampton, 

Virginia, Report No. 78-13, September 1978. 

32. Abarbanel, S.S. and Murman, E.M., "Stabi 1 ity of Two-Dimensional 

Hyperbolic Initial Boundary Value Problems for Explicit and Implicit 

Schemes," ICASE NASAILangley Research Center, Hampton, Virginia, Report 

NO. 81-29, August 1981. 

33. Trefethen, L.N. and Halpern, L., "We1 1-Posedness of One-way Wave 

Equations and Absorbing Boundary Condition," ICASE NASWLangley 

Research Center, Hampton, Virginia, Report No. 85-30, June 1985. 

34. Osher, S. and Tadmore, E., "On Convergence of Difference Approximations 

to Scalar Conservation Laws," ICASE NASWLangley Research Center, 

Hampton, Virginia, Report No. 85-28, May 1985. 

Goldberg, M. and Tadmor, E. "Scheme-Independent Stability Criteria for 

Difference Approximations of Hyperbolic Initial-Boundary Value 

Problems. II.," Math. Comp. Volume 36, April 1981, pp. 603-626. 

Osher, S., "Stability of Difference Approximations of Dissipative Type 

for Mixed Initial-Boundary Value Problem. I. ,I' Math. Comp., Volune 23, 

April 1969, pp. 335-340. 

35. 

36. 

37. Trefethen, L.N., "Instability of Difference Models for Hyperbolic 

Initial-Boundary Value Problems," Communications on Pure and Applied 

Math., Volume 37, May 1984, pp. 329-367. 

38. Moitra, A., Turkel, E., and Kumar, A., "Application of a Runge-Kutta 

Scheme for High-speed Inviscid Internal F1 ows," ICASE NASWLangley 

Research Center, Hampton, Virginia, Report No. 86-5, January 1986. 

39. Turkel, E., "Accuracy of Schemes with Nonuniforms Meshes for 

Compressible F1 uid F1 ows," ICASE NASWLan ley Research Center, Hampton, Virginia, Report No. 85-43, September 198 8 . 

46 
40. 

41. 

42. 

43. 

44. 

45. 

46. 

47. 

Turkel, E., Yaniv, S., and Landau, U., "Accuracy of Schemes for the 

Euler Equations with Nonuniform Meshes ,I' ICASE NASA/Langley Research 

Center, Hampton, Virginia, Report No. 85-59, December 1985. 

Jamesons, A., Schmidt, W., and Turkel, E., "Numerical Solutions of the 

Euler Equations by Finite Volume Methods Using Runge-Kutta Time 

Stepping Schemes ,'I AIAA 14th F1 uids and P1 asma Dynamics Conference, 

Paper 81-1259, June 1981. 

Moitra, A:, "Numerical Solution of the Euler Equations for High-speed, 

B1 ended WI ng-Body Configurations," AIAA 23rd Aerospace Sciences 

Meeting, Paper 85-0123, January 1985. 

Berger, M. J., "On Conservation at Grid Interfaces ,I' ICASE NASA/Langley 

Research Center, Hampton, Virginia, Report No. 84-43, September 1984. 

Rektorys, K. (Ed.), "Survey of Applicable Mathematics," p. 207, - MIT 

Press, Cambridge, Mass. 1969. 

Erickson, L.-E., "Transfinite Mesh Generation and Computer-Aided 

Analysis of Mesh Effects, Ph. D. Thesis, Department of Computer 

Sciences, Uppsala University, Uppsala, Sweden. 

Anderson, D.A., Tannehill, J.C., and Pletcher, R.H., "Computational 

F1 uid Mechanics and Heat Transfer," Hemisphere Publishing Corporation, 

New York, New York, 1984. 

Thompson, J.F., Warsi, Z.U.A., and Mastin, W.C., "Numerical Grid 

Generation," North-Holland, New York, New York. 

47