Numerical Grid Generation
Foundations and Applications
By: Joe E. Thompson, Z.U.A. Warsi and C. Wayne Mastin
X. CONFORMAL MAPPING
Innovations in conformal mapping continue to extend this classical technique to more complicated configurations, and surveys of the various techniques available are given in Ref. [7] and Ref. [1]. Some specific recommendations of techniques and tools are given in Ref. [7]. Conformal systems have advantage the of introducing the fewest additional terms in transformed partial differential equations. Considerable understanding of the theory of functions of a complex variable may be necessary for effective applications, though.
Although the complex variable techniques by which conformal transformations are usually generated are inherently twodimensional, certain more general cases can be treated by rotating or stacking twodimensional systems:
Systems can also be generated on curved surfaces, as has been done by cartographers, for stacking. Examples of the use of conformal mapping in the construction of threedimensional configurations are noted in Ref. [5].
A curvilinear coordinate system generated by a conformal mapping is very rigid in the sense that little control can be exerted over the distribution of the grid points. Conformal mappings also do not exist in three dimensions (except for trivial cases). Furthermore, the coordinate system tends to be more difficult to construct than when using algebraic or elliptic systems. In spite of these facts, conformal mappings continue to play a significant role in grid generation. A number of recent developments and applications of conformal transformations are noted in Ref. [1], [5], and [7].
The desirability of a coordinate system generated by a conformal transformation lies in the form of the transformed equations. For example, consider the diffusion equation

(1) 
Now and satisfy the CauchyRiemann equations

(2) 
or equivalently

(3) 
It follows that in this case g_{12}=g_{21}=0 and g_{11}=g_{22}=. Equation (1) can be written in curvilinear coordinates, using Eq. (III46), as

(4) 
where the Laplacian is defined in terms of the curvilinear coordinates. Therefore it is observed that the diffusion equation remains essentially unchanged. The only effect of the transformation is a change in the diffusion coefficient. Neumann boundary conditions are also unchanged in conformal coordinates. The boundary condition
where is normal to a =constant coordinate line, is expressed in curvilinear coordinates as
1. Construction by FiniteDifferences
The literature abounds with methods for constructing conformal mappings. As can be seen in the review article, Ref. [1], these methods may include the construction of SchwarzChristoffel transformations, the solution of integral equations, or expansions in terms of power series or Fourier series. Since this chapter is not intended to be a comprehensive treatment of conformal mapping, only the simple, yet frequently used, finite difference method based on elliptic systems is discussed here.
Consider the problem of conformally mapping the interior of the contour onto the interior of a rectangle. The Riemann Mapping Theorem states that such a mapping exists, and it also implies that the mapping is uniquely determined by specifying three real parameters. Suppose we wish to indicate four specific points on which are to map to the vertices of the rectangle. If the rectangle is fixed, then the problem is overdetermined and no conformal mapping exists. Therefore, the mapping must determine one of the dimensions of the rectangular region which we will now denote as the set
Rather than allow a rectangle with variable width, one can equivalently introduce the parameter M in (3) so that

(5) 
where
The mapping is no longer conformal, but the conformal mapping can be easily obtained by simply multiplying the coordinate by M. On the unit square the functions x and y now satisfy

(6) 
Two boundary conditions are needed in order to determine a unique solution for this elliptic system. One condition is derived from the equation of the boundary curve I which might be

(7) 
The other condition comes from applying the orthgonality equation, g_{12}=0. This condition also follows on eliminating the parameter M in Eq. (5). The implementation of the boundary conditions is done in the following order. First a boundary value for x or y is computed from the orthogonality constraint. If the boundary point lies along =0, then we may use

(8) 
A forward difference is used to approximate the derivatives and central differences for the derivatives. The same equations are used along = 1 with backward differences for the derivatives. Once an x or y value has been computed from Eq. (8), the other coordinate value is given by writing (7) in the form

(9) 
Although either equation in (8) could be used, it is advisable to choose either the first or second equation, depending on whether x_{ } or y_{ } has the largest absolute value. This avoids not only the possibility of division by zero but also problems with the solvability of the implicit equation (7). The same techniques are used along an =constant coordinate line. In this case the orthogonality constraint can be written as

(10) 
Now the parameter M must also be determined. It follows from Eq. (5) that

(10) 
An iterative algorithm is used to construct the mapping, and any algorithm which can be used for the elliptic systems in Chapter VI can also be used here. At each iteration a new set of boundary values for x and y are computed using Eq. (8)(10). There are two options in computing a value for M. Either a different value at each point can be computed from Eq. (11), or a constant value can be computed from a relation such as

(12) 
where 01. Eq. (12) is derived from the equation in Eq. (5) by integrating along an constant coordinate line. This same technique can also be used to derive an alternate formula for finding the boundary values at x and y. Along = 0, for example, we have
The constant M, called the conformal module of the region by complex analysts, has a simple geometric interpretation. From Eq. (11) it is noted that M is simply the aspect ratio of the grid cells. There exist highly accurate numerical methods for computing both M and the boundary values for x and y. If these values are computed first, then the system (6) can be solved by a direct elliptic solver.
The only control over the distribution of grid points with a nonformal mapping is by changing the points which map to the vertices of the rectangular region. However, most of the advantageous features are retained when the conformal mapping is combined with onedimensional stretching transformations. Thus we will consider a new set of computational variables, and , with and serving as intermediate variables defined by the onedimensional equations

(13) 
If x and y are solutions of Eq. (6), then in terms of the new computational variables,

(14) 

(15) 
In Eq. (15), the covariant metric tensor components are defined relative to the transformation from the physical x,y variables to the computational , variables.
The application of this transformation to the diffusion equation (1) results in the following transformed equation:

(16) 
Note that the coefficients of the and derivatives in Eq. (16) are functions of and , respectively. Therefore, only onedimensional arrays are needed to store these coefficients. It can be further noted that the steadystate equation (A_{t}=0) is a separable elliptic equation which can be solved using a direct elliptic solver.
In the above development the stretching functions given by Eq. (13) are used to control the grid point distribution. Clearly the derivatives of these functions must be nonvanishing, and these derivatives may as well be taken to be positive so that the orientation of the physical boundary is preserved. The function f() is a contraction mapping if f'() < 1 and an expansion mapping if f'() > 1. Therefore, relative to a conformal mapping, the =constant coordinate lines will be closer together where f'() < 1 and farther apart when f'() > 1. The same control over the =constant coordinate lines is exerted by the function h(). Several onedimensional functions are discussed in Chapter VIII.
On any particular boundary segment, say =0, it is in theory possible to match any desired distribution of grid points provided the correct stretching function f() can be determined. However, there is no known way of generating this stretching function so that Eq. (14), together with the boundary conditions given by Eq. (7) and g_{12}=0, will have a solution with the prescribed boundary values along =0. The solution to this problem lies in the implicit determination of the stretching function in the solution of the elliptic system. Suppose that h() = , so that Eq. (15) can be written as

(17) 
This equation allows Eq. (14) to be written as

(18) 
This quasilinear system can be solved with the orthogonality condition on all boundary components except = 0 where we will now impose the Dirichlet condition
where defines the desired distribution of grid points.
The ability to specify grid points along a boundary component extends the usefulness of conformal mappings. For example, one can assign coordinates around an airfoil and along the branch cut in a Ctype coordinate system so that the coordinate lines pass smoothly through the cut. In many segmented systems the grid points can be chosen so that coordinate lines pass smoothly from one subregion into the next. One disadvantage of this method is the reported slow convergence in the iterative solution of (18) for certain problems. An alternate method of achieving the same result would be to generate a conformal mapping from Eq. (6) and then use interpolation to redistribute the grid lines. Note that the interpolation scheme may affect the orthogonality of the coordinate system to some degree.
2. SchwarzChristoffel Transformation
Conformal mappings of circular disks or halfplanes onto polygonal regions are defined by the SchwarzChristoffel formula. Suppose the points _{1},_{2},....,_{n}, lie on the real axis of the plane. Then the mapping defined by

(19) 
transforms the upper half plane onto a polygonal region with interior angles of  _{i} = _{i}. However this is not exactly what is needed in most grid generation problems. Presumably one would be given a polygonal region with vertices z_{1}, z_{2},...,z_{n}. Thus the parameters A,B,_{1},_{2},...,_{n} must be determined so that the real axis maps onto the given polygon:
There are several numerical techniques for the approximation of the parameters in the SchwarzChristoffel transformations. Since a conformal mapping of a simplyconnected region has three degrees of freedom, three of these parameters must be given in order for the mapping to be uniquely determined. In certain infinite regions, the value of B can be calculated from the asymptotic behavior of the mapping function. We can also set z_{n}=_{n}=0, which implies that A=0 from Eq. (19). The remaining parameters to be determined are _{1},_{2},...,_{n1}. Alternately, as is commonly done in bounded regions, we can choose the values _{1},_{2},_{3} which are to map the points z_{1},z_{2},z_{3}. In this case the parameters to be determined are A,B,_{4},_{5},...,_{n}. The basic algorithm for determining the unknown parameters consists of computing the distances , using Eq. (19) and a quadrature formula to approximate the integral, and then iterating on the parameters until these distances are correct. Once these parameters in the transformation have been computed to the desired accuracy, the image of any point in the upper halfplane is formed by numerically evaluating the integral in Eq. (19).
SohwarzChristoffel transformations are not limited to regions with polygonal boundaries. They can be used in composition with other conformal mapping methods to map regions with curved boundaries onto various computational regions. For example, an integral equation method can be used to map a physical region with curved boundary components onto the unit disk, which can be easily transformed onto the upper halfplane. Now the upper halfplane can be mapped onto the computation region, which may consist of several rectangular blocks, by Eq. (19). There are also direct generalizations of the SchwarzChristoffel transformation for regions with curved boundaries. These are obtained by considering the limiting case of Eq. (19) as n>.
Recent extensions of the SchwarzChristoffel transformation to curved contours have made this procedure a powerful tool for treating complicated internal and other configurations. These improvements also lead to smoother metric coefficients for boundaries with slope discontinuities than in older methods for the SchwarzChristoffel transformation. This procedure for the SchwarzChristoffel transformation may also be more efficient than other conformal procedures involving an intermediate mapping of a nearcircle for mapping contours and circles in some eases. Several sources on the recent developments and applications of the SchwarzChristoffel transformation are cited in Ref. [1] and [5].
3. Construction from Integral Equations
Integral equations have played a major role in the solution of partial differential equations. Mathematicians have often resorted to integral equations when attempting to prove the existence and uniqueness of solutions. Numerical analysts turned to the socalled panel methods for solving partial differential equations in two and threedimensional regions. These mehtods replaced the partial differential equations by a set of integral equations and thereby reduced the dimension of the problem, since panel methods only involve boundary integrals. The application of integral equations depends on the availability of fundamental solutions of the partial differential equation. Therefore they are especially useful in the solution of Laplace's equation. Numerous solutions of Laplace's equation can be generated by determining the real and imaginary parts of analytic functions. As most conformal mappings can be reduced to the solution of boundaryvalue problems for Laplace's equation, it should come as no surprise that integral equations can be a valuable tool in the construction of conformal mappings. Only the basic integral equation method of Symm (cf. Ref. [1]) will be presented here. This method has proven to be robust, yet is easily derived and involves only the solution of a system of linear equations.
Suppose the simplyconnected region D, bounded by the contour , is to be <1. Let z=z_{o} be the point in D which maps to the origin =0. If the Dirichlet problem

(20) 
can be solved and the harmonic conjugate h of q can be found, then it can be directly verified that the analytic function

(21) 
maps onto <1. Due to the form of the series expansion for the exponential, it can also be shown that this function has a nonvanishing derivative, and hence the conformal mapping of D onto the unit disk is given by Eq. (21). We now turn to the problem of solving the boundary value problem in Eq. (20). Suppose there exists a solution of the form

(22) 
for z on . Regardless of the value of the function (), the function q(z) is harmonic on D. In order that q(z) satisfy the boundary condition, it is clear that we need to choose () such that, for z on ,

(23) 
This is then the integral equation for determining the unknown function (). The harmonic conjugate of is arg(z). Thus the function of h(z) can be expressed as

(24) 
Note that the function h(z) is only unique up to an addition constant. The addition of a constant to h(z) results in a rotation of the conformal mapping defined in Eq. (21).
The practicality of this method depends on the efficient solution of the integral equation in Eq. (23). In order to solve this equation numerically, divide into n intervals, _{j}, j=1,2,...,n and assume () has a constant value _{j}, on _{j}. Let z_{j} be a fixed point of _{j}. Now Eq. (23) can be approximated by the linear system of equations

(25) 
There are two alternatives in computing the coefficients in this system. If the _{j} are assumed to be straight lines, then the integrals can be calculated analytically. Otherwise, each integral must be computed numerically. Once these coefficients have been computed, the system can be solved to yield a step function which approximates the function (). The values of _{j} are now used to estimate the functions q(z) and h(z):

(26) 
Again the above integrals would, in general, be computed numerically. These values of q(z) and h(z) would be substituted in Eq. (21) to yield the image in the unit disk of any given point z in the region D.
This integral equation method is a very efficient and accurate method. However, it has one deficiency in regard to grid generation and the numerical solution of partial differential equations. The transformation which is constructed maps the physical region D onto the canonical region, which in this case is the unit disk. The unit disk could be the computational region, or it could be mapped onto a rectangular region by an auxiliary transformation. In any case, what is needed is the mapping from the unit disk onto the physical region. Therefore an interpolation scheme would be needed to approximate the inverse of the computed mapping.
It is sometimes more efficient to generate the final grid by solving the Laplace system numerically with Dirichlet boundary conditions from the conformal transformations, especially if a fast Poisson solver can be applied.
4. Elementary Complex Transformations
An extensive list of complex mappings is compiled in Ref. [44]. However, these mappings are only for regions with special boundary curves. If a strictly conformal transformation is not necessary, then these mappings may be used to create what are called nearly conformal mappings. For example, suppose an airfoil shape can be modeled as the image of a circle under the Joukowski transformation

(27) 
Under the inverse transformation, a given airfoil will map to a curve which is nearly circular. The region about the nearly circular curve can be mapped onto the region about a circular region by a simple algebraic transformation. One scheme for accomplishing this final mapping would be to divide each complex number on a given ray from the center by the modulus of the complex number on the curve. The composite mapping in this case would be a nearlyconformal mapping of the exterior of the airfoil onto the exterior of a circle. The inverse mapping, which could be explicitly defined, would define a nearlyorthogonal 0type grid about the airfoil.
Analytic functions are not only of value in mapping regions about airfoils, but are also helpful in the more general problem of generating grids in the neighborhood of boundary points with slope discontinuities. With most algebraic methods of grid generation, these slope discontinuities will propagate into the physical region resulting in nonsmooth grid lines and the associated increase in truncation error in the numerical solution of partial differential equations. The general idea can be conveyed with the following example. Suppose we have a region where the boundary has an interior angle of at the point z_{o}. Under the mapping

(28) 
the corner is eliminated. While this simple mapping may be useful in transforming the interior of a contour, the mapping of the exterior region would not be onetoone. The elimination of corners for regions surrounding a contour can be effected by applying the KarmanTrefftz mapping defined by

(29) 
where is the conjugate of z_{o}. The exponent a depends on the exterior angle and the region should be translated, if necessary, so that is an interior point of the contour.
This transformation may be applied sucoessively to eliminate any number of corners on the boundary of the physical region.
Elementary complex functions can therefore serve to precondition a region. Corners which are to map to sides of a computational rectangle can be eliminated. Conversely, rightangle corners can be formed at points of the physical reigon which are to map to vertices of the computational region thereby eliminating problems of extreme nonorthogonality.
The trend in treating more complicated regions is to break the mapping up into a sequence of more simple mappings. Contours, such as airfoils, are generally mapped to nearcircles by one or more simple transformations, and then the nearcircle is mapped to a circle by a series transformation, e.g., the Theodorsen procedure. It is necessary for convergence that the near circle be sufficiently near to being a circle. A series for the differential form is generally superior to the usual Theodorsen form for general bodies. This series appears in terms of arc length and surface angle, rather than the polar coordinates of the Theodorsen form which can lead to infinite derivatives and multiple values. The ordering of the points can break down in the Theodorsen form for closely spaced points also. The differential form is applicable, however, as long as there are no corners, even for twisted contours. In this and other series transformations, the differential form is usually more tolerant of odd shapes.
Multiplebody configurations can be treated by a sequence of transformations which map each body to a circle in succession, while maintaining previously established circles. Another procedure, invovles iteratively mapping each body to a circle with no special consideration of the others. This process generally requires only a few iterations to converge. Some recent applications are noted in Ref. [1] and [5].