- Email: [email protected]
Algebraic grid generation
137 ALGEBRAIC
GRID GENERATION
ROBERT E. SMITH NASA Langley Research Hampton, VA 23665
Center
ABSTRACT
Three methods dimensional
are described
physical
domains
The methods
domain.
and do not require
are based
They are simply
transfinite
interpolation,
technique.
The primary
of physical
require active
relatively computer
generation described, parameter boundary
The basic mathematical
topology
attention
generation,
and grid derivative
and are called
and the two-boundary
is that they provide Secondly,
the application
the methods structure
is given
and control
computational
or the use of complex
grid spacing.
with
two- and three-
functions
methods
Consequently,
in conjunction
and particular variable
method,
of the methods
few computations.
of grids.
equations
to as algebraic
the multisurface
advantage
in bounded
interpolation
of differential
referred
grid shape and physical
graphics
grids
grid in a rectangular
on mathematical
the solution
variables.
control
for transforming
into a uniform
requirements
of inter-
is advocated
of each method
to surface
function
explicit they
for rapid is
representation,
generation.
Physical
are presented.
NOMENCLATURE
;;z
vector
valued
representation
of surface points
tbz ,,
vector
valued
representation
of surface
e
function
E
magnitude
I
I
relating
variable
f.g,h
control
; -f-f FlcF2 G -1
vector
JtJ
normalized
spanning
between
of vectors
derivatives
arc length to the computational surfaces
tangent
to a spanning
function
functions valued
intermediate
representation vector
valued
integral
of interpolants
Jacobian
matrix,
h K
control parameter
L,N.M
number
inverse
of defining
the physical r,s,t
normalized
s
a set defining
of the physical representations
Jacobian
in the
domain
arc lengths points
on a surface
domain
matrix
in a grid concentration points
domain
of the physical
I,
function
J. and
K
directions
in
a set defining
points
in the computational
a set defining
points
in the physical
vector
valued
representation
orthogonality vector
magnitude
tangent
dependent
coordinates
physical
coordinate
blending
functions delta
domain
function
function
for linear
blending
function
for cubic
increments
curve
functions
blending
computational
linear
from the physical
physical
kronecket
of a surface
coefficient
to a piecewise
variables
domain
domain
interpolation interpolation
coordinates
for a uniform
computational
grid
Superscripts nth, mth partial
n.m
derivative
Subscripts I,JtK
index
for known points
k
index
for surfaces
R
index
in the physical
domain
INTRODUCTION The numerical
solution
of partial
geometries
and with varying
coordinate
systems
and physical equations conditions finding
requirements.
and refine
coordinate
generation," examined.
accuracy
and in this paper
Transfinite is a highly
of surface
described
algebraic definition5
is applied
can complicate
the application
for
both geometric the basic
of boundary
regions.
The process
of
representations
is called
grid generation
methods are
the multisurface
formulas
positions
"grid
method,
in terms of homotopic
and/or
derivatives
of the techniques.
interpolation4
generalized
interpolation
in terms of point
elemants
irragular
the need
that reflect
interpolation,
Interpolation
about
has created
in critical
in discrete
are transfinite
equations
transformations
three algebraic
technique.
and constraints
are the essential
of methods
Coordinate
solution
The methods
scales
transformations
but they can simplify
transformations
and the two-boundary mappings3
characteristic
and associated
of motion,l#'
differential
through
by Cordon
grid generation pioneered a series
and Hall in the early method.
by Steven
Coons.
of univariate
1970's
It is an outgrowth Transfinite
interpolations
where
139
blending functions and the associated parameters derivatives) determine a grid.
(point position and/or
Sriksson6 and Riezi and Sriksson' have adapted
the original transfinite interpolation formulation to use only exterior boundary descriptions and derivatives of certain boundaries.
They have also incorporated
exponentials into the blending functions to concentrate the grid near an exterior boundary.
The developers of the GIM codeSI
tion for grid generation.
use transfinite interpola-
They define boundaries in terms of algebraic
geometric formulas (linear segments, circular arc, conic, etc.) and use linear blending functions for the interior grid computation. The multisurface method 10,ll developed by Peter Eiseman provides formulas for grid definition based on grid descriptions of two boundary surfaces and an arbitrary number of intermediate control Surfaces.
Choosing interpolants
(defined similar to blending functions) and the placement of the control surfaces determines grid shape and spacing.
The multisurface method has been
used by Eiseman in numerous applicationsl2R13 but most notably for computing grids about turbine cascades.12 The two-boundary techniqueln2'14 described by this author is based on the description of two exterior boundaries and the application of either linear or hermite cubic polynomial interpolation to compute the interior grid.
For
cubic interpolation, surface derivatives combined with magnitude coefficients control the orthogonality of the grid at and near the boundaries.
lCowalski,15
applying the two-boundary technique , extended the derivative magnitude coefficients to a functional form for variable orthogonality control. Four additional topics are discussed.
They are surface paramaterization,
grid spacing control, grid topology, and grid computation.
These topics
compliment the basic; mathematical structure of the algebraic grid generation methods and should be considered in their application.
An introduction to
boundary-fitted coordinate systems , which sets the stage for grid generation, precedes the description of the algebraic methods.
BOUNDARY-FITTED COORDINATE TSANSFOSMATIONS The motivation for discrete coordinate transformations or grid generation is the numerical solution of partial differential equations.
Numerical solu-
tions are obtained by either finite difference or finite element techniques, however, the emphasis here is directed at finite difference methods.
Normally,
equations of motion are derived relative to a Cartesian coordinate system. When boundary conditions must be applied on irregular subdomains of the Cartesian coordinate system, and when there are regions within the subdomain
140 that have rapidly varying solutions, it is desirable to transform the equations
to a more appropriate coordinate system.
An ideal coordinate system for
obtaining numerical solutions is a boundary-fitted coordinate system where the physical boundaries of an irregular subdomain transform into exterior boundaries of a rectangular region, and where regions of rapid change are amplified.
If
the bounded subdomain of the Cartesian coordinate system is called the physical domain, then the rectangular coordinate system where a solution is obtained is called the computational domain (Fig. 1).
5
A transformation between the two
L!-+z
LY __-_+
5
Fig. 1. Computational
domain--physical domain
domains is a unique single valued functional relation. symbolically by letting and
5,
17, and
5 = S(xrY,s),
5
x,
Y, and
z
be coordinates in the computational domain, then
rl= Q(X,Y,Z),
5 = S(X,Y,Z),
x = x(5.r1,5), Y = y(S,n,<), z = Z(S,ll,S). The bounds of the computational domain are defined by
o
This is represented
be coordinates in the physical domain
and conversely
o
v
141
Fig. A uniform
grid
2. Computational
(Fig. 2) is superimposed
grid
onto the computational
domain by letting
AC = constantl, Aq = constant2, AC = constant3. Given the functional
the uniform
relations
grid in the computational
grid in the physical
In order to transform
the equations
respect
to the independent
partial
derivatives
if
v, and
u,
equations to
x,
5,
rl, and
w
y, and 5
z
variables
with respect are velocities
of motion,
domain
is transformed
to a corresponding
domain. of motion, x,
Y, and
x,
are transformed
5,
Y, and
then the first derivatives
of
That is
Q, and 2
u,
to first derivatives
by chain differentiation.
derivatives
with
must be transformed
2
to the variables in the
partial
5.
For example,
directions v, and with
w
to
in the with respect
respect
to
142
The matrix
is the Jacobian
5 =
matrix
S(X~Y,Z)~
are known,
rl = rl(x,~,z),
and
then the Jacobian
matrix
the functions. 5
of
Jacobian
can be obtained
5 =
x,
Y, and
s
If the functional
however,
x = xK..rlr5),
Y = y(C,17,5) and
found by differentiating
to explicitly
to determine
matrix
by differentiating
relations
S(X,Y,Z)
can be directly
It is not necessary,
as functions
inverse
of the transformation.
the functions
z = s(5,n,5)
the Jacobian
know
5,
matrix.
rl, and The
143
with respect to
5,
II, and
5.
With these derivatives
J = Transposed of Cofactor (J-1)
,
IJ-? where
IJ-11
IS the Jacobian determinate.
Thus
Higher derivative analysis can be pursued in a similar
provided
/J-11 # 0.
fashion.
For more information on the transformation of partial differential
equations the reader is referred to Reference 11. It is rare to find algebraic expressions of the computational coordinates as functions of the physical coordinates.
The preferred approach is to
express the physical domain as a function of the computational domain and differentiate the physical grid with respect to the computational grid. very important that derivative evaluation be performed and
It is
144 incorporated motion
into the finite difference
in such a manner
References
ALGEBRAIC
17 and 18 address
GRID GENERATION
In this section
three
algebraic
is to outline
techniques
can be compared
for common
control
Transfinite Probably
follows
Intermediate
and individual
and the physical functions
which
are
so that the merit.
domain enhance
the
interpolation and comprehensive
for application
is presented their
physical
section.
features
domain
techniques
structure
may also be postulated.
the most recent
interpolation dynamics
spacing
grid generation
the mathematical
salient
upon the computational
in the previous
of grid
of
are not created.
METHODS similar
The objective
presented
of the equations
errors
this subject.
discussed.
Each case is based
approximation
that geometrically-induced
by Rizzi
and Eriksson.'l
A transformation
format.
domain
description
such as those arising
is a vector-valued
of tramfinite
in computational
This description
from the computational
fluid
generally domain
to the
function
(1)
where
O
The first formulation method"
tional
is based
domain
S
of transfinite
on knowing
C
=
5 t
interpolation,
a sparsely-organized K=M
I&T' 'IIJK' 5
which
we call the "point
set of points
in the COmputa-
and the corresponding
point
set in
145
the physical
domain
Sp =
( xIJK, yIJK,
K=M J=N )I=L zIJK
Note that the set K=l
SC
is
-
J=l I=1 not the uniform
computational
at the intersection each point valued
of three perpendicular
is at the intersection
representation
Each point
grid.
planes,
of three surfaces
on each surface
Fig. 3. Tramfinite
in the computational
domain
and in the physical (Fig. 3).
is
domain
The vector-
is given by
interpolation--point
description
(2a)
= gJ(c,r;,, J=l...N.
(2b)
= ;K(c,rl,, K=l...M.
(2c)
146
Intarpolhtion between the points in the physical domain is performed by defining a set of blending functions:
a1m;
I=l*..L,
BJm;
hl...N,
YK(S)’
K=l. ..M,
with the conditions
aI
= 6
8,$)
= 6JQ,
6
6
If.
=
0,
I
IL'
!z=l...L
IG.1..
.N
t’ R,
= 1, I = II, I!2
JP. = 0, J # k,
6
6 = 1, J = R, Jp. 6
6
KI1 = 0, K # 1,
KR
f
1,
x = 2.
The transfinite interpolation method is the application of the recursive algorithm L (3a)
147
(3bl J=l
(3c)
Each step of the algorithm is a univariate interpolation in one of three possible directions, or the steps can be combined into a single equation. if it is assumed that
Also,
~(~,rl,S) is continuous, the order of the interpolation
direction is not important.
Obviously, a large quantity of geometric informaDeriving appropriate blending functions is
tion is required to define a grid.
the key element and it can vary from one grid problem to another. Eriksson' uses only the outer boundary surfaces (Fig. 4) and out of surface derivatives at certain boundaries to define an interior grid.
This is reason-
able since normally a great deal of geometric information is known at bounding surfaces, but not always away from'them.
Eriksson's presentation is as follows
and is referred to as the "outer surface" method.
n*
%Lrl,O) a?
-
N = 0.1.2... Fig. 4. Transfinite interpolation--outer surface description
148
11=1,2 n=O,l...Q
a=1,2 n=O,l...R
A set of blending
functions
is defined
crp(5)
&=1,2,
n=O,l...P
pm
+1,2,
n=O,l...Q
Ypm
k=l,2,
n=O,l...R
by
with the conditions
amY;“) (5)
a? =6k16nm
and where
the
Aij = 0,
6
functions
i # j,
are defined
"ij = 1,
(i = j).
by
149
The transfinite
interpolation
2 ~lCS,rl.r;)=
algorithm
becomes
P
z1 9.=1 n=O
(4a)
a(") +l,S) II
(4b)
The boundary
sets are
K=M J=N
K=M
K=M
J=l
J=l
J=l
K=M ll=2
K=M
K=M I=L K=l ' 1=l
P.=2 J=N
derivatives
functions
are required
functions
is critical
J=l I=1
I=L
J=l I=1
also, outward
J=N I=L
J=N
at certain
of these boundaries
for this formulation. to the successful
as well
Again the choice
application
.
as the blending
of the blending
of the method.
150 The multisurface
method
The multisurface generating boundary
method
coordinates
surface
$(c,r;).
zltE,F)
and
by Peter
Eiseman
an inner boundary
An arbitrary
are introduced
~2&L$_lcS#r, between
developed
between
;f,t6,C)
humber
to control
(Fig. 5).
Each
is a procedure
surface
sl(<,<)
of internal
and an outer
surfaces
the coordinate surface
for
representation
representation
is such
that
I-: ,
,a. ..___
-___
Tangent a spannj
_____
,’
,’
::
I’_
Fig. 5. The multisurface The physical
domain
can be written
method
as
(5)
where
0(5(1,
o
o
151 The variable n is the independent variable spanning between surfaces.
With
this introduction the description of the multisurface method generally follows that presented in Reference 16. It is assumed that the set of surfaces described above are ordered from bounding surface to bounding surface, and for a fixed corresponding point on each surface.
5
and
5
there is a
The intermediate surfaces are not
coordinate surfaces, but instead are surfaces which are used to establish a field of tangent vectors to the coordinate curve spanning across the surfaces. For the time being, it is assumed that the bounding surfaces are coordinate surfaces.
A smooth interpolation connecting the bounding surfaces results in
a smooth vector field of tangent directions but with unspecified magnitudes. A unique vector field of tangents is obtained by correctly choosing magnitudes which on integration fit precisely the bounding surfaces.
This is demonstrated
with the vector field of tangents given by
1
- zktc,C)
k=l..N, (Fig. 6).
(6)
Y
tr, x
z
Fig. 6. Tangents to a piecewise linear curve and a partition of the spanning variable from the computational domain The coefficients Rk are scalars which determine the magnitude of the vectors but not the direction.
Using the independent variable
direction, a partition
'I1 < n2... < nN_l
n
for the spanning
can be specified in correspondence
with the tangents in Eq. (61. The partitioned variable can be used to represent the tangents as discrete vector-valued functions which map into
Zk(S,C).
given by
The first derivative of
5Lrl.C)
with respect to
'lk n
is
152
N-l
+ +M,
(7)
= c k=l
where
and
6 =0 kP.
k#R,
6 =lk=R. kQ The interpolants
$k('l) are defined exactly like the blending functions in
Eq. (3) but here they are used to describe a derivative function and multiply a tangent vector field.
Integrating Eq. (7) with an initial
rl and
$(C.C.)
yields
N-l
c
(8)
k=l
where
ri GkUl) =
r
Q,(x) dx.
are chosen so that each ) = 1, then the If the magnitudes EkGk('N-l Ek evaluation of Eq. (8) at nN_l reduces to S,(c,<). This allows Eq. (8) to be expressed
N-l
c k=l
Gk(n)
Gk($1)
i
~,+,K.~) - ~kK,c)
1
which is referred to as Eiseman's general multisurface transformation.
(9)
153 The basic
ingredients
of the multisurface
'I1 < '12 ..- < '1N_1' the interpolents to be polynomials
'k surfaces
of degree
is of degree
N + 1.
N
$k in
method
are the partition
and the surfaces
zk&5).
11, the curve connecting
In a systematic
Choosing
the bounding
fashion
N-l
Jl,cm =
I-I
(rl - i-Ii).
i=l ifk
An example '12 = 1,
is a three
surface
$1 = 1 - rl, and
2 GlVl) = rl - +
~(S,rlL)=
transformation
(N - 1 = 2)
and with
rll = 0,
$2 = 11 then
2 , G2(Q) = %
,
1Q5,5) + L-
1
1 - r1(2 - rl)
r1(2 - rl) -
1
n2~,(5,<) + nZs’,&m
where Blvl) = 1 - 2rl + n2,
B,(rl) = 211 - 2u2,
Comparing
the multisurface
the multisurface and allows
method
interpolation
that blending
functions
method
requires
with transfinite interpolants
in only one coordinate can be derived
interpolation
q,(q)
direction.
from Eiseman
(Eq. (3)),
and one set of surfaces, It is apparent
transformation
formulas
154
starting with interpolants.
It is important to remember that the most difficult
aspect of algebraic grid generation is the determination of functions (blending functions interpolants, etc.) which control a grid.
The emphasis in the multi-
surface development is on deriving interpolants which provide satisfactory control.
The two-boundary technique The two-boundary technique has been described by Smith1'2 and Smith and Weigel.13
The technique has common characteristics with Eriksson's formulation
of transfinite interpolation where position and derivatives on exterior boundaries along with blending functions are used to define the physical domain. For the two-boundary technique blending functions are specified to be linear and cubic polynomials as described by Coons.5
These blending functions can
also be derived from Eiseman's two-surface definition and a special modification 16 Later control functions are incorporated of the four-surface definition. to further enhance grid spacing control. The technique is based on defining two nonintersecting surfaces and
z,(c,c) (Fig. 7) where
~,K,r;) =
Y1(W
and
Fig. 7. The two-boundary technique
155
The physical domain is expressed
Explicit forms of the two-boundary technique are linear and hermite cubic interpolation.
The linear form is
2 ~CS,rl,r;)= c k=l
(10)
where
xlm
= 1 -
9,
x,m = 17. The cubic formulation is
(11)
where !-$Ol) = 2n3 - 3n2 + 1,
lJ2ul) = -2n3 + 312, (12) I.1301)= n3 - 2112+ 1,
U*ul)
=
O
o3- n2,
156 the cross product
of surface
derivatives
or normal
derivatives
at the boundaries
is given by
+ i
J:
f
aZ
ayk
+E.5)
-@E,S)
ayk
a=k
-pG5)
The constants boundaries.
Tk
control
For nonzero
the magnitudes Tk
orthogonal
at the boundaries
the effect
of orthogonality
If the magnitudes unsatisfactory 5
(Tk(c,5))
further
are too large,
(Fig. 8). which
NO orthogonality magnitude
variable
The key ingredients are two nonintersecting normal
magnitude
are control section
deals with
surfaces.
bounding
which
surface
is
Tk
forces grid.
and is of
5
and
over the domain.
Unsatisfactory orthogonality magnitude
8. Grid orthogonality
for the two-boundary
functions.
functions
double-valued
Tk of orthogonality
effect
of the
of the physical
to be a function
Satisfactory orthogonality magnitude
Fig.
derivatives
the magnitudes
into the interior
allows
(13)
from this formulation
Increasing
the grid becomes
Kowalski
allows
-@E’S)
of the normal
a grid resulting zk(E,r;).
, k=l,Z.
surfaces,
technique, normal
as it is presented
magnitude
constants
It is later shown that additional govern
the spacing
representation
which
of a grid.
here,
or
ingredients
The following
is used to define
the bounding
157 SURFACE
REPRESENTATION
Surface
AND PARAMETERIZATION
representation
tion methods.
Normally,
is an important the initial
aspect
description
of the algebraic of a surface
grid genera-
is in terms of an
K=M I=L organized
point
set
S =
(Fig. 9).
It is desirable
to
K=l 1=l find a functional
representation
x(r,t)
s(r,t)
=
_
!_ y(r,t)
z(r,t)
with two independent independent process 1.
construct
putational 3. cedure
interpolate
S
with normalized
S, and related
(i.e.,
5
problems
and
5).
to A
is:
arc length or approximate
Fig. 10. Arc length parameterization
functions
relating
the arc lengths
to the corn-
(Fig. 10); and
the data set
(Fig. 11).
contains
domain
(Fig. 10);
Single valued
coordinates
t, which
for many grid generation
representation
such as bicubic
variables
and
the data set
arc length
Fig. 9. Surface
r
from the computational
that is recommended parameterize
normalized
2.
variables
variables
s
with a bidirectional
splines with the arc lengths being
interpolation
pro-
the independent
158
Fig.
Approximate
11. Bi-directional
arc lengths
are computed
interpolation
from the set
S
by l/2
+
rIK = rI-lK
- %-lKJ2
+ (YIK
-
YI-1K)2
1
2
+
(ZIK
+
('IK - 'IK_1j2
- =I-1K)
1
l/2
t
IK = %K-1
tll
+
=
I
= l...L.
K
= l...M.
r
+
(yIK - yIK_lj2
= 0,
rll
Normalized
IK - xIK_1)2
approximate
rIK IK=<'
t
arc lengths
tIK IK=<'
are
O
IK-
0 < tIK 2 1.
K=M After
forming I
sz =
t
the sets
,
Sx = {xIK,rIKrtIK):::
Sy = {y,..rIK,tIKif
I=1
\K=M three bidirectional
interpolations
and
I=1 can be performed
159
for
x(r.t),
and
y(r,t)
are defined
unit interval.
constructing
has the effect
of controlling
surface.
Each surface
variables
to the parametric
called
"control
The intermediate
z(r,tl.
on the unit interval,
and likewise
single-valued
functions"
variables.
variables
and
functions
the location
can have different
5
<
r = f(c)
of an arbitrary functions
and
and
t
on the
t = g (Cl
grid point on the
relating
The functions
and are discussed
r
are defined
the computational
f(5)
in more detail
and
g(C)
are
at a later point.
lB?IFOBMITY In the previous variables
describing
computational
The variable
In a similar between
This is particularly
surfaces
other
B
to defining
unless
some control
applying
control
reference
replaces
the normalized
along the space curve.
of grid points One approach
relative technique.
between
and has no physical
5,
5, and
Tk
domain. yields
along a space curve
(Fig. 12).
Fig. 13. A uniform distribution of grid points with respect to arc length along a spanning curve
it is often desirable
(Fig. 13).
for spanning
entity
the shape of the space curve,
function
the
from the computational
domain
the curve is specified.
distribution
boundaries compute
along
manner
for the two-boundary
for a fixed
in the physical
to parametric
can be paramaterized
functions
a coordinate
Fig. 12. Natural distribution of grid points along a spanning curve
(grid points)
but more complex
surfaces
desirable
the variable
(10) coordinates
In addition
are related
(Eq. (12)) is a mathematical
than it represents
discretizing
from Eq.
spanning
variables
rJ used in the cubic blending
the boundary
Uniformly
computational
surfaces.
variable
to arc length.
meaning
section
n
in the blending
to initially along
of points
is fixed
functions.
specify
Before
a uniform
the space curve connecting
for obtaining
arc length or approximate This establishes
a distribution
This distribution
a uniform
distribution
normalized
arc length
an empirical
relation
between
the is to (5) the
160 variable
*
and the variable
single-valued yields
interpolation
uniform
the blending
n.
for
distributions
functions
Uniformly
discretizing
0, and substituting
of grid points
can be redefined
s, performing
into the blending
along the curve.
in terms of
s
a
functions
Alternately,
where
now
s = e(n)
lJl(l) = 2s3 - 3s2 + 1,
2 U2(11) = 2s3 + 3s ,
P3(‘1) =
s3-
P4(‘1) =
s3 - s2,
s =
s = e(n).
0,
This procedure additional
is complex
because
interpolation.
to the uniform valued
2s2 + s,
distribution
function
it has two steps and the necessity
Control
of the grid spacing
can be accomplished
on the unit interval
distribution
by constructing
for an relative
a single
such that
s = h[e(n)]-
The function
GRID SPACING
h[e(n)]
of a grid in the phyiscal
the computational
"piecewise reader
coordinates
or surface
local control"
for the spanning
is referred
through
to References
approach
is the construction functions
direction.
is primarily
Eiseman
presents
functions
which
Control
in two dimensions
by how
functions,
what he calls
of interpolants
11 and 16 for this approach.
constraints.
technique
affected
into the blending
the derivation
of control
or surface
the two-boundary
domain
are incorporated
constraints.
blending using
function
CONTROL
The spacing
interpolants,
is a control
and the
Another
are embedded
functions and with
in the
are demonstrated
cubic blending
functions. The relationship for the two-boundary
between
the computational
technique
domain
in two dimensions
and the physical
is given by
domain
and
161
x(5.n)
= xl(rl)U1b)
+ x2(r2)P2W
+ T1 %1dr
113
)U (s)+T2
dy
$(r 2
2
)U (s), 4 (14)
+ ~,(r,)k~,(s) - Tl %r dr
y(S,rl) = ~l(rl)!~lW
113
1~ (s) - T2 $r2)U4W, 2
and
Lll(sI= 2s3 - 3s2 + 1,
1-12(s) = -2s3 + 3s2,
U3(S) = s3 - 2s2
+
St
lJ4(s) = s3 - s2,
where
Xl (rl)'Y1(rl)
position on the first boundary as a function of normalized arc length along the boundary position on the second boundary as a function of normalized arc length along the boundary
dxl -&r 1
dyl 1 ,-(rl) 1 drl
first derivative along the first boundary with respect to normalized arc length along the boundary first derivative along the second boundary with respect to normalized arc length along the boundary normal derivative magnitudes for the respective boundaries
Tl'T2 5 = f,(S)
normalized arc length along the first boundary normalized arc length along the second boundary
r2 = f,(S) s = h[e(n)]
arc length along the grid curves connecting the two boundaries
E.11 0(5<1 o
@ I
coordinates from the computational domain
162 Uniformly discretizing
5
and
rl and given the other quantities described
above a corresponding grid is generated in the physical domain from Wq. (14). For grid
curves
connecting the two boundaries (Fig. 141, their relation-
ship and spacing relative to their neighboring grid curves is based on position, derivatives, and derivative magnitudes at the two boundaries, and the blending functions. the
Given that the blending functions are the same for all grid curves,
spacing between the curves is only a function of boundary information.
The boundary positions and derivatives are a function of normalized arc lengths which are Ln turn written as functions of the computational coordinate is the functions
f,(6)
between grid curves.
and
f,(c)
5.
It
that ultimately control the spaciqg
When there is relatively low slope in these functions
(Fig. 14), there is concentration of
grid curves, and when there is relatively
high slope the grid curves are dispersed
(Fig. 15).
In a similar manner
1
1 91
Fig. 14. Spafisg between neighboring grid curves
),I
Fig. 15. Effect of control functions
163 the grid points along a grid curve are distributed by the blending functions. A control function
h[eUl)]
relating
D
to the normalized arc length
deterknines the final grid point distribution along the grid curve (Fig. 16).
Pig. 16. Control of grid points along grid curves The functions
fl(~),f,t~)
and
are called control functions.
h[e(Vl
They should be single valued, smooth, and have smooth derivatives.
Another
condition is that the functions are defined on the unit square (Fig. 17).
OS1
rs
1
t
090
5
n
5
1,o
Fig. 17. Domain for the definition of control functions The control functions can be analytic functions such as
164 A where
the parameter
first boundary. where
control
In general,
spline
in Reference
SIDE BOUNDARY
the two primary as previously
Another
functions
approach
of grid curves
are restrictive for arbitrary
on the unit square.
near the
relative control
This approach
to
is the is
15.
technique
boundaries.
described,
interpolation.
+,nL)
functions
CONSTRAINTS
The two-boundary
(Fig. lea).
the concentration
analytic
can be applied.
use of smoothing described
controls
K
can be constrained This
and then applying
An example
intersecting
by applying
the technique
the recursive
with one side boundary
The two-boundary
= %+,(5)
by boundaries
is accomplished
,s,K))
technique
is performed
,+,(F;)
,-F&'(f.(5)rgr
formulas
constraint
of transfinite
is demonstrated
with the formulation
,s,(C)) ,h(rl)). to)), :,l(f,(SLg,(O)L
h(Q)) ith constraint
+2
s1
t nt
(f,(S,).f,t5) .h(n))
Fig.
lea. step one in the two-boundary
The second
step is from the transfinite
linear blending
$(5,nn5)
function
= +'W
(Fig.
Fig. 18b. Side boundary constraint
technique interpolation
formulation
18b)
+ (1 - 5) 8;(f1KLf2WW 1
- ~~~~~f,~~,.g,~0,,,~~~f2~5)'42~0~~,h~~~~
1
.
with a
165
GRID GENERATION TOPOLOGY The algebraic grid generation techniques that are presented are defined with the assumption that a uniform rectangular computational domain transforms into a physical domain.
Also, exterior boundaries of the computational
domain transform into boundaries on the physical domain.
Consequently the
topology of the physical domain strongly influences how a grid generation technique is applied.
It is obvious that single six-sided box (computational
domain) or a square in two dimensions is not going to transform into all physical domains.
Further, in certain cases, transformations can only be made
by introducing singularities.
Problems most often arise when there are
closed boundaries and in this section some topological considerations are described. For boundaries in two dimensions, there are two primary types of physical domains that transform from a square computational domain.
They are O-type
domains and C-type domains, and are more commonly referred to as O-grids and C-grids.
Several two-dimensional domains are described schematically in
Figure 19 with corresponding boundary numbers in the computational and physica domains.
Also, multiple computational and physical domains can be coupled
with a common boundary.
A resulting grid as well as grid derivatives should
be continuous across a common boundary.
If there is a discontinuity, special
consideration should be taken in the finite difference procedure for the solution of the equations of motion.
2 computational domain
physical domain
4
3
1
Fig. 19a. Simple O-grid
166
4
6
physical domain
computational domain
6
5
5
1
3
2
I
I Fig. 19b. Simple C-grid
6 physical domain
computational domain
a
11
2
4
13
5
I Fig.19c. Hultiple,body O-grid
8 computational domain
9
I 10
4m_
Fig. 19d. Multiple body C-grid
5
6
3
2
7 1
167 I
I
8 ‘7
9’
computational 12 domain 11
i 10 i 4 6
5
1,2,3
I
I
Fig. lge. Combination domain C-grid
2
I I
computational 3
6
2
6
domain 4j7
8
I I
1
I I
5 .~.__
-__.
10
11
1.
12
9
Fig. lgf. L-shape domain For closed boundaries in three dimensions two suitable types of domains are O-O domains and C-O domains.
The topologies associated with closed
boundary domains in three dimensions is more complex than for two dimensions. A primary reason is that a planar surface from the computational domain will not transform into a closed three-dimensional surface without introducing 6 singularities. Rizzi and Eriksson extensively discuss the problems of three-dimensional closed surface topology and associated singularities and the reader is referred to Reference 6. A final note on topology and singularities is that every effort should be made to place an unavoidable singularity in a region where there is little change in the basic equations of motion.
Near a singularity there are large
changes in the derivatives of the computational coordinates with respect to
168 the physical points
coordinates.
These
can lead to inaccuracy
ference
procedure
large changes
and/or possibly
for the solution
between
neighboring
instability
of the equations
grid
in a finite dif-
of motion.
GRID COMPUTATION Algebraic data.
grid generation
These
data."
data can be divided
Fixed data describes
variable
data describes
control
functions.
cursor
Variable
grid derivatives
tory grid characteristics The algebraic ment because given
rate between because grid.
the computation What
the computer
An aside point
graphics
data which
whereas
surfaces
is ideally
control
and "variable
surfaces,
control
of input
suited
of for this
a grid can be modified
terminal,
the resulting
and the data again modified
grid and
until
satisfac-
are achieved. methods
is explicit is necessary,
and graphics
of the large number
work well
can be worked
and the resulting
control
in an interactive
with no iteration however,
terminal
of line segments
for the algebraic
the grid characteristics grid points)
computer
grid generation
grid solution.
"fixed data"
such as internal
observed,
a large quantity
such as bounding
input at a graphics
visually
require
into two types:
quantities
or numeric
generally
quantities
Interactive
type of application. with
methods
necessary
be high
(9600 baud or greater)
that are discussed
out on a relative directly
for a
is that the communication
that must be displayed
methods
environ-
applied
for a is that
small grid to compute
(few a larger
grid.
DISCUSSION
AND CONCLUSIONS
The three algebraic basically Blending
functions,
geometric
procedures
govern
to where
singularities
singularities
procedure
along with physical
of a uniform
the most
of the physical
there are closed boundaries,
in the finite difference
functions
The methods
are
characteristics.
rectangular
are relatively
extensive
simple
computational
of generality.
procedure,
is the topology
be taken relative
grid.
that are described
common
and do not require
and they have a high degree
grid generation where
or control
they are explicit
Next to the computational
several
the transformation
grid into a physical
to understand,
techniques
with
interpolants,
constraints
computational
effort,
grid generation
interpolation
exist
important domain.
consideration
are introduced.
Care should
and the corresponding
for the solution
for
For three dimensions
of the equations
effect of
169 Grid topology
motion.
importance Another means
is not extensively
consideration
is that grid derivatives
that each step in a grid generation Wiqqles
result.
discussed
in this paper,
but its
is emphasized.
is one step propagate
method
must be smooth. must produce
This
a smooth
into the next step and finally
into
the grid. A final point
is that the use of interactive
generation
is highly
generation
is truely
control
is adaptive,
advantageous.
and require
interactive applications
software
is an important
graphics
for grid
when grid
and the grid
in the control
grid generation
few computations,
This implies
of motion
intervention
Since algebraic
relatively
environment.
human
computer
the time is reached
coupled with the equations instantaneous
in the next best approach. explicit
Until
methods
process are
they work very well in an
that the development aspect of algebraic
of computer
grid generation.
REFERENCES 1.
2. 3. 4.
5. 6.
7.
8.
9.
10. 11. 12. 13.
Smith, R. E., Two-Boundary Grid Generation for the Solution of the ThreeDimensional Navier-Stokes Equations, Ph.D. Dissertation, Old Dominion University, May 1981. Smith, R. E., "Two-Boundary Grid Generation for the Solution of the ThreeDimensional Navier-Stokes Equations," NASA TM-83123, May 1981. Dugundji, James, Topology, Allyn and Bacon Inc., 1968. Gordon, W. J. and Hall, C. A., "Construction of Curvilinear Coordinate Systems and Application to Mesh Generation," International Journal for Numerical Methods in Engineering, Vol. 7, 461-477, 1973. Coons, S. A., "Surfaces for Computer Aided Design of Space Forms," MAC TR-41, Contract No. AF-33(6000-42859) MIT, June 1967. Eriksson, Lars-Erik, "Three-Dimensional Spline-Generated Coordinate Transformations for Grids Around Wing-Body Configurations," Numerical Grid Generation Techniques, NASA CP 2166, 1980. Rizzi, A. and Eriksson, L. E., "Transfinite Mesh Generation and Damped Euler Equation Algorithm for Transonic Flow Around Wing-Body Configurations,' AIAA 5th Computational Fluid Dynamics Conference, June 1981, Palo Alto, California. Grid Generation Anderson, P. G. and Spradley, L. W., "Finite-Difference by Multivariate Blending Function Interpolation," Numerical Grid Generation Techniques, NASA CP 2166, 1980. Spradley, L. W. Stalnaker, J. F., and Ratliff, A. W., "Computation of Three-Dimensional Viscous Flows with the Navier-Stokes Equations," AIAA Paper 80-1348, AIAA 13th Fluid and Plasma Dynamics Conference, Snowmass, Colorado. Eisemen, P. R., "A Multi-Surface Method of Coordinate Generation," Journal of Computational Physics, Vol. 33, No. 1, October 1979. Eiseman. P. R.. "Geometric Methods in Computational Fluid Dynamics," ICASE Report No. 81-11, April 1980. _ Eiseman, P. R., "A Coordinate System for a Viscous Transonic Cascade Analysis," Journal of Computational Physics, Vol. 29, 1978. Coordinates About Wings," Fourth AIAA Eiseman, P. R., "Three-Dimensional Computer Fluid Dynamics Conference, Williamsburg, Virginia, July 1979.
170
14. Smith, R. E. and W&gel, B. L., "Analytic and Approximate Boundary-Fitted Coordinate Systems for Fluid Flow Simulation," AIAA Paper 80-0192, AIAA 18th Aerospace Sciences Meeting, January 1980, Pasadena, California. 15. Kowalski, E. J., "Boundary-Fitted Coordinate Systems for Arbitrary Computational Regions," Numerical Grid Generation Techniques, NASA CP 2166, 1980. 16. Eiseman, P. R. and Smith, R. E., "Mesh Generation Using Algebraic Techniques," Numerical Grid Generation Techniques, NASA CP 2166, 1980. 17. Steger, J. L., Implicit Finite Difference Simulation of Flow About Arbitrary Geometries with Applications to Airfoils," AIAA Paper 77-665, AIAA 10th Fluid and Plasma Dynamics Conference, Albuquerque, New Mexico, June 1977. 18. Iiindman, R. G., "Geometrically Induced Errors and their Relationship to the Form of the Governing Equations and the Treatment of Generalized Mappings," AIAA Paper 81-1008, AIAA 5th Computational Fluid Dynamics Conference, Palo Alto, California, June 1981. 19. Smith, R. E., Kndlinski, R. A., and Everton, E. L., "A Grid Spacing Control Technique for Algebraic Grid Generation Methods," AIAA Paper 82-0226, AIAA 20th Aerospace Sciences Meeting, January 1982, Orlando, Florida.
GRID GENERATION
ROBERT E. SMITH NASA Langley Research Hampton, VA 23665
Center
ABSTRACT
Three methods dimensional
are described
physical
domains
The methods
domain.
and do not require
are based
They are simply
transfinite
interpolation,
technique.
The primary
of physical
require active
relatively computer
generation described, parameter boundary
The basic mathematical
topology
attention
generation,
and grid derivative
and are called
and the two-boundary
is that they provide Secondly,
the application
the methods structure
is given
and control
computational
or the use of complex
grid spacing.
with
two- and three-
functions
methods
Consequently,
in conjunction
and particular variable
method,
of the methods
few computations.
of grids.
equations
to as algebraic
the multisurface
advantage
in bounded
interpolation
of differential
referred
grid shape and physical
graphics
grids
grid in a rectangular
on mathematical
the solution
variables.
control
for transforming
into a uniform
requirements
of inter-
is advocated
of each method
to surface
function
explicit they
for rapid is
representation,
generation.
Physical
are presented.
NOMENCLATURE
;;z
vector
valued
representation
of surface points
tbz ,,
vector
valued
representation
of surface
e
function
E
magnitude
I
I
relating
variable
f.g,h
control
; -f-f FlcF2 G -1
vector
JtJ
normalized
spanning
between
of vectors
derivatives
arc length to the computational surfaces
tangent
to a spanning
function
functions valued
intermediate
representation vector
valued
integral
of interpolants
Jacobian
matrix,
h K
control parameter
L,N.M
number
inverse
of defining
the physical r,s,t
normalized
s
a set defining
of the physical representations
Jacobian
in the
domain
arc lengths points
on a surface
domain
matrix
in a grid concentration points
domain
of the physical
I,
function
J. and
K
directions
in
a set defining
points
in the computational
a set defining
points
in the physical
vector
valued
representation
orthogonality vector
magnitude
tangent
dependent
coordinates
physical
coordinate
blending
functions delta
domain
function
function
for linear
blending
function
for cubic
increments
curve
functions
blending
computational
linear
from the physical
physical
kronecket
of a surface
coefficient
to a piecewise
variables
domain
domain
interpolation interpolation
coordinates
for a uniform
computational
grid
Superscripts nth, mth partial
n.m
derivative
Subscripts I,JtK
index
for known points
k
index
for surfaces
R
index
in the physical
domain
INTRODUCTION The numerical
solution
of partial
geometries
and with varying
coordinate
systems
and physical equations conditions finding
requirements.
and refine
coordinate
generation," examined.
accuracy
and in this paper
Transfinite is a highly
of surface
described
algebraic definition5
is applied
can complicate
the application
for
both geometric the basic
of boundary
regions.
The process
of
representations
is called
grid generation
methods are
the multisurface
formulas
positions
"grid
method,
in terms of homotopic
and/or
derivatives
of the techniques.
interpolation4
generalized
interpolation
in terms of point
elemants
irragular
the need
that reflect
interpolation,
Interpolation
about
has created
in critical
in discrete
are transfinite
equations
transformations
three algebraic
technique.
and constraints
are the essential
of methods
Coordinate
solution
The methods
scales
transformations
but they can simplify
transformations
and the two-boundary mappings3
characteristic
and associated
of motion,l#'
differential
through
by Cordon
grid generation pioneered a series
and Hall in the early method.
by Steven
Coons.
of univariate
1970's
It is an outgrowth Transfinite
interpolations
where
139
blending functions and the associated parameters derivatives) determine a grid.
(point position and/or
Sriksson6 and Riezi and Sriksson' have adapted
the original transfinite interpolation formulation to use only exterior boundary descriptions and derivatives of certain boundaries.
They have also incorporated
exponentials into the blending functions to concentrate the grid near an exterior boundary.
The developers of the GIM codeSI
tion for grid generation.
use transfinite interpola-
They define boundaries in terms of algebraic
geometric formulas (linear segments, circular arc, conic, etc.) and use linear blending functions for the interior grid computation. The multisurface method 10,ll developed by Peter Eiseman provides formulas for grid definition based on grid descriptions of two boundary surfaces and an arbitrary number of intermediate control Surfaces.
Choosing interpolants
(defined similar to blending functions) and the placement of the control surfaces determines grid shape and spacing.
The multisurface method has been
used by Eiseman in numerous applicationsl2R13 but most notably for computing grids about turbine cascades.12 The two-boundary techniqueln2'14 described by this author is based on the description of two exterior boundaries and the application of either linear or hermite cubic polynomial interpolation to compute the interior grid.
For
cubic interpolation, surface derivatives combined with magnitude coefficients control the orthogonality of the grid at and near the boundaries.
lCowalski,15
applying the two-boundary technique , extended the derivative magnitude coefficients to a functional form for variable orthogonality control. Four additional topics are discussed.
They are surface paramaterization,
grid spacing control, grid topology, and grid computation.
These topics
compliment the basic; mathematical structure of the algebraic grid generation methods and should be considered in their application.
An introduction to
boundary-fitted coordinate systems , which sets the stage for grid generation, precedes the description of the algebraic methods.
BOUNDARY-FITTED COORDINATE TSANSFOSMATIONS The motivation for discrete coordinate transformations or grid generation is the numerical solution of partial differential equations.
Numerical solu-
tions are obtained by either finite difference or finite element techniques, however, the emphasis here is directed at finite difference methods.
Normally,
equations of motion are derived relative to a Cartesian coordinate system. When boundary conditions must be applied on irregular subdomains of the Cartesian coordinate system, and when there are regions within the subdomain
140 that have rapidly varying solutions, it is desirable to transform the equations
to a more appropriate coordinate system.
An ideal coordinate system for
obtaining numerical solutions is a boundary-fitted coordinate system where the physical boundaries of an irregular subdomain transform into exterior boundaries of a rectangular region, and where regions of rapid change are amplified.
If
the bounded subdomain of the Cartesian coordinate system is called the physical domain, then the rectangular coordinate system where a solution is obtained is called the computational domain (Fig. 1).
5
A transformation between the two
L!-+z
LY __-_+
5
Fig. 1. Computational
domain--physical domain
domains is a unique single valued functional relation. symbolically by letting and
5,
17, and
5 = S(xrY,s),
5
x,
Y, and
z
be coordinates in the computational domain, then
rl= Q(X,Y,Z),
5 = S(X,Y,Z),
x = x(5.r1,5), Y = y(S,n,<), z = Z(S,ll,S). The bounds of the computational domain are defined by
o
This is represented
be coordinates in the physical domain
and conversely
o
v
141
Fig. A uniform
grid
2. Computational
(Fig. 2) is superimposed
grid
onto the computational
domain by letting
AC = constantl, Aq = constant2, AC = constant3. Given the functional
the uniform
relations
grid in the computational
grid in the physical
In order to transform
the equations
respect
to the independent
partial
derivatives
if
v, and
u,
equations to
x,
5,
rl, and
w
y, and 5
z
variables
with respect are velocities
of motion,
domain
is transformed
to a corresponding
domain. of motion, x,
Y, and
x,
are transformed
5,
Y, and
then the first derivatives
of
That is
Q, and 2
u,
to first derivatives
by chain differentiation.
derivatives
with
must be transformed
2
to the variables in the
partial
5.
For example,
directions v, and with
w
to
in the with respect
respect
to
142
The matrix
is the Jacobian
5 =
matrix
S(X~Y,Z)~
are known,
rl = rl(x,~,z),
and
then the Jacobian
matrix
the functions. 5
of
Jacobian
can be obtained
5 =
x,
Y, and
s
If the functional
however,
x = xK..rlr5),
Y = y(C,17,5) and
found by differentiating
to explicitly
to determine
matrix
by differentiating
relations
S(X,Y,Z)
can be directly
It is not necessary,
as functions
inverse
of the transformation.
the functions
z = s(5,n,5)
the Jacobian
know
5,
matrix.
rl, and The
143
with respect to
5,
II, and
5.
With these derivatives
J = Transposed of Cofactor (J-1)
,
IJ-? where
IJ-11
IS the Jacobian determinate.
Thus
Higher derivative analysis can be pursued in a similar
provided
/J-11 # 0.
fashion.
For more information on the transformation of partial differential
equations the reader is referred to Reference 11. It is rare to find algebraic expressions of the computational coordinates as functions of the physical coordinates.
The preferred approach is to
express the physical domain as a function of the computational domain and differentiate the physical grid with respect to the computational grid. very important that derivative evaluation be performed and
It is
144 incorporated motion
into the finite difference
in such a manner
References
ALGEBRAIC
17 and 18 address
GRID GENERATION
In this section
three
algebraic
is to outline
techniques
can be compared
for common
control
Transfinite Probably
follows
Intermediate
and individual
and the physical functions
which
are
so that the merit.
domain enhance
the
interpolation and comprehensive
for application
is presented their
physical
section.
features
domain
techniques
structure
may also be postulated.
the most recent
interpolation dynamics
spacing
grid generation
the mathematical
salient
upon the computational
in the previous
of grid
of
are not created.
METHODS similar
The objective
presented
of the equations
errors
this subject.
discussed.
Each case is based
approximation
that geometrically-induced
by Rizzi
and Eriksson.'l
A transformation
format.
domain
description
such as those arising
is a vector-valued
of tramfinite
in computational
This description
from the computational
fluid
generally domain
to the
function
(1)
where
O
The first formulation method"
tional
is based
domain
S
of transfinite
on knowing
C
=
5 t
interpolation,
a sparsely-organized K=M
I&T' 'IIJK' 5
which
we call the "point
set of points
in the COmputa-
and the corresponding
point
set in
145
the physical
domain
Sp =
( xIJK, yIJK,
K=M J=N )I=L zIJK
Note that the set K=l
SC
is
-
J=l I=1 not the uniform
computational
at the intersection each point valued
of three perpendicular
is at the intersection
representation
Each point
grid.
planes,
of three surfaces
on each surface
Fig. 3. Tramfinite
in the computational
domain
and in the physical (Fig. 3).
is
domain
The vector-
is given by
interpolation--point
description
(2a)
= gJ(c,r;,, J=l...N.
(2b)
= ;K(c,rl,, K=l...M.
(2c)
146
Intarpolhtion between the points in the physical domain is performed by defining a set of blending functions:
a1m;
I=l*..L,
BJm;
hl...N,
YK(S)’
K=l. ..M,
with the conditions
aI
= 6
8,$)
= 6JQ,
6
6
If.
=
0,
I
IL'
!z=l...L
IG.1..
.N
t’ R,
= 1, I = II, I!2
JP. = 0, J # k,
6
6 = 1, J = R, Jp. 6
6
KI1 = 0, K # 1,
KR
f
1,
x = 2.
The transfinite interpolation method is the application of the recursive algorithm L (3a)
147
(3bl J=l
(3c)
Each step of the algorithm is a univariate interpolation in one of three possible directions, or the steps can be combined into a single equation. if it is assumed that
Also,
~(~,rl,S) is continuous, the order of the interpolation
direction is not important.
Obviously, a large quantity of geometric informaDeriving appropriate blending functions is
tion is required to define a grid.
the key element and it can vary from one grid problem to another. Eriksson' uses only the outer boundary surfaces (Fig. 4) and out of surface derivatives at certain boundaries to define an interior grid.
This is reason-
able since normally a great deal of geometric information is known at bounding surfaces, but not always away from'them.
Eriksson's presentation is as follows
and is referred to as the "outer surface" method.
n*
%Lrl,O) a?
-
N = 0.1.2... Fig. 4. Transfinite interpolation--outer surface description
148
11=1,2 n=O,l...Q
a=1,2 n=O,l...R
A set of blending
functions
is defined
crp(5)
&=1,2,
n=O,l...P
pm
+1,2,
n=O,l...Q
Ypm
k=l,2,
n=O,l...R
by
with the conditions
amY;“) (5)
a? =6k16nm
and where
the
Aij = 0,
6
functions
i # j,
are defined
"ij = 1,
(i = j).
by
149
The transfinite
interpolation
2 ~lCS,rl.r;)=
algorithm
becomes
P
z1 9.=1 n=O
(4a)
a(") +l,S) II
(4b)
The boundary
sets are
K=M J=N
K=M
K=M
J=l
J=l
J=l
K=M ll=2
K=M
K=M I=L K=l ' 1=l
P.=2 J=N
derivatives
functions
are required
functions
is critical
J=l I=1
I=L
J=l I=1
also, outward
J=N I=L
J=N
at certain
of these boundaries
for this formulation. to the successful
as well
Again the choice
application
.
as the blending
of the blending
of the method.
150 The multisurface
method
The multisurface generating boundary
method
coordinates
surface
$(c,r;).
zltE,F)
and
by Peter
Eiseman
an inner boundary
An arbitrary
are introduced
~2&L$_lcS#r, between
developed
between
;f,t6,C)
humber
to control
(Fig. 5).
Each
is a procedure
surface
sl(<,<)
of internal
and an outer
surfaces
the coordinate surface
for
representation
representation
is such
that
I-: ,
,a. ..___
-___
Tangent a spannj
_____
,’
,’
::
I’_
Fig. 5. The multisurface The physical
domain
can be written
method
as
(5)
where
0(5(1,
o
o
151 The variable n is the independent variable spanning between surfaces.
With
this introduction the description of the multisurface method generally follows that presented in Reference 16. It is assumed that the set of surfaces described above are ordered from bounding surface to bounding surface, and for a fixed corresponding point on each surface.
5
and
5
there is a
The intermediate surfaces are not
coordinate surfaces, but instead are surfaces which are used to establish a field of tangent vectors to the coordinate curve spanning across the surfaces. For the time being, it is assumed that the bounding surfaces are coordinate surfaces.
A smooth interpolation connecting the bounding surfaces results in
a smooth vector field of tangent directions but with unspecified magnitudes. A unique vector field of tangents is obtained by correctly choosing magnitudes which on integration fit precisely the bounding surfaces.
This is demonstrated
with the vector field of tangents given by
1
- zktc,C)
k=l..N, (Fig. 6).
(6)
Y
tr, x
z
Fig. 6. Tangents to a piecewise linear curve and a partition of the spanning variable from the computational domain The coefficients Rk are scalars which determine the magnitude of the vectors but not the direction.
Using the independent variable
direction, a partition
'I1 < n2... < nN_l
n
for the spanning
can be specified in correspondence
with the tangents in Eq. (61. The partitioned variable can be used to represent the tangents as discrete vector-valued functions which map into
Zk(S,C).
given by
The first derivative of
5Lrl.C)
with respect to
'lk n
is
152
N-l
+ +M,
(7)
= c k=l
where
and
6 =0 kP.
k#R,
6 =lk=R. kQ The interpolants
$k('l) are defined exactly like the blending functions in
Eq. (3) but here they are used to describe a derivative function and multiply a tangent vector field.
Integrating Eq. (7) with an initial
rl and
$(C.C.)
yields
N-l
c
(8)
k=l
where
ri GkUl) =
r
Q,(x) dx.
are chosen so that each ) = 1, then the If the magnitudes EkGk('N-l Ek evaluation of Eq. (8) at nN_l reduces to S,(c,<). This allows Eq. (8) to be expressed
N-l
c k=l
Gk(n)
Gk($1)
i
~,+,K.~) - ~kK,c)
1
which is referred to as Eiseman's general multisurface transformation.
(9)
153 The basic
ingredients
of the multisurface
'I1 < '12 ..- < '1N_1' the interpolents to be polynomials
'k surfaces
of degree
is of degree
N + 1.
N
$k in
method
are the partition
and the surfaces
zk&5).
11, the curve connecting
In a systematic
Choosing
the bounding
fashion
N-l
Jl,cm =
I-I
(rl - i-Ii).
i=l ifk
An example '12 = 1,
is a three
surface
$1 = 1 - rl, and
2 GlVl) = rl - +
~(S,rlL)=
transformation
(N - 1 = 2)
and with
rll = 0,
$2 = 11 then
2 , G2(Q) = %
,
1Q5,5) + L-
1
1 - r1(2 - rl)
r1(2 - rl) -
1
n2~,(5,<) + nZs’,&m
where Blvl) = 1 - 2rl + n2,
B,(rl) = 211 - 2u2,
Comparing
the multisurface
the multisurface and allows
method
interpolation
that blending
functions
method
requires
with transfinite interpolants
in only one coordinate can be derived
interpolation
q,(q)
direction.
from Eiseman
(Eq. (3)),
and one set of surfaces, It is apparent
transformation
formulas
154
starting with interpolants.
It is important to remember that the most difficult
aspect of algebraic grid generation is the determination of functions (blending functions interpolants, etc.) which control a grid.
The emphasis in the multi-
surface development is on deriving interpolants which provide satisfactory control.
The two-boundary technique The two-boundary technique has been described by Smith1'2 and Smith and Weigel.13
The technique has common characteristics with Eriksson's formulation
of transfinite interpolation where position and derivatives on exterior boundaries along with blending functions are used to define the physical domain. For the two-boundary technique blending functions are specified to be linear and cubic polynomials as described by Coons.5
These blending functions can
also be derived from Eiseman's two-surface definition and a special modification 16 Later control functions are incorporated of the four-surface definition. to further enhance grid spacing control. The technique is based on defining two nonintersecting surfaces and
z,(c,c) (Fig. 7) where
~,K,r;) =
Y1(W
and
Fig. 7. The two-boundary technique
155
The physical domain is expressed
Explicit forms of the two-boundary technique are linear and hermite cubic interpolation.
The linear form is
2 ~CS,rl,r;)= c k=l
(10)
where
xlm
= 1 -
9,
x,m = 17. The cubic formulation is
(11)
where !-$Ol) = 2n3 - 3n2 + 1,
lJ2ul) = -2n3 + 312, (12) I.1301)= n3 - 2112+ 1,
U*ul)
=
O
o3- n2,
156 the cross product
of surface
derivatives
or normal
derivatives
at the boundaries
is given by
+ i
J:
f
aZ
ayk
+E.5)
-@E,S)
ayk
a=k
-pG5)
The constants boundaries.
Tk
control
For nonzero
the magnitudes Tk
orthogonal
at the boundaries
the effect
of orthogonality
If the magnitudes unsatisfactory 5
(Tk(c,5))
further
are too large,
(Fig. 8). which
NO orthogonality magnitude
variable
The key ingredients are two nonintersecting normal
magnitude
are control section
deals with
surfaces.
bounding
which
surface
is
Tk
forces grid.
and is of
5
and
over the domain.
Unsatisfactory orthogonality magnitude
8. Grid orthogonality
for the two-boundary
functions.
functions
double-valued
Tk of orthogonality
effect
of the
of the physical
to be a function
Satisfactory orthogonality magnitude
Fig.
derivatives
the magnitudes
into the interior
allows
(13)
from this formulation
Increasing
the grid becomes
Kowalski
allows
-@E’S)
of the normal
a grid resulting zk(E,r;).
, k=l,Z.
surfaces,
technique, normal
as it is presented
magnitude
constants
It is later shown that additional govern
the spacing
representation
which
of a grid.
here,
or
ingredients
The following
is used to define
the bounding
157 SURFACE
REPRESENTATION
Surface
AND PARAMETERIZATION
representation
tion methods.
Normally,
is an important the initial
aspect
description
of the algebraic of a surface
grid genera-
is in terms of an
K=M I=L organized
point
set
S =
(Fig. 9).
It is desirable
to
K=l 1=l find a functional
representation
x(r,t)
s(r,t)
=
_
!_ y(r,t)
z(r,t)
with two independent independent process 1.
construct
putational 3. cedure
interpolate
S
with normalized
S, and related
(i.e.,
5
problems
and
5).
to A
is:
arc length or approximate
Fig. 10. Arc length parameterization
functions
relating
the arc lengths
to the corn-
(Fig. 10); and
the data set
(Fig. 11).
contains
domain
(Fig. 10);
Single valued
coordinates
t, which
for many grid generation
representation
such as bicubic
variables
and
the data set
arc length
Fig. 9. Surface
r
from the computational
that is recommended parameterize
normalized
2.
variables
variables
s
with a bidirectional
splines with the arc lengths being
interpolation
pro-
the independent
158
Fig.
Approximate
11. Bi-directional
arc lengths
are computed
interpolation
from the set
S
by l/2
+
rIK = rI-lK
- %-lKJ2
+ (YIK
-
YI-1K)2
1
2
+
(ZIK
+
('IK - 'IK_1j2
- =I-1K)
1
l/2
t
IK = %K-1
tll
+
=
I
= l...L.
K
= l...M.
r
+
(yIK - yIK_lj2
= 0,
rll
Normalized
IK - xIK_1)2
approximate
rIK IK=<'
t
arc lengths
tIK IK=<'
are
O
IK-
0 < tIK 2 1.
K=M After
forming I
sz =
t
the sets
,
Sx = {xIK,rIKrtIK):::
Sy = {y,..rIK,tIKif
I=1
\K=M three bidirectional
interpolations
and
I=1 can be performed
159
for
x(r.t),
and
y(r,t)
are defined
unit interval.
constructing
has the effect
of controlling
surface.
Each surface
variables
to the parametric
called
"control
The intermediate
z(r,tl.
on the unit interval,
and likewise
single-valued
functions"
variables.
variables
and
functions
the location
can have different
5
<
r = f(c)
of an arbitrary functions
and
and
t
on the
t = g (Cl
grid point on the
relating
The functions
and are discussed
r
are defined
the computational
f(5)
in more detail
and
g(C)
are
at a later point.
lB?IFOBMITY In the previous variables
describing
computational
The variable
In a similar between
This is particularly
surfaces
other
B
to defining
unless
some control
applying
control
reference
replaces
the normalized
along the space curve.
of grid points One approach
relative technique.
between
and has no physical
5,
5, and
Tk
domain. yields
along a space curve
(Fig. 12).
Fig. 13. A uniform distribution of grid points with respect to arc length along a spanning curve
it is often desirable
(Fig. 13).
for spanning
entity
the shape of the space curve,
function
the
from the computational
domain
the curve is specified.
distribution
boundaries compute
along
manner
for the two-boundary
for a fixed
in the physical
to parametric
can be paramaterized
functions
a coordinate
Fig. 12. Natural distribution of grid points along a spanning curve
(grid points)
but more complex
surfaces
desirable
the variable
(10) coordinates
In addition
are related
(Eq. (12)) is a mathematical
than it represents
discretizing
from Eq.
spanning
variables
rJ used in the cubic blending
the boundary
Uniformly
computational
surfaces.
variable
to arc length.
meaning
section
n
in the blending
to initially along
of points
is fixed
functions.
specify
Before
a uniform
the space curve connecting
for obtaining
arc length or approximate This establishes
a distribution
This distribution
a uniform
distribution
normalized
arc length
an empirical
relation
between
the is to (5) the
160 variable
*
and the variable
single-valued yields
interpolation
uniform
the blending
n.
for
distributions
functions
Uniformly
discretizing
0, and substituting
of grid points
can be redefined
s, performing
into the blending
along the curve.
in terms of
s
a
functions
Alternately,
where
now
s = e(n)
lJl(l) = 2s3 - 3s2 + 1,
2 U2(11) = 2s3 + 3s ,
P3(‘1) =
s3-
P4(‘1) =
s3 - s2,
s =
s = e(n).
0,
This procedure additional
is complex
because
interpolation.
to the uniform valued
2s2 + s,
distribution
function
it has two steps and the necessity
Control
of the grid spacing
can be accomplished
on the unit interval
distribution
by constructing
for an relative
a single
such that
s = h[e(n)]-
The function
GRID SPACING
h[e(n)]
of a grid in the phyiscal
the computational
"piecewise reader
coordinates
or surface
local control"
for the spanning
is referred
through
to References
approach
is the construction functions
direction.
is primarily
Eiseman
presents
functions
which
Control
in two dimensions
by how
functions,
what he calls
of interpolants
11 and 16 for this approach.
constraints.
technique
affected
into the blending
the derivation
of control
or surface
the two-boundary
domain
are incorporated
constraints.
blending using
function
CONTROL
The spacing
interpolants,
is a control
and the
Another
are embedded
functions and with
in the
are demonstrated
cubic blending
functions. The relationship for the two-boundary
between
the computational
technique
domain
in two dimensions
and the physical
is given by
domain
and
161
x(5.n)
= xl(rl)U1b)
+ x2(r2)P2W
+ T1 %1dr
113
)U (s)+T2
dy
$(r 2
2
)U (s), 4 (14)
+ ~,(r,)k~,(s) - Tl %r dr
y(S,rl) = ~l(rl)!~lW
113
1~ (s) - T2 $r2)U4W, 2
and
Lll(sI= 2s3 - 3s2 + 1,
1-12(s) = -2s3 + 3s2,
U3(S) = s3 - 2s2
+
St
lJ4(s) = s3 - s2,
where
Xl (rl)'Y1(rl)
position on the first boundary as a function of normalized arc length along the boundary position on the second boundary as a function of normalized arc length along the boundary
dxl -&r 1
dyl 1 ,-(rl) 1 drl
first derivative along the first boundary with respect to normalized arc length along the boundary first derivative along the second boundary with respect to normalized arc length along the boundary normal derivative magnitudes for the respective boundaries
Tl'T2 5 = f,(S)
normalized arc length along the first boundary normalized arc length along the second boundary
r2 = f,(S) s = h[e(n)]
arc length along the grid curves connecting the two boundaries
E.11 0(5<1 o
@ I
coordinates from the computational domain
162 Uniformly discretizing
5
and
rl and given the other quantities described
above a corresponding grid is generated in the physical domain from Wq. (14). For grid
curves
connecting the two boundaries (Fig. 141, their relation-
ship and spacing relative to their neighboring grid curves is based on position, derivatives, and derivative magnitudes at the two boundaries, and the blending functions. the
Given that the blending functions are the same for all grid curves,
spacing between the curves is only a function of boundary information.
The boundary positions and derivatives are a function of normalized arc lengths which are Ln turn written as functions of the computational coordinate is the functions
f,(6)
between grid curves.
and
f,(c)
5.
It
that ultimately control the spaciqg
When there is relatively low slope in these functions
(Fig. 14), there is concentration of
grid curves, and when there is relatively
high slope the grid curves are dispersed
(Fig. 15).
In a similar manner
1
1 91
Fig. 14. Spafisg between neighboring grid curves
),I
Fig. 15. Effect of control functions
163 the grid points along a grid curve are distributed by the blending functions. A control function
h[eUl)]
relating
D
to the normalized arc length
deterknines the final grid point distribution along the grid curve (Fig. 16).
Pig. 16. Control of grid points along grid curves The functions
fl(~),f,t~)
and
are called control functions.
h[e(Vl
They should be single valued, smooth, and have smooth derivatives.
Another
condition is that the functions are defined on the unit square (Fig. 17).
OS1
rs
1
t
090
5
n
5
1,o
Fig. 17. Domain for the definition of control functions The control functions can be analytic functions such as
164 A where
the parameter
first boundary. where
control
In general,
spline
in Reference
SIDE BOUNDARY
the two primary as previously
Another
functions
approach
of grid curves
are restrictive for arbitrary
on the unit square.
near the
relative control
This approach
to
is the is
15.
technique
boundaries.
described,
interpolation.
+,nL)
functions
CONSTRAINTS
The two-boundary
(Fig. lea).
the concentration
analytic
can be applied.
use of smoothing described
controls
K
can be constrained This
and then applying
An example
intersecting
by applying
the technique
the recursive
with one side boundary
The two-boundary
= %+,(5)
by boundaries
is accomplished
,s,K))
technique
is performed
,+,(F;)
,-F&'(f.(5)rgr
formulas
constraint
of transfinite
is demonstrated
with the formulation
,s,(C)) ,h(rl)). to)), :,l(f,(SLg,(O)L
h(Q)) ith constraint
+2
s1
t nt
(f,(S,).f,t5) .h(n))
Fig.
lea. step one in the two-boundary
The second
step is from the transfinite
linear blending
$(5,nn5)
function
= +'W
(Fig.
Fig. 18b. Side boundary constraint
technique interpolation
formulation
18b)
+ (1 - 5) 8;(f1KLf2WW 1
- ~~~~~f,~~,.g,~0,,,~~~f2~5)'42~0~~,h~~~~
1
.
with a
165
GRID GENERATION TOPOLOGY The algebraic grid generation techniques that are presented are defined with the assumption that a uniform rectangular computational domain transforms into a physical domain.
Also, exterior boundaries of the computational
domain transform into boundaries on the physical domain.
Consequently the
topology of the physical domain strongly influences how a grid generation technique is applied.
It is obvious that single six-sided box (computational
domain) or a square in two dimensions is not going to transform into all physical domains.
Further, in certain cases, transformations can only be made
by introducing singularities.
Problems most often arise when there are
closed boundaries and in this section some topological considerations are described. For boundaries in two dimensions, there are two primary types of physical domains that transform from a square computational domain.
They are O-type
domains and C-type domains, and are more commonly referred to as O-grids and C-grids.
Several two-dimensional domains are described schematically in
Figure 19 with corresponding boundary numbers in the computational and physica domains.
Also, multiple computational and physical domains can be coupled
with a common boundary.
A resulting grid as well as grid derivatives should
be continuous across a common boundary.
If there is a discontinuity, special
consideration should be taken in the finite difference procedure for the solution of the equations of motion.
2 computational domain
physical domain
4
3
1
Fig. 19a. Simple O-grid
166
4
6
physical domain
computational domain
6
5
5
1
3
2
I
I Fig. 19b. Simple C-grid
6 physical domain
computational domain
a
11
2
4
13
5
I Fig.19c. Hultiple,body O-grid
8 computational domain
9
I 10
4m_
Fig. 19d. Multiple body C-grid
5
6
3
2
7 1
167 I
I
8 ‘7
9’
computational 12 domain 11
i 10 i 4 6
5
1,2,3
I
I
Fig. lge. Combination domain C-grid
2
I I
computational 3
6
2
6
domain 4j7
8
I I
1
I I
5 .~.__
-__.
10
11
1.
12
9
Fig. lgf. L-shape domain For closed boundaries in three dimensions two suitable types of domains are O-O domains and C-O domains.
The topologies associated with closed
boundary domains in three dimensions is more complex than for two dimensions. A primary reason is that a planar surface from the computational domain will not transform into a closed three-dimensional surface without introducing 6 singularities. Rizzi and Eriksson extensively discuss the problems of three-dimensional closed surface topology and associated singularities and the reader is referred to Reference 6. A final note on topology and singularities is that every effort should be made to place an unavoidable singularity in a region where there is little change in the basic equations of motion.
Near a singularity there are large
changes in the derivatives of the computational coordinates with respect to
168 the physical points
coordinates.
These
can lead to inaccuracy
ference
procedure
large changes
and/or possibly
for the solution
between
neighboring
instability
of the equations
grid
in a finite dif-
of motion.
GRID COMPUTATION Algebraic data.
grid generation
These
data."
data can be divided
Fixed data describes
variable
data describes
control
functions.
cursor
Variable
grid derivatives
tory grid characteristics The algebraic ment because given
rate between because grid.
the computation What
the computer
An aside point
graphics
data which
whereas
surfaces
is ideally
control
and "variable
surfaces,
control
of input
suited
of for this
a grid can be modified
terminal,
the resulting
and the data again modified
grid and
until
satisfac-
are achieved. methods
is explicit is necessary,
and graphics
of the large number
work well
can be worked
and the resulting
control
in an interactive
with no iteration however,
terminal
of line segments
for the algebraic
the grid characteristics grid points)
computer
grid generation
grid solution.
"fixed data"
such as internal
observed,
a large quantity
such as bounding
input at a graphics
visually
require
into two types:
quantities
or numeric
generally
quantities
Interactive
type of application. with
methods
necessary
be high
(9600 baud or greater)
that are discussed
out on a relative directly
for a
is that the communication
that must be displayed
methods
environ-
applied
for a is that
small grid to compute
(few a larger
grid.
DISCUSSION
AND CONCLUSIONS
The three algebraic basically Blending
functions,
geometric
procedures
govern
to where
singularities
singularities
procedure
along with physical
of a uniform
the most
of the physical
there are closed boundaries,
in the finite difference
functions
The methods
are
characteristics.
rectangular
are relatively
extensive
simple
computational
of generality.
procedure,
is the topology
be taken relative
grid.
that are described
common
and do not require
and they have a high degree
grid generation where
or control
they are explicit
Next to the computational
several
the transformation
grid into a physical
to understand,
techniques
with
interpolants,
constraints
computational
effort,
grid generation
interpolation
exist
important domain.
consideration
are introduced.
Care should
and the corresponding
for the solution
for
For three dimensions
of the equations
effect of
169 Grid topology
motion.
importance Another means
is not extensively
consideration
is that grid derivatives
that each step in a grid generation Wiqqles
result.
discussed
in this paper,
but its
is emphasized.
is one step propagate
method
must be smooth. must produce
This
a smooth
into the next step and finally
into
the grid. A final point
is that the use of interactive
generation
is highly
generation
is truely
control
is adaptive,
advantageous.
and require
interactive applications
software
is an important
graphics
for grid
when grid
and the grid
in the control
grid generation
few computations,
This implies
of motion
intervention
Since algebraic
relatively
environment.
human
computer
the time is reached
coupled with the equations instantaneous
in the next best approach. explicit
Until
methods
process are
they work very well in an
that the development aspect of algebraic
of computer
grid generation.
REFERENCES 1.
2. 3. 4.
5. 6.
7.
8.
9.
10. 11. 12. 13.
Smith, R. E., Two-Boundary Grid Generation for the Solution of the ThreeDimensional Navier-Stokes Equations, Ph.D. Dissertation, Old Dominion University, May 1981. Smith, R. E., "Two-Boundary Grid Generation for the Solution of the ThreeDimensional Navier-Stokes Equations," NASA TM-83123, May 1981. Dugundji, James, Topology, Allyn and Bacon Inc., 1968. Gordon, W. J. and Hall, C. A., "Construction of Curvilinear Coordinate Systems and Application to Mesh Generation," International Journal for Numerical Methods in Engineering, Vol. 7, 461-477, 1973. Coons, S. A., "Surfaces for Computer Aided Design of Space Forms," MAC TR-41, Contract No. AF-33(6000-42859) MIT, June 1967. Eriksson, Lars-Erik, "Three-Dimensional Spline-Generated Coordinate Transformations for Grids Around Wing-Body Configurations," Numerical Grid Generation Techniques, NASA CP 2166, 1980. Rizzi, A. and Eriksson, L. E., "Transfinite Mesh Generation and Damped Euler Equation Algorithm for Transonic Flow Around Wing-Body Configurations,' AIAA 5th Computational Fluid Dynamics Conference, June 1981, Palo Alto, California. Grid Generation Anderson, P. G. and Spradley, L. W., "Finite-Difference by Multivariate Blending Function Interpolation," Numerical Grid Generation Techniques, NASA CP 2166, 1980. Spradley, L. W. Stalnaker, J. F., and Ratliff, A. W., "Computation of Three-Dimensional Viscous Flows with the Navier-Stokes Equations," AIAA Paper 80-1348, AIAA 13th Fluid and Plasma Dynamics Conference, Snowmass, Colorado. Eisemen, P. R., "A Multi-Surface Method of Coordinate Generation," Journal of Computational Physics, Vol. 33, No. 1, October 1979. Eiseman. P. R.. "Geometric Methods in Computational Fluid Dynamics," ICASE Report No. 81-11, April 1980. _ Eiseman, P. R., "A Coordinate System for a Viscous Transonic Cascade Analysis," Journal of Computational Physics, Vol. 29, 1978. Coordinates About Wings," Fourth AIAA Eiseman, P. R., "Three-Dimensional Computer Fluid Dynamics Conference, Williamsburg, Virginia, July 1979.
170
14. Smith, R. E. and W&gel, B. L., "Analytic and Approximate Boundary-Fitted Coordinate Systems for Fluid Flow Simulation," AIAA Paper 80-0192, AIAA 18th Aerospace Sciences Meeting, January 1980, Pasadena, California. 15. Kowalski, E. J., "Boundary-Fitted Coordinate Systems for Arbitrary Computational Regions," Numerical Grid Generation Techniques, NASA CP 2166, 1980. 16. Eiseman, P. R. and Smith, R. E., "Mesh Generation Using Algebraic Techniques," Numerical Grid Generation Techniques, NASA CP 2166, 1980. 17. Steger, J. L., Implicit Finite Difference Simulation of Flow About Arbitrary Geometries with Applications to Airfoils," AIAA Paper 77-665, AIAA 10th Fluid and Plasma Dynamics Conference, Albuquerque, New Mexico, June 1977. 18. Iiindman, R. G., "Geometrically Induced Errors and their Relationship to the Form of the Governing Equations and the Treatment of Generalized Mappings," AIAA Paper 81-1008, AIAA 5th Computational Fluid Dynamics Conference, Palo Alto, California, June 1981. 19. Smith, R. E., Kndlinski, R. A., and Everton, E. L., "A Grid Spacing Control Technique for Algebraic Grid Generation Methods," AIAA Paper 82-0226, AIAA 20th Aerospace Sciences Meeting, January 1982, Orlando, Florida.