rather than a mathematician and whose main interest is in applications. With the advent of parallel computers. this is critical to understanding a general 2D and 3D multi-block grid generator. Much emphasis is given to algorithms and data structures.R. New York. the International Space Course at the Technical University of Munich (1993). In addition. as will be outlined in the chapter on grid generation on parallel computers. The best way to understand theory is to develop tools and to generate examples. However. PDC. substantial progess has been achieved in automatic blocking (GridPro package). namely the Grid and GridPro [44] package. Computational Fluid Dynamics. for the International Space Course at the Technical University of Stuttgart 1995). White Pains. scalability. which are presented from a real working code. much more complex problems can be solved. Rome (1994).de is being built where further information will be available. This material is presented here for the Þrst time. Brussels (1992).Concepts of Grid Generation
1
1..e. and is missing from other presentations on the subject. a web address www. No claim is made that the packages presented here are the only or the best implementation of the concepts presented. presenting the key issues of the parallelizing strategy. and at the Von Karman Institute. envisaged for the second half of 1997. Material that has been prepared for presentations at the California Institute of Technology and the NASA Ames Research Center has also been partly included. The recent Workshops hely by the Aerothermodynamics Section of the European Space Agency’s Technology Center (ESTEC) (1994. entitled Grid Generation. a CD-Rom containing a collection of 2D and 3D grids can be obtained (handling fee) from the authors. The emphasis is on how to use and write a grid generation package. and communication.cle. Although the emphasis is on grid generation in Aerodynamics. 1995) has proved that grid generation is one of the pacing items in CFD. Multi-blocks may well be the key to achieve that quantum leap. In these notes Numerical Grid Generation is presented in a way that is understandable to the scientist and engineer who is more a programmer. Hence. To keep the number of pages for this article at a manageable size. However. at the Istituto per Applicacione del Calculo (IAC). domain decomposition. which offer the promise of a quantum leap in computing power for Computational Fluid Dynamics (CFD). These notes are the continuation of courses in grid generation that were given at Mississippi State University (1996). i. in cooperation with P. Some recent publications can be found on the World Wide Web under . The reader interested in these aspects should contact the authors. A discussion of the dependence of convergence speed of an implicit ßow solver on blocknumber for parallel systems wil also be given.
. The multi-block concept gets an additional motivation with the advent of parallel and distributed computing systems. the introduction to parallel computing and the section on interfacing multiblock grids to ßow solvers had to be omitted.. a certain amount of low level information is presented. and Parallel Computing. built by the authors. Furthermore. provided computational grids of sufÞcient quality can be generated. Eiseman.0 Overview and Status of Modern Grid Generation
These lecture notes are the less mathematical version and an excerpt from a forthcoming book being currently written by the Þrst author. the codes presented can be used for any application that can be modeled by block-structured grids.

In 3D. modeling a somewhat simpliÞed Orbiter geometry. The object–oriented grid generation approach is demonstrated in 2D with a 284 block grid. in general. Next. higher order solvers need overlaps of 2 points in each direction. The software to generate this overlap for multi-block grids is also discussed. explaining the 3D mapping procedure. First. the grid construction process is explained and demonstrated for several 2D examples. and a mono and a 7-block grid are generated. the questions of grid quality and grid adaptation are addressed. The last vehicle that is presented is the Space Shuttle Orbiter where surface grid generation is given for a 4-block grid. using both the Grid and GridPro packages. Along with the examples. the ”Hyperboloid Flare in the F4 windtunnel”. using the tools provided by Grid . In the next step. a multi-block grid for the double ellipsoid is generated to demonstrate the removal of singularities. object–oriented grid generation in exemplied by the Electre body in the F4 windtunnel. we start with a simple mono-block grid for the double ellipsoid. it is shown how the 94-block Euler grid is converted to a 147-block Navier-Stokes grid. grids generated are slope continuous.2 Since grid generation is a means to solve problems in CFD and related Þelds.
. In the last step. After that. namely a Shuttle without body ßap. a chapter on interfacing the Þnal grid to the Euler or Navier-Stokes solver is provided. viz. In 3D. Although. The rationale for the choice of that topology is carefully explained. surface grid generation is demonstrated for the Hermes Space Plane. the surface grid and volume grid for a 94-block Euler mesh of a Shuttle with body ßap are presented.

Covariant and contravariant vectors are related by the metric coefÞcients (see
. y and z. It is assumed that a one-to-one mapping A M exists with x´ξ1 ξm µ :
¼ x1´ξ1
xn
ξm µ ξm
½
´
ξ1
µ
(2.0
Methods of Differential Geometry in Numerical Grid Generation
In the following a derivation of the most widely used formulas for nonorthogonal curvilinear 1 xn coordinate systems ( Fig.Concepts of Differential Geometry
3
2.2)
The tangent vector ek points in the direction of the respective coordinate line.1: Curvilinear nonorthogonal Coordinate System. 2. The contravariant base vectors e point i ˆ in the tangential directions but are not unit vectors.1) is given as needed in numerical grid generation. whereas in the transformed space θ and ψ form a rectangle. The coordinates x1 x2 and x3 correspond to the usual Cartesian coordinates x.3)
The ei are called contravariant base vectors and are orthogonal to the respective covariant vectors for i j.1)
As an example. we consider the surface of a sphere where Rm R2 and ξ1 θ π 2 θ π 2. ξ2 ψ 0 ψ 2π. These base vectors are called covariant base vectors. Vectors ei are Cartesian unit vectors. k ∂ξk 1´1µm
ek :
(2. Variables x 1 m are arbitrary curvilinear coordinates. In the physical space. we have the surface of a sphere. A second set of base vectors is deÞned by ei ¡ e j δij
(2.
e3 e2 e1 e3 e1 e2
Figure 2. The denote Cartesian coordinates and variables ξ ξ approach taken follows [9]. The tangent vectors or base vectors at a point P ¾ M are deÞned by ∂x .

6)
In order to measure the distance between neighboring points. contravariant components of a vector are found by parallel projections onto the axes. According to its transformation behavior. a vector is called contravariant or covariant. Let vibe the i components of a vector in the coordinate system described by the x an let vi be the components in i the system ξ (α) A vector is a contravariant vector if its components transform in the same way as the coordinate differentials: ∂xi dxi ∂ξi dξ j (2.9)
In Cartesian coordinates there is no difference between covariant and contravariant components since there is no difference between covariant and contravariant base vectors. A physical vector can either be represented by contravariant or covariant components vi ei v je j
v
(2. the Þrst fundamental form is introduced:
gi j :
ei ¡ e j
(2. whereas covariant components ar obtained by orthogonal projection. The components of the inverse matrix are found from gi j gik δk j
(2. In two dimensions.7)
The gi j are also called components of the metric tensor.4)
where the summation convention is employed. For the above example the two tangent vectors eθ and eψ are obtained by differentiating each of the functions x´θ ψµ y´θ ψµ and z´θ ψµ with respect to either θ or ψ.5) ∂xi vi ∂ξ j v j (β) A vector is covariant if it transforms in the same way as the gradient of a scalar function: φ´x1
∂φ ∂xi
xn µ
∂φ ∂ξ j ∂ξ j ∂φ ∂ξ j ∂xi ∂xi ∂ξ j ∂ξ j vi ∂xi v j
φ´ξ1
ξm µ (2.4 below).8)
The distance ds between two neighboring points is given by
Ô
ds
´gi j dξi dξ j µ
(2. Therefore the
.

into equation (2. one obtains the computational useful form ∂ξi ∂2 xl ∂xl ∂ξ j ∂ξk
Γijk
(2. not tensor components. One can therefore write Γkj dξ j ek i
dei
(2. They are. Taking the scalar product with e .10)
In order to Þnd the transformation rules of derivatives of scalars. it can be represented by the system of base vectors ek . along with equation (2. the Christoffel symbols of the Þrst and second kinds are introduced. which follows directly from their transformation behaviour.15)
. The relationship between the Γkj and the glm is found in the following way. Insertion of the deÞnii tion of the metric components glm in equation (2.12). the Christoffels symbols vanish. Further.7).2). Since dei is a vector. These symbols are also called i k Christoffels symbols of the second kind.12)
where a comma denotes partial differentiation.12) and interchanging indices leads to 1 ij g 2 ∂g jl ∂ξk
Γijk
· ∂gklj ∂gklj ∂ξ ∂ξ
(2. This changes the base vector ei by dei . dei is proportional to dξ j . one obtains directly from equation 2.11)
where the symbols Γkj are only coefÞcients of proportionality. the deÞnition of the tangent vector.Concepts of Differential Geometry
5
matrix gi j (the ˆ denotes the Cartesian system) is the unit matrix.13)
If the base vectors are independent of position. ˆ The components of gi j in any other coordinate system can be directly calculated using the chain rule: ∂xk ∂xl ˆ ˆ gkl ˆ ∂xi ∂x j
gi j
(2. Suppose that coordinate ξj is changed by an amount of dξ j . vectors and tensors. however.14)
If we use equation (2.11
Γkj i
ei j ¡ ek
(2. The Christoffels symbols of the Þrst kind are deÞned as gil Γljk
Γi jk :
(2.

6 The knowledge of the grid point distribution.17)
In curved space. From equation.21)
Raising the indices in equation (2.20)
where J is the Jacobian of the transformations. Forming the determinant from the components of equation (2.14) yields 1 im ∂gmi g 2 ∂ξ j
Ô 1 ∂ g Ô g ∂ξ j
Γiji
´lnÔgµ j
(2. the nabla operator ∂ ∂ξi
∇:
ei
(2. For the calculation of the cross product the Levi-Civita tensor is introduced. and therefore
Ô gεiˆjk
εi jk
(2. allows the numerical calculation of the Christoffels symbols.16)
Ô where g is the determinant of the metric tensor.10). then.18)
is used. equation (2.22)
. i.21) with the metric. results in g J2 .
ˆ εi jk :
´1
1 0
for i 1 j 2 k 3 and all even permutations for all odd permutations if any two indices are the same
(2. In two dimensions with curvilinear coordinates ξ η and Cartesian coordinates x y. and hence ∂xl ∂xm ∂xn ˆ ˆ ˆ ˆ εlmn ∂xi ∂x j ∂xk
εi jk
ˆ J εi jk
(2. In the following. one Þnds
ξx ξy ηx ηy Ô g J
J 1
yη xη yξ xξ ´xξyη yξxηµ
(2. one obtains εi jk
Ô i jk ˆ gε
(2.e. that is g is the Jacobian of the transformation. upper and lower indices are the same and are summed over.0) we know the transformation law for covariant components. partial differentiation is replaced by covariant differentiation which takes into account the fact that base vectors themselves have non-vanishing spatial derivatives. If we contract the Christoffels symbols.(2.19)
where ˆ again denotes the Cartesian coordinate system.

In the past. A comparison between the impact of the numerical technique and the computational grid used. and making it possible to use the efÞcient techniques of Þnite differences and Þnite volumes for complex geometries. followed by quasi-automatic grid generation (via a special language) and Þnally produces and visualizes the ßow solution. Data structures and ObjectOriented-Programming (OOP) are not taught in engineering courses and thus their importance is not recognized in this Þeld.0. and Motif as well as on Open GL. The advent of high-speed computers has made ßow computations past 3D aircraft and space-
. reveals that in many cases grid effects are the dominant factor on the accuracy of the ßow solution. the majority of the engineering software is still in Fortran. coined PAW (Parallel Aerospace Workbench). For parallelization PVM (Parallel Virtual Machine) or MPI (Message Passing Interface) should be used (or a similar package).2 A Short History on Grid Generation
Numerical grid generation has the dual distinction of being the youngest science in the area of numerical simulation and one of the most interesting Þelds of numerical research. Tcl/Tk. provides modules for surface repair. the number of researchers working in grid generation is at least an order of magnitude smaller than the number of scientists active in CFD. With the advent of massively parallel systems and high end graphics workstations a new level of performance in aerodynamic simulation can be achieved by the development of an integrated package. ANSI C X11. and Mastin. First. Winslow [1]. Strange enough. Despite this fact. of Computational Physics [2].Numerical Grid Generation
9
3. Hence. the latest hardware is used together with the software concepts of the late Þfties. generating a hardware independent parallel ßow solver. partly interactive package that. Although conformal mappings were used in aerospace for airfoils. a modiÞcation of a numerical scheme in 1D is sufÞcient to justify a new paper. originating the Þeld of boundary conformed grids. Software development has to be based on standards like UNIX. In 1974. e.
3. a CFD paper is easier to publish since.1 General Remarks on Grid Generation Techniques
CFD of today is marked by the simulation of ßows past complex geometries and/or utilizing complex physics. As a consequence. algorithms and programming have been considered to be outside the engineering domain. authored by Thompson. many codes in industry are not state of the art. a general methodology for irregular 2D grids was Þrst presented in the work of A. publishing a grid generation paper involves a large amount of programming (time consuming) and algorithm development as well as computer graphics to visualize the grid.0. Second. Thames. PAW would provide the aerospace engineer with a comprehensive. accepting the CAD geometry from the design ofÞce.0 Introduction and Overview of Numerical Grid Generation Techniques
3. The reason for this is most likely that publishing a paper in CFD is much more rewarding. a paper entitled ”Automatic numerical generation of body Þtted curvilinear coordinate systems for Þelds containing any number of arbitrary two-dimensional bodies” appeared in the J. This paper can be considered a landmark paper.g.

with emphasis on the actual grid generation process as applied to complete geometries. the majority of structured grid generation codes now employs the multiblock concept. with the advent of massively parallel systems. not attainable by wind tunnel experiments. which is unstructured on the block level but results in a structured grid within a block. Þnally begin to realize that a new way of thinking is needed in software engineering.g. The interfacing to the ßow solver is done automatically and results can be visualized directly by the Plot3D package from NASA Ames or by any other package that supports this format. and Open Gl. e. generally coming from the design ofÞce in a CAD format. which will allow realistic simulations of turbulent ßows with real gas effects (high temperature) past complex 3D conÞgurations. To determine the block topology of a complex SD (solution domain) requires a certain effort. Only with this approach it can be hoped to reduce the size and complexity of source code. to write and maintain very large software packages and to implement them on arbitrary computer architectures. The purpose of these notes is to present a clear. vivid. or after a few iterations the solution may blow up. X-Windows.1 the complete modelization process to obtain a ßow solution is shown. The chart then shows the grid generation process. used to Fortran since the late Þfties. the development of generally applicable and easy to use 3D grid generation codes still presents a major challenge.1. This information has to be processed to convert it into a format that can be read by the modules of Grid . On the numerical side. a revolution in scientiÞc programming is taking place. It is assumed that a surface description of the vehicle is available. a new era of CFD (Computational Fluid Dynamics) has started. characterized by the concept of object-oriented-programming.??) is a major step forward in automatic topology deÞnition. With the advent of modern programming languages like C and C++. where each module generates output that can be read by the following modules. either because there is no convergence of the ßow solution. In practice. The stiffness of the equations
. No longer are scientists and engineers restricted to wind tunnel and free ßight experiments. problems may be encountered with the grid generated and the ßow solution may not be obtained.g. Engineers. Grid generation packages like GridPro. Most important are tools for data conversion to directly use CAD data and to automatically interface a grid to the ßow solver.. CATIA. Although nearly two decades have past since the work of Thompson et al. Grid generation has to embedded in the overall solution process that is depicted in Fig. The recent AZ–manager code (see Sec. This depends very much on the ßow physics and to a large extent on the experience of the user. the widely known Tecplot software from AMTEC Engineering. based on boundary conforming grids.3. is one of the important results of these calculations. A clear understanding of many of the aspects of subsonic to hypersonic ßows. Special emphasis is placed on algorithms and examples will be presented along with a discussion of data structures and advanced programming techniques. comprehensive treatment of modern numerical grid generation. These codes do not yet have reached a stage where they can be used as a black box. producing some kind of pipeline.10 craft as well as automobiles a reality. allowing full portability. In addition. e. 3. Therefore. In the case of structured grids it was soon recognized that branch cuts and slits or slabs did not offer sufÞcient geometrical ßexibility. and Grid code [20] are based exclusively on C. to fully use available hardware. topics nearly almost always neglected in the scientiÞc and engineering Þeld. A case with Ma of 25 and thermo-chemical nonequilibrium Navier-Stokes equations is of course much harder to solve than an Euler problem at Ma 2. Various conversion modules are available. In the chart of Fig. the development of high-speed computers as well as graphics workstations has allowed aerodynamicists to perform intricate calculations that were unthinkable a few decades ago.

some results for Euler and N-S calculations are presented.12 depends on the physics and also on the grid. A Navier-Stokes grid for a 3D vehicle with a cell aspect ratio of 106 in the boundary layer is not only much more difÞcult to generate but it is also much harder for the numerical scheme to converge to the physical solution in an efÞcient way.
. Although the emphasis of these lecture notes is on grid generation.

there may be hundreds of bodies in the ßowÞeld. or. An O-type topology for a 4 element airfoil may be adequate for Euler solutions.3 What Is a Good Grid
First. the discussion will be restricted to multi–block grids. It is felt that any kind of geometry can be gridded by this approach. The large number of blocks is needed because of complex topology. the answer is straightforward. However.
. and hence the accuracy of the results may deteriorate. Grids having a singularity. need grid line orthogonality at Þxed boundaries. necessitating a C-H type topology. namely that it is not possible to deciding by simply looking at a grid whether it is good or not. since no information can cross the singularity. for example. otherwise the physics cannot be adequately represented. but may be insufÞcient to resolve the wake resulting from viscous ßow. Their effect is that they reduce the order of the numerical scheme. This is the case for some turbulence models to determine the distance from the wall. the wrong drag or skin friction.g. However.g. because the solver is run on a MIMD (Multiple Instruction Multiple Data) massively parallel system that may comprise several thousand processors. Grid topology also has an inßuence on the physical results. However. resulting in a superior grid quality when compared with unstructured grids. or C2 ) º nonalignment of the grid with the ßow º insufÞcient resolution to resolve proper physical length scales º grid topology not well suited to sufÞciently cover the ßow physics º grid lacks special features needed by physical submodels º grid is not singularity free (e. The other topics are strongly coupled with the ßow physics. at the nose part of an aircraft. resulting in. the magnitude of these effects will depend on the underlying ßow physics. e. at stagnation point)
The Þrst three points are mostly independent of the ßow physics. may lead to substantially increased computing time because of a slowdown in convergence. It is known from numerical error investigations that locally 1D ßow (proper grid alignmnet) will give an improved numerical accuracy.g. in particular for viscous ßows. Physical submodels may.0. Second. On the block level the grid is completely unstructured. grid lines conforming to the boundary of a body will produce better accuracy. C1 . there exist several criteria that tell when grid quality is not good. e. for instance.
º high degree of skewness º abrupt changes in grid spacing º insufÞcient grid line continuity (C0 . boundary layer). Length scales have to be properly resolved (e.g. it should be kept in mind that grids of tens of thousands of blocks have to be generated (automatically). for example.Numerical Grid Generation
13
3.

with a monoblock mesh gridline topology is Þxed. In early grid generation it was attempted to always map the physical solution domain (SD) to a single rectangle or a single box in the computational domain (CD). A 2D grid around a circle which is mapped on a single rectangle necessarily has O-type topology. which provides the necessary geometrical ßexibility and the computational efÞciency for the Þnite volume or Þnite difference techniques. The fourblock grid shown in Sec. In addition. The larger number of blocks is needed to get a special grid line topology. Slope continuity across neighboring block boundaries is achieved by overlapping edges (2D) or faces (3D). unless additional slits (double valued line or surface) or slabs (blocks are cut out of the SD) are introduced. Hence. that is a code is needing only 2 or 3 ”for” loops (C language). A 94 block Euler grid for the Shuttle has been generated. The number of gridpoints is reduced substantially. the Joukowski airfoil proÞle. e. and some degree of unstructuredness has to be introduced. it became soon obvious that certain grid line conÞgurations could not be obtained. From differential geometry the concept of an atlas consisting of a number of charts is known. and to better represent the ßow physics (see discussion in Sec. For reasons of geometrical ßexibility it is mandatory that each block has its own local coordinate system. The main advantage of the structured approach. a singularity invariably leads to a clustering of grid points in an area where they are not needed. it would be advantageous if the grid line distribution would be similar to the streamline pattern. The overlap feature facilitates the construction of the ßow solver substantially and allows the direct parallelization of the code on massively parallel
. This line needs special treatment in the ßow solution. The set of charts covers the atlas where charts may be overlapping. has been lost. It has been observed that convergence rate is reduced. This directly leads to the multiblock concept.0. The grid is smooth across block boundaries. This always leads to a singular line. Therefore.14
3. The solution domain is subdivided into a set of blocks or segments (in the following the words block and segment are used interchangeably). which normally occurs in the nose region.4 Aspects of Multiblock Grid Generation
Structured grids use general curvilinear coordinates to produce a body Þtted mesh. Additional requirements with regard to grid uniformity and orthogonality cannot be matched. For multiply connected SDs. In addition. However.2. now the connectivity of the charts has to be determined. although special numerical schemes have been devised to alleviate this problem. with a small enough Re number.g. 3. Hence blocks can be rotated with respect to each other (Sec. The ßow solver ParNSS needs an overlap of 2 faces. For grid generation an overlap of exactly 1 row or 1 column is needed (2D). which includes the body ßap. computing time may be increased substantially. Each chart is mapped onto a single rectangle. One can start with a simple monoblock grid that wraps around the vehicle in an O-type fashion. If one considers for example the 2D ßow past an inÞnitely long cylinder. This has the advantage that boundaries can be exactly described and hence BCs can be accurately modeled. Since multiblock grids are unstructured on the block level. namely that one only has to deal with a rectangle or a box. information about block connectivity is needed along with the of each block. a procedure well known in complex function theory and analytic mapping of 2D domains. Furthermore. The conclusion is that this amount of structuredness is too rigid. no special discretization in the ßow solver across boundaries is needed. branch cuts had to be introduced.3). For a vehicle like the Space Shuttle a variety of grids can be constructed (see Sec. ??).??). ?? removes the singularity but otherwise retains the Otype structure.

one has to discern between physical boundary points on Þxed surfaces that can be used for computation of the right hand side of the Poisson equation and matching boundary points on overlap surfaces connecting neighboring blocks. Therefore communication among blocks follows an irregular pattern. e. However.Numerical Grid Generation
systems using message passing. this does not mean that the solution domain in the computational plane has a regular structure.
15
Each (curvilinear) block in the physical plane is mapped onto a Cartesian block in the computational plane (CP). That means that Þrst the right hand side is determined for each grid point on a surface and then the values of the control functions are interpolated into the interior of the solution domain. However. All grid point positions on the faces of a block must be known before the Poisson equations for this block can be solved to determine the grid point positions of the interior points.g. In this context a grid point is denoted as boundary point if it lies on one of the six faces of a block in the CP.
. The positions of the latter ones are not known a priori but are determined in the solution process (see above). one for each coordinate direction. A coordinate transformation of the governing physical equations and their respective boundary conditions from the physical plane to the computational plane is also required. with some additional smoothing. For the parallelization the important point is that there is no nearest neighbor relation for the blocks. simply because of the communication overhead. rather it may look fairly fragmented. The right hand side of the Poisson equations is used for grid point control and is determined from the speciÞed grid point distribution on the surfaces. lattice gauge theory problems. The actual solution domain on which the governing physical equations are solved is therefore a set of connected. A parallel architecture that is based on nearest neighbor communication. caused by random communication. The grid point distribution within each block is generated by the solution of a set of three Poisson equations. will not perform well for complex aerodynamic applications. regular blocks in the computational plane.

Since there are no coordinate lines in the UG. In computational ßuid dynamics (CFD) in general and in high speed ßows in particular many ßow situations are encountered where the ßow in the vicinity of the body is aligned with the surface. we wish to explain the motives leading to this approach. factoring in the direction of the plus and minus eigenvalues of the Jacobians) is employed.e.1 Elliptic Equations for 2D Grid Generation
In the following the elliptic PDEs used in the grid generation process will be derived. This is due to the more complicated data structure needed for UGs. resulting in two factors if LU decomposition (that is. The resolution of the BL leads to large anisotropies in the length scales in the directions along and off the body. In addition. In the past. the ßow solver based on the UG is substantially slower than for the SG. by the magnitude of the chemical production terms.1 Equations of Numerical Grid Generation
3. such a grid has to have coordinate lines. and by some authors of up to 10. Here some type of SG is indispensable. i. numerical diffusion can be reduced. better accuracy is achieved. orthogonality will increase the accuracy when algebraic turbulence models are employed. the use of the so-called thin layer approach. to describe the surface of the body a structured approach is better suited. reduces the computer time by about 30 %. Moreover. it cannot be completely unstructured. using the mathematical tools of the previous section. Thompson et al. A fairly complex procedure would be needed to artiÞcially construct these lines. the BL (Boundary Layer) must be resolved. The use of a SG. i. that is retaining the viscous terms only in the direction off the body. which are connected. allows the alignment of the grid in that direction.g. Clearly.( see e. This demands that the grid is closely wrapped around the body to describe the physics of the BL (some 32 layers are used in general for SGs or UGs).e. A BFG exactly matches curved boundaries. [2]) introduced the concept of a boundary Þtted grid (BFG) or structureg grid (SG). In order to invert the implicit operator. implicit schemes will be of advantage.1. and for complex SDs generally will consist by a set of so called blocks. First. that is a grid whose grid lines are aligned with the contours of the body. Furthermore. or in three factors if the coordinate directions are used. The latter have been mainly used in structural mechanics. have been given in the literature. extremely small time steps will be necessary. This is especially true in the case of hypersonic ßow because of the high kinetic energy. the difference in surface temperature for S¨ nger at cruising speed (around a Ma 5) is about 500 K depending on laminar or turbulent ßow calculations. factorization is generally used. Hence. Since the time-step size in an explicit scheme is governed by the smallest length scale or. for example. For the unstructured approach there is no direct way to perform this type of factorization. Thus. facilitating the implementation of BCs (Boundary Condition) and also increasing the numerical accuracy at boundaries. i. SGs can be made orthogonal at boundaries. About two decades ago. in the case of chemical reacting ßow. a SD may be covered by a set of hundreds or even thousands of blocks. In order to calculate heat loads for vehicles ßying below Ma 8. This behavior is not demanded by accuracy but to retain the stability of the scheme. In the approach taken by the authors. Factors of 3.16
3. there is a prevailing ßow direction.e. resulting in locally quasi 1D ßow. In the solution of the N-S (Navier-Stokes) equations. turbulence models have to be used. Depending on the real
. Second. Moreover. this simpliÞcation is not possible. Þnite differences and Þnite elements have been used extensively to solve computational problems.

which can be perfectly matched by adaptive grid alignment. coupled to the solution process. but within each block the grid is structured and slope continuity is provided across block boundaries (on the coarsest level). which has been employed in the Grid [27] and in the GridPro[35] packages. the alignment of the grid can result in a more accurate solution than randomly Þlling the space with an enormously large number of smaller and smaller tetrahedrons or hexahedrons. however. The highest degree of freedom of course is obtained in UGs. In many cases of practical interest. [40]. Therefore turbulent calculations are of high importance.It is therefore felt. that is the blocks itself can be considered as Þnite elements. the automatic zoning feature of GridProhas been used to generate a grid of several thousand blocks (see Sec. In addition. The majority of physical phenomena encountered in external and also in internal ßows exhibit certain well ordered structures.
∆ξ
P. For the derivation of the transformed grid generation equations one starts from the original Poisson equations and then the role of the dependent and independent variables is interchanged. A comparison of these two approaches is given by Dannenhoffer [2] where local enrichment gives somewhat better results. Recently a 94 block grid (Euler) and a 147 block grid (N-S) for the Space Shuttle have been generated. and the fact that ξ and η are coordinates themselves. Moreover. ??). Only if a very complex wave pattern evolves due to special physical phenomena. the coupling of SGs with UGs is possible as has been shown by Weatherill and Shaw et al. generating dozens of shock waves coming from an explosion. If a block is reÞned locally. ∆η
Q
(3. In general. we Þnd
j gik ´ξ i k Γik ξ
∆ξ
j
µ
gik ei k ¡ e1
(3.2)
. it is much more costly to use. A mixture between the boundary Þtted grid approach and the completely unstructured approach is the use of multiblock grids. that mesh redistribution and alignment is totally adequate for the major part of the ßow situations encountered in CFD.1)
Using (2. such as a bow shock or system of reßected shocks or some type of shock-shock interaction. this feature cannot be maintained. the transformation equation of the Laplacian. On the block level the grid is unstructured. Such a grid is called a hybrid grid. demonstrating the variability of this approach.Numerical Grid Generation
17
surface temperature encountered in ßight a totally different type of vehicle has to be designed since a cooled Titanium wall cannot withstand a temperature of 1300 K for a long period (about 20 minutes). the UG seems to be advantageous.34). for example Þne resolution of a bow shock or a canopy shock and in situations where shocks are reßected. However. for example. especially in aerodynamics. SG provide sophisticated means both for clustering and adaptation using redistribution or local enrichment techniques. Only a SG can provide the alignment along with the orthogonality at the boundary to accurately perform these calculations.

i. The derivation for the surface grid generation equations is similar as in the 2D or 3D case.e. ´u vµ values are the boundary points in the ´ξ ηµ grids in the computational plane (CP).11 need to be interchanged. To obtain the familiar form of the equations. the dependent and independent variables in Eqs. Then one iteration step is performed using these boundary data.18. To obtain the new coordinates Poisson equations are used. the surface grid generation equations take the form A curved surface can be deÞned by the triple
j
∆ξ ∆η
gik ξ i k Γik ξ gik
j
P Q
j η i k Γik η j
(3. Overrelaxation is
xnew
xold · ω´x xold µ 1
ω
2
(3. until a certain number of iterations has been performed or until the change in the solution of two successive iterations is smaller than a speciÞed bound. After that. z where u and v parametrize the surface. Using the ∇ operator of equation 2. i. Therefore.2 Elliptic Equations for Surface Grid Generation
x y ´u vµ. those newly iterated points which are part of the overlap are used to update the boundary points of neighboring blocks and the whole cycle starts again. Introducing the matrix M of the contravariant row vectors from the transformation ´u vµ ´ξ ηµ.12)
. 3. i. which is a mapping from R2 R3 . the metric is given by the transformation of the parameter space ´u vµ to physical space ´x y zµ.9)
where the notation J2 g. 2βi achieved by computing the new values from
g12 .1.e. αi j g22 .
3. and γi
j
g11 was used. ξu ηu ξv ηv
M
(3. iterating on all interior points.Numerical Grid Generation
1 2´αi
19 αi j ´xi·1 j xi 1 j µ βi j ´xi·1 j·1 xi 1 j·1 xi·1 j 1 · xi 1 j 1µ· γi j ´xi j·1 xi j 1 µ· 1 2Ji2 j Pi j ´xi·1 j xi 1 j µ· 1 2Ji2 j Qi j ´xi j·1 xi j 1 µ
j
xi
j
j
· γi j µ 1
Ó
(3.e ei ∂i ´x y zµ with ∂1 ∂u and ∂2 ∂v . and not with respect to ξ and η. i.10)
The solution process for a multiblock grid is achieved by updating the boundaries of each block. by receiving the proper data from neighboring blocks and by sending overlapping data to neighboring blocks. a transformation from ´u vµ to ´ξ ηµ.e. DeÞning a new coordinate system on such a surface only means another way of parametrization of the surface.11)
In this case the derivatives are with respect to u and v.

ei ∂i ´ξ η ζµ with ∂1 ∂ξ ∂1 ∂x etc. η and ζ are coordinates themselves.z) to ´ξ η ζµ. and k can be used to indicate the grid point position.38)
where Ψ is a scalar function. In more compact notation the elliptic generation equations read. The elliptic type of Eqs.18. y. z coordinate values describing the surfaces which form the boundary of the physical solution domain. In addition. z are now the dependent variables. Similar as in 2D. To perform the transformation the ∆-operator in general curvilinear coordinates is used: ∆Ψ
j gik Ψ i k Γik Ψ
j
(3. in order to produce an orthogonal grid in the Þrst layer of grid points off the surface. η j . η and ζ. Normally. However. The x. ∆ζ
R
(3.Numerical Grid Generation
23
3.3
Elliptic Equations for 3D Grid Generation
The following general coordinate transformation from the cartesian coordinate system. The co -and contravariant base vectors are given by ei ∂i ´x y zµ . η or ζ form the boundaries and x. The grid in the CP is uniform with grid spacings ∆ξ = ∆η = ∆ζ = 1.
(3. to the CD. Dirichlet BCs are used.(3. this set of equations is not solved on the complex SD. ζ ξ´x y zµ η´x y zµ ζ´x y zµ
(3.
. ξ y´ξ η ζµ. The determinant of the metric tensor.36)
∆ξ
P. The Poisson equations for the grid generation read: ξxx · ξyy · ξzz P ηxx · ηyy · ηzz Q ζxx · ζyy · ζzz R where P Q R are so-called control functions that depend on ξ.35)
Since there is a one-to-one correspondence between grid points of the SD and CD (except for singular points or singular lines). instead it is transformed to the CD. indices i. denoted by coordinates ´x y zµ.2. y. g. but since ξ. j. η.y. denoted by coordinates (ξ. ζk ) in the SD as functions of x y z. is given by the transformation from (x. ∆η
Q. is considered. the equations become nonlinear. Equations 3. η z´ξ η ζµ. a set of Poisson equations is used to determine the positions of (ξi . prescribing the points on the surface.1.37 are solved in the computational plane using Eq. allowing the points to move on the surface. but for adaptation purposes von Neumann BCs are sometimes used. are now used as boundary conditions to solve the transformed equations below.37)
where planes of constant ξ. The one to one transformation is given by (except for a Þnite (small) number of singularities): x y z x´ξ η ζµ. ζ). proper Boundary Conditions (BC) have to be speciÞed.36) is not altered.

A comparison of these two approaches is given by Dannenhoffer [2] where local enrichment gives somewhat better results. the use of the so-called thin layer approach. especially in aerodynamics. The resolution of the BL leads to large anisotropies in the length scales in the directions along and off the body. This demands that the grid is closely wrapped around the body to describe the physics of the BL (some 32 layers are used in general for SGs). i. by the magnitude of the chemical production terms. reduces computer time by about 30 %. i. facilitating the implementation of BCs (Boundary Condition) and also increasing the numerical. In the solution of the N-S (Navier-Stokes) equations. SGs provide sophisticated means both for clustering and adaptation using redistribution or local enrichment techniques.Numerical Grid Generation
25
3. and thus makes it easier to use a sophisticated impicit scheme. However. In order to invert the implicit operator. this simpliÞcation is not directly possible. Boundary Fitted Grids (BFG) have to have coordinate lines. In the present approach. and by some authors of up to 10. however. and for complex SDs will consist of a set of blocks. Moreover. Here some type of SG is indispensable. that mesh redistribution and alignment is totally adequate for the major part of the ßow situations encountered in external ßows. In general. In CFD in general. many situations are encountered for which the ßow in the vicinity of the body is aligned with the surface. and in high speed ßows in particular. Factors of 3. In many cases of practical interest.e.2
Grid Generation Concepts
3.e. to describe the surface of the body a structured approach is better suited. Moreover. Second. This is due to the more complicated data structure needed for UGs. i. The highest degree of freedom of course is obtained in UGs.2. it is much more costly to use. there is a prevailing ßow direction. This is especially true in the case of hypersonic ßow because of the high kinetic energy. the BL (Boundary Layer) must be resolved. in the case of chemical reacting ßow. allows the alignment of the grid.e. Since there are no coordinate lines in the UG. for example Þne resolution of a bow shock or a canopy shock and in situations where shocks are reßected. Hence. resulting in locally 1D ßow. An important point for the accuracy of the solution is the capability of grid point clustering and solution adaptation. Thus. but to retain the stability of the scheme.
.1 Computational Aspects of Multiblock Grids
As has been discussed previously. A BFG exactly matches curved boundaries. the alignment of the grid can result in a more accurate solution than randomly Þlling the space with an enormously large number of smaller and smaller tetrahedrons or hexahedrons. Numerical accuracy will increase further when algebraic turbulence models are employed using an almost orthogonal mesh. In addition. It is therefore felt. numerical diffusion can be reduced. the ßow solver based on the UG approach is substantially slower than for SGs. implicit schemes will be of advantage. that is retaining the viscous terms only in the direction off the body. The use of a structured grid (SG). they cannot be completely unstructured. Since the time-step size in an explicit scheme is governed by the smallest length scale or. have been given in the literature. a SG produces a regular matrix. SGs can be made orthogonal at boundaries and almost orthogonal within the SD. better accuracy is achieved. A fairly complex procedure would be needed to artiÞcially construct these lines. a SD may be covered by a set of hundreds or thousands of blocks. This behavior is not demanded by accuracy. extremely small time steps will be necessary.

In the CP. The origin of this coordinate system in the lower left front–corner (Fig. k in the ξ. The coordinate values itself are given by the proper grid point indices i. the coupling of SGs with UGs is possible as has been shown by Shaw et al.
ζ 6 3
2
4 ξ 1
5 η
Figure 3. or system of reßected shocks or some type of shock-shock interaction. and K. that is the same command Þle is used. Only if a very complex wave patterns evolve due to special physical phenomena. each each cube has its own right-handed coordinate system ´ξ η ζµ. coupled to the ßow solution. J. the UG has advantages. the η direction from front to back and the ζ direction from bottom to top. generating dozens of shock waves. which is followed in Grid []. The coordinate directions in the CP are denoted by ξ.2. On the block level the grid is completely unstructured but in each block the grid is structured and slope continuity is provided across block boundaries.
3.2 Description of the Standard-Cube
The following description is used for both the grid generator and the ßow solver.2: Standard cube in CP (computational plane). such as a bow shock. In addition. for example. [].2.2). That means that the values reach from 1 to I in the ξ direction. which can be perfectly matched by adaptive grid alignment. All computations are done on a standard–cube in the computational plane (CP) as shown in Fig. respectively. respectively. η and ζ directions. η and ζ and block dimensions are given by I. from 1 to J in the η
. Such a grid is called a hybrid grid.26 The large majority of physical phenomena encountered in external and internal ßows exhibit certain well ordered structures. The grid is uniform in the CP. Each cube has its own local coordinate system. A mixture between the boundary Þtted grid approach and the unstructured approach is the use of multiblock grids. j. 3. 3. where the ξ direction goes from left to right.

Thus.2 this corresponds to the bottom and top faces. 3. face 2 will be uniquely deÞned by describing it as a J (η) plane with a j value 1. which are numbered in the following way: 1 bottom. Arrows indicate orientation of faces. A simple notation of planes within a block can be achieved by specifying the normal direction along with the proper coordinate value in that direction. see Fig. face 6 is the ξ η plane with a ζ–value of K 1. respectively. 2 front. First. rules are given how to describe the matching of faces between neighboring blocks. Each grid point represents an integer coordinate value in the computational plane. A standard box has six faces. having ξ–values 0 and I 1. Faces 3 and 4 then denote the left and right right boundaries of the standard box in Fig.3. which means face 2 is deÞned as the ξ ζ plane with a η value of 0 and face 5 is the ξ ζ plane with a η value J 1. face 4 is deÞned by the pair ´I I 1µ.2.e. Each block has its own coordinate system (right handed). plane η 1 to 2 and plane ξ 1 to 3.Numerical Grid Generation
27
direction. Thus faces 4 and 6 are matching with the orientations as shown. followed by the J and K directions. Both are η ζ planes. Face 1 is the ξ η plane with a ζ–value of 0. orientation of block 2 is obtained by rotation of 3 2 about the η–axis – rotations are positive in a counter clockwise sense – and a subsequent rotation of 3 2π about the new ζ–axis. Suppose blocks 1 and 2 are oriented as shown. 1) where the Þrst value is the direction of the normal vector and the second value is the plane index . it is shown how the orientation of a face of a block is determined. arrows are drawn in the direction of increasing coordinate values. i. Second. 3 left. The rule is that the lower coordinate varies Þrst and thereby the orientation is uniquely determined. The orientation of faces between neighboring blocks is determined as follows. The rule is that plane ζ 1 corresponds to 1. In the following the matching of blocks is outlined.5. 5 back. This means the determination of the proper orientation values between the two neighboring faces. determined from the rules
. This notation is also required in the visualization module. Grid points are stored in such away that the I direction is treated Þrst. Face 2 is the front plane and face 5 the back plane. For example. This implies that K planes are stored in sequence. 6 top. and from 1 to K in the ζ direction. 4 right. In Fig.
Y Z
6 4 2 X
physical domai
computational domain
Figure 3. To determine the orientation of a face. which will be numbered in the same way as a dye. The input description for the Grid control (topology) Þle as used in Grid is shown.3. respectively.3: Mapping of a block from SD to CP. For example. by the pair (J.

values 1 to 8 have to be provided. Grid visualization is based on X-Windows and Motif. Grid is a collection of ANSI Croutines and is modeled following the Unix toolbox concept. the package can be implemented on nearly any type of computer. Grid is a general grid generation package for multiblock 2D. The grid visualization module Xivis is based on X11. e.g. Modules are build from reusable subroutines and can be used in arbitrary order.000
.30 shown in the previous Þgure. including Unix PCs. and 3D solution domains that may be of arbitrary shape. there are very efÞcient means for grid point clustering. while a block of about 8.2. Its name is derived from using the star as a wildcard. 3. because the code cannot anticipate in which orientation the distribution is going to be used. and Grid use a combination of algebraic and elliptic grid generation methods.3 The Grid Grid Generation Toolbox
In the following the two codes developed by the authors Grid and GridProwill be presented. consists of about 1500 lines in C.000 points Euler grid for the Hermes space plane: 20 iterations with SOR were used for the complete spacecraft. and provides special features. To speed up convergence. which stands for the collection of the modules that make up Grid . Rotations are denoted by integer 0. to obtain a Navier-Stokes grid from an Euler grid.
3. Grid . if control faces are speciÞed. For the Hermes grid generation an Euler grid of 150. Thus cases 1 and case 7 in Fig. such as grid generation on surfaces.000 points in a few minutes on a PC. Grid point distribution on the block faces is frozen. and 3 marking the corresponding orientations. or Unix station to an SGI Power Challenge etc. that will be used to calculate the control functions for the current face of the block under consideration. Grid can be used in every area where structured grids are needed. surface. All major workstations are supported. e. fairly small in size.000 points was converted into a Navier-Stokes grid of 300. The grid generation modules Grid . allowing this system to be run on virtually any type of computer. which is shown by the fact that several large have been generated and visualized on a PC. Modules for surface grid generation. namely how to match face 4 of block 1 to face 6 of block 2. However. from a PC under Linux. including grid enrichment. 3D grid generation as well as for adaptation are provided. 3D blocks can be interactively cut out of the SD. The same is valid for elements in group 2. we have determined the input needed for Fig. Other codes can be found in the references. the point distribution on the face of a certain block.5. ranging from a PC to a Cray. where elliptic techniques are mainly used for smoothing purposes. 3. each designed for a special task. i. Grid generation is very fast and efÞcient..7 are obtained by rotating case 1 by π 2 and then do the mirroring. This technique was successfully used when Grid generated the 150. that is. 1.g. All cases in group 1 can be obtained by rotating about an angle of 0. 1 2π or 3 2π. Grid is a collection of C routines. all grids presented have been generated by these tools. 2. Therefore. An algebraic-elliptic generation technique is used with convergence acceleration. In addition. Thus. the 3D multiblock grid generation module. The package comprises interactive tools for surface grid generation and for grid enrichment. The control information for a 3D case therefore contains one of the 8 orientation numbers. Solaris.e. This is clear. Grid . while points in the interior are iterated. The code automatically recognizes if the orientation between two faces needs a mirroring.

because it would be completely impractical if each node read the entire input. Each block has its own local coordinate system. For the description of the matching of neighboring blocks in 2D or 3D control structures Òcntrl2d or Òcntrl3d are needed. which are embedded in 2. Òend. In order to identify the object type. In case of a plane. and Grid . which redirects the input to Þle Þlename until an endof-Þle is encountered. Then a Þle reference is made to the respective Þle containing the plane or volume data. object Òplane2d.
.
3. If a multiblock grid with a larger number of blocks (see Sec. the user can reference to a Þle containing that object. the use of a separate Þle for the boundary data of each block is the only possible way. Grid . Furthermore. three coordinate values are required. All grids are slope continuous across block boundaries. which in general is the control (command) Þle. therefore the object name Òvol3d. A volume is part of 3-dimensional space. Òplane3d. All graphics functions are running under X-Windows so the code is fully portable. ??) is to be generated. there are three internal object types. Òvol3d.Numerical Grid Generation
31
points was cut out around the winglets and iterated 100 times to improve grid quality and then was inserted back.or 3-dimensional space. Òplane2d. especially for the large input Þles needed in Grid . which allow a much more compact programming.4
Input for Grid Generation in 2D and 3D
The Grid package expects geometrical objects as input. two coordinate values are needed while for a surface in 3D. object types. describing the connectivity and orientation with respect to neighboring blocks. the actual data must be proceeded by an identiÞer. That is.z). that means. The ÒÞle command can be nested. which are used only inside a module. Òplane3d or Òvol3d as there are blocks. Grid input consists of a control Þle. Òcntrl2d. Special data types are used. For the parallel grid generator this feature has to be used. Hence block boundaries are not visible. The above mentioned ÒÞle command is a feature which makes the input more readable. each point is given by a pair (x.y. Òline3d. The widely used Plot3d (NASA Ames) grid Þle format is supported. Grid description is described by a set of keywords used in the control Þles and the coordinates Þles. The meaning and usage of these objects is as follows. File formats can be either in ASCII or binary. Lines (curves) are needed as boundaries of planes and surfaces. namely Òerror. so the user never has to specify any array dimensions. One additional command is ÒÞle Þlename.and output Þles: Òline2d.2. It is recommended to write the control information in one Þle and to create as many separate objects of type Òplane2d. Òcntrl3d.y) or by a triple (x. object Òplane3d). Instead of placing an object of type Òplane3d or Òvol3d directly after the control information of a block. reading is resumed by the calling Þle. which is automatically detected by the input routines. Volume data are used in a restart case. After that. which in general comprises a large number of coordinate values to describe the geometry of the boundary surface or the volume grid. In Grid dynamic storage allocation is used. Òdigit. The following object types can be speciÞed in Grid input. This concerns the command Þles of modules Grid . Control functions for grid adaptation can also be speciÞd. The same is true for surfaces.

The elliptic grid generation equations as well as the ßow equations are solved in the transformed or Computational Space (Plane) (see Sec . followed by the dimensions of the block (ξ η directions)11 10. an interblock boundary. The Þrst integer gives the number of the edge. In Physical Space (also denoted as the Physical Plane) coordinates are described by (x y z). The command Þle describes the blocking of the solution domain.cmd and the geometry data (edge description) is stored in rect.cmd) has to be speciÞed. In what follows. 3. each block can have its own local coordinate system and these systems can be rotated with respect to each other.4. The next four lines describe the edges of the rectangle. It should be noted that two different coordinate systems are used. The Þrst line of the command Þle contains the keyword cntrl2D that simply denotes a 2D SD. However. edge 2 is north.. the following rule must be obeyed (for three dimensions see later chapter): Rule: With respect to the coordinate axes of the block the numbering of the edges has to be such that: edge 1 is east. The grid point distribution can be generated interactively. sufÞcient for the rectangle example. the second one gives the number of points in the η ´J µ direction (number of rows).e. A side of type 1 is used for interpolating an initial grid. A second input rule determines the orientation of the edges with respect to the axes of the coordinate system. First the so called control Þle or command Þle (extenion . The four edges of a block are numbered 1.). followed by the type of the side. and K..2. i. In Sec. a description of the simple format of the command Þle is given. and edge 4 points south. For each point on this edge the corresponding
.2. It also contains the references to those Þles that contain the surface grid.4.32 3.1 Rectangle Grid Example
In this section we perform the simple task of generating a monoblock grid for a rectangle. to enable the code to perform the automatic detection of rotated blocks. and 4. 3. edge 3 is west. The 1 in next line denotes the block number. In this example.2 we generate a 6 block grid for a diamond shaped body. In this plane coordinate directions are denoted by I. For complex geometries. using the tools (see Chapter Grid Generation Tools) Poli and Spline2D. J. curvilinear coordinates (ξ η ζ) are needed whose directions are conforming to the shape of the conÞguration. In Grid the following rule applies (respectively in three dimensions): Rule: In specifying the dimensions of a block the Þrst number denotes the number of grid points in the ξ ´I µ direction (number of columns). a 2D problem is considered only. the surface grid simply consists of the grid point distribution on the four edges (sides) of the rectangle. namely the block connectivity information and the grid point distribution on the physical boundary of the solution domain (SD). Since this is the Þrst example. The complete command Þle format is described in section ??. a type 2 or 3 indicates a matching boundary. this simple example allows us to introduce most of the concepts needed for advanced grid generation. A type 0 or 1 describes a Þxed boundary. there are only Þxed edges. However. In this example it is assumed that the command (connectivity) information is stored in the Þle rect. In order to generate a grid two types of informa tion have to be provided.lin. 2. hence sides are either of type 0 or 1. In 2D the rotation is determined automatically by the code. Therefore. If a grid comprises several blocks.

making it possible to generate grids with hundreds of thousands of points on a workstation. The last line.10. this line has the form object type for example
. If this edge is Þxed. several edges can be used. Since the interpolation is linear. the sequence of the edge coordinates in Þle Þle rect. the corresponding opposite side in the neighboring block is investigated. the reader should refer back to the example in Fig.4.lin the Þrst line indicates that data are of type line 2d. an interpolated grid may be partly overwritten in this process. The proper type of these boundaries is therefore Òline2d. Þle rect. Having read this information. If an edge is not neighbored. Each data Þle begins with a line indicating the object type.2
Diamond Shape 6 Block Grid Example
To illustrate the multiblock concept. In that case.10.8 is shown in Fig. the algorithm continues through all blocks until a Þxed side is encountered. a simple 6 block grid for a diamond shape is constructed. The next value gives the number of coordinate pairs (x.lin must be consistent with the order of the edges in Þle rect. 3. allowing complete freedom. the Þrst line of this Þle starts with Òcntrl2d. Since the example is 2D.
Ò cntrl2d 1 11 10 110000 200000 300000 400000 Þle rect.lin
In Þle Þle rect. 3. To obtain a good initial grid. The next two integers are needed to describe the block connectivity.lin is a reference to the geometry Þle that contains the coordinate values of the 4 edges of the rectangle in exactly the same order as the edges occur in the command Þle. specifying the block and edge number of the neighboring block.10. In the following a complete description of the topology of multiblock grids is given. the quality of the initial grid depends on the choice of the edge where the interpolation process starts from. 3. Integers 5 and 6 describe the so called control functions (see next Chapter) used for grid point clustering. y) that form this line.8.cmd. a linear interpolation between the two corresponding Þxed edges is performed. see Fig. these numbers are 0. For example. If this edge is of type matching. 3. however. The meaning of the control information is explained below.Numerical Grid Generation
33
point on the opposite edge is known. The control Þle information for the grid of Fig. 3.
3. In Fig. that is are part of a 2D curve. In 3D.2. Moreover. The interpolated grid speeds up the grid generation process by a factor of 5 . It should be noted that all physiacl coordinates are given with respect to a Cartesian coordinate system that holds for all blocks. If the edge is matching.11 the corresponding coordinate values are speciÞed for all Þxed (physical) boundaries. it would have been possible to describe edge 2 Þrst and the edge 1 in the command Þle.

g.
Figure 3.8: A 6 block grid for diamond shaped body. resolving a boundary layer. e. This type of grid line conÞguration cannot be obtained by a mono-block grid.9: Grid line can also be clustered to match the physics of the ßow.Numerical Grid Generation
35
Figure 3.
.

J) plane3d IJ x(1. The object speciÞcation is valid until the next control line is encountered or if the end of the current input Þle is read. command information and geometrical information of objects may be mixed. .. . All control lines. one can start with the control information of type Òcntrl3d to specify the block. which can not be identiÞed are converted to the internal object type error. After that. a line may be given plane3d to tell Grid that surface information for this block follows or a line Þle Þlename could be given to indicate that surface data for the next block are stored in Þle with name Þlename. x(I) y(I) Object Òline2d indicates the boundary line of a 2D SD.. denoting the boundary points. In control Þles.. the input format for the different object types will now be explained in detail.. To be more speciÞc. object speciÞc information is expected. .. x(I..1) y(1.g.1) y(1.1) x(1.e.J) y(I. Objects plane2d and plane3d are of the following form.38
Òcntrl2d
After this control line. plane2d IJ x(1. .1) z(1...2) y(1.1) .J) z(I.. . .J) y(I.. In short notation an object of type Òline2dhas the following form: line2d I x(1) y(1) . x(I.1) y(I.2) . x(I.J)
. The value of I speciÞes the number of data points.. One can write any number of objects (in this example Òcntrl3d) after this control line as needed. .1) . followed by I pairs of x and y values.

1.1) x(1. A similar format is used for the control information in 3D. J are the number of grid points in the respective directions. The next four lines describe the four edges of the block. Hence..1) z(I. 2 north. cntrl3d nos IJK s1 st nb ns nr cb cs cr s2 st nb ns nr cb cs cr s3 st nb ns nr cb cs cr s4 st nb ns nr cb cs cr s5 st nb ns nr cb cs cr s6 st nb ns nr cb cs cr Again.1. The only restriction is that the same order is used when the boundary data are read.1) y(1. The distance of the Þrst layer of grid points of a Þxed boundary (type 6).1) ... s1 to s4 denote the side-number where 1 is east.1.J..1. which must have the same number of grid points as the current side. x(I. A side of type 2 is a matching side (overlap). . In this case. It is not necessary to set s1 to 1 etc. 1 is a Þxed side which is used for computing the initial algebraic grid. so the K planes (I x J coordinates each) are stored consecutively.1) . the corresponding values for nb and ns have to be given where nb is the number of the neighboring block and ns the number of the matching side of this block.J.1) z(1..K) The control information for 2d reads cntrl2d nos IJ s1 st nb ns cb cs s2 st nb ns cb cs s3 st nb ns cb cs s4 st nb ns cb cs where nos is the block number and I. these values should be set to zero. The edge control information can be in any order. If st is 0 or 1. st is the side-type. J and K are the dimensions in x.1) z(I.J.2. these values should reference a Þxed side. and a combination of distance control and orthogonality (type 7) can also be chosen.1. 3 west.2.2.J. 0 means Þxed side. followed by the J and K indices..J. x(I.1) z(1. I.1) y(I.Numerical Grid Generation
39
In Òvol3d the I-index is running fastest. vol3d IJK x(1. The prescribed distance values needed for side type 6 and 7 are speciÞed at the end of the input control Þle. .1.1) y(1. .K) z(I. The 2D code allows for orthogonal grids at a Þxed boundary (type 5).K) y(I.2) z(1. re-
.J. before the Poisson equations are used.1. nos denotes the block number.2) . cb and ct denote the control block and control side from which the distribution of boundary points is to be taken to calculate the controlfunctions for the right hand side of the Poisson equations.y and z-direction.1.1.2) y(1.1) y(I.1) x(1. and 4 south. x(I.

so for each side there is one line with side-speciÞc information. these values are 0. cs. In the 2D case. Variables must be in the following range. where a 1 denotes the side used for initialization. the control values at the opposite boundary are assumed to be zero. The same is true for the control side information cb.40 spectively. rotation value nr is necessary.
¯ s1 s6: [1. In the 3D case a value between 1 and 8 must be given. which is detected by the code. namely the sides are either parallel (+1) or antiparallel (-1).6] ¯ cr: [1. and no control functions are given at the opposite Þxed side. there are only 2 possibilities.
.N] (N is total number of blocks) ¯ ns: [1.6] ¯ nr: [0.N] ¯ cs: [1. The control values within a block are determined by linear interpolation.6] ¯ st: [0. s1 to s6 are the side numbers as described for the standard block.3]
side number
side type neighboring block
¯ nb: [1. If no clustering is desired. Each block has 6 sides. In addition to the neighboring block nb and the neighboring side ns. st again is the side type.8]
neighboring side of block nb rotation needed to orient current side to neighboring side block number for control information side number for control information orientation for control information
For a Þxed side control information can be used from other sides (faces) or from other blocks and the orientation of these sides (faces) has to be speciÞed (cr).3] ¯ cb: [1. and cr. If control functions are speciÞed only on one side. Only integer values can be assumed.

The local grid reÞnement can be achieved without inßuencing the far Þeld grid. 3.14. It is therefore mandatory to localize the grid line distribution. The number of grid lines can be controlled within the block. To this end the clamp clip technique is used. In the remaining solution domain this requirement only causes memory overhead and reduces the convergence speed. If the grid is reÞned along the physical boundary. see Figures 3.
.3 Local Grid Clustering Using Clamp Technique
Along Þxed walls a large number of grid lines is required in order to simulate boundary layers.12: So called clamp technique to localize grid line distribution. the grid lines can be closed in clamp blocks.13: Clamp 1 at hyperboloid ßare.Numerical Grid Generation
41
Figure 3.
3. Using clamp clips.12.13 and 3. The principle of it is to build a closed block system connected to the physical boundary. see Fig. The real power of this technique is demonstrated in the Hyperboloid Flare in F4 Windtunnel grid.
Figure 3. a reÞnement at the outer boundary is also obtained.

.42
Figure 3.14: Clamp 2 at hyperboloid ßare.

3. Fig. generating a huge amount of Þnite volumes in areas that don’t need a Þne resolution. doted lines are block boundaries. 3. As a more complex example of the above strategy a ”hyperboloid ßare in the F4 windtunnel” grid is presented. this reÞnement extends into the far Þeld and thus causes a substantial computational overhead.
3.15 a conventional topology is shown to capture the boundary layers both for the windtunnel walls and the hyperboloid ßare.12 was used. Second. To this end the ”clamp” technique described in Fig. identiÞcation of single objects is not possible.3. Introducing clamp clips to locations near physical boundaries keeps the grid lines local. solid lines are grid lines.15: A conventional topology for a hyperboloid ßare in a windtunnel. The Þrst disadvantage is that all gridlines are extending into the farÞeld. One of the major constraints in the generation of this grid was to contain the high degree of grid line clustering around the hyperboloid ßare close to the body. Numbers denote block numbers. Although the topology allows a reÞnement of the grid from an Euler to a N–S grid.16 shows the local topology of the hyperboloid ßare. 3.Numerical Grid Generation
43
Figure 3.1 Hyperboloid Flare in F4 Windtunnel Grid
In Fig.
.

. This topology is one part of the topology of 284 block grid.44
Figure 3.18.16: Topology of 36 block hyperboloid ßare. 3. see Fig.

Generally. Common to all adaptation techniques is the speciÞcation of a weight function. giving rise to various adaptation techniques. concentrating grid points in areas of larger solution variation.1)
For example the following expressions could be chosen as weight functions: w 1 · β ∂Φ .Grid Adaptation
47
4. In many ßows of interest. from the Ma number distribution. If there are no charges in the SD. 4.2)
w ds
cdξ
(4. the weight function w is a linear combination of the Mi 1 · c1 M1 · c2 M2 ·
w
· cnMn
(4. takes the form
w sξ or in differential form
c
(4.4 can be interpreted in several ways. 4. that is the rhs of the Laplace equation is different from
. When employing elliptic PDEs for grid generation. Adaptive grid methods monitor the solution and adjust the grid dynamically. denoted by M1 M2 Mn . these regions have a regular structure so that SGs are well suited for adaptation. the electrostatic potential is obtained by solving Laplace’s equation.0
Grid Adaptation Techniques
For the adaptation of a grid.2 can be differentiated again.g. Introducing charges into the SD. the weight function can be used to calculate the so called control functions P´ξ ηµ and Q´ξ ηµ (2D) that form the rhs of the elliptic PDEs. ds is arc length. using the gradient or w 1 · β κ where κ denotes the curvature of the Φ. Assuming a constant c. written as Þrst order differential equation. ∂ξ Φ is the so called monitor surface (explained below) and could be constructed.4)
Eqs.2 and 4.3)
where where c constant or c c´ξµ. The Mi are nonnegative functions. and ξ denotes the coordinate of the computational domain. e. leading to the second order differential equation
´w sξµξ
0
(4. The effect of P and Q can be described using an example from electrostatics. The following discussion is restricted to grid point movement. The Þrst adaptation technique presented is concerned with the usage of PDEs. The basic formula for the redistribution of grid points is the so called equidistribution statement for the weight function w that. the Þrst order differential equation of Eq. which are then called Poisson equations. obtained from the solution features of the ßow. knowledge about the ßow solution is needed to obtain the desired grid resolution in all regions that exhibit large variations.

algebraic adaptation is always done on a curve by curve basis. The MS is a piecewise linear (smooth) surface. and a ßow solution is available. obtained by interpreting Eq. the new surface grid on the MS is projected back down to the physical domain to yield an improved grid that more accurately represents the ßow physics and thus reduces the numerical error. it can be the linear combination of one or more of the variables that determine the physics of the problem. the Laplace operator would equidistribute the grid points in the interior. The amount of clustering is given by the weight function. p. The basic features of the adaptation algorithm are given below where it is assumed that a grid already exists in the physical SD. in order to achieve the same clustering in the interior of the SD as on the physical boundaries. It is constructed by the user. 4. ∆η Q (4. namely the equilibrium statement (see below) that can be interpreted as a differential equation. Ma number could be used. w.
The two techniques mentioned above are derived from the same principle.0.2 in a discrete way.5)
Interchanging the role of dependent and independent variables leads to
.48 zero (Poisson equation). e. Adaptation is performed by operating on a surface grid lying on the MS. Without controlfunctions. Regardless of the dimension. leading to a Poisson equation or an integral form. the MS. º after that. if the SD is a 2D region in the ´x yµ plane. is a curved surface in 3D space that is directly above the 2D region. To capture a shock.1 Adaptation by Controlfunctions
Let us reconsider the grid generation equations ∆ξ P. denoted by Φ´x yµ. If the areas where to cluster the grid are known a priori. concentrates the equipotential lines in the vicinity of the charges (negative) or repels the potential lines (positive charge). Examples of monitor surfaces will be given in the next section. The two techniques are described in more detail in the subsequent sections. the boundary point distribution could be speciÞed accordingly and P and Q could be computed from this distribution.
4. for example. based on the arc length of the coordinate curves. to cluster points in the boundary layer.g. The second technique which is widely used is of algebraic nature. which may depend on the gradient or the curvature of the monitor surface. combined with the pressure.
º the initial surface grid is obtained by projecting (lifting) the grid points of the physical domain up to the monitor surface. º surface grid points are subsequently repositioned on the MS by equally distributing the weight function. resulting in an algebraic approach. A monitor surface (MS) (monitor curve in 1D) is positioned above the physical SD..

or an η-line (in 2D). For example. The control functions can be used to adapt the grid lines to better represent certain solution features or to match the distribution of the boundary points in the interior. The latter is performed only once. x and y must be replaced by arc length s.e. 4. where the non-linearity appears in the expression for the metric coefÞcients. ξ varies) and Q from a boundary which is formed from a ξ-line.4 and dividing by w. one can assume that each grid point is coupled to its four nearest neighbors by springs having equal spring constants. the effects of the boundary grid point distribution would be smoothed out within a few grid points and therefore would not be felt in the interior. xξ yηη yη
P:
g11
Q:
g22
(4. however. in aerodynamics. since xξξ . we Þnd from equations g11 xξξ g22 yηη
xξ P. The Þrst process is a dynamic one. yields sξξ
wξ s w ξ
. If the distribution of boundary points is equidistant. all partial derivatives in the other direction vanish.8)
Expanding now Eq. Since a boundary is either a ξ. the opposite behavior is desirable. equations can be solved for P and Q. we obtain the transformed Laplace equation. yη can be calculated from the boundary point distribution. i.e. xξ as well as yηη . To visualize the distribution of the grid lines generated by Laplace’s equation. to choose a dense distribution in the direction off the convex surface to capture the BL. however.Grid Adaptation
g22 xξξ 2g12 xξη · g11 xηη · g´xξ P · xη Qµ g22 yξξ 2g12 yξη · g11 yηη · g´yξ P · yη Qµ
49
0 0
(4. when a storm surge model is used for the calculation of the water level in a bay (convex area). If we determine P from a boundary formed by an η-line (i. If the boundary is not a straight line. For that reason.
yη Q. One of the simplest ways to achieve the desired clustering could be to change the distribution of the boundary points representing the physics. Since the Laplacian produces an equidistribution in the interior. Let us assume that ξ only depends on x and η only depends on y. high resolution near the shore is required or. The equilibrium of such a spring system gives the distribution of the coordinate points. the BL for a convex surface might have to be modeled.6)
These equations are quasilinear. In principle. one observes the following fact: A convex boundary repels grid lines. a Poisson equation is needed and control functions P and Q have to be determined.
(4. there is no difference between these two approaches. before the computation is started.7)
Also. yielding xξξ . If we set the control functions P and Q to zero. Very often. whereas a concave boundary attracts grid lines.

The numerical values of P and Q in the interior are determined by linear interpoloation between two opposite Þxed sides.58953
Figure 4.23004
1.91325
1. i. the weight function w can be directly incorporated in the grid generation equations.e.67494 2.27274 2.2 illustrate the grid adaptation by using control functions. 4. Using the relation for P (same holds for Q). The use of these control functions guarantees that the speciÞed boundary point distribution is obtained in the interior of the SD. arc length s replaces x and y.and η-lines as on the boundary.63224
1. It can be seen that the adapted grids possess a high degree of smoothness [45].Q w
g22
wη w
(4.5 Numerical Parameters: CFL=10 (Implicit Solution) S-W part=10
-4
2.50
Supersonic Inlet: Euler Computation Stationary Solution
200
19x129 Adaptive Grid
150
100
50
0
0
50
100
150
200
250
300
350
400
450
500
550
Mach Number Contour
(computation with restart file)
Physical Parameters: Ma=3 o AoA=-0.
Comparison with the Þrst equation in Eqs.9)
or when resolved for w w c1 exp P dξ g11
where c1 is the constant of integration. For curvilinear boundaries. 4.1 and 4.
. Figs. adapted by control functions.1: 1 block adptive grid for supersonic inlet. this procedure yields the same grid point distribution along ξ.7 leads to
P
g11
wξ .

Examples of monitor surfaces are given in the following three Þgures. [33]. Their extension to higher . The initial grid in the physical SD is uniform. First. Grid points are uniformly distributed on the monitor surface so that arc length spacing is constant. The basic idea of these algorithms is taken from the papers [32].dimensions is straightforward.0. so that they can be used for time dependent problems. Two algorithms will be presented that are extremely fast. When projected down back to the physical solution domain.3: One–dimenisional monitor surface (MS). The following notations are used in one dimension. Lifting up the grid points produces the grid on the monitor surface.4: Repostioned grid on the monitor surface. and [31].
. this results in a clustering according to the gradient of the monitor surface. or 2D grid has already been generated.
4.2
Algebraic Adaptation Algorithms
For the algebraic adaptation it is assumed that a 1D. The algorithm adapts this grid according to the speciÞed set of weight functions that are generally obtained from the physical solution. too.52
ψ
χ
Figure 4. The variable s denotes arc lehgth on the MS
ψ
dsi
χ
Figure 4. consider the one-dimensional problem of grid point clustering.

We Þrst present the redistribution approach as described in [32] and [33]. there is a clustering of grid points where wi is large.5: Repositioning of grid points on the monitor surface according to the magnitude of curvature.g. because there is a one-to-one mapping between the grid points of the physical SD and those on the MS. i
º grid points in the computational domain are denoted by ξi . on the MS. resulting in a clustering of grid points to regions where curvature κ 0. and N is the total number of grid points (or Þnite volumes).
º the weight function depends on arc length. i. one could directly set ΦI : Ma´xI µ -Þltering of Ma values may i i
. the algorithm of Hsu et al. and in the computational SD are given by x s. i i for example. º half point values (cell faces) are averaged.Grid Adaptation
53
ψ
χ
Figure 4. After that.. The one-dimensional version of the algorithms employs the equidistribution statement in the form:
si·1 si
wds
c´ξi·1 ξi µ
where c is determined from the requirement of constant arc length.e. x could also be chosen as the independent variable. sI and the repositioned i i grid points which are given by xi si where i 2´1µN 1. i. [31] is outlined that uses the equidistribution statement with a variable c´ξµ. e. resulting in a grid ´sI µ.e.g. If. Lift up Lift the initial grid ´xI µ from the physical SD up to the MS. wi·1 :
2
1 2
´wi · wi·1µ
Large w result in small ∆si .
I º there are two sets of grid points: the initial grid points denoted by x . Of course. respectively. e. ∆ξi
1 is chosen. 1. and ξ. 2. Ma is used. wi w´si µ. setting c constant. In both algorithms the original arc length of the MS and curve remain constant. MS Determine the monitor surface from the ßow solution. use the Ma number distribution. I º the initial and repostioned points on the monitor surface are also denoted as P and Pi .
º the continuous variables in the physical SD. In general.

Since the left and right sR and boundary values. Inversion The value F ´s j µ is used to search the table from step 4. sN sI N 1 requiring that the ξ j values have uniform grid spacing in the computational domain (as usual) with ∆ξi 1.g. 4.54 be needed . superscript I has been ommitted for notational simplicity) s1 to s on the MS. returning an index m´ jµ for which F ´sI µ F ´s j µ F ´sI ·1 µ j 2´1µN 1 m m Linear interpolation within interval ´m m · 1µ is utilized and yields αj Insertion of α j directly results in sj sI · α j ´sI ·1 sI µ j m m m 2´1µN 1 (4. are known where s1 sI sL . one obtains from Eq. 6.13) F ´s j µ F ´sI ·1 µ m I F ´sm·1 µ F ´sI µ m
(4. namely F ´s j µ j 1 F ´sN µ N 1 j 2´1µN 1 (4. Build table From step 3 a table of values ´sI F ´sI µµ i 1´1µN is obtained. F ´s j µ has to be inverted. Integrate Determine the weight function ( e. using the Euclidean norm. Find F ´s j µ Determine the F values for the repositioning of grid points.11)
It is evident that in order to obtain the new sj values. It should be noted i i that the grid sI is not distributed as required by the equidistribution statement and also would i not lead to a uniform grid in the compuatational domain when projected back and therefore has to be repositioned. denoted as sL and sR . c is obtained by integrating from s1 to sN F ´sN µ c ξN ξ1 where F ´sµ :
s s1
wdl
The values of F ´sµ and F ´sN µ are related by the formula F ´sµ ξ ξ1 F ´sN µ ξN ξ1 (4. 3.11 the implicit condition for the repositioning of the s j values. magnitude of gradient of Ma number) and integrate the equidistribution statement from (all values are taken from the initial grid.10)
Note: In practice the integration is done numerically using the trapezoidal quadrature rule. 4.resulting in a table ´xI ΦI sI µ i 2´1µN 1 where arc length is computed by i i i piecewise linear approximation.12)
. 5.

w 1· κn .Grid Adaptation
7. w´s´ξµµ is needed. The clustering of grid points is along a given curve. Repositioning The new values in the physical domain are found by projecting back down xj xI · α j ´xI ·1 xI µ j m m m 2´1µN 1
55
Using a more compact notation Pj : Pj
´x j y j µ results in PI · α j ´PI ·1 PI µ j 2´1µN 1 m m m
8. λ is chosen to satisfy this constraint. F ´sµ had to be inverted to obtain the new arc length values. η and ζ. Adaptation in 3D is performed by independenty looping over all gridlines in directions ξ. that is s´ξµ has to be inverted to provide ξ ξ´sµ. ci : ∆sI .14)
The constant λ has been introduced to avoid an additional constraint among the c that follows i from the assumption of constant arc length. the algorithm is used as presented. Hence. In that case.e. weight function w has to be determined as a function of ξ and not of arc length s. Elliptic grid generation is used for a few smoothing iterations (3–4). it is used in the weight function. It should be noted that s is now the arc length of the curve along which the points are to be redistrituted. namely into normal and geodesic curvature. Smooth Relax the repositioned grid. i 2´1µN 1 i
. as a consequence. denoted by κ and κg . The algorithm of Hsu does not use a monitor surface. Another possibility is to divide Eq. Since κ represents changes n of the MS itself. but works directly on the curve s s´ξµ. They are determined from the initial grid point distribution. step 7 of the above algorithm has to be replaced by steps 5 to 7 as described in the algorithm of Hsu below. The curvature for this curve can be splitted into two parts. i. n The normal curvature is a measure of the rate at which the direction of the tangent leaves the MS while the geodesic curvature is the curvature of the curve in the MS. The grid was adapted to capture a moving shock and a strong circular gradient (the details are given in [33]). In practical computations. It should be noted that the equilibrium statement in step 3 was integrated such that the resulting integral could be solved for ξ and. Values for ci are not yet known. requiring a MS in 3D space. 4.6 the result of an adaptation is depicted. In Fig. If this approach is to be applied in clustering points along a curve instead along the x axis.
∆si wi
λci
(4. c´ξµ. 4. The algorithm of Hsu uses c denoted by points Pi . that is there is no direct coupling of the points in a plane or within a block. Therefore. The MS algorithm can be extended to a 2D physical SD. To each gridline in the physical SD there is a corresponding curve in the MS.3 by w so that the integration directly gives the values for s. however.

6: Adaptation of a grid using the monitor surface technique to capture a moving shock together with a strong circular gradient (vortex). The initial grid has uniform grid spacing.
.56
Figure 4.

constructed in that way. ¯ The request that total arc length remains Þxed. the weight functions w should be Þltered. is shown. For example. Then the Þne grid in the BL remains unchanged. For adaptation of the oblique shock η lines (ξ = constant) are adapted the using pressure gradient or Ma number gradient. For all of these lines wi 1 is chosen. This can be easily checked by comparing yI ·1 yI with xI ·1 xI .g. For 3D a loop over all coordinate directions is needed. ¯ Adaptation works on a curve by curve basis. while preserving the Þne grid in the BL. the role of x and y is interchanged. ¯ The adaptation is extremely fast. one can adapt the grid to an oblique shock. e. Suppose that grid lines in the BL are ξ lines (η = constant). Arc length has no physical meaning.58
¯ If the straight line in 7 has a slope of more than π 4. On colorplate E a solution adapted grid for a monoblock cone. may be changed. Otherwise kinks in the repositioned grid may occur. The grid is adapted to capture both the oblique shock and the BL. No points from the BL are moved to the oblique shock region. The Þnal grid depicts a very sharp oblique shock and also retains the Þne BL grid.
. i averaged. m m m m ¯ In order to have a smooth grid point distribution. ¯ Adaptation can be performed in a step by step fashion. It can be used for moving meshes.

geometries of high complexity are now of interest as well as very large meshes. which comes close to the ideal situation described above. Consequently. complex turbine geometries. is presented. For example.0 A Grid Generation Meta Language
5. for example. a separation of topology and geometry is not possible. The user does not have to input any surface grids.1 Topology Input Language
With the parallel computers of today substantially more complex ßuid ßow problems can be tackled. computations of up to 30 million grid points have been performed. for instance. surface and volume grids are generated in the same run without any user interaction. the Gridpro code allows the deÞnition of analytic surfaces that are built in or can be described by the user in a C function type style.2. The user has to perform tens of thousands of mouse clicks with no or little reusability of his input. and speciÞes the Þlenames used for the geometry description of the conÞguration to be gridded. it can be reused for a whole class of applications. see below. generating a hierarchy of increasingly complex grid objects. In addition. in turn. This grid generation process can be termed hands off grid generation. The surface can be described as a set of patches (quadrilaterals) or can be given in triangular form. or ßows including combustion have been computed in industry. Complete aircraft conÞgurations. An aircraft. portable. a preprocessor is used that accepts surface deÞnitions following the NASA IGES CFD standard [?] and converts all surfaces into triangular surfaces. Clearly. In addition. TIL allows the construction of complex grids by combining predeÞned objects along with operators for placement of these objects. These surface deÞnitions are the interface to the CAD data. A variety of surface deÞnitions can be used. grid generation codes have to be capable to handle this new class of application. based on the ANSI–C syntax that allows the construction of objects. has a certain topology. In the following a methodology.0. these objects may be composed of
. internally only triangular surfaces are used. since it strongly inßuences the resulting grid line conÞguration. rotated and multiplied. This approaoch has the major advantage that it is fully reusable. Once the topology has been described. 5.Grid Generation Language
59
5. the grid around an engine could be an object (also referred to as component). the basic engine object would have to be duplicated and positioned accordingly. allowing the construction of objects composed of other objects where. and that highly complex grids can be built in a step by step fashion from the bottom up. A compiler type grid generation language has been built. that is. the language should be hierarchical. Convential grid generation techniques derived from CAD systems that interactively work on the CAD data to generate Þrst the surface grid and then the volume grid are not useful for these large and complex grids. The topology deÞnition consumes a certain amount of work. The user provides a (small) input Þle that describes the so called TIL code to built the wireframe model. Moreover. Since an aircraft or spacecraft generally has more than one engine located at different positions of its structure. In general. for instance see Sec. but different geometry data describe different aircraft types. To this end a completely different grid generation approach will be presented.0. These features could be used to build an application speciÞc data base that can be used by the design engineer to quickly generate the grids needed. That is. One step further would be the deÞnition of objects that can be translated.

e. No claims are made that TIL is the only (or the best) implementation of the concepts discussed. It denotes a major deviation from the current interactive blocking approach and offers substantial advantages in both the complexity of the grids that can be generated and the human effort needed to obtain a high quality complex grid. The following Þgures show 2D and 3D grids that have been generated with the automatic zoning (blocking) approach. an aircraft library or a library for automobiles. The surface deÞnition Þles are speciÞed outside of the GridPro environment. The surfaces used are deÞned analytically or are given by a CAD system. The toplogy of a complex grid can be organized into components or objects that may be grouped. In order to generate a grid The following steps have to be followed in the grid generation process.
. A schedule Þle is needed (extension . to obtain a certain convergence level (the user should keep in mind.fra) contain the topology deÞnitions. i.sch) to specify directions for the numerical scheme of the grid generator. but is believed that it is a major step toward a new level of performance in grid generation. however. that there is nothing like a converged grid).1: Navier-Stokes grid for a four-element airfoil. e. a turbomachinery library. The GridPro [35] package reduces the amount of information needed by an order of magnitude. TIL has been used to generate the grids thar are presented below.60 more basic objects etc.
Figure 5. The examples a presented demonstrate the versatility of the approach and show the high quality of the grids generated. comprising 79 blocks. The Topology Input Language [35] has been devised with these features in mind. The versatility and relative ease of use (comparable in difÞculty in mastering LATEX) will be demonstrated by presenting several examples along with their TIL code.g. In this way a library can be built for different technical areas. The TIL Þles (extension . The grid generator needs three different Þles as input.The Þrst layer of grid points off the airfoil contour is spaced on the order of 10 6 based on chord length. in particular when used for parallel computing.

Grid Generation Language

61

Figure 5.2: The Þgure shows the block structure of the four element airfoil generated by GridPro.

Figure 5.3: Grid for a T joint that has numerous industrial applications.

62

Figure 5.4: Sphere in a torus. The input for this grid is presented in the following tables.

Figure 5.5: Complete 3D grid for a generic aircraft with ßaps, constructed from analytical surfaces. However, the topology is exactly the same as for a real aircraft.

Grid Generation Language

63

Figure 5.6: The picture shows a blowup of the engine region of the generic aircraft.Future TSTO or SSTO vehicles will exhibit a similar complex geometry, necessitating both the modelization of internal and external ßows.

Figure 5.7:

64
Figure 5.8:
.

see for example [37] and is not in GridPro. providing a certain smoothness. The labeling only serves to determine the grid topology.Grid Generation Language
65
Step 0 In a preprocessing step. and sOUT.2]). selecting the proper topology may be more effective regarding grid adaptation than the use of additional control functions. cOUT. º The validity of component names (or numbers) is restricted to the component in which they have been deÞned. In other words. Step 1 In the Þrst step. º TIL provides special data types. The major advantage is that TIL can be used to generate simple components Þrst that can the be assembled to form complex grids. a component library can be built up. variables of that data type can be declared. it is stressed.2 . the user designs the block topology. Step 2 In the second step. as a general rule. This is a creative process that cannot be entirely automated.
º Any grid generated by GridPro is comprised by components. In the following some general explanations about TIL are given.
. In the deÞnition of a function its name is followed by the formal list of parameters. e. namely cIN. that the user has to label the surfaces (curves in 2D) that form the boundary of the solution domain (SD). but should not be positioned too far from it. Thus. Any surface is identiÞed by its corners (or vertices). For the proper surface format. The letters c and s stand for corner and surface. This also holds for variable names. sIN. respectively while IN and OUT specify data to be imported or exported. TIL is more speciÞc than a computer language since it also speciÞes the direction of the data transfer.g. º Components in TIL are similar to functions in ANSI C or C++.5.2). all surfaces have to be prepared. next. Again. Table [5. This procedure is called carpeting. the vertex need not be directly on the surface. unless speciÞc rules are built into the grid generator. presented below.g. In the formal parameter list. IN indicates that variables are of input type and OUT denotes output variables. the expertise of the aerospace engineer is needed. The concepts function and component are used interchangeably. Therefore. The blocking depends on the ßow physics and. The choice of the topology itself leads to a clustering of grid points. as can be seen from Tables [5. Approximate corner coordinates should be speciÞed. the user has to specify a list of neighboring (or connecting) corners that may be on other surfaces. A realization of a component is achieved by using a function. which are labeled in increasing order. see below. The component is called using the values of the actual parameter list. º Components itself may comprise more basic components. In the following some of the TIL rules are presented that will be needed to understand the torus-sphere example. that is. º Each function contains a body in which the necessary deÞnitions and statements occur (see. The user has complete freedom in the deÞnition of the surfaces. In addition. enclosed in parentheses (see e. a list of all neighboring vertices has to be speciÞed. Table 5. How the SD is subdivided into surfaces is up to the user.4]. all surfaces (curves in 2D) that form the physical boundaries of the solution domain have to be labeled. The code will do the necessary projections and will also construct a grid on each surface. that is labels for the corners (or vertices) have to be provided.

4.
The input data for the sphere in a torus example. The next INPUT statement reads in component ball. In Table 5.4 that these corner points have to be known in component ball since there is a connection between these cross section corners and the corresponding corners on the ball surface. the 8 corners having been made available by the cOUT statement in the torus argument list serve as input for the ball component. 5. The argument list of ball contains a cIN data type for 8 corners. 1:1 or 1:3. only once.5. each deÞned as a component of type tsec (see below). 4.. º In connecting corners i on component I and j on component J.4) as well as corners 1 4 of component 2 (again a square) of component torus are made available to function main via the cOUT statement. labeled surface 2. The range parameter is given by . the coordinates of the 8 corner points (or vertices) of that cube are connected to corners c 1 .e. In Table 5. containing the coordinates of corners labeled 2.5.5). forming a 4D hypercube topology. presented in Fig. are given in the following tables. The elements of this array can be accessed by. e.2 the topology input for the sphere is depicted.5). numbered surface 1 of type ellipsoid (ellipse in 2D) with 3 numbers denoting the lengths of the major axes in the x.8 used in cIN refers to the numbering in function main(). The following INPUT statement reads in component torus. y-axis goes to the right. which consists of 4 cross sections. The sphere has a radius of 0. and shifts its origin by (0 0 1. Hence.y.
. In the argument list it is speciÞed that the analaytic torus surface is an input surface to component torus. It should be noted that 1:1. and z coordinate directions. if one speciÞes cOUT(2 4. the same corner point on the ball surface is connected to its 3 neighboring corner points on the cube. º Variable names are restricted to the block in which they are deÞned. The central circle of the torus is in the y-z plane (x-axis points up. follwed by the component number and the component name with the associated list of the actual arguments. and 6 in the respective component. Thus. z-axis points away from the reader) with its center at the origin and a radius of 1. The Þrst line after BEGIN speciÞes a surface. 5. The torus surface is labeled surface 1. the second one references the array index. while the cOUT data type uses the numbering of component torus. This distinction is important in order to address the correct variables.h. º In referencing corners there exist several possibilities. the corresponding IN or OUT array is initialized with 0 (for example see Table 5. For example. 0. In addition. 1.6) an array cOUT[1] to cOUT[4] is allocat ed. 5. The Þrst integer gives the component number of the child function.g. It is clear from Fig.25 and its center is at (0. In this context the £ does not denote a multiplication. The array is allocated in the calling function (or parent function).. The INPUT statement is used to read in predeÞned components. The topology of the ball component is equivalent to a cube. it sufÞces to either specify the connectivity in I or in J. i. º Arrays begin with index 1. The bINt() function is the parent function to all other functions (components). and that corners 1 4 of component 1 (these are the 4 inner points of a cross section forming the square that can be seen in Fig. The Þrst line after BEGIN of function main deÞnes the torus surface analytically. 5. This is in contrast to ANSI C where arrays start with index 0. its deÞntion is contained in the Þle surf/torus.3).1 the components comprising the torus-sphere grid are read in.66
º If a negative integer value is speciÞed for the In or OUT data types. The cross section of the torus has a radius of 0.c 8 that belong to two cross sections of the torus component.

which is reßected in its argument list (see Table 5.cIN (2:1.h".
The following 3 numbers specify the initial ccordinates of the corner.8)). -L. The third value is the number of grid points.8). Vertices can be identiÞed by integer values or by variable names. 5.8)).4).3 gives the TIl input for the torus component. The Þrst two variables specify the corners (vertices) of the edge. The Þnal coordinates of the corners are computed by the program. used for the transformation of the original cross section coordinates..(-8). A point on the surface is determined by angles Θ and Φ that are in the y-z plane and in the cross section plane.4. The x command excludes pairs of corner points that otherwise would be falsely connected.cIN (-8).cOUT (1.68
COMPONENT torus() BEGIN s 1 -analyt "surf/torus. The topology of the torus is obtained by cutting it in the cross section plane.707 0 0 0.. This cut results in 2 cross sections (see Fig.. 5.4 for visualization. lists all corner points to which the current corner is connected..707 0)). named tsec (Table 5.3: Topology Þle for the torus component.8)). # Input tsec(torus_section) with transformation # ((1 0 0) (0. x f 1:1 3:1 1:2 3:2 1:3 3:3 1:4 3:4.707 0) (0 0. The torus surface is given by the formula
..8). END Table 5. These coordinate values only serve to provide an initial solution for the grid generation process.. The link option. one of the sphere surfaces (4 corner points). Therefore. The surface of a torus is generated by rotating a circle about an axis that is in the plane of the circle without intersecting it.707 0 0 -0.707 0)* tsec(sIN (1). x f 1:5 3:5 1:6 3:6 1:7 3:7 1:8 3:8. respectively.cIN (3:1.8).707 0)* tsec(sIN (1). INPUT 4 (1 0 0 0 0. The argument list in The command g (see Table 5. while the other 2 are the left and right cross sections connecting to the sphere surfaces. or to another cross section (8 corner points).cOUT (1.(-8). Component tsec represents such a cross section of the torus and thus is a more basic object (again see Fig. Again. reference is made to Fig. The 9 numers after the INPUT statement form a matrix. and to the corner points above and below the sphere of another cross section..707 0 0 0.8)). which is the case for the cross sections at the cut.4). Each cross section contains 8 corners and its corners are either connected to 2 other cross sections ( 8 corner points for each cross section). The -s option followed by an integer indicates the surface that it belongs to.cIN (1:1.4). 5.cOUT (1.(-8).707 0 0 -0. Table 5. # Define surface INPUT 1 (1 0 0 0 0.cOUT (1. It is important to note that the torus component itself comprises 4 components. each tsec component needs an input of 16 corner points.(1:1.. INPUT 2 (1 0 0 0 -0.4) speciÞes the number of grid points to be distributed along an edge.8).707 0)* tsec(sIN (1).707 0)* tsec(sIN (1). INPUT 3 (1 0 0 0 -0.

shown in Figs. ??. A 6 block grid is obtained. Titan is the only moon in the solar system possessing an atmosphere (mainly nitrogen). In order to compute the microaerodynamics caused by the sensors the proper grid has to be generated. GCMS and sensors) is designed. In order to ensure that the laser sensor will function properly dust particles must not be convected to any of the lens surfaces of the optical instruments [46].0. Therefore.9. During the two hour descent measurements of the composition of the atomsphere will be performed by several sensors located at the windward side of tyhe space probe.
.
5. comprising 6 blocks.Grid Generation Language
71
5. However. it is important to note that each of the more complex grids is generated by modifying the TIL code of its predecessor. does not contain the small sensors that are on the windward side of the probe. Thus the Þrst topology is a grid that corresponds to a box in a box.2 TIL Code for the Cassini-Huygens ConÞguration
In the following we present the TIL code to generate a 3D grid for the Cassini–Huygens space probe. A sequence of grids of increasing geometrical complexity has been generated. With increasing comlexity the number of blocks increases as well. Cassini–Huygens is a joint NASA–ESA project and will be launched in 1997. The combination of all objects results in the Þnal grid. comprises 462 blocks. 5. After a ßight time of 8 years. The corner connection to be initialized results in the wire frame topology. shown in Fig.2. To generate this grid the following TIL code has been written.10 and ??.0.1
TIL Code for 6 Block Cassini–Huygens Grid
The general approch for constructing the Cassini–Huygens grids of increasing complexity was to Þrst produce an initial mesh for the plain space probe without any instruments which is of relatively simple shape equivalent to a box. The simplest version. extensive numerical simulations have been performed to investigate this problem. a topology for the Huygens body (without any additional elements such as SSP. This grid has a box–in–box structure: the outer box illustrates the far Þeld and the interior one is the Huygens body. The Þnal grid. the planet Saturn will be reached and the Huygens probe will separate from the Cassini orbiter and enter the Titan atmosphere. The reÞnement of the grid is achieved by adding other elements designed as different objects. As the Þrst step in generating the Huygens grid. Saturn’s largest moon [?]. This topology describes the spherical far Þeld and the body to be simulated. 5. shown in Fig. modeling the SSP and GCMS sensors.

Right: Topology side view. It is bounded by a large spherical far Þeld. the Huygens body and GCMS are clamped by quads. The ratio of the far Þeld radius and the Huygens radius is about 20. The far Þeld. Left: Topology top view.Grid Generation Language
75
Figure 5.11: Illustration of the block topology for the 24 block Huygens grid.
. in which the Huygens space probe is embedded.
Farfield
Farfield
GCMS
GCMS Huygens
Huygens
Figure 5. Again clamps are used for the 3D topological design.10: The 6 block Huygens grid is depicted.

one may generate an intersecting hyperquad with the same topology.
. a grid without folded lines is generated. shown in Fig.12: 24 block grid for Huygens space probe including SSP. the grid is skewed and has folded lines. 5. Based on this initial condition. In the second case.14. For instance. 5. The positions of these points must be given as spatial coordinates. 5.76
Figure 5. shown in Fig.15 an example is given to demonstrate the inßuence of the initial wireframe coordinates on the grid quality.14 and 5.13. there is no skewness in the conÞguration. but different initial coordinates. The accuracy of coordinate information plays an important role for the grid quality. and these positions are used as initial conditions for the grid generator. In the Þrst case. 5.15. shown in Fig.
ordering of points which represent the vertices of the wireframe structure. In Figs.

13: Topological description for a two–loop hyperquad is depicted. the topology remains unchanged while wireframe corners are different.14: The wireframe should reßect the actual geometry. In both cases. Some attention has to be given to the placement of the wireframe points. Although wireframe coordinates are not Þxed.Grid Generation Language
77
A
B b
A
B b
a d D
a d C D
c
c C
Figure 5. a better initial solution will lead to faster convergence and inproved grid quality.
Figure 5.
. It is obvious that in the second case an overlap is encountered and grid lines will be folded.

15: The grid quality is inßuenced by the initial coordinate values of the wireframe corners A distorted wireframe may cause grid skewness and folded lines.
.78
Figure 5.

indicating that it is now the active axis. topology deÞnition means constructing a wireframe model by placing wireframe points. Activating a Surface Surface deÞnition: AZ–manager has some CAD capabililty.
. On the Command Panel under topic Rotation two view options world and eye can be selected. The user can deÞne surfaces in 2D (lines) as well as in 3D (surfaces). The area in the box will be enlarged. RM: Drag the box that is appearing to desired size and then release mouse button.Automatic Blocking
79
6. green the y-axis. AZ–manager that is now described. MB: Rotation (see description below). is such an interactive topology generator. the axis of rotation lies always in the view (screen) plane. The axis becomes a thin line. Under eye view. and activating a surface as well as performing actions on it. the axis of rotation. Moving the mouse in a certain direction in the view plane.e. Surfaces are numbered consecutively. and assigning them to Þxed surfaces.0 Tools for Automatic Topology DeÞnition
The next major step forward in automatic blocking is to provide a tool that by taking graphical input automatically produces the TIL code described in Sec. This feature is particularly useful for external ßows wherethe outer boundary is not of Þxed shape. Under world view. and purple the z-axis. º Constructing the set of wireframe points. Once the topology deÞnition has been done. Yellow denotes the x-axis. Mouse Usage All three mouse buttons are used. To activate an axis. Topology The topology deÞnition process comprises four stages:
º Generating a geometry description or using existing CAD data. linking these points. º Linking corresponding wireframe points. In this context.0. rotating a 3D objects. It should be noted that this axis of rotation is invisible. The topology generation process comprises the following stages. º Assigning wireframe points to surfaces. automatically and uniquely deÞnes the axis of rotation being perpendicular to the direction of motion.1.
Rotation (Orientation of Coordinate System) There are two different modes of rotation. The reader should refer to this section when actively using AZ–manager.
In the following general rules for AZ–manager usage will be outlined. AZ–manager has miniCAD capabilities allowing the user to built his own geometry. move the mouse cursor over the axis such that the rotation symbol occurs and press LB. all grids are produced automatically starting GridPro from within AZ–manager and visualizing the grid while it is being computed. because only then the full meaning of the rules will become clear. concerning mouse usage. a rotation is about the active coordinate axis. i. 5. LB: Drag to translate the object across the screen. It is also described how ro read in CAD data as the general means for surface description of 3D bodies.

translated and rotated to any desired position. The active surface. special care has to be taken by the grid generator to have grid lines close to the box corners. By using or scrolling is done and the corresponding surface is selected.e. The Huygens surface data is in the Þle HuygensSurf. The surface turns blue and then resize by dragging the surface handles. However. is selceted by pressing surf on the command panel (the command panel is the right side of the screen. A second Þle with extension conn is needed by GridPro. The CP can be moved by the user. Stepping through the pictures. the active surface is marked by blue handles. he just has to remember that it is needed by GridPro. a doubly connected region. A body or conÞguartion comprises a set of surfaces. i. Surface Geometry Description In example 3 the surface of the Cassini–Huygens Space Probe is described by a set of 6 rectangular patches. The user needs to have the possibility to scroll through these surfaces and to perform certain actions on them.dat. on which actions are performed. while button move governs the display style during the mouse movement. that have been extracted from CAD data. Only one surface can active at a time. Its usage is described in Example 3.dat. The following set of pictures serve as a step by step introduction on how to actually generate the grid. Since wireframe points are interactively placed.
6. Pressing the hand button. the clustering is performed automatically in GridPro based on boundary curvature. When surfaces are generated. One should experiment with the options available to fully understand their meaning. Cutting Plane (CP) To place wireframe points in 3D space. given in pointwise form. because of a distorted scaling. results will not be identical with the pictures. Fit Objects on Scree When surfaces are read in or constructed they may not Þt on the screen. A clustering algorithm is needed to achieve that. they may not be visible. one selects the current static display style.
. The user does not have to specify any parameters. termed the cutting plane. Setting the Display Style of Objects There are various possibilties of displaying objects on the srceen. one should select lines. using the button stop of the STY LE option on the command panel. All wireframe points generated are automatically positioned on this plane. describing the connectivity of the patches. For the 2D examples below. containing a multitude of commands). that is. Several other IGES formats are also supported. Using a conversion program. Only the active surface is affected by user operations. the Þle HuygensSurf.1 Example Topology 1: A Circle in a Box
The following example constructs a 4 block grid as a 2D projection of a rod (cylinder) in a box. a plane has to be deÞned. the reader should make sure that his results are similar to those depicted in the Þgures. since the outer boundary is formed by a box. i. The format of this Þle is of no concern to the user. Dragging these handles resizes the active surface. Gridlines wrap around the rod in an O-type fashion. they might have to be scaled. Press unzoom as many times as necessary to Þt objects on screen.0.80 starting from 0. Turn hand on and make the surface active.conn is automatically generated.e.

Figure 6.. invoke the automatic blocking manager with command az. the blocking toplogy has to be speciÞed.1: In this Þrst example the grid for a circle in a box has to be generated.
.180 s4
s1
s3
s5 Y 50.
general
topo
surf
grid
3-dim 2-dim .2: To start with the example.50 r=20 s2 0. The AZ–manager window appears. Switch to 2–dim on the menubar. Second.. First the user has to construct the geometry (physical boundaries) as shown in the Þgure. as shown in Þgure.Automatic Blocking
81
160.0 X
Figure 6.

3: AZ–manager has mini CAD capability that is used to construct the boundary of the solution domain..
.4: new text needed AZ–manager has mini CAD capability that is used to construct the boundary of the solution domain..82
general
topo
surf load: file load:-plane . To this end click surf on the menubar and select load:–plane from this menu.
surf id: type: center: normal: orient: n-space: ok
cancel
Figure 6.
grid
Figure 6. To this end click surf on the menubar and select load:–plane from this menu.

.Automatic Blocking
83
surf id: type: center: normal: orient: n-space:
0 -plane 0 0 0 1 0 0 +side
cancel
ok
Figure 6. input contour data as indicated. The Þrst surface to be generatedd. Note that surfaces are consecutively numbered (surf id).6: Input plane information as shown. The next surface is the south side of the box. It is of type textbfplane that actually is a line in 2D.
surf id: type: center: normal: orient: n-space:
1 -plane 0 0 0 0 1 0 +side
cancel
ok
Figure 6. is the west side of the box.5: To construct the outer box.

Note that surface orientation has to be reversed by specifying –side.
. Press unzoom to make a box Þt on screen. This plane is the north side of the box. Note that surface orientation has to be reversed by specifying –side. This surface is the east side of the box.7: Input plane information as shown.84
surf id: type: center: normal: orient: n−space:
2 −plane 160 0 0 1 0 0 −side
cancel
ok
Figure 6.8: Input plane information as shown.
surf id: type: center: normal: orient: n−space:
3 −plane 0 180 0 0 1 0 −side
cancel
ok
Figure 6.

. In order to resize it.e. Position the mouse curser over the handle and drag it in the indicated direction. i. Next. activate the proper surface (blue color) – in this example the proper edge of the box.86
general topo surf SHOW: axis y 1
"hand"
surf
cut
CUT-P: hand 2 x CURRENT: surf 3
Figure 6.10: A surface may not have the desired size. Then activate hand from the CUT–P option on the command panel. deselect cut (if it is on. the radio button is pressed) of the SHOW option on the Command Panel.

Automatic Blocking

87

8

7

4

3

1 5 a

2 b 6

Figure 6.11: Activate the other box surfaces in the same way and drag them as shown in Þgure a. The dragging shoul be done such that the box contour becomes visible. The circle should be visible as well. Now wireframe points can be placed. Wireframe point don’t have to ly on a surface, but they should be close to the surface to which they are going to be assigned. Wireframe points are created by pressing ”c” on the keyboard and clicking at the respective positions with LB. Make sure that Corner Creation Mode appears in the upper left corner of the window. Finish point positioning as depicted in Þgure b.

88

Figure 6.12: Link the corners by pressing ”e” on the keyboard and click each pair of points with LB. Make sure that Link Creation Mode appears in the upper left corner of the window. After completing the linkage process, a wireframe model should like as depicted in the Þgure.

CURRENT surf 0(-1)

Figure 6.13: Activate surface 0. Press ”s” on the keyboard and click the two points, marked by ”X”. If points are selected, their color turns from yellow to white. That is the corners are assigned to surface 0.

as shown in Þgure.18: Since GridPro is 3D internally. A red line between the corresponding diagnoals indicates the succesful performance. The Þrst block to be revoned is formed by the 4 points which build diagnoals. Select Ggrid start to generate the grid. Perform the same action for points marked by crosses.
general
topo TIL read TIL save Ggrid start
Confirm
Current Topology Is Complete.
.Automatic Blocking
91
Figure 6. it is necessary to explicitly remove two blocks. Press ”f” and click two points.19: Save the topology Þle (TIL code) by selecting topo and clicking TIL save. marked by triangles. Launch Ggrid ? ok cancel
Figure 6.

At the right side of the menu window. The generation process is now suspended.21: Select panel=T from the menubar and click Grid.
grid load new
panel=T Topolog Grid Surface
Figure 6.92
general
topo TIL read TIL save Ggrid start Ggrid stop
Figure 6. Click shell by Make Shell at the right low side of the menu window.
. The grid appears. Select grid from the menubar and click load new. select STYLE and click lines under stop.tmp” can be read. The grid appears as wireframe model.20: Select topo from the menubar and click Ggrid stop. the grid Þle ”blk.

0
X
Figure 6.Automatic Blocking
93
6.22: Constrct geometry in the same way as described in example 1.180 s3 110.
Figure 6. would be possible.130 r=20 s0 s5
s2 s4 Y 50.23: It is suggested to place the wirefrace points in a row by row fashion. Other toplogies. of course.50 r=20 s1 0.0.
.2 Example Topology 2: Two Circles in a Box
160.

The outer 8 wireframe points should be placed on a sphere with radius 5000.dat". Select surf and click load:–ellip.
6. the inner 8 points are located on the Huygens surface.Automatic Blocking
SET DIMENSION 2 SET GRIDEN 10 COMPONENT Huygens () BEGIN s 1 -linear "HuygensSurf. at the same time.34: The topology to be generated for Huygens is of type box in a box. The Huygens surface is displaced on the srceen and can be rotated as described in ??. Construct a circle with radius 5000. END
99
Table 6.0.
Y z=-5000
inner box
X
Z z=5000
Figure 6.3 Example Topology 3: Cassini–Huygens Space Probe
In this example we will do our Þrst 3D problem and.2: This TIL code is read in using topo TIL read as in former examples).dat contains the Huygens surface data as described in ??. will use surface description data coming from a CAD system (see above).
. The Þle HuygensSurf.

In addition. being portable to any software or hardware environment. º Gather/Scatter operations on index vectors. º A collective operation. This can be achieved by writing code for an abstract machine.
The abstract machine has the following Þle system methods:
º Collective synchronized Open of a Þle.0 Parallelization Strategy for Complex Geometries
7. º Multiple write operations with a single shared Þle pointer. also implemented as pipeline operations on RISC processors such as the i860. with its particular processors and perhaps a vendor-speciÞc parallel Þle system. º Multiple atomic seek-and-read operations with a shared Þle pointer. if needed by a CFD code. The abstract machine speciÞed above is well suited for unstructured or structured mesh computations. yet taking advantage of the speciÞc machine conÞguration to maximum advantage. then implementing the abstract machine on speciÞc hardware platforms. º Blocking receives of messages chosen by tag or by source processor.
Some parallel processing systems support these operations very efÞciently. The abstract machine has a communication system with these methods:
º Sending locally blocking tagged messages. In writing cost effective parallel software whose quality can be ensured. º A system message buffer of given Þxed size. one has to separate the methods of the application from the details of implementation on a particular parallel architecture.102
7.0. For example the vector operations on a Cray or IBM RS/6000 processor.
. When a processor does not directly support one of the methods. and other such operations. it can be implemented in software.1 The Abstract Parallel Machine
The following philosphy is basic to the parallelizing approach of multiblock codes. called ”reduce” in MPI for sums over all the processors. the abstract machine has a processor with these BLAS1 vector methods:
º Adding vectors and multiplying by a scalar (saxpy/daxpy operations).

3 Encapsulation of Message Passing
Message passing routines are called indirectly.
7. These pieces of code are assembled in a communication library. host-node . In order to minimize code maintenance only a parallel version of the code is kept. making the mapping of several subdomains (blocks) to a single processor unnecessary. because a lenghty optimization stage might be required.0. A single processor machine is simply considered as a one processor parallel machine.not recommended) and I/O modes (host-only I/O. for instance. EfÞcient coding: In parallel programming a number of operations.Parallel Grid Generation
103
7. load balancing is achieved by mapping any arbitrary number of blocks to a single processor. node local I/O.2 General Strategy
Parallelization is achieved by domain decomposition. The decompositions will be performed statically via pre-processing of the grids rather than dynamically during the ßow solutions.0. but leads to a much more ßexible approach. etc. Portability: Encapsulation of message passing routines helps to reduce the effort of porting a parallel application to different message passing environments. using a set of encapsulated routines that in turn call PVM routines. Source code: Encapsulation allows to keep only one source code both for sequential and par-
. Since PVM and MPI lack this type of operations they are realized by the application itself. this may not be the case for very large grids. The differences can also be hidden in the interface library. fast parallel I/O hardware. Since for structured grids the grid topology (that is the number and size of blocks needed to obtain a certain gridline conÞguration) is dictated by the physics. The parallelization of I/O can be very different with respect to the programming models (SPMD. but the communication libraries of Intel and Ncube are directly supported as well as the IBM parallel environment. to use simulated annealing to minimize the cost function. This is adequate for the current problems However.) supported by the parallel machines. All codes are implemented under PVM or MPI. such as global reduction operations (”global sum”) are required quite frequently. This demands a slightly more sophisticated message passing scheme. One might surmise that for unstructured grids a domain decomposition algorithm will result in a load balanced application. However. Do-loop level parallelization will neither be economical nor effective. future developments may include unsteady calculations on moving and/or solution adaptive grids and thus any modiÞcations made will keep these future developments in mind while not allowing them to dominate the initial parallelization effort.

104 allel machines. Maintenance and further development: Encapsulation keeps message passing routines local. Thus, software maintenance and further development will be facilitated. Common message passing subset: The table below lists the basic message passing operations which are available in most of the message passing environments. These routines are functionally equivalent although the semantics of the arguments and return valves may differ in detail. Portability can be highly increased by restricting oneself to use only operations included in the common subset for implementing the interface routines.

º Since each processor of the parallel machine takes one or more blocks, there may not be enough blocks to run the problem on a massively parallel machine. We have tools to automatically split the blocks to allow the utilization of more processors. º The blocks may have very different sizes, so that the blocks must be distributed to the processors to produce a reasonable load balance. We have tools to solve this bin-packing problem by an algorithm discussed below.

An extremely simple message-passing model is ismplemented, consisting of only send and receive. The simplicity of this model implies easy portability. For an elementary Laplace solver on a square grid, each gridpoint takes the average of its four neighbors, requiring 5 ßops, and communicates 1 ßoating-point number for each gridpoint on the boundary. The Grid grid generator is a more sophisticted elliptic solver, where about 75 ßops occur per internal grid point, while grid coordinates are exchanged across boundaries. Ourßow solver,ParNSS, in contrast does a great deal of calculation per grid point, while the amount of communication is still rather smalli. Thus we may expect any implicit ßow solverto be highly efÞcient in terms of communication, as shown by the timing results below. When the complexity of the physics increases, with turbulence models and chemistry, we expect the efÞciency to get even better. This is why a ßow solveris a viable parallel program even when running on a workstation cluster with slow communication (Ethernet).

Parallel Machines and CFD For the kinds of applications that we are considering in this report, we have identiÞed four major issues concerning parallelism, whether on workstation clusters or massively parallel machines. Load balancing As discussed above, the number of blocks in the grid must be equal to or larger than the number of processors. We wish to distribute the blocks to processors to produce an almost equal number of gridpoints per processor; this is equivalent to minimizing the maximum number of gridpoints per processor. We have used the following simple algorithm to achieve this.

Parallel Grid Generation

105

The largest unassigned block is assigned to the processor with the smallest total number of gridpoints already assigned to it, then the next largest block, and so on until all blocks have been assigned. Given the distribution of blocks to processors, there is a maximum achievable parallel efÞciency, since the calculation proceeds at the pace of the slowest processor, i.e. the one with the maximum number of gridpoints to process. This peak efÞciency is the ratio of the average to the maximum of the number of gridpoints per processor, which directly proceeds from the standard deÞnition of parallel efÞciency. Convergence For convergence acceleration a block-implicit solution scheme is used, so that with a monoblock grid, the solution process is completely implicit, and when blocks are small, distant points become decoupled. Increasing the number of processors means that the number of blocks must increase, which in turn may affect the convergence properties of the solver. It should be noted [?] that any physical ßuid has a Þnite information propagation speed, so that a fully implicit scheme may be neither necessary nor desirable. Performance It is important to establish the maximum achievable performance of the code on the current generation of supercomputers. Work is in progress with the Intel Paragon machine and a Cray YMP. Scalability Massively parallel processing is only useful for large problems. For a ßow solver, we wish to determine how many processors may be effectively utilized for a given problem size – since we may not always run problems of 10 000 000 gridpoints.

7.0.4 Aspects of Software Engineering: C versus Fortran Only a few remarks shall serve to highlight the discussion. Fortran has been the programming language of science and engineering since its inception in the late Þfties. It has been updated in 1977 and recently new versions called Fortran90 (or HPF) have been devised. On the other hand, ANSI-C is a quasi de-facto standard that is available on virtually any hardware platform. In the following, the most important features of a modern high-level language are listed, mandatory to write portable, safe, robust, efÞcient, and reusable code that can also be maintained.

º F77: does not possess modern constructs º HP Fortran and F90: some modern features available º C: all features available, except classes ; dynamic storage allocation limited º C++: all features available Recommendation : C++ language of choice for the 90’s, but not all of its advanced features are needed for CFD

7.0.5 Numerical Solution Strategy: A Tangled Web It is important to note that the succesful solution of the parallel ßow equations can only be performed by a Triad numerical solution procedure. Numerical Triad is the concept of using grid generation, domain decomposition, and the numerical solution scheme itself. Each of the three Triad elements has its own unique contribution in the numerical solution process. However, in the past, these topics were considered mainly separately and their close interrelationship has not been fully recognized. In fact, it is not clear which of the three topics will have the major contribution to the accurate and efÞcient solution of the ßow equations. While it is generally accepted that grid quality has an inßuence on the overall accuracy of the solution, the solution dynamic adaptation process leads to an intimate coupling of numerical scheme and adaptation process, i.e. the solution scheme is modiÞed by this coupling as well as the grid generation process. When domain decomposition is used, it may produce a large number of independent blocks (or subdomains). Within each subdomain a block-implicit solution technique is used, leading to a decoupling of grid points. Each domain can be considered to be completely independent of its neighboring domains, parallelism simply being achieved by introducing a new boundary condition, denoted as inter-block or inter-domain boundary condition. Updating these boundary points is done by message–passing. It should be noted that exactly the same approach is used when the code is run in serial mode, except that no messages have to be sent to other processors. Instead, the updating is performed by simply linking the receive

1). or a point with it. a full coupling of all points in the solution domain would be unphysical. A major question arises in how the decomposition process affects the convergence rate of the implicit scheme. The update process via boundaries therefore should be sufÞcient. the N-S equations can be considered hyperbolic. Hence. This guarantees that the numerical solution is independent of the block topology. because of the Þnite propagation speed. needed to compute the mixed derivatives the viscous terms. Thus. 26 messages would have to be sent (instead of 6 for updating the faces) to fully update the boundaries of this block.
buffer of a block to its corresponding neighboring send buffer. Instead. and is therefore not desirable and not needed. To continue the discussion of convergence speed it should be remembered that for steady state computations implicit techniques converge faster than fully explicit schemes. Instead.Parallel Grid Generation
107
h
Figure 7. However. 7. but may have larger error coefÞcients. To retain the numerical order across block (domain) boundaries. A total of 26 messages would be needed to update values along diagonals. In all other cases. This would lead to an unacceptable large number of messages. the missing information is constructed by Þnite difference formulas that have the same order of truncation error. Hence. parallelizing a multi-block code neither demands rewriting the code nor changing its structure. an edge. in this case the Newton method will converge quadratically. This only occurs in the boundary layer when a steady state has been reached or has almost been reached. only block faces are updated (maximal 6 messages) and values along diagonals are approximated by a Þnite difference stencil. it should be noted that the N-S equations are not elliptic. an overlap of two points in each coordinate direction has been implemented. The small cube in the middle is surrounded by 26 blocks that share a face. It would be uneconomical to send these diagonal values by message passing. First. The only restriction comes from the computation of ßow variables along the diagonals on a face of a block (see Fig.1: Flow variables are needed along the diagonals to compute mixed second derivatives for viscous terms. Imagine a set of 27 cubes with edge length h 3 assembled into a large cube of edge length h.
. unless the time derivative is omitted and inertia terms are neglected (Stokes equations). since the initial solution is close to the Þnal solution.

Points marked by a cross are used for inviscid ßux computation.
. Diagonal points (circles) are needed to compute the mixed derivatives in the viscous ßuxes.108
3 3 2
4 4 1
1
Figure 7. Care has to be taken when a face vanishes and 3 lines intersect.2: The Þgure shows the computational stencil.

Because of their mutual interaction these approaches must not be separated. decomposing the solution domain should result in a convergence speed up. combined with a Conjugate-Gradient technique for convergence acceleration within a Newton-Iteration. the scheme is fully explicit and hence computationally less efÞcient than the fully implicit scheme. On the other hand. the concept of numerical solution procedure is much more general than devising a single numerical scheme for discretizing the N-S equations. although boundary values are dynamically updated. In the preconditioning process used for the ConjugateGradient technique. it should be evident that only a combination of grid generation scheme. since the inversion of a set of small matrices is faster than the inversion of the single large matrix. an optimal decomposition topology must exist that most likely depends on the ßow physics and the type of implict solution process. derived from the discretized N-S equations. the eigenvalue spectrum may be compressed. Having smaller matrices the condition number should not increase. clearly demonstrating the convergence speed up. Thus. domain decomposition may have a direct inßuence on the convergence speed of the numerical scheme. because the resulting matrices are smaller. using physical reasoning it is concluded that in general the condition number should decrease. the basis of the numerical solution technique is the Newton method. In this paper. However. since a full coupling is not required by the physics. if the decomposition leads to a blocksize of 1 point per block.
. a number of numerical experiments has been performed with the ParNSS code [?]. In other words. Block numbers have been varied from 5 to 140. Only the implementation of this interconnectedness in a parallel solver will lead to the optimal design tool.Parallel Grid Generation
109
The former are generally more computationally efÞcient. It is shown in [?] that this ratio is a measure of the convergence speed for generalized conjugate residual algorithms. using an otherwise identical grid. no theory has been developed. domain decomposition may be used to decrease the condition number (ratio of largest to smallest eigenvalues) of the matrix forming the left hand side. From these remarks. However. So far. in particular for meshes with large variations in grid spacing. Therefore. and domain decomposition approach will result in an effective. numerical solution procedure. Second. general numerical solution strategy for the parallel N-S equations on complex geometries.