Mesh generation

From CFD-Wiki

Introduction

The partial differential equations that governs fluid flow and heat transfer are not usually amenable to analytical solutions, except for very simple cases. Therefore, in order to analyze fluid flows, flow domains are split into smaller subdomains (made up of geometric primitives like hexahedra and tatrahedra in 3D, and quadrilaterals and triangles in 2D) and discretized governing equations are solved inside each of these portions of the domain. Typically, one of three methods is used to solve the approximate version of the system of equations: finite volumes, finite elements, or finite differences. Care must taken to ensure proper continuity of solution across the common interfaces between two subdomains, so that the approximate solutions inside various portions can be put together to give a complete picture of fluid flow in the entire domain. Each of these portions of the domain are known as elements or cells, and the collection of all elements is known as mesh or grid. The origin of the term mesh (or grid) goes back to early days of CFD when most analyses were 2D in nature. For 2D analysis, a domain split into elements resembles a wire mesh, hence the name.

An example of a 2D analysis domain (flow over a backward facing step) and its mesh are shown in pictures below.

The process of obtaining an appropriate mesh (or grid) is termed mesh generation (or grid generation), and has long been considered a bottleneck in the analysis process due to the lack of a fully automatic mesh generation procedure. Specialized software progams have been developed for the purpose of mesh and grid generation, and access to a good software package and expertise in using this software are vital to the success of a modeling effort.

Mesh classification

As CFD has developed, better algorithms and more computational power has become available to CFD analysts, resulting in diverse solver techniques. One of the direct results of this development has been the expansion of available mesh elements and mesh connectivity (how cells are connected to one another). The elements in a mesh can be classified in various ways - the easiest is based upon the dimension and type of the elements. Common elements in 2D are triangles or rectangles, and common elements in 3D are tetrahedra or bricks. The most basic form of mesh classification is based upon the connectivity of the mesh: structured or unstructured.

Structured Meshes

A structured mesh is characterized by regular connectivity that can be expressed as a two or three dimensional array. This restricts the element choices to quadrilaterals in 2D or hexahedra in 3D. The above example mesh is a structured mesh, as we could store the mesh connectivity in a 40 by 12 array. The regularity of the connectivity allows us to conserve space since neighborhood relationships are defined by the storage arrangement. Additional classification can be made upon whether the mesh is conformal or not.

Unstructured Meshes

An unstructured mesh is characterized by irregular connectivity is not readily expressed as a two or three dimensional array in computer memory. This allows for any possible element that a solver might be able to use. Compared to structured meshes, the storage requirements for an unstructured mesh can be substantially larger since the neighborhood connectivity must be explicitly stored.

Hybrid Meshes

A hybrid mesh is a mesh that contains structured portions and unstructured portions. Note that this definition requires knowledge of how the mesh is stored (and used). There is disagreement as to the correct application of the terms "hybrid" and "mixed." The term "mixed" is usually applied to meshes that contain elements associated with structured meshes and elements associated with unstructured meshes (presumably stored in an unstructured fashion).

Algorithms

Since there are many different flow solvers with many different mesh requirements, there are also many different mesh generation algorithms. These are best classified in terms of the mesh that is generated: structured, unstructured, or hybrid.

Structured Algorithms

Many of the algorithms for the generation of structured meshes are descendents of "numerical grid generation" algoritms, in which a differential equation is solved to determine the nodal placement of the grid.

Algebraic Grid Generation

The simplest way to obtain a grid would be to specify the grid coordinates $\vec{x}$ as the result of some vector function, or

$\vec{x}=\vec{x}\left(\vec{\xi}\right),$

where $\vec{\xi}$ is the "index" vector, sometimes referred to as a computational coordinate. For our purposes here the entries of the computational coordinate will range from zero to a maximum. If such a function can be found for a given geometry, then the actual generation of gridpoints is straightforward. The problem, however, is that the determination of the function is not neccessarily that easy. In practice, it is sometimes easier to add an itermediate parametric space, denoted by $\vec{s}$ , in between the physical space representation of the grid and the computational space representation of the grid:

$\vec{x}=\vec{x}\left(\vec{s}\left(\vec{\xi}\right)\right).$

The entries in the computational coordinate are taken from the unit interval. This representation can help simplify matters, especially in the one dimensional case.

Many mesh generation systems (both structured and unstructured) require the generation of boundary grids before interior cells can be generated. This is an area in which algebraic grid generation is ideal - typically, we want to specify boundary edge point distributions quickly, with a minimum of complexity, and a high degree of repeatability. These functions are often referred to as stretching functions, with hyperbolic trigonometric functions such as the hyperbolic tangent as a popular choice. A simple one-parameter hyperbolic tangent stretching function is defined by

$s\left(\xi\right) = 1 + \frac{\tanh\left[\delta\left(\xi/I-1\right)\right]}{\tanh\left(\delta\right)},$

where $\delta$ is the stretching factor and $\xi\in[0,I]$ . This function partitions the unit interval and allows the specification of a single location. This sort of distribution is good for wall-normal grid distribution in viscous flows. This distribution is due to Vinokur [REF?]. Vinokur's procedure for the determination of the proper stretching factor to obtain desired spacings uses the derivatives of the stretching functions. Suppose we wish for our first grid spacing to be $\Delta s$ . This can be taken to mean that $s(1)=\Delta s$ or that $ds/d\xi (1) = \Delta s$ . Vinokur's procedure guarantees the latter.

A related double-sided stretching function (that gives symmetric spacings about $\xi=I/2$ ) is given by

$u\left(\xi\right) = \frac{1}{2}\left[1 + \frac{\tanh\left[\delta\left(\xi/I-1/2\right)\right]}{\tanh\left(\delta/2\right)}\right].$

This function is good for duct flows, such as turbulent channel flow. In situations in which different grid spacings are desired, a stretching function can be constructed that has specified spacings at both ends: $\Delta s_1$ and $\Delta s_2$ . Vinokur gives such a function, first defining

$A=\frac{\sqrt{\Delta s_2}}{\sqrt{\Delta s_1}},\ \ B=\frac{1}{I\sqrt{\Delta s_2\Delta s_1}}.$

The stretching factor $\delta$ is found from the solution of the transcendental equation

$\frac{\sinh\left(\delta\right)}{\delta} = B.$

The final grid distribution is then given by

$s\left(\xi\right)=\frac{u\left(\xi\right)}{A + (1-A)u\left(\xi\right)}.$

Again, Vinokur's procedure ensures that the derivative conditions $ds/d\xi (1) = \Delta s_1$ and $ds/d\xi (I) = \Delta s_2$ , and not the grid spacings obtained in the direct evaluation of the stretching function. To compute the actual gridpoints in practice, we need to express the grid edge in terms of a parameter. For example, a line segment between two points $\vec{x}_1$ and $\vec{x}_2$ can be expressed as

$\vec{x}(s)=\vec{x}_1 + \left(\vec{x}_2-\vec{x}_1\right)s.$

We can then use any of the above streching functions to generate a properly spaced grid on the given line segment. Other parametric forms are available as well. Of particular interest is the cubic Bézier curve, which allows the specification of direction and location at both endpoints and can be written as

$\vec{x}(s)=\vec{P}_0 (1-s)^3 + 3\vec{P}_1 s (1-s)^2 + 3 \vec{P}_2 s^2 (1-s) + \vec{P}_3 s^3,$

where the $P_i$ 's are the control points.

For the generation of interior cells, algebraic techniques are also available, most usually in the form of interpolation between boundary faces.

Elliptic Grid Generation

The oldest numerical grid generation techniques are based upon the solution of elliptic PDE's. Typically, a Poission-type equation is solved given the boundary grid distribution to generate interior nodal points. The solution domain is often topologically equivalent to a cube in 3D and a square in 2D. Consider the solution domain shown below with the indicated boundary resolution.

The simplest technique we could use here would be a solution of the Laplace equation using the standard second order finite difference stencil. This approach simplies into a form that is easily solved using Jacobi or Gauss-Seidel iterative techniques:

$\nabla\vec{x}=0\$

with Dirichlet boundary conditions is discretized as

$\ \vec{x}_{i,j} = \frac{\vec{x}_{i+1,j}+\vec{x}_{i-1,j}+\vec{x}_{i,j+1}+\vec{x}_{i,j-1}}{4}$

An initial grid is shown below in (a), and the resulting grid (after iterations) is shown in (b).

Note that the grid spacing near the curved section increases and then decreases as we move left to right, and that grid lines near the left and right boundaries are not very orthogonal. These issues are reasons that production grid generation techniques are usually more complicated. The addition of control functions allows for better grid clustering properties, which will be necessary for viscous flow simulations.

Grid Marching Methods

Unstructured Algorithms

It is difficult make general statements about unstructured mesh generation algorithms because the most prominent methods are very different in nature. The most popular family of algorithms are those based upon Delaunay triangulation, but other methods, such as quadtree/octree approaches are also used.

Delaunay Methods

Many of the commonly used unstructured mesh generation techniques are based upon the properties of the Delaunay triangulation and its dual, the Voronoi diagram. Given a set of points in a plane, a Delaunay triangulation of these points is the set of triangles such that no point is inside the circumcircle of a triangle. The triangulation is unique if no three points are on the same line and no four points are on the same circle. A similar definition holds for higher dimension, with tetrahedra replacing triangles in 3D.

Quadtree/Octree Methods

This article is a stub, a short article which needs to be improved. You can help by expanding it.