FINAL REPORT 

In response to the requirements of 

NCC 2-515 

EFFICIENT GRID GENERATION 

Principal Investigator: Ryoichi Seki 

Period of Study: August 1, 1987 - July 31, 1989 

(NASA-CB-18 SO 34 ) 

Final Report, 1 Aug,,1987 - 31 Jul, 1989 

EPF ICXBBT 631 D GEM BR AT IOB 

(California State Uaiv,i, 126 p CSCL 12A 

d 89- 274 34 

Unclas 

G3/64 0224988 

CSUN 

Northridge, CA 91330 
Contents 

Preface 

List of Figures 

Abstract 

... 

111 

vi 

... vlll 

1 Introduction 1 

2 Classification of Partial Differential Equations 9 

3 Transformation Metrics 11 

4 Elliptic Grid Generation 13 

5 Hyperbolic Grid Generation 23 

5.1 Numerical technique ......................... 27 

6 Combining Elliptic and Hyperbolic Equations 36 

6.1 Solution Procedure .......................... 37 

6.2 Testing Procedure ........................... 39 

6.3 The Unmodified Algorithm ...................... 40 

6.4 Floating the Boundary ........................ 41 

6.5 Using Compressed-Cell Initial Grids ................. 42 

6.6 Curtailing Outer Boundary Movement within a Cavity ...... 43 

7 Grid Breakdown 94 

8 Derivation of the Navier-Stokes Equations 99 

9 Results and Conclusions 107 

-- iv 

PRECEDING PAGE BLANK NOT FiLMED 

3F ,a 
10 Notation and Definitions 

3&Xoo*A-- O*"i/i.J 

113 

i13 

V 
1 

2 

3 

4- 7 

Portion of boundary points and solution points in a 2-dimensional 

rectangular coordinate system with Ax and Ay constant. Boundary conforms to a coordinate line. .................. 

Variance in Ax and Ay to increase local solution accuracy and 

resolution with resultant additional number of solution points. .. 

Consequent change in Ax and Ay near boundary that doesn't conform to a coordinate line.. ...................... 

11" 

6 

7 

8 

The elliptic generator at work .................................. 19-22 

8-14 The hyperbolic generator at work on bodies of various concavities 

.............................................................. 29-36 

15-19 The combined elliptic-hyperbolic generator with a fixed outer 

boundary .................................................... 46-50 

20-33 The combined elliptic-hyperbolic generator with a floated outer 

boundary .................................................... 51-64 

34-45 The combined elliptic-hyperbolic generator with a floated outer 

boundary, using cell-volume growth. .......................... 65-76 

46 The combined elliptic-hyperbolic generator with a floated outer 

boundary, using cell-volume growth. Failure of technique. ....... .. 77 

47-62 The combined elliptic-hyperbolic generator with a floated outer 

boundary, using cell-volume growth. Forming a cut. .......... 78-93 

63 Portion of boundary and first row of grid cells next to boundary. 

Here k=2 ............................... 95 

64 Grid breakdown at a strong convexity. ............................. 98 

vi 
ORIGINAL PAGE 1% 

OF POOR QUALITY 

65 Fixed control volume showing area vector x and the three surface 

,.. stresses at one oi sac: EWO sun’aces :-uA~--, .. .._v ). ... ,i,. ...... vj 

66-69 Generating a grid about two bodies ........................ 110-112 

70-73 The relation of the parameter cam to concavity ................. 117 

vii 
ABSTRACT 

EFFICIENT GRID GENERATION 

by 

Bill Gutierrez 

Master of Science in Physics 

Because the governing equations in fluid dynamics contain partial differentials 

and are too difficult in most cases to solve analyticdy, these differentials are 

generally replaced by finite difference terms. These terms contain terms in the 

solution at nearby states. This procedure discretizes the field into a finite number 

of states. These states, when plotted, form a grid, or mesh, of points. It is at 

these states, or field points, that the solution is found. 

The optimum choice of states, the x,y,z coordinate values, minimizes error and 

computational time. But the process of finding these states is made more difficult 

by complex boundaries, and by the need to control step size differences between 

the states, that is, the need to control the spacing of field points. 

One solution technique uses a different set of state variables, which define a 

different coordinate system, to generate the grid more easily. A new method, developed by Dr. Joseph Steger, combines elliptic and hyperbolic partial differential 

equations into a mapping function between the the physical andr’computational 

coordinate systems. This system of equations offers more control than either equa- 

... 

Vlll 
tion provides alone. The author has modified Steger’s algorithm in order to allow 

bodies with stronger concavities to be used, offering the possibility of generating 

a single grid about multiple bodies. Work has also been done on identifying areas where grid breakdown occurs. This work was supported by NASA under the 

NASA Minority Graduate Researcher’s program. 

-_ ix 
. -. .a . _. . >..:e5 >: ..: "1 .<__ - I. . . L.--L ,I,,-! " >.I >-.":- - llL. &-,I_'. \A -..---I-.*.. 

The governing equations in fluid dynamics, as in many other fields, contain terms 

in partial differentials. In the numerical solution of these equations the differentials are generally replaced by finite difference terms derived from taylor series 

expansions. The solution is then found from given auxiliary data. For instance, 

if u is a function of n variables 

u=u(x1,22, ..., Zn) (1) 

and the equation relating u to these variables is a partial differential equation, 

each m'th partial derivative of u can be expressed as a sum of terms in u, and an 

error term 

where aj and bj are constants and O(Ax)P is the error term. It is important to 

note that, similar to the exact differential on the left, all other x variables are 

held constant while xi varies to produce the different u values. O(Ax)P, the order 

of error, is a measure of difference between the exact differential and the finite 

difference term. The smaller Ax is, the smaller the error. Because Ax is generally 

less than one, higher values of p will also decrease this error. 

A more specific example of the foregoing is a finite difference equation for the 

first partial derivative of u at x; = xi. Dropping notation of x variables held 

constant, equation (2) becomes in this case 

1 
here 

Using equation (2), the solution of u at any x; = z: requires the solution of u 

at neighboring values of zi. The solution for a given set of (z1,x2,. . . x,) values 

is thus dependent on the solution in some neighborhood of these x values. But 

therefore these neighboring solutions are dependent on their neighbors’ solutions. 

For each x; this dependency is generally continuous between one or two z; values 

at which u is definitely known: the auxiliary data. Therefore, the solution process 

starts near values of the x variables for which u is known. One technique, a solution 

by “lines”, is outlined as follows. One of the x variables, say x1 is stepped through 

a set of values. At each value a finite difference equation equal, to some order of 

accuracy, to the original partial differential equation, is set up. The resultant set 

of equations is then solved simultaneously. z2 may then be varied a small amount, 

Ax2, and the process of running through x1 values is repeated to solve another 

set of equations. When 22 reaches its maximum x3 can then be varied a small 

amount, Az3, and the previous processes repeated. 

The fact that u changes only by changing one x variable makes it convenient 

to view the problem within an n-dimensional coordinate system where each 2; is 

represented by an axis. Within this system the set of x values in equation (2) lies 

on a coordinate line’ of the system. In a two-dimensional system, where z = 21 

and y = 22, u, is calculated from values of x and u along a y coordinate line; uy 

‘a coordinate line is described by allowing x, to vary while keeping all other coordinate 

variables constant. a coordinate surface is described by keeping xi constant while varying the 

other variables. 

2 
is found from values of y (and u) along an x coordinate line. figure 1 shows the 

s.~.b,A ._l_ _- u*-w~~ &..L-- . -. . - --- . '- * __- -..-G .a-.,-.c-iibb~er. air3wu ue 

obvious from the figure that the field is now broken up into a set of points. It is 

at these points that the solution is obtained. 

To facilitate numerical computations it is advantageous to keep Az and Ay 

constant, for the following reasons. 

0 The derivation of the above finite difference equation, (3), is based on a 

constant Ax, allowing for a higher order of accuracy and simplicity of form 

than would otherwise be possible. 

0 A constant A leads to the symmetrical array of points in figure 1. These 

points can then be indexed as an array, with points in the x direction being 

indexed by j, and points in the y direction being indexed by k. When one 

boundary value of x , is indexed as x(1), and the other boundary as ~(jmaz), 

the following simple relations can then hold 

~(j) = ~(1) + (j - 1) x AX (5) 

x(jmax) = x(1) + (jmaz - 1) x Ax (6) 

Where j = 1,2,. . . jmax. A similar situation exists for y where the index 

runs from one to some maximum, kmax. One can then say that 

Further, since to each u there associates an x and a y, they can be indexed 

as 
Any point in the field can then be identified by its j,k indices. 

There are three major reasons why the A’s are not usually constant. 

0 To reduce errors due to O(Az)” and increase resolution of u, Az can be 

made smaller. 

0 To increase efficiency Ax can be made larger in areas where u is changing 

more slowly. 

0 Varying y may vary the boundary values of x, so that equation (6) holds 

only if Ax is changed. 

Figure 2 shows the result of varying Ax within the field, in the area labeled 

“increased resolution”. Because of the constancy requirement for a finite difference expression, this new point requires additional points along the x and y 

coordinate lines intersecting this new point. u must be solved for at these new 

points, increasing computational time. 

When the boundary does not lie on a coordinate line, as in figure 3, unless the 

boundary is a linear function of y and x, it is not possible to keep Ax’s or Ay’s 

constant in the field. Although interpolation can be employed or different finite 

difference terms used, accuracy tends to decrease, while computational time tends 

to increase. 

One solution is to define a new coordinate system, a computational system, in 

which 

0 All boundaries coincide with coordinate lines (or surfaces in 3-d). 

4 
.* 

ORlGtJ?<X. FAX 1s 

OF POOR QUALITY 

0 A minimum of coordinate lines are used to to describe the boundq. 

0 If the new coordinate system’s variables are denoted by (,q then u = u((, q), 

-7 --- .. - 

and 

0 the finite difference approximations contain terms in ( and 7, and A( and 

AT0 A( and Aq, the point spacing on the boundary and interior, are constant. 

0 The mapping between coordinate systems must be one-to-one so that 

0 To ensure nonsingularity between corresponding points the Jacobian determinate (J) must be nonzero. 

0 When [,q are the coordinates of a Cartesian system the grid is a rectangular 

mesh of points. 

The requirements lend themselves rather well to using partial differential equations as the mapping functions. Existence of 2nd order partials also ensures 

smoothness and continuity. Although it is possible to use conformal mapping 

techniques to map between two-dimensional coordinate systems, conformal mapping does not extend to three-dimensional systems, whereas partial differential 

equations can map between these systems. 

5 
Figure 1: Portion of boundary points and solution points in a 2-dimensional rectangular coordinate system with Ax and Ay constant. Boundary conforms to a 

coordinate line. 

6 
.. 

Figure 2: Variance in Ax and Ay to increase local solution accuracy and resolution 

with resultant additional number of solution points. 

7 
0 

Figure 1: Consequent change in Az and Ay near boundary that doesn’t conform 

to a coordinate line. 

8 
.. 

2 Classification of Partial Differential Equations 

The general form of a second-order partial differential equation in two variables is 

where u = u(2,y) and a,b,c,d,e,f, and g are functions of x and y for a linear 

partial differential equation. If these coefficients contain first partials of u (u, or 

u,) the equation is termed quasi-linear. 

By an appropriate transformation from 2, y to (, 7 (13) transforms into one 

of three canonical equations, depending on the value of b2 - 4sc. 

For b2 - 4ac > 0 (13) is termed hyperbolic and the canonical equation is 

Ut( - urn + ..’ = 0 (14) 

For b2 - 4ac = 0 (13) is termed parabolic and the canonical equation is 

ut< + = 0 (15) 

For b2 - 4ac < 0 (13) is termed elliptic and the canonical equation is 

Where the - - refers to terms in 2, y, u, ut, and u,,. Each type of equation defines 

what auxiliary data in u, or one of its derivatives, is required for the problem to 

be “well-posed” for solution. This also directly influences what solution procedure 

may be used. If both x and y are spatial functions the auxiliary data is a boundary 

condition. If x or y represents time the auxiliary data represents initial conditions. 

and possibly final conditions. Both hyperbolic and parabolic equations march out 

from the auxiliary data to an unspecified boundary (or final ) condition of u. 

For hyperbolic equations, however, the solution at a point within the domain of 

9 
integration is a function only of a certain portion of the given data. For parabolic 

eqctLkions, ;&e soiution at a point is depeuaent on the entire c~miiary uaiit. iUe 

elliptic equation requires that the auxiliary data in u surround the domain of 

integration. That is, u for the entire boundary of the domain of integration must 

be given. The solution is then found for u within the domain. 

10 
3 Transformation Metrics 

.. i... -_I > . .- I _- -- _I -. I. 3. 

Given a general curvilinear coordinate system denoted by coordinates and 77, 

with a one-to-one point mapping between it and a Cartesian system denoted by x 

and y, that is 

as well as 

then by the chain rule, 

Substituting the values for dx and dy of equations (22) and (23) into equation 

(21) produces 

dt = Sx(qdt + xsdrl) + Sy(Y€dS + Ysdrl). (24) 

Rearanging , 

11 
Since ( is not a function of 77 the terms within the second set of parentheses 

must equal zero and so the terms in the first set of parnetheses must equal 1. The 

resulting set of equations 

produces 

%l 

-277 

tx = - IJI 

IJI ty = 

where 

I JI = qY?l - XvYt 

is called the Jacobian determinant. The same procedure can be applied to 77 and 

the results are 

- 3 

IJI rlY - 

3 

IJI r1x = 

12 
-\ 

4 Elliptic Grid Generation .. -~ 

The simplest, and most commonly used partial differential equation for grid generation is Laplace's equation [3],[2] 

(30) 2 Vu~ux,+Uyy=O 

a homogeneous elliptic equation. 

Replacing u in turn by the transformation coordinate variables ( and 17 produces a system of equations 

v2e=( 5x2 + tyy )=($ 

7x2 + rlyy 

where 

When (,q are considered as the independent variables, equation (31) can be 

written, using the transformation metrics, as 

where 

13 
-- .I -. v - Using centrd--differenc; iorms fcir ;the first. and-seccnk a'eiivatives (see "ai;: 

pendix), noting that A( = AT = 1, and using indexing notation, one finite 

difference form of equation (32) is 

where 

After clearing the denominators and expanding, equation (33) becomes 

Gathering all 2 terms in k to the left side, this now becomes 

14 
-- 

I” 

where 

This system can be solved by the “line” solution technique mentioned earlier: 

by keeping k constant (starting at k=l) an equation for each Rj,l (from j=l to -. 

j=jmax) is generated. The resultant set of equations is then solved simultaneously 

for i. The index variable k is then incremented and the process, of generating a 

set of equations for all j, is repeated. The overall process of incrementing k and 

solving a set of equations continues until a final k=kmax value is reached. Because 

the values of all x and y are initially unknown, except at boundary points, initial 

guesses for them must be made. These initial x and y values are then used to 

generate the values for A,B, and G. As new vdues of ij,k ( 2j,k and Yj,k ) are 

found, they are used to solve the succeeding set of equations for the next set of 3 

values. 

equation (38) is a system of 2 equations, one for each unknown coordinate 

variable x, and y. Substituting x for each instance of the vector R’ the equation 

becomes 

The complete set of equations, for constant k, from j=l to jmax, can be put 

in matrix form as 

15 
, -- 

-. 

B2.k 

32 jmoz- 1 ,k 

3 0 Ajmoz-2,k Bjmaz-l,k 

where each Gzj,k equals the right-hand-side of equation (39). A similar set of 

equations are solved for y at each k level. 

Having run through all j and k values an improved set of x,y values is obtained. 

This improved set can then be used to obtain an even better set of values by 

repeating the procedure just described. This iterative process stops when at some 

iteration the change in x and y from a previous iteration is smaller than some 

given amount. 

The number of iterations can be reduced significantly by using a technique 

known as successive line over-relaxation, [6] which is essentially a forecasting technique. The latest values of x and y are compared with their previous values and 

some fraction of their difference is added to the latest value, to predict a better 

value for the z and y’s. For example, after solving the set of equations along a 

line using Thomas’ algorithm, and obtaining new values for x(j,k) at each k level, 

where n refers to the iteration number, the following modification can be applied. 

where n references the last iteration of x and n-1 references the previous value. 

w controls the amount of change made. When w lies between one and two the 

technique is known as over-relaxation. When w is between zero and one the 

technique is known as under-relaxation. Under-relaxation is useful when values 

16 
lof the dependent variable oscillate for a given set of values of the independent 

variables, that is, oscillation with respect to iteration 

The main advantage of Laplace's equation in grid generation is that it serves 

to smooth grid lines. This is because it contains within it the maximum principle: 

there are no maximums or minimums of u within the domain of integration. That 

is, u attains its highest or lowest values on the boundary of the domain, u(zo, go) is 

an average of u in the neighborhood of zo and yo. This serves to reduce disparities 

in u for neighboring points. This can be shown by noting first that, within the 

domain of integration, the dependent variable cannot have a maximum or minimum. For there to be an extremum the first partial derivatives of the dependent 

variable must each go to zero; the second partials must have the same sign and 

each must be nonzero. This would mean u,, + uw # 0 which is not the case here. 

Since u does not have an extremum within the boundary it must be bracketed by 

values of u on the boundary. 

Figure 4 is an initial guess for a grid about a body. This grid was generated 

by placing points evenly on a circular outer boundary and joining them to points 

on the inner boundary, the body. Then points were evenly spaced along these 

connecting lines to generate the nodes. Figure 5 shows the results of 100 iterations 

of the elliptic solver. There has been a general smoothing of lines throughout. 

Figure 6 is a detail of the initial grid in figure 4 showing grid breakdown near one 

of the convex points. Figure 7 is this same area after 100 iterations of the elliptic 

generator 

The main disadvantage to using this system of equations is that there is no 

control over grid point spacing. For each given set of boundary values there is 

17 
.- 

only one find distribution of interior points. An alternative is to use the nonhomogeneous form of the elliptic equation, known as Poisson’s equation 

V2[ = p 

v2q = q 

In such equations p and q can be used to modify the strong tendency of the 

Laplacian to equalize all cell volumes [7]. Another alternative takes the fully 

converged grid and re-spaces the points along ( lines in order to cluster points 

near the body. [4] 

18 
Initial Grid 

Figure 4: tau=.5, cam=0.8, Initial Grid 

file: s5t8.00~~ 

19 
.- 

The Effect of the Elliptic Solver 

Figure 5: tau=.5, cam=0.8, pass 100 

file: s5t8.10ps 

20 
Initial Grid 

Figure 6: tau=.5, cam=0.8, detail of initial grid of figure 4 

file: s5 t8.00ptps 

21 
.- 

The Effect of the Elliptic Solver 

Figure 7: tau=.5, cam=0.8, pass 100, detail of final grid of figure 5 

file: s5t8. lOptps . 22 
.- 

-- 

I 5 Hyperbolic Grid Generation 

Steger [8] developed a grid generation scheme based on a set of equations which 

together form a hyperbolic system. These equations are expressions of constraints 

made on the grid generating system. These constraints are as follows: 

0 Orthogonality is explicitly specified, 

0 Grid cell volume (i.e. area in two dimensions), V, is user defined. Since, in 

the numerical implementation A( = Aq = 1, 

Thus the Jacobian should be nonsingular and approximately equal to the 

physical cell volume, AzAy. 

The resultant system is thus represented by 

Jxvx +tyvy = 0 

txvy -tyvx = J 

or, after a change of variables, again using the transformation metrics, 
.- 

I --. 

Equations (49) and (50) comprise a system of nonlinear partial differential 

equations with initial data in 7. Through a process known as local linearization, 

shown below, it is found that this system is hyperbolic and thus can be marched 

in 7 from initial data along the body, 7 = 0. 

If the substitutions z = xo + i and y = yo + @ are made, where xo and yo 

represent a known nearby state, and L and @ are of the order of Ax and Ay, 

local linearization can be performed by making these substitutions for x and y in 

equations (49) and (50). Starting with the first equation, 

Similarly, the remaining terms in the two equations become 

Using these substitutions equation (49) becomes 
or, in matrix notation 

Doing the same for equation (50) 

(Y;q + +r/- - - (y,"., + + - ZiY!) = v (61) 

Since by equation (50), (zpy," - z:y,") = Vo the equation becomes 

+ $Yv - yp, 0 - Z9y( 0 = v + vo 

or 

(4 -4( ;) +( -Y; x;) (;) =v+vo (63) 

€ 9 

Combining equations (60) and (63) 

which can be represented as 

A& + BIT,, = f. 

25 
-- 

If B, whose determinant is the Jacobian transformation, is nonzero, all cell 

volumes are non-zero, and its inverse, E', exists. Equation (64) can be then be 

transformed into 

It has been shown [l] that this type of equation is hyperbolic and marches in 

the 17 direction if B-lA has real, distinct eigenvalues, or if B-lA is a symmetric 

matrix. Expanding B-lA and performing the indicated multiplication produces 

which is symmetric, fulfilling the latter requirement. If A represents the eigenvalues the eigenvalue equation is 

or 

Thus, since Vo, the inverse Jacobian, is non-zero the eigenvalues are real, and 

distinct when (z!z: t $#) is non-zero 

26 
5.1 Numerical technique 

The numerical solution of equation (66) can be carried out using the line inversion 

mentioned earlier in the introduction and in the section on the elliptic generator. 

The main difference here is that only one pass is performed in generating the 

grid. Unlike the elliptic generating system, which improved an initial grid, and 

therefore was dependent on the quality of that grid for convergence to a solution, 

the hyperbolic system generates a grid from the known boundary conditions. 

Replacing the partial derivatives in equation (65) by finite difference operators 

produces 

where all terms are second order accurate except for the first order replacement 

of &. 

Since, in the numerical implementation, A( = A7 = 1, the expansion of 

equation (71) is 

Rearranging to gather terms involving IC + 1 on one side produces 

-+ -6 

a j,k dj- 1 ,k+ 1 + bj,k Rj,k+ lcj,kRj+ 1 ,k + 1 = RHS (73) 

where 

B-'A 

aj,$ = -- 2 (74) 

27 
bj,k = 1 

B-~A 

2 

- cj,k - 

RHS = B-lf +A!j,k 

The coefficients are all evaluated at the same index values of j and k, where 

x and y are known. The partials with respect to ( are straightforward central 

differences. The 7 partials can then be obtained by solving for them in equations 

(49) and (50) 

xc v 

(78) 

(79) 

Equation (73) generates the set of equations in j, from j=l to j=jmax. These 

equations are solved, as before, using Thomas' algorithm. 

The results of using this algorithm, designated hyg2d, can be seen in figures 8 

through 11 for bodies of increasing concavity2. Note that as the bodies increase 

in concavity the radial grid lines start to cluster, and finally, the grid breaks 

down. figure 12 and 13 show the marching out of the grid just before and at this 

breakdown. Figure 14 is a detail view of one point of the body. 

lSee the section Notations and Definitions and figures 69-72 for more information on concavity 

and how it is meant here. 

28 
The Effect of Camber (concavity) on the Hyperbolic Solver 

Figure 8: tau=.S, cam= 0 

file: hSOy0.ps 

29 
The Effect of Camber (concavity) on the Hyperbolic Solver 

Figure 9: tau=.5, cam= 0.40 

file: h50y40.p~ 

30 
The Effect of Camber (concavity) on the Hyperbolic Solver 

Figure 10: tau=.5, cam= 0.80 

file: h50y80.ps 

31 
The Effect of Camber (concavity) on the Hyperbolic Solver 

Figure 11 : tau=.5, cam= 1 .O 

file: h5ylO.p~ 

32 
... 

Marching Out of Hyperbolic Generator 

Figure 12: tau=.5, cam= 1 .O, just before cavi ty-caused breakdown 

file: h5y 1O.k 15ps 

33 
-- 

Marching Out of Hyperbolic Generator 

Figure 13: tau=.5, cam=l.O, breakdown due to strong concavity 

file: hSylO.kl6ps 

34 
The Effect of Camber (concavity) on the Hyperbolic Solver 

Figure 14: tau=.5, cam= 1.0 detail 

file: h5ylO.ptps 

35 
6 Combining Elliptic and Hyperbolic Equations 

Each of the grid generation schemes discussed, the elliptic and the hyperbolic, 

has a profile of advantages and disadvantages. The elliptic, through its smoothing 

action is able to generate a grid about more complex bodies than the hyperbolic. 

But it requires an initial grid to start from; it must iterate to convergence for 

a solution; Unless source terms are added, there is little control over point distribution; these source terms are not easy to use. The hyperbolic, on the other 

hand, creates a completed grid in the order of time it takes the elliptic to make 

one pass. This particular hyperbolic scheme has, additionally, orthogonality constraints, and user-control of grid cell volumes, providing a strong control over 

point distribution. Its chief disadvantage is that unlike the elliptic it breaks down 

sooner in areas of strong concavities or convexities. 

A scheme that uses the best of both would be better than either. Steger [5] 

has developed such a grid generation scheme. His procedure is to add a scalar 

multiple of the elliptic system, equation (32), to the right side of the hyperbolic 

system, equation (66), producing 

The control parameter p can then be used to emphasize the hyperbolic or elliptic system. For small p the hyperbolic dominates and for large p the elliptic 

dominates. 

36 
.- 

Viewed as a hyperbolic system the constraint of orthogonality is now relaxed 

by the strength of the elliptic terms when they sum to nonzero. Loosening this 

constraint helps avoid grid breakdown. In the volume equation the elliptic terms 

serve to smooth differences between cell volumes. Viewed as an elliptic system, the 

Laplace equation has been transformed into a Poisson equation with the hyperbolic terms serving as the source terms. These source terms allow concentration 

of 6, and 7 lines. Which side dominates is controlled by p. When it equals zero 

the system is purely hyperbolic. As it increases fiom zero the right side of the 

equations become more dominant until p is over about 3,000, when the system is 

effectively elliptic. This parameter can be made a function of spatial position, in 

order to, for instance cluster lines near the body. It can also be made a function 

of gradients in the flow field so that, in an area of high gradients, p can be set 

low, turning on the hyperbolic’s volume control, to increase point density in this 

area. In areas of grid breakdown p can be set high, turning on the elliptic, to 

repair the breakdown. 

6.1 Solution Procedure 

Replacing the partial derivatives in equation (6) by finite difference operators 

produces 

where all terms are second order accurate except for the first order replacement 

of &. 

37 
-- 

Since, in the numerical implementation, A( = Aq = 1, the expansion of 

equation (81) is 

Rearranging to gather terms involving k on one side produces 

where 

and all values in equations (84) - (88) not specifically indexed have index values 

of j and k, that is, they are evaluated at the points x = xj,k, y = yj,k. 

Equation (83) can now be solved similar to the elliptic equations. Starting with 

k=l, letting the j index run from 1 to some jmax, and keeping the k index constant, 

a set of equations are generated. The left-hand coefficients form a tridiagonal 

matrix, each element of which is a two-by-two array. The set of equations is 

38 
-- 

solved using a modified form of Thomas’ algorithm. The index variable k runs 

from 1 to some number kmax and therefore kmax sets of equations must be solved. 

Once these are solved the whole process is reiterated using the results of the last 

iteration, the resultant grid, as the source for generating new values for the righthand side and the coefficients, of equation (83). The iteration process halts when 

grid changes are below some threshold value. 

6.2 Testing Procedure 

Since the elliptic can repair breakdowns in the grid (due to the hyperbolic), and 

since breakdown is a function of the number and strength of convexities and 

concavities, it was decided to test this new solver by using bodies of increasing 

concavity and/or convexity. 

Although the hyperbolic generator breaks down in areas of either great concavity or great convexity, breakdowns due to the latter tend to be more localized 

and simpler to deal with. Uneven point spacing along the body on either side of 

the apex of the convexity causes the grid line radiating from the apex to swing 

over. This causes circumferential lines to cut through the body. But the number 

of nodes involved are generally limited to a few in either the circumferential (6) or 

radial (17) direction and can be handled by the fortified approach of Van Dalsem 

[91. 

The breakdowns caused by a concavity tend to be more global as the strength 

of the concavity, and the percent body surface the concavity takes up, increase. 

Therefore, it was decided to approach the problem by examining bodies containing 

concavities of increasing strength. 

39 
-. 

-- 

I 6.3 The Unmodified 

Implementation of the elliptic-like algorithm just described requires an initial grid. 

The main requirement placed on this grid is that its outer boundary not selfintersect or intersect with interior grid lines. The hyperbolic generator was used 

to generate the initial grid of figure 15. This generator offers the advantage of being 

closely related to the hyperbolic part of the combined generator. The similarity 

can be used to avoid problems such as “collapsing” grid lines, q lines that tend to 

fall back towards the body rather than away from it with succeeding iterations. 

This defect can occur when the initial grid has different spacing in the 7 direction 

than the hyperbolic component provides. Here p is set low near the body and 

increases geometrically in the radial direction. Thus the hyperbolic is enforced 

close to the body allowing clustering of points, while the elliptic is more in force 

deeper into the field out to the outer boundary, smoothing the sharper curvature 

of grid lines of the initial grids. Figure 16 shows the grid after 100 iterations. 

Note that because the outer boundary the clustering on it cannot be dealt with. 

An alternate initial grid, having a strong concavity and the final grid drawn 

from it are shown in figures 17 and 18. Figure 19 shows detail of the latter 

grid. at a convex point. This initial grid was generated by specifying a circular 

outer boundary and spacing the points on it evenly. Lines are drawn connecting 

these points with points on the body. Points are then evenly spaced along these 

connecting lines to form the grid. Here the outer boundary point spacing was 

fixed to avoid the clustering of the hyperbolic system. For high concavities it 

seems to generate a reasonable grid. Because the outer boundary point spacing 

is not a strong function of the inner boundary, their connecting lines (t lines) 

40 
-- 

sometimes have to make large sweeps. As shown in figure 19 orthogonality can 

break down. Other drawbacks are that there is no way of knowing which points 

are best connected to each other and, eventually a body shape is reached at which 

this algorithm breaks down. 

6.4 Floating the Boundary 

The major problem with the original algorithm is that it found a solution in 

an elliptic way. That is, the boundaries were kept fixed while the inner field 

points were adjusted. The first modification to the algorithm above, and used in 

subsequent modifications, was to float the boundary. This was achieved by turning 

the hyperbolic on at the outer boundary to generate a new outer boundary. This 

technique is more amenable in areas where there is no physical outer boundary, 

such as for airfoils, areas that use hyperbolic generators to form grids. 

Once the outer boundary is floated initial grids containing defects in their 

outer boundary, such as those generated by hyg2d, can be used. Such defects 

include boundaries that intersect other parts of the grid. The smoothing action 

of the elliptic, coupled with the boundary generating action of the hyperbolic can 

then repair defects in the initial grid. Figures 20 to 33 are a series of grids showing 

the effect of floating the boundary. 

The drawback to this latest technique is that at some point a body is reached 

whose concavity makes it too difficult to generate a grid about. The generator 

then fails to improve the grid. Even before this point is reached, though, the range 

of values p can take on narrows, making it difficult to choose a set of values for it. 

Here p is just a function of 7. For bodies with weak concavities it can take on any 

value at any 7. But as concavity grows grid improvement becomes very sensitive 

41 
to the due of p, to the point that it becomes very difficult to find dues for p 

that will enable a solution to converge. 

6.5 Using Compressed-Cell Initial Grids 

The problem with the previous technique was in the initial grids used. They 

contain too many defects and slow the convergence procedure. If cell volumes are 

extremely small the distance from inner to outer boundary makes intersections 

of grid lines more difficult. This led to modifying the algorithm to handle a 

different kind of initial grid, one with fewer defects. An initial grid of compacted 

cell volumes was used. During each iteration all cell volumes were then increased 

a fixed percent of their total volume. This had the benefit of starting with an 

initial grid which had very few, if any, defects and smoothing out defects as they 

developed, rather than starting with a grid containing many defects. Defects were 

generally isolated then, at or near the outer boundary. Each iteration was both 

a growth of the grid, somewhat like the hyperbolic, and also each iteration was a 

smoothing out in the interior, by the elliptic. Since the outer boundary is purely a 

product of the next lowest level, which was constantly being improved, this made 

for a very robust generator. Compared to the previous method this technique 

allows for a wider range of values for p given the same body as before. Figures 

34 to 45 are examples of this technique for bodies of high concavity. Cell volumes 

have been set to grow about fifteen percent per iteration. Figure 34 contains the 

initial grid but, since its cell volumes are compressed the cells can only be seen 

at at several magnifications (figure 35). Figures 36 to 40 show the grid expanding 

during successive passes of the combined solver. Figure 41 is a close-up of the the 

grid at the 200th pass. Figures 42 to 45 shows that the technique can be applied 

42 
to even deeper cavities. 

6.6 Curtailing Outer Boundary Movement within a Cavity 

There are two main drawbacks to using the technique just discussed as it stands. 

0 A point is reached in which the cavity is too deep for this technique to work. 

see figure 46 

.. a Although a grid with no breakdown is possible for strong conc iL,A.L.25, -I--. ":'e:? 

tends to be less control over point distribution within these areas. 

The technique in the previous section uses just the hyperbolic component to 

generate a new outer boundary with each iteration. The hyperbolic generates 

the points on this boundary, the line from data on the next lower section 

of the grid, the 7kmaz-l line. By the orthogonality constraint these points are 

normal to the qkmaz-1 line. As the cell volumes increase this causes boundary 

self-intersection within a cavity with some depth. Because the interior is derived 

primarily elliptically and therefore a function of both boundaries there is a tendency therefore for grid breakdown near the outer boundary. But, since the inner 

boundary is okay, there is, by the maximum principle less breakdown of these 

interior grid lines than the boundary line. At each iteration these areas of breakdown are improved by the elliptic through communication with inner grid points. 

These improved areas are then used to generate the latest outer boundary. The 

fact that the numerical technique starts at the inner boundary may be part of the 

reason why the elliptic smoothes these areas so well. 

43 
-- 

The result of the above tendencies is that within the cavity the outer boundary 

moves, at each iteration, largely in the 77 direction and out of the cavity. Thus in 

strong concavities the ratio of cell side lengths can be very high. See figure? And, 

while volume control normally gives control over point distribution, in this case it 

is overidden by the smoothing action of the elliptic coupled with the hyperbolic’s 

regeneration of the outer boundary, forcing the boundary out of the cavity. Within 

the cavity accuracy and resolution are thus lower in the 77 direction because of the 

resultant spacing. It was thoUg&L&&s& &a&: pbinkdensity wa.t+.lowesC inithe 

cavity, contributing to the increased cell aspect ratio, but analysis confirmed that 

for a given range of 7 point density is about the same within and without the 

cavity. This conforms with the fact that cell volumes had been functions of 77 

only. 

In and near an area of a concavity the gradients of the flow variables can be 

expected to be higher than average. Therefore A, and A, should actually be 

smaller than elsewhere to resolve the flow and provide a stable solution. There 

then appears to be a conflict between the need to maintain a small cell size and 

at the same time the need to cover the flow field out to several body diameters. 

It was therefore decided to modify how the outer boundary is generated within 

the cavity. A point on the boundary is chosen to serve as a hinge point. As the 

boundary points move away normal to the ~k,,,==-~ line their movement is modified 

to some specified degree by a forcing function that brings corresponding points, 

on each side of the hinge point, together. When a set of points meet they are no 

longer allowed to float. Their corresponding cell volumes are then held constant 

through later iterations. Figures 47 to 62 show the results of using this technique 

44 
-- 

on the grid that caused breakdown of the previous technique. There is The result 

is that within the cavity the boundary overlaps itself along a line of symmetry 

of the cavity and it is now possible to generate a grid for deeper cavities also a 

higher point density within the grid. 

45 
-- 

Using the Combined Solver with a Fixed Outer Boundary 

Figure 15: tau=.5, cam= 0.8, initial grid 

file: eh5y8.00~~ 

46 
Using the Combined Solver with a Fixed Outer Boundary 

Figure 16: tau=.5, cam= 0.8, pass number 100 

file: eh5y8. lops 

47 
.- 

The Effect of Camber (Concavity) on the Combined 

Elliptic-Hyperbolic Grid Generator 

Figure 17: tau=.S, cam=2.0, Initial Grid 

file : s5t20.pOpstmp 

48 
.- 

The Effect of Camber (Concavity) on the Combined 

Elliptic-Hyperbolic Grid Generator 

Figure 18: tau=.5, cam=2.0, pass 100 

file: sehSt20.1 Ops 

49 
.- 

Figure 19: detail of figure 18 

file :seh5t20.1 Optps 

50 
Grid Improvement by Floating the outer Boundarv 

Figure 20: tau=.5, cam=l.O, initial grid 

fie: ehf510.00ps 

51 
-- 

Grid Impmvement by Floating the outer Boundary 

Figure 21: tau=.5, cam=l.O, pass number 1 

file: ehf510.01ps 

52 
Grid Improvement by Floating the outer Boundary 

Figure 22: tau=.5, cam=l.O, pass number 2 

file: ehf5 10.02ps 

53 
.- 

Grid Improvement by Floating the outer Boundary 

Figure 23: tau=.5, cam=l .O, pass number 3 

fie: ehf5 10.03~s 

54 
-. 

Grid Improvement by Floating the outer Boundaro 

Figure 24: tau=.5, cam=l.O, pass number 4 

file: ehf5 10.04~s 

55 
Grid Improvement by Fhting the outer Bwndary 

Figure 25: tau=.5, cam=l .O, pass number 5 

fde: ehf5 10.05~s 

56 
Figure 26: tau=& cam=l .O, pass number 6 

file: ehf5 10.06~s 

57 
Grid Improvement by Floa€hg the outer Bmndarv 

Figure 27: tau=.5, cam=l .O, pass number 7 

file: ehf510.07~~ 

58 
Figure 28: tau=.5, cam=l .O, pass number 8 

file: ehf510.08~~ 

59 
.- 

&id LFnprovemen$ by Floslting the outer Bsundary 

Figure 29: tau=.5, cm=l .O, pass number 9 

file: ehf5 10.09~~ 

60 
Grid Improvement by Floating the outer Boundam 

Figure 30: tau=.5, cam=l.O, pass number 10 

file: ehf5 10.1 ps 

61 
.. 

Figure 31: tau=.5, cam=l.O, pass number 20 

file: ehf5 10.2~s 

62 
Grid Improvement by Floating the outer Boundary 

Figure 32: tau=.$ cam=l.O, pass number 50 

file: ehf510.5ps 

63 
Grid Improvement by Floating the outer Boundary 

Figure 33: tau=.5, cam=l.O, pass number 100 

file: ehf5 10.1 Ops 

64 
.- 

..- 

Initial Grid with Compressed Cells 

Figure 34: tau=.5, cam=2.0, full figure 

file: hg5y20.ps, initial dseta=le-5, final dseta=le-3 

65 
Figure 35: tau=.5, cam=2.0, detail of initial grid near convex point 

file: hg5y20.ptps initial dseta=le-5 final dseta=le-3 
1. 

Y 

T 

Figure 36: tau=.5, cam=2.0, 40th pass 

file: fg520.4elxps 

61 
. , . . . .... I 'C. 

- . . , .,,. , ._, . . :.,.i e5 

6~ poow QUALITY, 

Y 

t 

Figure 37: tau=.5, cam=2.0, 50th pass 

file: hg520.5elxps 

68 
-- 

Figure 38: tau=.5, cam=2.0, 80th pass 

file: fgS20.8elxps 

69 
-- 

Figure 39: tau=.5, cam=2.0, 100th pass 

file: fg520.1Oelxps 

70 
-- 

Cell Volume Growth Coupled with floating Boundary 

Figure 40: tau=.5, cam=2.0, 200th pass 

file: fg520.2Oelxps 

71 
.- 

I-- 

,. 

Cell Volume Growth Couded with floatin(= Bolandaw 

Figure 41: tau=.5, cam=2.0, 200th pass, detail near convex area 

file: fg520.2Oelxpts 

72 
-- 

IY 

t 

Figure 42: tau=.5, cam=4.0, Full Figure 

file: hg5y40.p~ initial dseta=le-5 final dseta=le-3 

73 
Cell V 0119pvm- taatjarg Boundary 

T 

Figure 43: tau=.5, cam=4.0, 50th pass 

file: fg540.5~~ 

74 
.. 

Cell Volume Growth Coupled with Floating Boundary 

Figure 44: tau=.5, cam=4.0, 100th pass 

file: fg540.10~~ 

75 
.- 

Cell Volume Growth Coupled with Floating Boundary 

Figure 45: tau=.S, cam=4.0, 200th pass 

file: fg540.20~~ 

76 
Ce€l Vekmre Gpowth Coupled wkh Floating Boundary 

Failure of Technique Due to Depth of Cavity 

T 

Figure 46: tau=.5, cam=6.0, 500th pass 

file: fg560.50~~ 

77 
~ - -- 

Forming a Cut Ina-Cavity-, -- . -. 

Initial Grid 

Figure 47: tau=& cam=6.0, initial grid 

file: fgc56ops 

78 
.- ORIGINAL PZ?"Z ;s 

OF POOR QUALITY 

Fomi??z a Cut In a cavity 

Figure 48: tau=.5, cam=6.0, initial grid, detail of cavity 

file: fgc560.00ptps 

I 

79 
Y 

t ‘1 

I 

Figure 49: tau=.5, cam=6.0, pass 20 

file: fgc560.2~~ 

80 
.. 

Figure 50: tau=.$ cam=6.0, pass 20 , detail of cavity 

file: fgc560.2ptps 

81 
.. ORIGINAL PXX IS 

OF POOR QUALITY 

Y 

t 

Figure 51: tau=.5, cam=6.0, pass 30 

file: fgc560.3~~ 

I 

82 
Forming a Cut In a Cavity 

Figure 52: tau=.5, cam=6.0, pass 30 , detail of cavity 

file: fgc560.3ptps 

83 
Forming a Cut In a Cavity 

Y 

1’ 

Figure 53: tau=.5, cam=6.0, pass 40 

file: fgc560.4~~ 

84 
-. 

-- 

Forming a Cut In a Cavity 

Figure 54 tau=.5, cam=6.0, pass 40 , detail of cavity 

fi 1 e: fgc560.4~ tps 

85 
T 

Forming a Cut In a Cavity 

*X 

Figure 55: tau=.5, cam=6.0, pass 50 

file: fgc560.5~~ 

86 
-- 

Forming a Cut In a Cavity 

, 

Figure 56: tau=.5, cam=6.0, pass 50 , detail of cavity 

file: fgc560.5ptps 

87 
-- 

Fixming a Gut In a cavity 

-- 

Figure 57: tau=.5, cam=6.0, pass 80 

file: fgc560.8~~ 

88 
-. 

Formin? a Cut In a Cavity 

Figure 58: tau=.5, cam=6.0, pass 80 , detail of cavity 

file: fgc560.8ptps 

89 
-. 

Figure 59: tau=.5, cam=6.0, pass 100 

file: fgc560.1 Ops 

90 
-- 

Figure 60: tau=.5, cam=6.0, pass 100 , detail of cavity 

file: fgc560.1 Optps 

91 
Figure 61: tau=.5, cam=6.0, pass 300 

file: fgc560.30~~ 

92 
Figure 62: tau=.5, cam=6.0, pass 300 , detail of cavity 

file: fgc560.30ptps 

93 
Since it is required that the mapping between real space and computational space 

be one-to-one, that is 

R' = ii(C) 

and 

c' = c'(R') 

The intersection of grid line segments, other than at nodes, would violate this 

requirement for the following reason. The intersection of these grid lines is an 

intersection of lines connecting points which have adjoining index values. For 

instance one line segment may join P(j,k) with P(j+l,k) the other may join P(i,l) 

with P(i,l+l). Since t and q are linear functions of these indices and the indices 

are different in each set, the range of values of t and q are different between one 

set of points and the other set. Thus at the point of intersetion in the r space 

coordinate system two values of 6 exist for the one value of i, violating the 

one-to-one requirement. 

Since the only grid line that can be assumed satisfactory is a given boundary 

line, a check for line intersection should most logically start there and work towards 

the other boundary. But a more elementary check should be to first check whether 

the grid points nearest the boundary, and therefore the lines connecting them to 

the body, are outside the body, as they should be. If k indexes points in the radial 

or q direction, and j indexes points in the circumferential or ( direction, then these 

first interior grid points can be noted as P(j,Ic) = P(j,2) G joins the 

I 

94 
body at the point Pj,l which is between the points Pj*-,1. 

Figure 63: Portion of boundary and first row of grid cells next to boundary. Here 

k=2 

If one defines following vectors 

then, for j increasing clockwise around the the body and k increasing radially 

away, it can be shown that if either of the following sets of conditions are met I 

and 

r'tx?2>0 

or 

forF1 x F2 < 0 

either 

?lXri>O 

95 
-- 

or 

GXF2>0 (954 

then 6 and Pj,2 is on the wrong side of the boundary, the previous 77 line. If 

all the points at the k=2 line pass this test, the second test should be run before 

checking points at the next k level. This first test compares each point in the 

grid with points having neighboring j and k values. Thus, it uses relatively little 

computational time per point, but each point must be checked. 

If the first test checks out the second test should be to check whether the k 

line intersects a previous k line. Figure 64 shows an example of this defect and the 

main reason to check for it. For strong convexities, when distances between the 

apex and its two adjoining points on the boundary are significantly different, the 

hyperbolic generator can cause the radial line emerging from the apex, the line 

(the radial line near the bottom in figure 64), to swing far enough over such that 

intersections occur. This can be checked by testing the 7 line segment at each 

point for intersection with the next lower 77 line. This check, like the previous 

one, though it checks each point in the grid, makes only one or two checks for 

intersection per point. 

A third breakdown, differing from the previous two in that there is no relation 

between the j values of the points compared, occurs when different parts of a 

boundary (the kmax line), or boundaries of different grids, intersect. The simplest 

procedure is to test each boundary line segment for intersection with all other 

boundary line segments. An improvement on this procedure is the following. Since 

intersection occurs only for lines whose x and y values overlap for some portion 

of the lines, computational time can be reduced, if the the check for intersection 

96 
.- 

is made only after the simpler checks for x and y overlap are made. This can be 

done by first sorting either x or y into a'monotonic sequence and examining the 

new arrangement of the old j values. When successive j values differ by more than 

one, indicating overlap, the other overlap check can be made. In other words, if 

the original sequence is x(j), the new sequence can be denoted as x(j(i)). Then if 

j(i) and j(i+l) differ by more than one, there is overlap of x values of two separate 

line segments and the next test, to check for overlap of y values is made. If there 

is y overlap the check for intersection can be made. Because only points on the 

outer boundary are checked, the check for boundary intersection takes about the 

same order of magnitude of time as the previous two procedures the check for 

breakdown is of the same order of time as the previous two. 

97 
Grid Breakdown Due to Strong Convexity 

and Unequal Spacing Along Body 

Figure 64: tau=.5, cam=2.0, detail of hg5y20.a~ used as initial grid 

file : hg5y20.ptdps, bold line indicates body/boundary 

98 
8 Derivation of the Navier-Stokes Equations 

The governing equations in fluid dynamics are based on the three conservation 

laws of mass, momentum, and energy. These equations can be derived in the 

following way [lo], [6] using a control-volume of fixed dimensions within which 

the three laws are analyzed. 

AX 

Figure 65: Fixed control volume showing area vector A and the three surface 

stresses at one of the two surfaces normal to the x axis. 

For the mass, or continuity, equation the change in mass per unit time, dmpt 

within the control volume (vol), is balanced by the net sum of mass entering and 

leaving the volume in unit time. Since density, p, equals dm/dvoZ, the rate of 

change in mass within the volume is 
(97) 

d $ p dvol 

at - 

, 

= jzdvol. 

.. IThe mass entering the volume per unit time can be represented as 

The minus sign is required because Ais directed away from the body. Combining equations (98) and (99) produces 

JZdvol = - iurface p U' ''0 

Using Gauss' theorem, J 8 di = J' $ ddvol, equation (100) becomes 

J -duo1 = - J f (p Z)dvol, at 

which leads to the point form, since integration is the same for all terms, 

-+V.(pii)=O + -# 

at 

The first term on the left exists when mass density within the control volume 

is unsteady. The second term indicates the rate of mass per unit volume building 

up within the control volume due to differences in density and/or flow rates at 

the various surfaces. When the first term is zero, indicating no change in mass 

density at a point, the second term must also be zero. 

The momentum equation is a vector equation, containing three scalar equations, one for each direction. It is derived from Newton's second law: the rate of 

change in momentum for a body is equal to the forces acting on the body. For 

the ked control volume, with mass flowing through it, the rate of change of momentum (of the mass within it), is equal to the net rate momentum is convected 

100 
I- across its entire surface by the flow, plus the sum of the forces acting on the mass 

within it 

where P; is the momentum within the control volume, P, is the momentum 

convected across the surface, Fb represents body forces, such as gravity, and F, 

represents surface forces such as pressure and friction. 

-L -b 

-+ 

Since dP = .ii dm, because ii may vary within the control volume, and dm = 

pdvol, the rate of change of momentum within the volume is 

The momentum convected across the volume surface, per unit time, (covering 

a distance dr) is found by multiplying the differential of momentum at the surface with the dot product of the flow velocity ii and the surface element i, and 

integrating over the whole surface 

dt dt 

dt 

The last change is by Gauss' theorem. As before, the negative sign in front of 

the integral is due to the direction of the area element. 

101 
When these momentum terms are on one side of the equation they then equal 

the sum of the forces acting on the mass in the control volume: body forces, such 

as gravity, and surface forces, or stresses, such as pressure and friction. 

For the control volume the body forces are integrated over the entire volume. 

If & is the total of body forces on the mass m in the volume then & = @b/m is 

the force per unit mass. And since density is defined as p dm/dvoZ, p& is the 

force per unit volume. Then the total force due to body forces is 

(110) 

Surfaces forces can be divided into two distint types: normal stresses and 

tangential stresses. Normal stresses, denoted as r3i = cii, are pressure forces and 

act normal to the surface in question. The subscript indicates both the surface 

being acted upon and the direction of the force, which is normal to the surface. 

So Zz acts on a surface area dA lying on a y-z plane (the line normal to a y-z 

plane lies in the x direction). 

Tangential forces are shearing, or friction, stresses, such as those associated 

with viscosity. At a surface there are two independent tangential directions, therefore an indication of this must be made in denoting the symbol for this stress: 

7f'k = r;& The first subscript, like that for the normal stress, denotes the plane 

the surface lies on. The second subscript indicates the direction of the force. So 

the two tangential forces acting at the y-z surface dAi are Fsy = rzyi acting in 

the y direction, and Fzz = r,,i,acting in the z direction. 

The normal and shearing stresses form a stress tensor which is symmetric 

I 

102 
ORIGINAL PAGE I$ 

OF POOR QUALITY 

The surfaces forces are integrated over the entire surface. Since the surface 

stresses form a tensor, its transformation by Gauss' theorem produces the following equation for the sum of the surfaces forces ga acting on the control volume 

where 

-* dSi j v.s=- i, j = 1,2,3. dx j 

The momentum equation is then, after removing the integration signs, 

and the x component is 

Since u is the negative of pressure this equation becomes 

103 
The equations for the y and z components of momentum are found in a similar 

way. ox Tz7J rxz \ 

The energy equation derives from the first law of thermodynamics 

dE = dQ - dW (119) 

In the case of a control volume with mass passing through it, the rate E;, the 

energy within the control volume, is changing, equals the rate energy, E,, is being 

convected in by the flow, plus the rate heat, Qconducted, is conducted in (ignoring 

heat radiated in for flows of sufficient speed), plus the rate work Wa is being done 

on the mass within the volume due to body forces, plus the rate work, Wdp, is 

being done on the mass due to surface pressure forces, plus the rate work, W,,, is 

being done on the mass due to surface viscous forces 

The energy per unit mass is p(e + 1/2u2 + ep), with e, 1/2u2 and ep representing 

internal, kinetic, and potential energy respectively; u2 = ii - ii. Similar to the 

previous equations for mass and momentum, the energy terms are replaced by 

p(ei + 1/2u2 + epi)dvoZ aEi 

- at = lol 

- dEc = -f p(e+1/2u2+ep)G-di dt area 

= dol p(ec + 1/2u2 + e,)Zdvol. 

By Fourier's law, heat conduction per unit area is q = -kdT/dn, where k is a 

proportionality constant, T is temperature, n is the direction normal to the area, 

and the negative sign indicates that heat flux is from higher to lower temperatures. 

104 
In vector form, and integrating over the entire surface of the control volume the 

heat term in equation (120) becomes 

dT -* 

-IC-& &A dQ conduct ed 

dt = -.Lea an 

dT 

area an 

= f k-ii-dz 

= lei $ (k-ii)dvol dT 

dn 

= lo, ’? (ka T)dvol. 

The work wb done by the body forces is found by integrating these forces over 

the entire volume. Since W = SdW = Fdz and F = Jpfdvol, then Pf.b is the 

body force per unit volume in the zi direction and the energy due to this body 

force is 

This leads to 

- -- 

dt dt 

Alternatively, this energy term (in fact all the energy 

(133) 

terms associated with the 

momentum equation ) can be obtained from the momentum equation by taking 

105 
ORIGINAL PAGE IS 

OF POOR QUALITY 

the dot product of the momentum change due to the body forces with the velocity 

k&oraf-tfue am- : -?a: :.,,,, e.- 

Applied to the surface force term in the momentum equation, the technique 

produces 

and that due to the x component of change is 

The energy equation is then, after dropping the integration signs to get the 

differential form, and moving the convection energy term to the left, 

106 
9 Results and Conclusions 

Steger's original algorithm for solving the combined elliptic-hyperbolic equations, 

using fixed boundaries, was shown to work well on certain bodies, for certain 

initial grids. As concavity increased, the initial grids generated by the hyperbolic 

generator, hyg2d, broke down, necessitating the use of some other initial-grid 

generator. This latter generator, not being a function of the body, spaced points 

on the boundary irrespective of body shape, causing loss of orthogonality in some 

areas. 

The author then made modifications to the original algorithm'to increase the 

range of bodies that the combined solver can handle. The first modification was 

to float the outer boundary. That is, at each iteration, after the interior of the grid 

was improved, the hyperbolic component, using the improved grid, generated a 

new and generally improved outer boundary. This allowed the use of initial grids 

with defects in the outer boundary. But, though an increased body complexity 

could now be dealt with, still, for strong enough concavities this algorithm breaks 

down. Additionally, as concavity increased, p, the parameter used to control the 

amount of elliptic or hyperbolic, became more limited in the values it could take 

on. Only values within this range would allow for convergence to solution. Finding 

these values became more difficult as cavities grew in depth. 

A second modification allows the use of a different kind of initial grid, one 

with few or no defects. This is a grid of compressed grid cells where the length 

of cells in the 17 direction is very small compared to the length in the 6 direction. 

By starting with such a defect-free grid and allowing grid cells to grow with each 

iteration, combined with the floating outer boundary it is possible to generate a 

107 
.- 

grid about more complex (higher concavity) bodies. Additionally p can be more 

arbitrarily picked, affording more variability between hyperbolic and elliptic. 

In the last modification to the algorithm, boundary movement within a cavity 

is modified. Instead of being generated by the hyperbolic alone it is modified to 

some degree by a forcing function that moves boundary points about a a chosen 

boundary point so that corresponding boundary points merge. This causes a cut 

to form within the cavity, made up of the outer boundary. This modification is 

not complete. There still remains the problem of forcing the point at which the 

outer boundary joins itself to move far enough from the body such that the entire 

relevant flow field is covered. 

When this is done a grid about a two-dimensional body with deep cavities 

may be achieved. Then, it should also be possible to generate a grid about multiple bodies. By attaching them with imaginary lines, such as to make one inner 

boundary, the cavities so formed may be manageable. This idea is demonstrated 

in figures 66 through 69. Figure 66 shows two bodies connected by a line, which 

forms a continuous inner boundary that outlines the two bodies. Figure 67 is a 

grid generated about this new body. Figures 68 and 69 are detailed views of the 

area surrounding the connection line. The main point here is that there is little 

if any breakdown of the grid, no intersections are noticeable. 

Several methods of identifying grid breakdown have been described. Although 

the logic has been coded, and the first two techniques are in the grid generator 

program ehgrd, there remains yet the problem of using this information to control 

p in order to eliminate the breakdown. One possibility is to simply turn the 

elliptic on (set p high) wherever grid breakdown is detected. This, though, may 

108 
hamper using p for resolution of flow variables. 

109 
Drawing One Grid About Two bodies by Connecting 

Them with a Line to Form One Body 

Y 1' 

Figure 66: Airfoil and slat connected by a line 

file : fl4fp2.p~ 

Y 

T 

Figure 67: 100th pass of the combined solver; boundary is 

floated and grid cell volumes increase 5 percent per iteration 

File : fl4f3.10~~ 

110 
Drawing One Grid About Two bodies by Connecting 

Them with a Line to Form One Body 

Figure 68: 100th pass, detail of final grid in figure 66 

File : f 140.1 Optaps 

111 
Drawing One Grid About Two bodies by Connecting 

Them with a Line to Form One Body 

Figure 69: 100th pass, detail of final grid in figure 66 

File : fl4f3.lOptbps 

112 
10 Notation and Definitions 

uxx 

Linear Finite Difference Operators 

Operat or Symbol Difference Represent ation 

Forward 

(Indexed) 

A,u = U(Z + AX) - U(Z) 

n,u = U(Z;+1) - u(2;) 

%+I - Ui - 

(1404 

(140b) 

(140c) 

Backward vz 

(Indexed) 

(141a) 

(141b) 

(141c) 

Central 

(Indexed) 

Zzu = u(z + Az) - u(z - Az) 

&u = U(Z;+1) - u(z;-1) 

- Ui+1 - %-I 

(142a) 

(142b) 

(142c) 

Central 6, 6,u = U(Z + fAz) - u(z - ~Az) (143) 

Central 6: 6Zu = hz(6,u) (144a) 

=U(Z + Az) - ~u(z) + U(I - Az) (144b) 

(Indexed) 6:u = u(z;+1) - 2u(z;) + U(Zi-1) (144c) 

- u;+1 - 2u; + Ui-1 (144d) 

2nd order accurate, 1st partial derivative 

113 
u, = 

u, = 

6,u 

- + O(AX)~ 2Ax 

6& 

2Ax 

2nd order accurate, 2nd partial derivative 

if u = u(z,y), and if only a finite number of x and y values are required, as is 

generally the case in numerical analysis, these values can be indexed as 

x = z(1) f 21 

Y = Y(l) = Yz 

u = u(1) G Ul 

I = 1,2,. . .Zmax. (149) 

Additionally, if the spacing in x and y values is constant a rectangular grid of 

points is formed, allowing double indexing of u 

- x = x(m) = x, 

- Y = dn) = Yn 

u = u(m,n) G u,,, 

m = 1,2,. . .mmax. 

n = 1,2, ... nmax 

114 
-. 

But if this is not the case a transformation of variables may be performed that 

allows and q such a rectangular grid to exist When z = z(t,q) and y = y(t,q), 

can be indexed as 

t = t(j) tj j = 1,2,. . . jmax, (151) 

7 = q(k) qk k = 1,2,. . Amax, (152) 

This leads to 

t 

If u is some derivative with respect to t or q the finite difference equivalent 

can be stated as 

imax 

i=imin 

i=imin 

Now if u happens to be a derivative with respect to x or y the chain rule can 

be used to transform into t,q variables so that 

ux = qtx + U$lx 

115 
Oj2lGlNAL PAGE IS 

OF POOR QUALITY 

Also by use of the metric transformation equations (, and 7, can be changed to 

the u's neighboring u(j,k), so that if v, = u, v, = v,(j,k), and, for a given value 

of j and k, jo,Ico 

imz 

v8(j,k) = aivrjo+i,k) 

where ai is some rational constant 

i=imin 

and i may take on negative values. 

Convexity and Concavity 

(159) 

Convexity refers to the convex portions of a body. Concavity refers to the 

relative depth and curvature of the concave portions of a body. The algorithm 

used to generate the two-dimensional bodies used in this paper controlled the 

body shape by a parameter named cam. This paramter when set to zero, as in 

figure 70, produces a body with no concave portion. As cam increases the bodies 

produced have what can be called stronger concavities as shown in figures 71, 72 

and 73. 

116 
The Relation of the Parameter Cam to Ccnc:~+ . 'i I 

- Pzirameter tau = 0.5 -. .-... - 

V 

L--+X 

a 

Figure 70: cam= 0.0 

LX 

Figure 72: cam= 0.80 

V 

LX 

Figure 

a 

71: cam= 0.40 

A 

i, 

Figure 73 : cam= 1.0 

117 
. ORiGiWiL. PAX 88 

OF POOR QUALITY 

[l] Garabedian P. R.: Partial Differential Equations, John Wiley, New York, 

1964. 

[2] Thompson, J. F., Warsi, Z. U. A., Masten, C. W.: Numerical Grid Generation 

Foundations and Applications, North-Holland, New York, 19885. 

[3] Thompson, J. F., Thames, F. C., Masten, C. W.: Automatic Numerical Generation of Body-Fitted Curvilinear Coordinate System for Field Containing 

Any Number of Arbitrary Two-Dimensional Bodies. J. Comp. Physics, vol. 

15, no. 3, July 1974, pp 299-319 

[4] Sorenson, R. L., Steger, J. L.: Simplified Clustering of Nonorthogonal Grids 

Generated by Elliptic Partial Differential Equations NASA TM-73252, 1977. 

[5] Steger, J. L.: Generation of Three-Dimensional Body-Fitted Grids by Solving 

Hyperbolic Partial Differential Equations NASA TM-101069, 1989. 

[6] Anderson,D. A., Tannehill, J. C., Pletcher R. H.: Computational Fluid Mechanics, Hemisphere Publishing Corporation, New York, 1984. 

[7] Sorenson, R. L.: A Computer Program to Generate Two-Dimensional Grids 

About Airfoils and Other Shapes by the Use of Poisson's Equation, NASA 

TM-81198, 1980. 

[8] Steger, J. L.: Generation of Body-Fitted Coordinates Using Hyperbolic Partial Differential Equations. SIAM J. Sci. Stat. Comput., vol 1, no 4, December 

1980. 

118 
[9] Van Dalsem, W. R., Steger, J. L.: The Fortified Navier-Stokes Approach, 

proceedings of the Workshhop on Computational Fluid Dynamics, Univ. of 

California, Davis, June 17-18, 1986. 

[lo] Schlicting Herman : Boundary Layer Theory, McGraw-Hill Book Company, 

New York, 1979. 

119