text:
Main Handbook of grid generation
Handbook of grid generationJoe F Thompson, B K Soni, N P Weatherill (eds.)
Block-Structured Grids -- Unstructured Grids -- Surface Definition -- Adaptation and Quality
You can write a book review and share your experiences. Other readers will always be interested in your opinion of the books you've read. Whether you've loved the book or not, if you give your honest and detailed thoughts then people will find new books that are right for them.
©1999 CRC Press LLC Acquiring Editor: Project Editor: Marketing Manager: Cover design: Manufacturing Manager: B. Stern Sylvia Wood J. Stark Dawn Boyd Carol Slatter Library of Congress Cataloging-in-Publication Data Thompson, Joe F. Handbook of Grid Generation / Joe F. Thompson, Bharat Soni, Nigel Weatherill, editors. p. cm. Includes bibliographical references and index. ISBN 0-8493-2687-7 (alk. paper) 1. Numerical grid generation (Numerical analysis) I. Thompson, Joe F. II. Soni, B.K. III. Weatherill, N.P. QA377.H3183 1998 519.4--dc21 98-34260 CIP This book contains information obtained from authentic and highly regarded sources. Reprinted material is quoted with permission, and sources are indicated. A wide variety of references are listed. Reasonable efforts have been made to publish reliable data and information, but the author and the publisher cannot assume responsibility for the validity of all materials or for the consequences of their use. Neither this book nor any part may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, microfilming, and recording, or by any information storage or retrieval system, without prior permission in writing from the publisher. All rights reserved. Authorization to photocopy items for internal or personal use, or the personal or internal use of specific clients, may be granted by CRC Press LLC, provided that $.50 per page photocopied is paid directly to Copyright Clearance Center, 27 Congress Street, Salem, MA 01970 USA. The fee code for users of the Transactional Reporting Service is ISBN 0-8493-2687-7/99/$0.00+$.50. The fee is subject to change without notice. For organizations that have been granted a photocopy license by the CCC, a separate system of payment has been arranged. The consent of CRC Press LLC does not extend to copying for general distribution, for promotion, for creating new works, or for resale. Specific permission must be obtained in writing from CRC Press LLC for such copying. Direct all inquiries to CRC Press LLC, 2000 Corporate Blvd., N.W., Boca Raton, Florida 33431. Trademark Notice: Product or corporate names may be trademarks or registered trademarks, and are used only for identification and explanation, without intent to infringe. © 1999 by CRC Press LLC No claim to original U.S. Government works International Standard Book Number 0-8493-2687-7 Library of Congress Card Number 98-34260 Printed in the United States of America 1 2 3 4 5 6 7 8 9 0 Printed on acid-free paper ©1999 CRC Press LLC Foreword Grid (mesh) generation is, of course, only a means to an end: a necessary tool in the computational simulation of physical field phenomena and processes. (The terms grid and mesh are used interchangeably, with identical meaning, throughout this handbook.) And grid generation is, unfortunately from a technology standpoint, still something of an art, as well as a science. Mathematics provides the essential foundation for moving the grid generation process from a user-intensive craft to an automated system. But there is both art and science in the design of the mathematics for — not of — grid generation systems, since there are no inherent laws (equations) of grid generation to be discovered. The grid generation process is not unique; rather it must be designed. There are, however, criteria of optimization that can serve to guide this design. The grid generation process has matured now to the point where the use of developed codes — freeware and commercial — is generally to be recommended over the construction of grid generation codes by end users doing computational field simulation. Some understanding of the process of grid generation — and its underlying principles, mathematics, and technology — is important, however, for informed and effective use of these developed systems. And there are always extensions and enhancements to be made to meet new occasions, especially in coupling the grid with the solution process thereon. This handbook is designed to provide essential grid generation technology for practice, with sufficient detail and development for general understanding by the informed practitioner. Complete details for the grid generation specialist are left to the sources cited. A basic introduction to the fundamental concepts and approaches is provided by Chapter l, which covers the state of practice in the entire field in a very broad sweep. An even more basic introduction for those with little familiarity with the subject is given by the Preface that precedes this first chapter. Appendixes provide information on a number of available grid generation codes, both commercial and freeware, and give some representative and illustrative grid configurations. The grid generation process in general proceeds from first defining the boundary geometry as discussed in Part III. Points are distributed on the curves that form the edges of boundary sections. A surface grid is then generated on the boundary surface, and finally, a volume grid is generated in the field. Chapter 13, although directed at structured grids, gives a general overview of the entire grid generation process and the fundamental choices and considerations involved from the standpoint of the user. Chapter 2, though also largely directed at structured grids, covers essential mathematical elements from tensor analysis and differential geometry relevant to the entire subject, particularly the aspects of curve and surfaces. The other chapters of this handbook cover the various aspects of grid generation in some detail, but still from the standpoint of practice, with citations of relevant sources for more complete discussion of the underlying technology. The chapters are grouped into four parts: structured grids, unstructured grids, surface definition, and adaptation/quality. An introduction to each part provides a road map through the material of the chapters. A source of fundamentals on structured grid generation is the 1985 textbook of Thompson, Warsi, and Mastin, now out of print but accessible on the Web at www.erc.msstate.edu. A recent comprehensive text of both structured and unstructured grids is that of Carey 1997 from Taylor and Francis publishers. ©1999 CRC Press LLC The first step in generating a grid is, of course, to acquire and input the boundary data. This boundary data may be in the form of output from a CAD system, or may simply be sets of boundary points acquired from drawings. CAD boundary data are generally in the form of some parametric description of boundary curves and surfaces, typically consisting of multiple segments for which assembly and some adjustments may be required. Point boundary data may be in the form of 1D arrays of points describing boundary curves and 2D arrays for boundary surfaces, or could be an unorganized cloud of points on a surface. In the latter case, conversion to some surface tessellation or parametric description is required. These initial steps of boundary definition are common in general to both structured and unstructured grid generation. And, unfortunately, considerable human intervention may be necessary in this setup phase of the process. The setup of the boundary definition from the CAD approach is discussed in general in Chapter 13, while details of application, together with procedures for boundary curve and surface parametric representations, are covered in Part III. There is then the fundamental choice of whether to use a structured or unstructured grid. Structured grids are covered in Part I, and unstructured grids are covered in Part II. The next step with either type of grid is the generation of the corresponding type of grid on the boundary surfaces — preceded, of course, by a distribution on points on the curves that form the edges of these surfaces. This surface grid generation is covered in Chapters 9 and 19 for structured and unstructured grids, respectively. Finally, the quality of the grid, with relation to the accuracy of the numerical solution being done on the grid, and the adaptation of the grid to improve that accuracy are covered in Part IV. Grid generation is still under active research and development, particularly in regard to automation, adaptation, and hybrid combinations. This handbook is therefore necessarily a snapshot in time, especially in these areas, but much of the material has matured now, and this collection should be of enduring value as a source and reference. Bharat K. Soni Joe F. Thompson Nigel P. Weatherill Starkville, MS, and Swansea, Wales, UK ©1999 CRC Press LLC Contributors Michael J. Aftosmis Gerald Farin Olivier-Pierre Jacquotte NASA Ames Research Center Moffett Field, CA Arizona State University Tempe, AZ Research Directorate (DRET) Paris, France Timothy J. Baker David R. Ferguson Brian A. Jean Princeton University Princeton, NJ The Boeing Company Seattle, WA U.S. Army Corps of Engineers Waterways Experiment Station Vicksburg, MS Mark W. Beall Luca Formaggia Rensselaer Polytechnic Institute Troy, NY Ecole Polytechnique Federale de Lausanne Lausanne, Switzerland Yannis Kallinderis Timothy Gatzke O.B. Khairullina The Boeing Company St. Louis, MO Urals Branch of the Russian Academy of Sciences Ekaterinburg, Russia Marsha J. Berger Courant Institute New York University William M. Chan MCAT, Inc. at NASA Ames Research Center Moffett Field Paul-Louis George Zheming Cheng Bernd Hammann Program Development Corporation White Plains, NY University of California at Davis Davis, CA Hugues L. de Cougny Rensselaer Polytechnic Institute Troy, NY Luís Eça Technical University of Lisbon Lisbon, Portugal Peter R. Eiseman Program Development Corporation White Plains, NY Austin L. Evans NASA Lewis Research Center Cleveland, OH ©1999 CRC Press LLC INRIA Le Chesnay Cedex, France O. Hassan University of Wales Swansea Swansea, UK Jochem Häuser CLE Salzgitter Bad Salzgitter, Germany Frédéric Hecht INRIA Le Chesnay Cedex, France Sergey A. Ivanenko Computer Center of the Russian Academy of Sciences Moscow, Russia University of Texas Austin, TX Ahmed Khamayseh Los Alamos National Laboratory Los Alamos, NM Andrew Kuprat Los Alamos National Laboratory Los Alamos, NM Kelly R. Laflin North Carolina State University Raleigh, NC Kunwoo Lee Seoul National University Seoul, Korea David L. Marcum Mississippi State University Starkville, MS C. Wayne Mastin Nichols Research Corporation Vicksburg, MS D. Scott McRae E. J. Probert Joe F. Thompson North Carolina State University Raleigh, NC University of Wales Swansea Swansea, UK Mississippi State University Starkville, MS Robert L. Meakin Anshuman Razdan O.V. Ushakova Army Aeroflightdynamics Directorate (AMCOM) Moffett Field, CA Arizona State University Tempe, AZ Urals Branch of the Russian Academy of Sciences Ekaterinburg, Russia John E. Melton Robert Schneiders MAGMA Giessereitechnologie GmbH Aachen, Germany Zahir U.A. Warsi NASA Ames Research Center Moffett Field, CA David P. Miller Jonathon A. Shaw Nigel P. Weatherill NASA Lewis Research Center Cleveland, OH Aircraft Research Association Bedford, U.K. University of Wales Swansea Swansea, UK K. Morgan A.F. Sidorov Yang Xia University of Wales Swansea Swansea, UK Urals Branch of the Russian Academy of Sciences Ekaterinburg, Russia CLE Salzgitter Bad Germany Robert M. O'Bara Rensselaer Polytechnic Institute Troy, NY Sangkun Park Mark S. Shephard Rensselaer Polytechnic Institute Troy, NY Information Technology R&D Center Seoul, Korea Robert E. Smith J. Peraire Bharat K. Soni Massachusetts Institute of Technology Cambridge, MA Mississippi State University Starkville, MS J. Peiró Imperial College London, UK ©1999 CRC Press LLC NASA Langley Research Center Hampton, VA Stefan P. Spekreijse National Aerospace Laboratory (NLR) Emmeloord, The Netherlands Mississippi State University Starkville, MS Tzu-Yi Yu Chaoyang University of Technology Wufeng, Taiwan Paul A. Zegeling University of Utrecht Utrecht, The Netherlands Acknowledgments Grid (mesh) generation is truly a worldwide active research area of computation science, and this handbook is the work of individual authors from around the world. It has been a distinct pleasure, and an opportunity for professional enhancement, to work with these dedicated researchers in the course of the preparation of this book over the past two years. The material comes from universities, industry, and government laboratories in 10 countries in North America, Europe, and Asia. And we three are from three different countries of origin, though we have collaborated for years. The attention to quality that has been the norm in the authoring of these many chapters has made our editing responsibility a straightforward process. These chapters should serve well to present the current state of the art in grid generation to practitioners, researchers, and students. The assembly and editing of the material for this handbook from all over the world via the Internet has been a rewarding experience in its own right, and speaks well for the potential for worldwide collaborative efforts in research. Our thanks go to Mississippi State University and the University of Wales Swansea for the encouragement and support of our efforts to produce this handbook. Specifically at Mississippi State, the work of Roger Smith in administering the electronic communication is to be noted, as are the efforts of Alisha Davis, who handled the word processing. Bob Stern of CRC Press has been great to work with and appreciation is due to him for recognizing the need for this handbook and for his editorial guidance and assistance throughout its preparation. His efforts, and those of Sylvia Wood, Suzanne Lassandro and Dawn Mesa, also at CRC, have made this a pleasant process. We naturally are especially grateful for the support of our wives, Purnima, Emilie, Barbara, and our families in this and all our efforts. And finally, Mississippi and Wales — two great places to live and work. Bharat K. Soni Joe F. Thompson Nigel P. Weatherill Author/Editors ©1999 CRC Press LLC Preface: An Elementary Introduction Joe F. Thompson, Bharat K. Soni, and Nigel P. Weatherill This first section is an elementary introduction provided for those with little familiarity with grid (mesh) generation in order to establish a base from which the technical development of the chapters in this handbook can proceed. (The terms grid and mesh are used interchangeably throughout with identical meaning.) The intent is not to introduce numerical solution procedures, but rather to introduce the idea of using numerically generated grid (mesh) systems as a foundation of such solutions. P-1 Discretizations The numerical solution of partial differential equations (PDEs) requires first the discretization of the equations to replace the continuous differential equations with a system of simultaneous algebraic difference equations. There are several bifurcations on the way to the construction of the solution process, the first of which concerns whether to represent the differential equations at discrete points or over discrete cells. The discretization is accomplished by covering the solution field with discrete points that can, of course, be connected in various manners to form a network of discrete cells. The choice lies in whether to represent the differential equations at the discrete points or over the discrete cells. P-1.1 Point Discretization In the former case (called finite difference), the derivatives in the PDEs are represented at the points by algebraic difference expressions obtained by performing Taylor series expansions of the solution variables at several neighbors of the point of evaluation. This amounts to taking the solution to be represented by polynomials between the points. This can be unrealistic if the solution varies too strongly between the points. One remedy is, of course, to use more points so that the spacing between points is reduced. This, however, can be expensive, since there will then be more points at which the equations must be evaluated. This is exacerbated if the points are equally spaced and strong variations in the solution occur over scattered regions of the field, since numerous points will be wasted in regions of small variation. An alternative, of course, is to make the points unequally spaced. P-1.2 Cell Discretization The other possibility of this first bifurcation is to return the PDEs to their more fundamental integral form and then to represent the integrals over discrete cells. Here there is yet another bifurcation — whether to represent the solution variables over the cell in terms of selected functions and then to integrate ©1999 CRC Press LLC FIGURE 1 these functions analytically over the volume (finite element), or to balance the fluxes through the cell sides (finite volume). The finite element approach itself comes in two basic forms: the variational, where the PDEs are replaced by a more fundamental integral variational principle (from which they arise through the calculus of variations), or the weighted residual (Galerkin) approach, in which the PDEs are multiplied by certain functions and then integrated over the cell. In the finite volume approach, the solution variables are considered to be constant within a cell, and the fluxes through the cell sides (which separate discontinuous solution values) are best calculated with a procedure that represents the dissolution of such a discontinuity during the time step (Riemann solver). P-2 Curvilinear (Structured) Grids The finite difference approach, using the discrete points, is associated historically with rectangular Cartesian grids, since such a regular lattice structure provides easy identification of neighboring points to be used in the representation of derivatives, while the finite element approach has always been, by the nature of its construction on discrete cells of general shape, considered well suited for irregular regions, since a network of such cells can be made to fill any arbitrarily shaped region and each cell is an entity unto itself, the representation being on a cell, not across cells. P-2.1 Boundary-Fitted Grids The finite difference method is not, however, limited to rectangular grids and has long been applied on other readily available analytical coordinate systems (cylindrical, spherical, elliptical, etc.) that still form a regular lattice. albeit curvilinear, that allows easy identification of neighboring points. These special curvilinear coordinate systems are all orthogonal, as are the rectangular Cartesian systems, and they also can exactly cover special regions (e.g., cylindrical coordinates covering the annular region between two concentric circles) in the same way that a Cartesian grid fills a rectangular region. The cardinal feature in each case is that some coordinate line is coincident with each portion of the boundary. In fact, these curvilinear systems can be considered to be logically rectangular, and from a programming standpoint are no different, conceptually, from the Cartesian system. Thus, for example, the cylindrical grid in Figure 1, where the radial coordinate r varies from r1 on the inner boundary to r2 on the outer and the azimuthal coordinate θ varies from 0 to 2π, can be diagrammed logically as shown in Figure 2. ©1999 CRC Press LLC FIGURE 2 The continuity of the azimuthal coordinate can be represented by defining extra "phantom" columns to the left of 0 and to the right of 2π and setting values on each phantom column equal to those on the corresponding "real" columns inside of 2π and 0, respectively. This latter, logically rectangular, view of the cylindrical grid is the one used in programming anyway, and without being told of the cylindrical configuration, a programmer would not realize any difference here from programming in Cartesian coordinates — there would simply be a different set of equations to be programmed on what is still a logically rectangular grid, e.g., the Laplacian on a Cartesian grid (with ξ = x and η = y), 2 ∇ f = f + fηη ξξ becomes (with ξ = θ and η = r) ∇2 f = f fη ξξ + fηη + η η2 on a cylindrical grid. The key point here is that in the logical (i.e., programming) sense there is really no essential difference between Cartesian grids and the cylindrical systems: both can be programmed as nested loops; the equations simply are different. Another key point is that the cylindrical grid fits the boundary of a cylindrical region just as the Cartesian grid fits the boundary of a rectangular region. This allows boundary conditions to be represented in the same manner in each case also (see Figure 3). By contrast, the use of a Cartesian grid on a cylindrical region requires a stair-stepped boundary and necessitates interpolation in the application of boundary conditions (Figure 4) — the proverbial square peg in a round hole. P-2.2 Block Structure (The Sponge Analogy) The best way to visualize the correspondence of a curvilinear grid in the physical field with a logically rectangular grid in the computational field is through the sponge analogy. Consider a rectangular sponge within which an equally spaced Cartesian grid has been drawn. Now wrap the sponge around a circular cylinder and connect the two ends of the sponge together. Clearly the original Cartesian grid in the sponge now has become a curvilinear grid fitted to the cylinder. But the rectangular logical form of the grid lattice is still preserved, and a programmer could still operate in the logically underformed sponge in constructing the loop and the difference expressions, simply having been given different equations to program. The correspondence of "phantom" points just outside one of the connected faces of the sponge ©1999 CRC Press LLC FIGURE 3 FIGURE 4 with "real" points just inside the face to which it is connected is clear — this is simply the correspondence of 370° with 10° in a cylindrical system. Such a sponge could just as well be around a cylinder of noncircular cross section, regardless of the cross-sectional shape. To carry the analogy further, the sponge could, in principle, be wrapped around a body of any shape, or could be expanded and compressed to fill any region (e.g., expanding to fill a sphere), again producing a curvilinear grid filling the region and having the same correspondence to a logically rectangular grid (Figure 5). The programmer need not know, in fact, what has been done to the sponge. It is also clear from this analogy that the sponge could deform in time; i.e., the curvilinear grid could move in physical space, while the logically rectangular grid could still be considered fixed in computational space (image the sponge filling a beating heart). Again, the programmer need not be told that the boundaries are moving, but simply again be given a different set of equations that will include a transformation of the time derivatives as well. It is not hard to see, however, that for some boundary shapes the sponge may have to be so greatly deformed that the curvilinear grid will be so highly skewed and twisted that it is not usable in a numerical solution. The solution to this problem is to use not one, but rather a group of sponges to fill the physical field. Each sponge has its own logically rectangular grid that deforms to a curvilinear grid when the sponge is put in place in the field. Each sponge now abuts with adjacent sponges, and the correspondence ©1999 CRC Press LLC FIGURE 5 across an interface is analogous to that across the two connected faces of the single sponge in the cylindrical case above — here it is simply that the "phantom" points just outside one sponge correspond to "real" points just inside a different sponge. Block-structured grid codes are based on this multiple-sponge analogy, with the physical field being filled with a group of grid blocks with correspondence of grid lines, and in fact complete continuity, across the interfaces between blocks. This approach has been carried to a high degree of application in the aerospace industry (cf. Chapter 13), with complete aircraft configurations being treated with a hundred or so blocks. Current grid generation systems seek to make the setup of this block structure both graphical and easy for the user. The ultimate goal is to automate the process (cf. Chapter 10). 2.3 Grid Generation Approaches With these obvious advantages of specialized curvilinear coordinate systems fitted to the boundaries of cylindrical, spherical, elliptical, and other analytic regions, it has been natural to use grids based on these systems for finite difference solutions on such regions. In the late 1960s the visual analogy between potential solutions (electrostatic potential, potential flow, etc.) that are solutions of Laplace's equation, ∇2φ = 0, and curvilinear grids led to the idea of generating grid lines in arbitrary regions as the solution of Laplace's equation. Thus, whereas potential flow is described in terms of a stream function ψ and a velocity potential φ that are orthogonal and satisfy ∇2ψ = 0, ∇2φ = 0 (Figure 6), a curvilinear grid could be generated by solving the system ∇2ξ = 0, ∇2η = 0 with η a constant on the upper and lower boundaries in the above region, while ξ is constant on the left and right boundaries (Figure 7). Here again, for purposes of programming, the grid forms a logically rectangular lattice (Figure 8). The problem of generating a curvilinear grid to fit any arbitrary region thus becomes a boundary value problem — the generation of interior values for the curvilinear coordinates from specified values on the boundary of the region (cf. Chapter 4). In order to set this up, note that we have for the boundary value problem the generation of interior values of the curvilinear coordinates ξ and η from specified constant values on opposing boundaries (Figure 9). Clearly ξ and η must vary monotonically and over the same range over the boundaries on which they are not specified, else the grid would overlap on itself. Thus, on the lower and upper boundaries, ξ here must vary monotonically from ξ 1 on the left to ξ 2 on the right. Similarly, on the left and right boundaries, ©1999 CRC Press LLC FIGURE 6 FIGURE 7 FIGURE 8 η must vary monotonically from η1 at the bottom to η2 at the top. The next question is what this variation should be. This is, in fact, up to the user. Ultimately, the discrete grid will be constructed by plotting lines of constant ξ and lines of constant η at equal intervals of each, with the size of the interval determined by the number of grid lines desired. Thus, if there are to be 10 grid lines running from side to side between the top and bottom of the region, 10 points would be selected on the left and right sides — with their locations being up to the user. Once these points are located, η can be said to assume, at the 10 points on each side, 10 values at equal intervals between its top and bottom values, η1 and η2. With this specification on the sides, the curvilinear coordinate η is thus specified on the entire boundary of ©1999 CRC Press LLC FIGURE 9 FIGURE 10 FIGURE 11 the region, and its interior values can be determined as a boundary value problem. A similar specification of ξ on the bottom and top boundaries by placing points on these boundaries sets up the determination of ξ in the interior from its boundary values. Now the problem can be considered a boundary value problem in the physical field for the curvilinear coordinates ξ and η (Figure 10) or can be considered a boundary value problem in the logical field for the Cartesian coordinates, x and y (Figure 11). Note that the boundary points are by nature equally spaced on the boundary of the logical field regardless of the distribution on the boundaries of the physical field. Continuing the potential analogy, the curvilinear grid can be generated by solving the system ∇2ξ = 0, ∇2η = 0, in the first case, or by solving the transformation of these equations (transformation relations are covered in Chapter 2), in the ©1999 CRC Press LLC FIGURE 12 αxξξ − 2 βxξη + γxηη = 0 αyξξ − 2 βyξη + γyηη = 0 α = xη2 + yη2 γ = xξ2 + yξ2 β = xξ xη + yξ yη second case. Although the equation set is longer in the second case, the solution region is rectangular, and the differencing can be done on a uniformly spaced rectangular grid. This is, therefore, the preferred approach. Note that the placing of points in any desired distribution on the boundary of the physical region, where x and y are the independent variables, amounts to setting (x,y) values at equally spaced points on the rectangular boundary of the logical field, where ξ and η are the independent variables. This is the case regardless of the shape of the physical boundary. This boundary value problem for the curvilinear grid can be generalized beyond the analogy with potential solutions, and in fact is in no way tied to the Laplace equation. The simplest approach is to generate the interior values by interpolation from the boundary values — a process called algebraic grid generation (cf. Chapter 3). There are several variants of this process. Thus for the region considered above, a grid could be generated by interpolating linearly between corresponding points on the top and bottom boundaries (Figure 12). Note that the point distributions on the side boundaries have no effect here. Alternatively, the interpolation could be between pairs of points on the side boundaries (Figure 13). The second case is, however, obviously unusable since the grid overlaps the boundary. Here the lack of influence from the points on the bottom boundary is disastrous. Another alternative is transfinite interpolation in which the interpolation is done in one (either) direction as above, but then the resulting error on the two sides not involved is interpolated in the other direction and subtracted from the first result. This procedure includes effects from all of the boundary and consequently matches the point distribution that is set on the entire boundary. This is the preferred approach, and it provides a framework for placing any one-dimensional interpolation into a multipledimensional form. It is possible to include any type of interpolation, such as cubic, which gives orthogonality at the boundaries, in the transfinite interpolation format. ©1999 CRC Press LLC FIGURE 13 It is still possible in some cases for the grid to overlap the boundaries with transfinite interpolation, and there is no control over the skewness of the grid. This gives incentive to now return to the grids generated from solving the Laplace equation. The Laplace equation is, by its very nature, a smoother, tending to average values at points with those at neighboring points. It can be shown from the calculus of variations, in fact, that grids generated from the Laplace equation are the smoothest possible. There arises, however, the need to concentrate coordinate lines in certain areas of anticipated strong solution variation, such as near solid walls in viscous flow. This can be accomplished by departing from the Laplace equation and designing a partial differential equation system for grid generation: designing because, unlike physics, there are no laws governing grid generation waiting to be discovered. The first approach to this, historically, was the obvious: simply replace the Laplace equation with Poisson equations ∇2ξ = P, ∇2η = Q and leave the control functions on the right-hand sides to be specified by the user (with appeal to Urania, the muse of science, for guidance). This does in fact work but the approach has evolved over the years, guided both by logical intuition and the calculus of variations, to use a similar set of equations but with a somewhat different right-hand side. Also, the user has been relieved of the responsibility for specifying the control functions, which are now generally evaluated automatically by the code from the boundary point distributions that are set by the user (cf. Chapter 4). These functions may also be adjusted by the code to achieve orthogonality at the boundary and/or to reduce the grid skewness or otherwise improve the grid quality (cf. Chapter 6). Algebraic grid generation, based on transfinite interpolation, is typically used to provide an initial solution to start an iterative solution of the partial differential equation for this elliptic grid generation system that provides a smoother grid, but with selective concentration of lines, and is less likely to result in overlapping of the boundary. This elliptic grid generation has an analogy to stretching a membrane attached to the boundaries (cf. Chapter 33) Grid lines inscribed on the underformed membrane move in space as the membrane is selectively stretched, but the region between the boundaries is always covered by the grid. Another form of grid generation from partial differential equations has an analogy with the waves emanating from a stone tossed into a pool This hyperbolic grid generation uses a set of hyperbolic equations, rather than the Poisson equation, to grow an orthogonal grid outward from a boundary (cf. Chapter 5). This approach is, in fact, faster than the elliptic grid generation, since no iterative solution is involved, but it is not possible to fit a specified outer boundary. Hyperbolic grid generation is thus limited in its use to open regions. As with the elliptic system, it is possible to control the spacing of the grid lines, and the orthogonality helps prevent skewness. ©1999 CRC Press LLC The control of grid line spacing can be extended to dynamically couple the grid generation system with the physical solution to be performed on the grid in order to resolve developing gradients in the solution wherever such variations appear in the field (cf. Chapter 34 and 35). With such adaptive grids, certain solution variables, such as pressure or temperature, are made to feed back to the control functions in the grid generations system to adjust the grid before the next cycle of the physical solution algorithm on the grid. P-2.4 Variations Structured grids today are typically generated and applied in the block-structured form described above with the multiple-sponge analogy. A variation is the chimera (from the monster of Greek mythology, composed of disparate parts) approach in which separate grids are generated about various boundary components, e.g., bodies in the field, and these separate grids are simply overlaid on a background grid and perhaps on each other in a hierarchy (cf. Chapter 11). The physical solution process on this composite grid proceeds with values being transferred between grids by interpolation. This approach has a number of advantages: (1) simplicity in grid generation since the various grids are generated separately, (2) bodies can be added to, or taken out of, the field easily, (3) bodies can carry their grids when moving relative to the background (think of simulating the kicking of a field goal with the ball and its grid tumbling end over end), (4) the separate grids can be used selectively to concentrate points in regions of developing gradients that may be in motion. The disadvantages are the complexity of setup (but this is being attacked in new code development) and the necessity for the interpolation between grids. Another approach of interest is the hybrid combination with separate structured grids over the various boundaries, connected by unstructured grids (cf. Chapter 23). There is great incentive to use structured grids over boundaries in viscous flow simulation because the boundary layer requires very small spacing out from the wall, resulting either in very long skewed triangular cells or a prohibitively and unnecessarily large number of small cells when unstructured grids are used. This hybrid approach is less well developed but can be expected to receive more attention. P-2.5 Transformation The use of numerically generated nonorthogonal curvilinear grids in the numerical solution of PDEs is not, in principle, any more difficult than using Cartesian grids: the differencing and solution techniques can be the same; there are simply more terms in the equations. For instance, the first derivative fx could be represented in difference form on a Cartesian grid as ( f x )ij = fi +1, j − fi −1, j 2∇x or if the spacing is not uniform, though the grid is still rectangular, by ( f x )ij = ©1999 CRC Press LLC fi +1, j − fi −1, j xi +1, j − xi −1, j To use a curvilinear grid, this derivative is transformed so that the curvilinear coordinate (ξ,η) rather than the Cartesian coordinate x,y, are the independent coordinates. Thus fx = ( xξ fη − xη fξ ) J where J = xξ yη – xη yξ is the Jacobian of the transformation and represents the cell volume. This then could be written in a difference form, taking ∆ξ and ∆η to be unity without loss of generality, using (f ) ξ ij ( fη )ij ( ) ( ) 1 i +1, j − i −1, j 2 1 fi +1, j − fi, j −1 = 2 = with analogous expressions for xξ , xη , yξ , yη. Movement of the grid, either to follow moving boundaries or to dynamically adapt to developing solution gradients, is not really a complication, since the time derivative can also be transformed as ( ft )r = ( ft )ξ − ( f x x + f y y˙ ) where the time derivative on the left is taken at a fixed position in space, i.e., is the time derivative appearing in the PDEs while the one on the right is that seen by a particular grid point moving with a speed ( x˙ , y˙ ). The spatial derivatives (fx , fy ) are transformed as was discussed above. There is no need to interpolate solution values from the grid at one time step to the displaced grid at the next time step, since that transfer is accomplished by the grid speed terms ( x˙ , y˙ ) in the above transformation relation. The straightforwardness of the use of curvilinear grids is further illustrated by the appearance of the generic convection–diffusion equations; ft + ∇ ⋅ (uf ) + ∇ ⋅ (v∇f ) + S = O where u is the velocity, v is a diffusion coefficient, and S is a source term, after transformation: 3 ( i =1 ) 3 3 ( ) At + Σ U i + v∇ 2ξ i Aξ i + Σ Σ g ij vAξ j i =1 j =1 3 ξi + A Σ a i ⋅ uξ i + S = 0 i =1 where now the time derivative is understood to be that seen by a certain (moving) grid point. Here the elements of the contravariant metric tensor g ij are given by g ij = a i ⋅ a j ©1999 CRC Press LLC where the ai are the contravariant base vectors (which are simply normals to the cell sides): ( ) a i = a j × ak / g (i, j,k cyclic) with the ai the covariant base vectors (tangents to the coordinate lines): ai = rξ i g is the Jacobian of the transformation (the where r is the Cartesian coordinate of a grid point, and cell volume): g = a1 ⋅ ( a2 × a3 ) Also, the contravariant velocity (normal to the cell sides) is U i = a ⋅ (u − r ) where u is the fluid velocity and r is the velocity of the moving grid. For comparison, the Cartesian grid formulation is 3 3 ( ) 3 At + Σ ui Ax i + Σ Σ δ ij vAx j i =1 i =1 j =1 3 xi + A Σ (ui ) x + S = 0 i =1 i The formulation has thus been complicated by the curvilinear grid only in the sense that the coefficient ui has been replaced by the coefficient U i + v(∇2ξ i ), and the Kronecker delta in the double summation has been replaced by g ij (thus expanding that summation from three terms to nine terms), and through the insertion of variable coefficients in the last summation. When it is considered that the transformed equation is to be solved on a fixed rectangular field with a uniform square grid, while the original equation would have to be solved on a field with moving curved boundaries, the advantages of using the curvilinear systems are clear. These advantages are further evidenced by consideration of boundary conditions. In general, boundary conditions for the example being treated would be of the form αA + βn ⋅ (u∇A) = γ where n is the unit normal to the boundary and α, β, and γ are specified. These conditions transform to αA + β ©1999 CRC Press LLC v g ii 3 Σ gij Aξ j =1 j =γ for a boundary on which ξ i is constant. For comparison, the original boundary conditions can be written in the form 3 αA + βv Σ n j Ax j = γ i =1 The transformed boundary conditions thus have the same form as the original conditions, but with the coefficient nj replaced by g ij/ g ii . The important simplification is the fact that the boundary to which the transformed conditions are applied is fixed and flat (coincident with a curvilinear coordinate surface). This permits a discrete representation of the derivative Aξ j along the transformed boundary without the need for interpolation. By contrast, the derivative Ax j in the original conditions cannot be discretized along the physical boundary without interpolation since the boundary is curved and may be in motion. Although the transformed equation clearly contains more terms, the differencing is the same as on a rectangular grid, i.e., it is done on the logically rectangular computational lattice, and the solution field is logically rectangular. Note that it is not necessary to discover and implement a transformation for each new boundary shape — rather the above formulation applies for all, simply with different values of (x, y, z) at the grid points. The transformed PDE can also be expressed in conservative form as ( 3 3 g A + Σ g U i A + v Σ g ij Aξ j + gS = 0 t i = 1 i =1 ξ i ) for use in the finite volume approach. For more information on transformations, see Chapter 2. P-3 P-3.1 Unstructured Grids Connectivities and Data Structures The basic difference between structured and unstructured grids lies in the form of the data structure which most appropriately describes the grid. A structured grid of quadrilaterals consists of a set of coordinates and connectivities that naturally map into elements of a matrix. Neighboring points in a mesh in the physical space are the neighboring elements in the mesh matrix (Figure 14). Thus, for example, a two-dimensional array x(i,j) can be used to store the x-coordinates of points in a 2D grid. The index i can be chosen to describe the position of points in one direction, while j describes the position of points in the other direction. Hence, in this way, the indices i and j represent the two families of curvilinear lines. These ideas naturally extend to three dimensions. For an unstructured mesh the points cannot be represented in such a manner and additional information has to be provided. For any particular point, the connection with other points must be defined explicitly in the connectivity matrix (Figure 15). ©1999 CRC Press LLC FIGURE 14 FIGURE 15 A typical form of data format for an unstructured grid in two dimensions is Number of Points, Number of Elements x1, y1 x2, y2 x3, y3 … n1, n2, n3 n4, n5, n6 n7, n8, n9 … where (x1, y1) are the coordinates of point i, and ni, 1=1,N are the point numbers with, for example, the triad (n1, n2, n3) forming a triangle. Other forms of connectivity matrices are equally valid, for example, connections can be based upon edges. The real advantage of the unstructured mesh is, however, because the points and connectivities ©1999 CRC Press LLC do not possess any global structure. It is possible, therefore, to add and delete nodes and elements as the geometry requires or, in a flow adaptivity scheme, as flow gradients or errors evolve. Hence the unstructured approach is ideally suited for the discretization of complicated geometrical domains and complex flowfield features. However, the lack of any global directional features in an unstructured grid makes the application of line sweep solution algorithms more difficult to apply than on structured grids. P-3.2 Grid Generation Approaches In contrast to the generation of structured grids, algorithms to construct unstructured grids are frequently based upon geometrical ideas. There are now many techniques available, many of which are described within this Handbook. For this elementary overview it is not appropriate to discuss details but to comment on general procedures. P-3.2.1 Triangle and Tetrahedra Creation by Delaunay Triangulation The Delaunay approach to unstructured grid generation is now popular. The basic concepts go back as far as Dirichlet, who in a paper in 1850 discussed the basic geometrical concept. Dirichlet proposed a method whereby a given domain could be systematically decomposed into a set of packed convex polygons. Given two points in the plane, P and Q, the perpendicular bisector of the line joining the two points subdivides the plane into two regions, V and W. The region V is closer to P than it is to Q. Extending these ideas, it is clear that for a given set of points in the plane, the regions Vi are territories that can be assigned to each point so that Vi represents the space closer to Pi than to any other point in the set. This geometrical construction of tiles is known as the Dirichlet tessellation. This tessellation of a closed domain results in a set of non-overlapping convex polygons, called Voronoï regions, covering the entire domain. From this description, it is apparent that in two dimensions, the territorial boundary that forms a side of a Voronoï polygon must be midway between the two points it separates and is thus a segment of the perpendicular bisector of the line joining these two points. If all point pairs that have some segment of a boundary in common are joined by straight lines, the result is a triangulation of the convex hull of the set of points Pi. This triangulation is known as the Delaunay triangulation. Equivalent constructions can be defined in higher dimensions. In three dimensions, the territorial boundary that forms a face of a Voronoï polyhedron is equidistant between the two points it separates. If all point pairs that have a common face in the Voronoï construction are connected, then a set of tetrahedra is formed that covers the convex hull of the data points. For the number of points which may be required in grid for computational analysis, it might appear that the above procedure would be difficult and computationally expensive to construct. However, there are several algorithms that can form the construction in a very efficient manner. These are discussed at length in Chapters 1, 16 and 20. The approach is very flexible in that it can automatically create grids with the minimum of user interaction for arbitrary geometries. P-3.2.2 Triangle and Tetrahedra Creation by the Advancing Front Method A grid generation technique based on the simultaneous point generation and connection is the advancing front method. Unlike the Delaunay approach, advancing front methods are not based on any geometrical criteria. They encompass the logical procedure of starting with a boundary grid of edges, in two dimensions, triangular faces, in three dimensions, and creating a point and constructing an element. Slowing the initial boundary advances into the domain until the domain is filled with elements. The placing of ©1999 CRC Press LLC points within the domain is, like the Delaunay approach, controlled by a combination of a background mesh and sources that provides the required data to ensure adequate resolution of the domain. The algorithms that generate grids in this way are based on fast geometrical search routines. Details are to be found in Chapter 17. It is possible to combine techniques from both the Delaunay and the Advancing Front methods to produce effective grid generation procedures – a sort of combination that tries to utilize the advantages of both approaches. Chapter 18 discusses one such approach. The Delaunay triangulation produces elements that are isotropic in nature. Although the Advancing Front method can produce elements with stretching, it cannot produce high quality meshes with stretching factors applicable to some problems, such as high Reynolds number viscous flows. Hence, it is necessary to augment the standard procedures outlined above. In general, this is done by introducing a mapping that ensures that regular isotropic grids can be generated but once mapped back to the physical space are distorted in a well defined manner to give appropriate element stretching. Such a method is described in detail in Chapter 20. P-3.2.3 Unstructured Grids of Quadrilaterals and Hexahedra The preference of some developers for quadrilateral or hexahedral element based unstructured meshes has resulted in effort devoted to the generation of such meshes. In two dimensions, it is possible to modify the Advancing Front algorithm to construct quadrilaterals, although the additional complexity in extending this approach to three dimensions has not yet been overcome for practical geometries. An alternative approach that has seen some success is that of "paving." This approach relies upon iteratively layering or paving rows of elements in the interior of a region. As rows overlap or coincide they are carefully connected together. It is fair to conclude that almost without exception the methods for the construction of unstructured hexahedral based grids are heuristic in nature, requiring considerable effort to include the many possible geometrical occurrences. Chapter 21 discusses in detail aspects of this kind of grid generation. P-3.2.4 Surface Mesh Generation The generation of unstructured grids on surfaces is, in itself, one of the most difficult and yet important aspects of mesh generation in three dimensions. The surface mesh influences the field mesh close to the boundary. Surface meshes have the same requirement for smoothness and continuity as the field meshes for which they act as boundary conditions, but in addition, they are required to conform to the geometry surfaces, including lines of intersection and must accurately resolve regions of high curvature. The approach usually taken to generate grids on surfaces is to represent the geometry in parametric coordinates. A parametric representation of a surface is straightforward to construct and provides a description of a surface in terms of two parametric coordinates. This is of particular importance, since the generation of a mesh on a surface then involves using grid generation techniques developed for two space dimensions. A full description of these procedures is given in Chapter 19. P-3.3 Grid Adaptation Techniques To resolve features of a solution field accurately it is, in general, necessary to introduce grid adaptivity techniques. Adaptivity is based on the equidistribution of errors principle, namely, wi dsi = constant ©1999 CRC Press LLC where wi is the error or activity indicator at node i and dsi is the local grid point spacing at node i. Central to adaptivity techniques and the satisfaction of this equidistribution principle is to define an appropriate indicator wi. Adaptivity criteria are based on an assessment of the error in the solution of the governing equations or are constructed to detect features of the field. These estimators are intimately connected to the analysis equations to be solved. For example, some of the main features of a solution of the Euler equations can be shock waves, stagnation points and vortices, and any indicator should accurately identify these flow characteristics. However, for the Navier-Stokes equations, it is important not only to refine the mesh in order to capture these features but, in addition, to adequately resolve viscous dominated phenomena such as the boundary layers. Hence, it seems likely that, certainly in the near future, adaptivity criteria will be a combination of measures, each dependent on some aspects of the flow and, in turn, on the flow equations. There is also an extensive choice of criteria based on error analysis. Such measures include, a comparison of computational stencils of different orders of magnitude, comparison of the same derivatives on different meshes, e.g., Richardson extrapolation, and resort to classical error estimation theory. No generally applicable theory exists for errors associated with hyperbolic equations, hence, to date combinations of rather ad hoc methods have been used. Once an adaptivity criterion has been established, the equidistribution principle is achieved through a variety of methods, including point enrichment, point derefinement, node movement and remeshing, or combinations of these. For more information on grid adaption techniques, see Chapter 35. P-3.3.1 Grid Refinement Grid refinement, or h-refinement, involves the addition of points into regions where adaptation is required. Such a procedure clearly provides additional resolution at the expense of increasing the number of points in the computation. Grid refinement on unstructured grids is readily implemented. The addition of a point or points involves a local reconnection of the elements, and the resulting grid has the same form as the initial grid. Hence, the same solver can be used on the enriched grid as was used on the initial grid. It is important that the adaptivity criteria resolve both the discontinuous features of the solution (i.e. shock waves, contacts) and the smooth features as the number of grid points are increased. A desirable feature of any adaptive method to ensure convergence is that the local cell size goes to zero in the limit of an infinite number of mesh points. Grid refinement on a structured or multiblock grids is not so straightforward. The addition of points will, in general, break the regular array of points. The resulting distributed grid points no longer naturally fit into the elements of an array. Furthermore, some points will not "conform" to the grid in that they have a different number of connections to other points. Hence grid refinement on structured grids requires a modification to the basic data structure and also the existence of so-called non-conforming nodes requires modifications to the solver. Clearly, point enrichment on structured grids is not as natural a process as the method applied on unstructured grids and hence is not so widely employed. Work has been undertaken to implement point enrichment on structured grids and the results demonstrate the benefits to be gained from the additional effort in modifications to the data structure and the solve. P-3.3.2 Grid Movement Grid movement satisfies the equidistribution principle through the migration of points from regions of low activity into regions of high activity. The number of nodes in this case remains fixed. Traditionally, algorithms to move points involve some optimization principle. Typically, expressions for smoothness, ©1999 CRC Press LLC orthogonality and weighting according to the analysis field or errors are constructed and then an optimization is performed such that movement can be driven by a weight function, but not at the expense of loss of smoothness and orthogonality. Such methods are in general, applicable to both structured and unstructured grids. An alternative approach is to use a weighted Laplacian function. Such a formulation is often used to smooth grids, and of course the formal version of the formulation is used as the elliptic grid generator presented earlier. P-3.3.3 Combinations of Node Movement, Point Enrichment and Derefinement An optimum approach to adaptation is to combine node movement and point enrichment with derefinement. These procedures should be implemented in a dynamic way, i.e., applied at regular intervals within the simulation. Such an approach also provides the possibility of using movement and enrichment to independently capture different features of the analysis. P-3.3.4 Grid Remeshing One method of adaptation which, to date, has been primarily used on unstructured grids, is adaptive remeshing. As already indicated, unstructured meshes can be generated using the concept of a background mesh. For an initial mesh, this is usually some very coarse triangulation that covers the domain and on which the spatial distribution is consistent with the given geometry. For adaptive remeshing, the solution achieved on an initial mesh is used to define the local point spacing on the background mesh which was itself the initial mesh used for the simulation. The mesh is regenerated using the new point spacing on the background mesh. Such an approach can result in a second adapted mesh that contains fewer points than that contained in the initial mesh. However, there is the overhead of regeneration of the mesh which in three dimensions can be considerable. Nevertheless, impressive demonstrations of its use have been published. ©1999 CRC Press LLC Contents Foreword Contributors Acknowledgments Preface: An Elementary Introduction Joe F. Thompson, Bharat K. Soni, and Nigel P. Weatherill 1 Fundamental Concepts and Approaches Joe F. Thompson and Nigel P. Weatherill PART I Block-Structured Grids Introduction to Structured Grids Joe F. Thompson 2 Mathematics of Space and Surface Grid Generation Zahir U.A. Warsi 3 Transfinite Interpolation (TFI) Generation Systems Robert E. Smith 4 Elliptic Generation Systems 5 Hyperbolic Methods for Surface and Field Grid Generation Stefan P. Spekreijse William M. Chan 6 Boundary Orthogonality in Elliptic Grid Generation Ahmed Khamayseh, Andrew Kuprat, and C. Wayne Mastin 7 Orthogonal Generation Systems 8 Harmonic Mappings 9 Surface Grid Generation Systems 10 Luís Eça Sergey A. Ivanenko Ahmed Khamayseh and Andrew Kuprat A New Approach to Automated Multiblock Decomposition for Grid Generation: A Hypercube++ Approach Sangkun Park and Kunwoo Lee ©1999 CRC Press LLC 11 12 Composite Overset Structured Grids Robert L. Meakin Parallel Multiblock Structured Grids Jochem Häuser, Peter R. Eiseman, Yang Xia, and Zheming Cheng 13 Block-Structured Applications PART II Timothy Gatzke Unstructured Grids Introduction to Unstructured Grids Nigel P. Weatherill 14 Data Structures for Unstructured Mesh Generation 15 Automatic Grid Generation Using Spatially Based Trees Luca Formaggia Mark S. Shephard, Hugues L. de Cougny, Robert M. O'Bara, and Mark W. Beall 16 Delaunay–Voronoï Methods 17 Advancing Front Grid Generation 18 Unstructured Grid Generation Using Automatic Point Insertion and Local Reconnection David L. Marcum 19 Surface Grid Generation 20 Nonisotropic Grids 21 Quadrilateral and Hexahedral Element Meshes 22 Timothy J. Baker J. Peraire, J. Peiró, and K. Morgan J. Peiró Paul Louis George and Frédéric Hecht Adaptive Cartesian Mesh Generation Robert Schneiders Michael J. Aftosmis, Marsha J. Berger, and John E. Melton 23 Hybrid Grids 24 Parallel Unstructured Grid Generation Jonathon A. Shaw Hugues L. de Cougny and Mark S. Shephard 25 Hybrid Grids and Their Applications 26 Unstructured Grids: Procedures and Applications ©1999 CRC Press LLC Yannis Kallinderis Nigel P. Weatherill PART III Surface Definition Introduction to Surface Definition Bharat K. Soni 27 Spline Geometry: A Numerical Analysis View 28 Computer-Aided Geometric Design 29 Computer-Aided Geometric Design Techniques for Surface Grid Generation Bernd Hamann, Brian Jean, and Anshuman Razdan 30 NURBS in Structured Grid Generation 31 David R. Ferguson Gerald Farin Tzu-Yi Yu and Bharat K. Soni NASA IGES And NASA-IGES NURBS Only Standard Austin L. Evans and David P. Miller PART IV Adaptation and Quality Introduction to Adaptation and Quality Bharat K. Soni 32 Truncation Error on Structured Grids 33 Grid Optimization Methods for Quality Improvement and Adaptation Olivier-Pierre Jacquotte 34 Dynamic Grid Adaptation and Grid Quality C. Wayne Mastin D. Scott McRae and Kelly R. Laflin 35 Grid Control and Adaptation 36 Variational Methods of Construction of Optimal Grids O. Hassan and E. J. Probert A. F. Sidorov, and O. V. Ushakova 37 Moving Grid Techniques Paul A. Zegeling Appendix A: Grid Software and Configurations Appendix B: Grid Configurations ©1999 CRC Press LLC Bharat K. Soni Bharat K. Soni O. B. Khairullina, I Block-Structured Grids Joe F. Thompson Introduction to Structured Grids The grid generation process, in general, proceeds from first defining the boundary geometry as discussed in Part III. Then points are distributed on the curves that form the edges of boundary sections. A surface grid is then generated on the boundary surface, and finally a volume grid is generated in the field. Chapter 13 gives a general overview of the entire grid generation process and the fundamental choices and considerations involved from the standpoint of the user. The underlying essential mathematics of structured grid generation, including essential concepts from differential geometry and tensor analysis, is collected in Chapter 2. The mathematical constructs explained in this chapter are utilized throughout the chapters of this handbook. The distribution of points on boundary curves (edges of boundary surfaces) is commonly done through several distribution functions as described in Section 3.6 of Chapter 3. (The mathematics of curves is covered in Section 2.3 of Chapter 2.) These functions have been adopted over time as providing point distributions that comply with certain constraints that must be applied in order to control error that can be introduced into the solution by the grid if the spacing changes too rapidly, as discussed in Chapter 32 of Part IV. Structured grids can be generated algebraically or as the solution of PDEs. Algebraic grid generation is simply some form of interpolation from boundary points — the variants just use different kinds of interpolation. The most fundamental and versatile form — and now commonly incorporated in grid generation codes — is TFI (transfinite interpolation), which is introduced in Section 1.3.5 of Chapter 1 and described in Chapter 3. The basic equations of TFI are given in Section 3.4 of Chapter 3, and the specific equations for application with and without orthogonality at the boundaries are given in Section 3.5. Algebraic grid generation based on TFI is the fastest procedure for structured grids, and is also commonly used to generate an initial grid in generation systems based on PDEs. Grids generated algebraically can, however, have some problems with smoothness and may overlap strongly convex portions of boundaries. ©1999 CRC Press LLC Generation systems based on PDEs can produce smoother grids with fewer problems with boundary overlap. Such generation systems are therefore often used to smooth algebraic grids. Since grid generation is essentially a boundary-value problem, grids can be generated from point distributions on boundaries by solving elliptic PDEs in the field. The smoothness properties and extremum principles inherent in some such PDE systems can serve to produce smooth grids without boundary overlap. The PDE solution is generally one by iteration, and therefore elliptic grid generation is not as fast as algebraic grid generation. The elliptic PDEs for grid generation are not unique, of course, but must be designed. This design has converged over the years to the elliptic system given in Section 1.3.3 of Chapter 1, which forms the basis for most grid generation codes today. This formulation incorporates control functions that are determined from the boundary point distribution to control the grid line spacing and orientation in the field to be compatible with that on the boundary. Procedures for the determination of these control functions in grid codes have evolved in time to the forms noted in this section of Chapter 1, which can accomplish boundary orthogonality through iterative adjustment during the generation process. A more recent and general formulation, with a sounder basis for evaluation of the control functions, is given here in Chapter 4: for 2D in Section 4.2 and for 3D in Section 4.4. This iterative solution of the elliptic system is often done by SOR, but a Picard iteration is given in Section 4.2.2 of Chapter 4, and a conjugate gradient solution is given in Section 12.10.4 of Chapter 13, in connection with parallel implementation. The generation of a grid on a boundary surface is a necessary prelude to the generation of a volume grid, and this is generally done by representing the boundary surface parametrically by NURBS or another spline formulation, and then generating the grid in parameter space either algebraically or using PDEs. This is perfectly analogous to 2D grid generation except that surface curvature terms appear in the PDEs. With the generation system operating in parameter space, the resulting grid is guaranteed to lie on the boundary surface. The parametric representation of the boundary surface is covered in Chapter 29, utilizing the underlying curve and surface constructs given in Chapter 28. Other aspects of surface generation are covered in the other chapters in Part III, and the mathematical foundations are given in Section 2.4 and in Section 2.5.2 of Chapter 2. Algebraic surface grid generation is simply the application of TFI to generate values of the surface parameters on the surface from the values set on the edges of the boundary surface by the grid point distribution on those edges, as covered in Section 9.2 of Chapter 9. Elliptic surface grid generation operates with the PDEs formulated in terms of the surface parameters, and surface curvature terms appearing in the PDEs (see Section 2.5.2 of Chapter 2). A commonly applied procedure is given in Section 9.3 of Chapter 9, and a more recent and general procedure is given in Section 4.3 of Chapter 4. Hyperbolic surface grid generation is covered in Section 5.3 of Chapter 5. It is generally advantageous, in view of such things as boundary layer phenomena and turbulence models, to have the grid orthogonal to boundaries even though orthogonality is not imposed in the field. This is commonly done through iterative adjustment of the control functions as described in Chapter 6: in Section 6.2 for 2D grids, Section 6.3 for surface grids, and Section 6.4 for volume grids. Another procedure in 2D, also using the control functions, is given in Section 4.2 of Chapter 4. An alternative approach to grid generation via PDEs is to use a hyperbolic generation system rather than an elliptic. Elliptic equations admit boundary conditions, i.e., grid point distributions, on all boundaries of a region. Hyperbolic systems, however, can take boundary conditions only on a portion of the boundary. Therefore, while elliptic grid generation systems can produce a grid in the entire volume from point distributions of the entire boundary, hyperbolic systems generate the grid by marching outward from a portion of the boundary. Hyperbolic grid generation systems therefore cannot be used to generate a grid in the entirety of a volume defined by a complete boundary. Chapter 5 covers hyperbolic grid generation in volumes in Section 6.2 and on surfaces in Section 6.3. Structured grids are not generally made orthogonal, although orthogonality at boundaries is often incorporated, as has been noted above. In fact, 3D orthogonality is not, in general, possible without imposing certain conditions on the grids on the boundary surfaces. And even in 2D, orthogonality imposes severe restrictions on the grid distribution. Transformed PDEs, however, take a much simpler ©1999 CRC Press LLC form on orthogonal grids, providing some incentive for their use when feasible — with relatively simple boundary configurations and physical problems without strong localized gradients. Chapter 7 covers orthogonal grid generation systems. As has been noted, PDEs for grid generation are designed, not discovered. Considerable research has gone into this topic, leading to generally standard elliptic (Chapter 4) and hyperbolic (Chapter 5) grid generation systems. The underlying theory of harmonic mappings provides a framework for the development of elliptic grid generation systems, and this topic is treated in some depth in Chapter 9. This theoretical base also leads to the formulation of adaptive grid systems, also covered in this chapter. Adaptive grids are most fundamentally formulated from variational principles, and this is covered in Chapter 36 of Part IV. Adaptive grids and grid quality are covered in the chapters of Part IV. A strong and versatile alternative to block-structured grids is the overset grid approach (originally called chimera, after the composite monster of Greek mythology). With this approach, individual structured grids are generated around separate boundary components, e.g., bodies, and these separate grids simply overlap each other in some hierarchy. Data is transferred between overlapping grids by interpolation. The overset grid approach is covered here in Chapter 11. The grid generation involved is typically done by hyperbolic generation systems, described in Chapter 5. The mathematics and technology of structured grid generation have matured now so that the techniques covered in Part I can be expected to be of enduring utility. The block structure is versatile, and serves as the foundation for efficient solutions because of its inherently simple data structure. Construction of the block configuration by hand, even with graphically interactive tools, is very labor intensive, however, as noted in Chapter 13. Automation of the block structure, rather than graphical interaction, is the goal, and this is an area of active research and development (Section 21.2 of Chapter 21 is relevant here). A very promising recent approach is included in Chapter 11. Finally, operation on parallel processors is essential now, and the block structure provides a natural means of domain decomposition, as covered in Sections 12.8–12.10 of Chapter 12. The operation of the block structure is discussed in Sections 12.2–12.6 of Chapter 12. Chapter 12 also covers a script-based meta-language approach to structured grid generation in Section 12.7. Although most available grid generation systems have departed from the script-based approach in favor of graphical interaction, the script-based approach has definite advantages in design cycles. ©1999 CRC Press LLC 1 Fundamental Concepts and Approaches 1.1 1.2 1.3 Introduction Mesh Generation Considerations Structured Grids Composite Grids • Block-Structured Grids • Elliptic Systems • Hyperbolic System • Algebraic System • Adaptive Grid Schemes 1.4 Joe F. Thompson Nigel P. Weatherill Unstructured Grids The Delaunay Triangulation • Point Creation • Other Unstructured Grid Techniques • Unstructured Grid Generation on Surfaces • Adaptation on Unstructured Grids • Summary 1.1 Introduction This introductory chapter uses fluid mechanics as an example field problem for reference; the applicability of the concepts discussed is, however, not in any way limited to this area. Fluid mechanics is described by nonlinear equations, which cannot generally be solved analytically, but which have been solved using various approximate methods including expansion and perturbation methods, sundry particle and vortex tracing methods, collocation and integral methods, and finite difference, finite volume, and finite element methods. Generally the finite difference, finite volume, and finite element discretization methods have been the most successful, but to use them it is necessary to discretize the field using a grid (mesh). (The terms grid and mesh are used interchangeably throughout with identical meaning.) The mesh can be structured or unstructured, but it must be generated under some of the various constraints described below, which can often be difficult to satisfy completely. In fact, at present it can take orders of magnitude more person-hours to construct the grid than it does to construct and analyze the physical solution on the grid. This is especially true now that solution codes of wide applicability are becoming available. Computational fluid dynamics (CFD) is a prime example, and grid generation has been cited repeatedly as a major pacing item (cf. Thompson [1996]). The same is true for other areas of computational field simulation. The proceedings of the several international conferences on grid generation (Thompson [1982], Hauser and Taylor [1986], Sengupta, et al. [1988], Arcilla, et al. [1991], Eiseman, et al. [1994], Soni et al. [1996]) as well as those of the NASA conferences (Smith [1980], Smith [1992], Choo [1995]) provide numerous illustrations of application to CFD and some other fields. A recent comprehensive text is Carey [1997]. ©1999 CRC Press LLC 1.2 Mesh Generation Considerations The generated mesh must be sufficiently dense that the numerical approximation is an accurate one, but it cannot be so dense that the solution is impractical to obtain. Generally, the grid spacing should be smoothly and sufficiently refined to resolve changes in the gradients of the solution. If the grid is also boundary-conforming and curvilinear, the application of boundary conditions is simplified. Boundaryconforming curvilinear grids may also allow the use of various approximate equations such as boundarylayer equations. The grid should also be constructed with computational efficiency in mind. The accuracy of a numerical approximation can also be impaired, if a grid changes discontinuously or is too skewed. Various computers often require well-organized data, and memory requirements can grow to impractical limits unless the data is organized well. Finally, the choice of a grid should not lead to overly complex computer codes. A mesh is a set of points distributed over a calculation field for a numerical solution of a set of partial differential equations (PDEs). This set may be structured, e.g., formed by the intersections of curvilinear coordinate surfaces, or unstructured, i.e., with no relation to coordinate directions. In the first case the points form quadrilateral cells in 2D, or hexahedral cells in 3D (with nonplanar sides). The unstructured mesh generally consists of triangles and tetrahedra in 2D and 3D, respectively, in its most basic form, but may be made of hexahedra or elements of any shape in general. The structured grid can be generated algebraically by interpolation from boundaries, e.g., transfinite interpolation, or by solving a set of partial differential equations in the region. An entire subject, complete with textbook (Thompson, Warsi, and Mastin [1985], now on the Web at www.erc.msstate.edu), has developed around the generation of structured grids having the fundamental characteristic that some curvilinear coordinate surface is coincident with each boundary segment, i.e., boundary conforming. A later text is Knupp and Steinberg [1993]. Castillo [1991] provides a compilation of mathematical aspects as well. Structured grid generation is also covered in the recent text of Carey [1997]. Several earlier surveys of the field are still useful for basic understanding (Thompson, Warsi, and Mastin [1982], Thompson [1984], Thompson [1985], Eiseman [1985]). Structured boundary-conforming meshes have been widely applied in computational fluid dynamics. Basically, the algebraic generation systems (Chapter 3) are faster, but the grids generated from partial differential equations are generally smoother. The hyperbolic generation systems (Chapter 5) are faster than the elliptic systems, but are more limited in the geometries that can be treated. The elliptic systems (Chapter 4) are the most generally applicable with complicated boundary configurations, but transfinite interpolation is also effective in a composite grid framework. The generation of unstructured meshes can be done by tessellation of a point distribution that could be random but is more likely to have been produced by some ordered procedure. This tessellation is not unique, and involves some type of nearest-neighbor search, such as the Delaunay triangulation (Chapter 16). Other approaches are the advancing front procedure (Chapter 17) and the finite octree method (Chapter 15). The recent text of Carey [1997] covers unstructured grid generation as well as structured grid generation. An earlier text on unstructured grids is George [1991]. General configurations can conceivably be treated with either type of mesh, and hybrid combinations (Chapter 23) are also possible, using individual structured meshes near boundaries, with these subregions being connected by an unstructured mesh. Still another approach is overlaid grids (chimera) (Chapter 11), in which separate boundary-conforming structured grids are generated for each component of a complex configuration, and data is communicated between the various component grids by interpolation. Of particular importance is the development of dynamically adaptive meshes (Chapters 33–36) coupled with the physical solution. In this mode the mesh is locally refined by the selective addition of points, and/or is moved to concentrate points, in order to resolve developing gradients in the physical solution on the mesh. Both of these approaches have seen considerable development and show much promise in particular areas. Implementation of solution algorithms on structured, unstructured, and overlaid grids places differing requirements on the algorithms. Various conflicts arise between the grid and solution procedures in ©1999 CRC Press LLC regard to requirements and ease of operation. In particular, unstructured grids require a much more complex solution data structure, but are more easily generated and adapted. Structured grids provide a more natural representation of normal derivative boundary conditions and allow more straightforward approximations based on prevailing directions, e.g., parallel or normal to a boundary or flow direction. The structure also leads to a much simpler data set construction, and allows the use of directional time splitting and flux representations. On the other hand, unstructured grids can be much more readily imagined for complicated boundary configurations. 1.3 Structured Grids 1.3.1 Composite Grids Structured grid generation had its roots in the U.S. in the work of Winslow and Crowley at Lawrence Livermore National Lab in the late 1960s (Winslow [1967]), and in Russia from Godunov and Prokopov [1967] at about the same time. (There is also that enigmatic Biblical reference to the "four corners of the earth," once thought to proclaim a flat earth but now seen to be a prophesy of structured grid generation.) Another very fundamental component was the work of Bill Gordon at Drexel on transfinite interpolation for the automotive industry, introduced to the emerging grid generation community at the grid conference in Nashville in 1982 (Gordon and Thiel [1982]). 184.108.40.206 Terminology The use of composite grids has been the key to the treatment of general 3D configurations with structured grids. Here in general, composite refers to the fact that the physical region is divided into subregions, within each of which a structured grid is generated. These subgrids may be patched together at common interfaces, may be overlaid, or may be connected by an unstructured grid. Considerable confusion has arisen in regard to terminology for composite grids, making it difficult to immediately classify papers on the subject. Composite grids in which the subgrids share common interfaces are referred to as block, patched, embedded, or zonal grids in the literature. The use of the first two of these terms is fairly consistent with this type of grid (patched comes from the common interfaces, block from the logically rectangular structure), but the last two are sometimes also applied to overlaid grids. Overlaid (overset) grids are often called chimera grids after the composite monster of Greek mythology. Unfortunately, the common interface grids can also be said to overlap, since they typically use surrounding layers of points to achieve continuity. Embedded grids can be almost anything, and the term is probably best avoided. The use of zonal comes mostly from CFD applications where the suggestion of applying different solution equations sets in different flow regions is made. Perhaps block or patched would be best for the common interface grids, chimera for the overlaid (avoiding overlaid) grids, and hybrid for the structured–unstructured combinations. 184.108.40.206 Forms With this terminology adopted, the block (or patched) grids may be completely continuous at the interfaces, have slope or line continuity, or be discontinuous (sharing a common interface but not common points thereon). (Block seems to cover all of these possibilities, but patched is being stretched a bit in the latter case.) Complete continuity is achieved through a surrounding layer of (image, phantom) points at which values are kept equal to those at corresponding object points inside an adjacent block. This requires a data indexing procedure to link the blocks across the interfaces. With complete continuity, the interface is not fixed (not even in shape), but is determined in the course of the solution. This type of interface necessitates an elliptic generation system. Slope continuity requires that the grid generation procedure incorporate some control over the intersection angle at boundaries (usually, but not necessarily, orthogonality), as can be done through Hermite interpolation in algebraic generation systems or through iterative adjustment of the control functions in elliptic systems. In this case the points on the interface are fixed, and the subgrids are generated independently, except for the use of the common interface ©1999 CRC Press LLC points and a common (presumably orthogonal) angle of intersection with the interface. The PDE coding construction is greatly simplified with either complete or slope continuity, since then no algorithm modifications are necessary at the interfaces. The chimera (overlaid) grids are composed of completely independent component grids that overlap a background grid, other component grids and/or other component boundary elements, creating holes in the component grids. This requires flagging procedures to locate grid points that lie out of the field of computation, but such holes can be handled even in tridiagonal solvers by placing ones at the corresponding positions on the matrix diagonal and all zeros off the diagonal. These overlaid grids also require interpolation to transfer data between grids, and that subject is the principal focus of effort in regard to the use of this type of composite grid. The hybrid structured–unstructured grids avoid this interpolation by replacing the overlaid region with an unstructured grid connecting logically rectangular structured component grids. This can require modification of solution codes, however. 1.3.2 Block-Structured Grids Block-structured grids opened the door to real-world CFD in the late 1980s, and many real applications are still based on these grids (see Chapters 12 and 13). The idea appears in the proceedings of the grid conference in Nashville in 1982 (Rubbert and Lee [1982]), but it was Weatherill and Forsey [1984], and Miki and Takagi [1984] that attracted attention to the block-structured approach. Today's structured grid codes are based on this approach. Although the grid is logically rectangular within each block, the blocks fit together in an unstructured manner. Block-structured generation systems that maintain complete continuity across block interfaces allow difference representations to be applied on the block interfaces as in the rest of the field. Complete continuity across block interfaces in the field is accomplished by treating the interface in the manner of a branch cut, with correspondence between "phantom" points outside one block with "real" points inside the adjacent block. The curvilinear grid system can be constructed simply by setting values in a rectangular array of position vectors, rijk (i = 1, 2,..., I j = 1, 2,..., J k = 1, 2,..., K ) and identifying the indices i, j, k with the three curvilinear coordinates. The position vector r is a threevector giving the values of the x, y, and z Cartesian coordinates of a grid point. Since all increments in the curvilinear coordinates cancel out of the transformation relations for derivative operators, there is no loss of generality in defining the discretization to be on integer values of these coordinates. Fundamental to this curvilinear coordinate system is the coincidence of some coordinate surface with each segment of the boundary of the physical region, in the same manner that surfaces of constant radius coincide with the inner and outer boundary segments of the region between two concentric spheres filled with a polar coordinate system. This is accomplished by placing a two-dimensional array of points on a physical boundary segment and entering these values into the array rijk of position vectors, with one index constant, e.g., in , rijk with i from 1 to I and j from 1 to J. The curvilinear coordinate k is thus constant on this physical boundary segment. With values set on the sides of the rectangular array of position vectors in this manner, the generation of the grid is accomplished by determining the values in its interior, e.g., by interpolation or a PDE solution. The set of values rijk then forms the nodes of a curvilinear coordinate system filling the physical region. A physical region bounded by six generally curved sides can thus be considered to have been transformed to a logically rectangular computational region on which the curvilinear coordinates are the independent variables. Although in principle it is possible to establish a correspondence between any physical region and a single empty logically rectangular block for general three-dimensional configurations, the resulting grid is likely to be much too skewed and irregular to be usable when the boundary geometry is complicated. A better approach with complicated physical boundaries is to segment the physical region into contiguous ©1999 CRC Press LLC subregions, each bounded by six curved sides (four in 2D) — each of which transforms to a logically rectangular block in the computational region. Each subregion has its own curvilinear coordinate system irrespective of that in the adjacent subregions (see Figure 13.5). This then allows both the grid generation and numerical solutions on the grid to be constructed to operate in a logically rectangular computational region, regardless of the shape or complexity of the full physical region. The full region is treated by performing the solution operations in all of the logically rectangular computational blocks. With the block-structured framework, partial differential equation solution procedures written to operate on logically rectangular regions can be incorporated into a code for general configurations in a straightforward manner, since the code only needs to treat a rectangular block. The entire physical field then can be treated in a loop over all the blocks. Transformation relations for partial differential equations are covered in Chapter 2 and in Thompson, Warsi, and Mastin [1985], on the Web. Discretization error related to the grid is covered in Chapter 32. The evaluation and control of grid quality (Chapter 33) is an ongoing area of active research. The generally curved surfaces bounding the subregions in the physical region form internal interfaces across which information must be transferred, i.e., from the sides of one logically rectangular computational block being paired with another on the same, or different, block, since both correspond to the same physical surface. Grid lines at the interfaces may meet with complete continuity, with or without slope continuity, or may not meet at all. Complete continuity of grid lines across the interface requires that the interface be treated as a branch cut on which the generation system is solved just as it is in the interior of blocks. The interface locations are then not fixed, but are determined by the grid generation system. This is most easily handled in coding by providing an extra layer of points surrounding each block. Here the grid points on an interface of one block are coincident in physical space with those on another interface of the same or another block, and also the grid points on the surrounding layer outside the first interface are coincident with those just inside the other interface, and vice versa. This coincidence can be maintained during the course of an iterative solution of an elliptic generation system by setting the values on the surrounding layers equal to those at the corresponding interior points after each iteration. All the blocks are thus iterated to convergence together, so that the entire composite grid is generated at once. The same procedure is followed by PDE solution codes on the block-structured grid. The construction of codes for complicated regions is greatly simplified by the block structure since, with the use of the surrounding layer of points on each block, a PDE code is only required basically to operate on a logically rectangular computational region. The necessary correspondence of points on the surrounding layers (image points) with interior points (object points) is set up by the grid code and made available to the PDE code. 1.3.3 Elliptic Systems Elliptic grid generation is treated in detail in Chapter 4. This section provides an overview of the technology as applied in the EAGLE system (Thompson [1987]), as an example of the technology applied in several current grid generation codes. 184.108.40.206 Generation System An elliptic grid generation system used in many codes is 3 3 3 Σ Σ g mn rξ mξ n + nΣ=1 g nn Pn rξ n = 0 m =1 n =1 where the gmn are the elements of the contravariant metric tensor, g mn = ∇ξ m ⋅ ∇ξ n ©1999 CRC Press LLC (1.1) These elements are more conveniently expressed for computation in terms of the elements of the covariant metric tensor, gmn, gmn = rξ m ⋅ rξ n which can be calculated directly. Thus ( g mn = gik g jl − gil g jk )/ g ( m, i, j ) cyclic, (n, k, l) cyclic where g, the square of the Jacobian, is given by ( g = det gmn = rξ 1 ⋅ rξ 2 × rξ 2 ) In these relations, r is the Cartesian position vector of a grid point (r = ix + jy + kz) and the ξ i (i = 1,2,3) are the three curvilinear coordinates. The Pn (n = 1,2,3) are the control functions that serve to control the spacing and orientation of the grid lines in the field. The first and second coordinate derivatives are normally calculated using second-order central differences. One-sided differences dependent on the sign of the control function Pn (backward for Pn < 0 and forward for Pn > 0) are useful to enhance convergence with very strong control functions. The control functions are evaluated either directly from the initial algebraic grid and then smoothed, or by interpolation from the boundary-point distributions. 126.96.36.199 Control Functions The now-standard procedure in block-structured systems is to first generate surface grids on block faces — both boundary and in-field block interfaces — from point distributions placed on the face edges by distribution functions. Then volume grids are generated within the blocks. In both this surface and volume grid generation, the first step is normally TFI, to be followed by elliptic generation with control functions interpolated into the field in accordance with boundary point distribution and surface curvature. The three components of the elliptic grid generation system, Eq. 1.1, provide a set of three equations that can be solved simultaneously at each point for the three control functions, Pn (n = 1,2,3), with the derivatives here represented by central differences. This produces control functions that will reproduce the algebraic grid from the elliptic system in a single iteration, of course. Thus, evaluation of the control functions in this manner would be of trivial interest except that these control functions can be smoothed before being used in the elliptic generation system. This smoothing is done by replacing the control function at each point with the average of the four neighbors in the two curvilinear directions (one in 2D) other than that of the function. Thus Pi is smoothed in the ξ j and ξ k directions, where i, j, k are cyclic. No smoothing is done in the direction of the function because to do so would smooth the spacing distribution. An algebraic grid is generated by transfinite interpolation (Chapter 3) from the boundary point distribution, to serve as the starting solution for the iterative solution of the elliptic system. With the boundary point distribution set from the hyperbolic sine or tangent functions, which have been shown to give reduced truncation error (Chapters 3 and 32), this algebraic grid has a good spacing distribution but may have slope breaks propagated from corners into the field. The use of smoothed control functions evaluated from the algebraic grid produces a smooth grid that retains essentially the spacing of the algebraic grid. The elliptic generation system can be solved by point SOR iteration using a field of locally optimum acceleration parameters. These optimum parameters make the solution robust and capable of convergence with strong control functions. ©1999 CRC Press LLC Control functions can also be evaluated on the boundaries using the specified boundary point distribution in the generation system, with certain necessary assumptions (orthogonality at the boundary) to eliminate some terms, and then can be interpolated from the boundaries in this manner. More general regions can, however, be treated by interpolating elements of the control functions separately. Thus control functions on a line on which ξ n varies can be expressed as Pn = An = + Sn ρn (1.2) where An is the logarithmic derivative of the arc length, Sn is the arc length spacing, and Ρn the radius of curvature of the surface on which ξ n is constant. The arc length spacing, Sn, and the arc length contribution, An, to the control function can be interpolated into the interior of the block from the four sides on which they are known by twodimensional interpolation. The radius of curvature, ρn, however is interpolated one-dimensionally between the two surfaces on which it is evaluated. The control function is finally formed by adding the arc length spacing divided by the radius of curvature to the arc length contribution according to Eq. 1.2. This procedure allows very general regions with widely varying boundary curvature to be treated. A more general construction of the control functions is given in Chapter 4. 126.96.36.199 Boundary Orthogonality The standard approach used to achieve orthogonality and specified off-boundary spacing on boundaries has been the iterative adjustment of control functions in elliptic generation systems, first introduced by Sorenson in the GRAPE code in the 1980s (Sorenson [1989]). Various modifications of this basic concept have been introduced in several codes, and the general approach is now common (see Chapter 6). A second-order elliptic generation system allows either the point locations on the boundary or the coordinate line slope at the boundary to be specified, but not both. It is possible, however, to iteratively adjust the control functions in the generation system until not only a specified line slope but also the spacing of the first coordinate surface off the boundary is achieved, with the point locations on the boundary specified. These relations can be applied on each successive coordinate surface off the boundary, with the off-surface spacing determined by a hyperbolic sine distribution from the spacing specified at the boundary. The control function increments are attenuated away from the boundary, and contributions are accumulated from all orthogonal boundary sections. Since the iterative adjustment of the control functions is a feedback loop, it is necessary to limit the acceleration parameters for stability. This allows the basic control function structure evaluated from the algebraic grid, or from the boundary-point distributions, to be retained, and thus relieves the iterative process from the need to establish this basic form of the control functions. The extent of the orthogonality into the field can also be controlled. This orthogonality feature is also applicable on specified grid surfaces within the field, allowing grid surfaces in the field to be kept fixed while retaining continuity of slope of the grid lines crossing the surface. This is quite useful in controlling the skewness of grid lines in some troublesome areas. Alternatively, boundary orthogonality can be achieved through Neumann boundary conditions, which allow the boundary points to move over a surface spline. The boundary point locations by Newton iteration on the spline to be at the foot of normals to the adjacent field points. This is done as follows: The Neumann point on the section currently closest to the field point R is first located. This is done by sweeping the section in ever expanding squares centered on the Neumann point. (These squares are actually limited by the section edges, of course, and hence, may become rectangles.) Next the quadrant about this closest point above which the field point lies is determined by comparing the dot products of the distance vector (from the closest point to the field point) with the tangent vectors (rξ , rη ) to the two grid lines on the section. The quadrant is identified by the signs of these two dot products. The Neumann boundary point in question, r, is then moved to the foot of the normal from the adjacent field point to the surface. This position is found as the solution of the nonlinear system ©1999 CRC Press LLC rξ ⋅ (R − r) = 0, rη ⋅ (R − r) = 0 (1.3) by Newton iteration. The location of the closest current boundary point and the examination of dot products described above has determined the surface cell, i.e., the quadrant, on which this solution lies so that the iteration can be confined to a single cell. Provision can also be made for extrapolated zero-curvature boundary conditions and for mirrorimage reflective boundary conditions on symmetry planes. 188.8.131.52 Surface Grids In the case of a curved surface, the surface is splined and the surface grid is generated in terms of surface parametric coordinates. The generation of a grid on a general surface (Chapter 9) is a two-dimensional grid problem in its own right, which can also be done either by interpolation or a PDE solution. In general, this is a 2D boundary value problem on a curved surface, i.e., the determination of the locations of points on the surface from specified distributions of points on the four edges of the surface. This is best approached through the use of surface parametric coordinates, whereby the surface is first defined by a 2D array of points, rmn' e.g., a set of cross sections. The surface is then splined, and the spline coordinates (u,v; surface parametric coordinates) are then made the dependent variables for the interpolation or PDE generation system. The generation of the surface grid can then be accomplished by first specifying the boundary points in the array rij on the four edges of the surface grid; converting these Cartesian coordinate values to spline coordinate values (uij, vij) on the edges; then determining the interior values in the arrays uij and vij from the edge values by interpolation or PDE solution; and finally converting these spline values to Cartesian coordinates rij (see Figure 9.1). 1.3.4 Hyperbolic System Elliptic generation systems operate throughout the entirety of a region, while hyperbolic systems move outward from boundaries. An alternate approach to boundary orthogonality and spacing is to incorporate a hyperbolic generation system near the boundary, transitioning to an elliptic system in the far field. It is also possible to base a grid generation system on hyperbolic partial differential equations, rather than elliptic equations (Chapter 5). In this case the grid is generated by numerically solving a hyperbolic system, marching in the direction of one curvilinear coordinate between two boundary curves in two dimensions, or between two boundary surfaces in three dimensions. The hyperbolic system, however, allows only one boundary to be specified, and is therefore of interest only for use in calculation on physically unbounded regions where the precise location of a computational outer boundary is not important. The hyperbolic grid generation system has the advantage of being generally faster than elliptic generation systems but, as just noted, is applicable only to certain configurations. Hyperbolic generation systems can be used to generate orthogonal grids. In two dimensions the condition of orthogonality is simply g12 = 0. If either the cell area g or the cell diagonal length (squared), g11 + g22, is a specified function of the curvilinear coordinates, i.e., g = F(ξ , η) or g11 + g22 = F(ξ , η) then the system consisting of g12 = 0 and either of the two equations just above is hyperbolic. Since the system is hyperbolic, a noniterative marching solution can be constructed proceeding in one coordinate direction away from a specified boundary. The cell volume distribution in the field can be controlled by the specified control function F. One form of this specification is as follows: Let points be distributed on a circle having a perimeter equal to that of the specified boundary at the same arc length distribution as on that boundary. Then specify a radial distribution of concentric circles about this circle according