网格生成mesh genaration(英文)

e cc p Mesh Generation Mark Filipiak Edinburgh Parallel Computing Centre The University of Edinburgh Version 1.0 November 1996 Edinburgh Parallel Computing Centre 3 Table of Contents 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1 Discretisation and mesh type . . . . . . . . . . . . . . . . . . . 6 1.2 Mesh characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 Structured meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1 Boundary-fitted meshes . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Problem solution on curved grids . . . . . . . . . . . . . . 10 2.3 Boundary fitting grids on a single block . . . . . . . . . 11 2.4 Algebraic grid generation: interpolation . . . . . . . . . 11 2.4.1 Transfinite interpolation (TFI). . . . . . . . . . . . 11 2.5 PDE grid generation . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.5.1 Elliptic grid generation. . . . . . . . . . . . . . . . . . 14 2.6 Implementation in 3D . . . . . . . . . . . . . . . . . . . . . . . . 16 2.7 Other methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.7.1 Hyperbolic grid generation . . . . . . . . . . . . . . 17 2.7.2 Parabolic grid generation. . . . . . . . . . . . . . . . 17 2.8 Multiblock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.8.1 C, O, H grids . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.8.2 Multiblock . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3 Unstructured Meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1 Mesh requirements for the Finite Element Method23 3.2 Mesh generation methods . . . . . . . . . . . . . . . . . . . . . 24 3.2.1 Decomposition and mapping . . . . . . . . . . . . 24 3.2.2 Grid based methods . . . . . . . . . . . . . . . . . . . . 25 3.2.3 Advancing front . . . . . . . . . . . . . . . . . . . . . . . 26 3.2.4 Delaunay triangulation . . . . . . . . . . . . . . . . . 31 3.2.5 Other methods. . . . . . . . . . . . . . . . . . . . . . . . . 35 3.2.6 Smoothing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4 Adaptive Meshing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.1 Adaptive meshes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.2 Parallel mesh generation . . . . . . . . . . . . . . . . . . . . . . 38 5 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 Edinburgh Parallel Computing Centre 5 1 Introduction Continuous physical systems, such as the airflow around an aircraft, the stress con- centration in a dam, the electric field in an integrated circuit, or the concentration of reactants in a chemical reactor, are generally modelled using partial differential equa- tions. To perform simulations of these systems on a computer, these continuum equa- tions need to be discretised, resulting in a finite number of points in space (and time) at which variables such as velocity, density, electric field are calculated. The usual methods of discretisation, finite differences, finite volumes and finite elements, use neighbouring points to calculate derivatives, and so there is the concept of a mesh or grid on which the computation is performed. There are two mesh types, characterised by the connectivity of the points. Structured meshes have a regular connectivity, which means that each point has the same number of neighbours (for some grids a small number of points will have a different number of neighbours). Unstructured meshes have irregular connectivity: each point can have a different number of neighbours. Figure 1 gives an example of each type of grid. In some cases part of the grid is structured and part unstructured (e.g., in vis- cous flows where the boundary layer could be structured and the rest of the flow unstructured). Figure 1: Structured mesh (left) and unstructured mesh (right). In the rest of this chapter, the various discretisation methods are described, with their mesh requirements. Chapter 2 describes the methods used to generate structured meshes in simple domains (algebraic and elliptic methods) and the extension to com- plex domains using multiblock. Chapter 3 describes methods used to generate unstructured meshes, concentrating on the two main methods for producing triangu- lar/tetrahedral meshes: advancing front and Delaunay triangulation. Chapter 4 gives a very brief introduction to adaptive meshing. 6 Technology Watch Report 1.1 Discretisation and mesh type The main discretisation methods are finite differences [8], finite volumes (which is equivalent to finite differences) [9] and finite elements [10]. To illustrate the methods, we consider the conservation form of the convection equation where is the density, is the velocity, and is a source term. The term is the flux of . The finite difference formulation approximates the derivatives using neighbouring points, e.g., for a regular, rectangular grid with spacing in the direction . Although irregular grids can be used for finite differences, regular grids are invariably used. These lead to simple difference schemes and efficient solvers (on vector machines, for example). There can be problems with the finite difference scheme at coordinate singularities for certain regular grids (e.g., spherical polar coordinates) [12] In the finite volume formulation, the physical space is split up into small volumes and the partial differential equations integrated over each of these volumes: Then the variables are approximated by their average values in each volume, and the fluxes through the surfaces of each volume are approximated as function of the varia- bles in neighbouring volumes. Finite volume discretisation can use both regular and irregular meshes. In an irregu- lar mesh the fluxes through the surfaces are still well defined. The finite element method also splits the spaces up into small volumes, the ‘elements’. In each element, the variables and fluxes are approximated using weighting func- tions. The computational variables are the coefficients of these weighting (or ‘shape’) functions. Finite element methods generally use irregular meshes: there is no special advantage in using regular meshes. 1.2 Mesh characteristics For all types of meshes, there are certain characteristics that we want to control. • The local density of points. High density gives more accuracy, but computation takes longer. This leads to adaptive meshing methods , treated very briefly in t∂ ∂ρ ρU( )∇•+ S= ρ U S ρU ρ h x x∂ ∂ρ 1 h--- ρ xn 1+( ) ρ xn( )–[ ]≈ V td d ρ Ωd V ∫ ρU( ) n• dΓ ∂V ∫°+ S Ωd V ∫= Edinburgh Parallel Computing Centre 7 Chapter 4. • The smoothness of the point distribution. Large variations in grid density or shape can cause numerical diffusion or anti-diffusion and dispersion or refrac- tion of waves. This can lead to inaccurate results or instability [12]. • The shape of the grid volumes. For instance, boundary layers in fluid flow re- quire a grid that is very compressed normal to the flow direction. In the finite element method using triangular elements the maximum angle must be bound- ed strictly below to prove convergence of the method as the element size is reduced [11]. For simple domains, the choice between regular or irregular meshes is governed mainly by the discretisation method. However, for complex domains (e.g., the inside of a hydro-electric turbine) irregular meshes are preferred to regular meshes. Irregu- lar mesh generation (at least for triangular or tetrahedral elements) can be fully auto- matic and fast. Regular mesh generation requires the domain to be split up into simple blocks which are then meshed automatically. This block decomposition is at best semi-automatic and can require man-months of user effort. pi Edinburgh Parallel Computing Centre 9 2 Structured meshes This chapter begins with a discussion of boundary-fitted grids and the discretisation of PDEs on them, then deals with the main methods for grid generation of simple domains (using algebraic or differential equation methods) and then explains the multiblock concept used for more complicated domains. References [12] and [13] are detailed expositions of structured mesh generation. 2.1 Boundary-fitted meshes Structured meshes are characterised by regular connectivity, i.e., the points of the grid can be indexed (by 2 indices in 2D, 3 indices in 3D) and the neighbours of each point can calculated rather than looked up (e.g., the neighbours of the point are at , , etc.). Meshes on a rectangular domain are trivial to generate (though come care needs to be taken in the discretisation at convex corners) and struc- tured mesh generation techniques concentrate on meshing domains with irregular boundaries, e.g., the flow in blood vessels, or the deformation, flow and heat transfer in metal being formed in dies. Generally the meshes are generated so that they fit the boundaries, with one coordinate surface forming (part of) the boundary. This gives accurate solutions near the boundary and enables the use of fast and accurate solvers. For fluid flow these grids allow the easy application of turbulence models, which usu- ally require the grid to be aligned with the boundary. The alternative is to use a rec- tangular grid which is clipped at the boundary, with local grid refinement near sharp features on the boundary (‘Cartesian grids’). This will reduce the truncation order at the boundary and will require the mesh cells to be clipped at the boundary, increasing the complexity of the solver (even for a finite volume code). Cartesian grid generation is very fast and has been used for Euler aerodynamics [14]; it does not appear to be applicable to Navier-Stokes aerodynamics. The most common method of generating boundary-fitting grids is to have one contin- uous grid that fits to all the boundaries. The effect is to fit a contiguous set of rectan- gular computational domains to a physical domain with curved boundaries (Figure 2). Figure 2: Two-block boundary-fitted grid It is difficult to fit complex domains with one mapping from a rectangular computa- tional domain without generating excessively skewed grids. To get round this prob- i j,( ) i 1 j,+( ) i 1 j,–( ) 10 Technology Watch Report lem the domain is split up into blocks and each block is gridded, with some continuity requirements at the block interfaces; this is multiblock. The decomposition of the domain into blocks is usually done manually using CAD techniques and is slow. An alternative to continuous boundary-fitted grids with multiple blocks is to use a boundary fitting grid near each boundary, and simple rectangular grid in the interior, and interpolate between them. These are called overset or chimera grids (Figure 3)[15]. Figure 3: Overset grids This type of grid is easier to generate than a multiblock grid since each grid is local and doesn’t need to match to others. The individual grids will generally be of high quality (low distortion). However, the interpolation can be difficult, especially with more than two grids overlapping, and it adds to the solver time (increasing it by 10%- 30%). The overlapping grids cannot be too different in resolution and this can cause problems with the grids required for boundary layers in viscous flow simulations. Chimera grids are very useful for moving boundaries (e.g., helicopter blades) or mul- tiple boundaries (e.g., particles in a fluid). Most of the grid remains fixed but the inter- polation changes as the grids move with the boundaries. Chimera grids do have certain advantages, and the recent work [16] on conservative interpolation methods have increased their usefulness. However, the bulk of struc- tured grid generation is based on multiblock type grids, and we will concentrate on these for the rest of this chapter. 2.2 Problem solution on curved grids Once the grid is generated, the physical problem has to be discretised and solved on this grid. It is possible to calculate finite difference or finite volume formulas using the physical grid directly, but this will reduce the truncation order of the FD/FV scheme since the grids are generally non-uniform. The preferred method is to transform the equations used to model the problem into computational space. Since the physical grid is defined as some transformation from the rectangular computational grid, this process is straightforward. The resulting equations in the computational space will include the effects of the coordinate transformation, e.g., time derivatives of moving grids. The finite differencing can be done on the transformed equations, on the com- putational grid. Since the computational grid is usually completely regular, high order methods can be used, and the truncation order is preserved if the resulting FD/ FV equations are transformed back to the physical grid. (Flow solvers can be devel- oped to use chimera grids, with the extra book-keeping required to interpolate between each grid.) Although the order of FD/FV methods can be preserved, there are other effects of cur- vilinear grids on the accuracy of the solution. Edinburgh Parallel Computing Centre 11 • Grid size variation gives numerical diffusion (or anti-diffusion leading to insta- bility). The effect is worst with small grid size and large solution gradients. • Non-orthogonality, or skewness, of the (physical) grid increases the coefficient of the truncation error, but doesn’t change the order. Skewness of up to are acceptable in the centre of the grid, but one-sided differences are used at bound- aries, where orthogonality should be maintained. These limits on skewness are the main reason that multiblock methods are used. • In a multiblock scheme, the corner junctions between blocks are usually points with non-standard connectivity (e.g., in 2D, 5 grid lines from a point rather than 4) and these need special treatment (e.g., for finite volume discretisation these points need to be discretised in physical space). 2.3 Boundary fitting grids on a single block For simple domains, a single grid can be used without leading to excessive skewness. The computational space is a rectangle or cuboid, or at least has a rectangular bound- ary, with a regular, rectangular grid. We then need to define a 1-1 mapping from the computational space to the physical space. The methods currently used are • algebraic grid generation, where the grid is interpolated from the physical boundaries, • methods where a PDE for the grid point coordinates is solved, • variational grid generation, where some functional, e.g., smoothness, is maxim- ised/minimised (variational grid generation is not treated in this report, see [13] for details), • other methods, (not treated in this report), e.g., conformal mapping [12], [13]. 2.4 Algebraic grid generation: interpolation The method used in the algebraic method of grid generation is transfinite interpolation (TFI). 2.4.1 Transfinite interpolation (TFI) Working in 2D (the extension to 3D is straightforward), take the computational domain as the unit square with coordinates and , and the physical domain as some region with coordinates and . To generate the grid in the physical space, we would create a grid in the unit square and map this into the physical domain. There are two requirements for this mapping (Figure 4): • It must be 1-1, and • the boundaries of the computational space must map into the boundaries of the physical space, 45° 0 1,[ ] 0 1,[ ]× s t x y 12 Technology Watch Report but otherwise the mapping can be arbitrary, although it may not produce a good grid. Figure 4: Mapping of boundaries TFI is one kind of mapping, where the physical coordinates, treated as a function of the computational coordinates, are interpolated from their values on the boundaries of the computational domain. Since the data is given at a non-denumerable number of points, the interpolation is called transfinite. The 2D interpolation is constructed as a linear combination of 2 1D interpolations (also called projections) and their product. First define the blending functions , and , that will be used to interpolate in each direction. interpolates values from the boundary into the domain, interpolates values from the boundary into the domain, and similarly for , in the direction. The require- ments on and are and similarly for , in the direction. The simplest blending function is linear, giving linear interpolation The 1D transfinite interpolations (projections) are formed as follows (for the coordi- nate) s t 10 0 1 Computational domain Physical domain φ0 φ1 θ0 θ1 φ0 s 0= φ1 s 1= θ0 θ1 t φ0 s( ) φ1 s( ) φ0 0( ) 1 φ0 1( ), 0= = φ1 0( ) 0 φ1 1( ), 1= = θ0 θ1 t φ0 s( ) 1 s–= φ1 s( ) s= x P s x[ ] s t,( ) φ0 s( ) x 0 t,( ) φ1 s( ) x 1 t,( )+= Pt x[ ] s t,( ) θ0 t( ) x s 0,( ) θ1 t( ) x s 1,( )+= Edinburgh Parallel Computing Centre 13 and the product projection is which is the finite interpolant for the values of at the four corners of the computa- tional domain. The 2D transfinite interpolation is then which interpolates to the entire boundary. To form the grid in the physical space, this interpolant is used to map the points of a regular grid in the computational space (Fig- ure 5). Figure 5: Transfinite interpolation It is possible to extend the TFI to so that it interpolates to several coordinate lines, not just the boundaries. This is useful to control the grid density and also may be required to ensure a 1-1 mapping (i.e, a valid grid). The blending functions would then be, for example, cubic splines. It is also possible to match normal derivatives, e.g. , at the boundaries. This allows grid slope continuity between the blocks of a multiblock grid. There are some problems with TFI. The mapping will propagate boundary singulari- ties (corners) into the interior of the domain, which is not ideal for fluid flow simula- P s Pt x[ ] s t,( ) φ0 s( ) θ0 t( ) x 0 0,( ) φ1 s( ) θ0 t( ) x 1 0,( )+= φ+ 0 s( ) θ1 t( ) x 0 1,( ) φ1 s( ) θ1 t( ) x 1 1,( )+ x P s Pt⊕( ) x[ ] Ps x[ ] Pt x[ ] PsPt x[ ]–+= 10 0 1 s t t s dx ds⁄ 14 Technology Watch Report tions. A more serious problem is that if the mapping is not 1-1 then overspill occurs (Figure 6). Figure 6: Overspill This can be corrected be re-parameterising the boundaries, or adding constraint lines (surfaces in 3D) inside the domain. However, TFI grid generation is very fast and is the only competitive method for 3D. It is possible to use an PDE based generator to refine a grid produced using TFI. 2.5 PDE grid generation The idea behind these methods is to define PDEs for the physical space coordinates in terms of the computational space coordinates, and then solve these equations on a grid in computational space to create a grid in physical space. The main type of PDE used is elliptic, which is suitable for domains with closed boundaries (for unbounded domains, a fictitious boundary at large distances is used). The alternatives are hyper- bolic and parabolic equations, which are used for unbounded domains. 2.5.1 Elliptic grid generation These methods start with an PDE for the computational space coordinates and in terms of the physical space coordinates and (the extension to 3D is straight- forward). The simplest example is the Laplace equation 10 0 1 t s s t ξ1 ξ2 x 1 x 2 ∇2ξi 0= Edinburgh Parallel Computing Centre 15 Parts of the physical boundary, corresponding to coordinate lines, are then ‘equipo- tential surfaces’ (Figure 7). Figure 7: Elliptic grid generation boundaries In 2D, the extremum principle guarantees that the mapping is 1-1, and the Laplace operator is smoothing so that boundary slope discontinuities are