Tetrahedral Mesh Generator For CFD Simulation of ...

1 downloads 0 Views 943KB Size Report
conforming boundary recovery was also successfully done where .... First the swapping process to recover missing faces as much as possible and secondly by ...
Proceedings of ICEE 2009 3rd International Conference on Energy and Environment, 7-8 December 2009, Malacca, Malaysia

Tetrahedral Mesh Generator For CFD Simulation of Complex Geometry Ng Yee Luon

Mohd. Zamri Yusoff, Norshah Hafeez Shuaib

OYL Research & Development Sdn. Bhd., Sungai Buloh, Selangor, MALAYSIA [email protected]

Centre for Advanced Computational Engineering, College of Engineering, Universiti Tenaga Nasional, Kajang, MALAYSIA [email protected], [email protected]

used as a guide to decide where should the next vertex be inserted. In general, in order to ensure acceptable quality, the vertex should be inserted as far as possible from other vertices [6]. This will reduce the likelihood of generating bad elements with short edges.

Abstract— Computational Fluid Dynamics (CFD) simulation is becoming an important aspect in the design and analysis of engineering equipment. CFD allows engineers to reduce losses related to fluid flows and hence increase the performance. The first part of CFD simulation is building suitable mesh for the geometry to be analyzed. Tetrahedral mesh is widely used in meshing complex domain due to its flexibility. One of the common techniques used to generate tetrahedral mesh is Delaunay triangulation. Although Delaunay triangulation can always generate tetrahedral elements that fill any domain, the quality of the elements is not guaranteed. This is particularly important in CFD simulations whereby the accuracy of results and convergence rate will be highly affected by the mesh quality. This paper describes the development of a tetrahedral mesh generator which is capable of generating sufficiently high quality tetrahedral cells in arbitrary complex geometry. A few test cases are presented which show that slivers are eliminated and the overall quality is good. The combination of the constrained and conforming boundary recovery was also successfully done where all the missing faces are rebuilt.

In this paper, a complete three dimensional tetrahedral mesh generation algorithm based on Delaunay triangulation is presented. This includes generation of the tetrahedral elements, recovery of the boundary surfaces, removal of the elements outside the domain and methods to improve the tetrahedral elements quality. In addition to the above, a surface mesh generator which can be used to generate input for the tetrahedral volume mesh generator is also described. A few test cases are presented which show that slivers are eliminated and the overall quality is good. II.

Keywords-component; Mesh, Tetrahedral, CFD, Delaunay

I.

INTRODUCTION

Meshing is the process of discretizing problem domain into many sub domains before the numerical calculation can be performed. Among different types of mesh elements, tetrahedral is the most commonly used due to its flexibility and ability in filling any arbitrary three dimensional (3D) complex geometries. Many techniques are available to generate tetrahedral meshes such as Advancing front e.g. [1], Octree e.g. [2] and Delaunay e.g. [3]. Delaunay triangulation based algorithm is widely used due to its effectiveness. A set of triangles is a Delaunay triangle when none of the triangle having vertex inside their circum-circle. This characteristic avoids generation of long and skinny triangle making them good in numerical simulation.

A. BASE Searching Insertion of a new vertex, P, into a set of Delaunay triangles/tetrahedral will result some of the cells violating the Delaunay criteria. The affected tetrahedral elements need to be removed and replaced with new tetrahedral elements that include point P. To start the process, the location of P with respect to the initial tetrahedral set needed to be determined. The tetrahedron which contains the point P is called BASE [8].

Delaunay triangulation for two or three dimensional spaces can be achieved by using Bowyer-Watson algorithm which was introduced independently by Bowyer [4] and Watson [5]. In the Bowyer-Watson algorithm, each point is processed one at a time and any elements that violate the empty circle properties are deleted. Triangles/tetrahedral generated can be

978-1-4244-5145-6/09/$26.00 ©2009 IEEE

DELAUNAY TRIANGULATION

The Delaunay triangulation process with the BowyerWatson algorithm was started with five large tetrahedron that enclosed all the input surface vertices. The point is inserted successively one by one into the domain and the region affected by the new vertex is determined. This is followed by Delaunay triangulation of the affected domain which contains the newly introduced vertex. The process can be divided into three stages which are BASE searching, Determination of the CORE and Triangulation of the CORE. The following sections will briefly describe the three stages involved. More details description can be found in [7]

The search is started from the last constructed tetrahedron and travel from cell to cell until the cell that contains P is

330

found. The search moves from one cell to the next by simply comparing the volume vectors of the tetrahedron formed by face to be crossed, F, and point P and neighboring tetrahedron containing F. If the vectors are of the same sign, face F is crossed, bringing the search closer to P. Figure 1 shows the search path to the BASE in 2D.

In the present work, the missing segments and faces are recovered separately. All the input edges are checked against the generated edges and a new vertex is added at the center of the missing edge. This process is repeated until all missing segments are recovered. For the faces, all missing faces were first identified by comparing all input faces against the generated tetrahedralization data. The missing input faces are then flagged. The recovery process is divided into two procedures. First the swapping process to recover missing faces as much as possible and secondly by adding vertex at the missing face. Vertex is added at the circum-center when the circum-center lies inside the triangle or at the intersecting point between the missing triangle and the edge crossing it. B. Domain Recovery Since the Delaunay triangulation start by few large cells that contain entire input vertices, unwanted cells are also generated. They are removed before refinement process. Some of the tetrahedral which lie outside the domain can be easily identified since they use at least one of the non input vertices. The problem arises when the domain shape is having internal boundary. They cannot be identified through node index since all their nodes are input vertices. Following [9], the point in solid test is used to identify tetrahedral which are outside the domain.

Figure 1 : Search path to BASE. Only highlighted elements are checked. The rest are skipped

B. Determination of the CORE CORE is the cavity that is left behind after the tetrahedral that violate the Delaunay criteria are deleted. The CORE can be easily identified when the BASE is known. The search can start from the adjacent cells of the BASE, when they are found to be part of the CORE, their adjacent cells will be checked. The “branch” search process is continued until no new tetrahedral is part of the CORE. The boundary of the CORE is the faces remain after the tetrahedral is removed. The new tetrahedral set is constructed simply by connecting new point to the boundary face of the CORE.

IV.

A. Quality Measurement In numerical simulation, the shape of the elements affects the accuracy and stability of the simulation. For tetrahedral, large angle will give large interpolation error, which will lead to instability of the simulation; while small angle tends to produce ill condition matrix [10]. Hence the refinement scheme must be able to avoid the creation of these two types of elements. Quality measurement used must be able to detect such elements.

C. Triangulation of the CORE New tetrahedral are created by connecting vertex P to all the CORE boundary faces. This is followed by updating of the adjacency data. In order to ensure fast searching, the adjacency data for CORE boundary is calculated before the actual construction. The common edge of each face pair is also recorded. During actual construction, after a CORE face is used, its adjacent faces immediately appear in “next to use” list. The adjacency data can be updated without the need of search and match since the next three cells built must be the neighbors. If the neighbor’s faces are already used, only the adjacency data are updated. III.

MESH QUALITY IMPROVEMENT

In 2D, radius to shortest edge ratio is a good quality measurement, since minimizing this ratio will also maximize the minimum angle of the triangle. In 3D, the radius to shortest edge ratio in 3D can detect most of the bad quality tetrahedral, except for sliver. The creation of the sliver is the natural output of the Delaunay triangulation rather than numerical error [8]. Even well spaced input cannot prevent sliver [11, 12]. In the present work, radius to shortest edge ratio is used as the main measure of quality of the mesh. In order to detect sliver, dihedral angle of all elements are also measured.

BOUNDARY AND DOMAIN RECOVERY

After input vertices are tetrahedralized, some of the tetrahedral generated are outside the domain and they may not conform to the input boundary. The recovery of these boundary faces is necessary in order to maintain the validity of the domain.

B. Vertex Addition After the initial tetrahedralization and domain recovery, the tetrahedral are refined by inserting vertex at circum-centre of bad quality cells. The process is repeated until all tetrahedral are within the bound of the quality. If the maximum radius to shortest edge allowed is set to 1, the algorithm is guaranteed to terminate when a new vertex is inserted at the circum-centre of the elements as the edge of the new tetrahedron is at least of the same length as the radius of the deleted tetrahedron [6,13]. No

A. Boundary Recovery In 2D domain, the missing segment can easily be recovered through recursive splitting procedure. Since the Delaunay triangulation always connects a vertex to its nearest neighbor, the entire length of the missing segment can be recovered once the vertices spacing along the edge is sufficiently small [6]. But in 3D, a missing face may contain one or more missing edges.

331

edge shorter than pre-existing edge can be created. So, the algorithm will terminate when it run out of space to add new vertex.

C

C B

A D

“Auto sizing” is also included based on the input surface. Here, each vertex is assigned a weight based on the length of the edges using the vertex. Later, the maximum of the tetrahedral is the average weight of the tetrahedral vertices. Hence the tetrahedron using small input surface has small final size. A weight of the new vertex added is slowly increased by the expansion rate set by user.

D

C

A

C. Swapping After the vertex insertion process, the only bad tetrahedral lefts are those having high “sliverness”. They are improved by applying swapping. Swapping is the process of changing connectivity of the cell while keeping the vertex location. Since sliver cannot be detected through radius to shortest edge ratio, after vertex insertion process, the quality parameter used is dihedral angle. Both minimum and maximum dihedral angle is recorded to keep the current quality state of the tetrahedron. The mesh improvement scheme only target elements which have quality worst than a specified angle.

The more general swaps used for any number of tetrahedral larger than four are edge removal and multi-face removal. They are the reverse of each other. Edge removal replaces common edge of a group of tetrahedral with a few faces created by non common nodes of those elements. Multi-face removal, as the name suggests, remove all the faces between two common nodes and replace with an edge connecting them. Figure 4 shows these two operations, where common edge AB is replaced with 3 internal faces in edge removal process and in multi-face removal these three internal faces are deleted and edge AB is restored.

A

Multiface Removal

D

B

C

Swap A 3-2

A

E

Common edge

Swap 2-3

Swapping process is local and it only changes small number of elements around the targeted region. All dihedral angles after transformation are calculated and compared with the value before. The transformation is only accepted when the worst case after swapping is better than before. This is to ensure that the mesh can only be improved and never become worst.

D

D

C

C A A

B

Figure 4 : Multi-face removal and edge removal

B

B E

A Edge removal

E

C

B D

Figure 3 : Swap 2-2 and swap 4-4

Most of the sliver can be removed by applying 3-2 swap. The diagonal edge of the sliver is used as reference and the two other tetrahedral that are using this edge are found. Later, this common edge is swapped, which effectively removed the sliver. The swap 3-2 can also be reversed to become swap 2-3. In swap 2-3, common face of the two tetrahedral is destroyed and replaced with common edge that are used by three elements. Cheng et al. [11] removed most of the sliver through this method. Figure 2 shows the swap 3-2 and swap 2-3 operations.

D

C

A B

D

E

B

A

B

B

Figure 2 : Swap 3-2 and swap 2-3 Operation

Other swaps are swap 4-4 and swap 2-2. Swap 4-4 is used to change the common edge connectivity of the 4 tetrahedral without changing the number of elements. Swap 2-2 is very similar to swap 4-4 but working on the common edge that lies on the boundary and shared by two tetrahedral. Figure 3 shows these two operations.

D. Smoothing Smoothing is the opposite of swapping. It changes the position of the vertices without touching their connectivity. For triangular and tetrahedral elements, the popular smoothing method is Laplacian smoothing, which moves the vertex toward the centroid of the polygon formed by surrounding vertices [10]. With the connectivity table, the surrounding nodes can easily be found and the process is local as in swapping. The drawback of Laplacian smoothing is that its quality function of the tetrahedron is not a smooth function of

332

the vertex coordinate. Hence the quality after smoothing may be worst than before. Pierre et al. [12] used variational tetrahedral meshing to optimize vertex position by minimizing the energy function, but their mesh generated is not conform to input boundary. Klingner et al. [10] used non-smooth optimized smoothing to counter this problem and extends it to boundary smoothing for flat surface. In the present work, Laplacian smoothing is used due to its speed and ease of implementation. A similar approach as in swapping is used to ensure that elements quality can only improved after the smoothing process.

A. Stanford Bunny The model of Stanford Bunny, as shown in Figure 6, is repaired and re-meshed (closing a few holes around the original model). The input obtained from [14] consists of 7878 vertices and 15752 triangles. The growth rate, called GROW, is set to 1.5. After input vertices are tetrahedralized, some of the tetrahedral generated are outside the domain and they may not conform to the input boundary. Figure 7 shows the mesh generated.

The swapping-smoothing procedure is repeated iteratively. This means that after all tetrahedral are swapped, smoothing is applied to the remaining bad elements that failed to be improved by swapping. It was found that, this approach successfully improve many tetrahedral that cannot be improved by swapping or smoothing alone. After smoothing, even though the element quality is still far from ideal, it allowed the swapping operation to be applied on tetrahedron that failed to swap during the last iteration. V.

e2

e2

e1

e1

e2 e1

e2

e2 e1

e1

Edge e1-e2 Recovered

Figure 5 : Segment recovering by swap

SURFACE MESHING

A triangular surface mesher is incorporated into the current volume mesher. In order to allow 2D mesher to triangulate a surface in 3D space, the surface must first be converted to 2D plane. A transformation matrix is calculated and used to rotate the surface to XY plane. The transformation matrix is constructed using axis of rotation calculated from initial surface normal vector and the normal vector of the target plane. All vertices belonging to the surface are multiplied by the transformation matrix to eliminate the Z element. After surface triangulation process is completed, the inverse of the transformation matrix is used to rotate the surface back to its original plane. The 2D triangulation process is exactly the same as in 3D. In 2D, the radius to shortest edge ratio alone is sufficient as the quality measure for the triangles generated since all thin and skinny triangles have large radius to shortest edge ratio. Hence, minimizing the ratio will automatically ensure good angle bound [6].

Figure 6 : Stanford Bunny [14]

65377 elements were generated in 11 seconds. The worst element has minimum and maximum dihedral angles of 16.47ْ and 146.42ْ respectively with skewness of 0.77. Only 59 elements or 0.090% are having skewness larger than 0.6.

The complete constrained boundary recovery is possible in 2D problem. Only series of edge swapping is applied to recover the missing edge, without the need of any new vertex addition. When an input edge is detected as missing, all edges that are crossed by this missing edge are recorded and swapped one by one. If any swapping process resulted in invalid elements, the edge is initially skipped and only returned after all other missing edges are swapped. Figure 5 demonstrates an example of edge recovery by swapping in 2D problems. VI.

RESULTS AND DISCUSSION

Two models are tested in order to analyze the performance of the tetrahedral volume mesher. The first model is the famous Stanford Bunny which is widely used as computer graphic test case [14]. The second model is a typical public transport interchange (PTI) in Hong Kong [15]. The model is used to investigate the pollution distribution inside the PTI. Figure 7 : Internal mesh of Bunny with 1.5 GROW

333

The data in Table I show rapid increase in elements generated as the GROW was decreases. Only a slight increase of element number observed when GROW was reduced from 3.0 to 1.5. But rapid increased is observed when GROW is further reduced to 1.2. The closer GROW to 1, the more density of the volume mesh is affected by the density of the surface mesh. But, even if GROW is set to be very large, the radius to shortest edge ratio criteria will still reduce GROW to certain level. The initial mesh will contains large numbers of bad quality elements which will then improved in the mesh refinement stages through vertex insertion process until the quality is acceptable. With smaller GROW, the mesh is already of good quality even before the refinement process. The worst quality element for all three cases is located at the same region, as shown in Figure 11. The bad input surface triangle caused difficulty for the mesh refinement via swapping or smoothing process. This is because the angle of the tetrahedron is bound by the angle of the boundary triangle.

Different values of GROW are tested for the same model. Number of elements generated and their quality are compared. Figures 8 to 10 shows the cutting view of the model for GROW of 1.5, 1.2 and 3.0 respectively. Their relevant data is shows in Table I. TABLE I.

GROW Element num. Gen. time (s) Worst element (min/max ) % > 0.6 skew

DATA FOR DIFFERENT GROW

1.2 84181 10.0 13.96/145.32

1.5 65377 9.8 16.47/146.4

3.0 62805 9.7 14.37/143.7

0.051 %

0.090 %

0.10 %

Figure 8 : Cutting view for 1.5 GROW

Figure 11 : Worst element is affected by surface mesh

B. Public Transport Interchange This model is the typical public transport interchange (PTI) in Hong Kong [15]. The model consists of 8 buses, 7 columns and a wall. The node density along the edges is set to 0.5. The surface mesh are generated by the surface mesher. Figure 12 shows the domain in wire frame. Figure 13 and 14 shows the top and bottom view of the surface mesh generated.

Figure 9 : Cutting view for 1.2 GROW

The triangles are smaller near the edge due to the higher node density of the edge as input by user. Further away from the edge, the triangle starts to grow into a larger size. If the space is small, for example within the gap between bus and wall, the elements will remain small in order to effectively resolve the gap. With GROW set to 1.5, 27620 surface triangles and 135256 tetrahedral are generated. 0.102% of the elements have skewness larger than 0.6. The worst element has minimum angle of 24.26ْ and maximum angle of 146.24ْ. The total mesh generating time is 22.2s, where 2s is surface mesh generation. Figure 15 shows the cutting view of the element generated.

Figure 10 : Cutting view for 3.0 GROW

334

COLUMN

Slight adjustment of GROW can greatly increase or decrease the mesh size (and number), especially when the value is closer to 1.00. This is due to GROW controlling mesh density of both surface mesh and volume mesh. Comparing to fixed surface mesh input (which is not affected by GROW), GROW has more impact on the final mesh size for the model where surface mesh is generated by embedded surface mesher. Figure 16 shows the surface mesh of the same PTI model with GROW set to 1.1. Figure 17 shows the cutting view of the mesh generated. No large elements, as in Figure 15, is visible now.

BUS

WALL Figure 12 : PTI domain [15]

Figure 16 : PTI surface mesh with GROW 1.1 Figure 13 : PTI surface mesh top view

Figure 17 : Cutting view of the PTI with GROW 1.1

The number of surface triangles generated is 41872. Number of tetrahedral is 352889, which is more than double the element generated when GROW set to 1.5. 0.0246% elements have skewness larger than 0.6, with worst element having minimum and maximum angles of 19.24ْ and 145.90ْ respectively. The total mesh generation time is 38.8s with 2.437s is surface mesh generation time. Both models did not contain any sliver. But this is due to the input surface mesh which did not contain small angle and the surface is smooth and of good quality. Reports shows that most of the sliver survived is at the boundary, especially near small angle [6,10,11].

Figure 14 :PTI surface mesh bottom view

VII. CONCLUSION A complete 3D tetrahedral volume mesh generator based on the Delaunay triangulation is presented. It was found that the iterative swap-smooth approach, although simple, is very

Figure 15 : Cutting view of the PTI with GROW 1.5

335

effective in removing bad quality tetrahedral. The quality of the mesh generated is generally good.

[7]

ACKNOWLEDGMENT

[8]

The work described in this paper is funded by Ministry of Science, Technology and Innovation (MOSTI), Malaysia under ScienceFund Grant No. 01-02-03-SF0130. The financial support from MOSTI which enable the author to present the work is greatly acknowledged.

[9]

[10]

[11]

REFERENCES [1] [2]

[3]

[4] [5] [6]

R. Lohner, “Progress in Grid Generation via Advancing Front Technique,” Engineering with Computers, 12, 186-210, 1996. M.S. Shephard and M.K. Georges, “Three Dimensional Mesh Generation by Finite Octree Technique”, International Journal for Numerical Methods in Engineering, 32, 709-749, 1991 S. Rebay, “Efficient Unstructured Mesh Generation by Means of Delaunay Triangulation and Bowyer-Watson Algorithm”, Journal of Computational Physics, 106, 125-138, 1993. A. Bowyer, Computing Dirichlet tessellations, The Computer Journal, 24(2):162–166, 1981. D. F. Watson, Computing the n-dimensional tessellation with application to Voronoi polytopes, The Computer Journal, 24(2):167–172, 1981. J.R. Shewchuk, Tetrahedral Mesh Generation by Delaunay Refinement, Proceedings of the Fourteenth Annual Symposium on Computational Geometry, 86-95, 1998.

[12]

[13]

[14] [15]

336

Y.L. Ng, M.Z. Yusoff M. Z, N.H. Shuaib, Development of Improved Three Dimensional Unstructured Tetrahedral Mesh Generator, Proceedings of World Academy of Science, Engineering and Technology, Vol 59, November 2009. H. Borouchaki and S.H. Lo, Fast Delaunay Triangulation in Three Dimensions, Comput. Methods Appl. Mech. Engrg. 128: 153-167, 1995. J.R. Shewchuk, Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator. In First Workshop on Applied Computational Geometry. ACM, 124-133, 1996. B.M. Klingner, J.R. Shewchuk, Aggressive Tetrahedral Mesh Improvement, Proceedings of the Sixteenth International Meshing Roundtable. 3-23, 2007. S.W. Cheng, T.K. Dey, H. Edelsbrunner, M.A. Facello, S.H. Teng, Sliver Exudation, J.ACM, 47, 883-904, 2000. A. Pierre, S.D. Cohen, Y. Mariette, D. Mathieu, Variational Tetrahedral Meshing, Special issue on Proceedings of SIGGRAPH, ACM Transactions on Graphics 24: 617-625, 2005. J.R. Shewchuk, Theoretically Guaranteed Delaunay Mesh Generation – In practice, Short course, Fourteenth International Meshing Roundtable, San Diego, 2005. The Stanford 3D Scanning Laboratory at http://graphics.stanford.edu/data/3Dscanrep/ Z. Lin, F. Jiang, T.T. Chow, C.F. Tsang and W.Z. Lu, CFD Analysis of Ventilation Effectiveness in a Public Transport Interchange, Building and Environment, Vol. 41, pp.254-261, 2006.