UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The impact of mesh regularity on errors Fan, Hongliang 2017

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2017_may_fan_hongliang.pdf [ 5.36MB ]
JSON: 24-1.0343046.json
JSON-LD: 24-1.0343046-ld.json
RDF/XML (Pretty): 24-1.0343046-rdf.xml
RDF/JSON: 24-1.0343046-rdf.json
Turtle: 24-1.0343046-turtle.txt
N-Triples: 24-1.0343046-rdf-ntriples.txt
Original Record: 24-1.0343046-source.json
Full Text

Full Text

The Impact of Mesh Regularity onErrorsbyHongliang FanB.Eng., Beihang University, 2014A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)March 2017c© Hongliang Fan 2017AbstractMesh quality plays an important role in improving the accuracy of the numerical simulation.There are different quality metrics for specific numerical cases. A regular mesh consisting ofthe equilateral triangles is one of them and is expected to improve the error performance. In thisstudy, Engwirda’s frontal-Delaunay scheme and Marcum’s advancing front local reconnectionscheme are described along with the conventional Delaunay triangulation. They are shownto improve the mesh regularity effectively. Even though several numerical test cases showthat more regular meshes barely improve the error performance, the time cost in the solverof regular meshes is smaller than the Delaunay mesh. The time cost decrease in the solverpays off the additional cost in the mesh generation stage. For simple test cases, more regularmeshes obtain lower errors than conventional Delaunay meshes with similar time costs. Formore complicated cases, the improvement in errors is small but regular meshes can save time,especially for a high order solver. Generally speaking, a regular mesh does not improve theerror performance as much as we expect, but it is worth generating.iiPrefaceThe research ideas and methods explored in this thesis are the fruits of a close working re-lationship between Dr. Carl Ollivier-Gooch and Hongliang Fan. The implementation of themethods, the data analysis, and the manuscript preparation were done by Hongliang Fan withinvaluable guidance from Dr. Carl Ollivier-Gooch throughout the process.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.1 Delaunay mesh refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.3 Finite-volume solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253 Frontal-Delaunay triangulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 30ivTable of Contents4 Advancing-front local reconnection method . . . . . . . . . . . . . . . . . . . . 405 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 495.1 Poisson equation on a square domain . . . . . . . . . . . . . . . . . . . . . . 495.2 Advection-diffusion equation in a channel . . . . . . . . . . . . . . . . . . . . 595.3 Potential flow around a circular cylinder . . . . . . . . . . . . . . . . . . . . . 665.4 NACA-0012 airfoil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 736 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84vList of Tables5.1 Mesh Generation Time Cost . . . . . . . . . . . . . . . . . . . . . . . . . . . 575.2 Time cost in mesh generation and second order solver . . . . . . . . . . . . . . 665.3 Time cost in mesh generation and fourth order solver . . . . . . . . . . . . . . 665.4 Time cost in mesh generation (seconds) . . . . . . . . . . . . . . . . . . . . . 735.5 Fourth order solver time cost . . . . . . . . . . . . . . . . . . . . . . . . . . . 735.6 Comparison of time cost in mesh generation and the fourth order solver (α = 2◦) 775.7 Time cost in the fourth order solver . . . . . . . . . . . . . . . . . . . . . . . . 78viList of Figures1.1 Meshes around a circular cylinder (Near View) . . . . . . . . . . . . . . . . . 21.2 Simple 1-D example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.1 Delaunay Triangulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.2 Delaunay Triangulation breaks the boundary . . . . . . . . . . . . . . . . . . . 92.3 Constrained Delaunay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.4 Diametral circle (dotted line) and diametral lens (dashed line) . . . . . . . . . . 112.5 Local feature size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.6 Boundary Encroachment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.7 Watson Insertion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.8 Flow chart for refinement algorithm . . . . . . . . . . . . . . . . . . . . . . . 162.9 Boundary discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.10 Constrained Delaunay trianglulation before boundary recovery . . . . . . . . . 182.11 Boundary discretization based on TVT split . . . . . . . . . . . . . . . . . . . 192.12 New boundary discretization (after second insertion) . . . . . . . . . . . . . . 212.13 New boundary discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.14 Uniform-degree edge swapping . . . . . . . . . . . . . . . . . . . . . . . . . . 222.15 The minimum sine value φ(x,y) . . . . . . . . . . . . . . . . . . . . . . . . . 242.16 Contours of φ(x,y) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242.17 The active set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.18 Neighboring CV stencil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28viiList of Figures3.1 Delaunay triangulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.2 Advancing-front illustration . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.3 Engwirda’s frontal-Delaunay triangulation . . . . . . . . . . . . . . . . . . . . 333.4 shape-optimal vertex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.5 Size-optimal vertex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.6 Choosing between the three kinds of vertices . . . . . . . . . . . . . . . . . . 353.7 Illustration of Engwirda’s frontal-Delauany method . . . . . . . . . . . . . . . 394.1 Marcum’s original scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 414.2 AFLR mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 434.3 Layer index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 434.4 Four potential vertices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.5 Illustration of modified AFLR scheme . . . . . . . . . . . . . . . . . . . . . . 455.1 Meshes generated on a square domain . . . . . . . . . . . . . . . . . . . . . . 505.2 The process of Engwirda’s scheme . . . . . . . . . . . . . . . . . . . . . . . . 515.3 The process of AFLR scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 525.4 Angle distribution for meshes on square domain . . . . . . . . . . . . . . . . . 535.5 Sinusoidal manufactured solution of the Poisson equation . . . . . . . . . . . . 545.6 Mesh regularity, truncation error, and discretization error (second order) . . . . 555.7 Order of Accuracy for Second Order Discretization . . . . . . . . . . . . . . . 565.8 Order of Accuracy for the Fourth Order Discretization . . . . . . . . . . . . . 565.9 Mesh generation time on a square domain . . . . . . . . . . . . . . . . . . . . 585.10 Time cost of the mesh generation and the solver . . . . . . . . . . . . . . . . . 585.11 Meshes in a channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 605.12 Angle distribution of meshes in a channel . . . . . . . . . . . . . . . . . . . . 615.13 Bad AFLR mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 615.14 AFLR scheme failure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62viiiList of Figures5.15 Near view of the transition area . . . . . . . . . . . . . . . . . . . . . . . . . . 625.16 Exact solution for advection-diffusion equation . . . . . . . . . . . . . . . . . 645.17 Order of Accuracy for second order discretization . . . . . . . . . . . . . . . . 645.18 Order of Accuracy for fourth order discretization . . . . . . . . . . . . . . . . 655.19 Meshes around a circular cylinder . . . . . . . . . . . . . . . . . . . . . . . . 675.20 Angle distribution of meshes around a circular cylinder . . . . . . . . . . . . . 685.21 Exact solution for the non-dimensional pressure around a circular cylinder . . . 705.22 Truncation error for the energy equation . . . . . . . . . . . . . . . . . . . . . 715.23 Discretization error for the pressure . . . . . . . . . . . . . . . . . . . . . . . 715.24 Time cost for both mesh generator and solver (fourth order) . . . . . . . . . . . 725.25 Regularity for meshes around NACA-0012 airfoil . . . . . . . . . . . . . . . . 755.26 Angle distribution of meshes around NACA-0012 airfoil . . . . . . . . . . . . 765.27 NACA-0012 airfoil with angle of attack 2◦ . . . . . . . . . . . . . . . . . . . . 765.28 NACA-0012 airfoil with angle of attack 4◦ . . . . . . . . . . . . . . . . . . . . 775.29 Time cost of mesh generation and solver in total . . . . . . . . . . . . . . . . . 785.30 Time cost in the fourth order solver for AoA 4◦ . . . . . . . . . . . . . . . . . 795.31 Time cost in the second order solver for AoA 2◦ . . . . . . . . . . . . . . . . . 80ixList of Algorithms2.1 Initial Queuing algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.2 Quality Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.3 Uniform degree swapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.1 Engwirda’s frontal-Delaunay scheme . . . . . . . . . . . . . . . . . . . . . . . 363.2 Engwirda’s frontal-Delaunay Queuing scheme . . . . . . . . . . . . . . . . . . 363.3 Quality Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37xAcknowledgementsI sincerely thank my research supervisor Dr. Carl Ollivier-Gooch for his patience, encour-agement, and guidance in helping me to accomplish my graduate study. He is an awesomesupervisor and mentor.Also I would like to thank Gary Yan, Daniel Zaide, and the rest of my colleagues in theANSLab research group for their suggestions and advice, which are inspiring and valuable tomy study.xiDedicationI am grateful for my parents’ unconditional love and support. Their encouragement and under-standing make me stronger.xiiChapter 1IntroductionComputational Fluid Dynamics, or CFD, is one of the three major methods in the study of flu-ids, along with theoretical and experimental analysis. CFD simulates the physical phenomenain the flow of fluids numerically [21]. These phenomena, such as shockwaves and the turbu-lence around a wing, are governed by the Navier-Stokes equations. In most cases, the modelsof interest are non-linear and have no analytical solutions. CFD allows an alternative for theunderstanding of complex flows via numerical simulations instead of physical experiments.Compared to experimental and theoretical fluid analysis, CFD is young and dramatically grow-ing [2]. In the 1960s, digital computers began to be available, which made possible numericalsimulations. In the following 50 years, CFD has grown with the development of more pow-erful computers and algorithms. Now, CFD has been applied to multiple disciplines such asfluid mechanics, electrical and electronic engineering, meteorology, biomedical engineeringand more.The complete CFD process contains three major steps: pre-processing, flow solution andpost-processing.Pre-Processing: This stage prepares all the information the flow solver needs. Given a partic-ular problem, we choose a sufficient and representative model of the physics. Differentmodels have different approximations. The approximation introduces physical modelingerror, which means the flow of interest is not analyzed perfectly.The geometry of the computational domain usually is defined by computer aided de-sign (CAD) software. Then the domain is divided into finite non-overlapping cells orcontrol volumes. This division process is known as mesh generation. It is also the focus1Chapter 1. Introduction(a) Structured mesh (b) Unstructured meshFigure 1.1: Meshes around a circular cylinder (Near View)of this study. Generally speaking, the more cells there are in the domain, the more ac-curate the solution will be. However, finer meshes also require higher computational re-sources. Besides the size of the mesh, the connectivity and the shape of the cells also haveimpacts on the accuracy. Based on how the cells are connected, as illustrated in Fig. 1.1,meshes can be classified as structured and unstructured meshes. The connectivity of thestructured mesh is predictable. As shown in Fig. 1.1a, the cell (i, j) always has adjacentcells marked as (i± 1, j± 1). This regular connectivity makes structured meshes easierto deal with in the flow solver. It works well in simple geometries. However, when itcomes to complicated geometries, it is not always possible to generate structured meshesautomatically. Therefore, unstructured meshes are more widely used in CFD for compli-cated geometries. Though the connectivity must be stored explicitly in the flow solver,which makes the codes more complicated and challenging to write, current high-orderaccurate algorithms make it worthwhile to adopt unstructured meshes on complicated21.1. Motivationgeometries. The mesh generation software used in the study is GRUMMP (Generationand Refinement of Unstructured, Mixed-Element Meshes in Parallel) [27].Flow Solution: In the pre-processing stage, we choose the mathematical governing equa-tions which, in fluid dynamics, are partial differential equations with some algebraicconstraints. To solve these equations numerically, we need to discretize the mathemat-ical equations into algebraic equations within each cell. After decades of development,the current three major families of discretization schemes are finite-difference [21, 1],finite-element [30], and finite-volume schemes [37]. The finite-volume scheme is usedin this study because of its advantages both in conservation of physics quantities and itseasy application on unstructured meshes [21]. The algebraic equations will be solvednumerically instead of the original mathematical equations. The approximation in thediscretization both in space and time introduces truncation errors. The solver stops whenthe residual of the discrete equation is close to zero. The difference between the numer-ical solution and the exact solution is known as discretization error. In this study, we useour in-house code ANSLib (Advanced Numerical Simulation Library) [26] as the flowsolver.Post-Processing: This stage is the visualization and analysis of the numerical solution ob-tained from the flow solution stage, which includes geometry and mesh display, vectorand contour plots, error estimation on the variables of interest and more. Visualization isusually completed with the help of plotting tools like Paraview [3], MATLAB [23], andgnuplot [38].1.1 MotivationAs mentioned in the previous section, the accuracy of CFD is evaluated by how we control thephysical modeling error and the numerical error. The latter is the target in this study. Usually,the more cells in the mesh, the more accurate we expect the solution to be. However, finer31.1. MotivationFigure 1.2: Simple 1-D examplemeshes introduce larger costs. The cell shape also influences the solution accuracy. We takethe one-dimensional case as a simple example to illustrate the concept. The governing equationfor steady heat conduction isd2Tdx2= 0. (1.1)Integrating the governing equation in cell i and applying Gauss’s theorem, we get1hi(dTdx∣∣∣∣i+ 12− dTdx∣∣∣∣i− 12)= 0. (1.2)The solution in cell i is represented by a reconstruction polynomialT (x) = Ti+(x− xi) dTdx∣∣∣∣i+(x− xi)22d2Tdx2∣∣∣∣i+(x− xi)36d3Tdx3∣∣∣∣i+ · · · , (1.3)where xi indicates the location of the cell reference point. As shown in Fig. 1.2, if we choosethe cell center of the control volume (CV) i as the cell reference point and integrate the recon-struction over the cells, the cell averages are,T i =1hiˆCViT (x)dx = T |i+h2i24d2Tdx2∣∣∣∣i+O(h4i )41.1. MotivationT i+1 =1hi+1ˆCVi+1T (x)dx= T |i+(12hi+12hi+1)dTdx∣∣∣∣i+(18h2i +14hihi+1+16h2i+1)d2Tdx2∣∣∣∣i+(148h3i +116h2i hi+1+112hih2i+1+124h3i+1)d3Tdx3∣∣∣∣i+O(h4i )T i−1 =1hi−1ˆCVi−1T (x)dx= T |i+(−12hi− 12hi−1)dTdx∣∣∣∣i+(18h2i +14hihi−1+16h2i−1)d2Tdx2∣∣∣∣i+(− 148h3i −116h2i hi−1−112hih2i−1−124h3i−1)d3Tdx3∣∣∣∣i+O(h4i ).The first order surface gradients are,dTdx∣∣∣∣i+ 12=T i+1−T ihi+1+hi2(1.4)dTdx∣∣∣∣i− 12=T i−T i−1hi−1+hi2. (1.5)Substitute Eqn. 1.4 and Eqn. 1.5 into Eqn. 1.3,1hi(dTdx∣∣∣∣i+ 12− dTdx∣∣∣∣i− 12)=d2Tdx2+hi+1−2hi+hi−13hid2Tdx2+ (1.6)−hihi−1+hihi+1−h2i−1+h2i+112hid3Tdx3+O(h2).Comparing Eqn. 1.6 to Eqn. 1.1, we see that our discretization is not exact. This error, calledthe truncation error, is zeroth-order for arbitrary choices of cell sizes. On the other hand, if themesh is regular, or hi+1 = hi = hi−1, the truncation error is second-order. Smooth variations ofh also lead to second-order truncation error. Juretic´ and Gosman proved that the cell shape alsoaffects the truncation error in 2-D and 3-D cases [19]. The prediction of error in the solution,or the discretization error, and the relations between truncation and discretization errors are51.1. Motivationcomplex [32, 7].Katz et al. [20] showed that mesh regularity affects truncation and discretization errors.The truncation error obtained from the regular mesh is higher order than those from perturbedmeshes. As for the solution error, the advantage in the order of accuracy is only observed ininviscid fluxes. The viscous fluxes are not sensitive to mesh regularity in the sense of having aneffect on solution error. Diskin et al. [8] also found that the relations between mesh character-istics and the solution accuracy are complicated. Mesh irregularities affect gradient errors anddiscretization errors in different ways. The viscosity affects the sensitivity to mesh regularity,as does the solution reconstruction method. For isotropic meshes, mesh regularity does notaffect gradient error much if the unweighted least squares method is used. However, if usingthe Green-Gauss method, the regular quadrilateral mesh offers a higher order of accuracy thanothers in his test. In addition, Jalali et al [16, 18] and Yan et al [39] investigated the accuracy ofdifferent discretization schemes on a wide range of meshes from perfect to perturbed meshesfor cell-centered and vertex-centered control volumes respectively. The degradation on the or-der of accuracy for truncation errors on perturbed meshes is observed in their research. Theyfound that a particular discretization scheme may minimize the truncation error for a certainkind of mesh.Most research on mesh regularity and errors are applied for perfect meshes and their per-turbations for simple geometric domains. This study focuses on how regular meshes can beautomatically generated for more complicated geometries, and how that affects the truncationand discretization errors. In particular, this thesis seeks to answer the following questions:1. How regular a triangular mesh can we generate for realistic cases?2. What is the additional cost to generate such a mesh?3. Considering that the cell shape affects the accuracy greatly, can we expect small dis-cretization error using unstructured regular meshes?4. Can a coarse regular mesh offer the same accuracy as a finer irregular mesh with lower61.2. Outlinecosts in the flow solver?1.2 OutlineIn the next chapter, we will recall some background information on both mesh generationtechniques and unstructured mesh finite volume schemes. The Delaunay mesh refinementalgorithm is discussed. We will also introduce post-processing techniques for mesh generationsuch as smoothing and swapping.In Chapter 3 and 4, the methodology to generate more regular meshes in this study isdescribed. Two mesh generation techniques are used in the study: Engwirda’s frontal-Delaunayalgorithm and Marcum’s advancing front local reconnection algorithm. Both are proven toimprove the mesh regularity dramatically. Several changes are made in these two schemes towork on our current GRUMMP framework.In Chapter 5, several numerical test cases are applied. In each case, the meshes generatedby different mesh generation techniques are shown and discussed first. The comparison of themesh regularity shows the advantage of new mesh generation techniques. The truncation errorand discretization error are also visualized to illustrate the effects of the mesh regularity. Thecost is shown at last to discuss whether it is worth it to generate more regular meshes.In Chapter 6, the thesis concludes with the summary of this research.7Chapter 2BackgroundDelaunay mesh refinement is chosen as a base-line mesh generation scheme because of itstheoretical guarantees and its practical performance. It will be introduced in Section 2.1. Post-processing techniques such as edge swapping and vertex smoothing, which help improve themesh quality, will be discussed in Section 2.2. Section 2.3 provides a brief review of the finite-volume method.2.1 Delaunay mesh refinementIn two dimensions, a triangulation is a set of triangles that connect all given vertices. Trianglescannot intersect each other and all triangles collectively must fill the polygonal boundary of thedomain. A Delaunay triangulation is the triangulation in which all triangles have empty cir-cumcircles, which means that there are no vertices lying inside circumcircles1 but the verticesare allowed to be located on circumcircles. Fig. 2.1 shows a simple example of a Delaunaytriangulation. All the circumcircles (dashed lines) are empty. If no four vertices lie on thesame circle, the Delaunay triangulation is guaranteed to be a unique triangulation. The proofcan be found in Reference [6, 9]. However, a Delaunay triangulation may break the domainboundary in some cases. As shown in Fig. 2.2, the boundary is marked with bold lines. Be-cause one boundary point fails to meet the empty circumcircle criterion (left), the Delaunaytriangulation (right) of these vertices changes the domain boundary. To solve this problem, theconstrained Delaunay triangulation was proposed [5]. Unlike the Delaunay triangulation which1The circumcircle is a triangle’s circumscribed circle. It is the unique circle that passes through all the tri-angle’s three vertices. The center of the circumcircle is called the circumcenter, and its radius is called thecircumradius.82.1. Delaunay mesh refinementFigure 2.1: Delaunay TriangulationFigure 2.2: Delaunay Triangulation breaks the boundary92.1. Delaunay mesh refinement(a) PSLG (Eight points and one segment) (b) Constrained Delaunay Triangulation (BoundaryProtected)Figure 2.3: Constrained Delaunaystarts with a set of vertices, the constrained Delaunay triangulation starts with a planar straightline graph (PSLG). The PSLG contains both vertices and segments, as shown in Fig. 2.3a. Ev-ery endpoint of the segments must be among the vertices in this PSLG. The segments are onlypermitted to intersect at endpoints. In a constrained Delaunay triangulation, the segments inPSLG serve as edges of the triangles. Only vertices that are visible to a triangle are consideredwhen checking empty circumcircles. Visibility means that no input segment stands betweenthe vertex and the triangle. By setting boundary edges as segments in the PSLG, we can pro-tect them from breaking by the Delaunay triangulation. Fig. 2.3b illustrates how a constrainedDelaunay triangulation conserves boundaries. The red vertex is located inside the circumcircleof a triangle. However, the segment in the PSLG (marked as bold line) separates the vertexfrom that triangle. Therefore the red vertex is not visible to the triangle whose circumcirclecontains it. This is not regarded as a violation of the empty circumcircle rule so the boundaryis conserved.The algorithm used in GRUMMP to generate guaranteed-quality triangular meshes [4, 28]is based on Ruppert’s scheme [33]. Ruppert’s scheme starts with the constrained Delaunay tri-angulation. The Delaunay triangulation maximizes the minimum angles in meshes but usually,bad quality triangles still exist. The quality of the triangulation is improved by inserting pointsinto the constrained Delaunay triangulation. Ruppert’s scheme inserts circumcenters of the bad102.1. Delaunay mesh refinement30° 30°Figure 2.4: Diametral circle (dotted line) and diametral lens (dashed line)quality triangles into the mesh unless the circumcenters encroach boundaries. Encroachmenthappens when a point is located inside a circle with a boundary edge as its diameter. Ruppertshowed that the algorithm terminates with the minimum angle in the mesh larger than 20.7◦.Shewchuk showed that when using diametral lenses instead of diametral circles to determineencroachment, the minimum angle in the mesh is greater than 25.7◦ [34]. The difference be-tween diametral lenses and diametral circles is illustrated in Fig. 2.4.To generate a guaranteed-quality triangular mesh, both element size and shape need to beconsidered. The size is controlled by length scales. If the circumradius of a triangle is lessthan the average length scales of its three vertices, that triangle satisfies the size constraint.GRUMMP uses a user controlled refinement parameter R to determine how small the elementsare. Small triangles are usually not required everywhere in the mesh. In the far field or theareas where variables of interest tend to be constant, triangles are permitted to be much larger.Therefore, a grading parameter G is introduced to control how fast the edge length grows overdistance. The length scale (LS) for point p is defined with these two parameters,LS(p) = min(lfs(p)R, minneighbours qiLS(qi)+1G|qi− p|). (2.1)112.1. Delaunay mesh refinementFigure 2.5: Local feature sizeqi denotes the neighbour points to point p. In Fig. 2.5, local feature sizes (lfs) for points (redpoints) are the radii of those circles. The lfs is defined as the radius of the smallest circlecentered at point p that touches two disjoint parts (vertices or edges) of the domain boundary.A large value of R reduces the target edge length while a larger value of G leads to a slowerincrease in element size over the same distance. As mentioned before, the minimum angle inconstrained Delaunay triangulation is guaranteed to be larger than 25.7◦ if the diametral lens isused for encroachment. A triangle is regarded as good in shape when the minimum angle of itis larger than 25.7◦. In GRUMMP, we use the normalized shortest edge length to circumradiusratio to evaluate the shape. The shape quality ρ isρ =`m√3r.The quality varies from 0 to 1. An equilateral triangle enjoy the maximum shape quality 1and a degenerate triangle has minimum shape quality 0. According to the law of sines, theminimum edge length `m, the minimum angle αm, and the circumradius r have the relationship122.1. Delaunay mesh refinementbelow,`msin(αm)= 2r.Therefore the shape quality can also be expressed asρ = 2sin(αm)/√3. (2.2)When the diametrial lens is used for encroachment, the minimum shape quality allowed in themesh is:ρ =2sin(25.7◦)√3≈ 0.5.If the diametral circle is used instead, the minimum shape quality is about 0.4. Only thosetriangles which satisfy both size and shape constraints are considered as good quality elements.The process of inserting vertices into an initial constrained Delaunay triangulation to im-prove the mesh quality is called mesh refinement. In GRUMMP, all bad triangles in the meshare queued initially. The triangle with the worst quality comes at the top of the queue. Theinsertion procedure grabs the top triangle in the queue and calculates the circumcenter of thatbad triangle as the location for a new vertex to insert. If the newly inserted vertex encroachesboundaries, the vertex is rejected and instead, the encroached boundary edge is divided in half,as shown in Fig. 2.6. If not, we use Watson’s method [31] to insert the vertex into the mesh.All triangles whose circumcircles contain that newly inserted vertex are deleted, leaving a hullwith that new vertex in it. Connecting faces of the hull with the newly inserted vertex formsseveral new triangles. New triangles with bad quality need to be queued for later splitting whilegood ones do not. The refinement procedure terminates when the queue is empty. Fig. 2.7 illus-trates Watson insertion when a new vertex does not encroach boundaries. The gray triangle inFig. 2.7a is the bad one to be refined, whose circumcenter is marked as red. The triangles withbold lines are those whose circumcircles contain that red vertex. The hull and the triangulationafter Watson insertion are shown in Fig. 2.7b and Fig. 2.7c.The priority queue often plays an important role in the mesh quality. In GRUMMP, badly132.1. Delaunay mesh refinement(a) The newly inserted vertex encroaches the bound-ary (bold line)(b) The boundary is split in halfFigure 2.6: Boundary Encroachment(a) The circumcenter of the bad tri-angle is calculated(b) Deleting the triangles leaving thehull(c) New triangles generatedFigure 2.7: Watson Insertion142.1. Delaunay mesh refinementAlgorithm 2.1 Initial Queuing algorithmi = 0while i <number of triangles in the mesh doif Triangle Ti is still in the mesh and fails quality test thenInput Ti with priority from quality test into the queueend ifi = i+1end whileAlgorithm 2.2 Quality TestCalculate the circumradius RCalculate the average length scale l¯ of three vertices in TFind the shortest edge length lminif R > l¯ thenQuality =−Rl¯elseQuality = lmin√3Rif Quality > WorstShapeAllowable thenQuality = NeverQueueend ifend ifsized elements, the ones with the large circumradius to averaged length scale ratio, come beforebadly shaped elements, largest first. If triangles satisfy the size constraint, the shape qualityis considered. The triangle with the smallest shape quality comes at the top of the queue.Algorithm 2.1 and Algorithm 2.2 illustrate how we deal with the priority queue in GRUMMP.Triangles with “NeverQueue” quality will not added to the queue. Other details such as howto handle small angles in the domain can be found in Reference [4]. The flow chart for thisscheme is shown in Fig. Delaunay mesh refinement Constrained Delaunay triangulation with well represented boundaries   Queue all cells  Pick the first cell in the queue  Does the circumcenter encroach boundaries? Split the boundary and connect new triangles Watson insertion Queue new triangles Empty queue? The End Yes No Yes No Figure 2.8: Flow chart for refinement algorithm162.1. Delaunay mesh refinementMesh generation for domains whose boundaries are straight lines is straightforward. Bound-aries can be simply split into edges with the user-controlled size. The constrained Delaunaytriangulation with the boundary discretized serves as the starting triangulation of the refine-ment algorithm discussed before. However, the discretization of curved boundaries is morechallenging. As illustrated in Fig. 2.9, for high curvature curves, they should be discretizedinto more edges locally to capture local curvature features. However, the fine discretization isnot required for the low curvature parts on curves. The algorithm for discretizing both parts indifferent ways automatically is not straightforward. In GRUMMP, the technique is based on atotal variation of tangent angles (TVT). The details can be found in Reference [4]. The TVTboundary discretization scheme starts with a very coarse discretization. Fig. 2.10 shows thestarting boundary edges in blue lines while the actual geometry curve is marked as red. Afterqueueing all the starting boundary edges, pick the first one in the queue and split it if the edgelength is larger than the local length scale. If the local curvature is small, the edge is split at thegeometric average location. Otherwise the split location is calculated based on TVT so that thechange in normal direction is divided equally between the two new segments. The length scalefor a point p on a curved boundary is calculated by Eqn. 2.1, where local feature size shouldbe modified with curvature,lfsc(p) = min(ρ(p), lfs(p)).ρ denotes to the radius of curvature. The boundary discretization based on this approach isshown in Fig. 2.11. The TVT split scheme successfully uses more edges around the leadingedge to capture high curvature curve. The local feature size at the trailing edge is large in thestarting triangulation. Therefore there is not much for the TVT split scheme to do. The verycoarse discretization around the trailing edge is not good news because it may fail to capturerapidly changing solution features. This can be fixed by manually setting a small length scaleat the trailing edge but there is still a problem. Usually, the curvature is nearly constant aroundthe trailing edge so boundary edges are simply split in half. In this scheme, there are sometransition areas where the length of one edge is twice as large as the adjacent one. This makes172.1. Delaunay mesh refinementFigure 2.9: Boundary discretizationFigure 2.10: Constrained Delaunay trianglulation before boundary recovery182.1. Delaunay mesh refinementFigure 2.11: Boundary discretization based on TVT splitit difficult to generate equilateral triangles.The improved boundary discretization algorithm is introduced to deal with trailing edgesof airfoils [40]. The new scheme starts at end point of an airfoil, which is usually the trailingedge. Given a parametric curve, C(t), the length at parameter t is marked as `(t). We denotetotal curve length as ¯`. The average length scale is easily calculated asLS(t) =¯`Rb,where Rb is the refinement parameter to decide how fine the boundary would be discretized.Approximately, it denotes the number of edges along the curve. Boundary grading parameterGb is also needed to determine the boundary edge length growth rate. Gb is the ratio of thelargest boundary edge length to the smallest one. A large value of Gb results in a small lengthscale at trailing edge. Based on the average length scale, we would expect length scales aresmaller near end points or parts with high curvature. For airfoils, these would be trailing edgesand leading edges. Within a certain fraction (In this study, α = 0.25) of total length distance to192.2. Post-processingend points, which denotes the region of trailing edges,min(`(t), ¯`− `(t))< `end = α ¯` α ∈ (0,0.5),smaller length scales are applied,LS(t) =12( ¯`Rb− 1Gb¯`Rb)(sin(pi`(t)`end−pi/2)+1)+1Gb¯`Rb.This also makes the edge length growth near the trailing edge more smooth to fix the problemthat one edge length is twice large as the adjacent one in the previous scheme. Near the leadingedge, if ρ(t)< ¯`Rb , it is regarded as a high curvature area. The smaller length scale for this areaisLS(t) = max(ρ(t),¯`Rb1Gb).Instead of dividing current edges recursively as used in the TVT boundary discretization scheme,the new algorithm discretizes the boundary by marching from one end to another. The startingparameter is 0. Insert the point P1 with a distance LS(0) = ¯`/(Rb ·Gb) from the end point P0on a curve. The parameter t1 of this distance is obtained from the parametric curve. Point P2 isinserted with a distance LS(t1) from P1. The first two insertions are shown in Fig. 2.12. Keepinserting points until the point to be inserted is too close to the other end point. In my study,distance less than 0.5LS(tmax) is regarded as the termination of the insertion. Fig. 2.13 shows asmooth boundary discretization with this new scheme. Both the leading edge and trailing edgeare treated properly.2.2 Post-processingAfter the refinement, when all triangles have good shapes and sizes, post-processing proceduressuch as swapping and smoothing are applied to improve mesh quality.202.2. Post-processingFigure 2.12: New boundary discretization (after second insertion)Figure 2.13: New boundary discretization212.2. Post-processingFigure 2.14: Uniform-degree edge swappingAlgorithm 2.3 Uniform degree swapping1: Queue all faces2: while The queue is not empty do3: Grab the first face in the queue and its two adjacent triangles4: Find degrees of four vertices and calculate their variance from six5: Calculate the variance after swapping6: Remove the face from the queue7: if The variance after swapping is smaller then8: Do swapping9: Add the new diagonal to the queue10: end if11: end while2.2.1 SwappingIn the triangular mesh, two adjacent triangles have four different vertices. If they form a convexquadrilateral, there are two possible configurations. Based on a certain criterion, if the otherconfiguration is better than the current one, it can be easily acquired by changing the diagonal.A regular mesh consists of equilateral triangles, where all vertices have six incident cells.Therefore in this study, we use uniform-degree edge swapping where the degree of a vertexis the number of incident cells as shown in Fig. 2.14. The total variance of degree for six isused as the criterion. The optimal variance is zero when all four vertices have six neighbours.The pseudo code for uniform degree swapping is shown in Algorithm 2.3. In Fig. 2.14, thevariance before swapping is 4. After the swapping, all vertex degrees are six and the varianceis 0, which is smaller than before swapping. Therefore the swap is performed.222.2. Post-processing2.2.2 SmoothingA number of smoothing methods have been used to improve mesh quality. These approachescan be classified as local and global smoothing techniques. A local method adjusts the geomet-ric location of one vertex at a time to achieve the optimal shape quality in a neighborhood. Theoverall mesh quality improves by applying the local smoothing for every adjustable vertex inthe mesh. It is efficient if each local adjustment is inexpensive. The global smoothing methodadjusts all vertex locations in the same time. It needs to solve an optimization problem as largeas the number of vertices in a mesh, which is computationally expensive.One of the commonly-used local smoothing techniques is Laplacian smoothing [11]. Itrelocates the free vertex to the arithmetic mean location of its incident vertices,xfree =∑i=Ni=0 xiNyfree =∑i=Ni=0 yiN.N refers to the number of incident vertices to the free vertex and x, y are the spatial coordinates.This method is computationally inexpensive, but it does not guarantee any improvement in theelement quality. Actually, it is possible to produce inverted elements which have negativevolume.In GRUMMP, we use the optimization-based local smoothing technique proposed by Fre-itag et al. [12]. Optimization techniques use functions and their gradients to find the positionwhere the functions obtain optimal values. Only a brief review of this method is shown hereand the details can be found in Reference [13, 14]. To get a regular triangular mesh, we con-sider the optimal location for a free vertex to be where the minimum sine value is maximized.The challenge is that the function to be optimized as shown in Fig. 2.15,φ(x,y) = mini∈Sfi(x,y) = mini∈Ssin(θi(x,y)),is only piecewise smooth and therefore is not differentiable everywhere. In this case, a simple232.2. Post-processing010.20.40.5(x,y)0.6y0.80 0.5x0-0.5-0.5- 2.15: The minimum sine value φ(x,y)Figure 2.16: Contours of φ(x,y)242.3. Finite-volume solver(a) Starting location (b) Two angles in the active setFigure 2.17: The active setregular hexagon topology as shown in Fig. 2.16 is used to facilitate the illustration. We denotethe active set A as the set of minimum angles when the point is located at (x,y). For example,in Fig. 2.17a, only one angle is in the active set, in other words, there is only one angle hasthe smallest sine value in this case. The first search step attempts to find the location wheremore than one angles are in the active set as shown in Fig. 2.17b. For this simple case, the twoactive angles set the trajectory straight up. The next step stops where the active set changes.The approximate route is shown in Fig. Finite-volume solverThe finite-volume solver starts with the conservative form of the governing equations, whichcan be cast into a generic form [2]∂U∂ t+∇ ·F = S (2.3)where U denotes the solution vector, F the flux vector, and S the source term.Integrating Eqn. 2.3 over an arbitrary control volume (CV) i and applying Gauss’s theorem,we obtain ¨CVi∂U∂ tdA+˛∂CViF · nˆds =¨CViSdA. (2.4)252.3. Finite-volume solverAssuming the discretized physical domain to be constant, Eqn. 2.4 can be further simplified asdU idt+1Ai˛∂CViF · nˆds = 1Ai¨CViSdA (2.5)where Ai denotes the area of the control volume, U i = 1Ai˜CViUdA the control volume average,and nˆ the outward unit normal vector.The general finite-volume numerical method can be summarized in the following stages [21]:1. Approximate U(x,y,z) in each control volume with a polynomial U i(x,y,z)=Ui+ ∂U i∂x (x−xi)+ ∂U i∂y (y− yi)+ ∂U i∂ z (z− zi)+ · · · . Given the value of U i for each control volume, per-form solution reconstruction to obtain the polynomial coefficients (Ui,∂U i∂x ,∂U i∂y ,∂U i∂ z , · · · ).Using this polynomial approximation of U , find U at the control volume boundary, andevaluate flux F at the boundary.2. Since there is a distinct approximation in each control volume, two distinct values ofthe flux are generally obtained at the boundary between two control volumes. Applysome strategy for resolving the discontinuity in the flux at the control volume boundaryto produce a single value of the flux.3. Integrate the flux along the control volume boundary.4. Advance the solution in time to obtain new values of U .In this section, the solution reconstruction is further discussed. The details of other algorithmssuch as flux evaluation and flux integration can be found in Reference [2, 29].A two-dimensional solution reconstruction is used as an illustration. The reconstructedsolution U˜i(x,y) is obtained by a Taylor expansionU˜i(x,y) =U |i+ ∂U∂x∣∣∣∣i(x− xi)+ ∂U∂y∣∣∣∣i(y− yi)+∂ 2U∂x2∣∣∣∣i(x− xi)22+∂ 2U∂x∂y∣∣∣∣i(x− xi)(y− yi)+ ∂2U∂y2∣∣∣∣i(y− yi)22+ · · · (2.6)262.3. Finite-volume solverwhere xi and yi denote the locations of the reference point of the control volume i.Conservation of the mean within the control volume i requires thatUi =¨CViU˜i(x,y)dA. (2.7)Substituting Eqn. 2.6 into Eqn. 2.7, we get thatU i =U |i+ ∂U∂x∣∣∣∣ixi+∂U∂y∣∣∣∣iyi+∂ 2U∂x2∣∣∣∣ix2i2+∂ 2U∂x∂y∣∣∣∣ixy+∂ 2U∂y2∣∣∣∣iy2i2+ · · · (2.8)where the geometric moments are represented asxmyni =1Ai¨Vi(x− xi)m(y− yi)ndA (2.9)Accuracy of the reconstruction is determined by the number of derivatives evaluated in thepolynomial. To make it k-exact, or (k+ 1)-order accurate, the polynomial has to be degree k.For instance, a second-order accurate reconstruction,U˜i(x,y) =U |i+ ∂U∂x∣∣∣∣i(x− xi)+ ∂U∂y∣∣∣∣i(y− yi)+O(∆x2,∆y2),requires U |i, ∂U∂x∣∣∣i, and ∂U∂y∣∣∣ito be evaluated.To compute these derivatives, we seek to minimize the error in predicting the mean valueof the function for neighbouring control volumes. The mean value for the control volume CVjis,1A j¨CV jU˜idA = U |i + ∂U∂x∣∣∣∣i1A jˆCV j(x − xi)dA + ∂U∂y∣∣∣∣i1A jˆCV j(y − yi)dA + · · · . (2.10)To avoid computing moments of each control volume, we replace x− xi with x− x j + x j− xi272.3. Finite-volume solver2nd2nd2nd2nd3rd3rd3rd3rd3rd3rd4th4th4th4th4th4th4th4th 4thFigure 2.18: Neighboring CV stenciland y− yi with y− y j + y j− yi. Eqn. 2.10 can be represented as1A j¨CV jU˜idA =U |i+ ∂U∂x∣∣∣∣ixˆi j +∂U∂y∣∣∣∣iyˆi j + · · · (2.11)withx̂myni j =1A j¨CV j((x− x j)+(x j− xi))m · ((y− y j)+(y j− yi))ndA=n∑l=0m∑k=0n!l!(n− l)!m!k!(m− k)!(x j− xi)k · (y j− yi)l · xm−kyn−l jThe minimum number of neighbouring control volumes in the stencil is equal to the numberof derivatives to be evaluated. In practice, we choose three neighbors for the second-orderaccuracy and nine for the third-order, as shown in Fig. 2.18. The least-squares system for282.3. Finite-volume solverevaluating derivatives is,1 x¯i y¯i x2i xyi y2i · · ·wi1 wi1x̂i1 wi1ŷi1 wi1x̂2i1 wi1x̂yi1 wi1ŷ2i1 · · ·wi2 wi2x̂i2 wi2ŷi2 wi2x̂2i2 wi2x̂yi2 wi2ŷ2i2 · · ·wi3 wi3x̂i3 wi3ŷi3 wi3x̂2i3 wi3x̂yi3 wi3ŷ2i3 · · ·.................. . . .wiN wiN x̂iN wiN ŷiN wiN x̂2iN wiN x̂yiN wiN ŷ2iN · · ·U∂U∂x∂U∂y∂ 2U2∂x2∂ 2U∂x∂y∂ 2U2∂y2...i=U¯iwi1U¯1wi2U¯2wi2U¯2...wiNU¯Nwhere geometric weights wi j are used to specify the relative importance of good prediction forneighboring control volumes. The weights are based on the distance between control volumereference points, wi j = 1|xi−x j| . The first row of this least-squares system is the mean con-straint as in Eqn. 2.8. It is eliminated by Gaussian elimination resulting in an unconstrainedleast-squares system which is then solved for every control volume using the Singular-Value-Decomposition (SVD) method. After the derivatives are evaluated, the flux evaluation andintegral [29, 36] are performed.29Chapter 3Frontal-Delaunay triangulationThe regularity of the mesh obtained from the Delaunay triangulation is not satisfying enough.Even with the help of swapping and smoothing, there are still lots of triangles that are far fromequilateral in the mesh, as shown in Fig. 3.1.Another competing two-dimensional unstructured mesh generation scheme, the advancing-front method, begins with a discretization of boundaries. The boundary edges form the initialfront. A particular edge from the front is selected as the base edge, and a new triangle is createdwith this base edge and a newly created vertex or an existing vertex. The new vertex is insertedto make an optimal triangle. After insertion, the base edge is removed from the front since itis obscured by the new triangle, and the newly generated edges are either added to or removedfrom the front based on their visibility. This procedure terminates when the front is closed.Fig. 3.2 gives a simple illustration of the advancing-front scheme. The front is marked by thesolid lines in the figures. As shown in Fig. 3.2a, the first triangle is generated. Two newlygenerated edges are added to the front while the boundary edge on that triangle is removedfrom the front because it is no longer visible for other edges in the front. Fig. 3.2b shows thefirst case when the insertion is not valid. The optimal triangle to be generated is marked asred but the new edges are intersected with the existing ones. Therefore the scheme connectsthe red edge in Fig. 3.2a with the nearby existing vertex, resulting in the blue triangle. Withthe insertion continuing, the optimal triangle may require insertion of a vertex too close to anexisting vertex as shown in Fig. 3.2c. If the optimal triangle is generated nevertheless, a veryshort edge will be the side-effect. Therefore, the nearby existing vertex rather than the optimalvertex is connected with the base edge, forming the blue triangle in Fig. 3.2d. The advancing-30Chapter 3. Frontal-Delaunay triangulationFigure 3.1: Delaunay triangulationfront scheme terminates when the front is empty. Fig. 3.2e shows what the front (solid lines)looks like after the final insertion. If we choose any edge from the front, the new optimaltriangle will either intersect with or be too close to the existing triangles. Therefore, the blueedge in Fig. 3.2f is applied to close the front.Frontal-Delaunay algorithms [31, 24, 25] are the hybridization of Delaunay-refinement andadvancing-front techniques, in which the Delaunay triangulation is used to define the topol-ogy of a mesh while new points are inserted in a manner consistent with the advancing-frontmethod. Engwirda’s frontal-Delaunay method [10, 9] is applied in this study because it im-proves the mesh quality compared to the Delaunay triangulation, as shown in Fig. 3.3.Engwirda’s frontal-Delaunay method is similar in execution to the classic Delaunay re-finement algorithm, but uses different queuing and point placement strategies. The priorityqueue for bad triangles is sorted by radius-edge ratio, the shape quality (Eqn. 2.2) introducedin Delaunay-refinement triangulation. Only bad triangles whose shortest edges are smallerthan 1.5 times the local feature sizes are added to the queue. This constraint ensures refine-31Chapter 3. Frontal-Delaunay triangulation(a) After the first insertion (b) Rejection due to intersections(c) Rejection due to being too close to existingvertex(d) Connecting existing vertex(e) The front after the last insertion (f) Final meshFigure 3.2: Advancing-front illustration32Chapter 3. Frontal-Delaunay triangulationFigure 3.3: Engwirda’s frontal-Delaunay triangulationment occurs in a frontal fashion, because the triangle to be refined is guaranteed to be adjacentto a boundary or to a triangle that meets the quality criterion. Two new types of vertices to beinserted are used by Engwirda: the shape-optimal vertex and the size-optimal vertex. Both areoff-centers proposed by Üngor [35], which are located on the Voronoi diagram of the corre-sponding Delaunay triangulation.The shape-optimal vertex c′, as shown in Fig. 3.4, is placed so that 4ABc′ is the largestisosceles triangle that satisfies the bound on the minimum allowable angle αmin. c denotes thecircumcenter and the gray line denotes the Voronoi edge associated with the base edge (theshortest edge of the bad triangle). The altitude of the new vertex to the base edge is calculatedash′ =||emin||21tan(αmin/2),where ||emin|| denotes the minimum edge length. By positioning the new vertex at c′, the newlygenerated triangle just satisfies the minimum shape constraint.The size-optimal vertex c′′ generates a triangle satisfying the local size constraint. The33Chapter 3. Frontal-Delaunay triangulationFigure 3.4: shape-optimal vertexFigure 3.5: Size-optimal vertexqueuing algorithm of Engwirda’s frontal-Delaunay scheme guarantees that only the triangleswith their shortest edges satisfying size constraints are considered. Therefore, when calculatingthe position of c′′, only the two newly created edges need to be compared with the local featuresizes. An iterative bisection method is used to find the position of the new vertex. Take edge e1in Fig. 3.5 as an example. Pick the circumcenter and the middle point of the base edge as thestarting points c0, c1. Then calculate the local length scales of these starting points by Eqn. 2.1.DefineSi = ||e1||− 12(lfs(v1)− lfs(ci)).34Chapter 3. Frontal-Delaunay triangulationv1v2cemin c''c'Figure 3.6: Choosing between the three kinds of verticesIf S0 and S1 have the same sign, the size-optimal point is farther than the circumcenter, andits location is irrelevant as we will see. Otherwise choose the middle point of edge c0c1 asc2. Determine S2 and choose the interval [c0,c2] or [c2,c1] so that there is a zero inside theinterval. Iterate this procedure until Sn < 10−4||emin||. Next, repeat this algorithm for the edgee2. Choose the point closer to the base edge between this two vertices as the size-optimal vertexc′′. All three edges in the newly generated triangle4v1v2c′′ will satisfy the size constraint.Overall, we compare three candidate vertices: the circumcenter c, the shape-optimal vertexc′, and the size-optimal vertex c′′. Each of them will form a new triangle with the same shortestedge emin as the bad triangle. Of c, c′, and c′′, the one closest to the shortest edge is chosenfor insertion. Fig. 3.6 shows an example of which vertex should be inserted. In this case, size-optimal vertex c′′ is closer to the shortest edge than shape-optimal c′. Therefore, by choosingc′′ as the off-center, both the shape and size constraints are satisfied. Since the off-center hasto be on the Voronoi edge, the circumcenter is chosen if it is closer to the base edge than theoff-center. Otherwise, one of the off-centers is chosen, which is the case in Fig. 3.6.35Chapter 3. Frontal-Delaunay triangulationWith the Steiner vertex calculation method in mind, the basic algorithm of Engwirda’sfrontal-Delaunay method is shown in Algorithm 3.1. Rejection for a vertex too close to exist-ing ones in original advancing-front scheme is not required here. A proper quality bound issufficient to deal with this problem. In this techniques, if an unwanted small triangle is gen-erated, it will be regared as good sized triangle as in Algorithm 3.3. It will not be added tothe queue. Even if it is a badly shaped triangle, usually the next splitting will make it pass thequality test with the help of a carefully chosen Steiner vertex. Therefore a vertex inserted closeto existing ones is permitted and will not blow up the algorithm globally.Algorithm 3.1 Engwirda’s frontal-Delaunay scheme1: Starting with a constrained Delaunay triangulation2: for Every triangle in the mesh do3: Queue triangle4: end for5: while The queue is not empty do6: Take the top triangle in the queue7: if Triangle is not in the mesh then8: Continue9: else10: Calculate the location of the new vertex11: Use Watson insertion to insert the new vertex12: Queue the newly generated triangles13: end if14: end whileAlgorithm 3.2 Engwirda’s frontal-Delaunay Queuing scheme1: procedure QUEUE(Triangle)2: if QualityTest(Triangle) == BAD then3: Find the minimum edge length of that triangle4: if Minimum edge length < 1.5× local feature size then5: Calculate the normalized shape quality6: Add the triangle into the queue with shape quality as priority7: end if8: end if9: end procedureA simple example is used to illustrate the procedure of Engwirda’s frontal-Delaunay method.The geometry is a 20×20 square, we set the local feature size to 5 everywhere at the begining.36Chapter 3. Frontal-Delaunay triangulationAlgorithm 3.3 Quality Test1: procedure QUALITY TEST(Triangle)2: if Average length of three edges on that triangle > 1.5× local feature size then3: BadSizeTriangle4: else5: GoodSizeTriangle6: end if7: if Shape quality of the triangle < minimum allowable value then8: BadShapeTriangle9: else10: GoodShapeTriangle11: end if12: if BadSizeTriangle or BadShapeTriangle then13: return BAD14: else15: return GOOD16: end if17: end procedureThe minimum angle allowed in the mesh is 33◦. Therefore, in this mesh, any triangle withminimum angle smaller than 33◦ or average edge length greater than 7.5 would be regardedas a bad triangle. Also, only the triangles whose minimum edge lengths are smaller than 7.5are allowed to be added to the queue. To simplify the illustration, the size-optimal vertex isalways chosen as the insertion point, which means the newly generated triangles are equilateraltriangles. The refinement scheme starts with a constrained Delaunay triangulation as shown inFig. 2.10. All the boundary edge lengths are 5. Several triangles in the mesh are colored redor blue. Red triangles denote that they are badly shaped triangles while blue triangles denotethey are badly sized triangles. At this stage, only the red triangles are added into the queuebecause the minimum edge lengths of the blue triangles are too large to be queued. All the redtriangles have the same shape quality. They are all on the top of the queue. In the next step, anytriangle could be the one to be refined. In Fig. 3.7b, the hatched triangle is the newly generatedone. Watson insertion is applied. After connecting the new vertex to all the edges on the emptyhull, several new triangles are generated. Only triangles in the updated queue are shown in thefigure. Darker red indicates it has worse shape quality, therefore it has higher priority to be37Chapter 3. Frontal-Delaunay triangulationrefined. After several insertion, we come to a stage where there are only a few triangles in thequeue, as shown in Fig. 3.7c. After the last insertion, as shown in Fig. 3.7d, all the trianglesin the mesh satisfy both the size constraint and shape constraint. The queue is empty and thescheme terminates. There are still several triangles that are far from equilateral. Increasing thequality bound will help (Fig, 3.7e). However, it comes with the risk of generating unnecessarysmall cells or even termination problems, as shown in Fig, 3.7f.38Chapter 3. Frontal-Delaunay triangulation(a) Constrained Delaunay triangulation (b) After the first insertion(c) The mesh just before final insertion (d) Final mesh(e) Increase the quality bound (f) Generating unnecessary small trianglesFigure 3.7: Illustration of Engwirda’s frontal-Delauany method39Chapter 4Advancing-front local reconnectionmethodMarcum’s advancing-front local reconnection (AFLR) method is also proven to generate highquality meshes [22]. In this chapter, a similar annulus domain as the one in Marcum’s paperis used as an example to illustrate this AFLR scheme. The inner circle is radius 1 while theouter circle is radius 4. The circles are divided into 16 edges, as shown in Fig. 4.1a. First, thelocal feature sizes (lfs) are calculated. To simplify the illustration, all the vertices on the innercircle have the same lfs 2pi/16≈ 0.4. The lfs on the outer circle is 8pi/16≈ 1.6. The gradientparameter g is set to 2.5 to make the lfs equation (Eqn. 2.1) valid. Set the center of both circlesat (0,0), then the lfs in this domain can be simplified aslfs(r,θ) = 0.4+ r−12.5 , 1≤ r ≤ 4 .In Marcum’s scheme, the active/off flags are applied to all the triangles to indicate whetherthey are good or bad in a size criterion. If the length of each edge on a triangle is less than the1.5 times average length scale of the two end vertices of that edge, the triangle is regarded asgood triangle and applied the off flag. The off flag means that this triangle is turned off andcurrently not to be refined. Otherwise, the active flag is applied. As shown in Fig, 4.1a, allthe boundary edges are exactly equal to the average length scales. As for the edges connectingtwo boundaries, the allowable lengths are 1.5× (1.6+0.4)2 = 1.5. Since the radius differenceof the two circles is 3, all these edges fail the size quality test. Therefore all the triangles inconstrained Delaunay triangulation are marked as active. The insertion points are calculated40Chapter 4. Advancing-front local reconnection method(a) Starting constrained Delaunay triangulation (b) The points to be inserted(c) Remove all triangles and reconnect (d) After local reconnection(e) Insertion vertices (f) Final meshFigure 4.1: Marcum’s original scheme41Chapter 4. Advancing-front local reconnection methodfor all the active triangles. In this scheme, the new point will form a equilateral triangle withthe shortest edge of the active triangle. In Fig. 4.1b, all the new vertices are shown in bluedots. Then, the triangles who contains the new vertices need to be removed and reconnectedwith the new vertices, which means every triangle for this case. Fig. 4.1c shows the mesh afterthe first insertion. A local reconnection scheme is applied to maximize the minimum angles inthe mesh, as shown in Fig. 4.1d. The active/off flags are applied to the new triangles. All theactive triangles are shadowed in Fig. 4.1d. All the vertices to be inserted are calculated onceagain. In Fig. 4.1e, the green dots are generated from edges closer to the center while red dotsare from the outer circle. Some of these are too close to each other. Therefore, a rejection testis applied. In Marcum’s scheme, if the distance between two vertices is smaller than 0.7 timesthe larger lfs between the two, the vertex with larger lfs is rejected. In Fig. 4.1e, all the greenvertices are rejected from the mesh. After the local reconnection, the final mesh is shown inFig. 4.1f.Because of the difficulty of implementing Marcum’s method into GRUMMP’s framework,a few necessary modifications are made. It gives comparable results as Marcum’s algorithm,as shown in Fig. 4.2 . Instead of only considering the size constraint in Marcum’s scheme, Ialso take the shape quality (Eqn. 2.2) into consideration. A triangle is set to be off, or in otherwords, marked as a good triangle only if it satisfies both constraints. To be consistent with theprevious two algorithms, the modified algorithm still uses Watson’s insertion to improve thelocal topology. The front and active/off flags, or quality tests, are updated after each insertion,rather than updating after all the points are inserted into the active elements as introduced inMarcum’s scheme. A new denotation, layer index, is introduced to make sure the refinement ismarching layer-wise as in Marcum’s method. The layer index indicates which layer the edge ison. Starting with the constrained Delaunay mesh, the inner boundary edges are marked as layerindex 0 and all other edges are marked as layer index ∞, as shown in Fig. 4.3a. When usingWatson insertion, an empty hull (Fig. 4.3b) is obtained. The minimum layer index is foundamong all the edges on the hull. In this circumstance, it is 0. After the insertion, the layer42Chapter 4. Advancing-front local reconnection methodFigure 4.2: AFLR meshindex of all newly created edges is the minimum layer index +1. As shown in Fig. 4.3c, thenewly inserted edges are all marked as layer index 1. The local reconnection is omitted in ourmethod because we are using Watson’s insertion, which guarantees a local max-min (maximizethe minimum angles) criterion. The elements in the queue are edges instead of triangles as inthe other schemes. The queue priority scheme will be discussed later.Another big difference from Marcum’s scheme is the calculator for vertices to be inserted.In Fig. 4.4, c, c′, and c′′ are the circumcentre, the shape-optimal vertex and the size-optimalLayer Index0(a) Layer index before the inser-tionLayer Index0(b) Empty hullLayer Index01(c) Layer index after the insertionFigure 4.3: Layer index43Chapter 4. Advancing-front local reconnection methodv1v2cemin c''c'ceFigure 4.4: Four potential verticesvertex respectively, as discussed in Chapter 3. ce denotes the vertex generating an equilateraltriangle with shortest edge emin. The vertex nearest to emin among the four vertices is chosenfor insertion.The brief procedure of the modified AFLR scheme is illustrated by the same annulus ex-ample. In this case, the inner circle is set as layer index 0 at the beginning. Other edges in themesh are set as layer index ∞. The triangles in Fig. 4.5 are colored by the layer index. Sincethe layer index is only applied to edges, the cell layer index is introduced as:Cell Layer Index =3maxi=1(LayerIndex(ei))+0.1×3mini=1(LayerIndex(ei)) . (4.1)This is shown in Fig. 4.5b. After the first insertion, the newly generated triangle has two newedges with layer index 1 and a layer index 0 boundary edge. Therefore, its cell layer index is1+0.1×0 = 1.The general procedure is outlined as follows:44Chapter 4. Advancing-front local reconnection method(a) Constrained Delaunay triangulation (b) The first insertion(c) First layer finished (d) Final meshFigure 4.5: Illustration of modified AFLR scheme45Chapter 4. Advancing-front local reconnection method1. Start with the constrained Delaunay triangulation (Fig. 4.5a).2. Initialize the data structure. For each triangle, if its each edge length is smaller than 1.5times the averaged lengthscale and its radius-edge ratio is larger than the threshold value,this triangle is treated as a good triangle. Otherwise, it is treated as bad.Example In the annulus case, all the triangles are bad. In Fig. 4.5a, the blue triangleshave minimum angles less than 29◦, which is the threshold value in this case. Be-sides that, all the triangles fail the size constraint test as discussed in Marcum’soriginal scheme.3. Queue boundary edges with bad adjacent triangles and interior edges which have a goodadjacent cell on one side and a bad adjacent cell on the other side.Example In Fig. 4.5a, only the boundary edges are queued because they all have badadjacent triangles. All the interior edges are shared by bad triangles, therefore theyare not queued. In Fig. 4.5b, after the first insertion, a good triangle is generated.The boundary edge on it is removed from the queue. Then the two newly generatedinterior edges are queued because now they have a good triangle on one side and abad one on the other side.4. The queue priority is calculated by the sum of layer index, edge length, and shape qualityof the adjacent bad triangle. The smaller the sum is, the higher priority the edge has.ExamplePriority = Layer Index+ ||e||+ ||emin||√3rr denotes the circumradius of the adjacent bad triangle. The layer index is an in-teger ranging in [0,∞]. The shape quality ||emin||√3rranges in [0,1] and ||e|| is usuallymuch smaller than 1 (not in this example). The layer index is dominant in priority.Therefore when comparing two similar edges, the one with a lower layer index has46Chapter 4. Advancing-front local reconnection methoda higher priority. This ensures the scheme works in a layer-wise marching pattern.When the layer index is the same, the shape quality is dominant. This means thatwhen comparing two edges on the same layer, the one adjacent to a worse trianglehas a higher priority, which is consistent to Engwirda’s frontal-Delaunay schemethat the worst shaped triangle is refined first. If these two values are exactly thesame, the shorter edge jumps ahead. The small edge or triangle may locate in thearea of interest, where we want regular triangles.5. Choose the first edge from the queue and find the bad triangle adjacent to it. Calculatefour kinds of points: the size-optimal vertex and shape-optimal vertex mentioned in theprevious algorithm, the point generating an exact equilateral triangle and the circumcen-ter. Choose the point that is nearest to the base edge.Example This is illustrated in Fig. 4.4.6. Reject the candidate point if the insertion generates an edge shorter than 0.7 times thebase edge length.Example This rejection is inherited from Marcum’s original scheme, preventing unnec-essary small triangles being generated.7. Use Watson insertion to insert the new point. The new edges created by the insertion areadded to the queue if they satisfy the criteria in step 3.8. Remove the edge from the queue regardless of whether the point is rejected or accepted.9. Repeat steps 5-8 until the queue is empty.Example Take the same annulus case, thanks to the layer index dominant priority, allthe boundary edges on the inner circle has the highest priority. The first layer isgenerated based on boundaries as shown in yellow segments in Fig. 4.5c. In this47Chapter 4. Advancing-front local reconnection methodcircumstance, these edges form a slightly larger circle, or front. They are queuedbecause sharing by good triangles and bad ones. Those edges between the innercircle and these yellow fronts are removed from the queue because they have goodtriangles on both sides. After layer-wise insertions, the scheme terminates underconditions like those in Fig. 4.5d. All triangles in the mesh are good ones, whichalso means all the edges are adjacent to good triangles. In other words, no moreedges are queued.48Chapter 5Numerical resultsThe mesh regularity, truncation error, discretization error, and time costs for several numericaltest cases are discussed in this chapter to answer the questions posed in Section Poisson equation on a square domain5.1.1 Mesh regularityThree different mesh refinement schemes are applied to a square domain, of which the edgelength is 1. Starting from the same constrained Delaunay triangulation, shown in Fig. 5.1a, thethree schemes generate different meshes. Delauany refinement, as shown in Fig. 5.1b, doesnot generate any obvious chunks of equilateral triangles. In Fig. 5.1c, Engwirda’s algorithmuses a different priority queue and different insertion points from Delaunay refinement. Onlythe triangles whose shortest edges satisfy the local length scales can be added to the queue.This ensures the newly created triangles are based on the existing good triangles. As shownin Fig. 5.2, Engwirda’s scheme introduces a “snake” pattern. It turns into the area where theworst triangles exist. It generates several large chunks of equilateral triangles in the middle.However, badly shaped triangles emerge when the scheme tries to fill the gaps between thesechunks because they are not perfectly aligned. Marcum’s AFLR scheme refines the mesh ina different way. It starts from the boundary and marches into the interior layer by layer. Asshown in Fig. 5.3, each layer consists of nearly equilateral triangles along the boundaries anda few badly shaped ones at the corners. In the end, the mesh is almost perfect except the crossshown in Fig. 5.1d.495.1. Poisson equation on a square domain(a) Starting Delaunay Triangulation (b) Delaunay Refinement(c) Engwirda frontal-Delaunay (d) Marcum’s AFLRFigure 5.1: Meshes generated on a square domain505.1. Poisson equation on a square domain(a) After 100 insertions (b) After 500 insertions(c) After 1000 insertions (d) After 1500 insertionsFigure 5.2: The process of Engwirda’s scheme515.1. Poisson equation on a square domain(a) After the first layer (b) After the fifth layerFigure 5.3: The process of AFLR schemeDelaunay refinement uses circumcenters as Steiner points. The badly sized triangles arerefined first. Inserting a circumcenter into a badly sized triangle does not guarantee the newtriangle satisfies the shape constraint. Therefore, a bad triangle may be refined several times tosatisfy the size and shape criterions. Because both the AFLR scheme and Engwirda’s frontal-Delaunay scheme select Steiner points more cleverly, the newly created triangle will satisfythe size and shape constraints. Generally speaking, the sizes of meshes obtained from AFLRscheme and frontal-Delaunay scheme are smaller than those from Delaunay refinement underthe same constraints. They both improve the regularity of the meshes, but which scheme gen-erates the most equilateral triangles for this case? As illustrated in Fig. 5.4, the AFLR meshhas about 70 percent of all the angles between 55◦ and 65◦, while Engwirda’s frontal-Delaunaymesh is only 5 percent behind. The traditional Delauany method generates the worst mesh forthis regularity comparison.5.1.2 Truncation error and discretization errorThe first test case for the relationship between errors and mesh regularity is solving the Poissonequation,∂ 2T∂x2+∂ 2T∂y2= S, (5.1)525.1. Poisson equation on a square domain20 30 40 50 60 70 80 90 100 110 120051015202530354045Angle(Degree)Percentage  DelaunayEngwirdaAFLRFigure 5.4: Angle distribution for meshes on square domainon a square domain. Substitute the manufactured solution shown in Fig. 5.5,T = pisin(pix)sin(piy)8,into the Poisson equation 5.1, then the source term S can be obtained easily,S =−pi3 sin(pix)sin(piy)4.A second order finite-volume discretization and least-squares reconstruction are used[36,15]. Setting the control volume averages of the manufactured solution as the initial condition,the flux integral at the first step,1ACV[−ˆCVSdA+˛∂CV(∇T ) ·~nds],535.1. Poisson equation on a square domainFigure 5.5: Sinusoidal manufactured solution of the Poisson equationis the truncation error. The converged numerical solution is used to compute the discretizationerror with the control volume averages of the manufactured solution. In Fig. 5.6, we canobserve that in the areas consisting of nearly equilateral triangles, the magnitudes of truncationerror and discretization error tend to be very near zero, relative to the error where the mesh isless regular. Even in the Delaunay mesh, where most triangles are far from equilateral ones,we can still find a few well shaped control volumes which have very small error magnitudes.In the meshes obtained from Engwirda’s frontal-Delaunay and Marcum’s AFLR schemes, thispattern is much more obvious. The areas marked as red in the left-most figures in Fig. 5.6indicate chunks of nearly equilateral triangles. In the same areas, we can expect nearly zerotruncation errors. However, the discretization error seems more sensitive to the mesh regularity.As we can see, nearly zero discretization error only occurs at the perfect equilateral triangles.Therefore when we compare the L2 norm of the errors, shown in Fig. 5.7, a large decreasein the magnitude of the truncation error can be found in more regular meshes, along with anincrease in the order of accuracy. For the discretization error, on the other hand, the order ofaccuracy for all three kinds of meshes tends to be the same. A small magnitude decrease canstill be observed in the two more regular meshes.545.1. Poisson equation on a square domain(a) Delaunay Mesh Regularity (b) Delaunay Truncation Error (c) Delaunay Discretization Error(d) Frontal-Delaunay Mesh Regularity (e) Frontal-Delaunay Truncation Error (f) Frontal-Delaunay Discretization Error(g) AFLR Mesh Regularity (h) AFLR Truncation Error (i) AFLR Discretization ErrorFigure 5.6: Mesh regularity, truncation error, and discretization error (second order)555.1. Poisson equation on a square domain102 103 104 10510−1100Number of Control VolumesL2 norm of Truncation Error  DelaunayEngwirdaAFLR0.5th Order(a) Truncation Error102 103 104 10510−610−510−410−310−2Number of Control VolumesL2 norm of Discretization Error  DelaunayEngwirdaAFLR2nd Order(b) Discretization ErrorFigure 5.7: Order of Accuracy for Second Order Discretization10 2 10 3 10 4 10 5Number of Control Volumes10 -510 -410 -310 -210 -1L2 Norm of Truncation ErrorDelaunayEngwirdaAFLR2.5th Order(a) Truncation Error10 2 10 3 10 4 10 5Number of Control Volumes10 -710 -610 -510 -410 -3L2 Norm of Discretization ErrorDelaunayEngwirdaAFLR3rd Order(b) Discretization ErrorFigure 5.8: Order of Accuracy for the Fourth Order DiscretizationThe fourth order discretization scheme [17, 29] is employed to show whether high orderreconstruction and flux integral will affect the relationship between errors and the mesh regu-larity. Fig. 5.8 shows the convergence when performing the fourth order discretization. For thetruncation error, the advantage of regular meshes for the order of accuracy disappears. They allconverge at 2.5th order. But a regular mesh still gives lower magnitude of the error. When itcomes to the discretization error, the AFLR mesh, which has most equilateral triangles, reducesthe error magnitude dramatically, and increases the order of accuracy. The size of Marcum’sAFLR mesh can be 1/16 of the Delaunay mesh while the error magnitude is about the same.565.1. Poisson equation on a square domainMethod Number of Vertices in mesh Total Time(s) Time per insertion(s)Delaunay 1915 0.957 4.997×10−4Engwirda 1597 0.965 6.043×10−4AFLR 1188 1.507 1.269×10−3Table 5.1: Mesh Generation Time Cost5.1.3 Time costThe two mesh generation schemes which generate more regular meshes take more time for eachinsertion because they calculate multiple possible Steiner points, and then reject some of them.However, they may require fewer cells to fill up the domain than the Delaunay triangulation.In Table 5.1, we see that Engwirda’s frontal-Delaunay mesh generation scheme only takes alittle longer than Delaunay triangulation in each insertion. With the help of fewer vertices inthe mesh, the time cost in mesh generation is nearly the same. Marcum’s AFLR scheme is thebest scheme in terms of regularity, but checking whether a candidate vertex is too close to theexisting ones takes a lot of time for each insertion. Even though it has the smallest mesh size,the time cost for AFLR is around twice that of the other two schemes. The difference in thetime cost between the two new schemes and Delauany triangulation grows larger as the size ofmesh increases, as shown in Fig. 5.9.Considering that the regularity of mesh can reduce the error magnitude, the two newschemes have smaller time costs for a given error when combining both the mesh generationcost and the solver cost together. Fig. 5.10 shows the combined time cost. For the second orderdiscretization, it is not worth waiting for the most regular AFLR mesh. In that case, Engwirda’sfrontal-Delaunay scheme is the optimal one, in that it reaches a certain error magnitude fastest.For the fourth order discretization, the advantage of Marcum’s AFLR scheme is obvious. Eventhough it takes much longer time in mesh generation, the smaller mesh size saves a lot of timein the solver and produces much lower error for a given target edge length in the mesh. Thispays off the large difference in mesh generation time cost when the mesh gets larger. In prac-tice, usually the mesh generation code is run once and the mesh can be solved several times. In575.1. Poisson equation on a square domain0 1000 2000 3000 4000 5000 6000 700001234567Number of vertices in meshMesh generation time (s)  DelaunayEngwirdaAFLRFigure 5.9: Mesh generation time on a square domain10−1 100 10110−610−510−410−310−2Wall clock time in secondsL2 Norm of Discretization Error  DelaunayEngwirdaAFLR(a) The 2nd Order10−1 100 101 10210−710−610−510−410−3Wall clock time in secondsL2 Norm of Discretization Error  DelaunayEngwirdaAFLR(b) The 4th OrderFigure 5.10: Time cost of the mesh generation and the solver585.2. Advection-diffusion equation in a channelthat case, the regular mesh is definitely worth waiting for.5.2 Advection-diffusion equation in a channelThe second test case is solving the advection-diffusion equation in a channel.5.2.1 Mesh regularityThe geometry is a rectangular channel with height H = 1 and length L = 3. As shown inFig. 5.11, similar results as meshes on a square domain can be observed. The mesh fromAFLR scheme has a strikingly large area of equilateral triangles. The only exceptions are theinlet, outlet, and a few cells in the middle. In this case, only the upper and lower boundaries areset to be layer index 0 at the beginning. The inlet and outlet have∞ layer index. The refinementis marched from upper and lower boundaries into the middle. Similar to the square case, theright angles make this scheme fail to generate perfect layers. This is why the inlet and outlethave bad cells adjacent to them. Engwirda’s frontal-Delaunay scheme does not improve thetriangles near the boundaries. However it has several large chunks of perfect triangles in theinterior domain. The Delaunay refinement still does not guarantee any equilateral triangles.The angle distribution shown in Fig. 5.12 also confirms the mesh qualities are similar to theprevious square domain case.The AFLR scheme has its limit. The finest mesh obtained from Marcum’s AFLR schemeis discussed here to show where the AFLR scheme fails. The equilateral triangles layers havethe limit. In Fig. 5.13, the triangles are colored by the cell layer index (Eqn. 4.1). In thisfigure, the maximum layer index is 59.3. There are about 50 layers in the mesh. The first20 layers are almost perfect. However, at the transient layer where the perfect pattern breaks,there are several triangles that barely passed the local quality test. They make the edge lengthsin the next layer vary a lot and result in the perfect layers terminating. If the minimumangle allowed in the mesh is set higher than 28◦, the hatched triangle shown in Fig. 5.14a is595.2. Advection-diffusion equation in a channel(a) Delaunay triangulation (1652 vertices)(b) Engwirda’s frontal-Delaunay (1401 vertices)(c) Marcum’s AFLR (1689 vertices)Figure 5.11: Meshes in a channel605.2. Advection-diffusion equation in a channel30 40 50 60 70 80 90 100 110051015202530354045AnglePercentage  DelaunayEngwirdaAFLRFigure 5.12: Angle distribution of meshes in a channelFigure 5.13: Bad AFLR mesh615.2. Advection-diffusion equation in a channel(a) The normal pattern(b) AFLR fails circumstanceFigure 5.14: AFLR scheme failureFigure 5.15: Near view of the transition area625.2. Advection-diffusion equation in a channelrefined. This generates another three almost equilateral triangles, keeping this layer completeand regular. However, if angles less than or equal to 28◦ are allowed in the mesh, the triangleto be refined would be the hatched one in Fig. 5.14b. This leads to two adjacent edges withvery different edge lengths. The triangles generated are also irregular. Because there is nosmoothing technique performed in the middle of refinement, in the next layer, the trianglesaround this area will be worse, as shown in Fig. 5.15. The point P is generated in the patternfrom Fig. 5.14b. Slight misallignment is made on the layer and after two more layers, theperfect pattern blows up. One obvious amendment for this case is to increase the quality boundto make each layer perfect. But not all geometries can be meshed perfectly. Increasing thequality bound is highly likely to introduce problems with robustness.5.2.2 Truncation error and discretization errorThe 2D advection-diffusion equation~u ·∇φ = α∇2φis solved in this case. A manufactured exact solution (shown in Fig. 5.16) is provided:φ =(eR1xR2eR2L− eR2xR1eR1LR2eR2L−R1eR1L)sin(piy),where R1 = H2α +√H24α2 +pi2 and R2 = H2α −√H24α2 +pi2. The coefficient of diffusion α is setas 0.01. The boundary conditions for both upper and lower boundaries are enforced by φ = 0,whereas at the outlet, the gradient is set as zero. For the inlet boundary, the boundary conditionis set as the exact solution φ = sin(piy). Notice that φ changes in the y-direction more thanthe x-direction. The AFLR mesh is more regular than Engwirda’s frontal-Delaunay mesh inthe y-direction, so we expect the truncation error on the AFLR mesh will be much less thanEngwirda’s mesh. Both more regular meshes can be expected to give smaller truncation error635.2. Advection-diffusion equation in a channelFigure 5.16: Exact solution for advection-diffusion equation102 103 104 105Number of Control Volumes10-410-310-210-1L 2 norm of Truncation errorDelaunayEngwirdaAFLR1st Order2nd Order(a) Truncation error102 103 104 105Number of Control Volumes10-510-410-310-210-1L 2 norm of Discretization errorDelaunayEngwirdaAFLR2.5th Order(b) Discretization errorFigure 5.17: Order of Accuracy for second order discretizationthan traditional Delaunay refinement.The convergence of truncation errors is shown in Fig. 5.17a. The Delaunay mesh convergesat first order, while the other two more regular meshes can achieve second-order for coarsemeshes. As the mesh becomes finer, they become about first order. The large truncation erroroccurs at the area where irregular triangles exist. When we refine the mesh, the truncation errorfor the perfect control volumes decreases at second order whereas the bad control volumes areonly about 0.25th order. On the very coarse mesh, the magnitudes of truncation error for regularand irregular control volumes are comparable, so they tend to get a higher order of accuracyglobally. However, when the mesh is refined, the truncation errors for badly shaped control645.2. Advection-diffusion equation in a channel102 103 104 105Number of Control Volumes10-610-510-410-3L 2 norm of Truncation errorDelaunayEngwirdaAFLR3rd Order(a) Truncation error102 103 104 105Number of Control Volumes10-810-710-610-510-410-3L 2 norm of Disretization errorDelaunayEngwirdaAFLR4th Order(b) Discretization errorFigure 5.18: Order of Accuracy for fourth order discretizationvolumes can be hundreds of times larger than the global average. This large error magnitudeis dominant, which lowers the global order of accuracy. For the discretization error, shown inFig. 5.17b, the improvement is not obvious. Similar to the previous Poisson equation test case,the relation between truncation error and discretization error is complicated. Small truncationerror does not guarantee small discretization error.Fourth order discretization is also applied for this case. Again, a decrease in the magnitudeof truncation error can be observed but for the discretization error, the improvement is small.5.2.3 Time costSince the discretization error is not improved as we expected, the time cost is compared tosee if it is worthwhile to generate regular meshes. In Table 5.2, the mesh generation timecost is comparable to the second order solving time cost. The AFLR mesh has very longmesh generation time and it is not worth waiting since it does not decrease time cost in solverand discretization errors. It takes nearly the same time to generate the Delaunay mesh andEngwirda’s frontal-Delaunay mesh. With the help of the smaller frontal-Delaunay mesh size,it takes a shorter time for Engwirda’s frontal-Delaunay mesh to converge. The decrease in timecost is not large but it does no harm to generate Engwirda’s frontal-Delaunay mesh. As for thefourth order solver shown in Table 5.3, the time cost decrease for Engwirda’s mesh is larger.655.3. Potential flow around a circular cylinderMesh Size Delaunay Engwirda AFLRMesh Solver Mesh Solver Mesh Solver~400 0.28 0.29 0.28 0.31 0.55 0.31~1600 0.75 0.81 0.71 0.63 2.03 0.87~6400 2.65 3.01 2.60 2.42 8.52 3.16~25000 10.28 13.59 10.86 11.20 33.00 14.68Table 5.2: Time cost in mesh generation and second order solverMesh Size Delaunay Engwirda AFLRMesh Solver Mesh Solver Mesh Solver~400 0.28 0.60 0.28 0.59 0.55 0.82~1600 0.75 2.60 0.71 1.89 2.03 2.54~6400 2.65 9.52 2.60 8.76 8.52 10.98~25000 10.28 43.02 10.86 32.18 33.00 39.79Table 5.3: Time cost in mesh generation and fourth order solverConsidering the slightly better error performance, Engwirda’s mesh is definitely recommendedfor high order solver in this test case.5.3 Potential flow around a circular cylinderThe incompressible flow is solved in ANSLib’s compressible flow solver with some modifica-tions.5.3.1 Mesh regularityThe domain is bounded by two circles. The inner circle denotes a circular cylinder whoseradius is 1. The wall boundary condition is set for the circular cylinder. The outer circledenotes the far field boundary. The radius is 50. Three different mesh generation schemesare applied to generate meshes in which cell sizes increase gradually with the distance fromthe circular cylinder. Fig. 5.19 shows that in the mesh obtained from the AFLR scheme,most triangles have minimum angles about 50◦. The circle is a good geometry for a layer-wisemarching method. Even though they are not nearly equilateral triangles, they grow in a smooth665.3. Potential flow around a circular cylinder(a) Delaunay mesh(b) Engwirda’s frontal-Delauany(c) Marcum’s AFLRFigure 5.19: Meshes around a circular cylinder675.3. Potential flow around a circular cylinder30 40 50 60 70 80 90 100 11005101520253035AnglePercentage  DelaunayEngwirdaAFLRFigure 5.20: Angle distribution of meshes around a circular cylindervariation. The AFLR scheme usually fails to generate beautiful triangles everywhere becausebad triangles are kept in the mesh as discussed in the previous section. In this case it onlyhappens in the far field. It will not affect the overall error improvement since the flow varieslittle in that area. Engwirda’s frontal-Delaunay scheme also generates many triangles whoseminimum angles are above 50◦. The Delaunay mesh still is the worst in terms of the regularity,as shown in Fig. Truncation error and discretization errorIn this case, the Euler equations are solved∂∂ tρρuρvE+∂∂xρuρu2+Pρuvu(E +P)+∂∂yρvρuvρv2+Pv(E +P)= 0, (5.2)685.3. Potential flow around a circular cylinderwhere the energy E = P/(γ − 1) + ρ(u2 + v2)/2. The analytic solution is available for theincompressible potential flow around a circular cylinder. We consider the case where the centerof the cylinder is (0,0). There is no circulation and the far field flow is in the x-direction. Thecomplex conjugate velocity W = u− iv is :W =V∞− V∞R2z2where z = x+ iy and V∞ is the free stream velocity. In our in-house code, the equations aresolved in non-dimensional form, so the non-dimensional Mach number M∞ = 0.3 is used asV∞. The radius R = 1. For this test case, the incompressible solution can be simplified toρ = 1W = M∞−M∞ · 1(x+ iy)2,u = real(W ) and v =−imag(W ). The pressure is calculated by Bernoulli’s equationp = p∗− 12ρ(u2+ v2),where p∗ = 1γ ·(1+ γ−12 M2∞) γγ−1 .The analytic solution shown in Fig. 5.21 is for an incompressible solver. To use this asa manufactured solution for compressible Euler equation 5.2, the source term calculated bysubstituting the analytic solution into Eqn. 5.2 needs to be added.∂∂ tρρuρvE+∂∂xρuρu2+Pρuvu(E +P)+∂∂yρvρuvρv2+Pv(E +P)=000S,695.3. Potential flow around a circular cylinderFigure 5.21: Exact solution for the non-dimensional pressure around a circular cylinderOnly the energy equation has a source term:S =∂∂x(u ·(γ pγ−1 +u2+ v22))+∂∂y(v ·(γ pγ−1 +u2+ v22)).With the help of the complex conjugate velocity W , the source term can be evaluated in S(x,y).The truncation error is estimated by the flux integral when setting the manufactured solu-tion as the initial solution. As shown in Fig. 5.22a, the magnitude of truncation error in theAFLR mesh is reduced compared to the other two meshes. The Engwirda’s frontal-Delaunaymesh also decreases the truncation error magnitude. As for the discretization error, shown inFig. 5.23a, Marcum’s AFLR method improves the magnitude of the error. All the three meshesachieve second-order accuracy on discretization error, which is consistent to the discretizationscheme in use. When the fourth order discretization scheme is applied, similar results can beobserved in Fig. 5.22b for the truncation errors. The discretization error shown in Fig. 5.23billustrate that the AFLR mesh improves both error magnitude and the order of accuracy. Inthe fourth order solver, the finest Delaunay mesh and Engwirda’s frontal-Delaunay mesh failto converge. This is because that we are solving the incompressible flow in a compressible705.3. Potential flow around a circular cylinder102 103 104 105Number of Control Volumes10-510-410-310-2L 2 norm of truncation errorDelaunayEngwirdaAFLR1.5th order(a) The second order solver102 103 104 105Number of Control Volumes10-510-410-310-2L 2 norm of truncation errorDelaunayEngwirdaAFLR2nd order(b) The fourth order solverFigure 5.22: Truncation error for the energy equation102 103 104 105Number of Control Volumes10-610-510-410-310-2L 2 norm of discretization errorDelaunayEngwirdaAFLR2nd order(a) The second order solver102 103 104Number of Control Volumes10-610-510-410-310-2L 2 norm of discretization errorDelaunayEngwirdaAFLR3rd order1st order(b) The fourth order solverFigure 5.23: Discretization error for the pressure715.3. Potential flow around a circular cylinder100 101 102Time cost in seconds10-510-410-310-2L 2 norm of discretization errorDelaunayEngwirdaAFLRFigure 5.24: Time cost for both mesh generator and solver (fourth order)flow solver. The solver tries to work against the source term. The lack of dissipation in thehigh order solver makes it more sensible to the errors. The AFLR mesh is the best in errormagnitude but the time cost must be checked since the Euler equations are harder to solve thanthe Poisson equation and the advection-diffusion equation.5.3.3 Time costThe advantage of the AFLR scheme in the mesh regularity also contributes to time savings. Asshown in Fig. 5.24, AFLR is definitely the best mesh generation scheme for this potential flowaround a cylinder case as measured by combined time to a given solution error. Note that themesh generation time cost for AFLR scheme is about twice as high as the other two (Table 5.4).But the smaller size and the mesh regularity help to reduce the time cost in the high order solver.The number of non-linear iterations required by the AFLR mesh is the smallest, which containsthe time demanding tasks such as the Jacobian calculation and factorization.For solving the potential flow around a circular cylinder, Marcum’s AFLR scheme is obvi-725.4. NACA-0012 airfoilMesh Size Delaunay-refinement Engwirda’s frontal-Delaunay Marcum’s AFLR~300 0.2687 0.2433 0.3647~1,000 0.7430 0.7026 1.0411~4,000 2.8018 2.9261 4.6033~20,000 12.5435 11.5693 20.3652Table 5.4: Time cost in mesh generation (seconds)Delaunay Engwirda AFLRMesh size 6405 5675 3840Non-linear iteration (NLI) 9 10 6Linear iteration (LI) 78 234 37Time cost (s) 65.90 67.17 26.24Table 5.5: Fourth order solver time costous the best choice in both reducing errors and saving time.5.4 NACA-0012 airfoilIn this test case, the compressible inviscid flow around the NACA-0012 airfoil is solved. TheMach number is 0.5 and different angles of attack (α) are applied.5.4.1 Mesh regularitySimilar to the cylinder case, Engwirda’s frontal-Delaunay scheme and Marcum’s AFLR schemeimprove the mesh regularity in different ways. As shown in Fig. 5.25b, in Engwirda’s mesh,large areas of dark red are always expected but at the same time, they lack continuity. A perfectequilateral triangle is not necessarily adjacent to another equilateral triangle or nearly equilat-eral triangle. On the other hand, in the mesh obtained from Marcum’s AFLR scheme, thefirst few layers around the boundary usually contain lots of equilateral triangles. As shown inFig. 5.25c, the areas around the leading edge and trailing edge are filled with nearly equilateraltriangles. They also show a continuous pattern thanks to the layer-wise marching pattern. Inthis case, the triangles around the airfoil in the AFLR meshes are not all nearly equilateral735.4. NACA-0012 airfoiltriangles. Several of them have the minimum angle of 45◦. As shown in Fig. 5.26, the AFLRscheme generates fewer nearly 60◦ angles than Engwirda’s frontal-Delaunay scheme. But bothof them are still quite good when compared with the Delaunay mesh.5.4.2 ErrorsThe exact solution of this case is not available, therefore the truncation error is omitted forthis case. In Fig. 5.27, the lift coefficient and drag coefficient are estimated when the angle ofattack is 2◦. Fig. 5.27a shows the convergence of the lift coefficient. The solid lines are for thesecond order discretization and dashed lines denote the fourth order discretization. For the sec-ond order discretization, the AFLR mesh, which has the best mesh quality around the airfoil,converges slightly faster than Engwirda’s mesh. The Delaunay mesh has a different conver-gence pattern, with a slightly different grid-converged value. For higher order discretization,the AFLR mesh and Engwirda’s frontal-Delaunay mesh perform consistently with second or-der scheme. As for the Delaunay mesh, the convergence pattern is similar to the other twomeshes. But it still gets a different value. As shown in Fig. 5.27b, the convergence of the dragcoefficient is almost exactly the same for all three meshes. Generally speaking, mesh regular-ity does not lead to significant improvement on lift coefficients and drag coefficients. For thesecond order scheme, mesh regularity helps to converge faster to an acceptable solution. Butwhen it comes to higher order, the improvement is not significant.The numerical tests are re-run for the same meshes when setting the angle of attack as 4◦.Similar results can be observed in Fig. 5.28. Again, mesh regularity does not help in this case.745.4. NACA-0012 airfoil(a) Delaunay mesh(b) Engwirda mesh(c) AFLR meshFigure 5.25: Regularity for meshes around NACA-0012 airfoil755.4. NACA-0012 airfoil30 40 50 60 70 80 90 100 11005101520253035AnglePercentage  DelaunayEngwirdaAFLRFigure 5.26: Angle distribution of meshes around NACA-0012 airfoil102 103 104 105Number of Control Volume0.2520.2540.2560.2580.260.2620.2640.2660.2680.27Lift coefficientDelaunayEngwirdaAFLR2nd Order4th Order(a) Lift coefficient102 103 104 105Number of Control Volume012345678Drag coefficient10-3DelaunayEngwirdaAFLR2nd Order4th Order(b) Drag coefficientFigure 5.27: NACA-0012 airfoil with angle of attack 2◦765.4. NACA-0012 airfoil102 103 104 105Number of Control Volume0.5050.510.5150.520.5250.530.5350.54Lift coefficientDelaunayEngwirdaAFLR2nd Order4th Order(a) Lift coefficient102 103 104 105Number of Control Volume23456789101112Drag coefficient10-3DelaunayEngwirdaAFLR2nd Order4th Order(b) Drag coefficientFigure 5.28: NACA-0012 airfoil with angle of attack 4◦5.4.3 Time costSince mesh regularity does not help to get an accurate solution on a coarser mesh as expected,it is doubtful whether it reduces the total time cost. As shown in Fig. 5.29, the improvementsin time cost for both second order and fourth order scheme are not obvious. Table 5.6 showsthe time cost for the total of mesh generation and solving with 2◦ angle of attack and thefourth order scheme. We can observe that the time cost in the solver is much higher than thatin the mesh generator. Therefore in this case, the additional time cost of Marcum’s AFLRscheme can be neglected. The advantage of the mesh regularity shortening the solving processis only observed for fine meshes. The time cost of Engwirda’s mesh is about half of the oneof Delaunay mesh. Marcum’s mesh also decreases time cost by about 25 percent. Whenrunning the cases several times, more regular meshes are recommended to reduce the timeMesh Size Delaunay Engwirda AFLRMesh Solver Mesh Solver Mesh Solver~1000 0.48 5.85 0.55 5.49 0.52 5.64~4000 1.77 25.54 1.85 22.40 4.01 23.83~16000 8.42 146.52 8.41 126.68 22.60 122.24~64000 30.66 1251.67 27.14 685.97 79.90 919.58Table 5.6: Comparison of time cost in mesh generation and the fourth order solver (α = 2◦)775.4. NACA-0012 airfoil100 101 102 103 104Time cost in secons0.2520.2540.2560.2580.260.2620.2640.2660.2680.27Lift coefficientDelaunayEngwirdaAFLR2nd Order4th Order(a) Time cost vs lift coefficient100 101 102 103 104Time cost in seconds012345678Drag coefficient10-3DelaunayEngwirdaAFLR2nd Order4th Order(b) Time cost vs drag coefficientFigure 5.29: Time cost of mesh generation and solver in totalDelaunay Engwirda AFLRMesh Size 54341 47613 51683NLI 22 15 16LI 156 134 177Time Cost 1251.67 686.97 919.58(a) AoA 2◦Delaunay Engwirda AFLRMesh Size 54341 47613 51683NLI 24 16 13LI 198 140 176Time Cost 1245.98 733.67 671.33(b) AoA 4◦Table 5.7: Time cost in the fourth order solver785.4. NACA-0012 airfoil~1000 ~4000 ~16000 ~64000Mesh Size0200400600800100012001400Time cost in fourth order solverDelaunayEngwirdaAFLRFigure 5.30: Time cost in the fourth order solver for AoA 4◦cost. For 4◦ angle of attack, similar results (Fig. 5.30) can be obtained: the AFLR mesh andEngwirda’s frontal-Delaunay mesh helps to decrease the time cost in the solver. Similar to thecircular cylinder case, mesh regularity helps reduce the mesh size and the number of non-lineariteration required. As shown in Table 5.7b, even though the AFLR mesh is a bit larger thanEngwirda’s frontal-Delaunay mesh, the smaller number of non-linear iterations required helpsto reduce the time cost. Again, the Delaunay mesh is worse than the other two in terms ofthe number of non-linear iterations. Notice that this improvement is only observed with a highorder discretization scheme. The decrease in time cost for the second order scheme is verylimited as shown in Fig. 5.31.795.4. NACA-0012 airfoil~1000 ~4000 ~16000 ~64000Mesh Size050100150200250Time cost in second order solverDelaunayEngwirdaAFLRFigure 5.31: Time cost in the second order solver for AoA 2◦80Chapter 6ConclusionsThe previous research on mesh regularity is completed on structured triangular meshes [8, 20].However, the impact of unstructured mesh regularity on errors is missing. In my thesis, Iincreased the unstructured mesh regularity with Engwirda’s frontal-Delaunay algorithm andMarcum’s AFLR algorithm. The impact of that on errors is a supplement to the previousresearch. Besides, the impact of mesh regularity on time costs in both mesh generation and theflow solver is also considered in this thesis, which is another missing part from the previousresearch. The unstructured mesh regularity reduces the time cost in the solver, especially in thehigh order solver, which is not expected before this study.Marcum’s AFLR algorithm works well in terms of the unstructured mesh regularity. Theoriginal scheme is not consistent with the Delaunay refinement framework, which brings dif-ficulties to implement this algorithm. In my thesis, I made several necessary modificationsto implement this algorithm in the Delaunay refinement framework. The layer index intro-duced in Chapter 4 works well to achieve the layer-wise pattern in the Delaunay refinementframework. The point selection method in Marcum’s algorithm is supplemented with the off-centres. The improved boundary discretization technique is introduced in Section 2.1. It isspecially designed for airfoils. It discretizes the boundary in a smooth way. This smoothnessof the boundary edge length helps the modified AFLR algorithm perform better in that theabrupt change of edge lengths will break the layer-wise pattern. With the help of these modi-fications, the modified AFLR algorithm works well and I can reproduce the similar results asthose from Marcum’s algorithm.After introducing the three different mesh generation schemes and several numerical test81Chapter 6. Conclusionscases, the questions posed in Section 1.1 can now be answered.1. How regular a triangular mesh can we generate for realistic cases?Depends on geometries. For simple geometries such as a square or annulus, almost70% of triangles in the most regular mesh we generated have minimum anglesabove 55◦. For more realistic cases like airfoils, the ratio decreases to 30%.2. What is the additional cost to generate such a mesh?Point calculation and rejection. In my implementation, both Engwirda’s frontal-Delaunayand Marcum’s AFLR scheme use similar vertex location calculators, which calcu-late two or three more possible vertex locations than Delaunay triangulation. Thesize-optimal vertices even require iterative bisection method to complete. Mar-cum’s AFLR method also adopts a rejection scheme. When it rejects a vertex tobe inserted, all the calculation time cost for this vertex is in vain. Both the pointlocation calculator and the rejection scheme increase the time costs at every step.The rejection scheme is dominant. Considering that the regular mesh size is usu-ally smaller than the Delaunay mesh, Engwirda’s frontal-Delaunay scheme, whichonly costs a little more because of the location calculator every step, can get similartime cost as Delauany mesh. The AFLR mesh costs much more than the other two,especially for very fine meshes.3. Considering that cell shape affects the accuracy greatly, can we expect small discretiza-tion error using unstructured regular meshes?Not really. In my numerical tests, regular meshes did not consistently produce lowererror.4. Can a coarse regular mesh offer the same accuracy as a finer irregular mesh with lowercosts in the flow solver?82Chapter 6. ConclusionsKind of. As answered in the previous question, the solution accuracy is not much im-proved by the regularity. It is not likely that a coarse regular mesh achieves thesame error magnitude as the double-sized irregular mesh. Although regularity failsto reduce discretization error, it still helps to decrease time costs in the solver. Thesmaller mesh sizes and faster convergence in the solver both contribute to this.However, for the second order solver, this improvement is not significant. Forhigher order solvers, the smaller time cost on the regular mesh is an advantage.Engwirda’s frontal-Deluany scheme and Marcum’s AFLR scheme both generate more regu-lar meshes than Delaunay meshes. For some simple geometries like an annulus or a square,AFLR scheme is able to generate almost perfect equilateral triangles everywhere. It is a goodchoice since it gives smallest truncation error and discretization error with reasonable addi-tional time cost in mesh generation. For more complicated numerical cases, especially for thehigh-order solver, Engwirda’s frontal-Delaunay is optimal because of its good regularity, gooderror performance, and decrease in time cost.83Bibliography[1] J. D. Anderson. Computational Fluid Dynamics: The Basics with Applications. McGraw-Hill Education, 1995.[2] John David Anderson and J Wendt. Computational Fluid Dynamics, volume 206.Springer, 1995.[3] Utkarsh Ayachit. The paraview guide: a parallel visualization application. 2015.[4] Charles Boivin and Carl Ollivier-Gooch. Guaranteed-quality triangular mesh generationfor domains with curved boundaries. International Journal for Numerical Methods inEngineering, 55(10):1185–1213, 2002.[5] L Paul Chew. Constrained Delaunay triangulations. Algorithmica, 4(1-4):97–108, 1989.[6] Boris Delaunay. Sur la sphere vide. Izv. Akad. Nauk SSSR, Otdelenie Matematicheskii iEstestvennyka Nauk, 7(793-800):1–2, 1934.[7] Boris Diskin and James L Thomas. Notes on accuracy of finite-volume discretizationschemes on irregular grids. Applied Numerical Mathematics, 60(3):224–226, 2010.[8] Boris Diskin and James L Thomas. Effects of mesh regularity on accuracy of finite-volume schemes. In 50th AIAA Aerospace Sciences Meeting, pages 2012–0609, 2012.[9] Darren Engwirda. Unstructured Tessellation & Mesh Generation. PhD thesis, Universityof Sydney, 2014.[10] Darren Engwirda and David Ivers. Face-centred Voronoi refinement for surface meshgeneration. Procedia Engineering, 82:8–20, 2014.84Bibliography[11] David A Field. Laplacian smoothing and Delaunay triangulations. Communications inApplied Numerical Methods, 4(6):709–712, 1988.[12] Lori Freitag, Mark Jones, and Paul Plassmann. A parallel algorithm for mesh smoothing.SIAM Journal on Scientific Computing, 20(6):2023–2040, 1999.[13] Lori A Freitag, Mark Jones, and Paul Plassmann. An efficient parallel algorithm formesh smoothing. In Proceedings of the Fourth International Meshing Roundtable, pages47–58, 1995.[14] Lori A Freitag and Carl Ollivier-Gooch. Tetrahedral mesh improvement using swap-ping and smoothing. International Journal for Numerical Methods in Engineering,40(21):3979–4002, 1997.[15] Alireza Jalali. Truncation error analysis of unstructured finite volume discretizationschemes. Master’s thesis, University of British Columbia, 2012.[16] Alireza Jalali and Carl Ollivier-Gooch. Accuracy assessment of finite volume discretiza-tions of convective fluxes on unstructured meshes. In 51st AIAA Aerospace SciencesMeeting including the New Horizons Forum and Aerospace Exposition, number 2013-0705.[17] Alireza Jalali and Carl Ollivier-Gooch. Accuracy assessment of finite volume discretiza-tions of diffusive fluxes on unstructured meshes. In 50th AIAA Aerospace Sciences Meet-ing including the New Horizons Forum and Aerospace Exposition, number 2012-0608.[18] Alireza Jalali, Mahkame Sharbatdar, and Carl Ollivier-Gooch. Accuracy analysis of un-structured finite volume discretization schemes for diffusive fluxes. Computers & Fluids,101:220–232, 2014.[19] F Juretic´ and AD Gosman. Error analysis of the finite-volume method with respect tomesh type. Numerical Heat Transfer, Part B: Fundamentals, 57(6):414–439, 2010.85Bibliography[20] Aaron Katz and Venkateswaran Sankaran. Mesh quality effects on the accuracy of CFDsolutions on unstructured meshes. Journal of Computational Physics, 230(20):7670–7686, 2011.[21] Harvard Lomax, Thomas H Pulliam, and David W Zingg. Fundamentals of Computa-tional Fluid Dynamics. Springer Science & Business Media, 2013.[22] David L Marcum and Nigel P Weatherill. Unstructured grid generation using iterativepoint insertion and local reconnection. AIAA journal, 33(9):1619–1625, 1995.[23] MATLAB. R2013b. The MathWorks Inc., Natick, Massachusetts, 2013.[24] Dimitri J Mavriplis. An advancing front Delaunay triangulation algorithm designed forrobustness. Journal of Computational Physics, 117(1):90–101, 1995.[25] J-D Müller, PL Roe, and H Deconinck. A frontal approach for internal node genera-tion in Delaunay triangulations. International Journal for Numerical Methods in Fluids,17(3):241–255, 1993.[26] Carl Ollivier-Gooch. A toolkit for numerical simulation of PDEs: I. Fundamentals ofgeneric finite-volume simulation. Computer Methods in Applied Mechanics and Engi-neering, 192(9):1147–1175, 2003.[27] Carl Ollivier-Gooch. GRUMMP version 0.7.0 user’s guide. University of BritishColumbia, 2016.[28] Carl Ollivier-Gooch and Charles Boivin. Guaranteed-quality simplical mesh generationwith cell size and grading control. Engineering with Computers, 17(3):269–286, 2001.[29] Carl Ollivier-Gooch and Michael Van Altena. A high-order-accurate unstructured meshfinite-volume scheme for the advection–diffusion equation. Journal of ComputationalPhysics, 181(2):729–752, 2002.86Bibliography[30] Olivier Pironneau. Finite Element Methods for Fluids. Wiley Chichester, 1989.[31] S Rebay. Efficient unstructured mesh generation by means of Delaunay triangulation andBowyer-Watson algorithm. Journal of Computational Physics, 106(1):125–138, 1993.[32] Christopher J Roy. Review of discretization error estimators in scientific computing.AIAA Paper, (2010-126).[33] Jim Ruppert. A Delaunay refinement algorithm for quality 2-dimensional mesh genera-tion. Journal of Algorithms, 18(3):548–585, 1995.[34] Jonathan R Shewchuk. Delaunay refinement mesh generation. Technical report, DTICDocument, 1997.[35] Alper Üngör. Off-centers: A new type of Steiner points for computing size-optimalquality-guaranteed Delaunay triangulations. In Latin American Symposium on Theoreti-cal Informatics, pages 152–161. Springer, 2004.[36] Michael Van Altena. High-order finite-volume discretisations for solving a modifiedadvection-diffusion problem on unstructured triangular meshes. PhD thesis, Universityof British Columbia, 1999.[37] H K Versteeg and W Malalasekera. An Introduction to Computational Fluid Dynamics,the Finite Volume Method. Longman Scientific & Technical, 1995.[38] Thomas Williams, Colin Kelley, and many others. Gnuplot 4.6: an interactive plottingprogram. http://gnuplot.sourceforge.net/, April 2013.[39] Gary Yan, Varun Prakash Puneria, Alireza Jalali, and Carl Ollivier-Gooch. Truncationand discretization error for diffusion schemes on unstructured meshes. AIAA SciTech.American Institute of Aeronautics and Astronautics, (2014-0478).[40] Daniel Zaide. Personal Email Communication.87


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items